text:
"""
    HubbardModelAttractive(lattice; params...)
    HubbardModelAttractive(L, dims; params...)
    HubbardModelAttractive(params::Dict)
    HubbardModelAttractive(params::NamedTuple)
    HubbardModelAttractive(; params...)

Defines an attractive (negative `U`) Hubbard model on a given (or derived)
`lattice`. If a linear system size `L` and dimensionality `dims` is given, the
`lattice` will be a Cubic lattice of fitting size.

Additional parameters (keyword arguments) include:
* `l::AbstractLattice = lattice`: The lattice the model uses. The keyword
argument takes precedence over the argument `lattice`.
* `U::Float64 = 1.0 > 0.0` is the absolute value of the Hubbard Interaction.
* `t::Float64 = 1.0` is the hopping strength.
* `mu::Float64` is the chemical potential.

Internally, a discrete Hubbard Stratonovich transformation (Hirsch
transformation) is used in the spin/magnetic channel to enable DQMC. The
resulting Hubbard Stratonovich fiels is real.
Furthermore, we use spin up/down symmetry to speed up the simulation. As a
result the greens matrix is of size (N, N) with N the number of sites, and the
element G[i, j] corresponds to the up-up and down-down element.
"""
@with_kw_noshow struct HubbardModelAttractive{LT<:AbstractLattice} <: HubbardModel
    # user optional
    mu::Float64 = 0.0
    U::Float64 = 1.0
    @assert U >= 0. "U must be positive."
    t::Float64 = 1.0
    l::LT

    # non-user fields
    flv::Int = 1
    # to avoid allocations (TODO always real?)
    IG::Vector{Float64} = Vector{Float64}(undef, length(l))
    G::Vector{Float64} = Vector{Float64}(undef, length(l))
end

HubbardModelAttractive(params::Dict{Symbol}) = HubbardModelAttractive(; params...)
HubbardModelAttractive(params::NamedTuple) = HubbardModelAttractive(; params...)
function HubbardModelAttractive(lattice::AbstractLattice; kwargs...)
    HubbardModelAttractive(l = lattice; kwargs...)
end
function HubbardModelAttractive(L, dims; kwargs...)
    l = choose_lattice(HubbardModelAttractive, dims, L)
    HubbardModelAttractive(l = l; kwargs...)
end

# cosmetics
import Base.summary
import Base.show
Base.summary(model::HubbardModelAttractive) = "attractive Hubbard model"
function Base.show(io::IO, model::HubbardModelAttractive)
    print(io, "attractive Hubbard model, $(length(model.l)) sites")
end
Base.show(io::IO, m::MIME"text/plain", model::HubbardModelAttractive) = print(io, model)

# Convenience
@inline parameters(m::HubbardModelAttractive) = (N = length(m.l), t = m.t, U = -m.U, mu = m.mu)

"""
Calculates the hopping matrix \$T_{i, j}\$ where \$i, j\$ are
site indices.

Note that since we have a time reversal symmetry relating spin-up
to spin-down we only consider one spin sector (one flavor) for the attractive
Hubbard model in the DQMC simulation.

This isn't a performance critical method as it is only used once before the
actual simulation.
"""
function hopping_matrix(mc::DQMC, m::HubbardModelAttractive{L}) where {L<:AbstractLattice}
    N = length(m.l)
    T = diagm(0 => fill(-m.mu, N))

    # Nearest neighbor hoppings
    @inbounds @views begin
        for (src, trg) in neighbors(m.l, Val(true))
            trg == -1 && continue
            T[trg, src] += -m.t
        end
    end

    return T
end

"""
Calculate the interaction matrix exponential `expV = exp(- power * delta_tau * V(slice))`
and store it in `result::Matrix`.

This is a performance critical method.
"""
@inline @bm function interaction_matrix_exp!(mc::DQMC, m::HubbardModelAttractive,
            result::Diagonal, conf::HubbardConf, slice::Int, power::Float64=1.)
    dtau = mc.parameters.delta_tau
    lambda = acosh(exp(0.5 * m.U * dtau))

    N = size(result, 1)
    @inbounds for i in 1:N
        result.diag[i] = exp(sign(power) * lambda * conf[i, slice])
    end
    nothing
end

@inline @bm function propose_local(
        mc::DQMC, m::HubbardModelAttractive, i::Int, slice::Int, conf::HubbardConf
    )
    # see for example dos Santos Introduction to quantum Monte-Carlo
    greens = mc.stack.greens
    dtau = mc.parameters.delta_tau
    lambda = acosh(exp(m.U * dtau/2))

    @inbounds ΔE_boson = -2. * lambda * conf[i, slice]
    γ = exp(ΔE_boson) - 1
    @inbounds detratio = (1 + γ * (1 - greens[i,i]))^2
    # squared because of two spin sectors

    return detratio, ΔE_boson, γ
end

@inline @bm function accept_local!(
        mc::DQMC, m::HubbardModelAttractive, i::Int, slice::Int, conf::HubbardConf,
        detratio, ΔE_boson, γ)
    greens = mc.stack.greens

    # Unoptimized Version
    # u = -greens[:, i]
    # u[i] += 1.
    # greens .-= kron(u * 1./(1 + gamma * u[i]), transpose(gamma * greens[i, :]))
    # conf[i, slice] .*= -1.

    # Optimized
    # copy! and `.=` allocate, this doesn't. Synced loop is marginally faster
    @turbo for j in eachindex(m.IG)
        m.IG[j] = -greens[j, i]
        m.G[j] = greens[i, j]
    end
    @inbounds m.IG[i] += 1.0
    # This is way faster for small systems and still ~33% faster at L = 15
    # Also no allocations here
    @inbounds x = γ / (1.0 + γ * m.IG[i])
    @turbo for k in eachindex(m.IG), l in eachindex(m.G)
        greens[k, l] -= m.IG[k] * x * m.G[l]
    end
    @inbounds conf[i, slice] *= -1
    nothing
end

@bm function propose_global_from_conf(mc::DQMC, m::HubbardModelAttractive, conf::AbstractArray)
    # G = calculate_greens(mc, current_slice(mc)-1, mc.s.greens_temp)
    # @assert G ≈ mc.s.greens
    # D = copy(mc.s.Dl)
    # new_greens = calculate_greens(mc, current_slice(mc)-1, mc.s.greens_temp, conf)

    # # Computing prod(D) is unstable at larger beta
    # detratio = 1.0
    # for i in eachindex(mc.s.Dl)
    #     detratio *= (D[i] / mc.s.Dl[i])^2
    # end
    # ΔE_Boson = energy_boson(mc, m, conf) - energy_boson(mc, m)
    # return detratio, ΔE_Boson, new_greens

    # I don't think we need this...
    @assert mc.stack.current_slice == 1
    @assert mc.stack.direction == 1

    # This should be just after calculating greens, so mc.s.Dl is from the UDT
    # decomposed G
    copyto!(mc.stack.tempvf, mc.stack.Dl)

    # -1?
    inv_det(mc, current_slice(mc)-1, conf)

    # This may help with stability
    detratio = 1.0
    for i in eachindex(mc.stack.tempvf)
        detratio *= mc.stack.tempvf[i] * mc.stack.Dr[i]
        # detratio *= mc.s.Dl[i] / D[i]
    end
    ΔE_Boson = energy_boson(mc, m, conf) - energy_boson(mc, m)
    # new_greens = finish_calculate_greens(
    #     mc.s.Ul, mc.s.Dl, mc.s.Tl, mc.s.Ur, mc.s.Dr, mc.s.Tr,
    #     mc.s.greens_temp, mc.s.pivot, mc.s.tempv
    # )
    # @info detratio^2
    return detratio^2, ΔE_Boson, nothing #new_greens
end

"""
Calculate energy contribution of the boson, i.e. Hubbard-Stratonovich/Hirsch field.
"""
@inline function energy_boson(mc::DQMC, m::HubbardModelAttractive, hsfield = conf(mc))
    dtau = mc.parameters.delta_tau
    lambda = acosh(exp(m.U * dtau/2))
    return lambda * sum(hsfield)
end

function greens(mc::DQMC, model::HubbardModelAttractive)
    G = greens!(mc)
    vcat(hcat(G, zeros(size(G))), hcat(zeros(size(G)), G))
end

function save_model(
        file::JLDFile,
        m::HubbardModelAttractive,
        entryname::String="Model"
    )
    write(file, entryname * "/VERSION", 1)
    write(file, entryname * "/tag", "HubbardModelAttractive")

    write(file, entryname * "/mu", m.mu)
    write(file, entryname * "/U", m.U)
    write(file, entryname * "/t", m.t)
    save_lattice(file, m.l, entryname * "/l")
    write(file, entryname * "/flv", m.flv)

    nothing
end

#     load_parameters(data, ::Type{<: DQMCParameters})
#
# Loads a DQMCParameters object from a given `data` dictionary produced by
# `JLD.load(filename)`.
function _load(data, ::Val{:HubbardModelAttractive})
    if !(data["VERSION"] == 1)
        throw(ErrorException("Failed to load HubbardModelAttractive version $(data["VERSION"])"))
    end

    l = _load(data["l"], to_tag(data["l"]))
    HubbardModelAttractive(
        mu = data["mu"],
        U = data["U"],
        t = data["t"],
        l = l,
        flv = data["flv"]
    )
end
to_tag(::Type{<: HubbardModelAttractive}) = Val(:HubbardModelAttractive)

################################################################################
### Measurement kernels
################################################################################
checkflavors(model::HubbardModelAttractive) = checkflavors(model, 1)

function cdc_kernel(mc, ::HubbardModelAttractive, ij::NTuple{2}, G::GreensMatrix)
    # spin up and down symmetric, so (i+N, i+N) = (i, i); (i+N, i) drops
    i, j = ij
    4 * dagger(G)[i, i] * dagger(G)[j, j] + 2 * dagger(G)[i, j] * G[i, j]
end
function cdc_kernel(mc, ::HubbardModelAttractive, ij::NTuple{2}, pg::NTuple{4})
    i, j = ij
    G00, G0l, Gl0, Gll = pg
    # spin up and down symmetric, so (i+N, i+N) = (i, i); (i+N, i) drops
    4 * dagger(Gll)[i, i] * dagger(G00)[j, j] + 2 * dagger(G0l)[i, j] * Gl0[i, j]
end

mx_kernel(mc, ::HubbardModelAttractive, i, G::GreensMatrix) = 0.0
my_kernel(mc, ::HubbardModelAttractive, i, G::GreensMatrix) = 0.0
mz_kernel(mc, ::HubbardModelAttractive, i, G::GreensMatrix) = 0.0

function sdc_x_kernel(mc, ::HubbardModelAttractive, ij::NTuple{2}, G::GreensMatrix)
    i, j = ij
    2 * dagger(G)[i, j] * G[i, j]
end
function sdc_y_kernel(mc, ::HubbardModelAttractive, ij::NTuple{2}, G::GreensMatrix)
    i, j = ij
    2 * dagger(G)[i, j] * G[i, j]
end
function sdc_z_kernel(mc, ::HubbardModelAttractive, ij::NTuple{2}, G::GreensMatrix)
    i, j = ij
    2 * dagger(G)[i, j] * G[i, j]
end

function sdc_x_kernel(mc, ::HubbardModelAttractive, ij::NTuple{2}, pg::NTuple{4})
    i, j = ij
    2 * dagger(pg[2])[i, j] * pg[3][i, j]
end
function sdc_y_kernel(mc, ::HubbardModelAttractive, ij::NTuple{2}, pg::NTuple{4})
    i, j = ij
    2 * dagger(pg[2])[i, j] * pg[3][i, j]
end
function sdc_z_kernel(mc, ::HubbardModelAttractive, ij::NTuple{2}, pg::NTuple{4})
    i, j = ij
    2 * dagger(pg[2])[i, j] * pg[3][i, j]
end

function pc_kernel(mc, ::HubbardModelAttractive, sites::NTuple{4}, G::GreensMatrix)
    src1, trg1, src2, trg2 = sites
    G[src1, src2] * G[trg1, trg2]
end
function pc_kernel(mc, ::HubbardModelAttractive, sites::NTuple{4}, pg::NTuple{4})
    src1, trg1, src2, trg2 = sites
    pg[3][src1, src2] * pg[3][trg1, trg2]
end
function pc_alt_kernel(mc, ::HubbardModelAttractive, sites::NTuple{4}, packed_greens::NTuple{4})
    src1, trg1, src2, trg2 = sites
	G00, G0l, Gl0, Gll = packed_greens
    dagger(G0l)[trg1, trg2] * dagger(G0l)[src1, src2]
end
function pc_ref_kernel(mc, ::HubbardModelAttractive, sites::NTuple{4}, packed_greens::NTuple{4})
    src1, trg1, src2, trg2 = sites
	G00, G0l, Gl0, Gll = packed_greens
    Gl0[src1, src2] * Gl0[trg1, trg2] +
    (I[trg2, trg1] - G0l[trg2, trg1]) * (I[src2, src1] - G0l[src2, src1])
end

function cc_kernel(mc, ::HubbardModelAttractive, sites::NTuple{4}, pg::NTuple{4})
    src1, trg1, src2, trg2 = sites
    G00, G0l, Gl0, Gll = pg
    T = mc.stack.hopping_matrix

    # up-up counts, down-down counts, mixed only on 11s or 22s
    s1 = src1; t1 = trg1
    s2 = src2; t2 = trg2
    output = -(
        4.0 * (
            T[t2, s2] * (I[s2, t2] - Gll[s2, t2]) -
            T[s2, t2] * (I[t2, s2] - Gll[t2, s2])
        ) * (
            T[t1, s1] * (I[s1, t1] - G00[s1, t1]) -
            T[s1, t1] * (I[t1, s1] - G00[t1, s1])
        ) +
        - 2.0 * T[t2, s2] * T[t1, s1] * G0l[s1, t2] * Gl0[s2, t1] +
        + 2.0 * T[t2, s2] * T[s1, t1] * G0l[t1, t2] * Gl0[s2, s1] +
        + 2.0 * T[s2, t2] * T[t1, s1] * G0l[s1, s2] * Gl0[t2, t1] +
        - 2.0 * T[s2, t2] * T[s1, t1] * G0l[t1, s2] * Gl0[t2, s1]
    )
    # output = 4.0 *
    #     (T[s1, t1] * Gll[t1, s1] - T[t1, s1] * Gll[s1, t1]) *
    #     (T[s2, t2] * G00[t2, s2] - T[t2, s2] * G00[s2, t2]) +
    #     + 2.0 * T[t1, s1] * T[t2, s2] * (- G0l[s2, t1]) * Gl0[s1, t2] +
    #     - 2.0 * T[s1, t1] * T[t2, s2] * (- G0l[s2, s1]) * Gl0[t1, t2] +
    #     - 2.0 * T[t1, s1] * T[s2, t2] * (- G0l[t2, t1]) * Gl0[s1, s2] +
    #     + 2.0 * T[s1, t1] * T[s2, t2] * (- G0l[t2, s1]) * Gl0[t1, s2] +

    output
end

function intE_kernel(mc, model::HubbardModelAttractive, G::GreensMatrix)
    # ⟨U (n↑ - 1/2)(n↓ - 1/2)⟩ = ...
    # = U [G↑↑ G↓↓ - G↓↑ G↑↓ - 0.5 G↑↑ - 0.5 G↓↓ + G↑↓ + 0.25]
    # = U [(G↑↑ - 1/2)(G↓↓ - 1/2) + G↑↓(1 + G↑↓)]
    # with up-up = down-down and up-down = 0
    - model.U * sum((diag(G.val) .- 0.5).^2)
