text:
# This file is a part of RRRMC.jl. License is MIT: http://github.com/carlobaldassi/RRRMC.jl/LICENCE.md

module SAT

using Random
using ExtractMacro
using ..Interface
using ..Common

using ...RRRMC # this is silly but it's required for correct cross-linking in docstrings, apparently

export GraphSAT

import ..Interface: energy, delta_energy, neighbors, allΔE,
                    delta_energy_residual, update_cache!, update_cache_residual!

function choose(N::Int, K::Int)
    out = fill(typemax(Int), K)
    @inbounds for k = 1:K
        out[k] = rand(1:(N-k+1))
        for l = 1:(k-1)
            if out[l] ≤ out[k]
                out[k] += 1
            end
        end
        for l = 1:(k-1)
            if out[l] > out[k]
                x = out[k]
                for j = k:-1:l+1
                    out[j] = out[j-1]
                end
                out[l] = x
                break
            end
        end
    end
    #@assert issorted(out)
    #@assert all(diff(out) .≠ 0)
    return out
end

function gen_randomKSAT(N::Integer, K::Integer, α::Real)
    N > 0 || throw(ArgumentError("N must be positive: $N"))
    K > 0 || throw(ArgumentError("K must be positive: $K"))
    α ≥ 0 || throw(ArgumentError("α must be non-negative: $α"))
    N ≥ K || throw(ArgumentError("N must not be less than K: $N < $K"))

    M = round(Int, α * N)
    A = Array{IVec}(undef, M)
    J = [BitArray(undef, K) for a = 1:M]
    for a = 1:M
        A[a] = choose(N, K)
        rand!(J[a])
    end
    return A, J
end

struct ClauseCache
    M::Int
    S::IVec  # S[a] = how many vars satisfy clause a
    I::IVec2 # I[a] = indices of the vars in S[a]
    ClauseCache(M::Integer, K::Integer) = new(M, zeros(Int, M), IVec[zeros(Int, K) for a = 1:M])
    ClauseCache(M::Integer, A::IVec2) = new(M, zeros(Int, M), IVec[zeros(Int, length(A[a])) for a = 1:M])
end

function clear!(cache::ClauseCache)
    @extract cache : S I
    fill!(S, 0)
    for Ia in I
        fill!(Ia, 0)
    end
    return cache
end

struct GraphSAT <: DiscrGraph{Int}
    N::Int
    M::Int
    K::Int
    A::IVec2
    J::Vector{BitVector}
    T::IVec2
    neighb::IVec2
    max_conn::Int
    cache::ClauseCache
    lfcache::LocalFields{Int}
    function GraphSAT(N::Integer, A::IVec2, J::Vector{BitVector})
        M = length(A)
        length(J) == M || throw(ArgumentError("Incompatible lengths of A and J: $M vs $(length(J))"))

        K = maximum(map(length, A))

        T = [Int[] for i = 1:N]
        for a = 1:M
            for i in A[a]
                push!(T[i], a)
            end
        end

        neighb = [Int[] for i = 1:N]
        for i in 1:N
            Ti = T[i]
            for a in Ti, j in A[a]
                if j ≠ i && j ∉ neighb[i]
                    push!(neighb[i], j)
                end
            end
        end

        # TODO: more input consistency checks?
        max_conn = maximum(map(length, T))
        cache = ClauseCache(M, A)
        lfcache = LocalFields{Int}(N)
        return new(N, M, K, A, J, T, neighb, max_conn, cache, lfcache)
    end
end

"""
  GraphSAT(N::Integer, α::Real, K::Integer)

A `DiscrGraph` implementing a random `K`-SAT graph with `N` spins and `αN` clauses.

The energy of the model is the number of violated clauses.
"""
function GraphSAT(N::Integer, K::Integer, α::Real)
    A, J = gen_randomKSAT(N, K, α)
    return GraphSAT(N, A, J)
end

function export_cnf(X::GraphSAT, filename::AbstractString)
    @extract X : N M A J
    open(filename, "w") do f
        println(f, "p cnf $N $M")
        for a = 1:M
            for (j,i) in zip(J[a],A[a])
                print(f, (2j-1) * i, " ")
            end
            println(f, "0")
        end
    end
end

function export_cnf(X::GraphSAT, filename::AbstractString, decimate::Vector{Int})
    @extract X : N M A=deepcopy(A) J=deepcopy(J) T=deepcopy(T)

    j = 1
    while j ≤ length(decimate)
        v = decimate[j]
        s = v > 0
        i = abs(v)
        for a in T[i]
            isempty(A[a]) && continue
            k = findfirst(==(i), A[a])
            @assert k ≢ nothing
            if J[a][k] == s
                empty!(A[a])
            else
                length(A[a]) > 1 || @error("contradiction")
                deleteat!(A[a], k)
                deleteat!(J[a], k)
                if length(A[a]) == 1
                    newv = A[a][1] * (2J[a][1]-1)
                    -newv ∈ decimate && @error("contradiction")
                    newv ∉ decimate && push!(decimate, newv)
                    empty!(A[a])
                end
            end
        end
        empty!(T[i])
        j += 1
    end

    nM = sum(map(x->!isempty(x),A)) + length(decimate)

    open(filename, "w") do f
        println(f, "p cnf $N $nM")
        for a = 1:M
            isempty(A[a]) && continue
            for (j,i) in zip(J[a],A[a])
                print(f, (2j-1) * i, " ")
            end
            println(f, "0")
        end
        for v in decimate
            println(f, "$v 0")
        end
    end
end

function energy(X::GraphSAT, C::Config)
    length(C) == X.N || throw(ArgumentError("different N: $(length(C)) $(X.N)"))
    @extract C : s
    @extract X : N M A J T cache lfcache
    clear!(cache)
    @extract cache : S I
    @extract lfcache : lfields lfields_last

    n = 0
    @inbounds for a = 1:M
        Ja = J[a]
        Aa = A[a]
        sat = 0
        Ia = I[a]
        for k = 1:length(Aa)
            i = Aa[k]
            si = s[i]
            Ji = Ja[k]
            Ji ⊻ si == 0 && (sat += 1; Ia[sat] = i)
        end
        S[a] = sat
        sat == 0 && (n += 1)
    end

    for i = 1:N
        Δ = 0
        @inbounds for a in T[i]
            Sa = S[a]
            Ia = I[a]
            if Sa == 1 && (Ia[1] == i)
                Δ += 1
            elseif S[a] == 0
                Δ -= 1
            end
        end
        lfields[i] = -Δ
    end

    lfcache.move_last = 0
    fill!(lfields_last, 0)

    return n
end

function delta_energy(X::GraphSAT, C::Config, move::Int)
    @assert X.N == C.N
    @assert 1 ≤ move ≤ C.N
    #@extract C : s
    @extract X : lfcache
    @extract lfcache : lfields

    @inbounds Δ = -lfields[move]
    return Δ

    # @inbounds begin
    #     Δ = 0.0
    #     Jx = J[move]
    #     σx = 2 * s[move] - 1
    #     Ax = A[move]
    #     for k = 1:length(Ax)
    #         y = Ax[k]
    #         σy = 2 * s[y] - 1
    #         Jxy = Jx[k]
    #         Δ += 2 * Jxy * σx * σy
    #     end
    # end
    # return Δ
end

