text:
abstract type Tri end
abstract type Tet end

struct eMesh{T1<:Union{Nothing,Tri},T2<:Union{Nothing,Tet}}
    point::Vector{SVector{3,Float64}}
    tri::Union{Nothing,Vector{SVector{3,Int64}}}
    tet::Union{Nothing,Vector{SVector{4,Int64}}}
    ϵ::Union{Nothing,Vector{Float64}}
    function eMesh( point::Vector{SVector{3,Float64}},
                    tri::Union{Nothing,Vector{SVector{3,Int64}}},
                    tet::Union{Nothing,Vector{SVector{4,Int64}}}=nothing,
                    ϵ::Union{Nothing,Vector{Float64}}=nothing)

        T1_ = ifelse(tri == nothing, Nothing, Tri)
        T2_ = ifelse(tet == nothing, Nothing, Tet)
        if T2_ == Tet
            @assert(isa(ϵ, Vector{Float64}))
            @assert(length(ϵ) == length(point), "length(ϵ) = $(length(ϵ)) but length(point) = $(length(point))")
            (length(ϵ) != 0) && @assert(0.0 < maximum(ϵ), "normalized penetration extent must be non-negative")
            (length(ϵ) != 0) && @assert(minimum(ϵ) == 0.0, "normalized penetration extent must be zero on the surface of the volume mesh")
            for k = 1:length(tet)
                (0.0 < volume(point[tet[k]])) || error("inverted tetrahedron")
            end
        else
            @assert(ϵ == nothing)
        end
        (T1_ == T2_ == Nothing) && error("a whole lot of nothing")
        return new{T1_,T2_}(point, tri, tet, ϵ)
    end
    function eMesh(hm::HomogenousMesh, tet::Union{Nothing,Vector{SVector{4,Int64}}}=nothing,
            ϵ::Union{Nothing,Vector{Float64}}=nothing)

        point = [SVector{3,Float64}(k) for k = hm.vertices]
        tri = [SVector{3,Int64}(k) for k = hm.faces]
        return eMesh(point, tri, tet, ϵ)
    end
    eMesh{Tri,Nothing}() = eMesh(Vector{SVector{3,Float64}}(), Vector{SVector{3,Int64}}(), nothing, nothing)
    eMesh{Nothing,Tet}() = eMesh(Vector{SVector{3,Float64}}(), nothing, Vector{SVector{4,Int64}}(), Vector{Float64}())
    eMesh{Tri,Tet}()     = eMesh(Vector{SVector{3,Float64}}(), Vector{SVector{3,Int64}}(), Vector{SVector{4,Int64}}(), Vector{Float64}())
end

as_tet_eMesh(e_mesh::eMesh{Tri,Tet}) = eMesh(e_mesh.point, nothing, e_mesh.tet, e_mesh.ϵ)
as_tet_eMesh(e_mesh::eMesh{Nothing,Tet}) = deepcopy(e_mesh)
as_tri_eMesh(e_mesh::eMesh{Tri,Tet}) = eMesh(e_mesh.point, e_mesh.tri, nothing, nothing)
as_tri_eMesh(e_mesh::eMesh{Tri,Nothing}) = deepcopy(e_mesh)

vertex_pos_for_tri_ind(eM::eMesh{Tri,T2}, k::Int64) where {T2} = eM.point[eM.tri[k]]
vertex_pos_for_tet_ind(eM::eMesh{T1,Tet}, k::Int64) where {T1} = eM.point[eM.tet[k]]

n_point(eM::eMesh) = length(eM.point)
n_tri(eM::eMesh) = length(eM.tri)
n_tet(eM::eMesh) = length(eM.tet)
get_tri(eM::eMesh) = eM.tri
get_tet(eM::eMesh) = eM.tet
get_point(eM::eMesh) = eM.point

function Base.empty!(e_mesh::eMesh{T1,T2}) where {T1,T2}
    empty!(e_mesh.point)
    (T1 == Nothing) || empty!(e_mesh.tri)
    (T2 == Nothing) || empty!(e_mesh.tet)
    (T2 == Nothing) || empty!(e_mesh.ϵ)
    return nothing
end

function Base.isempty(e_mesh::eMesh{T1,T2}) where {T1,T2}
    is_emp = isempty(e_mesh.point)
    if T1 == Tri
        is_emp = is_emp || isempty(e_mesh.tri)
    end
    if T2 == Tet
        is_emp = is_emp || isempty(e_mesh.tet)
        is_emp = is_emp || isempty(e_mesh.ϵ)
    end
    return is_emp
end

function Base.append!(eM_1::eMesh{T1,T2}, eM_2::eMesh{T1,T2}) where {T1,T2}
    n_1 = n_point(eM_1)
    if T1 != Nothing
        for k = eM_2.tri
            push!(eM_1.tri, k .+ n_1)
        end
    end
    if T2 != Nothing
        for k = eM_2.tet
            push!(eM_1.tet, k .+ n_1)
        end
        append!(eM_1.ϵ, eM_2.ϵ)
    end
    append!(eM_1.point, eM_2.point)
    return nothing
end

### VERIFICATION

function verify_mesh(eM::eMesh{T1,T2}) where {T1,T2}
    verify_mesh_triangles(eM)
    verify_mesh_tets(eM)
    return nothing
end

verify_mesh_tets(eM::eMesh{T1,Nothing}) where {T1} = nothing
function verify_mesh_tets(eM::eMesh{T1,Tet}) where {T1}
    length_tet = n_tet(eM)
    length_point = n_point(eM)
    length_ϵ = length(eM.ϵ)
    for k = 1:length_tet
        iΔ = eM.tet[k]
        all(1 .<= iΔ .<= length_point) || error("index in tet with sides $(iΔ) not within points")
    end
    (length_ϵ == length_point) || error("number of points not equal to number or ϵ")
    return nothing
end

verify_mesh_triangles(eM::eMesh{Nothing,T2}) where {T2} = nothing
function verify_mesh_triangles(eM::eMesh{Tri,T2}) where {T2}
    length_tri = n_tri(eM)
    length_point = n_point(eM)
    for k = 1:length_tri
        iΔ = eM.tri[k]
        all(1 .<= iΔ .<= length_point) || error("index in triangle with sides $(iΔ) not within points")
    end
    return nothing
end

### MESH MANIPULATION

function eMesh_transform!(e_mesh::eMesh{T1,T2}, extra_args...) where {T1,T2}
    dh = basic_dh(extra_args...)
    point = e_mesh.point
    for k = 1:n_point(e_mesh)
        point[k] = dh_vector_mul(dh, point[k])
    end
end

struct PointObj
    p::SVector{3,Float64}
    w::Float64
    b::Bool
    k::Int64
end

function crop_mesh(eM::eMesh{Tri,Nothing}, plane::SMatrix{1,4,Float64,4})
    function add_cropped_triangle!(eM_new::eMesh{Tri,Nothing}, iΔ::SVector{3,Int64})
        function plane_dot_2(eM_new::eMesh{Tri,Nothing}, plane::SMatrix{1,4,Float64,4}, k::Int64)
            p = get_point(eM_new)[k]
            w = dot(plane, onePad(p))
            b = -1.0e-12 < w
            return PointObj(p, w, b, k)
        end

        o1 = plane_dot_2(eM_new, plane, iΔ[1])
        o2 = plane_dot_2(eM_new, plane, iΔ[2])
        o3 = plane_dot_2(eM_new, plane, iΔ[3])
        s_bool = o1.b + o2.b + o3.b
        (s_bool == 3) && (push!(eM_new.tri, iΔ))  # all entire triangle
        ((s_bool == 0) || (s_bool == 3)) && (return nothing)  # exit function
        (o1.b == o2.b) && ((o1, o2, o3) = (o2, o3, o1))  # put 1 and 2 on opposite sides
        (o2.b == o3.b) && ((o1, o2, o3) = (o3, o1, o2))  # put 2 and 3 on opposite sides
        (o1.b == o2.b) || (o2.b == o3.b) && error("something is wrong")  # sanity check
        push!(eM_new.point, weightPoly(o1.p, o2.p, o1.w, o2.w))
        i12 = length(get_point(eM_new))
        push!(eM_new.point, weightPoly(o2.p, o3.p, o2.w, o3.w))
        i23 = length(get_point(eM_new))
        if o2.b  # add top of triangle
            push!(eM_new.tri, SVector{3,Int64}(i12, o2.k, i23))
        else  # add bottom quadralateral
            push!(eM_new.tri, SVector{3,Int64}(o1.k, i12, i23))
            push!(eM_new.tri, SVector{3,Int64}(o1.k, i23, o3.k))
        end
    end

    eM_new = eMesh{Tri,Nothing}()
    append!(eM_new.point, eM.point)
    for iΔ = eM.tri
        add_cropped_triangle!(eM_new, iΔ)
    end
    mesh_repair!(eM_new)
    return eM_new
end

function invert!(eM::eMesh{Tri,Nothing})
    for k = 1:n_tri(eM)
        tri_k = eM.tri[k]
        eM.tri[k] = SVector{3,Float64}(tri_k[3], tri_k[2], tri_k[1])
    end
end

### MESH REPAIR
function mesh_repair!(e_mesh::eMesh{T1,T2}) where {T1,T2}
    mesh_remove_unused_points!(e_mesh)
    mesh_inplace_rekey!(e_mesh)
    mesh_remove_unused_points!(e_mesh)
    return delete_triangles!(e_mesh)
end

function remove_degenerate!(eM::eMesh)
	area_or_vol(v::SVector{3,SVector{3,Float64}}) = area(v)
	area_or_vol(v::SVector{4,SVector{3,Float64}}) = volume(v)

	tol = 1.0e-6
	for vec_ind = (eM.tet, eM.tri)
		if vec_ind != nothing
			v = [area_or_vol(eM.point[vec_ind[k]]) for k = 1:length(vec_ind)]
			max_v = maximum(v)
			i_bad = findall(v .< max_v * tol)
			deleteat!(vec_ind, i_bad)
		end
	end
end

rekey!(v::Nothing, i::Vector{Int64}) = nothing
rekey!(v::Vector{SVector{N,Int64}}, i::Vector{Int64}) where {N} = replace!(x -> i[x], v)

function mesh_inplace_rekey(e_mesh::eMesh{T1,T2}) where {T1,T2}
    e_mesh_copy = deepcopy(e_mesh)
    mesh_inplace_rekey!(e_mesh_copy)
    return e_mesh_copy
end

function mesh_inplace_rekey!(e_mesh::eMesh{T1,T2}) where {T1,T2}
    # Expresses the mesh using the minumum number of points

    function shortest_side(e_mesh::eMesh)
        shortest_side(v::Nothing) = Inf
        shortest_side(ind::Vector{SVector{N,Int64}}) where {N} =  minimum([shortest_side(e_mesh.point[k]) for k = ind])
        function shortest_side(v::SVector{N,SVector{3,Float64}}) where {N}
            d_min = Inf
            for k = 1:N
                for kk = 1:(k-1)
                    d_min = min(d_min, norm(v[k] - v[kk]))
                end
            end
            return d_min
        end

        return min(shortest_side(e_mesh.tri), shortest_side(e_mesh.tet))
    end

    function inplace_rekey(point::Vector{SVector{3,Float64}}, min_side_length::Float64)
        balltree = BallTree(point; reorder = false)
        idxs = inrange(balltree, point, min_side_length * 0.499)  # Assume that all points closer to each other than hald a side length are duplicates
        return first.(sort!.(idxs))
    end

    if !isempty(e_mesh)  # reducing over an empty collection is not allowed
        min_side_length = shortest_side(e_mesh)
        new_key = inplace_rekey(e_mesh.point, min_side_length)
        rekey!(e_mesh.tri, new_key)
        rekey!(e_mesh.tet, new_key)
        mesh_remove_unused_points!(e_mesh)
        return new_key
    end
    return nothing
end

function mesh_remove_unused_points!(e_mesh::eMesh{T1,T2}) where {T1,T2}
    truth_vector = falses(n_point(e_mesh))
    (T1 != Nothing) && (for k = e_mesh.tri; truth_vector[k] = true; end)
    (T2 != Nothing) && (for k = e_mesh.tet; truth_vector[k] = true; end)
    new_key = cumsum(truth_vector)
    rekey!(e_mesh.tri, new_key)
    rekey!(e_mesh.tet, new_key)
    ind_truth = findall(truth_vector)
    point_new = e_mesh.point[ind_truth]
    empty!(e_mesh.point)
    append!(e_mesh.point, point_new)
    if T2 != Nothing
        ϵ_new = e_mesh.ϵ[ind_truth]
        empty!(e_mesh.ϵ)
        append!(e_mesh.ϵ, ϵ_new)
    end
    return nothing
end

delete_triangles!(e_mesh::eMesh{Nothing,T2}) where {T2} = -9999
function delete_triangles!(e_mesh::eMesh{Tri,T2}) where {T2}
    ### make the first index lowest
    for k = 1:n_tri(e_mesh)
        iΔ = e_mesh.tri[k]
        min_index = minimum(iΔ)
        (min_index == iΔ[1]) || (iΔ = SVector(iΔ[3], iΔ[1], iΔ[2]))
        (min_index == iΔ[1]) || (iΔ = SVector(iΔ[3], iΔ[1], iΔ[2]))
    end

    ### create a dictionary of key repition counts
    key_type = Vector{Tuple{Int64,SVector{3,Int64}}}
    dict_ind = Dict{SVector{3,Int64},key_type}()
    for k = 1:n_tri(e_mesh)
        iΔ = e_mesh.tri[k]
        sort_iΔ = sort(iΔ)
        !haskey(dict_ind, sort_iΔ) && (dict_ind[sort_iΔ] = key_type() )
        push!(dict_ind[sort_iΔ], (k, iΔ))
    end

    ### delete duplicates
    i_delete = Vector{Int64}()
    for (key_k, val_k) = dict_ind
        length_val_k = length(val_k)
        if length_val_k == 2
            val_k1 = val_k[1]
            val_k2 = val_k[2]
            (val_k1[2] == val_k2[2]) && error("non-opposing triangles")
            push!(i_delete, val_k1[1], val_k2[1])
        elseif 3 <= length_val_k
            error("something is wrong")
        end
    end
end