end


text:
module NurseAllocationFigures

using DataFrames
using CSV
using OrderedCollections: OrderedDict

using Distributions
using Dates

using Memoize
using ProgressMeter

using Gadfly
using Compose
import Cairo, Fontconfig

include("FiguresBase.jl")
using .FiguresBase

function compute_results_nurses(config, data, sent; use_rounding=false)
	N, T = data.N, data.T
	sent_null = zeros(Int, N, N, T)

	if use_rounding
		sent = round.(Int, sent)
	end

	model_title = titlecase(replace(config.experiment, "_" => " ")) * " Model"
	model_title = model_title == "Null Model" ? "No Transfers" : model_title

	active_nurses = compute_active_nurses(sent, data.nurses, use_rounding=use_rounding)
	active_nurses_nosent = compute_active_nurses(sent_null, data.nurses, use_rounding=use_rounding)

	shortage = compute_nurse_shortage(sent, data.nurses, data.demand, use_rounding=use_rounding)
	shortage_nosent = compute_nurse_shortage(sent_null, data.nurses, data.demand, use_rounding=use_rounding)

	load = compute_nurse_load(sent, data.nurses, data.demand, use_rounding=use_rounding)
	load_nosent = compute_nurse_load(sent_null, data.nurses, data.demand, use_rounding=use_rounding)

	results_table = DataFrame(
		date = permutedims(repeat(data.date_range, 1, N), (2,1))[:],
		location = repeat(data.locations, 1, T)[:],
		active_sent = active_nurses[:],
		active_nosent = active_nurses_nosent[:],
		shortage_sent = shortage[:],
		shortage_nosent = shortage_nosent[:],
		load_sent = load[:],
		load_nosent = load_nosent[:],
	)

	results = (
		sent = sent,
		active = active_nurses,
		active_nosent = active_nurses_nosent,
		shortage = shortage,
		shortage_nosent = shortage_nosent,
		load = load,
		load_nosent = load_nosent,
		table = results_table,
		model_title = model_title,
	)
	return results
end;

function compute_active_nurses(_sent, _nurses; use_rounding=false)
	N, _, T = size(_sent)
	_active = [(
		_nurses[i]
		- sum(_sent[i,:,1:t])
		+ sum(_sent[:,i,1:t])
	) for i in 1:N, t in 1:T]
	if use_rounding
		_active = round.(Int, _active)
	end
	return _active
end;

function compute_nurse_shortage(_sent, _nurses, _demand; use_rounding=false)
	return max.(0, compute_active_nurses(_sent, _nurses, use_rounding=use_rounding) .- _demand)
end;

function compute_nurse_load(_sent, _nurses, _demand; use_rounding=false)
	return compute_active_nurses(_sent, _nurses, use_rounding=use_rounding) ./ _demand
end;

function plot_metrics(config, data, results; bedtype=:all, display=true, save=true)
	N, T = data.N, data.T

	sent_out = sum(results.sent, dims=2)[:,1,:]
	sent_in = sum(results.sent, dims=1)[1,:,:]
	tfr = sent_out + sent_in

	tot_tfr = zeros(Int, data.N, data.N)
	for i in 1:data.N, j in (i+1):data.N
		q = cumsum(results.sent[i,j,:] - results.sent[j,i,:])
		z = maximum(abs.(q))
		tot_tfr[i,j] = z
	end

	s = sum(results.shortage)
	s_null = sum(results.shortage_nosent)
	s_red = (s_null - s) / s_null

	total_nurse_days = sum(results.active_nosent)

	if display
		println("Region: $(config.region)")
		println("Level: $(config.alloc_level)")
		println("Period: $(data.start_date) to $(data.end_date)")
		println("Bed type: $(bedtype)")
		println()
		printmetric("Number of days", T, int=true)
		printmetric("Number of hospitals", N, int=true)
		printmetric("Number of hospital-days", N*T, int=true)
		println()
		printmetric("Total nurses", sum(data.nurses))
		printmetric("Total nurse-days", sum(results.active_nosent))
		println()
		printmetric("Total sent", sum(results.sent))
		printmetric("Max sent by hospital-day", maximum(results.sent))
		printmetric("Max sent by day", maximum(sum(results.sent, dims=[1,2])))
		printmetric("Max sent by hospital", maximum(sum(results.sent, dims=[2,3])))
		println()
		printmetric("Number of transfers", sum(results.sent .> 0), int=true)
		printmetric("Percent of hospital-days with a transfer", mean(tfr .> 0), pct=true)
		printmetric("Percent of days with a transfer", mean(sum(tfr, dims=1) .> 0), pct=true)
		printmetric("Percent of hospitals with a transfer", mean(sum(tfr, dims=2) .> 0), pct=true)
		println()
		printmetric("Nurse shortage (without transfers)", s_null)
		printmetric("Nurse shortage (with transfers)", s)
		printmetric("Nurse shortage reduction", s_red, pct=true)
		println()
		printmetric("Percent shortage (without transfers)", s_null / total_nurse_days, pct=true)
		printmetric("Percent shortage (with transfers)", s / total_nurse_days, pct=true)
		printmetric("Shortage reduction", s_red, pct=true)
		println()
		printmetric("Number of hospital-days with a nurse shortage (without transfers)", sum(results.shortage_nosent .> 0))
		printmetric("Number of hospital-days with a nurse shortage (with transfers)", sum(results.shortage .> 0))
		printmetric("Percent of hospital-days with a nurse shortage (without transfers)", mean(results.shortage_nosent .> 0), pct=true)
		printmetric("Percent of hospital-days with a nurse shortage (with transfers)", mean(results.shortage .> 0), pct=true)
		println()
		printmetric("Maximum nurse shortage by hospital-day (without transfers)", maximum(results.shortage_nosent))
		printmetric("Maximum nurse shortage by hospital-day (with transfers)", maximum(results.shortage))
		println()
		printmetric("Maximum nurse load by hospital-day (without transfers)", maximum(skipbad(results.load_nosent)))
		printmetric("Maximum nurse load by hospital-day (with transfers)", maximum(skipbad(results.load)))
	end

	metrics = OrderedDict(
		"experiment date" => config.rundate,
		"experiment" => config.experiment,
		"region" => config.region,
		"allocation level" => config.alloc_level,
		"resource" => "nurses",
		"bed type" => bedtype,
		"start date" => data.start_date,
		"end date" => data.end_date,
		"number of days" => T,
		"number of hospitals" => N,
		"number of hospital-days" => N*T,
		"total nurses" => sum(data.nurses),
		"total nurse-days" => sum(results.active_nosent),
		"total nurses transferred" => sum(results.sent),
		"max sent by hospital-day" => maximum(results.sent),
		"max sent by day" => maximum(sum(results.sent, dims=[1,2])),
		"max sent by hospital" => maximum(sum(results.sent, dims=[2,3])),
		"number of transfers" => sum(results.sent .> 0),
		"proportion of hospital-days with a transfer" => mean(tfr .> 0),
		"proportion of days with a transfer" => mean(sum(tfr, dims=1) .> 0),
		"proportion of hospitals with a transfer" => mean(sum(tfr, dims=2) .> 0),
		"percent shortage (without transfers)" => (s_null / total_nurse_days),
		"percent shortage (with transfers)" => (s / total_nurse_days),
		"nurse shortage (without transfers)" => s_null,
		"nurse shortage (with transfers)" => s,
		"nurse shortage reduction" => s_red,
		"number of hospital-days with a nurse shortage (without transfers)" => sum(results.shortage_nosent .> 0),
		"number of hospital-days with a nurse shortage (with transfers)" => sum(results.shortage .> 0),
		"proportion of hospital-days with a nurse shortage (without transfers)" => mean(results.shortage_nosent .> 0),
		"proportion of hospital-days with a nurse shortage (with transfers)" => mean(results.shortage .> 0),
		"maximum nurse shortage by hospital-day (without transfers)" => maximum(results.shortage_nosent),
		"maximum nurse shortage by hospital-day (with transfers)" => maximum(results.shortage),
		"maximum nurse load by hospital-day (without transfers)" => maximum(skipbad(results.load_nosent)),
		"maximum nurse load by hospital-day (with transfers)" => maximum(skipbad(results.load)),
		"median non-zero transfer" => iszero(results.sent) ? 0 : median(nonzero(results.sent[:])),
		"mean non-zero transfer" => iszero(results.sent) ? 0 : mean(nonzero(results.sent[:])),
		"max non-zero transfer" => iszero(results.sent) ? 0 : maximum(nonzero(results.sent[:])),
	)
	metrics_df = DataFrame(metric=collect(keys(metrics)), value=collect(values(metrics)))
	metrics_str = metrics_totextable(metrics_df)

	config_metrics = OrderedDict(
		"experiment date" => config.rundate,
		"experiment" => config.experiment,
		"region" => config.region,
		"allocation level" => config.alloc_level,
		"resource" => "nurses",
		"bed type" => bedtype,
		"start date" => data.start_date,
		"end date" => data.end_date,
		"number of days" => T,
		"number of hospitals" => N,
		"number of hospital-days" => N*T,
	)
	config_metrics_df = DataFrame(metric=collect(keys(config_metrics)), value=collect(values(config_metrics)))
	config_metrics_str = metrics_totextable(config_metrics_df)

	roundmetric(x) = round(x, digits=1)
	to_percent(x) = string(roundmetric(x*100))*"%"
	metrics_fortable = OrderedDict(
		"experiment" => config.experiment,

		# result metrics

		"shortage" => roundmetric(s),
		"shortage reduction" => to_percent(s_red),

		"median non-zero shortage" => iszero(results.shortage) ? 0 : roundmetric(median(nonzero(results.shortage[:]))),
		"mean non-zero shortage" => iszero(results.shortage) ? 0 : roundmetric(mean(nonzero(results.shortage[:]))),
		"max non-zero shortage" => iszero(results.shortage) ? 0 : roundmetric(maximum(nonzero(results.shortage[:]))),

		"median load" => to_percent(median(skipbad(results.load[:]))),
		"mean load" => to_percent(mean(skipbad(results.load[:]))),
		"max load" => to_percent(maximum(skipbad(results.load[:]))),

		"number of hospital-days with a shortage" => roundmetric(sum(results.shortage .> 0)),
		"percent of hospital-days with a shortage" => to_percent(mean(results.shortage .> 0)),

		# transfer metrics

		"total nurse transfers" => roundmetric(sum(results.sent)),
		"percent of nurses transferred" => to_percent(sum(tot_tfr) / sum(data.nurses)),

		"median non-zero transfer" => iszero(results.sent) ? 0 : roundmetric(median(nonzero(results.sent[:]))),
		"mean non-zero transfer" => iszero(results.sent) ? 0 : roundmetric(mean(nonzero(results.sent[:]))),
		"max non-zero transfer" => iszero(results.sent) ? 0 : roundmetric(maximum(nonzero(results.sent[:]))),

		"percent of hospital-days with a transfer" => to_percent(mean(tfr .> 0)),
		"percent of days with a transfer" => to_percent(mean(sum(tfr, dims=1) .> 0)),
		"percent of hospitals with a transfer" => to_percent(mean(sum(tfr, dims=2) .> 0)),
	)
	metrics_table_df = DataFrame(metric=collect(keys(metrics_fortable)), value=collect(values(metrics_fortable)))
	metrics_table_str = metrics_totextable(metrics_table_df)

	if save
		bedtype_ext = (bedtype == :all) ? "" : (bedtype == :icu ? "_icu" : (bedtype == :ward ? "_ward" : ""))

		metrics_dir = "figures/$(config.region_abbrev)/$(config.rundate)/$(config.experiment)/metrics/"
		if !isdir(metrics_dir) mkpath(metrics_dir) end

		metrics_path_base = joinpath(metrics_dir, "metrics_$(config.experiment)$(bedtype_ext)")
		metrics_df |> CSV.write(metrics_path_base * ".csv")
		write(metrics_path_base * ".tex", metrics_str)

		metrics_table_path_base = joinpath(metrics_dir, "metrics_table_$(config.experiment)$(bedtype_ext)")
		metrics_table_df |> CSV.write(metrics_table_path_base * ".csv")
		write(metrics_table_path_base * ".tex", metrics_table_str)

		config_metrics_path_base = joinpath(metrics_dir, "config_$(config.experiment)$(bedtype_ext)")
		config_metrics_df |> CSV.write(config_metrics_path_base * ".csv")
		write(config_metrics_path_base * ".tex", config_metrics_str)
	end

	return
end;

function plot_shortage_distribution(config, data, results; bedtype=:all, display=true, save=true)
	colors = Scale.default_discrete_colors(4)[[1,3]]
	binwidth = 5
	shortage_nonzero = filter(x -> x > 0, results.shortage)
	shortage_nosent_nonzero = filter(x -> x > 0, results.shortage_nosent)
	nbins_sent = !isempty(shortage_nonzero) ? ceil(maximum(shortage_nonzero) / binwidth) : 1
	nbins_nosent = !isempty(shortage_nosent_nonzero) ? ceil(maximum(shortage_nosent_nonzero) / binwidth) : 1
	shortage_dist_plot = plot(
		layer(
			x=shortage_nonzero,
			Geom.histogram(
				bincount=nbins_sent,
				limits=(min=0, max=nbins_sent*binwidth)
			),
			style(default_color=colors[1], alphas=[0.7])
		),
		layer(
			x=shortage_nosent_nonzero,
			Geom.histogram(
				bincount=nbins_nosent,
				limits=(min=0, max=nbins_nosent*binwidth)
			),
			style(default_color=colors[2], alphas=[0.7])
		),
		Guide.title("Distribution of (Non-Zero) Shortage by Hospital-Day"),
		Guide.xlabel("Shortage"),
		Guide.ylabel("Count"),
		Guide.manual_color_key("", ["With Transfers", "Without Transfers"], colors, shape=[Shape.square], size=[1.6mm]),
		style(
			minor_label_font_size=18px,
			major_label_font_size=22px,
			minor_label_font="CMU Serif",
			major_label_font="CMU Serif",
			key_label_font_size=16px,
			key_title_font_size=18px,
			key_label_font="CMU Serif",
			key_title_font="CMU Serif",
			background_color=colorant"white",
		),
	)

	bedtype_ext = (bedtype == :all) ? "" : (bedtype == :icu ? "_icu" : (bedtype == :ward ? "_ward" : ""))

	if display
		shortage_dist_plot |> SVG(23cm, 10cm)
	end
	if save
		save_all(shortage_dist_plot, 23cm, 10cm, "shortage_dist"*bedtype_ext, config)
	end
end;