function update_cache!(X::GraphSAT, C::Config, move::Int)
    @assert X.N == C.N
    @assert 1 ≤ move ≤ C.N
    @extract C : N s
    @extract X : A T neighb cache lfcache

    @extract cache : S I
    @extract lfcache : lfields lfields_last move_last

    # if move_last == move
    #     @inbounds begin
    #         for j in neighbors[move]
    #             lfields[j], lfields_last[j] = lfields_last[j], lfields[j]
    #         end
    #         lfields[move] = -lfields[move]
    #         lfields_last[move] = -lfields_last[move]
    #     end
    #     return
    # end

    @inbounds for a in T[move]
        Sa = S[a]
        Ia = I[a]
        Aa = A[a]
        if Sa == 0
            S[a] = 1
            Ia[1] = move
            lfields[move] -= 2
            for j in Aa
                j == move && continue
                lfields[j] -= 1
            end
        else
            k = findfirst(==(move), Ia)
            if k ≢ nothing
                for l = k:Sa-1
                    Ia[l] = Ia[l+1]
                end
                Ia[Sa] = 0
                S[a] = Sa - 1
                if Sa == 1
                    lfields[move] += 2
                    for j in Aa
                        j == move && continue
                        lfields[j] += 1
                    end
                elseif Sa == 2
                    lfields[Ia[1]] -= 1
                end
            else
                if Sa == 1
                    lfields[Ia[1]] += 1
                end
                S[a] = Sa + 1
                Ia[Sa + 1] = move
            end
        end
    end

    # cache.move_last = move

    return
end

neighbors(X::GraphSAT, i::Int) = return X.neighb[i]

# TODO: improve
allΔE(X::GraphSAT) = tuple(0:X.max_conn...)

end


text:
using ThreeBodyDecay
using LaTeXStrings
using Plots

const mDstar = 2.00685 # 2 GeV
const ΓDstar = 30e-6 # 30 keV
tbs = let mD = 1.86483, mπ0 = 0.13498
    ThreeBodySystem(mD,mπ0,mD,mD+mDstar+0.01)
end

# amplitude

BWD(σ) = BW(σ,mDstar,100*ΓDstar)
A(σ3,σ1) = BWD(σ1) + BWD(σ3)
A(σ3,σ1,ν) = BWD(σ1) * wignerd(1, ν, 0, cosθ23(gσ2(σ3,σ1,tbs),σ3,tbs)) -
             BWD(σ3) * sum(wignerd(1, ν, λ, cos_plus_θhat3(σ1,gσ2(σ3,σ1,tbs),tbs)) *
                           wignerd(1, λ, 0, cosθ12(σ1,gσ2(σ3,σ1,tbs),tbs)) for λ=-1:1)

# intensity

I(σ3,σ1) = sum(abs2(A(σ3,σ1,ν)) for ν=-1:1)

let
    σ1v = LinRange(tbs.mthsq[1], tbs.sthsq[1],200)
    σ3v = LinRange(tbs.mthsq[3], tbs.sthsq[3],200)
    cal = [Kibble31(σ3,σ1,tbs) < 0.0 ? I(σ3,σ1) : NaN for σ3 in σ3v, σ1 in σ1v]
    heatmap(σ1v, σ3v, cal)
end

# sample density
density2d = getbinned2dDensity(
                (σ3,σ1) -> Kibble31(σ3,σ1,tbs) < 0.0 ? I(σ3,σ1) : 0.0,
                    (tbs.mthsq[1], tbs.sthsq[1]),
                    (tbs.mthsq[3], tbs.sthsq[3]), 100, 100);

let r = hcat([rand(density2d) for _ in 1:100000]...)
    histogram2d(r[2,:], r[1,:], bins=100)
end

ext: jl

max_stars_repo_name: mmikhasenko/ThreeBodyDecay

max_issues_repo_name: mmikhasenko/ThreeBodyDecay

max_forks_repo_name: mmikhasenko/ThreeBodyDecay


text:
module MSQR

using Yao
using QMPS
using Flux
using LinearAlgebra
using Statistics
using Requires

include("MPSSwapTest.jl")
include("MSQRtrain.jl")
include("SWAPtrain.jl")
@init @require CuYao="b48ca7a8-dd42-11e8-2b8e-1b7706800275" include("CUDA.jl")

end


text:
for op = (:+, :*, :&, :|, :$)
  @eval ($op)(a,b,c) = ($op)(($op)(a,b),c)
end

function g(x,y)
  return x * y
  x + y
end

cd("data") do
    open("outfile", "w") do f
        write(f, data)
    end
end

max_stars_repo_name: alejandro-medici/ace

max_issues_repo_name: GistBox/ace

max_forks_repo_name: GistBox/ace


text:
fusiontreedict(I) = FusionStyle(I) isa Abelian ? SingletonDict : FusionTreeDict

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

Perform an elementary braid (Artin generator) of neighbouring uncoupled indices `i` and
`i+1` on a fusion 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 fusion
tree with non-zero coefficient, namely `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′ = TupleTools.setindex(f.isdual, f.isdual[i], i+1)
    isdual′ = TupleTools.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′ = TupleTools.setindex(uncoupled, b, 1)
        uncoupled′ = TupleTools.setindex(uncoupled′, a, 2)
        if FusionStyle(I) isa Union{Abelian, SimpleNonAbelian}
            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 # DegenerateNonAbelian
            μ = 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′ = TupleTools.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′ = TupleTools.setindex(uncoupled, d, i)
    uncoupled′ = TupleTools.setindex(uncoupled′, b, i+1)
    if FusionStyle(I) isa Abelian
        inner′ = TupleTools.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 SimpleNonAbelian
        local newtrees
        for c′ in intersect(a ⊗ d, e ⊗ conj(b))
            coeff = oftype(oneT, if inv
                    conj(Rsymbol(d, c, e)*Fsymbol(d, a, b, e, c′, c))*Rsymbol(a, d, c′)
                else
                    Rsymbol(c, d, e)*conj(Fsymbol(d, a, b, e, c′, c)*Rsymbol(d, a, c′))
                end)
            iszero(coeff) && continue
            inner′ = TupleTools.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 # DegenerateNonAbelian
        local newtrees
        for c′ in intersect(a ⊗ d, e ⊗ conj(b))
            Rmat1 = inv ? Rsymbol(d, c, e)' : Rsymbol(c, d, e)
            Rmat2 = inv ? Rsymbol(a, d, c′) : Rsymbol(d, a, c′)'
            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′ = TupleTools.setindex(vertices, σ, i-1)
                    vertices′ = TupleTools.setindex(vertices′, λ, i)
                    inner′ = TupleTools.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

# braid fusion tree
"""
    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}
    TupleTools.isperm(p) || throw(ArgumentError("not a valid permutation: $p"))
    if FusionStyle(I) isa Abelian && BraidingStyle(I) isa SymmetricBraiding
        coeff = Rsymbol(one(I), one(I), one(I))
        for i = 1:N
            for j = 1:i-1
                if p[j] > p[i]
                    a, b = f.uncoupled[j], f.uncoupled[i]
                    coeff *= Rsymbol(a, b, first(a ⊗ b))
                end
            end
        end
        uncoupled′ = TupleTools._permute(f.uncoupled, p)
        coupled′ = f.coupled
        isdual′ = TupleTools._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 = TupleTools.setindex(levels, levels[s+1], s)
            levels = TupleTools.setindex(levels, l, s+1)
            trees, newtrees = newtrees, trees
            empty!(newtrees)
        end
        return trees
    end
end

# permute fusion tree
"""
    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

"""
    split(f::FusionTree{I, N}, M::Int)
    -> (::FusionTree{I, M}, ::FusionTree{I, N-M+1})

Split a fusion tree with the first M outgoing indices, and an incoming index corresponding
to the internal fusion tree index between outgoing indices N and N+1 of the original tree
`f`; and a second fusion tree whose first outgoing index is that same internal index. Its
remaining outgoing indices are the N-M outgoing indices of the original tree `f`, and also
the incoming index is the same. This is in the inverse of `insertat` in the sense that if
`f1, f2 = split(t, M) ⇒ f == insertat(f2, 1, f1)`.
"""
@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
        (FusionTree{I}((f.uncoupled[1],), f.uncoupled[1], (false,), (), ()), f)
    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 DegenerateNonAbelian
            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

"""
    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) == DegenerateNonAbelian()`, 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 DegenerateNonAbelian && μ === 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))
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 standard 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]

    uncoupled = TupleTools.deleteat(f1.uncoupled, i)
    coupled = f1.coupled
    isdual = TupleTools.deleteat(f1.isdual, i)
    if length(uncoupled) <= 2
        inner = ()
    else
        inner = TupleTools.deleteat(f1.innerlines, max(1, i-2))
    end
    if length(uncoupled) <= 1
        vertices = ()
    else
        vertices = TupleTools.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, 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′ = TupleTools.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, 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′ = TupleTools.insertafter(TupleTools.setindex(uncoupled, b, i), i, (c,))
    isdual′ = TupleTools.insertafter(TupleTools.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 Union{Abelian, SimpleNonAbelian}
        local newtrees
        for e in a ⊗ b
            coeff = conj(Fsymbol(a, b, c, d, e, e′))
            iszero(coeff) && continue
            inner′ = TupleTools.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′ = TupleTools.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′ = TupleTools.setindex(f1.vertices, ν, i-1)
                vertices′ = TupleTools.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, 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]
    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′ => coeff*coeff′)
                end
            end
        end
        return newtrees
    end