function sub_div_mesh(eM_ico::eMesh{Tri,T2}, n_div::Int64) where {T2}
    n_end(n::Int64) = div((n + 1) * n, 2)
    n_start(n::Int64) = 1 + n_end(n - 1)

    function sub_div_triangle(p::SVector{3,SVector{3,Float64}}, n_div::Int64)
        function sub_div_triangle_vert_index(n_div::Int64)
            i_tri = Vector{SVector{3,Int64}}()
            for k = 1:n_div
                for kk = 0:(k - 1)
                    i1 = n_start(k) + kk
                    i2 = i1 + k
                    i3 = i2 + 1
                    push!(i_tri, SVector{3,Int64}(i1, i2, i3))
                end
                for kk = 0:(k - 2)
                    i1 = n_start(k) + kk
                    i2 = i1 + k + 1
                    i3 = i2 - k
                    push!(i_tri, SVector{3,Int64}(i1, i2, i3))
                end
            end
            return i_tri
        end

        function get_new_point(n_vert::Int64, n_div::Int64, p::SVector{3,SVector{3,Float64}})
            function find_layer_1(i_point::Int64)
                i_end_layer = 1
                i_layer = 1
                while i_end_layer < i_point
                    i_layer += 1
                    i_end_layer += i_layer
                end
                norm_extent_1 = ifelse(i_point == 1, 0.0, (i_end_layer - i_point) / (i_layer - 1) )
                return i_layer, norm_extent_1
            end

            i_layer, norm_extent_1 = find_layer_1(n_vert)
            ϕ_1 = (n_div - i_layer + 1) / n_div
            ϕ_2 = (1 - ϕ_1) * norm_extent_1
            ϕ_3 = 1 - ϕ_1 - ϕ_2
            ϕ = SVector{3,Float64}(ϕ_1, ϕ_2, ϕ_3)
            return sum(p .* ϕ)
        end

        point = Vector{SVector{3,Float64}}()
        i_tri_div = sub_div_triangle_vert_index(n_div)
        for k = 1:n_end(n_div + 1)
            push!(point, get_new_point(k, n_div, p))
        end
        return eMesh(point, i_tri_div, nothing, nothing)
    end

    eM_ico_div = eMesh{Tri,Nothing}()
    for k = 1:n_tri(eM_ico)
        p = vertex_pos_for_tri_ind(eM_ico, k)
        append!(eM_ico_div, sub_div_triangle(p, n_div))
    end
    mesh_repair!(eM_ico_div)
    return eM_ico_div
end

### BASIC SHAPES
function output_eMesh_half_plane(plane_w::Float64=1.0, is_include_vis_sides::Bool=false)
    v1, v2, v3 = [SVector{3,Float64}(cos(theta), sin(theta), 0.0) for theta = (0.0, 2*pi/3, 4*pi/3)]
    v4 = SVector{3,Float64}(0.0, 0.0, -1.0) * plane_w
    point = [v1, v2, v3, v4]
    if is_include_vis_sides
        tri = [SVector{3,Int64}(1,2,3), SVector{3,Int64}(1,3,4),  SVector{3,Int64}(1,4,2), SVector{3,Int64}(2,4,3)]
    else
        tri = [SVector{3,Int64}(1,2,3)]
    end
    tet = [SVector{4,Int64}(4,1,2,3)]
    ϵ = [0.0, 0.0, 0.0, plane_w]
    return eMesh(point, tri, tet, ϵ)
end

function output_eMesh_sphere(rad::Union{Float64,SVector{3,Float64}}=1.0, n_div::Int64=4)
    function make_icosahedron()
        φ = Base.MathConstants.golden

        v = Vector{SVector{3,Float64}}()
        for s_1 = (-1.0, +1.0)
            for s_2 = (-1.0, +1.0)
                push!(v, SVector{3,Float64}(    0.0,     s_1, φ * s_2) )
                push!(v, SVector{3,Float64}(    s_1, φ * s_2,     0.0) )
                push!(v, SVector{3,Float64}(φ * s_2,     0.0,     s_1) )
            end
        end

        n_ico_vert = 12
        d = zeros(n_ico_vert, n_ico_vert)
        for k = 1:n_ico_vert
            for kk = 1:n_ico_vert
                d[k, kk] = norm(v[k] - v[kk])
            end
        end

        v_face_vert = Vector{SVector{3,Int64}}()
        b = d .== 2.0
        for i1 = 1:n_ico_vert
            for i2 = 1:n_ico_vert
                for i3 = 1:n_ico_vert
                    if (i1 < i2 < i3) && b[i1, i2] && b[i2, i3] && b[i1, i3]
                        p1 = v[i1]
                        p2 = v[i2]
                        p3 = v[i3]
                        n = triangleNormal(p1, p2, p3)
                        c = normalize(centroid(p1, p2, p3))
                        i_face = ifelse(n ≈ c, SVector{3,Int64}(i1, i2, i3), SVector{3,Int64}(i1, i3, i2))
                        push!(v_face_vert, i_face)
                    end
                end
            end
        end

        push!(v, SVector{3,Float64}(0,0,0))
        v_tet = Vector{SVector{4,Int64}}()
        for k = 1:length(v_face_vert)
            push!(v_tet, SVector{4,Int64}(13, v_face_vert[k]...))
        end
        ϵ = vcat(zeros(n_ico_vert), 1.0)
        return eMesh(v, v_face_vert, v_tet, ϵ)
    end

    function volumize_about(eM::eMesh{Tri,Nothing})
        i_tet = Vector{SVector{4,Int64}}()
        n_vert = n_point(eM)
        n_center = n_vert + 1
        for k = 1:n_tri(eM)
            push!(i_tet, SVector{4,Int64}(n_center, eM.tri[k]...))
        end
        ϵ = vcat(zeros(n_vert), 1.0)
        point = deepcopy(eM.point)
        push!(point, zeros(SVector{3,Float64}))
        return eMesh(point, eM.tri, i_tet, ϵ)
    end

    function project_to_sphere!(eM::eMesh)
        for k = 1:n_point(eM)
            eM.point[k] = normalize(eM.point[k])
        end
        return nothing
    end

    eM_ico = make_icosahedron()
    eM_ico_div = sub_div_mesh(eM_ico, n_div)
    project_to_sphere!(eM_ico_div)

    rad = ones(SVector{3,Float64}) .* rad
    eMesh_transform!(eM_ico_div, Diagonal(rad))

    return volumize_about(eM_ico_div)
end

function output_box_ind()
    oriented_box_faces = (
        SVector{4, Int64}(1,3,5,7),  # -x
        SVector{4, Int64}(2,6,4,8),  # +x
        SVector{4, Int64}(1,5,2,6),  # -y
        SVector{4, Int64}(3,4,7,8),  # +y
        SVector{4, Int64}(1,2,3,4),  # -z
        SVector{4, Int64}(5,7,6,8),  # +z
    )
    ϵ = zeros(Float64,8)
    push!(ϵ, 1.0)
    tri = Vector{SVector{3,Int64}}()
    tet = Vector{SVector{4,Int64}}()
    for k = 1:6
        bf_k = oriented_box_faces[k]
        push!(tri, bf_k[SVector{3,Int64}(1,3,4)])
        push!(tri, bf_k[SVector{3,Int64}(1,4,2)])
    end
    for k = 1:length(tri)
        tri_k = tri[k]
        push!(tet, SVector{4,Int64}(9, tri_k[1], tri_k[2], tri_k[3]))
    end
    return tri, tet, ϵ
end

function output_eMesh_box(r::Union{Float64,SVector{3,Float64}}=1.0, c::SVector{3,Float64}=zeros(SVector{3,Float64}))
    point = [
        SVector{3,Float64}(-1,-1,-1),
        SVector{3,Float64}(+1,-1,-1),
        SVector{3,Float64}(-1,+1,-1),
        SVector{3,Float64}(+1,+1,-1),
        SVector{3,Float64}(-1,-1,+1),
        SVector{3,Float64}(+1,-1,+1),
        SVector{3,Float64}(-1,+1,+1),
        SVector{3,Float64}(+1,+1,+1),
        SVector{3,Float64}( 0, 0, 0),
    ]
    tri, tet, ϵ = output_box_ind()
    e_mesh = eMesh(point, tri, tet, ϵ)
    r = r .* ones(SVector{3,Float64})
    eMesh_transform!(e_mesh, SMatrix{3,3,Float64,9}(r[1], 0, 0, 0, r[2], 0, 0, 0, r[3]))
    eMesh_transform!(e_mesh, c)
    return e_mesh
end

# function make_cone_points(rⁱ::Float64, rᵒ::Float64, h::Float64, k_slice::Int64, tot_slice::Int64)
#     polar_point(r::Float64, θ::Float64, z::Float64) = SVector{3,Float64}(r * cos(θ), r * sin(θ), z)
#
#     @assert(0 < rⁱ < rᵒ)
#     @assert(0 < h)
#     @assert(1 <= k_slice <= tot_slice)
#
#     v = Vector{SVector{3,Float64}}()
#     θ_space = LinRange{Float64}(0.0, 2*pi, tot_slice + 1)
#     θ_1 = θ_space[k_slice]
#     θ_2 = θ_space[k_slice + 1]
#     for k = 1:8
#         θ = ifelse(mod1(k, 4) <= 2, θ_1, θ_2)
#         r = ifelse(isodd(k), rⁱ, rᵒ)
#         z = ifelse(k <= 4, -0.5, +0.5) * h
#         push!(v, polar_point(r, θ, z))
#     end
#     push!(v, polar_point((rⁱ + rᵒ) / 2, (θ_1 + θ_2) / 2, 0.0))
# end
#
# function output_eMesh_slice(rⁱ::Float64, rᵒ::Float64, h::Float64, k_slice::Int64, tot_slice::Int64)
#     point = make_cone_points(rⁱ, rᵒ, h, k_slice, tot_slice)
#     tri, tet, ϵ = output_box_ind()
#     return eMesh(point, tri, tet, ϵ)
# end
#
# function output_eMesh_hole(rⁱ::Float64, rᵒ::Float64, h::Float64, tot_slice::Int64)
#     eM_cone = eMesh{Tri,Tet}()
#     for k = 1:tot_slice
#         append!(eM_cone, output_eMesh_slice(rⁱ, rᵒ, h, k, tot_slice))
#     end
#     return eM_cone
# end


text:
using Knet
# https://www.tensorflow.org/api_docs/python/tf/nn/sigmoid_cross_entropy_with_logits

function deeper_cut_combined_loss(scores, labels)
    part_detection_labels = labels[:, :, 1:global_num_joints, :]
    part_detection_scores = scores[:, :, 1:global_num_joints, :]
    part_detection_loss =
        sigmoid_cross_entropy_loss(part_detection_scores, part_detection_labels)

    score_channel_size = (scores|>size)[end-1]
    loc_ref_loss = 0
    if score_channel_size >= global_num_joints * 5
        loc_ref_labels = labels[:, :, global_num_joints+1:global_num_joints*3, :]
        loc_ref_scores = scores[:, :, global_num_joints+1:global_num_joints*3, :]

        loc_ref_labels_weights = labels[:, :, global_num_joints*3+1:global_num_joints*5, :]
        loc_ref_scores_weights = scores[:, :, global_num_joints*3+1:global_num_joints*5, :]

        loc_ref_weight_loss =
            sigmoid_cross_entropy_loss(loc_ref_scores_weights, loc_ref_labels_weights)

        loc_ref_loss =
            huber_loss(loc_ref_scores, loc_ref_labels; weights = loc_ref_scores_weights)

        loc_ref_loss = loc_ref_loss + loc_ref_weight_loss
    elseif score_channel_size == global_num_joints * 3
        loc_ref_labels = labels[:, :, global_num_joints+1:global_num_joints*3, :]
        loc_ref_scores = scores[:, :, global_num_joints+1:global_num_joints*3, :]
        loc_ref_loss = huber_loss(loc_ref_scores, loc_ref_labels)
    end

    return part_detection_loss + locref_loss_weight * loc_ref_loss
end

function sigmoid_cross_entropy_loss(scores, labels)
    # max.(y_hat, 0) .- y_hat .* y + log(1 .+ exp(-abs.(y_hat)))
    #   mask = Knet.atype()(zeros(Float32, size(labels)))
    #  mask[:, :, 1:14 ,:] .= 1
    # scores = scores .* mask
    # labels .* mask
    l = max.(0, scores) .- labels .* scores .+ log.(1 .+ exp.(-abs.(scores)))
    (sum(l) / length(l))
end

function huber_loss(scores, labels; weights = 1, beta = 1)
    diff = abs.(labels .- scores)
    cpu_diff = Array(diff)
    low_idx = findall(cpu_diff .< beta)
    high_idx = findall(cpu_diff .>= beta)

    # TODO: put weights into sigmoids
    sigmoided_weights = weights == 1 ? Knet.atype()(ones(size(scores))) : sigm.(weights)

    loss_sum = sum((diff[high_idx] .- 0.5 .* beta^2) .* sigmoided_weights[high_idx])
    loss_sum += sum((0.5 .* (diff[low_idx] .^ 2) .* sigmoided_weights[low_idx]))
    return sum(loss_sum) / length(diff)
end

#=
def huber_loss(labels, predictions, weight=1.0, k=1.0, scope=None):
    """Define a huber loss  https://en.wikipedia.org/wiki/Huber_loss
      tensor: tensor to regularize.
      k: value of k in the huber loss
      scope: Optional scope for op_scope.

    Huber loss:
    f(x) = if |x| <= k:
              0.5 * x^2
           else:
              k * |x| - 0.5 * k^2

    Returns:
      the L1 loss op.

    http://concise-bio.readthedocs.io/en/latest/_modules/concise/tf_helper.html
    """
    with ops.name_scope(scope, "absolute_difference",
                        [predictions, labels]) as scope:
        predictions.get_shape().assert_is_compatible_with(labels.get_shape())
        if weight is None:
            raise ValueError("`weight` cannot be None")
        predictions = math_ops.to_float(predictions)
        labels = math_ops.to_float(labels)
        diff = math_ops.subtract(predictions, labels)
        abs_diff = tf.abs(diff)
        losses = tf.where(abs_diff < k,
                          0.5 * tf.square(diff),
                          k * abs_diff - 0.5 * k ** 2)
        return tf.losses.compute_weighted_loss(losses, weight)

=#

ext: jl


text:
## Imported math functions

## need to import these
math_sympy_methods_base = (:sin, :cos, :tan, :sinh, :cosh, :tanh, :asin, :acos, :atan,
                           :asinh, :acosh, :atanh, :sec, :csc, :cot, :asec, :acsc, :acot,
                           :sech, :csch,
                           :coth, :acoth,
                           :log2, :log10, :log1p, :exponent, :exp, :exp2, :expm1,
                           :sqrt, :dawson,
                           :ceiling, :floor,
                           :factorial,
                           :gcd, :lcm,
                           :isqrt
)