function plot_load_distribution(config, data, results; bedtype=:all, display=true, save=true)
	colors = Scale.default_discrete_colors(4)[[1,3]]
	binwidth = 0.2
	load_sent = skipbad(results.load)
	load_nosent = skipbad(results.load_nosent)
	nbins_sent = ceil(maximum(load_sent) / binwidth)
	nbins_nosent = ceil(maximum(load_nosent) / binwidth)
	load_dist_plot = plot(
		layer(
			x=load_sent,
			Geom.histogram(
				bincount=nbins_sent,
				limits=(min=0, max=nbins_sent*binwidth)
			),
			style(default_color=colors[1], alphas=[0.7])
		),
		layer(
			x=load_nosent,
			Geom.histogram(
				bincount=nbins_nosent,
				limits=(min=0, max=nbins_nosent*binwidth)
			),
			style(default_color=colors[2], alphas=[0.7])
		),
		layer(xintercept=[1.0], Geom.vline(color="red", size=1.0mm), order=10),
		Guide.title("Distribution of Load by Hospital-Day"),
		Guide.xlabel("Load"),
		Guide.ylabel("Count"),
		Guide.manual_color_key("", ["With Transfers", "Without Transfers"], colors, shape=[Shape.square], size=[1.6mm]),
		Coord.Cartesian(xmax=5.0),
		style(
			minor_label_font_size=18px,
			major_label_font_size=22px,
			minor_label_font="CMU Serif",
			major_label_font="CMU Serif",
			key_label_font_size=16px,
			key_title_font_size=18px,
			key_label_font="CMU Serif",
			key_title_font="CMU Serif",
			background_color=colorant"white",
		),
	)

	bedtype_ext = (bedtype == :all) ? "" : (bedtype == :icu ? "_icu" : (bedtype == :ward ? "_ward" : ""))

	if display
		load_dist_plot |> SVG(23cm, 10cm)
	end
	if save
		save_all(load_dist_plot, 23cm, 10cm, "load_dist"*bedtype_ext, config);
	end
end;

function plot_maxload_distribution(config, data, results; bedtype=:all, display=true, save=true)
	colors = Scale.default_discrete_colors(4)[[1,3]]
	binwidth = 0.2
	maxload_sent = skipbad(maximum(results.load, dims=2))
	maxload_nosent = skipbad(maximum(results.load_nosent, dims=2))
	nbins_sent = ceil(maximum(maxload_sent) / binwidth)
	nbins_nosent = ceil(maximum(maxload_nosent) / binwidth)
	maxload_dist_plot = plot(
		layer(
			x=maxload_sent,
			Geom.histogram(
				bincount=nbins_sent,
				limits=(min=0, max=nbins_sent*binwidth)
			),
			style(default_color=colors[1], alphas=[0.7])
		),
		layer(
			x=maxload_nosent,
			Geom.histogram(
				bincount=nbins_nosent,
				limits=(min=0, max=nbins_nosent*binwidth)
			),
			style(default_color=colors[2], alphas=[0.7])
		),
		layer(xintercept=[1.0], Geom.vline(color="red", size=1.0mm), order=10),
		Guide.title("Distribution of Maximum Load by Hospital"),
		Guide.xlabel("Maximum Load"),
		Guide.ylabel("Count"),
		Guide.manual_color_key("", ["With Transfers", "Without Transfers"], colors, shape=[Shape.square], size=[1.6mm]),
		style(
			minor_label_font_size=18px,
			major_label_font_size=22px,
			minor_label_font="CMU Serif",
			major_label_font="CMU Serif",
			key_label_font_size=16px,
			key_title_font_size=18px,
			key_label_font="CMU Serif",
			key_title_font="CMU Serif",
			background_color=colorant"white",
		),
	)

	bedtype_ext = (bedtype == :all) ? "" : (bedtype == :icu ? "_icu" : (bedtype == :ward ? "_ward" : ""))

	if display
		maxload_dist_plot |> SVG(23cm, 10cm)
	end
	if save
		save_all(maxload_dist_plot, 23cm, 10cm, "maxload_dist"*bedtype_ext, config)
	end
end;

function plot_active_total(config, data, results; bedtype=:all, display=true, save=true)
	if bedtype == :icu
		plot_title = "Total COVID ICU Nurses in $(config.region)"
		fn = "active_total_icu"
	elseif bedtype == :ward
		plot_title = "Total COVID Ward Nurses in $(config.region)"
		fn = "active_total_ward"
	else
		plot_title = "Total COVID Nurses in $(config.region)"
		fn = "active_total"
	end

	overall_plt = Gadfly.plot(
		layer(
			x = data.date_range,
			y = sum(data.demand, dims=1)[:],
			Geom.point, Geom.line,
			style(default_color=colorant"red"),
		),
		layer(
			x = data.date_range,
			y = fill(sum(data.nurses), data.T),
			Geom.point, Geom.line,
		),
		Guide.xlabel("Date"), Guide.ylabel("Nurses"),
		Coord.Cartesian(ymin=0),
		Guide.xticks(ticks=DateTime(data.start_date):Day(14):DateTime(data.end_date)),
		Scale.x_continuous(labels=x -> Dates.format(x, dateformat"mm-dd")),
		Scale.y_continuous(format=:plain),
		Guide.title(plot_title),
		style(
			key_position=:none,
			minor_label_font_size=20px,
			major_label_font_size=25px,
			minor_label_font="CMU Serif",
			major_label_font="CMU Serif",
			background_color=colorant"white",
		),
	)

	key_labels = ["COVID Nurses", "COVID Nurse Demand"]
	key_colors = [Scale.default_discrete_colors(1)[1], colorant"red"]
	key_h = 0.5cm
	key = make_horizontal_key(key_labels, key_colors, height=key_h)

	overall_plt = compose(
		context(0,0,1,1),
		compose(
			context(0, 0, 1w, 1h - key_h),
			render(overall_plt),
		),
		compose(
			context(0, 1h - key_h, 1w, key_h),
			key,
		),
	)

	overall_plt = fill_background(overall_plt)

	if display
		overall_plt |> SVG(20cm, 12cm)
	end
	if save
		save_all(overall_plt, 20cm, 12cm, fn, config)
	end
end;

function plot_active(config, data, results; add_title=true, columns=3, display=true, save=true, subset=nothing, bedtype=:all)
	if subset == nothing
		locs = enumerate(data.locations)
	else
		ind = [findfirst(==(h), data.locations) for h in subset]
		locs = zip(ind, subset)
	end

	active_plot_list = []
	colors = Scale.default_discrete_colors(5)[[1,4]]
	for (h_idx, h_name) in locs
		plot_title = length(h_name) > 43 ? strip(h_name[1:40])*"..." : h_name
		plt = plot(
			layer(x=data.date_range, y=results.active_nosent[h_idx,:], Geom.point, Geom.line, style(default_color=colors[1])),
			layer(x=data.date_range, y=results.active[h_idx,:], Geom.point, Geom.line, style(default_color=colors[2])),
			layer(x=data.date_range, y=data.demand[h_idx,:], Geom.point, Geom.line, style(default_color=colorant"red")),
			Guide.xticks(ticks=DateTime(data.start_date):Day(14):DateTime(data.end_date)),
			Scale.x_continuous(labels=x -> Dates.format(x, dateformat"mm-dd")),
			Coord.Cartesian(ymin=0),
			Guide.xlabel("Date"), Guide.ylabel("Nurses"),
			Guide.title(plot_title),
			style(
				key_position=:none,
				minor_label_font_size=16px,
				major_label_font_size=18px,
				minor_label_font="CMU Serif",
				major_label_font="CMU Serif",
				background_color=colorant"white",
			),
		)
		push!(active_plot_list, plt)
	end

	@assert columns == 2 || columns == 3
	plot_width = Dict(2 => 16, 3 => 12)[columns]
	active_plots, active_plots_w, active_plots_h = plot_grid(active_plot_list, columns, plot_width=plot_width, return_plot=true, display=false)

	active_key_colors = vcat(colors, colorant"red")
	active_key = make_horizontal_key(["Nurses (Without Transfers)", "Nurses (With Transfers)", "Nurse Demand"], active_key_colors, height=0.5cm)

	if bedtype == :all
		plot_title = "COVID Nurses With and Without Transfers"
	elseif bedtype == :icu
		plot_title = "COVID Nurses in the ICU With and Without Transfers"
	elseif bedtype == :ward
		plot_title = "COVID Nurses in Wards With and Without Transfers"
	end

	active_plots, dy1 = add_key_top(active_key, active_plots)
	active_plots, dy2 = add_title ? add_title_top(plot_title, active_plots) : (active_plots, 0cm)
	active_plots = fill_background(active_plots)

	active_plots_w = active_plots_w
	active_plots_h = active_plots_h + dy1 + dy2

	if display
		active_plots |> SVG(active_plots_w, active_plots_h)
	end

	if save
		out_name = isnothing(subset) ? "active" : "active_subset"
		out_name = (bedtype == :all) ? out_name : (bedtype == :icu ? out_name*"_icu" : (bedtype == :ward ? out_name*"_ward" : out_name))
		save_all(active_plots, active_plots_w, active_plots_h, out_name, config)
	end
end;

function plot_demand(config, data; add_title=true, columns=3, display=true, save=true, subset=nothing, bedtype=:all)
	if subset == nothing
		locs = enumerate(data.locations)
	else
		ind = [findfirst(==(h), data.locations) for h in subset]
		locs = zip(ind, subset)
	end

	active_plot_list = []
	for (h_idx, h_name) in locs
		plot_title = length(h_name) > 43 ? strip(h_name[1:40])*"..." : h_name
		plt = plot(
			layer(x=data.date_range, y=data.demand[h_idx,:], Geom.point, Geom.line, style(default_color=colorant"red")),
			Guide.xticks(ticks=DateTime(data.start_date):Day(14):DateTime(data.end_date)),
			Scale.x_continuous(labels=x -> Dates.format(x, dateformat"mm-dd")),
			Coord.Cartesian(ymin=0),
			Guide.xlabel("Date"), Guide.ylabel("Nurses"),
			Guide.title(plot_title),
			style(
				key_position=:none,
				minor_label_font_size=16px,
				major_label_font_size=18px,
				minor_label_font="CMU Serif",
				major_label_font="CMU Serif",
				background_color=colorant"white",
			),
		)
		push!(active_plot_list, plt)
	end

	@assert columns == 2 || columns == 3
	plot_width = Dict(2 => 16, 3 => 12)[columns]
	active_plots, active_plots_w, active_plots_h = plot_grid(active_plot_list, columns, plot_width=plot_width, return_plot=true, display=false)

	if bedtype == :all
		plot_title = "COVID Nurse Demand"
	elseif bedtype == :icu
		plot_title = "COVID Nurse Demand in the ICU"
	elseif bedtype == :ward
		plot_title = "COVID Nurse Demand in Wards"
	end

	active_plots, dy = add_title ? add_title_top(plot_title, active_plots) : (active_plots, 0cm)
	active_plots = fill_background(active_plots)

	active_plots_w = active_plots_w
	active_plots_h = active_plots_h + dy

	if display
		active_plots |> SVG(active_plots_w, active_plots_h)
	end

	if save
		out_name = isnothing(subset) ? "demand" : "demand_subset"
		out_name = (bedtype == :all) ? out_name : (bedtype == :icu ? out_name*"_icu" : (bedtype == :ward ? out_name*"_ward" : out_name))
		save_all(active_plots, active_plots_w, active_plots_h, out_name, config)
	end
end;

function plot_load(config, data, results; display=true, save=true, subset=nothing, bedtype=:all)
	if subset == nothing
		results_table = results.table
		locations = data.locations
		nlocations = data.N
	else
		results_table = filter(row -> row.location in subset, results.table)
		locations = subset
		nlocations = length(subset)
	end

	## load plots ##
	load_plot_max = 3.0
	plt_nosent = Gadfly.plot(
		results_table, x=:date, y=:load_nosent, color=:location,
		Geom.line,
		Guide.xlabel("Date"), Guide.ylabel("Nurse Load"),
		Guide.title("COVID Nurse Load in $(config.region) - Without Transfers"),
		layer(ymin=[0.0], ymax=[1.0],           Geom.hband, alpha=[0.1], Theme(default_color="green"), order=-1),
		layer(ymin=[1.0], ymax=[load_plot_max], Geom.hband, alpha=[0.1], Theme(default_color="red"),   order=-1),
		Coord.cartesian(xmin=data.start_date+Day(1), xmax=data.end_date-Day(1), ymin=0, ymax=load_plot_max),
		Guide.xticks(ticks=DateTime(data.start_date):Day(14):DateTime(data.end_date)),
		Scale.x_continuous(labels=x -> Dates.format(x, dateformat"mm-dd")),
		style(
			key_position=:none,
			minor_label_font_size=20px,
			major_label_font_size=25px,
			minor_label_font="CMU Serif",
			major_label_font="CMU Serif",
		),
	)
	plt_sent = Gadfly.plot(
		results_table, x=:date, y=:load_sent, color=:location,
		Geom.line,
		Guide.xlabel("Date"), Guide.ylabel("Nurse Load"),
		Guide.title("COVID Nurse Load in $(config.region) - With Transfers"),
		layer(ymin=[0.0], ymax=[1.0],           Geom.hband, alpha=[0.1], Theme(default_color="green"), order=-1),
		layer(ymin=[1.0], ymax=[load_plot_max], Geom.hband, alpha=[0.1], Theme(default_color="red"),   order=-1),
		Coord.cartesian(xmin=data.start_date+Day(1), xmax=data.end_date-Day(1), ymin=0, ymax=load_plot_max),
		Guide.xticks(ticks=DateTime(data.start_date):Day(14):DateTime(data.end_date)),
		Scale.x_continuous(labels=x -> Dates.format(x, dateformat"mm-dd")),
		style(
			key_position=:none,
			minor_label_font_size=20px,
			major_label_font_size=25px,
			minor_label_font="CMU Serif",
			major_label_font="CMU Serif",
		),
	)

	## load key ##
	load_key_labels = map(h_name -> length(h_name) > 53 ? strip(h_name[1:50])*"..." : h_name, locations)
	load_key, load_key_h = make_key(load_key_labels, Scale.default_discrete_colors(nlocations), line_height=0.5cm, ncolumns=3, shape=:circle)

	load_plts_withkey = compose(
		context(0,0,1,1),
		compose(
			context(0, 0, 1w, 1h - load_key_h),
			compose(context(0, 0, 0.5w, 1cy), render(plt_nosent)),
			compose(context(0.5w, 0, 0.5w, 1cy), render(plt_sent)),
		),
		compose(
			context(0, 1h - load_key_h, 1w, load_key_h),
			load_key,
		),
	)

	load_plts_withkey = fill_background(load_plts_withkey)

	if display
		load_plts_withkey |> SVG(40cm, 14cm + load_key_h)
	end

	if save
		out_name = isnothing(subset) ? "load_bylocation" : "load_bylocation_subset"
		out_name = (bedtype == :all) ? out_name : (bedtype == :icu ? out_name*"_icu" : (bedtype == :ward ? out_name*"_ward" : out_name))
		save_all(load_plts_withkey, 40cm, 14cm + load_key_h, out_name, config)
	end
end;

function plot_sent_total(config, data, results; bedtype=:all, display=true, save=true)
	sent_plot = plot(
		x=data.date_range,
		y=sum(results.sent, dims=[1,2])[:],
		Geom.point, Geom.line,
		Guide.xticks(ticks=DateTime(data.start_date):Day(14):DateTime(data.end_date)),
		Scale.x_continuous(labels=x -> Dates.format(x, dateformat"mm-dd")),
		Coord.Cartesian(ymin=0),
		Guide.xlabel("Date"), Guide.ylabel("Nurse Transfers"),
		Guide.title("Total Nurse Transfers"),
		style(
			key_position=:none,
			minor_label_font_size=16px,
			major_label_font_size=18px,
			minor_label_font="CMU Serif",
			major_label_font="CMU Serif",
			background_color=colorant"white",
		),
	)

	bedtype_ext = (bedtype == :all) ? "" : (bedtype == :icu ? "_icu" : (bedtype == :ward ? "_ward" : ""))

	if display
		sent_plot |> SVG(20cm, 10cm)
	end
	if save
		save_all(sent_plot, 20cm, 10cm, "transfers_total"*bedtype_ext, config)
	end