end

# repartition double fusion tree
"""
    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}

Input is a double fusion 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 the tree by bending incoming to outgoing sectors (or vice versa) in order to
have `N` outgoing sectors.
"""
@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
    # DegenerateNonAbelian, 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
        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

# 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
#     # DegenerateNonAbelian, 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 FusionStyle(I) isa Union{Abelian, SimpleNonAbelian}
#         if N == N₁
#             return SingletonDict{Tuple{F1,F2},T}( (f1, f2) => coeff)
#         else
#             (f1′, f2′), coeff1 = first(N < N₁ ? _bendright(f1, f2) : _bendleft(f1, f2))
#             (f1′′, f2′′), coeff2 = first(_recursive_repartition(f1′, f2′, Val(N)))
#             return SingletonDict{Tuple{F1,F2},T}( (f1′′, f2′′) => coeff1*coeff2)
#         end
#     else
#         if N == N₁
#             return FusionTreeDict{Tuple{F1,F2},T}( (f1, f2) => coeff)
#         else
#             newtrees = FusionTreeDict{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))
#                     push!(newtrees, (f1′′, f2′′) => coeff1*coeff2)
#                 end
#             end
#             return newtrees
#         end
#     end
# end

# change to N₁ - 1, N₂ + 1
function _bendright(f1::FusionTree{I, N₁}, f2::FusionTree{I, N₂}) where {I<:Sector, N₁, N₂}
    # map final splitting vertex (a, b)<-c to fusion vertex a<-(c, dual(b))
    @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 Union{Abelian, SimpleNonAbelian}
        if FusionStyle(I) isa Abelian
            sqrtca = 1
        else
            sqrtca = sqrt(dim(c)/dim(a))
        end
        coeff = sqrtca*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 = sqrt(dim(c)/dim(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
# change to N₁ + 1, N₂ - 1
function _bendleft(f1::FusionTree{I}, f2::FusionTree{I}) where I
    # map final fusion vertex c<-(a, b) to splitting vertex (c, dual(b))<-a
    return fusiontreedict(I)((f1′, f2′) => conj(coeff) for
                                ((f2′, f1′), coeff) in _bendright(f2, f1))
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 pair 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 TupleTools.isperm((p1..., p2...))
    if FusionStyle(f1) isa Abelian &&
        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(sqrt(dim(u))*one(eltype(Fsymbol(u, u, u, u, u, u)))*one(eltype(Rsymbol(u, u, u))))
            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 double fusion 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 (`t1`) and incoming sectors (`t2`) respectively (with identical coupled sector
`t1.coupled == t2.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.
"""
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

# transpose double fusion tree
const transposecache = LRU{Any, Any}(; maxsize = 10^5)
const usetransposecache = Ref{Bool}(true)

"""
    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}

Input is a double fusion 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 (`t1`) and incoming sectors (`t2`) respectively (with identical coupled sector
`t1.coupled == t2.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.
"""
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))
    i1 = findfirst(==(1), p)
    @assert i1 !== nothing
    newtrees = repartition(f1, f2, N₁)
    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

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

# change to N₁ - 1, N₂ + 1
function _foldright(f1::FusionTree{I, N₁}, f2::FusionTree{I, N₂}) where {I<:Sector, N₁, N₂}
    # map first splitting vertex (a, b)<-c to fusion vertex b<-(dual(a), c)
    @assert N₁ > 0
    if FusionStyle(I) isa Abelian
        a = f1.uncoupled[1]
        isduala = f1.isdual[1]
        sqrtda = sqrt(dim(a))
        factor = frobeniusschur(a)
        if isduala
            factor = one(factor)
        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, (!f1.isdual[1], f2.isdual...))
        return fusiontreedict(I)((fl,fr)=>1)
    else
        a = f1.uncoupled[1]
        isduala = f1.isdual[1]
        sqrtda = sqrt(dim(a))
        factor = frobeniusschur(a)
        if isduala
            factor = one(factor)
        end
        c1 = dual(a)
        c2 = f1.coupled
        hasdegeneracies = FusionStyle(a) isa DegenerateNonAbelian
        local newtrees
        for c in c1 ⊗ c2
            N₁ == 1 && c != one(c) && continue
            for μ in (hasdegeneracies ? (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 * sqrtda * 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
# change to N₁ + 1, N₂ - 1
function _foldleft(f1::FusionTree{I}, f2::FusionTree{I}) where I
    # map first fusion vertex c<-(a, b) to splitting vertex (dual(a), c)<-b
    return fusiontreedict(I)((f1′, f2′) => conj(coeff) for
                                    ((f2′, f1′), coeff) in _foldright(f2, f1))
end

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

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


text:
using Distributions
using SpecialFunctions

function fim(d::Normal)
    μ,σ = params(d)
    [1 / σ^2 0.0; 0.0 2/σ^2]
end
function fim(d::Beta)
    α,β = params(d)
    [trigamma(α) - trigamma(α + β)  -trigamma(α + β);
     -trigamma(α + β)   trigamma(β) - trigamma(α + β)]
end
function fim(d::Gamma)
    α,θ = params(d);
    [polygamma(1,α) θ; θ  α*θ^2]
end

function fim(d::Product)
    BlockDiagonal([fim(marginal) for marginal in d.v])
end

max_stars_repo_name: ap-browning/FisherMetrics

max_issues_repo_name: ap-browning/FisherMetrics

max_forks_repo_name: ap-browning/FisherMetrics


text:
using CovidPOMDPs
using POMCPOW
using Test

macro isinferred(ex)
  quote
    try
      @inferred $(Expr(ex.head, esc.(ex.args)...))
      true
    catch err
      println(err)
      false
    end
  end
end

function validate_history(hist::SimHist)
    S = susceptible(hist)
    I = infected(hist)
    R = recovered(hist)
    N = hist.N
    T = hist.T

    @test all(≥(0), S)
    @test all(≥(0), I)
    @test all(≥(0), R)
    @test all(≥(0), hist.pos_test)

    # Population size preservation
    @test all(==(N), S[i] + I[i] + R[i] for i in 1:T)

    @test all( 0.0 .≤ getfield.(hist.actions,:testing_prop) .≤ 1.0)

    return nothing
end

@testset "Common" begin
    include("common.jl")
end

@testset verbose=true "Single Strain" begin
    include(joinpath("single_strain", "single_strain.jl"))
end

@testset verbose=true "Double Strain" begin
    include(joinpath("double_strain", "double_strain.jl"))
end


text:
#! /usr/bin/env julia
lines = map(x -> parse(Int, x), readlines())

i = 1
value = 0;
visited = BitSet()
while !(value in visited)
  push!(visited, value)
  value += lines[i]
  if i >= length(lines)
    i = 0
  end
  i += 1
end
println(value)

max_stars_repo_name: paolobueno/advent-of-julia

max_issues_repo_name: paolobueno/advent-of-julia

max_forks_repo_name: paolobueno/advent-of-julia


text:
n_obs(sem_fit) = n_obs(sem_fit.model)

n_obs(model::AbstractSemSingle) = n_obs(model.observed)

n_obs(model::SemEnsemble) = sum(n_obs.(model.sems))

#= n_obs(sem_fit::SemFit{Mi, So, St, Mo, O} where {Mi, So, St, Mo <: AbstractSemSingle, O}) =
    n_obs(
        sem_fit,
        sem_fit.model.observed,
        sem_fit.model.imply,
        sem_fit.model.diff,
        sem_fit.model.loss.functions...
        )

# RAM + SemML
n_obs(sem_fit::SemFit, obs::Union{SemObsCommon, SemObsMissing}, imp, diff, loss_ml) =
    n_obs(sem_fit.model)

n_obs(model::AbstractSemSingle) = n_obs(model.observed) =#


text:
function Base.show(io::IO, t::MIME"text/plain", x::Union{DecisionTrees.Tree,DecisionTrees.Leaf,DecisionTrees.Node,DecisionTrees.BaggedTree,DecisionTrees.BoostedTree})
    elt = string(typeof(x))
    println(io, elt)
    # println(io,"")
    show(io, t, fieldnames(typeof(x)))
    # println(io,"")
    # println(io,"")
    # println(io,"Use the dump() function to see the whole content of the object.")
end

function Base.show(io::IO, x::Union{DecisionTrees.Tree,DecisionTrees.Leaf,DecisionTrees.Node,DecisionTrees.BaggedTree,DecisionTrees.BoostedTree})
    elt = string(typeof(x))
    # println(io,"")
    println(io, elt)
    # println(io,"")
    show(io, fieldnames(typeof(x)))
    # println(io,"")
    # println(io,"")
    # println(io,"Use the dump() function to see the whole content of the object.")
end


text:
module TestGrouping
    using Base.Test
    using DataFrames

    df = DataFrame(a=rep(1:4, 2), b=rep(2:-1:1, 4), c=randn(8))
    #df[6, :a] = NA
    #df[7, :b] = NA

    cols = [:a, :b]

    sdf = sort(df, cols=cols)
    bdf = by(df, cols, df -> DataFrame(cmax = maximum(df[:c])))

    @test isequal(bdf[cols], unique(sdf[cols]))

    byf = by(df, :a, df -> DataFrame(bsum = sum(df[:b])))
end


text:
# NUMBER FORMAT OPTIONS
const Numtype = Float32
#const Numtype = Posit{16,2}
#const Numtype = Main.FiniteFloats.Finite16
#const Numtype = BigFloat
#setprecision(7)

# DOMAIN RESOLUTION AND RATIO
const nx = 100                  # number of grid cells in x-direction
const Lx = 1000e3               # length of the domain in x-direction
const L_ratio = 1               # Domain aspect ratio of Lx/Ly

# PHYSICAL CONSTANTS
const gravity = 10.             # gravitational acceleration
const water_depth = 500.        # layer thickness at rest
const ρ = 1e3                   # density
const ϕ = 30.                    # central latitue of the domain (for coriolis)
const ω = 2π/(24*3600)          # Earth's angular frequency [s^-1]
const R = 6.371e6               # Earth's radius [m]

# WIND FORCING OPTIONS
const wind_forcing = "channel"  # "channel", "double_gyre", "shear" or "none"
const Fx0 = 0.12                # wind stress strength [Pa], default 0.12

# BOTTOM TOPOGRAPHY OPTIONS
const topography_feature = "ridges" # "ridge", "seamount", "flat", "ridges"
const topofeat_height = 80.      # height of seamount
const topofeat_width = 200e3    # horizontal scale [m] of the seamount

# NEWTONIAN COOLING OPTIONS
const surface_forcing = false   # or true
const t_relax = 100.              # time scale of the interface_relaxation [days]
const η_refh = 5.               # height difference [m] of the interface relaxation profile
const η_refw = 100e3             # width [m] of the tangent used for the interface relaxation

# TIME STEPPING OPTIONS
const RKo = 4                   # Order of the RK time stepping scheme (3 or 4)
const cfl = 1.0                 # CFL number (1.0 recommended for RK4, 0.6 for RK3)
const Ndays = 100                # number of days to integrate for
const nstep_diff = 10           # diffusive part every nstep_diff time steps.
const nstep_advcor = 1         # advection and coriolis every nstep_advcor time steps.

# BOUNDARY CONDITION OPTIONS
const bc_x = "periodic"      # "periodic" or anything else for nonperiodic
const lbc = 1.                  # lateral boundary condition parameter
                                # 0 free-slip, 0<lbc<2 partial-slip, 2 no-slip

# MOMENTUM ADVECTION OPTIONS
const adv_scheme = "ArakawaHsu"   # "Sadourny" or "ArakawaHsu"

# BOTTOM FRICTION OPTIONS
const bottom_friction = "linear" # "linear" or "quadratic"
const drag = 1e-5               # bottom drag coefficient [dimensionless] for quadratic
const τdrag = 300.               # bottom drag coefficient [days] for linear

# DIFFUSION OPTIONS
const diffusion = "Constant"    # "Smagorinsky" or "Constant", biharmonic in both cases
const ν_const = 500             # [m^2/s] scaling constant for Constant biharmonic diffusion
const c_smag = 0.15             # Smagorinsky coefficient [dimensionless]

# TRACER ADVECTION
const tracer_advection = false   # "true" or "false"
const tracer_relaxation = false  # "true" or "false"
const injection_region = "rect"  # "west", "south" or "rect"
const sstrestart = false        # start from previous sst file
const Uadv = 0.5               # Velocity scale [m/s] for tracer advection
const SSTmax = 30.               # tracer (sea surface temperature) max for restoring
const SSTmin = 0.               # tracer (sea surface temperature) min for restoring
const τSST = 500.                # tracer restoring time scale [days]
const SSTw = 1000e3               # width [m] of the tangent used for the IC and interface relaxation
const SSTϕ = 0.5                # latitude/longitude ∈ [0,1] of sst edge

# OUTPUT OPTIONS
const output = true                # for nc output
const output_vars = ["u","sst"]
const output_dt = 6             # output time step in hours
#const outpath = "/network/aopp/chaos/pred/kloewer/julsdata/ssthr/"
const outpath = "/local/kloewer/julsdata/timest/"

# INITIAL CONDITIONS
const initial_cond = "rest"   # "rest" or "ncfile"
const initpath = "/network/aopp/chaos/pred/kloewer/julsdata/ssthr/"

const init_run_id = 2           # only for starting from ncfile

max_stars_repo_name: eschnett/juls

max_issues_repo_name: eschnett/juls

max_forks_repo_name: eschnett/juls


text:
module ExtensibleMCMC

    using StaticArrays, Distributions, Random, DataStructures, LinearAlgebra
    using Dates

    import Distributions: logpdf

    include("types.jl") # ✔

    const __PREVIOUS = Previous()
    const __PROPOSAL = Proposal()
    const __PRESTEP = PreMCMCStep()
    const __POSTSTEP = PostMCMCStep()

    include("utility_functions.jl") # ✔
    include("priors.jl")
    include("schedule.jl") # ✔
    include("chain_statistics.jl") # ✔
    _DIR = "transition_kernels"
    include(joinpath(_DIR, "random_walk.jl")) # ✔
    include(joinpath(_DIR, "adaptation.jl")) # ✔

    include("callbacks.jl") # ✔
    include("mcmc.jl") # ✔
    include("run.jl") # ✔
    include("updates.jl") # ✔/✗ (TODO implement conjugate updates, MALA, Hamiltionan Monte Carlo, and much more, this might have do be done after DiffusionMCMC.jl)
    include("workspaces.jl") # ✔
    _DIR = "example"
    include(joinpath(_DIR, "gsn_target.jl")) # ✔

    export MCMC
    export UniformRandomWalk, GaussianRandomWalk, GaussianRandomWalkMix
    export AdaptationUnifRW, HaarioTypeAdaptation, NoAdaptation
    export RandomWalkUpdate
    export GenericMCMCBackend
    export GenericChainStats
    export GsnTargetLaw
    export run!
    export get_decorators, isdecorator
    export ImproperPosPrior, ImproperPrior
    export SavingCallback, REPLCallback
end # module


text:
export writexml

####### WRITEXML ###########
writexml(node, f) = writexml(stdout, node, f)

function writexml(io::IO, doc::Document,f::String)
    f=open(f, "w")
    write(f,"""<?xml version="1.0" encoding="UTF-8"?> \n""")
    if hasroot(doc)
        f=writexml(io, root(doc), f,0)
    else
        @warn "No root found"
    end
    close(f)
end

function writexml(io::IO, n::TextNode,f::IOStream)
    write(f,n.content)
    return f
end

function writexml(io::IO, n::Node,f::IOStream, depth::Integer = 0)
    @warn "Unknown node type"
end

function writexml(io::IO, n::AttributeNode,f::IOStream)
   write(f,n.name*"=\""*n.value*"\"")
   return f
end

function writexml(io::IO, n::TextNode,f::IOStream, depth::Integer)
    write(f,"    "^depth*n.content * " \n")
    return f
end

function writexml(io::IO, parent::ElementNode,f::IOStream, depth::Integer = 0)
    tag = nodename(parent)
    write(f,"    "^depth*"<$tag")

    attrs = map(x -> x.name * "=\"$(x.value)\"", attributes(parent))
    attr_str = join(attrs, " ")
    if !isempty(attr_str)
        write(f," "*attr_str)
    end

    children = nodes(parent)
    len = length(children)

    if len == 0
        write(f,"/> \n")
    elseif len == 1 && istext(first(children))
        write(f,">")
        for n in children
            f=writexml(io, n,f)
        end
    else
        write(f,"> \n")
        for n in children
            f=writexml(io, n,f,depth + 1)
        end
        write(f,"    "^depth)
    end

    if len != 0
        write(f,"</$tag> \n")
    end

    return f
end


text:
using ImageInspector, MLDatasets, Plots

# gray images
x = FashionMNIST.traintensor(1);
plot(
    plot(image(x; flip = true); title = "flip = true"),
    plot(image(x; flip = false); title = "flip = false");
    axis = nothing,
    border = :none,
)

# color images
x1 = FashionMNIST.traintensor(1);
x2 = CIFAR10.traintensor(2);
plot(
    plot(image(x1)),
    plot(image(x2));
    axis = nothing,
    border = :none
)

# imagegrid
xs = FashionMNIST.traintensor(1:10);
plot(imagegrid(xs, 1:10; nrows = 2, sep = 2); axis = nothing, border = :none)

# imageplot
x = CIFAR10.traintensor(1:10);
imageplot(x, 1:10; nrows = 2, sep = 1)


text:
"""
    type Coord

An {x,y,z} coordinate type. Used throughout the ConstructiveSolidGeometry.jl package for speed.

# Constructors
* `Coord(x::Float64, y::Float64, z::Float64)`
"""
type Coord
    x::Float64
    y::Float64
    z::Float64
end

"""
    type Ray

A ray is defined by its origin and a unitized direction vector

# Constructors
* `Ray(origin::Coord, direction::Coord)`
"""
type Ray
    origin::Coord
    direction::Coord
end

"""
    abstract Surface

An abstract class that all surfaces (`Sphere`, `Plane`, `InfCylinder`) inherit from. Implementation of new shapes should inherit from `Surface`.
"""
abstract type Surface end

"""
    type Plane <: Surface

Defined by a point on the surface of the plane, its unit normal vector, and an optional boundary condition.

# Constructors
* `Plane(point::Coord, normal::Coord)`
* `Plane(point::Coord, normal::Coord, boundary::String)`

# Arguments
* `point::Coord`: Any point on the surface of the plane
* `normal::Coord`: A unit normal vector of the plane. Recommended to use `unitize(c::Coord)` if normalizing is needed.
* `boundary::String`: Optional boundary condition, defined as a `String`. Options are "transmission" (default), "vacuum", and "reflective".
"""
type Plane <: Surface
    point::Coord
    normal::Coord
    reflective::Bool
	vacuum::Bool
	Plane(point::Coord, normal::Coord, ref::Bool, vac::Bool) = new(point, normal, ref, vac)
	function Plane(point::Coord, normal::Coord, boundary::String)
		if boundary == "reflective"
			new(point, normal, true, false)
		elseif boundary == "vacuum"
			new(point, normal, false, true)
		else
			new(point, normal, false, false)
		end
	end
	Plane(point::Coord, normal::Coord) = new(point, normal, false, false)
end

"""
    type Cone <: Surface

Defined by the tip of the cone, its direction axis vector, the angle between the central axis and the cone surface, and an optional boundary condition.

# Constructors
* `Cone(tip::Coord, axis::Coord, theta::Float64)`
* `Cone(tip::Coord, axis::Coord, theta::Float64, boundary::String)`

# Arguments
* `tip::Coord`: The vertex (tip) of the cone
* `axis::Coord`: A unit vector representing the central axis of the cone. As the cone equation actually defines two cones eminating in a mirrored fashion from the tip, this direction vector also indicates which cone is the true cone. I.e., following the direction of the axis vector when starting from the tip should lead inside the cone you actually want.  Recommended to use `unitize(c::Coord)` if normalizing is needed.
* `theta::Float64`: The angle (in radians) between the central axis (must be between 0 and pi/2)
* `boundary::String`: Optional boundary condition, defined as a `String`. Options are \"transmission\" (default) or \"vacuum\".
"""
type Cone <: Surface
    tip::Coord
    axis::Coord
	theta::Float64
    reflective::Bool
	vacuum::Bool
	Cone(tip::Coord, axis::Coord, theta::Float64, ref::Bool, vac::Bool) = new(tip, axis, theta, ref, vac)
	function Cone(tip::Coord, axis::Coord, theta::Float64, boundary::String)
		if boundary == "reflective"
			new(tip, axis, theta, true, false)
		elseif boundary == "vacuum"
			new(tip, axis, theta, false, true)
		else
			new(tip, axis, theta, false, false)
		end
	end
	Cone(tip::Coord, axis::Coord, theta::Float64) = new(tip, axis, theta, false, false)
end

"""
    type Sphere <: Surface

Defined by the center of the sphere, its radius, and an optional boundary condition.

# Constructors
* `Sphere(center::Coord, radius::Float64)`
* `Sphere(center::Coord, radius::Float64, boundary::String)`

# Arguments
* `center::Coord`: The center of the sphere
* `radius::Float64`: The radius of the sphere
* `boundary::String`: Optional boundary condition, defined as a `String`. Options are \"transmission\" (default) or \"vacuum\".
"""
type Sphere <: Surface
    center::Coord
    radius::Float64
    reflective::Bool
	vacuum::Bool
	Sphere(c::Coord, r::Float64, ref::Bool, vac::Bool) = new(c, r, ref, vac)
	function Sphere(c::Coord, r::Float64, boundary::String)
		if boundary == "reflective"
			new(c, r, true, false)
		elseif boundary == "vacuum"
			new(c, r, false, true)
		else
			new(c, r, false, false)
		end
	end
	Sphere(c::Coord, r::Float64) = new(c, r, false, false)
end

"""
    type InfCylinder <: Surface

An arbitrary direction infinite cylinder defined by any point on its central axis, its radius, the unit normal direction of the cylinder, and an optional boundary condition. A finite cylinder can be generated by defining the intersection of an infinite cylinder and two planes.

# Constructors
* `InfCylinder(center::Coord, normal::Coord, radius::Float64)`
* `InfCylinder(center::Coord, normal::Coord, radius::Float64, boundary::String)`

# Arguments
* `center::Coord`: The center of the infinite cylinder
* `normal::Coord`: A unit normal direction vector of the cylinder (i.e., a vector along its central axis), Recommended to use `unitize(c::Coord)` if normalizing is needed.
* `radius::Float64`: The radius of the infinite cylinder
* `boundary::String`: Optional boundary condition, defined as a `String`. Options are \"transmission\" (default) or \"vacuum\".
"""
type InfCylinder <: Surface
    center::Coord
    normal::Coord
    radius::Float64
    reflective::Bool
	vacuum::Bool
	InfCylinder(c::Coord, n::Coord, r::Float64, ref::Bool, vac::Bool) = new(c, n, r, ref, vac)
	function InfCylinder(c::Coord, n::Coord, r::Float64, boundary::String)
		if boundary == "reflective"
			new(c, n, r, true, false)
		elseif boundary == "vacuum"
			new(c, n, r, false, true)
		else
			new(c, n, r, false, false)
		end
	end
	InfCylinder(c::Coord, n::Coord, r::Float64) = new(c, n, r, false, false)
end

"""
    type Box

An axis aligned box is defined by the minimum `Coord` and maximum `Coord` of the box. Note that a Box is only used by ConstructiveSolidGeometry.jl for bounding box purposes, and is not a valid surface to define CSG cells with. Instead, you must define all six planes of a box independently.

# Constructors
* `Box(min::Coord, max::Coord)`
"""
type Box
    lower_left::Coord
    upper_right::Coord
end

"""
    type Region

The volume that is defined by a surface and one of its halfspaces

# Constructors
* `Region(surface::Surface, halfspace::Int64)`

# Arguments
* `surface::Surface`: A `Sphere`, `Plane`, or `InfCylinder`
* `halfspace::Int64`: Either +1 or -1
"""
type Region
    surface::Surface
    halfspace::Int64
end

"""
    type Cell

Defined by an array of regions and the logical combination of those regions that define the cell

# Constructors
* `Cell(regions::Array{Region}, definition::Expr)`

# Arguments
* `regions::Array{Region}`: An array of regions that are used to define the cell
* `definition::Expr`: A logical expression that defines the volume of the cell. The intersection operator is ^, the union operator is |, and the complement operator is ~. Regions are defined by their integer indices in the regions array.
"""
type Cell
    regions::Array{Region}
    definition::Expr
end

"""
    type Geometry

The top level object that holds all the cells in the problem. This object contains all data regarding the geometry within a system.

# Constructors
* `Geometry(cells::Array{Cell}, bounding_box::Box)`

# Arguments
* `cells::Array{Cell}`: All cells inside the geometry. The cells must combine to fill the entire space of the bounding box. No two cells should overlap.
* `bounding_box::Box`: The bounding box around the problem.
"""
type Geometry
    cells::Array{Cell}
    bounding_box::Box
end

_p = Coord(0,0,0)
typeassert(_p, Coord)


text:
# Julia wrapper for header: cuda_device_runtime_api.h
# Automatically generated using Clang.jl

# Not fully automated, since Clang.jl can't handle the `extern "C"`, and we
# change the `ccall` invocation to refer to a manually-provided library (`extern ...`).
# We also provide 1.0-compatible definitions for functionality that is officially supported.

# TODO: checking

import Base.Sys: WORD_SIZE

function cudaDeviceGetAttribute(value, attr, device)
    ccall("extern cudaDeviceGetAttribute", llvmcall, cudaError_t,
          (Ptr{Cint}, cudaDeviceAttr, Cint),
          value, attr, device)
end

function cudaDeviceGetLimit(pValue, limit)
    ccall("extern cudaDeviceGetLimit", llvmcall, cudaError_t,
          (Ptr{Csize_t}, cudaLimit),
          pValue, limit)
end

function cudaDeviceGetCacheConfig(pCacheConfig)
    ccall("extern cudaDeviceGetCacheConfig", llvmcall, cudaError_t,
          (Ptr{cudaFuncCache},),
          pCacheConfig)
end

function cudaDeviceGetSharedMemConfig(pConfig)
    ccall("extern cudaDeviceGetSharedMemConfig", llvmcall, cudaError_t,
          (Ptr{cudaSharedMemConfig},),
          pConfig)
end

function cudaDeviceSynchronize()
    ccall("extern cudaDeviceSynchronize", llvmcall, cudaError_t, ())
end

function cudaGetLastError()
    ccall("extern cudaGetLastError", llvmcall, cudaError_t, ())
end

function cudaPeekAtLastError()
    ccall("extern cudaPeekAtLastError", llvmcall, cudaError_t, ())
end

function cudaGetErrorString(error)
    ccall("extern cudaGetErrorString", llvmcall, Cstring,
          (cudaError_t,),
          error)
end

function cudaGetErrorName(error)
    ccall("extern cudaGetErrorName", llvmcall, Cstring,
          (cudaError_t,),
          error)
end

function cudaGetDeviceCount(count)
    ccall("extern cudaGetDeviceCount", llvmcall, cudaError_t,
          (Ptr{Cint},),
          count)
end

function cudaGetDevice(device)
    ccall("extern cudaGetDevice", llvmcall, cudaError_t,
          (Ptr{Cint},),
          device)
end

function cudaStreamCreateWithFlags(pStream, flags)
    ccall("extern cudaStreamCreateWithFlags", llvmcall, cudaError_t,
        (Ptr{cudaStream_t}, UInt32),
        pStream, flags)
end

function cudaStreamDestroy(stream)
    ccall("extern cudaStreamDestroy", llvmcall, cudaError_t,
        (cudaStream_t,),
        stream)
end

function cudaStreamWaitEvent(stream, event, flags)
    ccall("extern cudaStreamWaitEvent", llvmcall, cudaError_t,
          (cudaStream_t, cudaEvent_t, UInt32),
          stream, event, flags)
end

function cudaStreamWaitEvent_ptsz(stream, event, flags)
    ccall("extern cudaStreamWaitEvent_ptsz", llvmcall, cudaError_t,
          (cudaStream_t, cudaEvent_t, UInt32),
          stream, event, flags)
end

function cudaEventCreateWithFlags(event, flags)
    ccall("extern cudaEventCreateWithFlags", llvmcall, cudaError_t,
          (Ptr{cudaEvent_t}, UInt32),
          event, flags)
end

function cudaEventRecord(event, stream)
    ccall("extern cudaEventRecord", llvmcall, cudaError_t,
          (cudaEvent_t, cudaStream_t),
          event, stream)
end

function cudaEventRecord_ptsz(event, stream)
    ccall("extern cudaEventRecord_ptsz", llvmcall, cudaError_t,
          (cudaEvent_t, cudaStream_t),
          event, stream)
end

function cudaEventDestroy(event)
    ccall("extern cudaEventDestroy", llvmcall, cudaError_t,
          (cudaEvent_t,),
          event)
end

function cudaFuncGetAttributes(attr, func)
    ccall("extern cudaFuncGetAttributes", llvmcall, cudaError_t,
          (Ptr{cudaFuncAttributes}, Ptr{Cvoid}),
          attr, func)
end

function cudaFree(devPtr)
    ccall("extern cudaFree", llvmcall, cudaError_t,
          (Ptr{Cvoid},),
          devPtr)
end

function cudaMalloc(devPtr, size)
    ccall("extern cudaMalloc", llvmcall, cudaError_t,
          (Ptr{Ptr{Cvoid}}, Csize_t),
          devPtr, size)
end

function cudaMemcpyAsync(dst, src, count, kind, stream)
    ccall("extern cudaMemcpyAsync", llvmcall, cudaError_t,
          (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t, cudaMemcpyKind, cudaStream_t),
          dst, src, count, kind, stream)
end

function cudaMemcpyAsync_ptsz(dst, src, count, kind, stream)
    ccall("extern cudaMemcpyAsync_ptsz", llvmcall, cudaError_t,
          (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t, cudaMemcpyKind, cudaStream_t),
          dst, src, count, kind, stream)
end

function cudaMemcpy2DAsync(dst, dpitch, src, spitch, width, height, kind, stream)
    ccall("extern cudaMemcpy2DAsync", llvmcall, cudaError_t,
          (Ptr{Cvoid}, Csize_t, Ptr{Cvoid}, Csize_t, Csize_t, Csize_t, cudaMemcpyKind,
           cudaStream_t),
          dst, dpitch, src, spitch, width, height, kind, stream)
end

function cudaMemcpy2DAsync_ptsz(dst, dpitch, src, spitch, width, height, kind, stream)
    ccall("extern cudaMemcpy2DAsync_ptsz", llvmcall, cudaError_t,
          (Ptr{Cvoid}, Csize_t, Ptr{Cvoid}, Csize_t, Csize_t, Csize_t, cudaMemcpyKind,
           cudaStream_t),
          dst, dpitch, src, spitch, width, height, kind, stream)
end

function cudaMemcpy3DAsync(p, stream)
    ccall("extern cudaMemcpy3DAsync", llvmcall, cudaError_t,
          (Ptr{cudaMemcpy3DParms}, cudaStream_t),
          p, stream)
end

function cudaMemcpy3DAsync_ptsz(p, stream)
    ccall("extern cudaMemcpy3DAsync_ptsz", llvmcall, cudaError_t,
          (Ptr{cudaMemcpy3DParms}, cudaStream_t),
          p, stream)
end

function cudaMemsetAsync(devPtr, value, count, stream)
    ccall("extern cudaMemsetAsync", llvmcall, cudaError_t,
          (Ptr{Cvoid}, Cint, Csize_t, cudaStream_t),
          devPtr, value, count, stream)
end

function cudaMemsetAsync_ptsz(devPtr, value, count, stream)
    ccall("extern cudaMemsetAsync_ptsz", llvmcall, cudaError_t,
          (Ptr{Cvoid}, Cint, Csize_t, cudaStream_t),
          devPtr, value, count, stream)
end

function cudaMemset2DAsync(devPtr, pitch, value, width, height, stream)
    ccall("extern cudaMemset2DAsync", llvmcall, cudaError_t,
          (Ptr{Cvoid}, Csize_t, Cint, Csize_t, Csize_t, cudaStream_t),
          devPtr, pitch, value, width, height, stream)
end

function cudaMemset2DAsync_ptsz(devPtr, pitch, value, width, height, stream)
    ccall("extern cudaMemset2DAsync_ptsz", llvmcall, cudaError_t,
          (Ptr{Cvoid}, Csize_t, Cint, Csize_t, Csize_t, cudaStream_t),
          devPtr, pitch, value, width, height, stream)
end

function cudaMemset3DAsync(pitchedDevPtr, value, extent, stream)
    ccall("extern cudaMemset3DAsync", llvmcall, cudaError_t,
          (cudaPitchedPtr, Cint, cudaExtent, cudaStream_t),
          pitchedDevPtr, value, extent, stream)
end

function cudaMemset3DAsync_ptsz(pitchedDevPtr, value, extent, stream)
    ccall("extern cudaMemset3DAsync_ptsz", llvmcall, cudaError_t,
          (cudaPitchedPtr, Cint, cudaExtent, cudaStream_t),
          pitchedDevPtr, value, extent, stream)
end

function cudaRuntimeGetVersion(runtimeVersion)
    ccall("extern cudaRuntimeGetVersion", llvmcall, cudaError_t,
          (Ptr{Cint},),
          runtimeVersion)
end

function cudaGetParameterBuffer(alignment, size)
    ccall("extern cudaGetParameterBuffer", llvmcall, Ptr{Cvoid},
          (Csize_t, Csize_t),
          alignment, size)
end

function cudaGetParameterBufferV2(func, gridDimension, blockDimension, sharedMemSize)
    ccall("extern cudaGetParameterBufferV2", llvmcall, Ptr{Cvoid},
        (Ptr{Cvoid}, dim3, dim3, UInt32),
        func, gridDimension, blockDimension, sharedMemSize)
end

function cudaLaunchDevice_ptsz(func, parameterBuffer, gridDimension, blockDimension,
                               sharedMemSize, stream)
    ccall("extern cudaLaunchDevice_ptsz", llvmcall, cudaError_t,
          (Ptr{Cvoid}, Ptr{Cvoid}, dim3, dim3, UInt32, cudaStream_t),
          func, parameterBuffer, gridDimension, blockDimension, sharedMemSize, stream)
end

function cudaLaunchDeviceV2_ptsz(parameterBuffer, stream)
    ccall("extern cudaLaunchDeviceV2_ptsz", llvmcall, cudaError_t,
          (Ptr{Cvoid}, cudaStream_t),
          parameterBuffer, stream)
end

function cudaLaunchDevice(func, parameterBuffer, gridDimension, blockDimension,
                          sharedMemSize, stream)
    ccall("extern cudaLaunchDevice", llvmcall, cudaError_t,
          (Ptr{Cvoid}, Ptr{Cvoid}, dim3, dim3, UInt32, cudaStream_t),
          func, parameterBuffer, gridDimension, blockDimension, sharedMemSize, stream)
end

function cudaLaunchDeviceV2(parameterBuffer, stream)
    ccall("extern cudaLaunchDeviceV2", llvmcall, cudaError_t,
        (Ptr{Cvoid}, cudaStream_t),
        parameterBuffer, stream)
end

function cudaOccupancyMaxActiveBlocksPerMultiprocessor(numBlocks, func, blockSize,
                                                       dynamicSmemSize)
    ccall("extern cudaOccupancyMaxActiveBlocksPerMultiprocessor", llvmcall, cudaError_t,
          (Ptr{Cint}, Ptr{Cvoid}, Cint, Csize_t),
          numBlocks, func, blockSize, dynamicSmemSize)
end

function cudaOccupancyMaxActiveBlocksPerMultiprocessorWithFlags(numBlocks, func, blockSize,
                                                                dynamicSmemSize, flags)
    ccall("extern cudaOccupancyMaxActiveBlocksPerMultiprocessorWithFlags", llvmcall, cudaError_t,
          (Ptr{Cint}, Ptr{Cvoid}, Cint, Csize_t, UInt32),
          numBlocks, func, blockSize, dynamicSmemSize, flags)
end

function cudaCGGetIntrinsicHandle(scope)
    ccall("extern cudaCGGetIntrinsicHandle", llvmcall, Culonglong,
        (cudaCGScope,),
        scope)
end

function cudaCGSynchronize(handle, flags)
    ccall("extern cudaCGSynchronize", llvmcall, cudaError_t,
        (Culonglong, UInt32),
        handle, flags)
end

function cudaCGSynchronizeGrid(handle, flags)
    ccall("extern cudaCGSynchronizeGrid", llvmcall, cudaError_t,
          (Culonglong, UInt32),
          handle, flags)
end

function cudaCGGetSize(numThreads, numGrids, handle)
    ccall("extern cudaCGGetSize", llvmcall, cudaError_t,
          (Ptr{UInt32}, Ptr{UInt32}, Culonglong),
          numThreads, numGrids, handle)
end

function cudaCGGetRank(threadRank, gridRank, handle)
    ccall("extern cudaCGGetRank", llvmcall, cudaError_t,
          (Ptr{UInt32}, Ptr{UInt32}, Culonglong),
          threadRank, gridRank, handle)
end


text:
# LOOPS as in python

animals = ["dog", "cat", "mouse"]

for animal in animals
    println("$animal is a mammal")
end

# with enumeration
for (i,animal) in enumerate(animals)
    println("$animal no. $i is a mammal")
end

verbs = ["plays","sleeps","eats"]
for (animal,does) in zip(animals,verbs)
    println("The $animal $does.")
end

max_stars_repo_name: milankl/misc

max_issues_repo_name: milankl/misc

max_forks_repo_name: milankl/misc


text:
export Catalog

"""
    Catalog( data, ssm)

Data type to store analogs and succesors observations
from the Space State model

- [Example Catalog](@ref)
"""
struct Catalog

    nt::Int64
    nv::Int64
    data::Array{Float64,2}
    analogs::Array{Float64,2}
    successors::Array{Float64,2}
    sources::StateSpaceModel

    function Catalog(data::Array{Float64,2}, ssm::StateSpaceModel)

        nv, nt = size(data)
        analogs = data[:, 1:end-ssm.dt_states]
        successors = data[:, ssm.dt_states+1:end]

        new(nt, nv, data, analogs, successors, ssm)

    end

end


text:
module Test

using sbp.Config: Tv
using Test

function taylor_poly(x::AbstractArray, p::Int64)
        y = x.^p / factorial(p)
        return y
end

function test_first_derivative(D::AbstractArray, xl::AbstractArray, xr::AbstractArray, p_int::Int64, p_bnd::Int64, ml::Int64, mr::Int64)

        @test isapprox(interior(D*xr.^0, ml, mr), interior(eps(Tv) .* xl, ml, mr), atol=1.0)
        for p in 1:p_int
                yr = taylor_poly(xr, p)
                yl = taylor_poly(xl, p - 1)
                @test isapprox(interior(D*yr, ml, mr), interior(yl, ml, mr))
        end

        @test isapprox(boundary(D*xr.^0, ml, mr), boundary(eps(Tv) .* xl, ml, mr), atol=1.0)
        for p in 1:p_bnd
                yr = taylor_poly(xr, p)
                yl = taylor_poly(xl, p - 1)
                @test isapprox(boundary(D*yr, ml, mr), boundary(yl, ml, mr))
        end

end

function test_interpolation(P::AbstractArray, xl::AbstractArray, xr::AbstractArray, p_int::Int64, p_bnd::Int64, ml::Int64, mr::Int64)

        for p in 0:p_int
                yr = taylor_poly(xr, p)
                yl = taylor_poly(xl, p)
                @test isapprox(interior(P*yr, ml, mr), interior(yl, ml, mr))
        end

        for p in 0:p_bnd
                yr = taylor_poly(xr, p)
                yl = taylor_poly(xl, p)
                @test isapprox(boundary(P*yr, ml, mr), boundary(yl, ml, mr))
        end

end

function test_quadrature(H::AbstractArray, x::AbstractArray, pv::Int64)

        for p in 0:pv
                xr = taylor_poly(x, p)
                xl = taylor_poly(x, p + 1)
                @test isapprox( sum(H*xr), xl[end] - xl[1])
        end
end

function test_first_derivative_sbp(Hp::AbstractArray,
                                   Dp::AbstractArray,
                                   Hm::AbstractArray,
                                   Dm::AbstractArray,
                                   B::AbstractArray)
        @test isapprox(Hp*Dp + (Hm*Dm)', B)
end

function test_first_derivative_sbp(H::AbstractArray,
                                   D::AbstractArray,
                                   B::AbstractArray)
        @test isapprox(H*D + (H*D)', B)
end

function test_interpolation_sbp(Hp::AbstractArray,
                                Dp::AbstractArray,
                                Hm::AbstractArray,
                                Dm::AbstractArray)
        @test isapprox(Hp*Dp,  (Hm*Dm)')
end

function interior(x::AbstractArray, ml::Int64, mr::Int64)
        return x[ml+1:size(x,1) - mr]
end

function boundary(x::AbstractArray, ml::Int64, mr::Int64)
        y = zeros(ml + mr)
        for i in 1:ml
                @inbounds y[i] = x[i]
        end
        m = size(x, 1)
        for i in 1:mr
                @inbounds y[i] = x[m - mr + i]
        end
        return y
end

end

max_stars_repo_name: ooreilly/sbpjl

max_issues_repo_name: ooreilly/sbpjl

max_forks_repo_name: ooreilly/sbpjl


text:
@testset "Maximum weight matching" begin
    g = complete_graph(3)
    w = [
        1 2 1
        1 1 1
        3 1 1
    ]
    match = maximum_weight_matching(g, w; optimizer=GLPK.Optimizer)
    @test match.mate[1] == 3
    @test match.weight ≈ 3

    g = complete_graph(3)
    w = zeros(3, 3)
    w[1, 2] = 1
    w[3, 2] = 1
    w[1, 3] = 1
    match = maximum_weight_matching(g, w; optimizer=GLPK.Optimizer)
    @test match.weight ≈ 1

    g = Graph(4)
    add_edge!(g, 1, 3)
    add_edge!(g, 1, 4)
    add_edge!(g, 2, 4)

    w = zeros(4, 4)
    w[1, 3] = 1
    w[1, 4] = 3
    w[2, 4] = 1

    match = maximum_weight_matching(g, w; optimizer=GLPK.Optimizer)
    @test match.weight ≈ 3
    @test match.mate[1] == 4
    @test match.mate[2] == -1
    @test match.mate[3] == -1
    @test match.mate[4] == 1

    g = Graph(4)
    add_edge!(g, 1, 2)
    add_edge!(g, 2, 3)
    add_edge!(g, 3, 1)
    add_edge!(g, 3, 4)
    match = maximum_weight_matching(g; optimizer=GLPK.Optimizer)
    @test match.weight ≈ 2
    @test match.mate[1] == 2
    @test match.mate[2] == 1
    @test match.mate[3] == 4
    @test match.mate[4] == 3

    w = zeros(4, 4)
    w[1, 2] = 1
    w[2, 3] = 1
    w[1, 3] = 1
    w[3, 4] = 1

    match = maximum_weight_matching(g, w; optimizer=GLPK.Optimizer)
    @test match.weight ≈ 2
    @test match.mate[1] == 2
    @test match.mate[2] == 1
    @test match.mate[3] == 4
    @test match.mate[4] == 3

    w = zeros(4, 4)
    w[1, 2] = 1
    w[2, 3] = 1
    w[1, 3] = 5
    w[3, 4] = 1

    match = maximum_weight_matching(g, w; optimizer=GLPK.Optimizer)
    @test match.weight ≈ 5
    @test match.mate[1] == 3
    @test match.mate[2] == -1
    @test match.mate[3] == 1
    @test match.mate[4] == -1
end

@testset "Maximum weight maximal matching" begin
    g = complete_bipartite_graph(2, 2)
    w = zeros(4, 4)
    w[1, 3] = 10.0
    w[1, 4] = 1.0
    w[2, 3] = 2.0
    w[2, 4] = 11.0
    match = maximum_weight_maximal_matching(g, w; optimizer=GLPK.Optimizer)
    @test match.weight ≈ 21
    @test match.mate[1] == 3
    @test match.mate[3] == 1
    @test match.mate[2] == 4
    @test match.mate[4] == 2

    g = complete_bipartite_graph(2, 4)
    w = zeros(6, 6)
    w[1, 3] = 10
    w[1, 4] = 0.5
    w[2, 3] = 11
    w[2, 4] = 1
    match = maximum_weight_maximal_matching(g, w; optimizer=GLPK.Optimizer)
    @test match.weight ≈ 11.5
    @test match.mate[1] == 4
    @test match.mate[4] == 1
    @test match.mate[2] == 3
    @test match.mate[3] == 2

    g = complete_bipartite_graph(2, 6)
    w = zeros(8, 8)
    w[1, 3] = 10
    w[1, 4] = 0.5
    w[2, 3] = 11
    w[2, 4] = 1
    w[2, 5] = -1
    w[2, 6] = -1
    match = maximum_weight_maximal_matching(g, w, 0; optimizer=GLPK.Optimizer)
    @test match.weight ≈ 11.5
    @test match.mate[1] == 4
    @test match.mate[4] == 1
    @test match.mate[2] == 3
    @test match.mate[3] == 2

    g = complete_bipartite_graph(4, 2)
    w = zeros(6, 6)
    w[3, 5] = 10
    w[3, 6] = 0.5
    w[2, 5] = 11
    w[1, 6] = 1
    w[1, 5] = -1

    match = maximum_weight_maximal_matching(g, w, 0; optimizer=GLPK.Optimizer)
    @test match.weight ≈ 12
    @test match.mate[1] == 6
    @test match.mate[2] == 5
    @test match.mate[3] == -1
    @test match.mate[4] == -1
    @test match.mate[5] == 2
    @test match.mate[6] == 1
end