math_object_methods_base = (:round,
                            )

## Export SymPy math functions and vectorize them
math_sympy_methods = (:radians2degrees, :degrees2radians,
                      :factorial2,
                      )

# hypot and atan2
hypot(x::Sym, y::Number) = sqrt(x^2 + y^2)
atan2(y::Sym, x::Number) = sympy_meth(:atan2, y, x)

## Log function handles arguments differently
log(x::Sym) = sympy_meth(:log, x)
log(b::Sym, x::Sym) = sympy_meth(:log, x, b)

### Trigonometry
Base.rad2deg(x::Sym) = radians2degrees(x)
Base.deg2rad(x::Sym) = degrees2radians(x)

## degree functions
for fn in (:cosd, :cotd, :cscd, :secd, :sind, :tand,
          :acosd, :acotd, :acscd, :asecd, :asind, :atand)

    rad_fn = string(fn)[1:end-1]
    @eval ($fn)(x::Sym) = sympy[Symbol($rad_fn)](x * Sym(sympy["pi"])/180)
end

for fn in (:cospi, :sinpi)
    rad_fn = string(fn)[1:end-2]
    @eval ($fn)(x::Sym) = sympy[Symbol($rad_fn)](x * Sym(sympy["pi"]))
end

## :asech, :acsch, :sinc, :cosc,
## These fail, so define from definitions
## http://mathworld.wolfram.com/InverseHyperbolicSecant.html
asech(z::Sym) = log(sqrt(1/z-1)*sqrt(1/z+1) + 1/z)
## http://mathworld.wolfram.com/InverseHyperbolicCosecant.html
acsch(z::Sym) = log(sqrt(1+1/z^2) + 1/z)
# Julia's sinc is defined to be zero at x=1
sinc(x::Sym) = piecewise((Sym(1), Eq(x, 0)), (sin(PI*x)/(PI*x), Gt(abs(x), 0)))
cosc(x::Sym) = diff(sinc(x))

# deprecate these when v0.4 support dropped in favor of `asech.(...)` form
asech(as::AbstractArray{Sym}) = map(asech, as)
acsch(as::AbstractArray{Sym}) = map(acsch, as)