end;

function plot_transfers(config, data, results; add_title=true, columns=3, display=true, save=true, subset=nothing, bedtype=:all)
	if subset == nothing
		locs = collect(enumerate(data.locations))
	else
		ind = [findfirst(==(h), data.locations) for h in subset]
		locs = zip(ind, subset)
	end

	transfer_plot_list = []
	colors = [colorant"seagreen4", colorant"firebrick1"]
	for (h_idx, h_name) in locs
		plot_title = length(h_name) > 43 ? strip(h_name[1:40])*"..." : h_name
		plt = plot(
			layer(x=data.date_range, y=sum(results.sent[:,h_idx,:], dims=1)[:], Geom.point, Geom.line, style(default_color=colors[1])),
			layer(x=data.date_range, y=sum(results.sent[h_idx,:,:], dims=1)[:], Geom.point, Geom.line, style(default_color=colors[2])),
			Guide.xticks(ticks=DateTime(data.start_date):Day(14):DateTime(data.end_date)),
			Scale.x_continuous(labels=x -> Dates.format(x, dateformat"mm-dd")),
			Coord.Cartesian(ymin=0),
			Guide.xlabel("Date"), Guide.ylabel("Nurses"),
			Guide.title(plot_title),
			style(
				key_position=:none,
				minor_label_font_size=16px,
				major_label_font_size=18px,
				minor_label_font="CMU Serif",
				major_label_font="CMU Serif",
				background_color=colorant"white",
			),
		)
		push!(transfer_plot_list, plt)
	end

	@assert columns == 2 || columns == 3
	plot_width = Dict(2 => 16, 3 => 12)[columns]
	transfer_plot, plt_w, plt_h = plot_grid(transfer_plot_list, columns, plot_width=plot_width, return_plot=true, display=false)

	transfer_key = make_horizontal_key(["Received", "Sent"], colors, height=0.5cm)

	transfer_plot, dy1 = add_key_top(transfer_key, transfer_plot)
	transfer_plot, dy2 = add_title ? add_title_top("Nurses Sent and Received by Hospital", transfer_plot) : (transfer_plot, 0cm)
	transfer_plot = fill_background(transfer_plot)

	transfer_plot_w = plt_w
	transfer_plot_h = plt_h + dy1 + dy2

	if display
		transfer_plot |> SVG(transfer_plot_w, transfer_plot_h)
	end

	if save
		out_name = isnothing(subset) ? "transfers" : "transfers_subset"
		out_name = (bedtype == :all) ? out_name : (bedtype == :icu ? out_name*"_icu" : (bedtype == :ward ? out_name*"_ward" : out_name))
		save_all(transfer_plot, transfer_plot_w, transfer_plot_h, out_name, config)
	end
end;

function plot_patients_per_nurse(config, combined_data; display=true, save=true)
	nbins = 16
	plt1 = plot(
		combined_data, x=:patients_per_nurse_icu,
		Geom.histogram(bincount=nbins),
		Guide.xlabel("Patient to Nurse Ratio"),
		Guide.ylabel("Frequency"),
		Guide.title("Patients per Nurse Distribution - ICU"),
		Coord.cartesian(xmin=0, xmax=4),
		style(
			minor_label_font_size=16px,
			major_label_font_size=18px,
			minor_label_font="CMU Serif",
			major_label_font="CMU Serif",
			background_color=colorant"white",
		),
	)
	plt2 = plot(
		combined_data, x=:patients_per_nurse_ward,
		Geom.histogram(bincount=nbins),
		Guide.xlabel("Patient to Nurse Ratio"),
		Guide.ylabel("Frequency"),
		Guide.title("Patients per Nurse Distribution - Ward"),
		Coord.cartesian(xmin=0, xmax=8),
		style(
			minor_label_font_size=16px,
			major_label_font_size=18px,
			minor_label_font="CMU Serif",
			major_label_font="CMU Serif",
			background_color=colorant"white",
		),
	)
	plt = hstack(plt1, plt2)

	if display
		plt |> SVG(28cm, 8cm)
	end
	if save
		config = merge(config, (experiment = "data",))
		save_all(plt, 28cm, 8cm, "patients_per_nurse", config)
	end
end;

end;

ext: jl

max_stars_repo_name: flixpar/covid-resource-allocation

max_issues_repo_name: flixpar/covid-resource-allocation

max_forks_repo_name: flixpar/covid-resource-allocation


text:
using Plots, Interact, Mux

app = @manipulate for angle=0:0.01:2pi
    angles = range(0, angle, length=100)
    x = cos.(angles)
    y = sin.(angles)
    plot(
        x, y,
        ratio=:equal,
        xlims=(-1.5, 1.5), ylims=(-1.5, 1.5),
        legend=:none,
        framestyle=:none)
end

WebIO.webio_serve(page("/", req -> app), 8085)


text:
function Smat(n, lam)
    if lam == 0
        e = ones(n)
        S = myspdiagm(spdiagm(0 => 0.5*e, 2 => -0.5*e), n)
        S[1] = 1
    else
        e = lam./(collect(0:n-1).+lam)
        S = myspdiagm(spdiagm(0 => e, 2 => -e), n)
    end
    return S
end

myspdiagm(M, N) = size(M, 1) !== N ? (return M[1:N, 1:N]) : (return M)


text:
using ArgParse
using PowerModelsSecurityConstrained

function scopf_main(parsed_args)
    ini_file = parsed_args["file"]
    scenario = parsed_args["scenario"]
    scoring_method = parsed_args["scoring-method"]
    time_limit_sol1 = parsed_args["time-limit"]
    time_limit_total = parsed_args["hard-time-limit"]
    distribute = parsed_args["distribute"]
    skip_solution2 = parsed_args["skip-solution2"]
    remove_solutions = parsed_args["remove-solutions"]
    gurobi = parsed_args["gurobi"]

    files, scenario_id = find_goc_files(ini_file, scenario_id=scenario)

    ini_dir = dirname(ini_file)
    output_dir = joinpath(ini_dir, scenario_id)

    info(LOGGER, "output dir: $(output_dir)")

    time_sol1_start = time()
    compute_solution1(files["con"], files["inl"], files["raw"], files["rop"], time_limit_sol1-60, scoring_method, "network name"; output_dir=output_dir, scenario_id=scenario_id, gurobi=gurobi)
    time_sol1 = time() - time_sol1_start

    if !skip_solution2
        compute_solution2(files["con"], files["inl"], files["raw"], files["rop"], trunc(Int, time_limit_total-time_sol1-60), scoring_method, "network name"; output_dir=output_dir, scenario_id=scenario_id)
    end

    if remove_solutions
        remove_solution_files(output_dir=output_dir)
        remove_detail_file(output_dir=output_dir)
    end
end

function parse_scopf_commandline()
    s = ArgParseSettings()

    @add_arg_table! s begin
        "--file", "-f"
            help = "the data initiation file (.ini)"
            required = true
        "--scenario", "-s"
            help = "the scenario to run (directory)"
            default = ""
        "--scoring-method", "-m"
            help = "the objective function (1,2,3,4)"
            arg_type = Int
            default = 2
        "--time-limit", "-t"
            help = "solution 1 runtime limit (seconds)"
            arg_type = Int
            default = 600000 #2700
        "--hard-time-limit"
            help = "solution 1 and 2 runtime limit (seconds)"
            arg_type = Int
            default = 600000 #34200
        "--distribute", "-d"
            help = "run on multiple processes"
            action = :store_true
        "--skip-solution2"
            help = "skip computation of solution2 file"
            action = :store_true
        "--remove-solutions"
            help = "delete the solution files after competition"
            action = :store_true
        "--gurobi", "-g"
            help = "use Gurobi for solving lp, qp and mip problems"
            action = :store_true
            #default = false
    end

    return parse_args(s)
end


text:
T=3
eps_exp=(1.)/(3.)
tsname="joBlock - basic"
@testset "$tsname" begin
for t=1:T # start test loop

    cn=rand(2:9,3)
    rm=rand(2:9,3)
    a1=rand(ComplexF64,rm[1],cn[1])
    A1=joMatrix(a1;DDT=ComplexF32,RDT=ComplexF64)
    a2=rand(ComplexF64,rm[1],cn[2])
    A2=joMatrix(a2;DDT=ComplexF32,RDT=ComplexF64)
    a3=rand(ComplexF64,rm[1],cn[3])
    A3=joMatrix(a3;DDT=ComplexF32,RDT=ComplexF64)
    b1=rand(ComplexF64,rm[2],cn[2])
    B1=joMatrix(b1;DDT=ComplexF32,RDT=ComplexF64)
    b2=rand(ComplexF64,rm[2],cn[3])
    B2=joMatrix(b2;DDT=ComplexF32,RDT=ComplexF64)
    b3=rand(ComplexF64,rm[2],cn[1])
    B3=joMatrix(b3;DDT=ComplexF32,RDT=ComplexF64)
    c1=rand(ComplexF64,rm[3],cn[3])
    C1=joMatrix(c1;DDT=ComplexF32,RDT=ComplexF64)
    c2=rand(ComplexF64,rm[3],cn[1])
    C2=joMatrix(c2;DDT=ComplexF32,RDT=ComplexF64)
    c3=rand(ComplexF64,rm[3],cn[2])
    C3=joMatrix(c3;DDT=ComplexF32,RDT=ComplexF64)
    bd=[a1 a2 a3; b1 b2 b3; c1 c2 c3]
    BD=joBlock([3,3,3],A1,A2,A3,B1,B2,B3,C1,C2,C3)
    BDm=BD.m
    BDn=BD.n

    verbose && println("$tsname ($BDm,$BDn)")
    @testset "$BDm x $BDn" begin
        vn=rand(ComplexF32,BDn);
        vm=rand(ComplexF64,BDm);
        mvn=rand(ComplexF32,BDn,2);
        mvm=rand(ComplexF64,BDm,2);
        @test norm(bd*vn-BD*vn)<eps(norm(bd*vn))^eps_exp
        @test norm(bd*mvn-BD*mvn)<eps(norm(bd*mvn))^eps_exp
        @test norm(transpose(bd)*vm-transpose(BD)*vm)<eps(norm(transpose(bd)*vm))^eps_exp
        @test norm(transpose(bd)*mvm-transpose(BD)*mvm)<eps(norm(transpose(bd)*mvm))^eps_exp
        @test norm(adjoint(bd)*vm-adjoint(BD)*vm)<eps(norm(adjoint(bd)*vm))^eps_exp
        @test norm(adjoint(bd)*mvm-adjoint(BD)*mvm)<eps(norm(adjoint(bd)*mvm))^eps_exp
        @test norm(conj(bd)*vn-conj(BD)*vn)<eps(norm(conj(bd)*vn))^eps_exp
        @test norm(conj(bd)*mvn-conj(BD)*mvn)<eps(norm(conj(bd)*mvn))^eps_exp
        @test norm((-bd)*vn-(-BD)*vn)<eps(norm((-bd)*vn))^eps_exp
        @test norm((-bd)*mvn-(-BD)*mvn)<eps(norm((-bd)*mvn))^eps_exp
    end

end # end test loop
end
tsname="joBlock - weighted"
@testset "$tsname" begin
for t=1:T # start test loop

    cn=rand(2:9,3)
    rm=rand(2:9,3)
    w=rand(ComplexF64,9)
    a1=rand(ComplexF64,rm[1],cn[1])
    A1=joMatrix(a1;DDT=ComplexF32,RDT=ComplexF64)
    a2=rand(ComplexF64,rm[1],cn[2])
    A2=joMatrix(a2;DDT=ComplexF32,RDT=ComplexF64)
    a3=rand(ComplexF64,rm[1],cn[3])
    A3=joMatrix(a3;DDT=ComplexF32,RDT=ComplexF64)
    b1=rand(ComplexF64,rm[2],cn[2])
    B1=joMatrix(b1;DDT=ComplexF32,RDT=ComplexF64)
    b2=rand(ComplexF64,rm[2],cn[3])
    B2=joMatrix(b2;DDT=ComplexF32,RDT=ComplexF64)
    b3=rand(ComplexF64,rm[2],cn[1])
    B3=joMatrix(b3;DDT=ComplexF32,RDT=ComplexF64)
    c1=rand(ComplexF64,rm[3],cn[3])
    C1=joMatrix(c1;DDT=ComplexF32,RDT=ComplexF64)
    c2=rand(ComplexF64,rm[3],cn[1])
    C2=joMatrix(c2;DDT=ComplexF32,RDT=ComplexF64)
    c3=rand(ComplexF64,rm[3],cn[2])
    C3=joMatrix(c3;DDT=ComplexF32,RDT=ComplexF64)
    bd=[w[1]*a1 w[2]*a2 w[3]*a3; w[4]*b1 w[5]*b2 w[6]*b3; w[7]*c1 w[8]*c2 w[9]*c3]
    BD=joBlock([3,3,3],A1,A2,A3,B1,B2,B3,C1,C2,C3;weights=w)
    BDm=BD.m
    BDn=BD.n

    verbose && println("$tsname ($BDm,$BDn)")
    @testset "$BDm x $BDn" begin
        vn=rand(ComplexF32,BDn);
        vm=rand(ComplexF64,BDm);
        mvn=rand(ComplexF32,BDn,2);
        mvm=rand(ComplexF64,BDm,2);
        @test norm(bd*vn-BD*vn)<eps(norm(bd*vn))^eps_exp
        @test norm(bd*mvn-BD*mvn)<eps(norm(bd*mvn))^eps_exp
        @test norm(transpose(bd)*vm-transpose(BD)*vm)<eps(norm(transpose(bd)*vm))^eps_exp
        @test norm(transpose(bd)*mvm-transpose(BD)*mvm)<eps(norm(transpose(bd)*mvm))^eps_exp
        @test norm(adjoint(bd)*vm-adjoint(BD)*vm)<eps(norm(adjoint(bd)*vm))^eps_exp
        @test norm(adjoint(bd)*mvm-adjoint(BD)*mvm)<eps(norm(adjoint(bd)*mvm))^eps_exp
        @test norm(conj(bd)*vn-conj(BD)*vn)<eps(norm(conj(bd)*vn))^eps_exp
        @test norm(conj(bd)*mvn-conj(BD)*mvn)<eps(norm(conj(bd)*mvn))^eps_exp
        @test norm((-bd)*vn-(-BD)*vn)<eps(norm((-bd)*vn))^eps_exp
        @test norm((-bd)*mvn-(-BD)*mvn)<eps(norm((-bd)*mvn))^eps_exp
    end

end # end test loop
end


text:
import JSON2

using Codex
using Test

@testset "Codex" begin
    include("gitlab.jl")
    include("utils.jl")
    include("stats.jl")
    include(joinpath("TransformExport", "transformexport.jl"))
    include(joinpath("Questionnaires", "questionnaires.jl"))
end


text:
using DataStructures
file = readlines("day_17.txt")

function get_neig(pos)
    return [(pos[1]+i,pos[2]+j,pos[3]+k) for i =-1:1 for j =-1:1 for k =-1:1
        if  !(i==0 && j==0 && k==0)]
end

cubes = DefaultDict(false)
recently_change = []

#reading input
for (index_i, line) in enumerate(file)
    for (index_j,word) in enumerate(line)
        if (word == '#')
            cubes[(index_i,index_j,0)] = true
            push!(recently_change,(index_i,index_j,0))
            append!(recently_change, get_neig((index_i,index_j,0)))
        end
    end
end

for cycle in 1:6
    new_cubes = DefaultDict(false)
    new_recently_change = []

    for pos in recently_change
        neig = get_neig(pos)
        sum_Act = sum([cubes[x] for x in neig])

        if cubes[pos] == true
            if sum_Act == 2 || sum_Act == 3
                new_cubes[pos] = true
                append!(new_recently_change,neig)
                push!(new_recently_change,pos)
            end
        else
            if sum_Act == 3
                new_cubes[pos] = true
                append!(new_recently_change,neig)
                push!(new_recently_change,pos)
            end
        end
    end

    global cubes = new_cubes
    global recently_change = new_recently_change
    println(cycle)
end

#heay, i could have made a better detect for possible active cubes, oh well
aa = 0
for (key, value) in cubes
    global aa+=value
end
print(aa)

#QUESTION 2 (the same thing but 3d)

file = readlines("day_17.txt")
function get_neig(pos)
    return [ (pos[1]+i,pos[2]+j,pos[3]+k,pos[4]+l)
        for i =-1:1 for j =-1:1 for k =-1:1 for l=-1:1
        if  !(i==0 && j==0 && k==0 && l==0)]
end

cubes = DefaultDict(false)

#reading input
for (index_i, line) in enumerate(file)
    for (index_j,word) in enumerate(line)
        if (word == '#')
            pos = (index_i,index_j,0,0)
            cubes[pos] = true
        end
    end
end

for cycle in 1:6
    new_cubes = DefaultDict(false)
    recently_change = DefaultDict(0)

    #count neig acc
    for (pos, value) in cubes
        if value == false
            print("AAAAAAAA")
        end
        neig = get_neig(pos)

        for i in neig
            recently_change[i] += value
        end
    end

    #find new cycle value
    for (pos, neig_acc) in recently_change
        curr_act = cubes[pos]
        new_act = false

        new_act |= (curr_act == true &&(neig_acc == 3 || neig_acc == 2))
        new_act |= (curr_act == false && neig_acc == 3 )

        if(new_act == true)
            new_cubes[pos] = true
        end
    end

    global cubes = new_cubes
    println(cycle)
end

aa = 0
for (key, value) in cubes
    global aa+=value
end
print(aa)


text:
module FockPackage

using LinearAlgebra: dot
using Printf: @printf, @sprintf
using ..Spatials: PID, AbstractBond, Bond, decompose
using ..DegreesOfFreedom: IID, Index, Internal, FilteredAttributes, IDFConfig, Table, OID, Operator, Operators, LaTeX, latexformat, coordpresent
using ..Terms: wildcard, constant, Subscript, Subscripts, subscriptexpr, Coupling, Couplings, @subscript, couplingcenters, Term, TermCouplings, TermAmplitude, TermModulate
using ...Interfaces: id, rank, kind
using ...Prerequisites: Float, delta, decimaltostr
using ...Mathematics.AlgebraOverFields: SimpleID, ID, Element
using ...Mathematics.VectorSpaces: VectorSpace, IsMultiIndexable, MultiIndexOrderStyle

import ..DegreesOfFreedom: script, otype, isHermitian
import ..Terms: couplingcenter, statistics, abbr, termfactor
import ...Interfaces: dims, inds, ⊗, ⋅, expand, expand!, permute

export ANNIHILATION, CREATION, MAJORANA, FID, FIndex, Fock
export usualfockindextotuple, nambufockindextotuple
export FockOperator, FOperator, BOperator, isnormalordered
export foptdefaultlatex, boptdefaultlatex
export FCID, FockCoupling
export σ⁰, σˣ, σʸ, σᶻ, σ⁺, σ⁻
export @fc_str, @σ⁰_str, @σˣ_str, @σʸ_str, @σᶻ_str, @σ⁺_str, @σ⁻_str
export Onsite, Hopping, Pairing
export hubbard, Hubbard
export interorbitalinterspin, InterOrbitalInterSpin
export interorbitalintraspin, InterOrbitalIntraSpin
export spinflip, SpinFlip
export pairhopping, PairHopping
export Coulomb
export FFockTerm, BFockTerm

"""
    ANNIHILATION

Indicate that the nambu index is ANNIHILATION.
"""
const ANNIHILATION = 1
"""
    CREATION

Indicate that the nambu index is CREATION.
"""
const CREATION = 2

"""
    MAJORANA

Indicate that the nambu index is MAJORANA.
"""
const MAJORANA = 0

"""
    FID <: IID

The Fock id.
"""
struct FID <: IID
    orbital::Int
    spin::Int
    nambu::Int
    FID(orbital::Int, spin::Int, nambu::Int) = (@assert nambu ∈ (0, 1, 2) "FID error: wrong input nambu($nambu)."; new(orbital, spin, nambu))
end
Base.fieldnames(::Type{FID}) = (:orbital, :spin, :nambu)

"""
    FID(; orbital::Int=1, spin::Int=1, nambu::Int=ANNIHILATION)

Create a Fock id.
"""
FID(; orbital::Int=1, spin::Int=1, nambu::Int=ANNIHILATION) = FID(orbital, spin, nambu)

"""
    adjoint(fid::FID) -> FID

Get the adjoint of a Fock id.
"""
Base.adjoint(fid::FID) = FID(fid.orbital, fid.spin, (fid.nambu == 0) ? 0 : 3-fid.nambu)

"""
    Fock <: Internal{FID}

The Fock interanl degrees of freedom.
"""
struct Fock <: Internal{FID}
    atom::Int
    norbital::Int
    nspin::Int
    nnambu::Int
    Fock(atom::Int, norbital::Int, nspin::Int, nnambu::Int) = (@assert nnambu ∈ (1, 2) "Fock error: wrong input nnambu($nnambu)."; new(atom, norbital, nspin, nnambu))
end
IsMultiIndexable(::Type{Fock}) = IsMultiIndexable(true)
MultiIndexOrderStyle(::Type{Fock}) = MultiIndexOrderStyle('C')
dims(fock::Fock) = (fock.norbital, fock.nspin, fock.nnambu)
inds(fid::FID, fock::Fock) = (fid.orbital, fid.spin, (fock.nnambu == 1) ? 1 : fid.nambu)
FID(inds::NTuple{3, Int}, fock::Fock) = FID(inds[1], inds[2], (fock.nnambu == 1) ? 0 : inds[3])

"""
    Fock(; atom::Int=1, norbital::Int=1, nspin::Int=2, nnambu::Int=2)

Construct a Fock degrees of freedom.
"""
Fock(; atom::Int=1, norbital::Int=1, nspin::Int=2, nnambu::Int=2) = Fock(atom, norbital, nspin, nnambu)

"""
    FIndex{S} <: Index{PID{S}, FID}

The Fock index.
"""
struct FIndex{S} <: Index{PID{S}, FID}
    scope::S
    site::Int
    orbital::Int
    spin::Int
    nambu::Int
    function FIndex(scope::S, site::Int, orbital::Int, spin::Int, nambu::Int) where S
        @assert nambu ∈ (0, 1, 2) "FIndex error: wrong input nambu($nambu)."
        new{S}(scope, site, orbital, spin, nambu)
    end
end
Base.fieldnames(::Type{<:FIndex}) = (:scope, :site, :orbital, :spin, :nambu)

"""
    union(::Type{P}, ::Type{FID}) where {P<:PID}

Get the union type of `PID` and `FID`.
"""
Base.union(::Type{P}, ::Type{FID}) where {P<:PID} = FIndex{fieldtype(P, :scope)}

"""
    script(oid::OID{<:FIndex}, ::Val{:site}) -> Int
    script(oid::OID{<:FIndex}, ::Val{:orbital}) -> Int
    script(oid::OID{<:FIndex}, ::Val{:spinint}) -> Int
    script(oid::OID{<:FIndex}, ::Val{:spinsym}) -> String
    script(oid::OID{<:FIndex}, ::Val{:nambu}) -> String

Get the required script of an Fock oid.
"""
script(oid::OID{<:FIndex}, ::Val{:site}) = oid.index.site
script(oid::OID{<:FIndex}, ::Val{:orbital}) = oid.index.orbital
script(oid::OID{<:FIndex}, ::Val{:spinint}) = oid.index.spin
script(oid::OID{<:FIndex}, ::Val{:spinsym}) = (oid.index.spin == 1) ? "↓" : (oid.index.spin == 2) ? "↑" : error("script error: spin($oid.index.spin) not in (1, 2).")
script(oid::OID{<:FIndex}, ::Val{:nambu}) = (oid.index.nambu == CREATION) ? "\\dagger" : ""

"""
    angle(id::OID{<:FIndex}, vectors::AbstractVector{<:AbstractVector{Float}}, values::AbstractVector{Float}) -> Complex{Float}

Get the twist phase corresponding to a Fock oid.
"""
function Base.angle(id::OID{<:FIndex}, vectors::AbstractVector{<:AbstractVector{Float}}, values::AbstractVector{Float})
    phase = (length(vectors) == 1) ? 2pi*dot(decompose(id.icoord, vectors[1]), values) :
            (length(vectors) == 2) ? 2pi*dot(decompose(id.icoord, vectors[1], vectors[2]), values) :
            (length(vectors) == 3) ? 2pi*dot(decompose(id.icoord, vectors[1], vectors[2], vectors[3]), values) :
            error("angle error: not supported number of input basis vectors.")
    (id.index.nambu == ANNIHILATION) ? phase : (id.index.nambu == CREATION) ? -phase : error("angle error: not supported Fock index.")
end

"""
    usualfockindextotuple

Indicate that the filtered attributes are `(:scope, :site, :orbital, :spin)` when converting a Fock index to tuple.
"""
const usualfockindextotuple = FilteredAttributes(:scope, :site, :orbital, :spin)
"""
    nambufockindextotuple

Indicate that the filtered attributes are `(:scope, :nambu, :site, :orbital, :spin)` when converting a Fock index to tuple.
"""
const nambufockindextotuple = FilteredAttributes(:scope, :nambu, :site, :orbital, :spin)

"""
    FockOperator{V<:Number, I<:ID{OID}} <: Operator{V, I}

Abstract type for all Fock operators.
"""
abstract type FockOperator{V<:Number, I<:ID{OID}} <: Operator{V, I} end

"""
    isnormalordered(opt::FockOperator) -> Bool

Judge whether a FockOperator is normal ordered.
"""
function isnormalordered(opt::FockOperator)
    flag = true
    for i = 1:rank(opt)
        flag && (opt.id[i].index.nambu == ANNIHILATION) && (flag = false)
        flag || (opt.id[i].index.nambu == CREATION) && return false
    end
    return true
end

"""
    FOperator(value::Number, id::ID{OID}=ID())

Fermionic Fock operator.
"""
struct FOperator{V<:Number, I<:ID{OID}} <: FockOperator{V, I}
    value::V
    id::I
    FOperator(value::Number) = new{typeof(value), Tuple{}}(value, ())
    FOperator(value::Number, id::Tuple{}) = new{typeof(value), Tuple{}}(value, ())
    FOperator(value::Number, id::ID{OID}) = new{typeof(value), typeof(id)}(value, id)
end

"""
    foptdefaultlatex

The default LaTeX format of the oids of a fermionic operator.
"""
const foptdefaultlatex = LaTeX{(:nambu,), (:site, :orbital, :spinsym)}('c')
latexformat(FOperator, foptdefaultlatex)

"""
    statistics(opt::FOperator) -> Char
    statistics(::Type{<:FOperator}) -> Char

Get the statistics of FOperator.
"""
statistics(opt::FOperator) = opt |> typeof |> statistics
statistics(::Type{<:FOperator}) = 'F'