## in Julia, not SymPy
cbrt(x::Sym) = x^(1//3)
Base.ceil(x::Sym) = ceiling(x)

functions_sympy_methods = (
                           :arg,
                           :conjugate,
                           :re,
                           :sign
                           )

## map Abs->abs, Max->max, Min->min
abs(ex::Sym, args...; kwargs...) = sympy_meth(:Abs, ex, args...; kwargs...)
Base.abs2(x::Sym) = re(x*conj(x))
Base.copysign(x::Sym, y::Sym) = abs(x)*sign(y)
Base.signbit(x::Sym) = x < 0
Base.flipsign(x::Sym, y) = signbit(y) ? -x : x
Base.real(x::Sym) = sympy_meth(:re, x)
Base.imag(x::Sym) = sympy_meth(:im, x)
Base.eps(::Type{Sym}) = zero(Sym)

## use SymPy Names here...
## XXX This pattern may be of general usage. XXX
## It allows fn(Any, Sym) or f(Sym, Any) to match

import Base: min, max
min(x::Sym, a) = sympy_meth(:Min, x, a)
min(a, x::Union{SA, Real}) where {SA <: Sym} = min(x,a)

max(x::Sym, a) = sympy_meth(:Max, x, a)
max(a, x::Union{SA, Real}) where {SA <: Sym} = max(x,a)

# SymPy names
Min(ex::Sym, ex1::Sym) = sympy_meth(:Min, ex, ex1)
Max(ex::Sym, ex1::Sym) = sympy_meth(:Max, ex, ex1)

## Calculus functions

## bring over for univariate function calls as well (not just expressions)

"""

Compute a symbolic limit lim_{x->c+} f(x)

By default, *right* limits are returned. The keyword argument `dir="-"` needs to be
specified for *left* limits.

The function `f(x)` can be expressed as a symbolic expression, or a univariate function

The values `x` and `c` can be expressed by two arguments `x, c` or a pair `x=>c`.

```
@vars x y
ex = 1/(x^2 - y^2)
limit(ex, x=>3)
limit(ex, x, 3)
limit(ex, x=>3, y=>2)  # first x then y not (x,y) -> (3,2)

fn(x) = x^x
limit(fn, 0)    # symbol not needed
```
"""
limit(ex::Sym, args...; kwargs...) = sympy_meth(:limit, ex, args...; kwargs...)
limit(ex::Sym, d::Pair; kwargs...) = limit(ex, d.first, d.second;kwargs...)
limit(ex::Sym, ds::Pair...) = reduce(limit, ex, ds)
limit(f::Function, c::Number=0; kwargs...) = (z = (symbols(gensym())); limit(f(z),z=>c;kwargs...))

export limit

function Base.diff(ex::Sym, args...; kwargs...)
    if funcname(ex) in map(string, relational_sympy_values)
        Eq(diff(lhs(ex), args...; kwargs...), diff(rhs(ex), args...; kwargs...))
    else
        sympy_meth(:diff, ex, args...; kwargs...)
    end
end

##

## diff for matrix doesn't handle vectors well, so we vectorize here
diff(exs::AbstractVector{Sym}, args...; kwargs...) = map(ex -> diff(ex, args...;kwargs...), exs)
diff(exs::AbstractMatrix{Sym}, args...; kwargs...) = map(ex -> diff(ex, args...;kwargs...), exs)

## find symbolic derivatives from a function
function diff(f::Function, k::Int=1; kwargs...)
    x = Sym("x")
    diff(f(x), x, k; kwargs...)
end

# set up derivative, call doit to implement
Derivative(ex::Sym, args...) = sympy_meth(:Derivative, ex, args...)
export(Derivative)
#### Piecewise functions

"""

Create a piecewise defined function.

To create conditions on the variable, the functions `Lt`, `Le`, `Eq`, `Ge`, and `Gt` can be used. For infix notation,
unicode operators can be used: `\\ll<tab>`, `\\leqq<tab>`, `\\Equal<tab>`, `\\geqq<tab>`, and `\\gg<tab>` (but *not* `\\ge<tab>` or `\\le<tab>`).

To combine terms, the unicode `\\vee<tab>` (for "or"), `\\wedge<tab>` (for "and") can be used

Examples:
```
x,a = symbols("x,a")
p = piecewise((1, x ≪ 1), (2, (Lt(1,x)) ∨ Lt(x,2)), (3, x ≫ 2)) ## using ∨ and ∧ for & and or
subs(p, x, 2) ## 2
p = piecewise((1, Lt(x, a)), (2, Ge(x,a)))  # same as piecewise((1,  x ≪ a), (2, x ≧ a))
subs(p, x, a - 1)
```

[Note: there is also an alias `Piecewise` for copy-n-pasting from python code, but despite the capital letter, this is not a constructor for a type.]
"""
function piecewise(args...)
    sympy_meth(:Piecewise, args...)
end
const Piecewise = piecewise

piecewise_fold(ex::Sym) = sympy_meth(:piecewise_fold, ex)

Base.ifelse(ex::Sym, a, b) = piecewise((a, ex), (b, true))

"""
Indicator expression: (Either `\\Chi[tab](x,a,b)` or `Indicator(x,a,b)`)

`Χ(x, a, b)` is `1` on `[a,b]` and 0 otherwise.

This is not a function taking `x`, but a symbolic expression of `x`.

"""
Χ(x::Sym, a=-oo, b=oo) = piecewise((1, Gt(x, a) ∧ Le(x, b)), (0,true))
Indicator(x::Sym, a=-oo, b=oo) = Χ(x, a, b)

import Base: &, |
(&)(a::Bool, b::Sym) = a & (b == SympyTRUE)
(|)(a::Bool, b::Sym) = a | (b == SympyTRUE)

export Indicator, Χ

### generic programming interface

## Handle ininf, and isnan by coercion to float
Base.isfinite(x::Sym) = isfinite(convert(Float64, x))
Base.isinf(x::Sym) = try isinf(convert(Float64, x)) catch e false end
Base.isnan(x::Sym) = try isnan(convert(Float64, x)) catch e false end

## we rename sympy.div -> polydiv
Base.div(x::Sym, y::SymOrNumber) = convert(Sym, sympy["floor"](x/convert(Sym,y)))
Base.rem(x::Sym, y::SymOrNumber) = x-Sym(y)*Sym(sympy["floor"](x/y))

## zero and one (zeros?)
Base.zero(x::Sym) = Sym(0)
Base.zero(::Type{Sym}) = Sym(0)

Base.one(x::Sym) = Sym(1)
Base.one(::Type{Sym}) = Sym(1)

## useful at times
Base.typemax(::Type{Sym}) = oo
Base.typemin(::Type{Sym}) = -oo

##################################################
## special numbers are initialized after compilation
if isdefined(PyCall,:PyNULL)
    pynull() = PyCall.PyNULL()
else
    pynull() = PyCall.PyObject()
end
global PI = Sym(pynull())
global E = Sym(pynull())
global IM = Sym(pynull())
global oo = Sym(pynull())

if isdefined(Base, :MathConstants)
    Base.convert(::Type{Sym}, x::Irrational{:π}) = PI
    Base.convert(::Type{Sym}, x::Irrational{:e}) = E
    Base.convert(::Type{Sym}, x::Irrational{:γ}) = Sym(sympy["EulerGamma"])
    Base.convert(::Type{Sym}, x::Irrational{:catalan}) = Sym(sympy["Catalan"])
    Base.convert(::Type{Sym}, x::Irrational{:φ}) = (1 + Sym(5)^(1//2))/2
else
    Base.convert(::Type{Sym}, x::Irrational{:π}) = PI
    Base.convert(::Type{Sym}, x::Irrational{:e}) = E
    Base.convert(::Type{Sym}, x::Irrational{:γ}) = Sym(sympy["EulerGamma"])
    Base.convert(::Type{Sym}, x::Irrational{:catalan}) = Sym(sympy["Catalan"])
    Base.convert(::Type{Sym}, x::Irrational{:φ}) = (1 + Sym(5)^(1//2))/2
end
function init_math()
    "PI is a symbolic  π. Using `julia`'s `pi` will give round off errors."
    copy!(PI.x,  sympy["pi"])

    "E is a symbolic  `e`. Using `julia`'s `e` will give round off errors."
    copy!(E.x, Sym(sympy["exp"](1)).x)

    "IM is a symbolic `im`"
    copy!(IM.x, sympy["I"])

    "oo is a symbolic infinity. Example: `integrate(exp(-x), x, 0, oo)`."
    copy!(oo.x, sympy["oo"])
end


text:
using Test
using StaticArrays: SVector
using QuantumLattices.Essentials.FockPackage
using QuantumLattices.Essentials.Spatials: Bond, Point, PID, rcoord, azimuthd
using QuantumLattices.Essentials.DegreesOfFreedom: Table, IDFConfig, OID, Operators, oidtype, otype, script, latexformat
using QuantumLattices.Essentials.Terms: Couplings, @subscript, statistics, abbr
using QuantumLattices.Interfaces: dims, inds, ⊗, ⋅, expand, permute
using QuantumLattices.Prerequisites: Float
using QuantumLattices.Mathematics.VectorSpaces: IsMultiIndexable, MultiIndexOrderStyle
using QuantumLattices.Mathematics.AlgebraOverFields: ID
import QuantumLattices.FockPackage: fockcouplingnambus

@testset "FID" begin
    fid = FID(orbital=1, spin=1)
    @test fid|>typeof|>fieldnames == (:orbital, :spin, :nambu)
    @test fid' == FID(1, 1, 2)
    @test fid'' == FID(1, 1, 1)

    fid = FID(1, 1, 0)
    @test fid' == fid
end

@testset "Fock" begin
    fock = Fock(norbital=1, nspin=2, nnambu=2)
    @test IsMultiIndexable(Fock) == IsMultiIndexable(true)
    @test MultiIndexOrderStyle(Fock) == MultiIndexOrderStyle('C')
    @test dims(fock) == (1, 2, 2)
    @test inds(FID(1, 1, 1), fock) == (1, 1, 1)
    @test FID((1, 1, 1), fock) == FID(1, 1, 1)
    @test fock|>collect == [FID(1, 1, 1), FID(1, 1, 2), FID(1, 2, 1), FID(1, 2, 2)]
end

@testset "FIndex" begin
    @test FIndex|>fieldnames == (:scope, :site, :orbital, :spin, :nambu)
    @test union(PID{Int}, FID) == FIndex{Int}
    @test angle(OID(FIndex(1, 1, 1, 1, 1), [0.0, 0.0], [1.0, 2.0], 1), [[1.0, 0.0], [0.0, 1.0]], [0.1, 0.0]) ≈ 2pi*0.1
    @test angle(OID(FIndex(1, 1, 1, 1, 2), [0.0, 0.0], [1.0, 2.0], 1), [[1.0, 0.0], [0.0, 1.0]], [0.0, 0.2]) ≈ -2pi*0.4
end

@testset "script" begin
    @test script(OID(FIndex('c', 1, 2, 1, 1)), Val(:site)) == 1
    @test script(OID(FIndex('c', 1, 2, 1, 1)), Val(:orbital)) == 2
    @test script(OID(FIndex('c', 1, 2, 3, 1)), Val(:spinint)) == 3
    @test script(OID(FIndex('c', 1, 2, 2, 1)), Val(:spinsym)) == "↑"
    @test script(OID(FIndex('c', 1, 2, 1, 1)), Val(:spinsym)) == "↓"
    @test script(OID(FIndex('c', 1, 2, 3, 1)), Val(:nambu)) == ""
    @test script(OID(FIndex('c', 1, 2, 3, 2)), Val(:nambu)) == "\\dagger"
end

@testset "oidtype" begin
    @test oidtype(FID, Point{2, PID{Int}}, Nothing, Val(true)) == OID{FIndex{Int}, SVector{2, Float}, SVector{2, Float}, Nothing}
    @test oidtype(FID, Point{2, PID{Int}}, Table, Val(true)) == OID{FIndex{Int}, SVector{2, Float}, SVector{2, Float}, Int}
    @test oidtype(FID, Point{2, PID{Int}}, Nothing, Val(false)) == OID{FIndex{Int}, Nothing, Nothing, Nothing}
    @test oidtype(FID, Point{2, PID{Int}}, Table, Val(false)) == OID{FIndex{Int}, Nothing, Nothing, Int}
end

@testset "FockOperator" begin
    @test latexformat(FOperator) == foptdefaultlatex

    opt = FOperator(1.0, (FIndex(1, 1, 1, 1, 2), FIndex(1, 1, 1, 1, 1)))
    @test opt|>isnormalordered
    @test opt|>statistics == opt|>typeof|>statistics == 'F'

    opt = FOperator(1.0, (FIndex(1, 2, 1, 1, 2), FIndex(1, 2, 1, 1, 1), FIndex(1, 1, 1, 2, 2), FIndex(1, 1, 1, 2, 1)))
    @test opt|>isnormalordered == false
    @test repr(opt) == "c^{\\dagger}_{2, 1, ↓}c^{}_{2, 1, ↓}c^{\\dagger}_{1, 1, ↑}c^{}_{1, 1, ↑}"

    opt1 = FOperator(1.5, (FIndex(1, 2, 1, 1, 2), FIndex(1, 2, 1, 1, 1)))
    opt2 = FOperator(2.0, (FIndex(1, 2, 1, 1, 1), FIndex(1, 2, 1, 1, 2)))
    @test opt1*opt2 == nothing

    opt1 = FOperator(1.5, (FIndex(1, 2, 1, 1, 2), FIndex(1, 2, 1, 1, 1)))
    opt2 = FOperator(2.0, (FIndex(1, 2, 1, 1, 2), FIndex(1, 2, 1, 1, 1)))
    @test opt1*opt2 == FOperator(3.0, (FIndex(1, 2, 1, 1, 2), FIndex(1, 2, 1, 1, 1), FIndex(1, 2, 1, 1, 2), FIndex(1, 2, 1, 1, 1)))

    @test latexformat(BOperator) == boptdefaultlatex

    opt = BOperator(1.0, (FIndex(1, 1, 1, 1, 2), FIndex(1, 1, 1, 1, 1)))
    @test opt|>statistics == opt|>typeof|>statistics == 'B'
    @test repr(opt) == "b^{\\dagger}_{1, 1, ↓}b^{}_{1, 1, ↓}"
end

@testset "permute" begin
    i1, i2 = OID(FIndex(1, 1, 1, 1, 2)), OID(FIndex(1, 1, 1, 1, 1))
    @test permute(FOperator, i1, i2) == (FOperator(1), FOperator(-1, ID(i2, i1)))
    @test permute(FOperator, i2, i1) == (FOperator(1), FOperator(-1, ID(i1, i2)))
    @test permute(BOperator, i1, i2) == (BOperator(1), BOperator(1, ID(i2, i1)))
    @test permute(BOperator, i2, i1) == (BOperator(-1), BOperator(1, ID(i1, i2)))

    i1, i2 = OID(FIndex(1, 1, 1, 2, 2)), OID(FIndex(1, 1, 1, 1, 1))
    @test permute(FOperator, i1, i2) == (FOperator(-1, ID(i2, i1)),)
    @test permute(FOperator, i2, i1) == (FOperator(-1, ID(i1, i2)),)
    @test permute(BOperator, i1, i2) == (BOperator(1, ID(i2, i1)),)
    @test permute(BOperator, i2, i1) == (BOperator(1, ID(i1, i2)),)
end

@testset "FCID" begin
    @test FCID(atom=1, orbital=1, spin=1, nambu=1)|>typeof|>fieldnames == (:center, :atom, :orbital, :spin, :nambu, :obsub, :spsub)
end

@testset "FockCoupling" begin
    @test FockCoupling{2}(1.0)|>string == "FockCoupling{2}(value=1.0)"
    @test FockCoupling{2}(1.0, atoms=(1, 1))|>string == "FockCoupling{2}(value=1.0, atoms=(1, 1))"
    @test FockCoupling{2}(1.0, atoms=(1, 1), spins=(1, 2))|>string == "FockCoupling{2}(value=1.0, atoms=(1, 1), spins=(1, 2))"
    @test FockCoupling{2}(2.0)|>repr == "2.0 {2}"

    fc1 = FockCoupling{2}(1.5, atoms=(2, 1), spins=(@subscript (x, 1)), centers=(1, 2))
    fc2 = FockCoupling{2}(2.0, atoms=(1, 2), orbitals=(@subscript (x, y) with x < y), centers=(1, 2))
    @test fc1|>repr == "1.5 sl(2, 1)⊗sp(x, 1)@(1, 2) with (*, $(fc1.id[1].spsub)) && (*, $(fc1.id[2].spsub))"
    @test fc2|>repr == "2.0 sl(1, 2)⊗ob(x, y)@(1, 2) with ($(fc2.id[1].obsub), *) && ($(fc2.id[2].obsub), *)"
    fc = fc1 * fc2
    @test fc|>repr == "3.0 sl(2, 1, 1, 2)⊗ob(*, *, x, y)⊗sp(x, 1, *, *)@(1, 2, 1, 2) with (*, $(fc1.id[1].spsub)) && (*, $(fc1.id[2].spsub)) && ($(fc2.id[1].obsub), *) && ($(fc2.id[2].obsub), *)"

    fc1 = FockCoupling{2}(1.5, spins=(@subscript (x, 1)), centers=(1, 2))
    fc2 = FockCoupling{2}(2.0, orbitals=(@subscript (x, y) with x < y), centers=(1, 2))
    fc = fc1 ⊗ fc2
    @test fc|>repr == "3.0 ob(x, y)⊗sp(x, 1)@(1, 2) with ($(fc2.id[1].obsub), $(fc1.id[2].spsub)) && ($(fc2.id[1].obsub), $(fc1.id[2].spsub))"

    fc1 = FockCoupling{2}(1.5, atoms=(2, 1), centers=(1, 2))
    fc2 = FockCoupling{2}(2.0, atoms=(1, 2), centers=(1, 2))
    fc = fc1 ⋅ fc2
    @test fc|>repr == "3.0 sl(2, 2)@(1, 2)"

    ex = expand(FockCoupling{2}(2.0, atoms=(1, 1)), PID(1, 1), Fock(atom=2, norbital=2, nspin=2, nnambu=2))
    @test collect(ex) == []

    ex = expand(FockCoupling{2}(2.0, atoms=(1, 2), orbitals=(1, 2)), (PID(1, 1), PID(1, 2)), (Fock(atom=1, norbital=2, nspin=2, nnambu=2), Fock(atom=2, norbital=2, nspin=2, nnambu=2)))
    @test IsMultiIndexable(typeof(ex)) == IsMultiIndexable(true)
    @test MultiIndexOrderStyle(typeof(ex)) == MultiIndexOrderStyle('C')
    @test dims(ex) == (1, 2)
    @test collect(ex) == [  (2.0, (FIndex(1, 1, 1, 1, 2), FIndex(1, 2, 2, 1, 1))),
                            (2.0, (FIndex(1, 1, 1, 2, 2), FIndex(1, 2, 2, 2, 1)))
    ]

    ex = expand(FockCoupling{4}(2.0, centers=(1, 1, 1, 1), spins=(2, 2, 1, 1)), PID(1, 1), Fock(atom=1, norbital=2, nspin=2, nnambu=2))
    @test dims(ex) == (2, 1)
    @test collect(ex) == [  (2.0, (FIndex(1, 1, 1, 2, 2), FIndex(1, 1, 1, 2, 1), FIndex(1, 1, 1, 1, 2), FIndex(1, 1, 1, 1, 1))),
                            (2.0, (FIndex(1, 1, 2, 2, 2), FIndex(1, 1, 2, 2, 1), FIndex(1, 1, 2, 1, 2), FIndex(1, 1, 2, 1, 1)))
    ]

    ex = expand(FockCoupling{4}(2.0, orbitals=(@subscript (α, α, β, β) with α < β), spins=(2, 1, 1, 2), nambus=(2, 2, 1, 1)), PID(1, 1), Fock(atom=1, norbital=3, nspin=2, nnambu=2))
    @test dims(ex) == (3, 1)
    @test collect(ex) == [  (2.0, (FIndex(1, 1, 1, 2, 2), FIndex(1, 1, 1, 1, 2), FIndex(1, 1, 2, 1, 1), FIndex(1, 1, 2, 2, 1))),
                            (2.0, (FIndex(1, 1, 1, 2, 2), FIndex(1, 1, 1, 1, 2), FIndex(1, 1, 3, 1, 1), FIndex(1, 1, 3, 2, 1))),
                            (2.0, (FIndex(1, 1, 2, 2, 2), FIndex(1, 1, 2, 1, 2), FIndex(1, 1, 3, 1, 1), FIndex(1, 1, 3, 2, 1)))
    ]

    fc1 = FockCoupling{2}(1.0, spins=(2, 2))
    fc2 = FockCoupling{2}(-1.0, spins=(1, 1))
    ex = expand(fc1*fc2, PID(1, 1), Fock(atom=1, norbital=2, nspin=2, nnambu=2))
    @test dims(ex) == (4, 1)
    @test collect(ex) == [  (-1.0, (FIndex(1, 1, 1, 2, 2), FIndex(1, 1, 1, 2, 1), FIndex(1, 1, 1, 1, 2), FIndex(1, 1, 1, 1, 1))),
                            (-1.0, (FIndex(1, 1, 1, 2, 2), FIndex(1, 1, 1, 2, 1), FIndex(1, 1, 2, 1, 2), FIndex(1, 1, 2, 1, 1))),
                            (-1.0, (FIndex(1, 1, 2, 2, 2), FIndex(1, 1, 2, 2, 1), FIndex(1, 1, 1, 1, 2), FIndex(1, 1, 1, 1, 1))),
                            (-1.0, (FIndex(1, 1, 2, 2, 2), FIndex(1, 1, 2, 2, 1), FIndex(1, 1, 2, 1, 2), FIndex(1, 1, 2, 1, 1)))
    ]
end

@testset "fockcouplingnambus" begin
    @test fockcouplingnambus(nothing, ('*', '*', '*', '*'), (1, 1, 2, 2)) == (0, 0, 2, 1)
    @test fockcouplingnambus(nothing, (0, 0, 1, 1), (1, 1, 2, 2)) == (0, 0, 1, 1)
    @test fockcouplingnambus(Val(:Pairing), ('*', '*'), (2, 2)) == (1, 1)
    @test fockcouplingnambus(Val(:Pairing), (1, 2), (2, 2)) == (1, 2)
end

@testset "σ⁰" begin
    @test σ⁰("sp") == FockCoupling{2}(1, spins=(1, 1))+FockCoupling{2}(1, spins=(2, 2))
    @test σ⁰("ob") == FockCoupling{2}(1, orbitals=(1, 1))+FockCoupling{2}(1, orbitals=(2, 2))
    @test σ⁰("sl") == FockCoupling{2}(1, atoms=(1, 1))+FockCoupling{2}(1, atoms=(2, 2))
    @test σ⁰("ph") == FockCoupling{2}(1, nambus=(1, 2))+FockCoupling{2}(1, nambus=(2, 1))
end

@testset "σˣ" begin
    @test σˣ("sp") == FockCoupling{2}(1, spins=(1, 2))+FockCoupling{2}(1, spins=(2, 1))
    @test σˣ("ob") == FockCoupling{2}(1, orbitals=(1, 2))+FockCoupling{2}(1, orbitals=(2, 1))
    @test σˣ("sl") == FockCoupling{2}(1, atoms=(1, 2))+FockCoupling{2}(1, atoms=(2, 1))
    @test σˣ("ph") == FockCoupling{2}(1, nambus=(1, 1))+FockCoupling{2}(1, nambus=(2, 2))
end

@testset "σʸ" begin
    @test σʸ("sp") == FockCoupling{2}(1im, spins=(1, 2))+FockCoupling{2}(-1im, spins=(2, 1))
    @test σʸ("ob") == FockCoupling{2}(1im, orbitals=(1, 2))+FockCoupling{2}(-1im, orbitals=(2, 1))
    @test σʸ("sl") == FockCoupling{2}(1im, atoms=(1, 2))+FockCoupling{2}(-1im, atoms=(2, 1))
    @test σʸ("ph") == FockCoupling{2}(1im, nambus=(1, 1))+FockCoupling{2}(-1im, nambus=(2, 2))
end

@testset "σᶻ" begin
    @test σᶻ("sp") == FockCoupling{2}(-1, spins=(1, 1))+FockCoupling{2}(1, spins=(2, 2))
    @test σᶻ("ob") == FockCoupling{2}(-1, orbitals=(1, 1))+FockCoupling{2}(1, orbitals=(2, 2))
    @test σᶻ("sl") == FockCoupling{2}(-1, atoms=(1, 1))+FockCoupling{2}(1, atoms=(2, 2))
    @test σᶻ("ph") == FockCoupling{2}(-1, nambus=(1, 2))+FockCoupling{2}(1, nambus=(2, 1))
end

@testset "σ⁺" begin
    @test σ⁺("sp") == Couplings(FockCoupling{2}(1, spins=(2, 1)))
    @test σ⁺("ob") == Couplings(FockCoupling{2}(1, orbitals=(2, 1)))
    @test σ⁺("sl") == Couplings(FockCoupling{2}(1, atoms=(2, 1)))
    @test σ⁺("ph") == Couplings(FockCoupling{2}(1, nambus=(2, 2)))
end

@testset "σ⁻" begin
    @test σ⁻("sp") == Couplings(FockCoupling{2}(1, spins=(1, 2)))
    @test σ⁻("ob") == Couplings(FockCoupling{2}(1, orbitals=(1, 2)))
    @test σ⁻("sl") == Couplings(FockCoupling{2}(1, atoms=(1, 2)))
    @test σ⁻("ph") == Couplings(FockCoupling{2}(1, nambus=(1, 1)))
end

@testset "fockcoupling" begin
    fc = fc"1.0 sl(1, 1, 1, 1)⊗ob(α, α, β, β)⊗sp(σ, γ, σ, γ)⊗ph(2, 1, 2, 1)@(1, 1, 2, 2) with (α<β, σ≠γ)"
    ob, sp = fc.obsubscripts[1].identifier, fc.spsubscripts[1].identifier
    @test repr(fc) == "1.0 sl(1, 1, 1, 1)⊗ob(α, α, β, β)⊗sp(σ, γ, σ, γ)⊗ph(2, 1, 2, 1)@(1, 1, 2, 2) with ($ob, $sp) && ($ob, $sp) && ($ob, $sp) && ($ob, $sp)"

    fc = fc"1.0 sl(1, 1, 1, 1)⊗ob(α, α, β, β)⊗sp(σ, γ, σ, γ)⊗ph(2, 1, 2, 1)@(1, 1, 2, 2)"
    ob, sp = fc.obsubscripts[1].identifier, fc.spsubscripts[1].identifier
    @test repr(fc) == "1.0 sl(1, 1, 1, 1)⊗ob(α, α, β, β)⊗sp(σ, γ, σ, γ)⊗ph(2, 1, 2, 1)@(1, 1, 2, 2) with ($ob, $sp) && ($ob, $sp) && ($ob, $sp) && ($ob, $sp)"

    fc = fc"1.0 sl(1, 1, 1, 1)⊗ob(α, α, β, β)⊗sp(σ, γ, σ, γ)⊗ph(2, 1, 2, 1)"
    ob, sp = fc.obsubscripts[1].identifier, fc.spsubscripts[1].identifier
    @test repr(fc) == "1.0 sl(1, 1, 1, 1)⊗ob(α, α, β, β)⊗sp(σ, γ, σ, γ)⊗ph(2, 1, 2, 1) with ($ob, $sp) && ($ob, $sp) && ($ob, $sp) && ($ob, $sp)"

    fc = fc"1.0 sl(1, 1, 1, 1)⊗ob(α, α, β, β)⊗sp(σ, γ, σ, γ)⊗ph(2, 1, 2, 1) with (α<β, σ≠γ)"
    ob, sp = fc.obsubscripts[1].identifier, fc.spsubscripts[1].identifier
    @test repr(fc) == "1.0 sl(1, 1, 1, 1)⊗ob(α, α, β, β)⊗sp(σ, γ, σ, γ)⊗ph(2, 1, 2, 1) with ($ob, $sp) && ($ob, $sp) && ($ob, $sp) && ($ob, $sp)"

    fc = fc"1.0 sl(1, 1, 1, 1)⊗ob(α, α, β, β)⊗sp(σ, γ, σ, γ)⊗ph(2, 1, 2, 1) with α<β, σ≠γ"
    ob, sp = fc.obsubscripts[1].identifier, fc.spsubscripts[1].identifier
    @test repr(fc) == "1.0 sl(1, 1, 1, 1)⊗ob(α, α, β, β)⊗sp(σ, γ, σ, γ)⊗ph(2, 1, 2, 1) with ($ob, $sp) && ($ob, $sp) && ($ob, $sp) && ($ob, $sp)"

    fc = fc"1.0 sl(1, 1, 1, 1)⊗ob(α, α, β, β)⊗sp(σ, γ, σ, γ)⊗ph(2, 1, 2, 1) with (α<β, *)"
    ob, sp = fc.obsubscripts[1].identifier, fc.spsubscripts[1].identifier
    @test repr(fc) == "1.0 sl(1, 1, 1, 1)⊗ob(α, α, β, β)⊗sp(σ, γ, σ, γ)⊗ph(2, 1, 2, 1) with ($ob, $sp) && ($ob, $sp) && ($ob, $sp) && ($ob, $sp)"

    fc = fc"1.0 sl(1, 1, 1, 1)⊗ob(α, α, β, β)⊗ph(2, 1, 2, 1) with α<β"
    ob, sp = fc.obsubscripts[1].identifier, fc.spsubscripts[1].identifier
    @test repr(fc) == "1.0 sl(1, 1, 1, 1)⊗ob(α, α, β, β)⊗ph(2, 1, 2, 1) with ($ob, $sp) && ($ob, $sp) && ($ob, $sp) && ($ob, $sp)"

    fc = fc"1.0 sl(1, 1, 1, 1)⊗ob(α, α, β, β)⊗sp(2, 1, 2, 1) with α<β"
    ob, sp = fc.obsubscripts[1].identifier, fc.spsubscripts[1].identifier
    @test repr(fc) == "1.0 sl(1, 1, 1, 1)⊗ob(α, α, β, β)⊗sp(2, 1, 2, 1) with ($ob, $sp) && ($ob, $sp) && ($ob, $sp) && ($ob, $sp)"

    fc = fc"1.0 sl(1, 1, 1, 1)⊗ob(1, 1, 1, 1)⊗ph(2, 1, 2, 1)"
    @test repr(fc) == "1.0 sl(1, 1, 1, 1)⊗ob(1, 1, 1, 1)⊗ph(2, 1, 2, 1)"

    fc = fc"1.0 sl(1, 1, 1, 1)⊗ph(2, 1, 2, 1)"
    @test repr(fc) == "1.0 sl(1, 1, 1, 1)⊗ph(2, 1, 2, 1)"

    fc = fc"1.0 @(1, 1, 1, 1)"
    @test repr(fc) == "1.0 @(1, 1, 1, 1)"

    fc = fc"1.0im {2}"
    @test repr(fc) == "1.0im {2}"
end

@testset "fockcouplings" begin
    @test σ⁰"sp" == σ⁰("sp") && σ⁰"sp@(1, 2)" == σ⁰("sp", centers=(1, 2))
    @test σ⁰"ob" == σ⁰("ob") && σ⁰"ob@(1, 2)" == σ⁰("ob", centers=(1, 2))
    @test σ⁰"sl" == σ⁰("sl") && σ⁰"sl@(1, 2)" == σ⁰("sl", centers=(1, 2))
    @test σ⁰"ph" == σ⁰("ph") && σ⁰"ph@(1, 2)" == σ⁰("ph", centers=(1, 2))

    @test σˣ"sp" == σˣ("sp") && σˣ"sp@(1, 2)" == σˣ("sp", centers=(1, 2))
    @test σˣ"ob" == σˣ("ob") && σˣ"ob@(1, 2)" == σˣ("ob", centers=(1, 2))
    @test σˣ"sl" == σˣ("sl") && σˣ"sl@(1, 2)" == σˣ("sl", centers=(1, 2))
    @test σˣ"ph" == σˣ("ph") && σˣ"ph@(1, 2)" == σˣ("ph", centers=(1, 2))

    @test σʸ"sp" == σʸ("sp") && σʸ"sp@(1, 2)" == σʸ("sp", centers=(1, 2))
    @test σʸ"ob" == σʸ("ob") && σʸ"ob@(1, 2)" == σʸ("ob", centers=(1, 2))
    @test σʸ"sl" == σʸ("sl") && σʸ"sl@(1, 2)" == σʸ("sl", centers=(1, 2))
    @test σʸ"ph" == σʸ("ph") && σʸ"ph@(1, 2)" == σʸ("ph", centers=(1, 2))

    @test σᶻ"sp" == σᶻ("sp") && σᶻ"sp@(1, 2)" == σᶻ("sp", centers=(1, 2))
    @test σᶻ"ob" == σᶻ("ob") && σᶻ"ob@(1, 2)" == σᶻ("ob", centers=(1, 2))
    @test σᶻ"sl" == σᶻ("sl") && σᶻ"sl@(1, 2)" == σᶻ("sl", centers=(1, 2))
    @test σᶻ"ph" == σᶻ("ph") && σᶻ"ph@(1, 2)" == σᶻ("ph", centers=(1, 2))

    @test σ⁺"sp" == σ⁺("sp") && σ⁺"sp@(1, 2)" == σ⁺("sp", centers=(1, 2))
    @test σ⁺"ob" == σ⁺("ob") && σ⁺"ob@(1, 2)" == σ⁺("ob", centers=(1, 2))
    @test σ⁺"sl" == σ⁺("sl") && σ⁺"sl@(1, 2)" == σ⁺("sl", centers=(1, 2))
    @test σ⁺"ph" == σ⁺("ph") && σ⁺"ph@(1, 2)" == σ⁺("ph", centers=(1, 2))

    @test σ⁻"sp" == σ⁻("sp") && σ⁻"sp@(1, 2)" == σ⁻("sp", centers=(1, 2))
    @test σ⁻"ob" == σ⁻("ob") && σ⁻"ob@(1, 2)" == σ⁻("ob", centers=(1, 2))
    @test σ⁻"sl" == σ⁻("sl") && σ⁻"sl@(1, 2)" == σ⁻("sl", centers=(1, 2))
    @test σ⁻"ph" == σ⁻("ph") && σ⁻"ph@(1, 2)" == σ⁻("ph", centers=(1, 2))
end

@testset "Onsite" begin
    term = Onsite{'F'}(:mu, 1.5)
    @test term|>abbr == :st
    @test otype(term|>typeof, OID{FIndex{Int}, SVector{2, Float}, SVector{2, Float}, Nothing}) == FOperator{Float, ID{OID{FIndex{Int}, SVector{2, Float}, SVector{2, Float}, Nothing}, 2}}
    @test otype(term|>typeof, OID{FIndex{Int}, SVector{2, Float}, SVector{2, Float}, Int}) == FOperator{Float, ID{OID{FIndex{Int}, SVector{2, Float}, SVector{2, Float}, Int}, 2}}

    term = Onsite{'B'}(:mu, 1.5, couplings=σˣ("sp")⊗σᶻ("ob"), modulate=true)
    @test term|>abbr == :st
    @test otype(term|>typeof, OID{FIndex{Int}, SVector{2, Float}, SVector{2, Float}, Nothing}) == BOperator{Float, ID{OID{FIndex{Int}, SVector{2, Float}, SVector{2, Float}, Nothing}, 2}}
    @test otype(term|>typeof, OID{FIndex{Int}, SVector{2, Float}, SVector{2, Float}, Int}) == BOperator{Float, ID{OID{FIndex{Int}, SVector{2, Float}, SVector{2, Float}, Int}, 2}}

    point = Point(PID('a', 1), (0.5, 0.5), (0.0, 0.0))
    config = IDFConfig{Fock}(pid->Fock(atom=pid.site%2, norbital=2, nspin=2, nnambu=2), [point.pid])
    table = Table(config, usualfockindextotuple)

    term = Onsite{'F'}(:mu, 1.5, couplings=σˣ("sp")⊗σᶻ("ob"), modulate=true)
    operators = Operators(  FOperator(+1.5, ID(OID(FIndex('a', 1, 2, 2, 2), [0.5, 0.5], [0.0, 0.0], 4), OID(FIndex('a', 1, 2, 1, 1), [0.5, 0.5], [0.0, 0.0], 3))),
                            FOperator(-1.5, ID(OID(FIndex('a', 1, 1, 1, 2), [0.5, 0.5], [0.0, 0.0], 1), OID(FIndex('a', 1, 1, 2, 1), [0.5, 0.5], [0.0, 0.0], 2)))
    )
    @test expand(term, point, config, table, true) == operators
    @test expand(term, point, config, table, false) == operators+operators'
    operators = Operators(  FOperator(+1.5, ID(OID(FIndex('a', 1, 2, 2, 2), [0.5, 0.5], [0.0, 0.0], nothing), OID(FIndex('a', 1, 2, 1, 1), [0.5, 0.5], [0.0, 0.0], nothing))),
                            FOperator(-1.5, ID(OID(FIndex('a', 1, 1, 1, 2), [0.5, 0.5], [0.0, 0.0], nothing), OID(FIndex('a', 1, 1, 2, 1), [0.5, 0.5], [0.0, 0.0], nothing)))
    )
    @test expand(term, point, config, nothing, true) == operators
    @test expand(term, point, config, nothing, false) == operators+operators'

    term = Onsite{'F'}(:mu, 1.5, couplings=σᶻ("sp")⊗σᶻ("ob"), modulate=true)
    operators = Operators(  FOperator(-0.75, ID(OID(FIndex('a', 1, 2, 1, 2), [0.5, 0.5], [0.0, 0.0], 3), OID(FIndex('a', 1, 2, 1, 1), [0.5, 0.5], [0.0, 0.0], 3))),
                        FOperator(-0.75, ID(OID(FIndex('a', 1, 1, 2, 2), [0.5, 0.5], [0.0, 0.0], 2), OID(FIndex('a', 1, 1, 2, 1), [0.5, 0.5], [0.0, 0.0], 2))),
                        FOperator(+0.75, ID(OID(FIndex('a', 1, 2, 2, 2), [0.5, 0.5], [0.0, 0.0], 4), OID(FIndex('a', 1, 2, 2, 1), [0.5, 0.5], [0.0, 0.0], 4))),
                        FOperator(+0.75, ID(OID(FIndex('a', 1, 1, 1, 2), [0.5, 0.5], [0.0, 0.0], 1), OID(FIndex('a', 1, 1, 1, 1), [0.5, 0.5], [0.0, 0.0], 1)))
    )
    @test expand(term, point, config, table, true) == operators
    @test expand(term, point, config, table, false) == operators+operators'
end

@testset "Hopping" begin
    bond = Bond(1, Point(PID('a', 1), (0.5, 0.5), (0.0, 0.0)), Point(PID('b', 2), (0.0, 0.0), (0.0, 0.0)))
    config = IDFConfig{Fock}(pid->Fock(atom=pid.site%2, norbital=2, nspin=2, nnambu=2), [bond.spoint.pid, bond.epoint.pid])
    table = Table(config, usualfockindextotuple)
    term = Hopping{'F'}(:t, 1.5, 1)
    operators = Operators(  FOperator(1.5, ID(OID(FIndex('b', 2, 2, 2, 2), [0.0, 0.0], [0.0, 0.0], 8), OID(FIndex('a', 1, 2, 2, 1), [0.5, 0.5], [0.0, 0.0], 4))),
                            FOperator(1.5, ID(OID(FIndex('b', 2, 2, 1, 2), [0.0, 0.0], [0.0, 0.0], 7), OID(FIndex('a', 1, 2, 1, 1), [0.5, 0.5], [0.0, 0.0], 3))),
                            FOperator(1.5, ID(OID(FIndex('b', 2, 1, 1, 2), [0.0, 0.0], [0.0, 0.0], 5), OID(FIndex('a', 1, 1, 1, 1), [0.5, 0.5], [0.0, 0.0], 1))),
                            FOperator(1.5, ID(OID(FIndex('b', 2, 1, 2, 2), [0.0, 0.0], [0.0, 0.0], 6), OID(FIndex('a', 1, 1, 2, 1), [0.5, 0.5], [0.0, 0.0], 2)))
    )
    @test term|>abbr == :hp
    @test expand(term, bond, config, table, true) == operators
    @test expand(term, bond, config, table, false) == operators+operators'
end

@testset "Pairing" begin
    bond = Bond(1, Point(PID('a', 1), (0.5, 0.5), (0.0, 0.0)), Point(PID('b', 2), (0.0, 0.0), (0.0, 0.0)))
    config = IDFConfig{Fock}(pid->Fock(atom=pid.site%2, norbital=1, nspin=2, nnambu=2), [bond.spoint.pid, bond.epoint.pid])
    table = Table(config, nambufockindextotuple)
    term = Pairing{'F'}(:Δ, 1.5, 1, couplings=FockCoupling{2}(spins=(2, 2)), amplitude=bond->(bond|>rcoord|>azimuthd ≈ 45 ? 1 : -1))
    operators = Operators(  FOperator(-1.5, ID(OID(FIndex('b', 2, 1, 2, 1), [0.0, 0.0], [0.0, 0.0], 6), OID(FIndex('a', 1, 1, 2, 1), [0.5, 0.5], [0.0, 0.0], 2))),
                            FOperator(+1.5, ID(OID(FIndex('a', 1, 1, 2, 1), [0.5, 0.5], [0.0, 0.0], 2), OID(FIndex('b', 2, 1, 2, 1), [0.0, 0.0], [0.0, 0.0], 6)))
    )
    @test term|>abbr == :pr
    @test expand(term, bond, config, table, true) == operators
    @test expand(term, bond, config, table, false) == operators+operators'

    point = Point(PID('a', 1), (0.5, 0.5), (0.0, 0.0))
    config = IDFConfig{Fock}(pid->Fock(atom=pid.site % 2, norbital=1, nspin=2, nnambu=2), [point.pid])
    table = Table(config, nambufockindextotuple)
    term = Pairing{'F'}(:Δ, 1.5, 0, couplings=FockCoupling{2}(spins=(2, 1))-FockCoupling{2}(spins=(1, 2)))
    operators = Operators(  FOperator(+1.5, ID(OID(FIndex('a', 1, 1, 2, 1), [0.5, 0.5], [0.0, 0.0], 2), OID(FIndex('a', 1, 1, 1, 1), [0.5, 0.5], [0.0, 0.0], 1))),
                            FOperator(-1.5, ID(OID(FIndex('a', 1, 1, 1, 1), [0.5, 0.5], [0.0, 0.0], 1), OID(FIndex('a', 1, 1, 2, 1), [0.5, 0.5], [0.0, 0.0], 2)))
    )
    @test term|>abbr == :pr
    @test expand(term, point, config, table, true) == operators
    @test expand(term, point, config, table, false) == operators+operators'
end

@testset "Hubbard" begin
    point = Point(PID('a', 1), (0.5, 0.5), (0.0, 0.0))
    config = IDFConfig{Fock}(pid->Fock(atom=pid.site%2, norbital=2, nspin=2, nnambu=2), [point.pid])
    table = Table(config, usualfockindextotuple)
    term = Hubbard{'F'}(:H, 2.5)
    operators = Operators(  FOperator(1.25, ID( OID(FIndex('a', 1, 1, 2, 2), [0.5, 0.5], [0.0, 0.0], 2), OID(FIndex('a', 1, 1, 2, 1), [0.5, 0.5], [0.0, 0.0], 2),
                                                OID(FIndex('a', 1, 1, 1, 2), [0.5, 0.5], [0.0, 0.0], 1), OID(FIndex('a', 1, 1, 1, 1), [0.5, 0.5], [0.0, 0.0], 1))),
                            FOperator(1.25, ID( OID(FIndex('a', 1, 2, 2, 2), [0.5, 0.5], [0.0, 0.0], 4), OID(FIndex('a', 1, 2, 2, 1), [0.5, 0.5], [0.0, 0.0], 4),
                                                OID(FIndex('a', 1, 2, 1, 2), [0.5, 0.5], [0.0, 0.0], 3), OID(FIndex('a', 1, 2, 1, 1), [0.5, 0.5], [0.0, 0.0], 3)))
    )
    @test term|>abbr == :hb
    @test expand(term, point, config, table, true) == operators
    @test expand(term, point, config, table, false) == operators*2
end

@testset "InterOrbitalInterSpin" begin
    point = Point(PID('a', 1), (0.5, 0.5), (0.0, 0.0))
    config = IDFConfig{Fock}(pid->Fock(atom=pid.site%2, norbital=2, nspin=2, nnambu=2), [point.pid])
    table = Table(config, usualfockindextotuple)
    term = InterOrbitalInterSpin{'F'}(:H, 2.5)
    operators = Operators(  FOperator(1.25, ID( OID(FIndex('a', 1, 1, 2, 2), [0.5, 0.5], [0.0, 0.0], 2), OID(FIndex('a', 1, 1, 2, 1), [0.5, 0.5], [0.0, 0.0], 2),
                                                OID(FIndex('a', 1, 2, 1, 2), [0.5, 0.5], [0.0, 0.0], 3), OID(FIndex('a', 1, 2, 1, 1), [0.5, 0.5], [0.0, 0.0], 3))),
                            FOperator(1.25, ID( OID(FIndex('a', 1, 1, 1, 2), [0.5, 0.5], [0.0, 0.0], 1), OID(FIndex('a', 1, 1, 1, 1), [0.5, 0.5], [0.0, 0.0], 1),
                                                OID(FIndex('a', 1, 2, 2, 2), [0.5, 0.5], [0.0, 0.0], 4), OID(FIndex('a', 1, 2, 2, 1), [0.5, 0.5], [0.0, 0.0], 4)))
    )
    @test term|>abbr == :nons
    @test expand(term, point, config, table, true) == operators
    @test expand(term, point, config, table, false) == operators*2
end

@testset "InterOrbitalIntraSpin" begin
    point = Point(PID('a', 1), (0.5, 0.5), (0.0, 0.0))
    config = IDFConfig{Fock}(pid->Fock(atom=pid.site%2, norbital=2, nspin=2, nnambu=2), [point.pid])
    table = Table(config, usualfockindextotuple)
    term = InterOrbitalIntraSpin{'F'}(:H, 2.5)
    operators = Operators(  FOperator(1.25, ID( OID(FIndex('a', 1, 1, 1, 2), [0.5, 0.5], [0.0, 0.0], 1), OID(FIndex('a', 1, 1, 1, 1), [0.5, 0.5], [0.0, 0.0], 1),
                                                OID(FIndex('a', 1, 2, 1, 2), [0.5, 0.5], [0.0, 0.0], 3), OID(FIndex('a', 1, 2, 1, 1), [0.5, 0.5], [0.0, 0.0], 3))),
                            FOperator(1.25, ID( OID(FIndex('a', 1, 1, 2, 2), [0.5, 0.5], [0.0, 0.0], 2), OID(FIndex('a', 1, 1, 2, 1), [0.5, 0.5], [0.0, 0.0], 2),
                                                OID(FIndex('a', 1, 2, 2, 2), [0.5, 0.5], [0.0, 0.0], 4), OID(FIndex('a', 1, 2, 2, 1), [0.5, 0.5], [0.0, 0.0], 4)))
    )
    @test term|>abbr == :noes
    @test expand(term, point, config, table, true) == operators
    @test expand(term, point, config, table, false) == operators*2
end

@testset "SpinFlip" begin
    point = Point(PID('a', 1), (0.5, 0.5), (0.0, 0.0))
    config = IDFConfig{Fock}(pid->Fock(atom=pid.site%2, norbital=2, nspin=2, nnambu=2), [point.pid])
    table = Table(config, usualfockindextotuple)
    term = SpinFlip{'F'}(:H, 2.5)
    operators = Operators(FOperator(2.5, ID(    OID(FIndex('a', 1, 1, 2, 2), [0.5, 0.5], [0.0, 0.0], 2), OID(FIndex('a', 1, 2, 1, 2), [0.5, 0.5], [0.0, 0.0], 3),
                                                OID(FIndex('a', 1, 1, 1, 1), [0.5, 0.5], [0.0, 0.0], 1), OID(FIndex('a', 1, 2, 2, 1), [0.5, 0.5], [0.0, 0.0], 4)))
    )
    @test term|>abbr == :sf
    @test expand(term, point, config, table, true) == operators
    @test expand(term, point, config, table, false) == operators+operators'
end

@testset "PairHopping" begin
    point = Point(PID('a', 1), (0.5, 0.5), (0.0, 0.0))
    config = IDFConfig{Fock}(pid->Fock(atom=pid.site%2, norbital=2, nspin=2, nnambu=2), [point.pid])
    table = Table(config, usualfockindextotuple)
    term = PairHopping{'F'}(:H, 2.5)
    operators = Operators(FOperator(2.5, ID(    OID(FIndex('a', 1, 1, 2, 2), [0.5, 0.5], [0.0, 0.0], 2), OID(FIndex('a', 1, 1, 1, 2), [0.5, 0.5], [0.0, 0.0], 1),
                                                OID(FIndex('a', 1, 2, 1, 1), [0.5, 0.5], [0.0, 0.0], 3), OID(FIndex('a', 1, 2, 2, 1), [0.5, 0.5], [0.0, 0.0], 4)))
    )
    @test term|>abbr == :ph
    @test expand(term, point, config, table, true) == operators
    @test expand(term, point, config, table, false) == operators+operators'
end

@testset "Coulomb" begin
    bond = Bond(1, Point(PID('a', 1), (0.5, 0.5), (0.0, 0.0)), Point(PID('b', 2), (0.0, 0.0), (0.0, 0.0)))
    config = IDFConfig{Fock}(pid->Fock(atom=pid.site%2, norbital=1, nspin=2, nnambu=2), [bond.spoint.pid, bond.epoint.pid])
    table = Table(config, usualfockindextotuple)

    term = Coulomb{'F'}(:V, 2.5, 1, couplings=σᶻ("sp")*σᶻ("sp"))
    operators = Operators(  FOperator(-1.25, ID(OID(FIndex('b', 2, 1, 1, 2), [0.0, 0.0], [0.0, 0.0], 3), OID(FIndex('b', 2, 1, 1, 1), [0.0, 0.0], [0.0, 0.0], 3),
                                                OID(FIndex('a', 1, 1, 2, 2), [0.5, 0.5], [0.0, 0.0], 2), OID(FIndex('a', 1, 1, 2, 1), [0.5, 0.5], [0.0, 0.0], 2))),
                            FOperator(+1.25, ID(OID(FIndex('b', 2, 1, 1, 2), [0.0, 0.0], [0.0, 0.0], 3), OID(FIndex('b', 2, 1, 1, 1), [0.0, 0.0], [0.0, 0.0], 3),
                                                OID(FIndex('a', 1, 1, 1, 2), [0.5, 0.5], [0.0, 0.0], 1), OID(FIndex('a', 1, 1, 1, 1), [0.5, 0.5], [0.0, 0.0], 1))),
                            FOperator(-1.25, ID(OID(FIndex('b', 2, 1, 2, 2), [0.0, 0.0], [0.0, 0.0], 4), OID(FIndex('b', 2, 1, 2, 1), [0.0, 0.0], [0.0, 0.0], 4),
                                                OID(FIndex('a', 1, 1, 1, 2), [0.5, 0.5], [0.0, 0.0], 1), OID(FIndex('a', 1, 1, 1, 1), [0.5, 0.5], [0.0, 0.0], 1))),
                            FOperator(+1.25, ID(OID(FIndex('b', 2, 1, 2, 2), [0.0, 0.0], [0.0, 0.0], 4), OID(FIndex('b', 2, 1, 2, 1), [0.0, 0.0], [0.0, 0.0], 4),
                                                OID(FIndex('a', 1, 1, 2, 2), [0.5, 0.5], [0.0, 0.0], 2), OID(FIndex('a', 1, 1, 2, 1), [0.5, 0.5], [0.0, 0.0], 2)))
    )
    @test term|>abbr == :cl
    @test expand(term, bond, config, table, true) == operators
    @test expand(term, bond, config, table, false) == operators+operators'

    term = Coulomb{'F'}(:V, 2.5, 1, couplings=σˣ("sp")*σᶻ("sp"))
    operators = Operators(  FOperator(-2.5, ID( OID(FIndex('b', 2, 1, 2, 2), [0.0, 0.0], [0.0, 0.0], 4), OID(FIndex('b', 2, 1, 1, 1), [0.0, 0.0], [0.0, 0.0], 3),
                                                OID(FIndex('a', 1, 1, 1, 2), [0.5, 0.5], [0.0, 0.0], 1), OID(FIndex('a', 1, 1, 1, 1), [0.5, 0.5], [0.0, 0.0], 1))),
                            FOperator(+2.5, ID( OID(FIndex('b', 2, 1, 1, 2), [0.0, 0.0], [0.0, 0.0], 3), OID(FIndex('b', 2, 1, 2, 1), [0.0, 0.0], [0.0, 0.0], 4),
                                                OID(FIndex('a', 1, 1, 2, 2), [0.5, 0.5], [0.0, 0.0], 2), OID(FIndex('a', 1, 1, 2, 1), [0.5, 0.5], [0.0, 0.0], 2))),
                            FOperator(+2.5, ID( OID(FIndex('b', 2, 1, 2, 2), [0.0, 0.0], [0.0, 0.0], 4), OID(FIndex('b', 2, 1, 1, 1), [0.0, 0.0], [0.0, 0.0], 3),
                                                OID(FIndex('a', 1, 1, 2, 2), [0.5, 0.5], [0.0, 0.0], 2), OID(FIndex('a', 1, 1, 2, 1), [0.5, 0.5], [0.0, 0.0], 2))),
                            FOperator(-2.5, ID( OID(FIndex('b', 2, 1, 1, 2), [0.0, 0.0], [0.0, 0.0], 3), OID(FIndex('b', 2, 1, 2, 1), [0.0, 0.0], [0.0, 0.0], 4),
                                                OID(FIndex('a', 1, 1, 1, 2), [0.5, 0.5], [0.0, 0.0], 1), OID(FIndex('a', 1, 1, 1, 1), [0.5, 0.5], [0.0, 0.0], 1)))
    )
    @test term|>abbr == :cl
    @test expand(term, bond, config, table, true) == operators
    @test expand(term, bond, config, table, false) == operators+operators'
end


text:
"""
    Episode iterators using https://github.com/JuliaML/Reinforce.jl/blob/master/src/episodes/iterators.jl as inspiration for api

"""

using Random

"""
    Episode(env, agent, maxn, rng)

This is a struct for managing the components of an episode iterator. You should only pass a reference to env and agent, while managing the reference seperately.

# Arguments
- `env::AbstractEnvironment`: The Environment (following RLCore interface)
- `agent::AbstractAgent`: The Agent (following RLCore interface)
- `maxn=0`: Max number of steps for the episode.
- `rng=nothing`: The Random Number Generator for the episode (can either be nothing or an AbstractRNG)

"""
mutable struct Episode{E<:AbstractEnvironment, A<:AbstractAgent, F<:AbstractFloat, AR}
    env::E
    agent::A
    total_reward::F # total reward of the episode

    niter::Int      # current step in this episode
    maxn::Int       # max steps in an episode - should be constant during an episode

    rng::AR
end

Episode(env, agent, maxn=0, rng=nothing) =
    Episode(env, agent, get_reward(env), 0, maxn, rng)

function _start!(ep::Episode)
    if ep.rng isa Nothing
        s = start!(ep.env)
        agent_ret = start!(ep.agent, s)

        ep.total_reward = zero(ep.total_reward)
        ep.niter = 1
        return (s, agent_ret), ep.niter
    else
        s = start!(ep.env, ep.rng)
        agent_ret = start!(ep.agent, s, ep.rng)

        ep.total_reward = zero(ep.total_reward)
        ep.niter = 1
        return (s, agent_ret), ep.niter
    end
end

_done(ep::Episode, state) =
  (ep.maxn != 0 && ep.niter >= ep.maxn) || is_terminal(ep.env)

function Base.iterate(ep::Episode, state = _start!(ep))
    _done(ep, state) && return nothing

    env = ep.env
    agent = ep.agent

    s = state[1][1]
    agent_ret = state[1][2]

    action = if agent_ret isa NamedTuple
        agent_ret.action
    else
        agent_ret
    end

    s′, r, t = if ep.rng isa Nothing
        step!(env, action)
    else
        step!(env, action, ep.rng)
    end

    agent_ret = if ep.rng isa Nothing
        step!(agent, s′, r, t)
    else
        step!(agent, s′, r, t, rng)
    end

    ep.total_reward += r
    ep.niter = state[2]

    (s, agent_ret, s′, r, t), ((s′, agent_ret), ep.niter + 1)
end

function run_episode!(f::Base.Callable, env::AbstractEnvironment, agent::AbstractAgent, args...)
  ep = Episode(env, agent, args...)
  for sars in ep
    f(sars)
  end
  ep.total_reward, ep.niter
end

# function episode!(env, agent, max_eps_steps, total_steps, max_steps; callback=nothing, rng::Random.AbstractRNG=nothing)
#     terminal = false

#     # Start Environment
#     s_t = if rng isa Nothing
#         start!(env)
#     else
#         start!(env, rng)
#     end

#     # Start Agent
#     agent_ret = if rng isa Nothing
#         start!(agent, s_t)
#     else
#         start!(agent, s_t, rng)
#     end

#     total_rew = 0.0
#     steps = 1
#     if !(callback isa Nothing)
#         callback(env, agent, (state=s_t,), agent_ret, steps)
#     end

#     if (total_steps+steps >= max_steps) # 5 Minutes of Gameplay = 18k steps.
#         terminal = true
#     end

#     while !terminal
#         # Figure out action!
#         action = if agent_ret isa NamedTuple
#             agent_ret.action
#         else
#             agent_ret
#         end

#         # Step Environment
#         s_tp1, rew, terminal = if rng isa Nothing
#             step!(env, action)
#         else
#             step!(env, action, rng)
#         end

#         # Step Agent
#         agent_ret = if rng isa Nothing
#             step!(agent, s_tp1, rew, terminal)
#         else
#             step!(agent, s_tp1, rew, terminal, rng)
#         end

#         total_rew += rew
#         steps += 1

#         if !(callback isa Nothing)
#             callback(env, agent, (state=s_tp1, rew=rew, terminal=terminal), agent_ret, steps)
#         end

#         if (total_steps+steps >= max_steps) || (steps >= max_eps_steps) # 5 Minutes of Gameplay = 18k steps.
#             break
#         end
#     end
#     return total_rew, steps
# end


text:
function parse_all(doc::Document, server)
    ps = CSTParser.ParseState(doc._content)
    StaticLint.clear_meta(getcst(doc))
    if endswith(doc._uri, ".jmd")
        doc.cst, ps = parse_jmd(ps, doc._content)
    else
        doc.cst, ps = CSTParser.parse(ps, true)
    end
    if doc.cst.typ === CSTParser.FileH
        doc.cst.val = doc.path
        doc.cst.ref = doc
    end
    ls_diags = Diagnostic[]
    if server.runlinter && doc._runlinter
        scopepass(getroot(doc))
        mark_errors(doc, ls_diags)
    end
    send(JSONRPC.Request{Val{Symbol("textDocument/publishDiagnostics")},PublishDiagnosticsParams}(nothing, PublishDiagnosticsParams(doc._uri, ls_diags)), server)
end

function mark_errors(doc, out = Diagnostic[])
    line_offsets = get_line_offsets(doc)
    errs = get_errors(doc.cst)
    n = length(errs)
    n == 0 && return out
    i = 1
    start = true
    offset = errs[i][1]

    r = Int[0, 0]
    pos = 0
    nlines = length(line_offsets)
    if offset > last(line_offsets)
        line = nlines
    else
        line = 1
        while line < nlines
            while line_offsets[line] <= offset < line_offsets[line + 1]
                ind = line_offsets[line]
                char = 0
                while offset > ind
                    ind = nextind(doc._content, ind)
                    char += 1
                end
                if start
                    r[1] = line
                    r[2] = char
                    offset += errs[i][2].span
                else
                    if errs[i][2].typ === CSTParser.IDENTIFIER
                        push!(out, Diagnostic(Range(r[1] - 1, r[2], line - 1, char), 2, "Julia", "Julia", "Missing reference: $(errs[i][2].val)", nothing))
                    elseif errs[i][2].typ === CSTParser.ErrorToken
                        push!(out, Diagnostic(Range(r[1] - 1, r[2], line - 1, char), 1, "Julia", "Julia", "Parsing error", nothing))
                    elseif errs[i][2].ref == StaticLint.IncorrectCallNargs
                        push!(out, Diagnostic(Range(r[1] - 1, r[2], line - 1, char), 2, "Julia", "Julia", "Incorrect number of arguments", nothing))
                    elseif errs[i][2].ref == StaticLint.IncorrectIterSpec
                        push!(out, Diagnostic(Range(r[1] - 1, r[2], line - 1, char), 2, "Julia", "Julia", "Incorrect specification for iterator", nothing))
                    elseif errs[i][2].ref == StaticLint.NothingEquality
                        push!(out, Diagnostic(Range(r[1] - 1, r[2], line - 1, char), 2, "Julia", "Julia", "Compare against `nothing` using === ", nothing))
                     elseif errs[i][2].ref == StaticLint.NothingNotEq
                        push!(out, Diagnostic(Range(r[1] - 1, r[2], line - 1, char), 2, "Julia", "Julia", "Compare against `nothing` using !== ", nothing))
                    end
                    i += 1
                    i>n && break
                    offset = errs[i][1]
                end
                start = !start
                offset = start ? errs[i][1] : errs[i][1] + errs[i][2].span
            end
            line += 1
        end
    end
    return out
end

function get_errors(x::EXPR, errs = Tuple{Int,EXPR}[], pos = 0)
    if x.typ === CSTParser.ErrorToken || (CSTParser.isidentifier(x) && !(StaticLint.hasref(x)) #= && !(x.parent isa EXPR && x.parent.typ == CSTParser.Quotenode) =#)
        if x.ref != StaticLint.NoReference
            push!(errs, (pos, x))
        end
    elseif x.ref isa StaticLint.LintCodes
        push!(errs, (pos, x))
    end
    if x.args !== nothing
        for i in 1:length(x.args)
            get_errors(x.args[i], errs, pos)
            pos += x.args[i].fullspan
        end
    end
    errs
end

function convert_diagnostic(h::LSDiagnostic{T}, doc::Document) where {T}
    rng = Range(doc, h.loc)
    code = 1
    Diagnostic(rng, code, string(T), string(typeof(h).name), isempty(h.message) ? string(T) : h.message, nothing)
end

function publish_diagnostics(doc::Document, server)
    ls_diags = map(doc.diagnostics) do diag
        convert_diagnostic(diag, doc)
    end
    publishDiagnosticsParams = PublishDiagnosticsParams(doc._uri, ls_diags)
    response =  JSONRPC.Request{Val{Symbol("textDocument/publishDiagnostics")},PublishDiagnosticsParams}(nothing, publishDiagnosticsParams)
    send(response, server)
end

function clear_diagnostics(uri::URI2, server)
    doc = server.documents[uri]
    empty!(doc.diagnostics)
    response =  JSONRPC.Request{Val{Symbol("textDocument/publishDiagnostics")},PublishDiagnosticsParams}(nothing, PublishDiagnosticsParams(doc._uri, Diagnostic[]))
    send(response, server)
end

function clear_diagnostics(server)
    for (uri, doc) in server.documents
        clear_diagnostics(uri, server)
    end
end


text:
using TimerOutputs
using CUDA

macro synctimeit(timer, label, body)
    if has_cuda()
        return :(
            if ($(esc(timer)).enabled)
                @timeit $(esc(timer)) $(esc(label)) CUDA.@sync $(esc(body))
            else
                $(esc(body))
            end
            )
    else
        return :(
            if ($(esc(timer)).enabled)
                @timeit $(esc(timer)) $(esc(label)) $(esc(body))
            else
                $(esc(body))
            end
            )
    end

end

macro withCUDA(expr)
    return has_cuda() ? :($(esc(expr))) : :(nothing)
end


text:
using Test
using GraphletCounting

using LightGraphs

@testset "Count triangles" begin
    expectedDG = SimpleDiGraph(6)
    add_edge!(expectedDG, 2, 1)
    add_edge!(expectedDG, 6, 1)
    add_edge!(expectedDG, 3, 1)
    add_edge!(expectedDG, 3, 4)
    add_edge!(expectedDG, 5, 4)
    sort_by_deg_or!(expectedDG)
    enumerate_wedges(expectedDG)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct ProximityPlacementGroupsApi <: SwaggerApi
    client::Swagger.Client
end

"""

Create or update a proximity placement group.
Param: resourceGroupName::String (required)
Param: proximityPlacementGroupName::String (required)
Param: parameters::ProximityPlacementGroup (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: ProximityPlacementGroup
"""
function _swaggerinternal_proximityPlacementGroupsCreateOrUpdate(_api::ProximityPlacementGroupsApi, resourceGroupName::String, proximityPlacementGroupName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PUT", ProximityPlacementGroup, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/proximityPlacementGroups/{proximityPlacementGroupName}", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "proximityPlacementGroupName", proximityPlacementGroupName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function proximityPlacementGroupsCreateOrUpdate(_api::ProximityPlacementGroupsApi, resourceGroupName::String, proximityPlacementGroupName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_proximityPlacementGroupsCreateOrUpdate(_api, resourceGroupName, proximityPlacementGroupName, parameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function proximityPlacementGroupsCreateOrUpdate(_api::ProximityPlacementGroupsApi, response_stream::Channel, resourceGroupName::String, proximityPlacementGroupName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_proximityPlacementGroupsCreateOrUpdate(_api, resourceGroupName, proximityPlacementGroupName, parameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Delete a proximity placement group.
Param: resourceGroupName::String (required)
Param: proximityPlacementGroupName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Nothing
"""
function _swaggerinternal_proximityPlacementGroupsDelete(_api::ProximityPlacementGroupsApi, resourceGroupName::String, proximityPlacementGroupName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "DELETE", Nothing, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/proximityPlacementGroups/{proximityPlacementGroupName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "proximityPlacementGroupName", proximityPlacementGroupName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function proximityPlacementGroupsDelete(_api::ProximityPlacementGroupsApi, resourceGroupName::String, proximityPlacementGroupName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_proximityPlacementGroupsDelete(_api, resourceGroupName, proximityPlacementGroupName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function proximityPlacementGroupsDelete(_api::ProximityPlacementGroupsApi, response_stream::Channel, resourceGroupName::String, proximityPlacementGroupName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_proximityPlacementGroupsDelete(_api, resourceGroupName, proximityPlacementGroupName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Retrieves information about a proximity placement group .
Param: resourceGroupName::String (required)
Param: proximityPlacementGroupName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Param: includeColocationStatus::String
Return: ProximityPlacementGroup
"""
function _swaggerinternal_proximityPlacementGroupsGet(_api::ProximityPlacementGroupsApi, resourceGroupName::String, proximityPlacementGroupName::String, api_version::String, subscriptionId::String; includeColocationStatus=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", ProximityPlacementGroup, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/proximityPlacementGroups/{proximityPlacementGroupName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "proximityPlacementGroupName", proximityPlacementGroupName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "includeColocationStatus", includeColocationStatus)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function proximityPlacementGroupsGet(_api::ProximityPlacementGroupsApi, resourceGroupName::String, proximityPlacementGroupName::String, api_version::String, subscriptionId::String; includeColocationStatus=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_proximityPlacementGroupsGet(_api, resourceGroupName, proximityPlacementGroupName, api_version, subscriptionId; includeColocationStatus=includeColocationStatus, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function proximityPlacementGroupsGet(_api::ProximityPlacementGroupsApi, response_stream::Channel, resourceGroupName::String, proximityPlacementGroupName::String, api_version::String, subscriptionId::String; includeColocationStatus=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_proximityPlacementGroupsGet(_api, resourceGroupName, proximityPlacementGroupName, api_version, subscriptionId; includeColocationStatus=includeColocationStatus, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Lists all proximity placement groups in a resource group.
Param: resourceGroupName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: ProximityPlacementGroupListResult
"""
function _swaggerinternal_proximityPlacementGroupsListByResourceGroup(_api::ProximityPlacementGroupsApi, resourceGroupName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", ProximityPlacementGroupListResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/proximityPlacementGroups", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function proximityPlacementGroupsListByResourceGroup(_api::ProximityPlacementGroupsApi, resourceGroupName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_proximityPlacementGroupsListByResourceGroup(_api, resourceGroupName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function proximityPlacementGroupsListByResourceGroup(_api::ProximityPlacementGroupsApi, response_stream::Channel, resourceGroupName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_proximityPlacementGroupsListByResourceGroup(_api, resourceGroupName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Lists all proximity placement groups in a subscription.
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: ProximityPlacementGroupListResult
"""
function _swaggerinternal_proximityPlacementGroupsListBySubscription(_api::ProximityPlacementGroupsApi, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", ProximityPlacementGroupListResult, "/subscriptions/{subscriptionId}/providers/Microsoft.Compute/proximityPlacementGroups", ["azure_auth"])
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function proximityPlacementGroupsListBySubscription(_api::ProximityPlacementGroupsApi, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_proximityPlacementGroupsListBySubscription(_api, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function proximityPlacementGroupsListBySubscription(_api::ProximityPlacementGroupsApi, response_stream::Channel, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_proximityPlacementGroupsListBySubscription(_api, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Update a proximity placement group.
Param: resourceGroupName::String (required)
Param: proximityPlacementGroupName::String (required)
Param: parameters::ProximityPlacementGroupUpdate (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: ProximityPlacementGroup
"""
function _swaggerinternal_proximityPlacementGroupsUpdate(_api::ProximityPlacementGroupsApi, resourceGroupName::String, proximityPlacementGroupName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PATCH", ProximityPlacementGroup, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/proximityPlacementGroups/{proximityPlacementGroupName}", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "proximityPlacementGroupName", proximityPlacementGroupName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function proximityPlacementGroupsUpdate(_api::ProximityPlacementGroupsApi, resourceGroupName::String, proximityPlacementGroupName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_proximityPlacementGroupsUpdate(_api, resourceGroupName, proximityPlacementGroupName, parameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function proximityPlacementGroupsUpdate(_api::ProximityPlacementGroupsApi, response_stream::Channel, resourceGroupName::String, proximityPlacementGroupName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_proximityPlacementGroupsUpdate(_api, resourceGroupName, proximityPlacementGroupName, parameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export proximityPlacementGroupsCreateOrUpdate, proximityPlacementGroupsDelete, proximityPlacementGroupsGet, proximityPlacementGroupsListByResourceGroup, proximityPlacementGroupsListBySubscription, proximityPlacementGroupsUpdate


text:
using JuMP
using Cbc
using SCIP
using Clp
d = zeros(6, 4);
sumd = zeros(6, 4);
orderquantity = [500, 800, 1500, 700, 900, 500];
order = [1,2,3,4,5,6]

ordertime = [3, 4, 5, 4, 5, 6]
orderitem = [1, 1, 1, 2, 2, 1]

d[3, 1] = 500
d[4, 1] = 800
d[5, 1] = 1500
d[4, 2] = 700
d[5, 2] = 900
d[6, 1] = 500

mps = Model(Cbc.Optimizer)
items = [1,2,3,4]
date = [1,2,3,4,5,6]
res = [1]
s_init = 0
alpha = [20, 30, 0, 0]
@variable(mps, x[t in date,i in items,] >= 0)
@variable(mps, s[t in date, i in items] >= 0)
@variable(mps, y[t in date, i in items], Bin)

@constraint(mps, dem_sat1[i in items], s_init + x[1, i] == d[1, i] + s[1, i])
@constraint(mps, dem_sat[t in date[2:end], i in items], s[ t - 1, i] + x[t, i] == d[t, i] + s[t, i])
@constraint(mps, vub[t in date, i in items], x[t, i] <= sum(d[k, i]  for k in date[t:end]) * y[t, i])
@constraint(mps, reslimit[t in date], sum(x[t,i] * alpha[i] for i in items) <= 3600 * 8)
@objective(mps, Min, sum(s[t, i] for i in items, t in date))
optimize!(mps)
x_val = value.(x)
x_val[2,1] = 460
x_val[3,1] = 840

sub1 = Model(Clp.Optimizer)
order1 = [1,2,3,6]
@variable(sub1, wx[t in date, p in order1] >= 0)
# @variable(sub1, wy[t in date, p in order1], Bin)
@constraint(sub1,  dem_order[p in order1], sum(wx[t, p] for t in 1:ordertime[p])  == orderquantity[p] )
# @constraint(sub1, xlimit[t in date, p in order1], wx[t, p] <= x_val[t,1] );
@constraint(sub1, xeq[t in date], sum(wx[t, p] for p in order1) == x_val[ t,1]);
@objective(sub1, Min, sum(wx[t, p] * (ordertime[p] - 1) for t in date, p in order1));
optimize!(sub1);
@show value.(wx);

sub2 = Model(Clp.Optimizer)
order2 = [4,5]
@variable(sub2, wx2[t in date, p in order2] >= 0)
@constraint(sub2 , dem_order2[p in order2], sum(wx2[t, p] for t in 1:ordertime[p]) == orderquantity[p])
@constraint(sub2, xeq2[t in date], sum(wx2[t, p] for p in order2) == x_val[t,2])

@objective(sub2, Min, sum(wx2[t, p] * (ordertime[p] - 1) for t in date, p in order2) )
optimize!(sub2)
@show value.(wx2)

sub0 =  Model(Clp.Optimizer)
@variable(sub0, wx0[t in date, p in order] >= 0)
@constraint(sub0, dem_order0[p in order], sum(wx0[t, p] for t in 1:ordertime[p]) == orderquantity[p])
@constraint(sub0, xeq01[t in date], sum(wx0[t, p] for p in order1) == x_val[t, 1])
@constraint(sub0, xeq02[t in date], sum(wx0[t, p] for p in order2) == x_val[t, 2])

@objective(sub0, Min, sum(wx0[t, p] * (ordertime[p] - 1) for t in date, p in order))
optimize!(sub0)
@show value.(wx0)

max_stars_repo_name: Behemoth-s/ProductPlanningByMIP

max_issues_repo_name: Behemoth-s/ProductPlanningByMIP

max_forks_repo_name: Behemoth-s/ProductPlanningByMIP


text:
module Day9
using DataStructures
export parse_input, part_1

function parse_input(input::String)::Tuple{Int, Int}
    m = match(r"(\d+) players; last marble is worth (\d+) points", input)
    return tuple(map(x -> parse(Int, x), m.captures)...)
end

function rotate!(deque::Deque, times::Int)
    if times < 0
        for i in 1:abs(times)
            pushfirst!(deque, pop!(deque))
        end
    else
        for i in 1:times
            push!(deque, popfirst!(deque))
        end
    end
end

function part_1(players, max_points)::Int
    marbles = Deque{Int}()
    push!(marbles, 0)
    points = zeros(Int, players)

    for value in 1:max_points
        if value % 23 == 0
            player = mod1(value, players)
            points[player] += value

            rotate!(marbles, -7)
            points[player] += pop!(marbles)
            rotate!(marbles, 1)
        else
            rotate!(marbles, 1)
            push!(marbles, value)
        end
    end

    return maximum(points)
end

end


text:
@testset "45.jump-game-ii.jl" begin
    @test jump45([2, 3, 1, 1, 4]) == 2
end


text:
using Flux
using ParameterSchedulers
using ParameterSchedulers:AbstractSchedule

struct Warmup <: AbstractSchedule
    eta::Float32
    warmup::Int
end

function (schedule::Warmup)(t::Int)
    schedule.eta * min(t^(-0.5), t * schedule.warmup^(-1.5))
end

Base.getindex(schedule::Warmup, t::Int) = schedule(t)
Base.iterate(schedule::Warmup, t::Int=1) = schedule(t), t + 1

Warmup(;η=0.001, w=10_000) = Warmup(η, w)

struct WarmupADAM
    adam::ADAM
    warmup::ScheduleIterator{Warmup,<:Any}
end

WarmupADAM(η, w, β) = WarmupADAM(
    ADAM(0.0, β),
    ScheduleIterator(Warmup(η, w))
)
WarmupADAM(;η=0.001, w=10_000, β=(0.9, 0.999)) = WarmupADAM(η, w, β)

function Flux.Optimise.apply!(opt::WarmupADAM, x, Δ)
    opt.adam.eta = next!(opt.warmup)
    return Flux.Optimise.apply!(opt.adam, x, Δ)
end


text:
using Jvuetify
using Test

@testset "Jvuetify.jl" begin
    # @test my_f(2,1) == 5
    # @test my_f(2,3) == 7
    @test true == true
end


text:
using BenchmarkTools
using ChainPartitioners
using Cthulhu
using SparseArrays
using JSON

include("matrices.jl")

function main(args)
    suite = BenchmarkGroup()

    suite["partition_stripe"] = BenchmarkGroup(["partition", "stripe"])
    suite["partition_plaid"] = BenchmarkGroup(["partition", "plaid"])
    suite["pack_stripe"] = BenchmarkGroup(["pack", "stripe"])

    work_model = (AffineWorkModel(0, 10, 1), "work_model")
    net_model = (AffineConnectivityModel(0, 10, 1, 100), "net_model")
    sym_model = (AffineMonotonizedSymmetricConnectivityModel(0, 0, 1, 100, 90), "sym_model")
    comm_model = (AffinePrimaryConnectivityModel(0, 10, 1, 0, 100), "comm_model")
    local_model = (AffineSecondaryConnectivityModel(0, 10, 1, 0, 100), "local_model")
    col_block_model = (ColumnBlockComponentCostModel{Int}(3, (w) -> 1 + w), "col_block_model")
    block_model = (BlockComponentCostModel{Int}(1, 3, (1, identity), (1, identity)), "block_model")

    for mtx in [
        "Pajek/GD99_c",
        "LPnetlib/lp_blend",
    ]
        A = SparseMatrixCSC(matrices[mtx])
        (m, n) = size(A)
        ϵ = 0.01
        for (f, f_key) = [
            (DynamicTotalChunker(net_model[1], 8), "DynamicTotalChunker(net_model, 8)"),
            (DynamicTotalChunker(col_block_model[1], 8), "DynamicTotalChunker(col_block_model, 8)"),
            (OverlapChunker(0.9, 8), "OverlapChunker(0.9, 8)"),
            (StrictChunker(8), "StrictChunker(8)"),
        ]
            suite["pack_stripe"]["pack_stripe($(mtx), $f_key)"] = @benchmarkable pack_stripe($A, $f)
        end
        for (Π, Π_key) in [
            (pack_stripe(A', EquiChunker(4)), "Π_spl"),
        ]
            for (f, f_key) = [
                (DynamicTotalChunker(block_model[1], 8), "DynamicTotalChunker(block_model, 8)"),
            ]
                suite["pack_stripe"]["pack_stripe($(mtx), $f_key, $Π_key)"] = @benchmarkable pack_stripe($A, $f, $Π)
            end
        end
        for K = [2, ceil(Int, n^(3/4))]
            mdls = [
                work_model,
                net_model,
            ]
            if issymmetric(pattern(A))
                append!(mdls, [
                    sym_model
                ])
            end
            for (f, f_key) = [
                (EquiSplitter(), "EquiSplitter()"),
                ((BisectIndexBottleneckSplitter(mdl), "BisectIndexBottleneckSplitter($mdl_key)") for (mdl, mdl_key) in mdls)...,
                ((BisectCostBottleneckSplitter(mdl, ϵ), "BisectCostBottleneckSplitter($mdl_key)") for (mdl, mdl_key) in mdls)...,
                ((LazyBisectCostBottleneckSplitter(mdl, ϵ), "LazyBisectCostBottleneckSplitter($mdl_key)") for (mdl, mdl_key) in mdls[2:end])...,
                (DynamicTotalSplitter(net_model[1]), "DynamicTotalSplitter(net_model)"),
            ]
                suite["partition_stripe"]["partition_stripe($(mtx), $K, $f_key)"] = @benchmarkable partition_stripe($A, $K, $f)
            end
            for (Π, Π_key, mdls) in [
                (partition_stripe(A', K, EquiSplitter()), "Π_spl", [local_model, comm_model]),
                (MapPartition(K, mod1.(1:m, K)), "Π_map", [comm_model,]),
            ]
                for (f, f_key) = [
                    ((BisectIndexBottleneckSplitter(mdl), "BisectIndexBottleneckSplitter($mdl_key)") for (mdl, mdl_key) in mdls)...,
                    ((BisectCostBottleneckSplitter(mdl, ϵ), "BisectCostBottleneckSplitter($mdl_key)") for (mdl, mdl_key) in mdls)...,
                    (LazyBisectCostBottleneckSplitter(comm_model[1], ϵ), "LazyBisectCostBottleneckSplitter(comm_model)"),
                    (GreedyBottleneckPartitioner(local_model[1]), "GreedyBottleneckPartitioner(local_model)"),
                    (MagneticPartitioner(), "MagneticPartitioner"),
                ]
                    suite["partition_stripe"]["partition_stripe($(mtx), $K, $f_key, $Π_key)"] = @benchmarkable partition_stripe($A, $K, $f, $Π)
                end
            end
        end
    end

    if length(args) == 1 && args[1] == "update"
        tune!(suite, verbose=true, seconds=60)
        BenchmarkTools.save(joinpath(@__DIR__, "params.json"), BenchmarkTools.params(suite))
        results = run(suite, verbose=true, seconds=60)
        BenchmarkTools.save(joinpath(@__DIR__, "benchmarks.json"), results)
        return
    elseif length(args) == 0
        params = BenchmarkTools.load(joinpath(@__DIR__, "params.json"))[1]
        BenchmarkTools.loadparams!(suite, params)
        results = run(suite, seconds=60, verbose=true)
        reference = BenchmarkTools.load(joinpath(@__DIR__, "benchmarks.json"))[1]
        judgements = judge(minimum(results), minimum(reference); time_tolerance = 0.1, memory_tolerance = 0.1)
        show(IOContext(stdout, :compact=>false, :limit=>false), regressions(judgements))
        show(IOContext(stdout, :compact=>false, :limit=>false), improvements(judgements))
        println()
        return
    elseif length(args) == 1
        key = args[1]
        params = BenchmarkTools.load(joinpath(@__DIR__, "params.json"))[1]
        BenchmarkTools.loadparams!(suite, params)
        results = run(suite[key], seconds=60, verbose=true)
        reference = BenchmarkTools.load(joinpath(@__DIR__, "benchmarks.json"))[1][key]
        judgements = judge(minimum(results), minimum(reference); time_tolerance = 0.1, memory_tolerance = 0.1)
        show(IOContext(stdout, :compact=>false, :limit=>false), regressions(judgements))
        #show(IOContext(stdout, :compact=>false, :limit=>false), improvements(judgements))
        println()
        return
    else
        throw(ArgumentError("unrecognized command line arguments"))
    end
end

main(ARGS)