"""
    permute(::Type{<:FOperator}, id1::OID{<:FIndex}, id2::OID{<:FIndex}, ::Any=nothing) -> Tuple{Vararg{FOperator}}

Permute two fermionic oid and get the result.
"""
function permute(::Type{<:FOperator}, id1::OID{<:FIndex}, id2::OID{<:FIndex}, ::Any=nothing)
    @assert (id1.index ≠ id2.index) || (id1.rcoord ≠ id2.rcoord) || (id1.icoord ≠ id2.icoord) "permute error: permuted ids should not be equal to each other."
    if (id1.index' == id2.index) && (id1.rcoord == id2.rcoord) && (id1.icoord == id2.icoord)
        return (FOperator(1), FOperator(-1, ID(id2, id1)))
    else
        return (FOperator(-1, ID(id2, id1)),)
    end
end

"""
    *(f1::FOperator, f2::FOperator) -> Union{Nothing, FOperator}

Get the multiplication of two fermionic Fock operators.
"""
function Base.:*(f1::FOperator, f2::FOperator)
    (rank(f1) > 0) && (rank(f2) > 0) && (f1.id[end] == f2.id[1]) && return nothing
    return invoke(*, Tuple{Element, Element}, f1, f2)
end

"""
    BOperator(value::Number, id::ID{OID}=ID())

Bosonic Fock operator.
"""
struct BOperator{V<:Number, I<:ID{OID}} <: FockOperator{V, I}
    value::V
    id::I
    BOperator(value::Number) = new{typeof(value), Tuple{}}(value, ())
    BOperator(value::Number, id::Tuple{}) = new{typeof(value), Tuple{}}(value, ())
    BOperator(value::Number, id::ID{OID}) = new{typeof(value), typeof(id)}(value, id)
end

"""
    boptdefaultlatex

The default LaTeX format of the oids of a bosonic operator.
"""
const boptdefaultlatex = LaTeX{(:nambu,), (:site, :orbital, :spinsym)}('b')
latexformat(BOperator, boptdefaultlatex)

"""
    statistics(opt::BOperator)
    statistics(::Type{<:BOperator})

Get the statistics of BOperator.
"""
statistics(opt::BOperator) = opt |> typeof |> statistics
statistics(::Type{<:BOperator}) = 'B'

"""
    permute(::Type{<:BOperator}, id1::OID{<:FIndex}, id2::OID{<:FIndex}, ::Any=nothing) -> Tuple{Vararg{BOperator}}

Permute two bosonic oid and get the result.
"""
function permute(::Type{<:BOperator}, id1::OID{<:FIndex}, id2::OID{<:FIndex}, ::Any=nothing)
    @assert (id1.index ≠ id2.index) || (id1.rcoord ≠ id2.rcoord) || (id1.icoord ≠ id2.icoord) "permute error: permuted ids should not be equal to each other."
    if (id1.index' == id2.index) && (id1.rcoord == id2.rcoord) && (id1.icoord == id2.icoord)
        if id1.index.nambu == CREATION
            return (BOperator(1), BOperator(1, ID(id2, id1)))
        else
            return (BOperator(-1), BOperator(1, ID(id2, id1)))
        end
    else
        return (BOperator(1, ID(id2, id1)),)
    end
end

"""
    FCID(; center=wildcard, atom=wildcard, orbital=wildcard, spin=wildcard, nambu=wildcard, obsub=wildcard, spsub=wildcard)

The id of a Fock coupling.
"""
struct FCID{C, A, O, S, N, OS, SS} <: SimpleID
    center::C
    atom::A
    orbital::O
    spin::S
    nambu::N
    obsub::OS
    spsub::SS
end
Base.fieldnames(::Type{<:FCID}) = (:center, :atom, :orbital, :spin, :nambu, :obsub, :spsub)
FCID(; center=wildcard, atom=wildcard, orbital=wildcard, spin=wildcard, nambu=wildcard, obsub=wildcard, spsub=wildcard) = FCID(center, atom, orbital, spin, nambu, obsub, spsub)

"""
    FockCoupling(value::Number, id::ID{FCID}, obsubscripts::Subscripts, spsubscripts::Subscripts)
    FockCoupling{N}(value::Number=1;
                    centers::Union{NTuple{N, Int}, Nothing}=nothing,
                    atoms::Union{NTuple{N, Int}, Nothing}=nothing,
                    orbitals::Union{NTuple{N, Int}, Subscript, Nothing}=nothing,
                    spins::Union{NTuple{N, Int}, Subscript, Nothing}=nothing,
                    nambus::Union{NTuple{N, Int}, Nothing}=nothing
                    ) where N

Fock coupling.
"""
struct FockCoupling{V, I<:ID{FCID}, OS<:Subscripts, SS<:Subscripts} <: Coupling{V, I}
    value::V
    id::I
    obsubscripts::OS
    spsubscripts::SS
    function FockCoupling(value::Number, id::ID{FCID}, obsubscripts::Subscripts, spsubscripts::Subscripts)
        new{typeof(value), typeof(id), typeof(obsubscripts), typeof(spsubscripts)}(value, id, obsubscripts, spsubscripts)
    end
end
function FockCoupling{N}(value::Number=1;
                        centers::Union{NTuple{N, Int}, Nothing}=nothing,
                        atoms::Union{NTuple{N, Int}, Nothing}=nothing,
                        orbitals::Union{NTuple{N, Int}, Subscript, Nothing}=nothing,
                        spins::Union{NTuple{N, Int}, Subscript, Nothing}=nothing,
                        nambus::Union{NTuple{N, Int}, Nothing}=nothing
                        ) where N
    (centers === nothing) && (centers = ntuple(i->wildcard, Val(N)))
    (atoms === nothing) && (atoms = ntuple(i->wildcard, Val(N)))
    obsubscript = isa(orbitals, NTuple{N, Int}) ? Subscript(orbitals) : isa(orbitals, Nothing) ? Subscript{N}() : orbitals
    spsubscript = isa(spins, NTuple{N, Int}) ? Subscript(spins) : isa(spins, Nothing) ? Subscript{N}() : spins
    orbitals = obsubscript.opattern
    spins = spsubscript.opattern
    (nambus === nothing) && (nambus = ntuple(i->wildcard, Val(N)))
    obsubs = ntuple(i->obsubscript.identifier, Val(N))
    spsubs = ntuple(i->spsubscript.identifier, Val(N))
    return FockCoupling(value, ID(FCID, centers, atoms, orbitals, spins, nambus, obsubs, spsubs), Subscripts(obsubscript), Subscripts(spsubscript))
end

"""
    show(io::IO, fc::FockCoupling)

Show a Fock coupling.
"""
function Base.show(io::IO, fc::FockCoupling)
    @printf io "FockCoupling{%s}(value=%s" rank(fc) decimaltostr(fc.value)
    cache = []
    for attrname in (:centers, :atoms, :orbitals, :spins, :nambus)
        attrvalue = getproperty(fc.id, attrname)
        any(attrvalue .≠ wildcard) && @printf io ", %s=(%s)" attrname join(attrvalue, ", ")
    end
    for attrname in (:obsubs, :spsubs)
        attrvalue = getproperty(fc.id, attrname)
        all(attrvalue .== wildcard) || all(attrvalue .== constant) || @printf io ", %s=(%s)" attrname join(attrvalue, ", ")
    end
    @printf io ")"
end

"""
    repr(fc::FockCoupling) -> String

Get the repr representation of a Fock coupling.
"""
function Base.repr(fc::FockCoupling)
    cache = []
    flag = true
    for (attrname, abbr) in zip((:atoms, :orbitals, :spins, :nambus), ("sl", "ob", "sp", "ph"))
        attrvalue = getproperty(fc.id, attrname)
        any(attrvalue .≠ wildcard) && push!(cache, @sprintf "%s(%s)" abbr join(attrvalue, ", "))
    end
    result = decimaltostr(fc.value)
    length(cache) > 0 && (flag = false; result = @sprintf "%s %s" result join(cache, "⊗"))
    any((centers = fc.id.centers) .≠ wildcard) && (flag = false; result = @sprintf "%s%s@(%s)" result ((length(cache) > 0) ? "" : " ") join(centers, ", "))
    obsubs, spsubs = fc.id.obsubs, fc.id.spsubs
    ((all(obsubs .== wildcard) || all(obsubs .== constant)) && (all(spsubs .== wildcard) || all(spsubs .== constant))) || (
        flag = false;
        result = @sprintf "%s with %s" result join(((@sprintf "(%s, %s)" obsub spsub) for (obsub, spsub) in zip(obsubs, spsubs)), " && ")
        )
    flag && (result = @sprintf "%s {%s}" result rank(fc))
    return result
end

"""
    *(fc1::FockCoupling, fc2::FockCoupling) -> FockCoupling

Get the multiplication between two Fock couplings.
"""
function Base.:*(fc1::FockCoupling, fc2::FockCoupling)
    FockCoupling(fc1.value*fc2.value, fc1.id*fc2.id, fc1.obsubscripts*fc2.obsubscripts, fc1.spsubscripts*fc2.spsubscripts)
end

"""
    ⊗(fc1::FockCoupling, fc2::FockCoupling) -> FockCoupling

Get the direct product between two Fock couplings.
"""
function ⊗(fc1::FockCoupling, fc2::FockCoupling)
    @assert fc1|>rank == fc2|>rank "⊗ error: dismatched rank."
    value, attrvalues, subscripts, N = fc1.value*fc2.value, [], [], rank(fc1)
    wildcards = NTuple{N, Char}(wildcard for i = 1:N)
    v1, v2 = fc1.id.centers, fc2.id.centers
    t1, t2 = (v1 == wildcards), (v2 == wildcards)
    push!(attrvalues, (t1 && t2) ? wildcards : t1 ? v2 : t2 ? v1 : (v1 == v2) ? v1 : error("⊗ error: dismatched centers ($v1 and $v2)."))
    for (i, attrname) in enumerate((:atoms, :orbitals, :spins, :nambus, :obsubs, :spsubs))
        v1, v2 = getproperty(fc1.id, attrname), getproperty(fc2.id, attrname)
        t1, t2 = (v1 == wildcards), (v2 == wildcards)
        if attrname ∉ (:obsubs, :spsubs)
            @assert t1 || t2 "⊗ error: dismatched $attrname ($v1 and $v2)."
        else
            @assert t1 || t2 "⊗ error: dismatched $attrname ($v1 and $v2)."
            (attrname == :obsubs) && push!(subscripts, t1 ? fc2.obsubscripts : fc1.obsubscripts)
            (attrname == :spsubs) && push!(subscripts, t1 ? fc2.spsubscripts : fc1.spsubscripts)
        end
        push!(attrvalues, t1 ? v2 : v1)
    end
    return FockCoupling(value, ID(FCID, attrvalues...), subscripts...)
end

"""
    ⋅(fc1::FockCoupling, fc2::FockCoupling) -> FockCoupling

Get the dot product of two rank-2 Fock couplings.

A rank-2 FockCoupling can be considered as a matrix acting on the sublattice, orbital, spin and nambu spaces.
The dot product here is defined as the multiplication between such matrices.
"""
function ⋅(fc1::FockCoupling, fc2::FockCoupling)
    @assert rank(fc1) == rank(fc2) == 2 "⋅ error: input fockcouplings must be rank-2."
    attrpairs = []
    value = fc1.value * fc2.value
    wildcards = (wildcard, wildcard)
    v1, v2 = fc1.id.centers, fc2.id.centers
    t1, t2 = (v1 == wildcards), (v2 == wildcards)
    push!(attrpairs, :centers=>((t1 && t2) ? nothing : t1 ? v2 : t2 ? v1 : (v1 == v2) ? v1 : error("⋅ error: dismatched centers ($v1 and $v2).")))
    for (i, attrname) in enumerate((:atoms, :orbitals, :spins, :nambus))
        v1, v2 = getproperty(fc1.id, attrname), getproperty(fc2.id, attrname)
        if isa(v1, NTuple{2, Int}) && isa(v2, NTuple{2, Int})
            value = value * delta(v1[2], v2[1])
            push!(attrpairs, attrname=>(v1[1], v2[2]))
        else
            @assert (v1 == wildcards) && (v2 == wildcards) "⋅ error: dismatched $attrname ($v1 and $v2)."
        end
    end
    return FockCoupling{2}(value; attrpairs...)
end

"""
    expand(fc::FockCoupling, pid::PID, fock::Fock, kind::Union{Val{K}, Nothing}=nothing) where K -> Union{FCExpand, Tuple{}}
    expand(fc::FockCoupling, pids::NTuple{R, PID}, focks::NTuple{R, Fock}, kind::Union{Val{K}, Nothing}=nothing) where {R, K} -> Union{FCExpand, Tuple{}}

Expand a Fock coupling with the given set of point ids and Fock degrees of freedom.
"""
expand(fc::FockCoupling, pid::PID, fock::Fock, kind::Union{Val{K}, Nothing}=nothing) where {K} = expand(fc, (pid,), (fock,), kind)
function expand(fc::FockCoupling, pids::NTuple{R, PID}, focks::NTuple{R, Fock}, kind::Union{Val{K}, Nothing}=nothing) where {R, K}
    centers = couplingcenters(typeof(fc), fc.id.centers, Val(R))
    rpids = NTuple{rank(fc), eltype(pids)}(pids[centers[i]] for i = 1:rank(fc))
    rfocks = NTuple{rank(fc), eltype(focks)}(focks[centers[i]] for i = 1:rank(fc))
    nambus = fockcouplingnambus(kind, fc.id.nambus, NTuple{rank(fc), Int}(rfocks[i].nnambu for i = 1:rank(fc)))
    for (i, atom) in enumerate(fc.id.atoms)
       isa(atom, Int) && (atom ≠ rfocks[i].atom) && return ()
    end
    obsbexpands = expand(fc.obsubscripts, NTuple{rank(fc), Int}(rfocks[i].norbital for i = 1:rank(fc))) |> collect
    spsbexpands = expand(fc.spsubscripts, NTuple{rank(fc), Int}(rfocks[i].nspin for i = 1:rank(fc))) |> collect
    return FCExpand(fc.value, rpids, obsbexpands, spsbexpands, nambus)
end
couplingcenter(::Type{<:FockCoupling}, i::Int, n::Int, ::Val{2}) = (i <= n/2) ? 1 : 2
@generated function fockcouplingnambus(::Union{Val{K}, Nothing}, nambus::NTuple{N, Any}, ranges::NTuple{N, Int}) where {N, K}
    exprs = [:(isa(nambus[$i], Int) ?
            (((ranges[$i] == 1) && (nambus[$i] == 0)) || ((ranges[$i] == 2) && (0 < nambus[$i] <= 2)) ? nambus[$i] : error("fockcouplingnambus error: nambu out of range.")) :
            ((ranges[$i] == 1) ? 0 : (($i)%2 == 1) ? CREATION : ANNIHILATION)
            ) for i = 1:N]
    return Expr(:tuple, exprs...)
end
struct FCExpand{V, N, S} <: VectorSpace{Tuple{V, NTuple{N, FIndex{S}}}}
    value::V
    pids::NTuple{N, PID{S}}
    obsbexpands::Vector{NTuple{N, Int}}
    spsbexpands::Vector{NTuple{N, Int}}
    nambus::NTuple{N, Int}
end
IsMultiIndexable(::Type{<:FCExpand}) = IsMultiIndexable(true)
MultiIndexOrderStyle(::Type{<:FCExpand}) = MultiIndexOrderStyle('C')
dims(fce::FCExpand) = (length(fce.obsbexpands), length(fce.spsbexpands))
@generated function Tuple(inds::NTuple{2, Int}, fce::FCExpand{V, N}) where {V, N}
    exprs = [:(FIndex(fce.pids[$i], FID(fce.obsbexpands[inds[1]][$i], fce.spsbexpands[inds[2]][$i], fce.nambus[$i]))) for i = 1:N]
    return Expr(:tuple, :(fce.value), Expr(:tuple, exprs...))
end

"""
    σ⁰(mode::String; centers::Union{NTuple{2, Int}, Nothing}=nothing) -> Couplings{ID{FCID, 2}, FockCoupling{Int, <:ID{FCID, 2}}}

The Pauli matrix σ⁰, which can act on the space of spins("sp"), orbitals("ob"), sublattices("sl") or particle-holes("ph").
"""
function σ⁰(mode::String; centers::Union{NTuple{2, Int}, Nothing}=nothing)
    @assert (mode == "sp") || (mode == "ob") || (mode == "sl") || (mode == "ph") "σ⁰ error: not supported mode($mode)."
    attrname = (mode == "sp") ? :spins : (mode == "ob") ? :orbitals : (mode == "sl") ? :atoms : :nambus
    (attrval1, attrval2) = (mode == "ph") ? ((ANNIHILATION, CREATION), (CREATION, ANNIHILATION)) : ((1, 1), (2, 2))
    FockCoupling{2}(1; attrname=>attrval1, :centers=>centers) + FockCoupling{2}(1; attrname=>attrval2, :centers=>centers)
end

"""
    σˣ(mode::String; centers::Union{NTuple{2, Int}, Nothing}=nothing) -> Couplings{ID{FCID, 2}, FockCoupling{Int, <:ID{FCID, 2}}}

The Pauli matrix σˣ, which can act on the space of spins("sp"), orbitals("ob"), sublattices("sl") or particle-holes("ph").
"""
function σˣ(mode::String; centers::Union{NTuple{2, Int}, Nothing} = nothing)
    @assert (mode == "sp") || (mode == "ob") || (mode == "sl") || (mode == "ph") "σˣ error: not supported mode($mode)."
    attrname = (mode == "sp") ? :spins : (mode == "ob") ? :orbitals : (mode == "sl") ? :atoms : :nambus
    (attrval1, attrval2) = (mode == "ph") ? ((ANNIHILATION, ANNIHILATION), (CREATION, CREATION)) : ((1, 2), (2, 1))
    FockCoupling{2}(1; attrname=>attrval1, :centers=>centers) + FockCoupling{2}(1; attrname=>attrval2, :centers=>centers)
end

"""
    σʸ(mode::String; centers::Union{NTuple{2, Int}, Nothing}=nothing) -> Couplings{ID{FCID, 2}, FockCoupling{Int, <:ID{FCID, 2}}}

The Pauli matrix σʸ, which can act on the space of spins("sp"), orbitals("ob"), sublattices("sl") or particle-holes("ph").
"""
function σʸ(mode::String; centers::Union{NTuple{2, Int}, Nothing} = nothing)
    @assert (mode == "sp") || (mode == "ob") || (mode == "sl") || (mode == "ph") "σʸ error: not supported mode($mode)."
    attrname = (mode == "sp") ? :spins : (mode == "ob") ? :orbitals : (mode == "sl") ? :atoms : :nambus
    (attrval1, attrval2) = (mode == "ph") ? ((ANNIHILATION, ANNIHILATION), (CREATION, CREATION)) : ((1, 2), (2, 1))
    FockCoupling{2}(1im; attrname=>attrval1, :centers=>centers) + FockCoupling{2}(-1im; attrname=>attrval2, :centers=>centers)
end

"""
    σᶻ(mode::String; centers::Union{NTuple{2, Int}, Nothing}=nothing) -> Couplings{ID{FCID, 2}, FockCoupling{Int, <:ID{FCID, 2}}}

The Pauli matrix σᶻ, which can act on the space of spins("sp"), orbitals("ob"), sublattices("sl") or particle-holes("ph").
"""
function σᶻ(mode::String; centers::Union{NTuple{2, Int}, Nothing} = nothing)
    @assert (mode == "sp") || (mode == "ob") || (mode == "sl") || (mode == "ph") "σᶻ error: not supported mode($mode)."
    attrname = (mode == "sp") ? :spins : (mode == "ob") ? :orbitals : (mode == "sl") ? :atoms : :nambus
    (attrval1, attrval2) = (mode == "ph") ? ((ANNIHILATION, CREATION), (CREATION, ANNIHILATION)) : ((1, 1), (2, 2))
    FockCoupling{2}(-1; attrname=>attrval1, :centers=>centers) + FockCoupling{2}(1; attrname=>attrval2, :centers=>centers)
end

"""
    σ⁺(mode::String; centers::Union{NTuple{2, Int}, Nothing}=nothing) -> Couplings{ID{FCID, 2}, FockCoupling{Int, <:ID{FCID, 2}}}

The Pauli matrix σ⁺, which can act on the space of spins("sp"), orbitals("ob"), sublattices("sl") or particle-holes("ph").
"""
function σ⁺(mode::String; centers::Union{NTuple{2, Int}, Nothing} = nothing)
    @assert (mode == "sp") || (mode == "ob") || (mode == "sl") || (mode == "ph") "σ⁺ error: not supported mode($mode)."
    attrname = (mode == "sp") ? :spins : (mode == "ob") ? :orbitals : (mode == "sl") ? :atoms : :nambus
    attrval = (mode == "ph") ? (CREATION, CREATION) : (2, 1)
    Couplings(FockCoupling{2}(1; attrname=>attrval, :centers=>centers))
end

"""
    σ⁻(mode::String; centers::Union{NTuple{2, Int}, Nothing}=nothing) -> Couplings{ID{FCID, 2}, FockCoupling{Int, <:ID{FCID, 2}}}

The Pauli matrix σ⁻, which can act on the space of spins("sp"), orbitals("ob"), sublattices("sl") or particle-holes("ph").
"""
function σ⁻(mode::String; centers::Union{NTuple{2, Int}, Nothing} = nothing)
    @assert (mode == "sp") || (mode == "ob") || (mode == "sl") || (mode == "ph") "σ⁻ error: not supported mode($mode)."
    attrname = (mode == "sp") ? :spins : (mode == "ob") ? :orbitals : (mode == "sl") ? :atoms : :nambus
    attrval = (mode == "ph") ? (ANNIHILATION, ANNIHILATION) : (1, 2)
    Couplings(FockCoupling{2}(1; attrname=>attrval, :centers=>centers))
end

"""
    fc"..." -> FockCoupling

Construct a FockCoupling from a literal string.
"""
macro fc_str(str)
    ps = split(str, " with ")
    conditions = (length(ps) == 2) ? fcconditions(ps[2]) : nothing
    pos=findfirst(x->x∈('@', '{', 's', 'o', 'p'), ps[1])
    @assert pos<length(ps[1]) "@fc_str error: wrong input pattern."
    ps = (ps[1][1:pos-1], replace(ps[1][pos:end], " "=>""))
    coeff = eval(Meta.parse(ps[1]))
    if (ps[2][1] == '{') && (ps[2][end] == '}')
        N = parse(Int, ps[2][end-1])
        return FockCoupling{N}(coeff)
    else
        attrpairs = []
        components = split(ps[2], '@')
        centers = (length(components) == 2) ? couplingcenters(components[2]) : nothing
        push!(attrpairs, :centers=>centers)
        N = (centers === nothing) ? nothing : length(centers)
        count = 0
        if length(components[1]) > 0
            for component in split(components[1], '⊗')
                attrname, attrvalue = fccomponent(component)
                (N === nothing) && (N = length(attrvalue))
                if isa(attrvalue, Expr)
                    @assert (attrname == :orbitals) || (attrname == :spins) "@fc_str error: wrong input pattern."
                    @assert N == length(attrvalue.args) "@fc_str error: dismatched ranks."
                    count += 1
                    condition = (conditions === nothing) ? :nothing : conditions[count]
                    push!(attrpairs, attrname=>eval(subscriptexpr(attrvalue, condition)))
                else
                    @assert N == length(attrvalue) "@fc_str error: dismatched ranks."
                    push!(attrpairs, attrname=>attrvalue)
                end
            end
        end
        return FockCoupling{N}(coeff; attrpairs...)
    end
end
function fcconditions(str::AbstractString)
    conditions = Meta.parse(str)
    conditions = (conditions.head == :tuple) ? conditions.args : [conditions]
    return [(condition == "*") ? :nothing : condition for condition in conditions]
end
function fccomponent(str::AbstractString)
    @assert (str[3] == '(') && (str[end] == ')') "fccomponent error: wrong input pattern."
    attrname = (str[1:2] == "sl") ? :atoms : (str[1:2] == "ob") ? :orbitals : (str[1:2] == "sp") ? :spins : (str[1:2] == "ph") ? :nambus : error("fccomponent error: wrong input pattern.")
    expr = Meta.parse(str[3:end])
    attrvalue = isa(expr, Expr) ? (all(isa(arg, Int) for arg in expr.args) ? Tuple(expr.args) : expr) : (expr,)
    return attrname => attrvalue
end

σᵅsplit(str::AbstractString) = (ps = split(replace(str, " "=>""), '@');  (length(ps) == 1) ? (ps[1], nothing) : (length(ps) == 2) ? (ps[1], couplingcenters(ps[2])) : "σᵅsplit error: wrong input pattern.")
σᵅname(mode::AbstractString) = (mode == "sp") ? :spins : (mode == "ob") ? :orbitals : (mode == "sl") ? :atoms : (mode == "ph") ? :nambus : error("σᵅname error: wrong input mode.")

"""
    σ⁰"sp"/σ⁰"sp@(c₁, c₂)" -> Couplings
    σ⁰"ob"/σ⁰"ob@(c₁, c₂)" -> Couplings
    σ⁰"sl"/σ⁰"sl@(c₁, c₂)" -> Couplings
    σ⁰"ph"/σ⁰"ph@(c₁, c₂)" -> Couplings

The Pauli matrix σ⁰, which can act on the space of spins("sp"), orbitals("ob"), sublattices("sl") or particle-holes("ph").
"""
macro σ⁰_str(str::String)
    mode, centers = σᵅsplit(str)
    attrname = σᵅname(mode)
    (attrval1, attrval2) = (mode == "ph") ? ((ANNIHILATION, CREATION), (CREATION, ANNIHILATION)) : ((1, 1), (2, 2))
    FockCoupling{2}(1; attrname=>attrval1, :centers=>centers) + FockCoupling{2}(1; attrname=>attrval2, :centers=>centers)
end

"""
    σˣ"sp"/σˣ"sp@(c₁, c₂)" -> Couplings
    σˣ"ob"/σˣ"ob@(c₁, c₂)" -> Couplings
    σˣ"sl"/σˣ"sl@(c₁, c₂)" -> Couplings
    σˣ"ph"/σˣ"ph@(c₁, c₂)" -> Couplings

The Pauli matrix σˣ, which can act on the space of spins("sp"), orbitals("ob"), sublattices("sl") or particle-holes("ph").
"""
macro σˣ_str(str::String)
    mode, centers = σᵅsplit(str)
    attrname = σᵅname(mode)
    (attrval1, attrval2) = (mode == "ph") ? ((ANNIHILATION, ANNIHILATION), (CREATION, CREATION)) : ((1, 2), (2, 1))
    FockCoupling{2}(1; attrname=>attrval1, :centers=>centers) + FockCoupling{2}(1; attrname=>attrval2, :centers=>centers)
end

"""
    σʸ"sp"/σʸ"sp@(c₁, c₂)" -> Couplings
    σʸ"ob"/σʸ"ob@(c₁, c₂)" -> Couplings
    σʸ"sl"/σʸ"sl@(c₁, c₂)" -> Couplings
    σʸ"ph"/σʸ"ph@(c₁, c₂)" -> Couplings

The Pauli matrix σʸ, which can act on the space of spins("sp"), orbitals("ob"), sublattices("sl") or particle-holes("ph").
"""
macro σʸ_str(str::String)
    mode, centers = σᵅsplit(str)
    attrname = σᵅname(mode)
    (attrval1, attrval2) = (mode == "ph") ? ((ANNIHILATION, ANNIHILATION), (CREATION, CREATION)) : ((1, 2), (2, 1))
    FockCoupling{2}(1im; attrname=>attrval1, :centers=>centers) + FockCoupling{2}(-1im; attrname=>attrval2, :centers=>centers)
end

"""
    σᶻ"sp"/σᶻ"sp@(c₁, c₂)" -> Couplings
    σᶻ"ob"/σᶻ"ob@(c₁, c₂)" -> Couplings
    σᶻ"sl"/σᶻ"sl@(c₁, c₂)" -> Couplings
    σᶻ"ph"/σᶻ"ph@(c₁, c₂)" -> Couplings

The Pauli matrix σᶻ, which can act on the space of spins("sp"), orbitals("ob"), sublattices("sl") or particle-holes("ph").
"""
macro σᶻ_str(str::String)
    mode, centers = σᵅsplit(str)
    attrname = σᵅname(mode)
    (attrval1, attrval2) = (mode == "ph") ? ((ANNIHILATION, CREATION), (CREATION, ANNIHILATION)) : ((1, 1), (2, 2))
    FockCoupling{2}(-1; attrname=>attrval1, :centers=>centers) + FockCoupling{2}(1; attrname=>attrval2, :centers=>centers)
end

"""
    σ⁺"sp"/σ⁺"sp@(c₁, c₂)" -> Couplings
    σ⁺"ob"/σ⁺"ob@(c₁, c₂)" -> Couplings
    σ⁺"sl"/σ⁺"sl@(c₁, c₂)" -> Couplings
    σ⁺"ph"/σ⁺"ph@(c₁, c₂)" -> Couplings

The Pauli matrix σ⁺, which can act on the space of spins("sp"), orbitals("ob"), sublattices("sl") or particle-holes("ph").
"""
macro σ⁺_str(str::String)
    mode, centers = σᵅsplit(str)
    attrname = σᵅname(mode)
    attrval = (mode == "ph") ? (CREATION, CREATION) : (2, 1)
    Couplings(FockCoupling{2}(1; attrname=>attrval, :centers=>centers))
end

"""
    σ⁻"sp"/σ⁻"sp@(c₁, c₂)" -> Couplings
    σ⁻"ob"/σ⁻"ob@(c₁, c₂)" -> Couplings
    σ⁻"sl"/σ⁻"sl@(c₁, c₂)" -> Couplings
    σ⁻"ph"/σ⁻"ph@(c₁, c₂)" -> Couplings

The Pauli matrix σ⁻, which can act on the space of spins("sp"), orbitals("ob"), sublattices("sl") or particle-holes("ph").
"""
macro σ⁻_str(str::String)
    mode, centers = σᵅsplit(str)
    attrname = σᵅname(mode)
    attrval = (mode == "ph") ? (ANNIHILATION, ANNIHILATION) : (1, 2)
    Couplings(FockCoupling{2}(1; attrname=>attrval, :centers=>centers))
end

"""
    Onsite{ST}( id::Symbol, value::Any;
                couplings::Union{Function, Coupling, Couplings, Nothing}=nothing,
                amplitude::Union{Function, Nothing}=nothing,
                modulate::Union{Function, Bool}=false,
                ) where {ST}

Onsite term.

Type alias for `Term{statistics, :Onsite, 2, id, V, Int, <:TermCouplings, <:TermAmplitude, <:Union{TermModulate, Nothing}}`.
"""
const Onsite{statistics, id, V, C<:TermCouplings, A<:TermAmplitude, M<:Union{TermModulate, Nothing}} = Term{statistics, :Onsite, 2, id, V, Int, C, A, M}
function Onsite{ST}(id::Symbol, value::Any;
                    couplings::Union{Function, Coupling, Couplings, Nothing}=nothing,
                    amplitude::Union{Function, Nothing}=nothing,
                    modulate::Union{Function, Bool}=false
                    ) where {ST}
    couplings = TermCouplings((couplings === nothing) ? FockCoupling{2}() : couplings)
    Term{ST, :Onsite, 2}(id, value, 0, couplings=couplings, amplitude=amplitude, modulate=modulate)
end
abbr(::Type{<:Onsite}) = :st
isHermitian(::Type{<:Onsite}) = nothing

"""
    Hopping{ST}(id::Symbol, value::Any, bondkind::Int=1;
                couplings::Union{Function, Coupling, Couplings, Nothing}=nothing,
                amplitude::Union{Function, Nothing}=nothing,
                modulate::Union{Function, Bool}=false,
                ) where {ST}

Hopping term.

Type alias for `Term{statistics, :Hopping, 2, id, V, Int, <:TermCouplings, <:TermAmplitude, <:Union{TermModulate, Nothing}}`.
"""
const Hopping{statistics, id, V, C<:TermCouplings, A<:TermAmplitude, M<:Union{TermModulate, Nothing}} = Term{statistics, :Hopping, 2, id, V, Int, C, A, M}
function Hopping{ST}(id::Symbol, value::Any, bondkind::Int=1;
                    couplings::Union{Function, Coupling, Couplings, Nothing}=nothing,
                    amplitude::Union{Function, Nothing}=nothing,
                    modulate::Union{Function, Bool}=false
                    ) where {ST}
    couplings = TermCouplings((couplings === nothing) ? FockCoupling{2}() : couplings)
    @assert bondkind ≠ 0 "Hopping error: input bondkind (neighbor) cannot be 0. Use `Onsite` instead."
    Term{ST, :Hopping, 2}(id, value, bondkind, couplings=couplings, amplitude=amplitude, modulate=modulate)
end
abbr(::Type{<:Hopping}) = :hp
isHermitian(::Type{<:Hopping}) = false

"""
    Pairing{ST}(id::Symbol, value::Any, bondkind::Int=0;
                couplings::Union{Function, Coupling, Couplings, Nothing}=nothing,
                amplitude::Union{Function, Nothing}=nothing,
                modulate::Union{Function, Bool}=false,
                ) where {ST}

Pairing term.

Type alias for `Term{statistics, :Pairing, 2, id, V, Int, <:TermCouplings, <:TermAmplitude, <:Union{TermModulate, Nothing}}`.
"""
const Pairing{statistics, id, V, C<:TermCouplings, A<:TermAmplitude, M<:Union{TermModulate, Nothing}} = Term{statistics, :Pairing, 2, id, V, Int, C, A, M}
function Pairing{ST}(id::Symbol, value::Any, bondkind::Int=0;
                    couplings::Union{Function, Coupling, Couplings, Nothing}=nothing,
                    amplitude::Union{Function, Nothing}=nothing,
                    modulate::Union{Function, Bool}=false
                    ) where {ST}
    couplings = TermCouplings((couplings === nothing) ? FockCoupling{2}() : couplings)
    Term{ST, :Pairing, 2}(id, value, bondkind, couplings=couplings, amplitude=amplitude, modulate=modulate)
end
abbr(::Type{<:Pairing}) = :pr
isHermitian(::Type{<:Pairing}) = false
function fockcouplingnambus(::Val{:Pairing}, nambus::NTuple{2, Any}, ranges::NTuple{2, Int})
    @assert ranges == (2, 2) "fockcouplingnambus error: ranges for Pairing terms must be (2, 2)."
    n1 = isa(nambus[1], Int) ? ((0 < nambus[1] <= 2) ? nambus[1] : error("fockcouplingnambus error: nambu out of range.")) : ANNIHILATION
    n2 = isa(nambus[2], Int) ? ((0 < nambus[2] <= 2) ? nambus[2] : error("fockcouplingnambus error: nambu out of range.")) : ANNIHILATION
    return (n1, n2)
end
function expand!(operators::Operators, term::Pairing, bond::AbstractBond, config::IDFConfig,
                table::Union{Table, Nothing}=nothing,
                half::Bool=false,
                coord::Union{Val{true}, Val{false}}=coordpresent)
    argtypes = Tuple{Operators, Term, AbstractBond, IDFConfig, Union{Table, Nothing}, Bool, Union{Val{true}, Val{false}}}
    invoke(expand!, argtypes, operators, term, bond, config, table, half, coord)
    isa(bond, Bond) && invoke(expand!, argtypes, operators, term, bond|>reverse, config, table, half, coord)
    return operators
end

const hubbard = FockCoupling{4}(spins=(2, 2, 1, 1), nambus=(CREATION, ANNIHILATION, CREATION, ANNIHILATION))
"""
    Hubbard{ST}(id::Symbol, value::Any; amplitude::Union{Function, Nothing}=nothing, modulate::Union{Function, Bool}=false) where ST

Hubbard term.

Type alias for `Term{statistics, :Hubbard, 4, id, V, Int, <:TermCouplings, <:TermAmplitude, <:Union{TermModulate, Nothing}}`.
"""
const Hubbard{statistics, id, V, C<:TermCouplings, A<:TermAmplitude, M<:Union{TermModulate, Nothing}} = Term{statistics, :Hubbard, 4, id, V, Int, C, A, M}
function Hubbard{ST}(id::Symbol, value::Any; amplitude::Union{Function, Nothing}=nothing, modulate::Union{Function, Bool}=false) where ST
    @assert value == value' "Hubbard error: only real values are allowed."
    Term{ST, :Hubbard, 4}(id, value, 0, couplings=hubbard, amplitude=amplitude, modulate=modulate)
end
abbr(::Type{<:Hubbard}) = :hb
isHermitian(::Type{<:Hubbard}) = true

const interorbitalinterspin = FockCoupling{4}(orbitals=(@subscript (α, α, β, β) with α < β), spins=(@subscript (σ₁, σ₁, σ₂, σ₂) with σ₁ ≠ σ₂), nambus=(CREATION, ANNIHILATION, CREATION, ANNIHILATION))
"""
    InterOrbitalInterSpin{ST}(id::Symbol, value::Any; amplitude::Union{Function, Nothing}=nothing, modulate::Union{Function, Bool}=false) where ST

Interorbital-interspin term.

Type alias for `Term{statistics, :InterOrbitalInterSpin, 4, id, V, Int, <:TermCouplings, <:TermAmplitude, <:Union{TermModulate, Nothing}}`.
"""
const InterOrbitalInterSpin{statistics, id, V, C<:TermCouplings, A<:TermAmplitude, M<:Union{TermModulate, Nothing}} = Term{statistics, :InterOrbitalInterSpin, 4, id, V, Int, C, A, M}
function InterOrbitalInterSpin{ST}( id::Symbol, value::Any; amplitude::Union{Function, Nothing}=nothing, modulate::Union{Function, Bool}=false) where ST
    @assert value == value' "InterOrbitalInterSpin error: only real values are allowed."
    Term{ST, :InterOrbitalInterSpin, 4}(id, value, 0, couplings=interorbitalinterspin, amplitude=amplitude, modulate=modulate)
end
abbr(::Type{<:InterOrbitalInterSpin}) = :nons
isHermitian(::Type{<:InterOrbitalInterSpin}) = true

const interorbitalintraspin = FockCoupling{4}(orbitals=(@subscript (α, α, β, β) with α < β), nambus=(CREATION, ANNIHILATION, CREATION, ANNIHILATION))
"""
    InterOrbitalIntraSpin{ST}(id::Symbol, value::Any; amplitude::Union{Function, Nothing}=nothing, modulate::Union{Function, Bool}=false) where ST

Interorbital-intraspin term.

Type alias for `Term{statistics, :InterOrbitalIntraSpin, 4, id, V, Int, <:TermCouplings, <:TermAmplitude, <:Union{TermModulate, Nothing}}`.
"""
const InterOrbitalIntraSpin{statistics, id, V, C<:TermCouplings, A<:TermAmplitude, M<:Union{TermModulate, Nothing}} = Term{statistics, :InterOrbitalIntraSpin, 4, id, V, Int, C, A, M}
function InterOrbitalIntraSpin{ST}( id::Symbol, value::Any; amplitude::Union{Function, Nothing}=nothing, modulate::Union{Function, Bool}=false) where ST
    @assert value == value' "InterOrbitalIntraSpin error: only real values are allowed."
    Term{ST, :InterOrbitalIntraSpin, 4}(id, value, 0, couplings=interorbitalintraspin, amplitude=amplitude, modulate=modulate)
end
abbr(::Type{<:InterOrbitalIntraSpin}) = :noes
isHermitian(::Type{<:InterOrbitalIntraSpin}) = true

const spinflip = FockCoupling{4}(orbitals=(@subscript (α, β, α, β) with α < β), spins=(2, 1, 1, 2), nambus=(CREATION, CREATION, ANNIHILATION, ANNIHILATION))
"""
    SpinFlip{ST}(id::Symbol, value::Any; amplitude::Union{Function, Nothing}=nothing, modulate::Union{Function, Bool}=false) where ST

Spin-flip term.

Type alias for `Term{statistics, :SpinFlip, 4, id, V, Int, <:TermCouplings, <:TermAmplitude, <:Union{TermModulate, Nothing}}`.
"""
const SpinFlip{statistics, id, V, C<:TermCouplings, A<:TermAmplitude, M<:Union{TermModulate, Nothing}} = Term{statistics, :SpinFlip, 4, id, V, Int, C, A, M}
function SpinFlip{ST}(id::Symbol, value::Any; amplitude::Union{Function, Nothing}=nothing, modulate::Union{Function, Bool}=false) where ST
    Term{ST, :SpinFlip, 4}(id, value, 0, couplings=spinflip, amplitude=amplitude, modulate=modulate)
end
abbr(::Type{<:SpinFlip}) = :sf
isHermitian(::Type{<:SpinFlip}) = false

const pairhopping = FockCoupling{4}(orbitals=(@subscript (α, α, β, β) with α < β), spins=(2, 1, 1, 2), nambus=(CREATION, CREATION, ANNIHILATION, ANNIHILATION))
"""
    PairHopping{ST}(id::Symbol, value::Any; amplitude::Union{Function, Nothing}=nothing, modulate::Union{Function, Bool}=false) where ST

Pair-hopping term.

Type alias for `Term{statistics, :PairHopping, 4, id, V, Int, <:TermCouplings, <:TermAmplitude, <:Union{TermModulate, Nothing}}`.
"""
const PairHopping{statistics, id, V, C<:TermCouplings, A<:TermAmplitude, M<:Union{TermModulate, Nothing}} = Term{statistics, :PairHopping, 4, id, V, Int, C, A, M}
function PairHopping{ST}(id::Symbol, value::Any; amplitude::Union{Function, Nothing}=nothing, modulate::Union{Function, Bool}=false) where ST
    Term{ST, :PairHopping, 4}(id, value, 0, couplings=pairhopping, amplitude=amplitude, modulate=modulate)
end
abbr(::Type{<:PairHopping}) = :ph
isHermitian(::Type{<:PairHopping}) = false

"""
    Coulomb{ST}(    id::Symbol, value::Any, bondkind::Int=1;
                    couplings::Union{Function, Coupling, Couplings, Nothing}=nothing,
                    amplitude::Union{Function, Nothing}=nothing,
                    modulate::Union{Function, Bool}=false
                    ) where ST

Coulomb term.

Type alias for `Term{statistics, :Coulomb, 4, id, V, Int, <:TermCouplings, <:TermAmplitude, <:Union{TermModulate, Nothing}}`.
"""
const Coulomb{statistics, id, V, C<:TermCouplings, A<:TermAmplitude, M<:Union{TermModulate, Nothing}} = Term{statistics, :Coulomb, 4, id, V, Int, C, A, M}
function Coulomb{ST}(id::Symbol, value::Any, bondkind::Int=1;
                    couplings::Union{Function, Coupling, Couplings, Nothing}=nothing,
                    amplitude::Union{Function, Nothing}=nothing,
                    modulate::Union{Function, Bool}=false
                    ) where ST
    couplings = TermCouplings((couplings === nothing) ? FockCoupling{2}()*FockCoupling{2}() : couplings)
    @assert bondkind ≠ 0 "Coulomb error: input bondkind (neighbor) cannot be 0. Use `Hubbard/InterOrbitalInterSpin/InterOrbitalIntraSpin/SpinFlip/PairHopping` instead."
    Term{ST, :Coulomb, 4}(id, value, bondkind, couplings=couplings, amplitude=amplitude, modulate=modulate)
end
abbr(::Type{<:Coulomb}) = :cl
isHermitian(::Type{<:Coulomb}) = nothing
termfactor(::Nothing, id::ID{OID, 4}, ::Val{:Coulomb}) = ((id[2]' == id[1]) && (id[4]' == id[3])) ? 2 : 1

"""
    FFockTerm

Fermionic Fock term types.
"""
const FFockTerm = Union{Onsite{'F'}, Hopping{'F'}, Pairing{'F'}, Hubbard{'F'}, InterOrbitalInterSpin{'F'}, InterOrbitalIntraSpin{'F'}, SpinFlip{'F'}, PairHopping{'F'}, Coulomb{'F'}}

"""
    BFockTerm

Bosonic Fock term types.
"""
const BFockTerm = Union{Onsite{'B'}, Hopping{'B'}, Pairing{'B'}, Hubbard{'B'}, InterOrbitalInterSpin{'B'}, InterOrbitalIntraSpin{'B'}, SpinFlip{'B'}, PairHopping{'B'}, Coulomb{'B'}}

"""
    otype(T::Type{<:FFockTerm}, I::Type{<:OID})
    otype(T::Type{<:BFockTerm}, I::Type{<:OID})

Get the operator type of a Fock term.
"""
otype(T::Type{<:FFockTerm}, I::Type{<:OID}) = FOperator{T|>valtype, ID{I, T|>rank}}
otype(T::Type{<:BFockTerm}, I::Type{<:OID}) = BOperator{T|>valtype, ID{I, T|>rank}}

end # module


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Desription
# ==============================================================================
#
#   Tests related to the kinematics of DCMs.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# File: ./src/dcm.jl
# ==================

# Functions: ddcm
# ---------------

@testset "Kinematics of DCMs (Float64)" begin
    T = Float64

    # Create a random DCM.
    Dba₀ = create_rotation_matrix(_rand_ang(), :Z) *
        create_rotation_matrix(_rand_ang(), :Y) *
        create_rotation_matrix(_rand_ang(), :X)

    # Create a random velocity vector.
    wba_a = @SVector randn(3)

    # Propagate the initial DCM using the sampled velocity vector.
    Δ   = 1e-7
    Dba = Dba₀
    num = 2_000_000

    for k in 1:num
        dDba = ddcm(Dba, Dba * wba_a)
        Dba  = Dba + dDba * Δ
        Dba  = orthonormalize(Dba)
    end

    # In the end, the vector aligned with `wba_a` must not change.
    v₀ = Dba₀ * wba_a
    v₁ = Dba * wba_a
    @test v₁ ≈ v₀ rtol = 1e-6
    @test eltype(v₀) === T
    @test eltype(v₁) === T

    # In the end, a vector perpendicular to `wba_a` must rotate the angle
    # compatible with the angular velocity and time of integration.
    aux   = @SVector randn(3)
    aux   = aux / norm(aux)
    vr_a  = wba_a × aux
    vr_a  = vr_a / norm(vr_a)
    vr_b₀ = Dba₀ * vr_a
    vr_b  = Dba  * vr_a
    θ     = acos(vr_b ⋅ vr_b₀)
    @test eltype(vr_b₀) === T
    @test eltype(vr_b) === T

    # Estimate θ based on the angular velocity.
    θest = mod(norm(wba_a) * Δ * num, 2π)
    θest > pi && (θest = 2π - θest)
    @test θ ≈ θest atol = 1e-6
end

@testset "Kinematics of DCMs (Float32)" begin
    T = Float32

    # Create a random DCM.
    Dba₀ = create_rotation_matrix(_rand_ang(T), :Z) *
        create_rotation_matrix(_rand_ang(T), :Y) *
        create_rotation_matrix(_rand_ang(T), :X)

    # Create a random velocity vector.
    wba_a = @SVector randn(T, 3)

    # Propagate the initial DCM using the sampled velocity vector.
    Δ   = 1f-7
    Dba = Dba₀
    num = 20_000

    for k in 1:num
        dDba = ddcm(Dba, Dba * wba_a)
        Dba  = Dba + dDba * Δ
        Dba  = orthonormalize(Dba)
    end

    # In the end, the vector aligned with `wba_a` must not change.
    v₀ = Dba₀ * wba_a
    v₁ = Dba * wba_a
    @test v₁ ≈ v₀ rtol = 1e-3
    @test eltype(v₀) === T
    @test eltype(v₁) === T

    # In the end, a vector perpendicular to `wba_a` must rotate the angle
    # compatible with the angular velocity and time of integration.
    aux   = @SVector randn(T, 3)
    aux   = aux / norm(aux)
    vr_a  = wba_a × aux
    vr_a  = vr_a / norm(vr_a)
    vr_b₀ = Dba₀ * vr_a
    vr_b  = Dba  * vr_a
    θ     = acos(vr_b ⋅ vr_b₀)
    @test eltype(vr_b₀) === T
    @test eltype(vr_b) === T

    # Estimate θ based on the angular velocity.
    θest = mod(norm(wba_a) * Δ * num, T(2π))
    θest > pi && (θest = T(2π) - θest)
    @test θ ≈ θest atol = 1e-3
end

@testset "Kinematics of DCMs (Errors)" begin
    @test_throws ArgumentError ddcm(DCM(I), [1, 2])
    @test_throws ArgumentError ddcm(DCM(I), [1, 2, 3, 4])
end