text:
using .GraphSignal, .GraphPartition, .BasisSpecification, LinearAlgebra

"""
    dmatrix = dvec2dmatrix(dvec::Matrix{Float64}, GP::GraphPart, BS::BasisSpec)

Given a vector of expansion coefficients, convert it to a matrix.

### Input Arguments
* `dvec::Matrix{Float64}`: a vector of expansion coefficients
* `GP::GraphPart`: a GraphPart object
* `BS::BasisSpec`: a BasisSpec object

### Output Argument
* `dmatrix::Array{Float64,3}`: a set of matrices of expansion coefficients
"""
function dvec2dmatrix(dvec::Matrix{Float64}, GP::GraphPart, BS::BasisSpec)

    #
    # 0. Preliminaries
    #
    # extract data
    # MATLAB: [levlist,levlengths] = ExtractData(BS,GP);
    #if isempty(BS.levlengths) && GP != nothing
    #    levlist2levlengths!(GP, BS)
    #end

    levlist = BS.levlist
    #levlengths = BS.levlengths

    # constants
    (N, jmax) = Base.size(GP.rs)
    N = N - 1
    fcols = Base.size(dvec, 2)

    # allocate space
    dmatrix = zeros(N, jmax, fcols)

    #
    # 1. Put the entries in the vector into the correct places in the matrix
    #
    # n = 1
    # for row = 1:length(levlist)
    #     dmatrix[n:(n + levlengths[row] - 1), levlist[row], :] =
    #         dvec[n:(n + levlengths[row] - 1), :]
    #     n += levlengths[row]
    # end
    for i = 1:length(levlist)
        dmatrix[levlist[i][1], levlist[i][2], :] = dvec[i,:]
    end
    #
    # 2. Return the `dmatrix` array
    #
    return dmatrix
end # of function dvec2dmatrix

"""
    (dvec, BS) = dmatrix2dvec(dmatrix::Array{Float64,3}, GP::GraphPart)

Given a matrix of expansion coefficients, convert it to a vector.

### Input Arguments
* `dmatrix::Array{Float64,3}`: matrices of expansion coefficients
* `GP::GraphPart`: an input GraphPart object
* `BS::BasisSpec`: an input BasisSpec object

### Outputs Arguments
* `dvec::Matrix{Float64}`: a vector of expansion coefficients
* `BS::BasisSpec`: an output BasisSpec object
"""
function dmatrix2dvec(dmatrix::Array{Float64,3}, GP::GraphPart)

    # constants
    (N, jmax, fcols) = Base.size(dmatrix)

    # THE BASIS IS NOT SPECIFIED ==> retain nonzero coeffs and specify the basis
    ## 0. Preliminaries
    # allocate/initialize
    dvec = dmatrix[:, jmax, :]
    levlist = jmax * ones(Int, N)

    ## 1. Make a vector of the nonzero basis entries
    for j = (jmax - 1):-1:1
        regioncount = countnz(GP.rs[:, j]) - 1
        for r = 1:regioncount
            indr = GP.rs[r, j]:(GP.rs[r + 1, j] - 1)
            if countnz(dmatrix[indr, j, :]) > 0
                dvec[indr, :] = dmatrix[indr, j, :]
                levlist[GP.rs[r, j]:(GP.rs[r + 1, j] - 1)] .= j
            end
        end
    end
    ## 2. Specify the corresponding basis and return things
    levlist = collect(enumerate(levlist))
    BS = BasisSpec(levlist)
    #levlist2levlengths!(GP, BS)
    return dvec, BS
end # of function dmatrix2dvec (with no BS input)

"""
    dvec = dmatrix2dvec(dmatrix::Array{Float64,3}, GP::GraphPart, BS::BasisSpec)

Given a matrix of expansion coefficients, convert it to a vector.
This function assumes that the input coefficient array `dmatrix` is in the coarse-to-fine format. If `BS.c2f == false`, then this function internally converts `dmatrix` into the fine-to-coarse format. Hence, if one supplies the f2c `dmatrix`, the results become wrong, and the subsequent procedure may result in error.

### Input Arguments
* `dmatrix::Array{Float64,3}`: matrices of expansion coefficients
* `GP::GraphPart`: an input GraphPart object
* `BS::BasisSpec`: an input BasisSpec object

### Outputs Arguments
* `dvec::Matrix{Float64}`: a vector of expansion coefficients
"""
function dmatrix2dvec(dmatrix::Array{Float64,3}, GP::GraphPart, BS::BasisSpec)

    # constants
    (N, jmax, fcols) = Base.size(dmatrix)

    # THE BASIS IS SPECIFIED ==> select the corresponding coefficients
    ## 0. Preliminaries
    # allocate space
    dvec = zeros(N, fcols)
    # if isempty(BS.levlengths)
    #     levlist2levlengths!(GP, BS)
    # else
    #     @warn("We assume that levlengths of the input BS were already computed properly prior to dmatrix2dvec.")
    # end
    levlist = BS.levlist
    #levlengths = BS.levlengths
    BSc2f = BS.c2f
    # put dmatrix in the fine-to-coarse arrangement, if necessary
    if !BSc2f # This assumes that dmatrix has not yet been arranged in f2c
    #if !BSc2f && isempty(GP.rsf2c)
        dmatrix = fine2coarse!(GP, dmatrix = dmatrix, coefp = true)
    end

    ## 1. Make a vector of the matrix entries specified by the BS object
    # n = 1
    # for row = 1:length(levlist)
    #     dvec[n:(n + levlengths[row] - 1), :] =
    #         dmatrix[n:(n + levlengths[row] - 1), levlist[row], :]
    #     n += levlengths[row]
    # end

    for i = 1:length(levlist)
        dvec[i,:] = dmatrix[levlist[i][1], levlist[i][2], :]
    end

    ## 2. Return the dvec
    return dvec
end # of dmatrix2dvec (with BS input)

"""
    levlist2levlengths!(GP::GraphPart, BS::BasisSpec)

Compute the levlengths info for a BasisSpec object

### Input Arguments
* `GP::GraphPart`: a GraphPart object
* `BS::BasisSpec`: BasisSpec object, without `levlengths`; after this function, the `levlengths` field is filled.
"""
function levlist2levlengths!(GP::GraphPart, BS::BasisSpec)

    # 0. Preliminaries
    levlist = BS.levlist
    rs = GP.rs
    # allocate space
    levlengths = zeros(typeof(rs[1]), length(levlist))

    # 1. Determine the length of each basis block specified by levlist
    if BS.c2f                   # coarse-to-fine case
        n = 0
        for row = 1:length(levlist)
            # find the regionstart(s) of the next region
            # MATLAB: IX = find( GP.rs(:,levlist(row))==n+1, 1, 'last' );
            IX = something(findlast(isequal(n+1), rs[:, levlist[row]]), 0)
            #IX = findlast( rs[:, levlist[row]], n + 1 )
            levlengths[row] = rs[IX+1, levlist[row]] - rs[IX, levlist[row]]
            n += levlengths[row]
        end
    else                        # fine-to-coarse case
        n = 0
        if isempty(GP.rsf2c)
            fine2coarse!(GP)
        end
        rsf2c = GP.rsf2c        # This must be here: GP.rsf2c cannot be null!
        for row = 1:length(levlist)
            IX = something(findlast(isequal(n+1), rsf2c[:, levlist[row]]), 0)
            #IX = findlast( rsf2c[:, levlist[row]], n + 1 )
            levlengths[row] = rsf2c[IX + 1, levlist[row]] - rsf2c[IX, levlist[row]]
            n += levlengths[row]
        end
    end

    # 2. Postprocessing
    # get rid of blocks with length 0
    # MATLAB: levlist(levlengths == 0) = [];
    #         levlengths(levlengths == 0) = [];
    #         BS = BasisSpec(levlist,levlengths,c2f,description);
    levlist = levlist[levlengths .!= 0] # . before != is important here!!
    levlengths = levlengths[levlengths .!= 0]
    BS.levlengths = levlengths

end # of function levlist2levlengths!

"""
    (levlistfull, levlengthsfull, transfull) = bsfull(GP::GraphPart, BS::BasisSpec, trans::Vector{Bool})

Given a BasisSpec object, return the full-length, redundant levlist, levlengths, and trans descriptions.

###  Input Arguments
* `GP::GraphPart`: an input GraphPart object
* `BS::BasisSpec`: an input BasisSpec object
* `trans::Vector{Bool}`: a specification of the transforms used for the HGLET-GHWT hybrid transform (default: null)
* `levlengthsp::Bool`: a flag to return levlengthsfull (default: false)
* `transp::Bool`: a flag to return transfull (default: false)

###  Output Arguments
* `levlistfull::Vector{UInt8}`: the full-length, redundant levels list description
* `levlengthsfull::Vector{UInt8}`: the full-length, redundant levels lengths description
* `transfull::Matrix{Bool}`: the full-length, redundant trans description
"""
function bsfull(GP::GraphPart, BS::BasisSpec;
                trans::Vector{Bool} = Vector{Bool}(undef, 0),
                levlengthsp::Bool = false, transp::Bool = false)

    ## 0. Preliminaries

    # extract data
    if isempty(BS.levlengths)
        levlist2levlengths!(GP, BS)
    end
    levlist = BS.levlist
    levlengths = BS.levlengths

    # allocate space
    N = Base.length(GP.ind)
    levlistfull = zeros(Int, N)
    # Assuming that the maximum value of levlist <= jmax \approx log2(N)
    # can be representable by `UInt8`, i.e., N < \approx 5.79E76, which is
    # reasonable. We cannot handle such a large N at this point.

    if levlengthsp
        levlengthsfull = zeros(Int, N)
    end

    if transp && !isempty(trans)
        cols = Base.size(trans, 2)
        # MATLAB: transfull = false(N, cols)
        transfull = falses(N, cols)
    else
        transfull = falses(0)   # transfull is empty.
    end

    ## 1. Fill out the redundant descriptions
    global idx = 0
    for row = 1:length(levlist)
        levlistfull[(idx + 1):(idx + levlengths[row])] .= levlist[row]
        if levlengthsp
            levlengthsfull[(idx + 1):(idx + levlengths[row])] .= levlengths[row]
        end
        if !isempty(transfull)
            transfull[(idx + 1):(idx + levlengths[row]), :] = repmat(trans[row, :], levlengths[row], 1)
        end
        global idx += levlengths[row]
    end

    ## 2. Prepare the returns
    if levlengthsp
        if transp
            return levlistfull, levlengthsfull, transfull
        else
            return levlistfull, levlengthsfull
        end
    else
        if transp
            return levlistfull, transfull
        else
            return levlistfull
        end
    end
end # of function bsfull

"""
    BS = bs_haar(GP::GraphPart)

Specify the Haar basis for a given graph partitioning

### Input Argument
* `GP::GraphPart`: an input GraphPart object

### Output Argument
* `BS::BasisSpec`: a BasisSpec object corresponding to the Haar basis
"""
# function bs_haar(GP::GraphPart)
#
#     # determine jmax
#     jmax = Base.size(GP.rs, 2)
#
#     # allocate space for levlist
#     levlist = zeros(Int, jmax)
#
#     # fill in levlist for the Haar basis
#     levlist[1] = jmax
#     for row = 2:jmax
#         levlist[row] = jmax + 2 - row
#     end
#
#     # make a BasisSpec object
#     BS = BasisSpec(levlist, c2f = false, description = "Haar basis")
#     # fill in the levlengths field of the BasisSpec object
#     levlist2levlengths!(GP, BS)
#     # return it
#     return BS
# end # of function bs_haar

function bs_haar(GP::GraphPart)
    if isempty(GP.tagf2c)
        fine2coarse!(GP)
    end
    # determine jmax
    jmax = Base.size(GP.rs, 2)

    # allocate space for levlist
    levlist = Vector{Tuple{Int,Int}}((findall(x->x==1, GP.tagf2c)))
    pushfirst!(levlist, (1, jmax))

    # make a BasisSpec object
    BS = BasisSpec(levlist, c2f = false, description = "Haar basis")
    # fill in the levlengths field of the BasisSpec object
    #levlist2levlengths!(GP, BS)
    # return it
    return BS
end # of function bs_haar

"""
    BS = bs_level(GP::GraphPart, j::Int, c2f::Bool = true)

Specify the basis corresponding to level j for a given graph partitioning

### Input Arguments
* `GP::GraphPart`: an input GraphPart object
* `j::Int`: the level to which the basis corresponds (`j = 0` is the global level)
* `c2f::Bool`: a flag for c2f or f2c (default: true, i.e., c2f)

### Output Argument
* `BS::BasisSpec`: a BasisSpec object corresponding to the level `j` basis
"""
function bs_level(GP::GraphPart, j::Int; c2f::Bool = true)

    # coarse-to-fine dictionary
    if c2f
        #rspointer = GP.rs[:, j + 1]
        bspec = "coarse-to-fine level $(j)"
    # fine-to-coarse dictionary
    else
        if isempty(GP.rsf2c)
            # if isempty(GP.tag) || isempty(GP.compinfo) # I don't think this part is necessary since GP should be partially filled.
            #    ghwt_core!(GP)
            # end
            fine2coarse!(GP)
        end
        #rspointer = GP.rsf2c[:, j + 1]
        bspec = "fine-to-coarse level $(j)"
    end

    # specify the level j basis
    #Nj = countnz(rspointer) - 1
    if isempty(GP.tag)          # otherwise, it doesn't work for HGLET!
        N = size(GP.rs)[1]-1
    else
        N = size(GP.tag, 1)
    end
    levlist = collect(enumerate((j + 1) * ones(Int,N)))
    BS = BasisSpec(levlist, c2f = c2f, description = bspec)

    # fill in the levlengths field of the BasisSpec object
    #levlist2levlengths!(GP, BS)

    # return it
    return BS
end # of bs_level

"""
    BS = bs_walsh(GP::GraphPart)

Specify the walsh basis corresponding to for a given graph partitioning

### Input Arguments
* `GP::GraphPart`: an input GraphPart object

### Output Argument
* `BS::BasisSpec`: a BasisSpec object corresponding to the walsh basis
"""
function bs_walsh(GP::GraphPart)
    BS = bs_level(GP, 0, c2f = true)
    BS.description = "Walsh basis"
    return BS
end # of bs_level

"""
function dvec_Threshold(dvec::Array{Float64}, SORH::String, keep::Float64,
                GP::GraphPart, BS::BasisSpec)

 Threshold HGLET / GHWT coefficients

### Input Arguments
*   `dvec::Array{Float64,2}`        the vector of expansion coefficients
*   `SORH::String`        use soft ('s') or hard ('h') thresholding
*   `keep::Float64`        a fraction between 0 and 1 which says how many coefficients should be kept
*   `GP::GraphPart`          a GraphPart object, used to identify scaling coefficients
*   `BS::BasisSpec`          a BasisSpec object, used to identify scaling coefficients

### Output Argument
*   `dvec_new::Vector{Float64}`        the thresholded expansion coefficients
*   `kept::Float64`        the number of coefficients kept
"""
function dvec_Threshold(dvec::Array{Float64,2}, SORH::String, keep::Float64,
                GP::GraphPart, BS::BasisSpec)
#
# Apply to 1-D only for now. Need to modify in the future
#
    dvec_new = deepcopy(dvec)
    # Obtain the number of coefficients kept
    if keep > 1 || keep < 0
      error("The input keep should between 0~1")
    else
      kept = Int64(round(keep*size(dvec,1)))
    end

    if SORH == "s" || SORH == "soft"
      ind = sortperm(abs.(dvec[:]), rev = true)
      T = abs(dvec[ind[kept+1]])

      if !isnothing(BS)
        #convert the tag matrix to a vector corresponding to BS
        tag = dmatrix2dvec(Array{Float64,3}(reshape(GP.tag,(size(GP.tag,1),size(GP.tag,2),1))),GP,BS)

        #soft-threshold the non-scaling coefficients
        ind = 1:size(dvec,1)
        ind = ind[(tag .> 0)[:]]
        dvec_new[ind] = sign.(dvec[ind]).*(abs.(dvec[ind]).-T).*(abs.(dvec[ind]).-T .>0)

      # if BS is not given, soft-threshold all coefficients
      else
        dvec_new = sign.(dvec).*(abs.(dvec).-T).*(abs.(dvec).-T.>0)
      end

      kept = count(!iszero, dvec_new)

    elseif SORH == "h" || SORH == "hard"
      ind = sortperm(abs.(dvec[:]), rev = true)
      dvec_new[ind[kept+1:end]] .= 0

      kept = count(!iszero, dvec_new)

    else

      kept = count(!iszero, dvec_new)

    end
    return dvec_new, kept
end


text:
# TODO: there should be a distinction between an object that will manage a full plot, vs a component of a plot.
# the PlotRecipe as currently implemented is more of a "custom component"
# a recipe should fully describe the plotting command(s) and call them, likewise for updating.
#   actually... maybe those should explicitly derive from PlottingObject???

abstract PlotRecipe

getRecipeXY(recipe::PlotRecipe) = Float64[], Float64[]
getRecipeArgs(recipe::PlotRecipe) = ()

plot(recipe::PlotRecipe, args...; kw...) = plot(getRecipeXY(recipe)..., args...; getRecipeArgs(recipe)..., kw...)
plot!(recipe::PlotRecipe, args...; kw...) = plot!(getRecipeXY(recipe)..., args...; getRecipeArgs(recipe)..., kw...)
plot!(plt::Plot, recipe::PlotRecipe, args...; kw...) = plot!(getRecipeXY(recipe)..., args...; getRecipeArgs(recipe)..., kw...)

num_series(x::AMat) = size(x,2)
num_series(x) = 1

_apply_recipe(d::Dict; kw...) = ()

# if it's not a recipe, just do nothing and return the args
function _apply_recipe(d::Dict, args...; issubplot=false, kw...)
    if issubplot && !haskey(d, :n) && !haskey(d, :layout)
        # put in a sensible default
        d[:n] = maximum(map(num_series, args))
    end
    args
end

# # -------------------------------------------------

# function rotate(x::Real, y::Real, θ::Real; center = (0,0))
#   cx = x - center[1]
#   cy = y - center[2]
#   xrot = cx * cos(θ) - cy * sin(θ)
#   yrot = cy * cos(θ) + cx * sin(θ)
#   xrot + center[1], yrot + center[2]
# end

# # -------------------------------------------------

# type EllipseRecipe <: PlotRecipe
#   w::Float64
#   h::Float64
#   x::Float64
#   y::Float64
#   θ::Float64
# end
# EllipseRecipe(w,h,x,y) = EllipseRecipe(w,h,x,y,0)

# # return x,y coords of a rotated ellipse, centered at the origin
# function rotatedEllipse(w, h, x, y, θ, rotθ)
#   # # coord before rotation
#   xpre = w * cos(θ)
#   ypre = h * sin(θ)

#   # rotate and translate
#   r = rotate(xpre, ypre, rotθ)
#   x + r[1], y + r[2]
# end

# function getRecipeXY(ep::EllipseRecipe)
#   x, y = unzip([rotatedEllipse(ep.w, ep.h, ep.x, ep.y, u, ep.θ) for u in linspace(0,2π,100)])
#   top = rotate(0, ep.h, ep.θ)
#   right = rotate(ep.w, 0, ep.θ)
#   linex = Float64[top[1], 0, right[1]] + ep.x
#   liney = Float64[top[2], 0, right[2]] + ep.y
#   Any[x, linex], Any[y, liney]
# end

# function getRecipeArgs(ep::EllipseRecipe)
#   [(:line, (3, [:dot :solid], [:red :blue], :path))]
# end

# # -------------------------------------------------

# "Correlation scatter matrix"
# function corrplot{T<:Real,S<:Real}(mat::AMat{T}, corrmat::AMat{S} = cor(mat);
#                                    colors = :redsblues,
#                                    labels = nothing, kw...)
#   m = size(mat,2)
#   centers = Float64[mean(extrema(mat[:,i])) for i in 1:m]

#   # might be a mistake?
#   @assert m <= 20
#   @assert size(corrmat) == (m,m)

#   # create a subplot grid, and a gradient from -1 to 1
#   p = subplot(rand(0,m^2); n=m^2, leg=false, grid=false, kw...)
#   cgrad = ColorGradient(colors, [-1,1])

#   # make all the plots
#   for i in 1:m
#     for j in 1:m
#       idx = p.layout[i,j]
#       plt = p.plts[idx]
#       if i==j
#         # histogram on diagonal
#         histogram!(plt, mat[:,i], c=:black)
#         i > 1 && plot!(plt, yticks = :none)
#       elseif i < j
#         # annotate correlation value in upper triangle
#         mi, mj = centers[i], centers[j]
#         plot!(plt, [mj], [mi],
#                    ann = (mj, mi, text(@sprintf("Corr:\n%0.3f", corrmat[i,j]), 15)),
#                    yticks=:none)
#       else
#         # scatter plots in lower triangle; color determined by correlation
#         c = RGBA(RGB(getColorZ(cgrad, corrmat[i,j])), 0.3)
#         scatter!(plt, mat[:,j], mat[:,i], w=0, ms=3, c=c, smooth=true)
#       end

#       if labels != nothing && length(labels) >= m
#         i == m && xlabel!(plt, string(labels[j]))
#         j == 1 && ylabel!(plt, string(labels[i]))
#       end
#     end
#   end

#   # link the axes
#   subplot!(p, link = (r,c) -> (true, r!=c))
# end

"Sparsity plot... heatmap of non-zero values of a matrix"
function spy{T<:Real}(y::AMat{T}; kw...)
  I,J,V = findnz(y)
  heatmap(J, I; leg=false, yflip=true, nbins=size(y), kw...)
end

ext: jl


text:
function gcam_go_techspecific(x, y)
    x * y
end


text:
module SelectiveInference

using LinearAlgebra
using IntervalSets
using Distributions
using Parameters
using Reexport

include("stepwise.jl")
@reexport using .Stepwise
include("../quadratic.jl")
using .Quadratic
include("../intervals.jl")
using .Intersection
include("../selective_p.jl")
using .MultiTruncatedDistributions
include("../ci.jl")
using .ConfidenceInterval

export SelectiveHypothesisTest, selective_p, selective_CI

const ε = 1e-4

struct TestStatistic{A,B,C,D}
    z_obs::A
    σ²_z::B
    𝐚::C
    𝐛::D
end
function parametric_representation(𝛈, 𝐲, 𝚺)
    z = 𝛈' * 𝐲
    σ²_z = 𝛈' * 𝚺 * 𝛈
    𝐛 = 𝚺 * 𝛈 / σ²_z
    𝐚 = 𝐲 - 𝐛 * 𝛈' * 𝐲
    @assert 𝐚 + 𝐛 * z ≈ 𝐲
    return TestStatistic(z, σ²_z, 𝐚, 𝐛)
end
"Make the hypothesis test problem: βⱼ = 0 or not for all j ∈ `alg_result.model.M`, where βⱼ is the coefficient of the j-th feature."
struct SelectiveHypothesisTest{A <: AbstractAlgorithm,B <: AbstractAlgResult,C,D,E,F}
    alg::A # used feature selection algorithm
    alg_result::B # the result of the algorithm `alg` for `𝐗`, `𝐲_obs`, `𝚺`
    test_stats::C # the information of test statistic for all hypotheses
    𝐗::D # used covariate matrix
    𝐲_obs::E # used response vector
    𝚺::F # used covariance matrix for response vector

    function SelectiveHypothesisTest(alg::S, alg_result::T, 𝐗::U, 𝐲_obs::V, 𝚺::W) where {S,T,U,V,W}
        M_obs = alg_result.model.M
        𝐗_obs = 𝐗[:, M_obs]
        𝚮 = 𝐗_obs / (𝐗_obs' * 𝐗_obs)
        test_stats = [parametric_representation(𝛈, 𝐲_obs, 𝚺) for 𝛈 in eachcol(𝚮)]
        new{S,T,typeof(test_stats),U,V,W}(alg, alg_result, test_stats, 𝐗, 𝐲_obs, 𝚺)
    end
end

function (::AIC)(𝐗, 𝚺⁻¹, test_stat::TestStatistic)
    𝐚, 𝐛 = test_stat.𝐚, test_stat.𝐛
    𝐀 = 𝚺⁻¹ - 𝚺⁻¹ * 𝐗 * pinv(𝐗' * 𝚺⁻¹ * 𝐗) * 𝐗' * 𝚺⁻¹
    return quadratic(𝐛' * 𝐀 * 𝐛, 2 * (𝐚' * 𝐀 * 𝐛), 𝐚' * 𝐀 * 𝐚 + 2 * size(𝐗, 2))
end
function (::BIC)(𝐗, 𝚺⁻¹, test_stat::TestStatistic)
    𝐚, 𝐛 = test_stat.𝐚, test_stat.𝐛
    𝐀 = 𝚺⁻¹ - 𝚺⁻¹ * 𝐗 * pinv(𝐗' * 𝚺⁻¹ * 𝐗) * 𝐗' * 𝚺⁻¹
    return quadratic(𝐛' * 𝐀 * 𝐛, 2 * (𝐚' * 𝐀 * 𝐛), 𝐚' * 𝐀 * 𝐚 + log(size(𝐗, 1)) * size(𝐗, 2))
end

# compute the region correspond to the selected model at `𝐲_z`
function region(𝐗, 𝐲_z, 𝚺,
    test_stat::TestStatistic,
    direction::Bidirection,
    criterion::AbstractCriterion;
    max_step=typemax(Int))
    alg_result = stepwise(𝐗, 𝐲_z, 𝚺, direction, criterion; max_step=max_step)
    history = alg_result.history

    d = size(𝐗, 2)
    model = Model(direction, d)
    𝚺⁻¹ = I / 𝚺
    Ω_z = [(-Inf)..Inf]

    for i in eachindex(history)
        @unpack M, Mᶜ = model
        # i-th selected feature
        idx = history[i]
        𝐗ᵢ = idx > 0 ? 𝐗[:, union(M, idx)] : 𝐗[:, setdiff(M, -idx)]
        ℓᵢ = criterion(𝐗ᵢ, 𝚺⁻¹, test_stat)

        # ℓᵢ is smaller than the value of previous model M
        i != 1 && (Ω_z = Ω_z ∩ (ℓᵢ ≤ criterion(𝐗[:, M], 𝚺⁻¹, test_stat)))
            # delete
        for j ∈ M
            Ω_z = Ω_z ∩ (ℓᵢ ≤ criterion(𝐗[:, setdiff(M, j)], 𝚺⁻¹, test_stat))
        end
        # add
        for j′ ∈ Mᶜ
            Ω_z = Ω_z ∩ (ℓᵢ ≤ criterion(𝐗[:, union(M, j′)], 𝚺⁻¹, test_stat))
        end

        update!(model, idx)
    end; @assert model == alg_result.model

    # stopping rule
    length(history) == max_step && return model, Ω_z
    @unpack M, Mᶜ = model
    minℓ = criterion(𝐗[:, M], 𝚺⁻¹, test_stat)
    ## delete
    for j ∈ M
        Ω_z = Ω_z ∩ (minℓ ≤ criterion(𝐗[:, setdiff(M, j)], 𝚺⁻¹, test_stat))
    end
    ## add
    for j′ ∈ Mᶜ
        Ω_z = Ω_z ∩ (minℓ ≤ criterion(𝐗[:, union(M, j′)], 𝚺⁻¹, test_stat))
    end

    return model, Ω_z
end
function region(𝐗, 𝐲_z, 𝚺,
    test_stat::TestStatistic,
    direction::Forward,
    criterion::AbstractCriterion;
    max_step=typemax(Int))
    alg_result = stepwise(𝐗, 𝐲_z, 𝚺, direction, criterion; max_step=max_step)
    history = alg_result.history

    d = size(𝐗, 2)
    model = Model(direction, d)
    𝚺⁻¹ = I / 𝚺
    Ω_z = [(-Inf)..Inf]

    for i in eachindex(history)
        @unpack M, Mᶜ = model
        # i-th selected feature
        idx = history[i]
        𝐗ᵢ = 𝐗[:, union(M, idx)]
        ℓᵢ = criterion(𝐗ᵢ, 𝚺⁻¹, test_stat)

        # ℓᵢ is smaller than the value of previous model M
        i != 1 && (Ω_z = Ω_z ∩ (ℓᵢ ≤ criterion(𝐗[:, M], 𝚺⁻¹, test_stat)))
        # add
        for j′ ∈ Mᶜ
            Ω_z = Ω_z ∩ (ℓᵢ ≤ criterion(𝐗[:, union(M, j′)], 𝚺⁻¹, test_stat))
        end

        update!(model, idx)
    end; @assert model == alg_result.model

    # stopping rule
    length(history) == max_step && return model, Ω_z
    @unpack M, Mᶜ = model
    minℓ = criterion(𝐗[:, M], 𝚺⁻¹, test_stat)
    ## add
    for j′ ∈ Mᶜ
        Ω_z = Ω_z ∩ (minℓ ≤ criterion(𝐗[:, union(M, j′)], 𝚺⁻¹, test_stat))
    end

    return model, Ω_z
end
function region(𝐗, 𝐲_z, 𝚺,
    test_stat::TestStatistic,
    direction::Backward,
    criterion::AbstractCriterion;
    max_step=typemax(Int))
    alg_result = stepwise(𝐗, 𝐲_z, 𝚺, direction, criterion; max_step=max_step)
    history = alg_result.history

    d = size(𝐗, 2)
    model = Model(direction, d)
    𝚺⁻¹ = I / 𝚺
    Ω_z = [(-Inf)..Inf]

    for i in eachindex(history)
        M = model.M
        # i-th selected feature
        idx = history[i]
        𝐗ᵢ = 𝐗[:, setdiff(M, -idx)]
        ℓᵢ = criterion(𝐗ᵢ, 𝚺⁻¹, test_stat)

        # ℓᵢ is smaller than the value of previous model M
        i != 1 && (Ω_z = Ω_z ∩ (ℓᵢ ≤ criterion(𝐗[:, M], 𝚺⁻¹, test_stat)))
            # delete
        for j ∈ M
            Ω_z = Ω_z ∩ (ℓᵢ ≤ criterion(𝐗[:, setdiff(M, j)], 𝚺⁻¹, test_stat))
        end

        update!(model, idx)
    end; @assert model == alg_result.model

    # stopping rule
    length(history) == max_step && return model, Ω_z
    M = model.M
    minℓ = criterion(𝐗[:, M], 𝚺⁻¹, test_stat)
    ## delete
    for j ∈ M
        Ω_z = Ω_z ∩ (minℓ ≤ criterion(𝐗[:, setdiff(M, j)], 𝚺⁻¹, test_stat))
    end

    return model, Ω_z
end
function region(𝐗, 𝐲_z, 𝚺, test_stat::TestStatistic, alg::StepwiseFeatureSelection)
    @unpack direction, criterion, max_step = alg
    region(𝐗, 𝐲_z, 𝚺, test_stat, direction, criterion; max_step=max_step)
end # wrapper of `region`

function compute_solution_path(𝐗, 𝚺, test_stat::TestStatistic, alg::AbstractAlgorithm;
    z_min=typemin(Int), z_max=typemax(Int))
    𝐚, 𝐛 = test_stat.𝐚, test_stat.𝐛
    z = float(z_min)
    𝐳 = [z]
    𝐌 = Model[]
    while z < z_max
        𝐲_z = 𝐚 + 𝐛 * z
        model, Ω_z = region(𝐗, 𝐲_z, 𝚺, test_stat, alg)
        for I_z ∈ Ω_z
            z ∈ I_z && (z = I_z.right + ε; break)
        end
        z < z_max ? push!(𝐳, z - ε) : push!(𝐳, z_max)
        push!(𝐌, model)
    end
    return 𝐳, 𝐌
end

function truncated_interval(model_obs, params...; kwargs...)
    𝐳, 𝐌 = compute_solution_path(params...; kwargs...)
        Z = ClosedInterval{eltype(𝐳)}[]
    for i in eachindex(𝐌)
        model_obs == 𝐌[i] && push!(Z, (𝐳[i]..(𝐳[i + 1])))
    end
    return Z
end

function selective_p(test_prob::SelectiveHypothesisTest)
    @unpack alg, alg_result, test_stats, 𝐗, 𝐲_obs, 𝚺 = test_prob
    model_obs = alg_result.model
    function p_value(test_stat::TestStatistic)
        z_obs, σ²_z = test_stat.z_obs, test_stat.σ²_z

        # compute Z = {z | selection_event}
        Z = truncated_interval(model_obs, 𝐗, 𝚺, test_stat, alg;
            z_min=-abs(z_obs) - 10 * √(σ²_z), z_max=abs(z_obs) + 10 * √(σ²_z))

        𝜋 = cdf(TruncatedDistribution(Normal(0, √(σ²_z)), Z), z_obs)
        return 2 * min(𝜋, 1 - 𝜋)
    end
    p_values = p_value.(test_stats)
    return p_values
end

function selective_CI(test_prob::SelectiveHypothesisTest; α=0.05)
    @unpack alg, alg_result, test_stats, 𝐗, 𝐲_obs, 𝚺 = test_prob
    model_obs = alg_result.model
    function CI(test_stat::TestStatistic)
        z_obs, σ²_z = test_stat.z_obs, test_stat.σ²_z

        # compute Z = {z | selection_event}
        Z = truncated_interval(model_obs, 𝐗, 𝚺, test_stat, alg;
            z_min=-abs(z_obs) - 10 * √(σ²_z), z_max=abs(z_obs) + 10 * √(σ²_z))

        return confidence_interval(z_obs, √(σ²_z), Z; α=α)
    end
    CIs = CI.(test_stats)
    return CIs
end

end # module


text:
using Test
using StaticArrays: SVector
using QuantumLattices.Essentials.QuantumSystems
using QuantumLattices.Essentials.QuantumSystems: heisenbergxyz, heisenbergpmz, gamma, dm
using QuantumLattices.Essentials.QuantumOperators: ID
using QuantumLattices.Essentials.Spatials: AbstractPID, PID, CPID, Point, Bond, rcoord, azimuthd
using QuantumLattices.Essentials.DegreesOfFreedom: Index, CompositeOID, OID, Hilbert, Operator, Operators, statistics, script, latexname, ishermitian
using QuantumLattices.Essentials.DegreesOfFreedom: IIDSpace, Coupling, Couplings, abbr, @subscript_str, @couplings, wildcard
using QuantumLattices.Interfaces: ⊗, ⋅, expand, permute, rank
using QuantumLattices.Prerequisites.Combinatorics: Permutations
using QuantumLattices.Prerequisites.VectorSpaces: shape, ndimshape

@testset "FID" begin
    fid = FID{:f}(orbital=1, spin=1)
    @test string(fid) == "FID{:f}(1, 1, 1)"
    @test statistics(fid) == statistics(typeof(fid)) == :f
    @test fid' == replace(fid, nambu=2)
    @test fid'' == replace(fid, nambu=1)

    fid = FID{:b}(1, 1, 0)
    @test string(fid) == "FID{:b}(1, 1, 0)"
    @test statistics(fid) == statistics(typeof(fid)) == :b
    @test fid' == fid
end

@testset "Fock" begin
    fock = Fock{:b}(norbital=1, nspin=2, nnambu=2)
    @test shape(fock) == (1:1, 1:2, 1:2)
    @test ndimshape(fock) == ndimshape(typeof(fock)) == 3
    @test CartesianIndex(FID{:b}(1, 1, 1), fock) == CartesianIndex(1, 1, 1)
    @test FID(CartesianIndex(1, 1, 1), fock) == FID{:b}(1, 1, 1)
    @test collect(fock) == [FID{:b}(1, 1, 1), FID{:b}(1, 2, 1), FID{:b}(1, 1, 2), FID{:b}(1, 2, 2)]
    @test statistics(fock) == statistics(typeof(fock)) == :b
    @test string(fock) == "Fock{:b}(norbital=1, nspin=2, nnambu=2)"

    @test summary(Fock{:b}(nspin=0, nnambu=1)) == "0-element Fock{:b}"
    @test summary(Fock{:b}(nspin=1, nnambu=1)) == "1-element Fock{:b}"
    @test summary(Fock{:f}(nspin=2, nnambu=1)) == "2-element Fock{:f}"

    @test match(FID{wildcard}, Fock{:f}) == match(FID{wildcard}, Fock{:b}) == true
    @test match(FID{:f}, Fock{:f}) == match(FID{:b}, Fock{:b}) == true
    @test match(FID{:b}, Fock{:f}) == match(FID{:f}, Fock{:b}) == false
end

@testset "latex" begin
    @test script(Val(:site), Index(PID(1), FID{:f}(2, 1, 1))) == 1
    @test script(Val(:orbital), Index(PID(1), FID{:f}(2, 1, 1))) == 2
    @test script(Val(:spint), Index(PID(1), FID{:f}(2, 3, 1))) == 3
    @test script(Val(:spsym), Index(PID(1), FID{:f}(2, 2, 1))) == "↑"
    @test script(Val(:spsym), Index(PID(1), FID{:f}(2, 1, 1))) == "↓"
    @test script(Val(:nambu), Index(PID(1), FID{:f}(2, 3, 1))) == ""
    @test script(Val(:nambu), Index(PID(1), FID{:f}(2, 3, 2))) == "\\dagger"

    @test latexname(Index{<:AbstractPID, <:FID{:f}}) == Symbol("Index{AbstractPID, FID{:f}}")
    @test latexname(CompositeOID{Index{<:AbstractPID, <:FID{:f}}}) == Symbol("CompositeOID{Index{AbstractPID, FID{:f}}}")
    @test latexname(Index{<:AbstractPID, <:FID{:b}}) == Symbol("Index{AbstractPID, FID{:b}}")
    @test latexname(CompositeOID{Index{<:AbstractPID, <:FID{:b}}}) == Symbol("CompositeOID{Index{AbstractPID, FID{:b}}}")
end

@testset "angle" begin
    @test angle(OID(Index(PID(1), FID{:f}(1, 1, 1)), [0.0, 0.0], [1.0, 2.0]), [[1.0, 0.0], [0.0, 1.0]], [0.1, 0.0]) ≈ 2pi*0.1
    @test angle(OID(Index(CPID(1, 1), FID{:f}(1, 1, 2)), [0.0, 0.0], [1.0, 2.0]), [[1.0, 0.0], [0.0, 1.0]], [0.0, 0.2]) ≈ -2pi*0.4
end

@testset "FockOperator" begin
    id₁ = OID(Index(PID(2), FID{:f}(1, 1, 2)), SVector(0.5, 0.0), SVector(0.0, 0.0))
    id₂ = OID(Index(PID(2), FID{:f}(1, 1, 1)), SVector(0.5, 0.0), SVector(0.0, 0.0))
    id₃ = OID(Index(PID(1), FID{:f}(1, 2, 2)), SVector(0.0, 0.0), SVector(0.0, 0.0))
    id₄ = OID(Index(PID(1), FID{:f}(1, 2, 1)), SVector(0.0, 0.0), SVector(0.0, 0.0))

    opt = Operator(1.0, ID(id₁, id₂))
    @test opt|>isnormalordered

    opt = Operator(1.0, ID(id₁, id₂, id₃, id₄))
    @test opt|>isnormalordered == false
    @test repr(opt) == "c^{\\dagger}_{2, 1, ↓}c^{}_{2, 1, ↓}c^{\\dagger}_{1, 1, ↑}c^{}_{1, 1, ↑}"

    op₁ = Operator(1.5, ID(id₁, id₂))
    op₂ = Operator(2.0, ID(id₂, id₁))
    @test op₁*op₂ == nothing

    op₁ = Operator(1.5, ID(id₁, id₂))
    op₂ = Operator(2.0, ID(id₁, id₂))
    @test op₁*op₂ == Operator(3.0, ID(id₁, id₂, id₁, id₂))

    @test permute(id₁, id₂) == (Operator(1), Operator(-1, ID(id₂, id₁)))
    @test permute(id₂, id₁) == (Operator(1), Operator(-1, ID(id₁, id₂)))
    @test permute(id₁, id₄) == (Operator(-1, ID(id₄, id₁)),)
    @test permute(id₄, id₁) == (Operator(-1, ID(id₁, id₄)),)

    id₁ = OID(Index(PID(2), FID{:b}(1, 1, 2)), SVector(0.5, 0.0), SVector(0.0, 0.0))
    id₂ = OID(Index(PID(2), FID{:b}(1, 1, 1)), SVector(0.5, 0.0), SVector(0.0, 0.0))
    id₃ = OID(Index(PID(1), FID{:b}(1, 2, 2)), SVector(0.0, 0.0), SVector(0.0, 0.0))
    id₄ = OID(Index(PID(1), FID{:b}(1, 2, 1)), SVector(0.0, 0.0), SVector(0.0, 0.0))

    opt = Operator(1.0, ID(id₁, id₂))
    @test repr(opt) == "b^{\\dagger}_{2, 1, ↓}b^{}_{2, 1, ↓}"

    @test permute(id₁, id₂) == (Operator(+1), Operator(1, ID(id₂, id₁)))
    @test permute(id₂, id₁) == (Operator(-1), Operator(1, ID(id₁, id₂)))
    @test permute(id₁, id₄) == (Operator(1, ID(id₄, id₁)),)
    @test permute(id₄, id₁) == (Operator(1, ID(id₁, id₄)),)
end

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

    fc₁ = FockCoupling{2}(1.5, spins=subscript"[x 1]")
    fc₂ = FockCoupling{2}(2.0, orbitals=subscript"[x y](x < y)")
    @test repr(fc₁) == "1.5 sp[x 1]"
    @test repr(fc₂) == "2.0 ob[x y](x < y)"
    fc = fc₁ * fc₂
    @test repr(fc) == "3.0 ob[* *]×[x y](x < y) ⊗ sp[x 1]×[* *]"

    fc₁ = FockCoupling{2}(1.5, spins=subscript"[x 1]")
    fc₂ = FockCoupling{2}(2.0, orbitals=subscript"[x y](x < y)")
    fc = fc₁ ⊗ fc₂
    @test repr(fc) == "3.0 ob[x y](x < y) ⊗ sp[x 1]"

    fc₁ = FockCoupling{2}(1.5, spins=(2, 1))
    fc₂ = FockCoupling{2}(2.0, spins=(1, 2))
    fc = fc₁ ⋅ fc₂
    @test repr(fc) == "3.0 sp[2 2]"

    fc = FockCoupling{2}(2.0, orbitals=(1, 2), nambus=(2, 1))
    bond = Bond(1, Point(CPID(1, 2), SVector(0.5), SVector(0.0)), Point(CPID(1, 1), SVector(0.0), SVector(0.0)))
    hilbert = Hilbert(pid=>Fock{:f}(norbital=2, nspin=2, nnambu=2) for pid in [bond.epoint.pid, bond.spoint.pid])
    ex = expand(fc, bond, hilbert, Val(:Hopping))
    @test collect(ex) == [
        (2.0, ID(OID(Index(CPID(1, 1), FID{:f}(1, 1, 2)), SVector(0.0), SVector(0.0)), OID(Index(CPID(1, 2), FID{:f}(2, 1, 1)), SVector(0.5), SVector(0.0)))),
        (2.0, ID(OID(Index(CPID(1, 1), FID{:f}(1, 2, 2)), SVector(0.0), SVector(0.0)), OID(Index(CPID(1, 2), FID{:f}(2, 2, 1)), SVector(0.5), SVector(0.0))))
    ]

    fc = FockCoupling{4}(2.0, spins=(2, 2, 1, 1), nambus=(2, 1, 2, 1))
    point = Point(PID(1), SVector(0.0), SVector(0.0))
    hilbert = Hilbert(point.pid=>Fock{:b}(norbital=2, nspin=2, nnambu=2))
    ex = expand(fc, point, hilbert, Val(:info))
    @test collect(ex) == [
        (2.0, ID(OID(Index(PID(1), FID{:b}(1, 2, 2)), SVector(0.0), SVector(0.0)),
                 OID(Index(PID(1), FID{:b}(1, 2, 1)), SVector(0.0), SVector(0.0)),
                 OID(Index(PID(1), FID{:b}(1, 1, 2)), SVector(0.0), SVector(0.0)),
                 OID(Index(PID(1), FID{:b}(1, 1, 1)), SVector(0.0), SVector(0.0))
                 )),
        (2.0, ID(OID(Index(PID(1), FID{:b}(2, 2, 2)), SVector(0.0), SVector(0.0)),
                 OID(Index(PID(1), FID{:b}(2, 2, 1)), SVector(0.0), SVector(0.0)),
                 OID(Index(PID(1), FID{:b}(2, 1, 2)), SVector(0.0), SVector(0.0)),
                 OID(Index(PID(1), FID{:b}(2, 1, 1)), SVector(0.0), SVector(0.0))
                 ))
    ]

    fc = FockCoupling{4}(2.0, orbitals=subscript"[α α β β](α < β)", spins=(2, 1, 1, 2), nambus=(2, 2, 1, 1))
    point = Point(PID(1), SVector(0.5), SVector(0.0))
    hilbert = Hilbert(point.pid=>Fock{:f}(norbital=3, nspin=2, nnambu=2))
    ex = expand(fc, point, hilbert, Val(:info))
    @test collect(ex) == [
        (2.0, ID(OID(Index(PID(1), FID{:f}(1, 2, 2)), SVector(0.5), SVector(0.0)),
                 OID(Index(PID(1), FID{:f}(1, 1, 2)), SVector(0.5), SVector(0.0)),
                 OID(Index(PID(1), FID{:f}(2, 1, 1)), SVector(0.5), SVector(0.0)),
                 OID(Index(PID(1), FID{:f}(2, 2, 1)), SVector(0.5), SVector(0.0))
                 )),
        (2.0, ID(OID(Index(PID(1), FID{:f}(1, 2, 2)), SVector(0.5), SVector(0.0)),
                 OID(Index(PID(1), FID{:f}(1, 1, 2)), SVector(0.5), SVector(0.0)),
                 OID(Index(PID(1), FID{:f}(3, 1, 1)), SVector(0.5), SVector(0.0)),
                 OID(Index(PID(1), FID{:f}(3, 2, 1)), SVector(0.5), SVector(0.0))
                 )),
        (2.0, ID(OID(Index(PID(1), FID{:f}(2, 2, 2)), SVector(0.5), SVector(0.0)),
                 OID(Index(PID(1), FID{:f}(2, 1, 2)), SVector(0.5), SVector(0.0)),
                 OID(Index(PID(1), FID{:f}(3, 1, 1)), SVector(0.5), SVector(0.0)),
                 OID(Index(PID(1), FID{:f}(3, 2, 1)), SVector(0.5), SVector(0.0))
                 ))
    ]

    fc₁ = FockCoupling{2}(+1.0, spins=(2, 2), nambus=(2, 1))
    fc₂ = FockCoupling{2}(-1.0, spins=(1, 1), nambus=(2, 1))
    point = Point(CPID(1, 1), SVector(0.0), SVector(0.0))
    hilbert = Hilbert(point.pid=>Fock{:f}(norbital=2, nspin=2, nnambu=2))
    ex = expand(fc₁*fc₂, point, hilbert, Val(:info))
    @test collect(ex) == [
        (-1.0, ID(OID(Index(CPID(1, 1), FID{:f}(1, 2, 2)), SVector(0.0), SVector(0.0)),
                  OID(Index(CPID(1, 1), FID{:f}(1, 2, 1)), SVector(0.0), SVector(0.0)),
                  OID(Index(CPID(1, 1), FID{:f}(1, 1, 2)), SVector(0.0), SVector(0.0)),
                  OID(Index(CPID(1, 1), FID{:f}(1, 1, 1)), SVector(0.0), SVector(0.0))
                  )),
        (-1.0, ID(OID(Index(CPID(1, 1), FID{:f}(2, 2, 2)), SVector(0.0), SVector(0.0)),
                  OID(Index(CPID(1, 1), FID{:f}(2, 2, 1)), SVector(0.0), SVector(0.0)),
                  OID(Index(CPID(1, 1), FID{:f}(1, 1, 2)), SVector(0.0), SVector(0.0)),
                  OID(Index(CPID(1, 1), FID{:f}(1, 1, 1)), SVector(0.0), SVector(0.0))
                  )),
        (-1.0, ID(OID(Index(CPID(1, 1), FID{:f}(1, 2, 2)), SVector(0.0), SVector(0.0)),
                  OID(Index(CPID(1, 1), FID{:f}(1, 2, 1)), SVector(0.0), SVector(0.0)),
                  OID(Index(CPID(1, 1), FID{:f}(2, 1, 2)), SVector(0.0), SVector(0.0)),
                  OID(Index(CPID(1, 1), FID{:f}(2, 1, 1)), SVector(0.0), SVector(0.0))
                  )),
        (-1.0, ID(OID(Index(CPID(1, 1), FID{:f}(2, 2, 2)), SVector(0.0), SVector(0.0)),
                  OID(Index(CPID(1, 1), FID{:f}(2, 2, 1)), SVector(0.0), SVector(0.0)),
                  OID(Index(CPID(1, 1), FID{:f}(2, 1, 2)), SVector(0.0), SVector(0.0)),
                  OID(Index(CPID(1, 1), FID{:f}(2, 1, 1)), SVector(0.0), SVector(0.0))
                  ))
    ]
end

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

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

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

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

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

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

@testset "fockcoupling" begin
    fc = fc"1.0 ob[α α β β](α < β) ⊗ sp[σ γ σ γ](σ ≠ γ) ⊗ ph[2 1 2 1]"
    @test repr(fc) == "1.0 ob[α α β β](α < β) ⊗ sp[σ γ σ γ](σ ≠ γ) ⊗ ph[2 1 2 1]"

    fc = fc"1.0 ob[α α β β] ⊗ sp[σ γ σ γ] ⊗ ph[2 1 2 1]"
    @test repr(fc) == "1.0 ob[α α β β] ⊗ sp[σ γ σ γ] ⊗ ph[2 1 2 1]"

    fc = fc"1.0 ob[α α β β](α < β) ⊗ sp[σ γ σ γ] ⊗ ph[2 1 2 1]"
    @test repr(fc) == "1.0 ob[α α β β](α < β) ⊗ sp[σ γ σ γ] ⊗ ph[2 1 2 1]"

    fc = fc"1.0 ob[α α β β](α < β) ⊗ ph[2 1 2 1]"
    @test repr(fc) == "1.0 ob[α α β β](α < β) ⊗ ph[2 1 2 1]"

    fc = fc"1.0 ob[α α β β](α < β) ⊗ sp[2 1 2 1]"
    @test repr(fc) == "1.0 ob[α α β β](α < β) ⊗ sp[2 1 2 1]"

    fc = fc"1.0 ob[1 1 1 1] ⊗ ph[2 1 2 1]"
    @test repr(fc) == "1.0 ob[1 1 1 1] ⊗ ph[2 1 2 1]"

    fc = fc"1.0 ph[2 1 2 1]"
    @test repr(fc) == "1.0 ph[2 1 2 1]"

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

@testset "Onsite" begin
    @test abbr(Onsite) == :st
    @test isnothing(ishermitian(Onsite))

    point = Point(PID(1), (0.5, 0.5), (0.0, 0.0))
    hilbert = Hilbert(point.pid=>Fock{:f}(norbital=2, nspin=2, nnambu=2))

    term = Onsite(:mu, 1.5, couplings=σˣ"sp"⊗σᶻ"ob", modulate=true)
    operators = Operators(
        Operator(+1.5, ID(OID(Index(PID(1), FID{:f}(2, 1, 2)), [0.5, 0.5], [0.0, 0.0]), OID(Index(PID(1), FID{:f}(2, 2, 1)), [0.5, 0.5], [0.0, 0.0]))),
        Operator(-1.5, ID(OID(Index(PID(1), FID{:f}(1, 1, 2)), [0.5, 0.5], [0.0, 0.0]), OID(Index(PID(1), FID{:f}(1, 2, 1)), [0.5, 0.5], [0.0, 0.0])))
    )
    @test expand(term, point, hilbert, half=true) == operators
    @test expand(term, point, hilbert, half=false) == operators+operators'

    term = Onsite(:mu, 1.5, couplings=σᶻ"sp"⊗σᶻ"ob", modulate=true)
    operators = Operators(
        Operator(-0.75, ID(OID(Index(PID(1), FID{:f}(2, 1, 2)), [0.5, 0.5], [0.0, 0.0]), OID(Index(PID(1), FID{:f}(2, 1, 1)), [0.5, 0.5], [0.0, 0.0]))),
        Operator(-0.75, ID(OID(Index(PID(1), FID{:f}(1, 2, 2)), [0.5, 0.5], [0.0, 0.0]), OID(Index(PID(1), FID{:f}(1, 2, 1)), [0.5, 0.5], [0.0, 0.0]))),
        Operator(+0.75, ID(OID(Index(PID(1), FID{:f}(2, 2, 2)), [0.5, 0.5], [0.0, 0.0]), OID(Index(PID(1), FID{:f}(2, 2, 1)), [0.5, 0.5], [0.0, 0.0]))),
        Operator(+0.75, ID(OID(Index(PID(1), FID{:f}(1, 1, 2)), [0.5, 0.5], [0.0, 0.0]), OID(Index(PID(1), FID{:f}(1, 1, 1)), [0.5, 0.5], [0.0, 0.0])))
    )
    @test expand(term, point, hilbert, half=true) == operators
    @test expand(term, point, hilbert, half=false) == operators+operators'
end

@testset "Hopping" begin
    bond = Bond(1, Point(CPID('a', 1), (0.5, 0.5), (0.0, 0.0)), Point(CPID('b', 2), (0.0, 0.0), (0.0, 0.0)))
    hilbert = Hilbert(pid=>Fock{:f}(norbital=2, nspin=2, nnambu=2) for pid in [bond.spoint.pid, bond.epoint.pid])
    term = Hopping(:t, 1.5, 1)
    operators = Operators(
        Operator(1.5, ID(OID(Index(CPID('b', 2), FID{:f}(2, 2, 2)), [0.0, 0.0], [0.0, 0.0]), OID(Index(CPID('a', 1), FID{:f}(2, 2, 1)), [0.5, 0.5], [0.0, 0.0]))),
        Operator(1.5, ID(OID(Index(CPID('b', 2), FID{:f}(2, 1, 2)), [0.0, 0.0], [0.0, 0.0]), OID(Index(CPID('a', 1), FID{:f}(2, 1, 1)), [0.5, 0.5], [0.0, 0.0]))),
        Operator(1.5, ID(OID(Index(CPID('b', 2), FID{:f}(1, 1, 2)), [0.0, 0.0], [0.0, 0.0]), OID(Index(CPID('a', 1), FID{:f}(1, 1, 1)), [0.5, 0.5], [0.0, 0.0]))),
        Operator(1.5, ID(OID(Index(CPID('b', 2), FID{:f}(1, 2, 2)), [0.0, 0.0], [0.0, 0.0]), OID(Index(CPID('a', 1), FID{:f}(1, 2, 1)), [0.5, 0.5], [0.0, 0.0])))
    )
    @test term|>abbr == term|>typeof|>abbr == :hp
    @test term|>ishermitian == term|>typeof|>ishermitian == false
    @test expand(term, bond, hilbert, half=true) == operators
    @test expand(term, bond, hilbert, half=false) == operators+operators'
end

@testset "Pairing" begin
    bond = Bond(1, Point(PID(1), (0.5, 0.5), (0.0, 0.0)), Point(PID(2), (0.0, 0.0), (0.0, 0.0)))
    hilbert = Hilbert(pid=>Fock{:f}(norbital=1, nspin=2, nnambu=2) for pid in [bond.spoint.pid, bond.epoint.pid])
    term = Pairing(:Δ, 1.5, 1, couplings=@couplings(FockCoupling{2}(spins=(2, 2))), amplitude=bond->(bond|>rcoord|>azimuthd ≈ 45 ? 1 : -1))
    operators = Operators(
        Operator(-1.5, ID(OID(Index(PID(2), FID{:f}(1, 2, 1)), [0.0, 0.0], [0.0, 0.0]), OID(Index(PID(1), FID{:f}(1, 2, 1)), [0.5, 0.5], [0.0, 0.0]))),
        Operator(+1.5, ID(OID(Index(PID(1), FID{:f}(1, 2, 1)), [0.5, 0.5], [0.0, 0.0]), OID(Index(PID(2), FID{:f}(1, 2, 1)), [0.0, 0.0], [0.0, 0.0])))
    )
    @test term|>abbr == term|>typeof|>abbr == :pr
    @test term|>ishermitian == term|>typeof|>ishermitian == false
    @test expand(term, bond, hilbert, half=true) == operators
    @test expand(term, bond, hilbert, half=false) == operators+operators'

    point = Point(CPID('a', 1), (0.5, 0.5), (0.0, 0.0))
    hilbert = Hilbert(point.pid=>Fock{:f}(norbital=1, nspin=2, nnambu=2))
    term = Pairing(:Δ, 1.5, 0, couplings=FockCoupling{2}(spins=(2, 1))-FockCoupling{2}(spins=(1, 2)))
    operators = Operators(
        Operator(+1.5, ID(OID(Index(CPID('a', 1), FID{:f}(1, 2, 1)), [0.5, 0.5], [0.0, 0.0]), OID(Index(CPID('a', 1), FID{:f}(1, 1, 1)), [0.5, 0.5], [0.0, 0.0]))),
        Operator(-1.5, ID(OID(Index(CPID('a', 1), FID{:f}(1, 1, 1)), [0.5, 0.5], [0.0, 0.0]), OID(Index(CPID('a', 1), FID{:f}(1, 2, 1)), [0.5, 0.5], [0.0, 0.0])))
    )
    @test term|>abbr == term|>typeof|>abbr == :pr
    @test expand(term, point, hilbert, half=true) == operators
    @test expand(term, point, hilbert, half=false) == operators+operators'
end

@testset "Hubbard" begin
    point = Point(PID(1), (0.5, 0.5), (0.0, 0.0))
    hilbert = Hilbert(point.pid=>Fock{:f}(norbital=2, nspin=2, nnambu=2))
    term = Hubbard(:H, 2.5)
    operators = Operators(
        Operator(1.25, ID(
            OID(Index(PID(1), FID{:f}(1, 2, 2)), [0.5, 0.5], [0.0, 0.0]),
            OID(Index(PID(1), FID{:f}(1, 2, 1)), [0.5, 0.5], [0.0, 0.0]),
            OID(Index(PID(1), FID{:f}(1, 1, 2)), [0.5, 0.5], [0.0, 0.0]),
            OID(Index(PID(1), FID{:f}(1, 1, 1)), [0.5, 0.5], [0.0, 0.0])
            )),
        Operator(1.25, ID(
            OID(Index(PID(1), FID{:f}(2, 2, 2)), [0.5, 0.5], [0.0, 0.0]),
            OID(Index(PID(1), FID{:f}(2, 2, 1)), [0.5, 0.5], [0.0, 0.0]),
            OID(Index(PID(1), FID{:f}(2, 1, 2)), [0.5, 0.5], [0.0, 0.0]),
            OID(Index(PID(1), FID{:f}(2, 1, 1)), [0.5, 0.5], [0.0, 0.0])
            ))
    )
    @test term|>abbr == term|>typeof|>abbr == :hb
    @test term|>ishermitian == term|>typeof|>ishermitian == true
    @test expand(term, point, hilbert, half=true) == operators
    @test expand(term, point, hilbert, half=false) == operators*2
end

@testset "InterOrbitalInterSpin" begin
    point = Point(PID(1), (0.5, 0.5), (0.0, 0.0))
    hilbert = Hilbert(point.pid=>Fock{:f}(norbital=2, nspin=2, nnambu=2))
    term = InterOrbitalInterSpin(:H, 2.5)
    operators = Operators(
        Operator(1.25, ID(
            OID(Index(PID(1), FID{:f}(1, 2, 2)), [0.5, 0.5], [0.0, 0.0]),
            OID(Index(PID(1), FID{:f}(1, 2, 1)), [0.5, 0.5], [0.0, 0.0]),
            OID(Index(PID(1), FID{:f}(2, 1, 2)), [0.5, 0.5], [0.0, 0.0]),
            OID(Index(PID(1), FID{:f}(2, 1, 1)), [0.5, 0.5], [0.0, 0.0])
            )),
        Operator(1.25, ID(
            OID(Index(PID(1), FID{:f}(1, 1, 2)), [0.5, 0.5], [0.0, 0.0]),
            OID(Index(PID(1), FID{:f}(1, 1, 1)), [0.5, 0.5], [0.0, 0.0]),
            OID(Index(PID(1), FID{:f}(2, 2, 2)), [0.5, 0.5], [0.0, 0.0]),
            OID(Index(PID(1), FID{:f}(2, 2, 1)), [0.5, 0.5], [0.0, 0.0])
            ))
    )
    @test term|>abbr == term|>typeof|>abbr == :nons
    @test term|>ishermitian == term|>typeof|>ishermitian == true
    @test expand(term, point, hilbert, half=true) == operators
    @test expand(term, point, hilbert, half=false) == operators*2
end

@testset "InterOrbitalIntraSpin" begin
    point = Point(PID(1), (0.5, 0.5), (0.0, 0.0))
    hilbert = Hilbert(point.pid=>Fock{:f}(norbital=2, nspin=2, nnambu=2))
    term = InterOrbitalIntraSpin(:H, 2.5)
    operators = Operators(
        Operator(1.25, ID(
            OID(Index(PID(1), FID{:f}(1, 1, 2)), [0.5, 0.5], [0.0, 0.0]),
            OID(Index(PID(1), FID{:f}(1, 1, 1)), [0.5, 0.5], [0.0, 0.0]),
            OID(Index(PID(1), FID{:f}(2, 1, 2)), [0.5, 0.5], [0.0, 0.0]),
            OID(Index(PID(1), FID{:f}(2, 1, 1)), [0.5, 0.5], [0.0, 0.0])
            )),
        Operator(1.25, ID(
            OID(Index(PID(1), FID{:f}(1, 2, 2)), [0.5, 0.5], [0.0, 0.0]),
            OID(Index(PID(1), FID{:f}(1, 2, 1)), [0.5, 0.5], [0.0, 0.0]),
            OID(Index(PID(1), FID{:f}(2, 2, 2)), [0.5, 0.5], [0.0, 0.0]),
            OID(Index(PID(1), FID{:f}(2, 2, 1)), [0.5, 0.5], [0.0, 0.0])
            ))
    )
    @test term|>abbr == term|>typeof|>abbr == :noes
    @test term|>ishermitian == term|>typeof|>ishermitian == true
    @test expand(term, point, hilbert, half=true) == operators
    @test expand(term, point, hilbert, half=false) == operators*2
end

@testset "SpinFlip" begin
    point = Point(PID(1), (0.5, 0.5), (0.0, 0.0))
    hilbert = Hilbert(point.pid=>Fock{:f}(norbital=2, nspin=2, nnambu=2))
    term = SpinFlip(:H, 2.5)
    operators = Operators(
        Operator(2.5, ID(
            OID(Index(PID(1), FID{:f}(1, 2, 2)), [0.5, 0.5], [0.0, 0.0]),
            OID(Index(PID(1), FID{:f}(2, 1, 2)), [0.5, 0.5], [0.0, 0.0]),
            OID(Index(PID(1), FID{:f}(1, 1, 1)), [0.5, 0.5], [0.0, 0.0]),
            OID(Index(PID(1), FID{:f}(2, 2, 1)), [0.5, 0.5], [0.0, 0.0])
            ))
    )
    @test term|>abbr == term|>typeof|>abbr == :sf
    @test term|>ishermitian == term|>typeof|>ishermitian == false
    @test expand(term, point, hilbert, half=true) == operators
    @test expand(term, point, hilbert, half=false) == operators+operators'
end

@testset "PairHopping" begin
    point = Point(PID(1), (0.5, 0.5), (0.0, 0.0))
    hilbert = Hilbert(point.pid=>Fock{:f}(norbital=2, nspin=2, nnambu=2))
    term = PairHopping(:H, 2.5)
    operators = Operators(
        Operator(2.5, ID(
            OID(Index(PID(1), FID{:f}(1, 2, 2)), [0.5, 0.5], [0.0, 0.0]),
            OID(Index(PID(1), FID{:f}(1, 1, 2)), [0.5, 0.5], [0.0, 0.0]),
            OID(Index(PID(1), FID{:f}(2, 1, 1)), [0.5, 0.5], [0.0, 0.0]),
            OID(Index(PID(1), FID{:f}(2, 2, 1)), [0.5, 0.5], [0.0, 0.0])
            ))
    )
    @test term|>abbr == term|>typeof|>abbr == :ph
    @test term|>ishermitian == term|>typeof|>ishermitian == false
    @test expand(term, point, hilbert, half=true) == operators
    @test expand(term, point, hilbert, half=false) == operators+operators'
end

@testset "Coulomb" begin
    bond = Bond(1, Point(PID(1), (0.5, 0.5), (0.0, 0.0)), Point(PID(2), (0.0, 0.0), (0.0, 0.0)))
    hilbert = Hilbert(pid=>Fock{:f}(norbital=1, nspin=2, nnambu=2) for pid in [bond.spoint.pid, bond.epoint.pid])

    term = Coulomb(:V, 2.5, 1, couplings=σᶻ"sp"*σᶻ"sp")
    operators = Operators(
        Operator(-1.25, ID(
            OID(Index(PID(2), FID{:f}(1, 1, 2)), [0.0, 0.0], [0.0, 0.0]),
            OID(Index(PID(2), FID{:f}(1, 1, 1)), [0.0, 0.0], [0.0, 0.0]),
            OID(Index(PID(1), FID{:f}(1, 2, 2)), [0.5, 0.5], [0.0, 0.0]),
            OID(Index(PID(1), FID{:f}(1, 2, 1)), [0.5, 0.5], [0.0, 0.0])
            )),
        Operator(+1.25, ID(
            OID(Index(PID(2), FID{:f}(1, 1, 2)), [0.0, 0.0], [0.0, 0.0]),
            OID(Index(PID(2), FID{:f}(1, 1, 1)), [0.0, 0.0], [0.0, 0.0]),
            OID(Index(PID(1), FID{:f}(1, 1, 2)), [0.5, 0.5], [0.0, 0.0]),
            OID(Index(PID(1), FID{:f}(1, 1, 1)), [0.5, 0.5], [0.0, 0.0])
            )),
        Operator(-1.25, ID(
            OID(Index(PID(2), FID{:f}(1, 2, 2)), [0.0, 0.0], [0.0, 0.0]),
            OID(Index(PID(2), FID{:f}(1, 2, 1)), [0.0, 0.0], [0.0, 0.0]),
            OID(Index(PID(1), FID{:f}(1, 1, 2)), [0.5, 0.5], [0.0, 0.0]),
            OID(Index(PID(1), FID{:f}(1, 1, 1)), [0.5, 0.5], [0.0, 0.0])
            )),
        Operator(+1.25, ID(
            OID(Index(PID(2), FID{:f}(1, 2, 2)), [0.0, 0.0], [0.0, 0.0]),
            OID(Index(PID(2), FID{:f}(1, 2, 1)), [0.0, 0.0], [0.0, 0.0]),
            OID(Index(PID(1), FID{:f}(1, 2, 2)), [0.5, 0.5], [0.0, 0.0]),
            OID(Index(PID(1), FID{:f}(1, 2, 1)), [0.5, 0.5], [0.0, 0.0])
            ))
    )
    @test term|>abbr == term|>typeof|>abbr == :cl
    @test isnothing(term|>ishermitian) && isnothing(term|>typeof|>ishermitian)
    @test expand(term, bond, hilbert, half=true) == operators
    @test expand(term, bond, hilbert, half=false) == operators+operators'

    term = Coulomb(:V, 2.5, 1, couplings=σˣ"sp"*σᶻ"sp")
    operators = Operators(
        Operator(-2.5, ID(
            OID(Index(PID(2), FID{:f}(1, 2, 2)), [0.0, 0.0], [0.0, 0.0]),
            OID(Index(PID(2), FID{:f}(1, 1, 1)), [0.0, 0.0], [0.0, 0.0]),
            OID(Index(PID(1), FID{:f}(1, 1, 2)), [0.5, 0.5], [0.0, 0.0]),
            OID(Index(PID(1), FID{:f}(1, 1, 1)), [0.5, 0.5], [0.0, 0.0])
            )),
        Operator(+2.5, ID(
            OID(Index(PID(2), FID{:f}(1, 1, 2)), [0.0, 0.0], [0.0, 0.0]),
            OID(Index(PID(2), FID{:f}(1, 2, 1)), [0.0, 0.0], [0.0, 0.0]),
            OID(Index(PID(1), FID{:f}(1, 2, 2)), [0.5, 0.5], [0.0, 0.0]),
            OID(Index(PID(1), FID{:f}(1, 2, 1)), [0.5, 0.5], [0.0, 0.0])
            )),
        Operator(+2.5, ID(
            OID(Index(PID(2), FID{:f}(1, 2, 2)), [0.0, 0.0], [0.0, 0.0]),
            OID(Index(PID(2), FID{:f}(1, 1, 1)), [0.0, 0.0], [0.0, 0.0]),
            OID(Index(PID(1), FID{:f}(1, 2, 2)), [0.5, 0.5], [0.0, 0.0]),
            OID(Index(PID(1), FID{:f}(1, 2, 1)), [0.5, 0.5], [0.0, 0.0])
            )),
        Operator(-2.5, ID(
            OID(Index(PID(2), FID{:f}(1, 1, 2)), [0.0, 0.0], [0.0, 0.0]),
            OID(Index(PID(2), FID{:f}(1, 2, 1)), [0.0, 0.0], [0.0, 0.0]),
            OID(Index(PID(1), FID{:f}(1, 1, 2)), [0.5, 0.5], [0.0, 0.0]),
            OID(Index(PID(1), FID{:f}(1, 1, 1)), [0.5, 0.5], [0.0, 0.0])
            ))
    )
    @test term|>abbr == term|>typeof|>abbr == :cl
    @test expand(term, bond, hilbert, half=true) == operators
    @test expand(term, bond, hilbert, half=false) == operators+operators'
end

@testset "SID" begin
    @test SID{1//2}('z', orbital=1)' == SID{1//2}('z', orbital=1)
    @test SID{1//2}('x', orbital=1)' == SID{1//2}('x', orbital=1)
    @test SID{1//2}('y', orbital=1)' == SID{1//2}('y', orbital=1)
    @test SID{1//2}('+', orbital=1)' == SID{1//2}('-', orbital=1)
    @test SID{1//2}('-', orbital=1)' == SID{1//2}('+', orbital=1)

    sid = SID{3//2}('x')
    @test string(sid) == "SID{3//2}(1, 'x')"
    @test replace(sid, tag='z') == SID{3//2}('z')
    @test totalspin(sid) == totalspin(typeof(sid)) == 3//2
    @test statistics(sid) == statistics(typeof(sid)) == :b
end

@testset "Matrix" begin
    @test isapprox(Matrix(SID{1//2}(1, 'z')), [[-0.5, 0.0] [0.0, 0.5]])
    @test isapprox(Matrix(SID{1//2}(1, 'x')), [[0.0, 0.5] [0.5, 0.0]])
    @test isapprox(Matrix(SID{1//2}(1, 'y')), [[0.0, -0.5im] [0.5im, 0.0]])
    @test isapprox(Matrix(SID{1//2}(1, '+')), [[0.0, 1.0] [0.0, 0.0]])
    @test isapprox(Matrix(SID{1//2}(1, '-')), [[0.0, 0.0] [1.0, 0.0]])

    @test isapprox(Matrix(SID{1}(1, 'z')), [[-1.0, 0.0, 0.0] [0.0, 0.0, 0.0] [0.0, 0.0, 1.0]])
    @test isapprox(Matrix(SID{1}(1, 'x')), [[0.0, √2/2, 0.0] [√2/2, 0.0, √2/2] [0.0, √2/2, 0.0]])
    @test isapprox(Matrix(SID{1}(1, 'y')), [[0.0, -√2im/2, 0.0] [√2im/2, 0.0, -√2im/2] [0.0, √2im/2, 0.0]])
    @test isapprox(Matrix(SID{1}(1, '+')), [[0.0, √2, 0.0] [0.0, 0.0, √2] [0.0, 0.0, 0.0]])
    @test isapprox(Matrix(SID{1}(1, '-')), [[0.0, 0.0, 0.0] [√2, 0.0, 0.0] [0.0, √2, 0.0]])
end

@testset "Spin" begin
    spin = Spin{1}(norbital=2)
    @test shape(spin) == (1:2, 1:5)
    @test ndimshape(spin) == ndimshape(typeof(spin)) == 2
    @test CartesianIndex(SID{1}(1, 'z'), spin) == CartesianIndex(1, 3)
    @test SID(CartesianIndex(1, 1), spin) == SID{1}(1, 'x')
    @test summary(spin) == "10-element Spin{1}"
    @test string(spin) == "Spin{1}(norbital=2)"
    @test totalspin(spin) == totalspin(typeof(spin)) == 1
    @test collect(spin) == [
        SID{1}(1, 'x'), SID{1}(2, 'x'), SID{1}(1, 'y'), SID{1}(2, 'y'), SID{1}(1, 'z'), SID{1}(2, 'z'),
        SID{1}(1, '+'), SID{1}(2, '+'), SID{1}(1, '-'), SID{1}(2, '-')
    ]
    @test shape(IIDSpace(SID{1//2}(:a, 'x'), Spin{1//2}(2))) == (1:2, 1:1)
    @test shape(IIDSpace(SID{1//2}(2, 'y'), Spin{1//2}(2))) == (2:2, 2:2)
    @test shape(IIDSpace(SID{1//2}(:a, wildcard), Spin{1//2}(2))) == (1:2, 1:5)
    @test shape(IIDSpace(SID{1//2}(2, wildcard), Spin{1//2}(2))) == (2:2, 1:5)

    @test match(SID{wildcard}, Spin{1//2}) == true
    @test match(SID{1//2}, Spin{1//2}) == true
    @test match(SID{1//2}, Spin{1}) == match(SID{1}, Spin{1//2}) == false
end

@testset "latex" begin
    index = Index(PID(1), SID{1//2}(2, 'z'))
    @test script(Val(:site), index) == 1
    @test script(Val(:orbital), index) == 2
    @test script(Val(:tag), index) == 'z'

    @test latexname(Index{<:AbstractPID, <:SID}) == Symbol("Index{AbstractPID, SID}")
    @test latexname(CompositeOID{<:Index{<:AbstractPID, <:SID}}) == Symbol("CompositeOID{Index{AbstractPID, SID}}")
end

@testset "SpinOperator" begin
    opt = Operator(1.0, ID(
        OID(Index(PID(1), SID{1//2}('+')), [0.0, 0.0], [0.0, 0.0]),
        OID(Index(PID(1), SID{1//2}('-')), [0.0, 0.0], [0.0, 0.0])
        ))
    @test opt' == Operator(1.0, ID(
        OID(Index(PID(1), SID{1//2}('+')), [0.0, 0.0], [0.0, 0.0]),
        OID(Index(PID(1), SID{1//2}('-')), [0.0, 0.0], [0.0, 0.0])
        ))
    @test ishermitian(opt)
    @test repr(opt) == "S^{+}_{1}S^{-}_{1}"
end

@testset "permute" begin
    soptrep(opt::Operator) = opt.value * prod([Matrix(opt.id[i].index.iid) for i = 1:rank(opt)])
    for S in (1//2, 1, 3//2)
        oids = [OID(Index(PID(1), SID{S}(2, tag)), [0.0, 0.0], [0.0, 0.0]) for tag in ('x', 'y', 'z', '+', '-')]
        for (id₁, id₂) in Permutations{2}(oids)
            left = soptrep(Operator(1, ID(id₁, id₂)))
            right = sum([soptrep(opt) for opt in permute(id₁, id₂)])
            @test isapprox(left, right)
        end
    end
    id₁ = OID(Index(PID(1), SID{1//2}(2, 'z')), [0.0, 0.0], [0.0, 0.0])
    id₂ = OID(Index(PID(2), SID{1//2}(2, 'z')), [0.0, 0.0], [0.0, 0.0])
    @test permute(id₁, id₂) == (Operator(1, ID(id₂, id₁)),)
end

@testset "SpinCoupling" begin
    @test string(SpinCoupling(1.0, ('+', '-'))) == "SpinCoupling(value=1.0, tags=S⁺S⁻)"
    @test string(SpinCoupling(1.0, ('z', 'z'))) == "SpinCoupling(value=1.0, tags=SᶻSᶻ)"
    @test string(SpinCoupling(1.0, ('-', '+'), orbitals=(1, 2))) == "SpinCoupling(value=1.0, tags=S⁻S⁺, orbitals=[1 2])"
    @test repr(SpinCoupling(2.0, ('x', 'y'))) == "2.0 SˣSʸ"

    sc₁ = SpinCoupling(1.5, ('+', '-'), orbitals=subscript"[α β](α > β)")
    sc₂ = SpinCoupling(2.0, ('+', '-'), orbitals=subscript"[α β](α < β)")
    @test repr(sc₁) == "1.5 S⁺S⁻ ob[α β](α > β)"
    @test repr(sc₂) == "2.0 S⁺S⁻ ob[α β](α < β)"

    sc = sc₁ * sc₂
    @test repr(sc) == "3.0 S⁺S⁻S⁺S⁻ ob[α β](α > β)×[α β](α < β)"

    sc = SpinCoupling(2.0, ('+', '-'), orbitals=(1, 2))
    bond = Bond(1, Point(CPID(1, 2), [0.5], [0.0]), Point(CPID(1, 1), [0.0], [0.0]))
    hilbert = Hilbert(pid=>Spin{1}(norbital=2) for pid in [bond.epoint.pid, bond.spoint.pid])
    ex = expand(sc, bond, hilbert, Val(:SpinTerm))
    @test collect(ex) == [(2.0, ID(
        OID(Index(CPID(1, 1), SID{1}(1, '+')), [0.0], [0.0]),
        OID(Index(CPID(1, 2), SID{1}(2, '-')), [0.5], [0.0])
        ))]

    sc = SpinCoupling(2.0, ('+', '-', '+', '-'), orbitals=subscript"[α α β β](α < β)")
    point = Point(PID(1), [0.0], [0.0])
    hilbert = Hilbert(point.pid=>Spin{1}(norbital=3))
    ex = expand(sc, point, hilbert, Val(:info))
    @test collect(ex) == [
        (2.0, ID(OID(Index(PID(1), SID{1}(1, '+')), [0.0], [0.0]),
                 OID(Index(PID(1), SID{1}(1, '-')), [0.0], [0.0]),
                 OID(Index(PID(1), SID{1}(2, '+')), [0.0], [0.0]),
                 OID(Index(PID(1), SID{1}(2, '-')), [0.0], [0.0])
                 )),
        (2.0, ID(OID(Index(PID(1), SID{1}(1, '+')), [0.0], [0.0]),
                 OID(Index(PID(1), SID{1}(1, '-')), [0.0], [0.0]),
                 OID(Index(PID(1), SID{1}(3, '+')), [0.0], [0.0]),
                 OID(Index(PID(1), SID{1}(3, '-')), [0.0], [0.0])
                 )),
        (2.0, ID(OID(Index(PID(1), SID{1}(2, '+')), [0.0], [0.0]),
                 OID(Index(PID(1), SID{1}(2, '-')), [0.0], [0.0]),
                 OID(Index(PID(1), SID{1}(3, '+')), [0.0], [0.0]),
                 OID(Index(PID(1), SID{1}(3, '-')), [0.0], [0.0])
                 ))
    ]
end

@testset "heisenberg" begin
    @test heisenberg"+-z ob[1 2]" == heisenberg"ob[1 2]" == heisenbergpmz(orbitals=(1, 2)) == Couplings(
        SpinCoupling(1//1, ('z', 'z'), orbitals=(1, 2)),
        SpinCoupling(1//2, ('+', '-'), orbitals=(1, 2)),
        SpinCoupling(1//2, ('-', '+'), orbitals=(1, 2))
    )
    @test heisenberg"" == heisenbergpmz() == Couplings(SpinCoupling(1//1, ('z', 'z')), SpinCoupling(1//2, ('+', '-')), SpinCoupling(1//2, ('-', '+')))
    @test heisenberg"xyz" == heisenbergxyz() == Couplings(SpinCoupling(1, ('x', 'x')), SpinCoupling(1, ('y', 'y')), SpinCoupling(1, ('z', 'z')))
end

@testset "ising" begin
    @test ising"x" == Couplings(SpinCoupling(1, ('x', 'x')))
    @test ising"y" == Couplings(SpinCoupling(1, ('y', 'y')))
    @test ising"z" == Couplings(SpinCoupling(1, ('z', 'z')))

    @test ising"x ob[α β](α < β)" == Couplings(SpinCoupling(1, ('x', 'x'), orbitals=subscript"[α β](α < β)"))
    @test ising"y ob[α β](α < β)" == Couplings(SpinCoupling(1, ('y', 'y'), orbitals=subscript"[α β](α < β)"))
    @test ising"z ob[α β](α < β)" == Couplings(SpinCoupling(1, ('z', 'z'), orbitals=subscript"[α β](α < β)"))
end

@testset "gamma" begin
    @test gamma"x ob[1 1]" == gamma('y', 'z', orbitals=(1, 1))== SpinCoupling(1, ('y', 'z'), orbitals=(1, 1)) + SpinCoupling(1, ('z', 'y'), orbitals=(1, 1))
    @test gamma"y" == gamma('z', 'x') == SpinCoupling(1, ('z', 'x')) + SpinCoupling(1, ('x', 'z'))
    @test gamma"z" == gamma('x', 'y') == SpinCoupling(1, ('x', 'y')) + SpinCoupling(1, ('y', 'x'))
end

@testset "dm" begin
    @test dm"x ob[1 1]" == dm('y', 'z', orbitals=(1, 1)) == SpinCoupling(1, ('y', 'z'), orbitals=(1, 1)) - SpinCoupling(1, ('z', 'y'), orbitals=(1, 1))
    @test dm"y" == dm('z', 'x') == SpinCoupling(1, ('z', 'x')) - SpinCoupling(1, ('x', 'z'))
    @test dm"z" == dm('x', 'y') == SpinCoupling(1, ('x', 'y')) - SpinCoupling(1, ('y', 'x'))
end

@testset "sᵅ" begin
    @test sˣ"ob[1]" == Couplings(SpinCoupling(1, ('x',), orbitals=(1,)))
    @test sʸ"" == Couplings(SpinCoupling(1, ('y',)))
    @test sᶻ"" == Couplings(SpinCoupling(1, ('z',)))
end

@testset "spincoupling" begin
    sc = sc"1.0 S⁺S⁻ ob[α β](α < β)"
    @test repr(sc) == "1.0 S⁺S⁻ ob[α β](α < β)"

    sc = sc"1.0 S⁺S⁻ ob[α β]"
    @test repr(sc) == "1.0 S⁺S⁻ ob[α β]"

    sc = sc"1.0 S⁺S⁻ ob[1 2]"
    @test repr(sc) == "1.0 S⁺S⁻ ob[1 2]"

    sc = sc"1.0 S⁺S⁻"
    @test repr(sc) == "1.0 S⁺S⁻"
end

@testset "SpinTerm" begin
    point = Point(PID(1), (0.5, 0.5), (0.0, 0.0))
    hilbert = Hilbert(point.pid=>Spin{1//2}(norbital=2))
    term = SpinTerm(:h, 1.5, 0, sᶻ"")
    operators = Operators(
        Operator(1.5, ID(OID(Index(PID(1), SID{1//2}(1, 'z')), [0.5, 0.5], [0.0, 0.0]))),
        Operator(1.5, ID(OID(Index(PID(1), SID{1//2}(2, 'z')), [0.5, 0.5], [0.0, 0.0])))
    )
    @test term|>abbr == term|>typeof|>abbr == :sp
    @test term|>ishermitian == term|>typeof|>ishermitian == true
    @test expand(term, point, hilbert) == operators

    bond = Bond(1, Point(CPID('a', 1), (0.0, 0.0), (0.0, 0.0)), Point(CPID('b', 1), (0.5, 0.5), (0.0, 0.0)))
    hilbert = Hilbert(pid=>Spin{1//2}(norbital=2) for pid in [bond.spoint.pid, bond.epoint.pid])
    term = SpinTerm(:J, 1.5, 1, heisenberg"")
    operators = Operators(
        Operator(1.50, ID(OID(Index(CPID('b', 1), SID{1//2}(2, 'z')), [0.5, 0.5], [0.0, 0.0]), OID(Index(CPID('a', 1), SID{1//2}(2, 'z')), [0.0, 0.0], [0.0, 0.0]))),
        Operator(0.75, ID(OID(Index(CPID('b', 1), SID{1//2}(2, '-')), [0.5, 0.5], [0.0, 0.0]), OID(Index(CPID('a', 1), SID{1//2}(2, '+')), [0.0, 0.0], [0.0, 0.0]))),
        Operator(0.75, ID(OID(Index(CPID('b', 1), SID{1//2}(1, '-')), [0.5, 0.5], [0.0, 0.0]), OID(Index(CPID('a', 1), SID{1//2}(1, '+')), [0.0, 0.0], [0.0, 0.0]))),
        Operator(0.75, ID(OID(Index(CPID('b', 1), SID{1//2}(1, '+')), [0.5, 0.5], [0.0, 0.0]), OID(Index(CPID('a', 1), SID{1//2}(1, '-')), [0.0, 0.0], [0.0, 0.0]))),
        Operator(1.50, ID(OID(Index(CPID('b', 1), SID{1//2}(1, 'z')), [0.5, 0.5], [0.0, 0.0]), OID(Index(CPID('a', 1), SID{1//2}(1, 'z')), [0.0, 0.0], [0.0, 0.0]))),
        Operator(0.75, ID(OID(Index(CPID('b', 1), SID{1//2}(2, '+')), [0.5, 0.5], [0.0, 0.0]), OID(Index(CPID('a', 1), SID{1//2}(2, '-')), [0.0, 0.0], [0.0, 0.0])))
    )
    @test expand(term, bond, hilbert) == operators
end

@testset "NID" begin
    @test NID('u')' == NID('u')
    @test NID('p')' == NID('p')

    nid = NID('p')
    @test statistics(nid) == statistics(nid) == :b
end

@testset "Phonon" begin
    pn = Phonon(3)
    @test shape(pn) == (1:2, 1:3)
    @test ndimshape(pn) == ndimshape(typeof(pn)) == 2
    for i = 1:length(pn)
        @test NID(CartesianIndex(pn[i], pn), pn) == pn[i]
    end
    @test collect(pn) == [NID('u', 'x'), NID('p', 'x'), NID('u', 'y'), NID('p', 'y'), NID('u', 'z'), NID('p', 'z')]

    @test shape(IIDSpace(NID('u'), Phonon(3))) == (1:1, 1:3)
    @test shape(IIDSpace(NID('u', 'x'), Phonon(3))) == (1:1, 1:1)
    @test shape(IIDSpace(NID('u', 'y'), Phonon(3))) == (1:1, 2:2)
    @test shape(IIDSpace(NID('u', 'z'), Phonon(3))) == (1:1, 3:3)

    @test shape(IIDSpace(NID('p'), Phonon(2))) == (2:2, 1:2)
    @test shape(IIDSpace(NID('p', 'x'), Phonon(3))) == (2:2, 1:1)
    @test shape(IIDSpace(NID('p', 'y'), Phonon(3))) == (2:2, 2:2)
    @test shape(IIDSpace(NID('p', 'z'), Phonon(3))) == (2:2, 3:3)
end

@testset "latex" begin
    index = Index(PID(1), NID('u', 'x'))
    @test script(Val(:BD), index, nlatex) == 'u'
    @test script(Val(:BD), OID(index, [0.0, 0.0], [0.0, 0.0]), nlatex) == 'u'
    @test script(Val(:site), index) == 1
    @test script(Val(:dir), index) == 'x'

    index = Index(PID(2), NID('p', 'y'))
    @test script(Val(:BD), index, nlatex) == 'p'
    @test script(Val(:BD), OID(index, [0.0, 0.0], [0.0, 0.0]), nlatex) == 'p'
    @test script(Val(:site), index) == 2
    @test script(Val(:dir), index) == 'y'

    latexname(Index{<:AbstractPID, <:NID}) == Symbol("Index{AbstractPID, NID}")
    latexname(CompositeOID{<:Index{<:AbstractPID, <:NID}}) == Symbol("CompositeOID{Index{AbstractPID, NID}}")
end

@testset "PhononOperator" begin
    opt = Operator(1.0, ID(
        OID(Index(PID(1), NID('p', 'x')), [0.0, 0.0], [0.0, 0.0]),
        OID(Index(PID(1), NID('p', 'x')), [0.0, 0.0], [0.0, 0.0])
        ))
    @test opt' == Operator(1.0, ID(
        OID(Index(PID(1), NID('p', 'x')), [0.0, 0.0], [0.0, 0.0]),
        OID(Index(PID(1), NID('p', 'x')), [0.0, 0.0], [0.0, 0.0])
        ))
    @test ishermitian(opt)
    @test repr(opt) == "(p^{x}_{1})^2"
end

@testset "permute" begin
    id₁ = OID(Index(PID(1), NID('u', 'x')), [0.0, 0.0], [0.0, 0.0])
    id₂ = OID(Index(PID(1), NID('p', 'x')), [0.0, 0.0], [0.0, 0.0])
    @test permute(id₁, id₂) == (Operator(+1im), Operator(1, ID(id₂, id₁)))
    @test permute(id₂, id₁) == (Operator(-1im), Operator(1, ID(id₁, id₂)))

    id₁ = OID(Index(PID(1), NID('u', 'x')), [0.0, 0.0], [0.0, 0.0])
    id₂ = OID(Index(PID(1), NID('u', 'x')), [0.0, 0.0], [0.0, 0.0])
    @test permute(id₁, id₂) == (Operator(1, ID(id₂, id₁)),)

    id₁ = OID(Index(PID(1), NID('u', 'x')), [0.0, 0.0], [0.0, 0.0])
    id₂ = OID(Index(PID(1), NID('p', 'y')), [0.0, 0.0], [0.0, 0.0])
    @test permute(id₁, id₂) == (Operator(1, ID(id₂, id₁)),)
end

@testset "PhononCoupling" begin
    pnc = PhononCoupling(1.0, ('p', 'p'))
    @test string(pnc) == "PhononCoupling(value=1.0, tags=[p p])"
    @test repr(pnc) == "1.0 [p p]"

    pnc = PhononCoupling(1.0, ('p', 'p'), dirs=('x', 'x'))
    @test string(pnc) == "PhononCoupling(value=1.0, tags=[p p], dirs=[x x])"
    @test repr(pnc) == "1.0 [p p] dr[x x]"

    pnc = PhononCoupling(2.0, ('p', 'p'))
    point = Point(PID(1), [0.5, 0.0], [0.0, 0.0])
    hilbert = Hilbert(point.pid=>Phonon(2))
    ex = expand(pnc, point, hilbert, Val(:PhononKinetic))
    @test collect(ex) == [
        (2.0, ID(OID(Index(PID(1), NID('p', 'x')), [0.5, 0.0], [0.0, 0.0]), OID(Index(PID(1), NID('p', 'x')), [0.5, 0.0], [0.0, 0.0]))),
        (2.0, ID(OID(Index(PID(1), NID('p', 'y')), [0.5, 0.0], [0.0, 0.0]), OID(Index(PID(1), NID('p', 'y')), [0.5, 0.0], [0.0, 0.0])))
    ]

    pnc = PhononCoupling(1.0, ('u', 'u'))
    bond = Bond(1, Point(PID(2), [0.5, 0.0], [0.0, 0.0]), Point(PID(1), [0.0, 0.0], [0.0, 0.0]))
    hilbert = Hilbert(pid=>Phonon(2) for pid in [bond.epoint.pid, bond.spoint.pid])
    ex = expand(pnc, bond, hilbert, Val(:PhononPotential))
    @test shape(ex) == (1:2, 1:3)
    @test collect(ex) ==[
        (+1.0, ID(OID(Index(PID(1), NID('u', 'x')), [0.0, 0.0], [0.0, 0.0]), OID(Index(PID(1), NID('u', 'x')), [0.0, 0.0], [0.0, 0.0]))),
        (+0.0, ID(OID(Index(PID(1), NID('u', 'y')), [0.0, 0.0], [0.0, 0.0]), OID(Index(PID(1), NID('u', 'y')), [0.0, 0.0], [0.0, 0.0]))),
        (-2.0, ID(OID(Index(PID(1), NID('u', 'x')), [0.0, 0.0], [0.0, 0.0]), OID(Index(PID(2), NID('u', 'x')), [0.5, 0.0], [0.0, 0.0]))),
        (+0.0, ID(OID(Index(PID(1), NID('u', 'y')), [0.0, 0.0], [0.0, 0.0]), OID(Index(PID(2), NID('u', 'y')), [0.5, 0.0], [0.0, 0.0]))),
        (+1.0, ID(OID(Index(PID(2), NID('u', 'x')), [0.5, 0.0], [0.0, 0.0]), OID(Index(PID(2), NID('u', 'x')), [0.5, 0.0], [0.0, 0.0]))),
        (+0.0, ID(OID(Index(PID(2), NID('u', 'y')), [0.5, 0.0], [0.0, 0.0]), OID(Index(PID(2), NID('u', 'y')), [0.5, 0.0], [0.0, 0.0])))
        ]

    @test kinetic"" == Couplings(PhononCoupling(1, ('p', 'p')))
    @test potential"" == Couplings(PhononCoupling(1, ('u', 'u')))
end

@testset "PhononKinetic" begin
    term = PhononKinetic(:T, 2.0)
    @test abbr(term) == abbr(typeof(term)) == :pnk
    @test ishermitian(term) == ishermitian(typeof(term)) == true

    point = Point(PID(1), [0.5, 0.0], [0.0, 0.0])
    hilbert = Hilbert(point.pid=>Phonon(2))
    operators = Operators(
        Operator(2.0, ID(OID(Index(PID(1), NID('p', 'x')), [0.5, 0.0], [0.0, 0.0]), OID(Index(PID(1), NID('p', 'x')), [0.5, 0.0], [0.0, 0.0]))),
        Operator(2.0, ID(OID(Index(PID(1), NID('p', 'y')), [0.5, 0.0], [0.0, 0.0]), OID(Index(PID(1), NID('p', 'y')), [0.5, 0.0], [0.0, 0.0])))
    )
    @test expand(term, point, hilbert) == operators
end

@testset "PhononPotential" begin
    term = PhononPotential(:V, 2.0, 1)
    @test abbr(term) == abbr(typeof(term)) == :pnp
    @test ishermitian(term) == ishermitian(typeof(term)) == true

    bond = Bond(1, Point(PID(2), [0.5, 0.0], [0.0, 0.0]), Point(PID(1), [0.0, 0.0], [0.0, 0.0]))
    hilbert = Hilbert(pid=>Phonon(2) for pid in [bond.epoint.pid, bond.spoint.pid])
    operators = Operators(
        Operator(+2.0, ID(OID(Index(PID(1), NID('u', 'x')), [0.0, 0.0], [0.0, 0.0]), OID(Index(PID(1), NID('u', 'x')), [0.0, 0.0], [0.0, 0.0]))),
        Operator(-4.0, ID(OID(Index(PID(1), NID('u', 'x')), [0.0, 0.0], [0.0, 0.0]), OID(Index(PID(2), NID('u', 'x')), [0.5, 0.0], [0.0, 0.0]))),
        Operator(+2.0, ID(OID(Index(PID(2), NID('u', 'x')), [0.5, 0.0], [0.0, 0.0]), OID(Index(PID(2), NID('u', 'x')), [0.5, 0.0], [0.0, 0.0]))),
    )
    @test expand(term, bond, hilbert) == operators

    bond = Bond(1, Point(PID(2), [0.0, 0.5], [0.0, 0.0]), Point(PID(1), [0.0, 0.0], [0.0, 0.0]))
    hilbert = Hilbert(pid=>Phonon(2) for pid in [bond.epoint.pid, bond.spoint.pid])
    operators = Operators(
        Operator(+2.0, ID(OID(Index(PID(1), NID('u', 'y')), [0.0, 0.0], [0.0, 0.0]), OID(Index(PID(1), NID('u', 'y')), [0.0, 0.0], [0.0, 0.0]))),
        Operator(-4.0, ID(OID(Index(PID(1), NID('u', 'y')), [0.0, 0.0], [0.0, 0.0]), OID(Index(PID(2), NID('u', 'y')), [0.0, 0.5], [0.0, 0.0]))),
        Operator(+2.0, ID(OID(Index(PID(2), NID('u', 'y')), [0.0, 0.5], [0.0, 0.0]), OID(Index(PID(2), NID('u', 'y')), [0.0, 0.5], [0.0, 0.0])))
    )
    @test expand(term, bond, hilbert) == operators

    bond = Bond(1, Point(PID(2), [0.5, 0.5], [0.0, 0.0]), Point(PID(1), [0.0, 0.0], [0.0, 0.0]))
    hilbert = Hilbert(pid=>Phonon(2) for pid in [bond.epoint.pid, bond.spoint.pid])
    operators = Operators(
        Operator(+1.0, ID(OID(Index(PID(1), NID('u', 'x')), [0.0, 0.0], [0.0, 0.0]), OID(Index(PID(1), NID('u', 'x')), [0.0, 0.0], [0.0, 0.0]))),
        Operator(+1.0, ID(OID(Index(PID(1), NID('u', 'y')), [0.0, 0.0], [0.0, 0.0]), OID(Index(PID(1), NID('u', 'y')), [0.0, 0.0], [0.0, 0.0]))),
        Operator(-2.0, ID(OID(Index(PID(1), NID('u', 'x')), [0.0, 0.0], [0.0, 0.0]), OID(Index(PID(2), NID('u', 'x')), [0.5, 0.5], [0.0, 0.0]))),
        Operator(-2.0, ID(OID(Index(PID(1), NID('u', 'y')), [0.0, 0.0], [0.0, 0.0]), OID(Index(PID(2), NID('u', 'y')), [0.5, 0.5], [0.0, 0.0]))),
        Operator(+1.0, ID(OID(Index(PID(2), NID('u', 'x')), [0.5, 0.5], [0.0, 0.0]), OID(Index(PID(2), NID('u', 'x')), [0.5, 0.5], [0.0, 0.0]))),
        Operator(+1.0, ID(OID(Index(PID(2), NID('u', 'y')), [0.5, 0.5], [0.0, 0.0]), OID(Index(PID(2), NID('u', 'y')), [0.5, 0.5], [0.0, 0.0])))
    )
    @test expand(term, bond, hilbert) ≈ operators
end

@testset "DMPhonon" begin
    @test dmphonon"" == Couplings(Coupling(1, ID(NID('u', wildcard), SID{wildcard}(1, wildcard))))

    term = DMPhonon(:dmp, 2.0, 1)
    @test abbr(term) == abbr(typeof(term)) == :dmp
    @test ishermitian(term) == ishermitian(typeof(term)) == false

    bond = Bond(1, Point(PID(2), [0.5, 0.5], [0.0, 0.0]), Point(PID(1), [0.0, 0.0], [0.0, 0.0]))
    hilbert = Hilbert(pid=>Phonon(2)⊗Spin{1//2}(1) for pid in [bond.epoint.pid, bond.spoint.pid])
    operators = Operators(
        Operator(√2/2, ID(OID(Index(PID(1), SID{1//2}(1, 'y')), [0.0, 0.0], [0.0, 0.0]), OID(Index(PID(2), NID('u', 'x')), [0.5, 0.5], [0.0, 0.0]))),
        Operator(√2/2, ID(OID(Index(PID(2), NID('u', 'x')), [0.5, 0.5], [0.0, 0.0]), OID(Index(PID(2), SID{1//2}(1, 'y')), [0.5, 0.5], [0.0, 0.0]))),
        Operator(√2/2, ID(OID(Index(PID(1), NID('u', 'x')), [0.0, 0.0], [0.0, 0.0]), OID(Index(PID(2), SID{1//2}(1, 'x')), [0.5, 0.5], [0.0, 0.0]))),
        Operator(√2/2, ID(OID(Index(PID(1), NID('u', 'x')), [0.0, 0.0], [0.0, 0.0]), OID(Index(PID(1), SID{1//2}(1, 'x')), [0.0, 0.0], [0.0, 0.0]))),
        Operator(-√2/2, ID(OID(Index(PID(1), SID{1//2}(1, 'x')), [0.0, 0.0], [0.0, 0.0]), OID(Index(PID(2), NID('u', 'x')), [0.5, 0.5], [0.0, 0.0]))),
        Operator(-√2/2, ID(OID(Index(PID(1), SID{1//2}(1, 'y')), [0.0, 0.0], [0.0, 0.0]), OID(Index(PID(1), NID('u', 'x')), [0.0, 0.0], [0.0, 0.0]))),
        Operator(-√2/2, ID(OID(Index(PID(1), SID{1//2}(1, 'x')), [0.0, 0.0], [0.0, 0.0]), OID(Index(PID(1), NID('u', 'y')), [0.0, 0.0], [0.0, 0.0]))),
        Operator(√2/2, ID(OID(Index(PID(2), NID('u', 'x')), [0.5, 0.5], [0.0, 0.0]), OID(Index(PID(1), SID{1//2}(1, 'y')), [0.0, 0.0], [0.0, 0.0]))),
        Operator(√2/2, ID(OID(Index(PID(1), NID('u', 'y')), [0.0, 0.0], [0.0, 0.0]), OID(Index(PID(1), SID{1//2}(1, 'y')), [0.0, 0.0], [0.0, 0.0]))),
        Operator(-√2/2, ID(OID(Index(PID(1), NID('u', 'y')), [0.0, 0.0], [0.0, 0.0]), OID(Index(PID(1), SID{1//2}(1, 'x')), [0.0, 0.0], [0.0, 0.0]))),
        Operator(√2/2, ID(OID(Index(PID(2), SID{1//2}(1, 'y')), [0.5, 0.5], [0.0, 0.0]), OID(Index(PID(1), NID('u', 'y')), [0.0, 0.0], [0.0, 0.0]))),
        Operator(-√2/2, ID(OID(Index(PID(2), SID{1//2}(1, 'x')), [0.5, 0.5], [0.0, 0.0]), OID(Index(PID(2), NID('u', 'x')), [0.5, 0.5], [0.0, 0.0]))),
        Operator(√2/2, ID(OID(Index(PID(2), NID('u', 'y')), [0.5, 0.5], [0.0, 0.0]), OID(Index(PID(2), SID{1//2}(1, 'x')), [0.5, 0.5], [0.0, 0.0]))),
        Operator(√2/2, ID(OID(Index(PID(1), SID{1//2}(1, 'y')), [0.0, 0.0], [0.0, 0.0]), OID(Index(PID(1), NID('u', 'y')), [0.0, 0.0], [0.0, 0.0]))),
        Operator(-√2/2, ID(OID(Index(PID(1), NID('u', 'x')), [0.0, 0.0], [0.0, 0.0]), OID(Index(PID(2), SID{1//2}(1, 'y')), [0.5, 0.5], [0.0, 0.0]))),
        Operator(-√2/2, ID(OID(Index(PID(2), SID{1//2}(1, 'y')), [0.5, 0.5], [0.0, 0.0]), OID(Index(PID(1), NID('u', 'x')), [0.0, 0.0], [0.0, 0.0]))),
        Operator(-√2/2, ID(OID(Index(PID(2), NID('u', 'y')), [0.5, 0.5], [0.0, 0.0]), OID(Index(PID(1), SID{1//2}(1, 'y')), [0.0, 0.0], [0.0, 0.0]))),
        Operator(√2/2, ID(OID(Index(PID(2), SID{1//2}(1, 'y')), [0.5, 0.5], [0.0, 0.0]), OID(Index(PID(2), NID('u', 'x')), [0.5, 0.5], [0.0, 0.0]))),
        Operator(√2/2, ID(OID(Index(PID(1), NID('u', 'y')), [0.0, 0.0], [0.0, 0.0]), OID(Index(PID(2), SID{1//2}(1, 'y')), [0.5, 0.5], [0.0, 0.0]))),
        Operator(√2/2, ID(OID(Index(PID(2), NID('u', 'y')), [0.5, 0.5], [0.0, 0.0]), OID(Index(PID(1), SID{1//2}(1, 'x')), [0.0, 0.0], [0.0, 0.0]))),
        Operator(-√2/2, ID(OID(Index(PID(2), SID{1//2}(1, 'x')), [0.5, 0.5], [0.0, 0.0]), OID(Index(PID(1), NID('u', 'y')), [0.0, 0.0], [0.0, 0.0]))),
        Operator(-√2/2, ID(OID(Index(PID(1), SID{1//2}(1, 'y')), [0.0, 0.0], [0.0, 0.0]), OID(Index(PID(2), NID('u', 'y')), [0.5, 0.5], [0.0, 0.0]))),
        Operator(√2/2, ID(OID(Index(PID(2), SID{1//2}(1, 'x')), [0.5, 0.5], [0.0, 0.0]), OID(Index(PID(2), NID('u', 'y')), [0.5, 0.5], [0.0, 0.0]))),
        Operator(√2/2, ID(OID(Index(PID(2), SID{1//2}(1, 'x')), [0.5, 0.5], [0.0, 0.0]), OID(Index(PID(1), NID('u', 'x')), [0.0, 0.0], [0.0, 0.0]))),
        Operator(-√2/2, ID(OID(Index(PID(1), NID('u', 'x')), [0.0, 0.0], [0.0, 0.0]), OID(Index(PID(1), SID{1//2}(1, 'y')), [0.0, 0.0], [0.0, 0.0]))),
        Operator(-√2/2, ID(OID(Index(PID(2), NID('u', 'x')), [0.5, 0.5], [0.0, 0.0]), OID(Index(PID(2), SID{1//2}(1, 'x')), [0.5, 0.5], [0.0, 0.0]))),
        Operator(-√2/2, ID(OID(Index(PID(2), SID{1//2}(1, 'y')), [0.5, 0.5], [0.0, 0.0]), OID(Index(PID(2), NID('u', 'y')), [0.5, 0.5], [0.0, 0.0]))),
        Operator(√2/2, ID(OID(Index(PID(1), SID{1//2}(1, 'x')), [0.0, 0.0], [0.0, 0.0]), OID(Index(PID(2), NID('u', 'y')), [0.5, 0.5], [0.0, 0.0]))),
        Operator(-√2/2, ID(OID(Index(PID(2), NID('u', 'y')), [0.5, 0.5], [0.0, 0.0]), OID(Index(PID(2), SID{1//2}(1, 'y')), [0.5, 0.5], [0.0, 0.0]))),
        Operator(√2/2, ID(OID(Index(PID(1), SID{1//2}(1, 'x')), [0.0, 0.0], [0.0, 0.0]), OID(Index(PID(1), NID('u', 'x')), [0.0, 0.0], [0.0, 0.0]))),
        Operator(-√2/2, ID(OID(Index(PID(2), NID('u', 'x')), [0.5, 0.5], [0.0, 0.0]), OID(Index(PID(1), SID{1//2}(1, 'x')), [0.0, 0.0], [0.0, 0.0]))),
        Operator(-√2/2, ID(OID(Index(PID(1), NID('u', 'y')), [0.0, 0.0], [0.0, 0.0]), OID(Index(PID(2), SID{1//2}(1, 'x')), [0.5, 0.5], [0.0, 0.0])))
    )
    @test expand(term, bond, hilbert) ≈ operators
end


text:
"""
    reproject(input_data, output_projection; shape_out = nothing, order = 1, hdu_in = 1, hdu_out = 1)

Reprojects image data to a new projection using interpolation.

# Arguments
- `input_data`: Image data which is being reprojected.
                It can be an ImageHDU, FITS object, name of a FITS file or a tuple of image matrix and WCSTransform.
- `output_projection`: Frame in which data is reprojected.
                       Frame can be taken from WCSTransform object, ImageHDU, FITS or name of FITS file.
- `shape_out`: Shape of image after reprojection.
- `order`: Order of interpolation.
           0: Nearest-neighbor
           1: Linear
           2: Quadratic
- `hdu_in`: Used to specify HDU number when giving input as FITS or name of FITS file.
- `hud_out:` Used to specify HDU number when giving output projection as FITS or name of FITS file.
"""
function reproject(input_data, output_projection; shape_out = nothing, order::Int = 1, hdu_in::Int = 1, hdu_out::Int = 1)
    if input_data isa ImageHDU || input_data isa Tuple{AbstractArray, WCSTransform}
        array_in, wcs_out = parse_input_data(input_data)
    else
        array_in, wcs_out = parse_input_data(input_data, hdu_in)
    end

    if output_projection isa FITS || output_projection isa String
        wcs_in, shape_out = parse_output_projection(output_projection, hdu_out)
    else
        wcs_in, shape_out = parse_output_projection(output_projection, shape_out)
    end

    type_in = wcs_to_celestial_frame(wcs_in)
    type_out = wcs_to_celestial_frame(wcs_out)

    if type_in == type_out && shape_out === nothing
        return array_in
    end

    img_out = fill(NaN, shape_out)
    array_in = pad_edges(array_in)
    itp = interpolator(array_in, order)
    shape_in = size(array_in)

    for i in 1:shape_out[1]
        for j in 1:shape_out[2]
            pix_coord_in = [float(i), float(j)]
            world_coord_in = pix_to_world(wcs_in, pix_coord_in)

            if type_in == "ICRS"
                coord_in = ICRSCoords(deg2rad(world_coord_in[1]), deg2rad(world_coord_in[2]))
            elseif type_in == "Gal"
                coord_in = GalCoords(deg2rad(world_coord_in[1]), deg2rad(world_coord_in[2]))
            elseif type_in == "FK5"
                coord_in = FK5Coords{wcs_in.equinox}(deg2rad(world_coord_in[1]), deg2rad(world_coord_in[2]))
            else
                throw(ArgumentError("Unsupported output WCS coordinate type"))
            end

            if type_out == "ICRS"
                coord_out = convert(ICRSCoords, coord_in)
            elseif type_out == "Gal"
                coord_out = convert(GalCoords, coord_in)
            elseif type_out == "FK5"
                coord_out = convert(FK5Coords{wcs_out.equinox}, coord_in)
            else
                throw(ArgumentError("Unsupported input WCS coordinate type"))
            end

            pix_coord_out = world_to_pix(wcs_out, [rad2deg(lon(coord_out)), rad2deg(lat(coord_out))])

            if 0.5 <= pix_coord_out[1] <= shape_in[1] - 1 && 0.5 <= pix_coord_out[2] <= shape_in[2] - 1
                img_out[i,j] = itp(pix_coord_out[1] + 1, pix_coord_out[2] + 1)
            end
        end
    end

    return img_out, (!isnan).(img_out)
end

"""
    interpolator(array_in, order::Int)

Returns an interpolator with the given array and order of interpolation.
"""
function interpolator(array_in::AbstractArray, order::Int)
   if order == 0
        itp = interpolate(array_in, BSpline(Constant()))
    elseif order == 1
        itp = interpolate(array_in, BSpline(Linear()))
    else
        itp = interpolate(array_in, BSpline(Quadratic(InPlace(OnCell()))))
    end

    return itp
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module TestGivens

using Test, LinearAlgebra, Random
using LinearAlgebra: rmul!, lmul!, Givens

# Test givens rotations
@testset for elty in (Float32, Float64, ComplexF32, ComplexF64)
    if elty <: Real
        raw_A = convert(Matrix{elty}, randn(10,10))
    else
        raw_A = convert(Matrix{elty}, complex.(randn(10,10),randn(10,10)))
    end
    @testset for A in (raw_A, view(raw_A, 1:10, 1:10))
        Ac = copy(A)
        R = LinearAlgebra.Rotation(LinearAlgebra.Givens{elty}[])
        for j = 1:8
            for i = j+2:10
                G, _ = givens(A, j+1, i, j)
                lmul!(G, A)
                rmul!(A, adjoint(G))
                lmul!(G, R)

                @test lmul!(G,Matrix{elty}(I, 10, 10)) == [G[i,j] for i=1:10,j=1:10]

                @testset "transposes" begin
                    @test G'*G*Matrix(elty(1)I, 10, 10)   ≈ Matrix(I, 10, 10)
                    @test (G*Matrix(elty(1)I, 10, 10))*G' ≈ Matrix(I, 10, 10)
                    @test copy(R')*(R*Matrix(elty(1)I, 10, 10)) ≈ Matrix(I, 10, 10)
                    @test_throws ErrorException transpose(G)
                    @test_throws ErrorException transpose(R)
                end
            end
        end
        @test_throws ArgumentError givens(A, 3, 3, 2)
        @test_throws ArgumentError givens(one(elty),zero(elty),2,2)
        G, _ = givens(one(elty),zero(elty),11,12)
        @test_throws DimensionMismatch lmul!(G, A)
        @test_throws DimensionMismatch rmul!(A, adjoint(G))
        @test abs.(A) ≈ abs.(hessenberg(Ac).H)
        @test opnorm(R*Matrix{elty}(I, 10, 10)) ≈ one(elty)

        I10 = Matrix{elty}(I, 10, 10)
        G, _ = givens(one(elty),zero(elty),9,10)
        @test (G*I10)' * (G*I10) ≈ I10
        K, _ = givens(zero(elty),one(elty),9,10)
        @test (K*I10)' * (K*I10) ≈ I10

        @testset "Givens * vectors" begin
            if isa(A, Array)
                x = A[:, 1]
            else
                x = view(A, 1:10, 1)
            end
            G, r = givens(x[2], x[4], 2, 4)
            @test (G*x)[2] ≈ r
            @test abs((G*x)[4]) < eps(real(elty))
            @inferred givens(x[2], x[4], 2, 4)

            G, r = givens(x, 2, 4)
            @test (G*x)[2] ≈ r
            @test abs((G*x)[4]) < eps(real(elty))
            @inferred givens(x, 2, 4)

            G, r = givens(x, 4, 2)
            @test (G*x)[4] ≈ r
            @test abs((G*x)[2]) < eps(real(elty))
        end
    end
end

# 36430
# dimensional correctness:
const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test")
isdefined(Main, :Furlongs) || @eval Main include(joinpath($(BASE_TEST_PATH), "testhelpers", "Furlongs.jl"))
using .Main.Furlongs

@testset "testing dimensions with Furlongs" begin
    @test_throws MethodError givens(Furlong(1.0), Furlong(2.0), 1, 2)
end

const TNumber = Union{Float64,ComplexF64}
struct MockUnitful{T<:TNumber} <: Number
    data::T
    MockUnitful(data::T) where T<:TNumber = new{T}(data)
end
import Base: *, /, one, oneunit
*(a::MockUnitful{T}, b::T) where T<:TNumber = MockUnitful(a.data * b)
*(a::T, b::MockUnitful{T}) where T<:TNumber = MockUnitful(a * b.data)
*(a::MockUnitful{T}, b::MockUnitful{T}) where T<:TNumber = MockUnitful(a.data * b.data)
/(a::MockUnitful{T}, b::MockUnitful{T}) where T<:TNumber = a.data / b.data
one(::Type{<:MockUnitful{T}}) where T = one(T)
oneunit(::Type{<:MockUnitful{T}}) where T = MockUnitful(one(T))

@testset "unitful givens rotation unitful $T " for T in (Float64, ComplexF64)
    g, r = givens(MockUnitful(T(3)), MockUnitful(T(4)), 1, 2)
    @test g.c ≈ 3/5
    @test g.s ≈ 4/5
    @test r.data ≈ 5.0
end

end # module TestGivens

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
using RKHSRegularization
using Test

@testset "RKHSRegularization.jl" begin
    # Write your tests here.
end

max_stars_repo_name: RoyCCWang/RKHSRegularization

max_issues_repo_name: RoyCCWang/RKHSRegularization

max_forks_repo_name: RoyCCWang/RKHSRegularization


text:
function main()

    io = open("krypton.txt", "w");
    B1 = 26102.88e-8
    B2 = 56946.82e-8
    λ1 = 14.36*10^(-9)
    λ2 = 89.81*10^(-9)
    C1 = 2.01e-6*(1e-6)^2
    C2 = 10.043e-3*(1e-6)^2

    p0 = 1
    k0(λ) = 2π/λ
    a = 14.9e-6
    T = 273+20
    T0 = 273
    λ0 = 800e-9
    δ(λ) = T0 / T * (B1 * λ^2/(λ^2 - C1) + B2 * λ^2 / (λ^2 - C2))
    n(λ, p, u) = sqrt(1 + δ(λ)*p/p0 - (u / k0(λ) / a)^2 )

    e11 = 2.405
    e12 = 5.52
    e13 = 8.6573
    e23 = 10.1735
    e33 = 11.6198
    e43 = 13.0152

    ps = range(0, 12, length=1000)

    for p in ps
        res13 = (n(λ0, p, e11) - n(λ0/3, p, e13)) * 10^4
        res23 = (n(λ0, p, e11) - n(λ0/3, p, e23)) * 10^4
        res33 = (n(λ0, p, e11) - n(λ0/3, p, e33)) * 10^4
        write(io, "$p \t $res13 \t $res23 \t $res33 \n")
    end
    close(io)
end
main()


text:
@testset "Configure using $method" for method in [
    "function arguments", "environment variables", "toml file"
]
    user_id = "user1"
    api_key = "12345"

    if method == "function arguments"
        config = configure(user_id=user_id, api_key=api_key; banyanconfig_path="tempfile.toml")
    elseif method == "environment variables"
        ENV["BANYAN_USER_ID"] = user_id
        ENV["BANYAN_API_KEY"] = api_key
        config = configure(; banyanconfig_path="tempfile.toml")
    elseif method == "toml file"
        delete!(ENV, "BANYAN_USER_ID")
        delete!(ENV, "BANYAN_API_KEY")
        config = configure(; banyanconfig_path="tempfile.toml")
    end
    println(config)
    @test config["banyan"]["user_id"] == user_id && config["banyan"]["api_key"] == api_key

    rm("tempfile.toml", force=true)
end

max_stars_repo_name: banyan-team/banyan-julia

max_issues_repo_name: banyan-team/banyan-julia

max_forks_repo_name: banyan-team/banyan-julia


text:
function get_dm(auth::TwAuth, options::Dict)
  Requests.json(get_oauth(auth, "https://api.twitter.com/1.1/direct_messages.json", options))
end

function get_dm_sent(auth::TwAuth, options::Dict)
  Requests.json(get_oauth(auth, "https://api.twitter.com/1.1/direct_messages/sent.json", options))
end

function get_dm_show(auth::TwAuth, options::Dict)
  Requests.json(get_oauth(auth, "https://api.twitter.com/1.1/direct_messages/show.json", options))
end

max_stars_repo_name: nick-paul/JuliaTwitter

max_issues_repo_name: nick-paul/JuliaTwitter

max_forks_repo_name: nick-paul/JuliaTwitter


text:
struct FBPAFFT{F<:AbstractCTFilter} <: AbstractFBP
    filter::F
end

FBPAFFT() = FBPAFFT(RamLak())
FBPAFFT(f::Function) = FBPAFFT(CTFilter(f))

_alg_method(::FBPAFFT) = IsIRadonDiag()
@_defiradonalgfn FBPAFFT fbpa_fft

function fbpa_fft end

@_defiradonfn fbpa_fft begin
    xys = Vector{NTuple{2,T}}(undef, rows * cols)
    δx::T = width(xs)
    δy::T = width(ys)
    h::T = hypot(δx, δy)
    λ::T = (nd - 1) / h
    κ::T = λ / ν
    @inbounds @simd for k ∈ eachindex(xys)
        x = xs[(k - 1) ÷ rows + 1] * κ
        y = ys[(k - 1) % rows + 1] * κ
        xys[k] = x, y
    end
    p = _iradon_progress(length(tomog), progress)
    Threads.@threads for k ∈ eachindex(xys)
        x, y = xys[k]
        @inbounds tomog[k] = sum(eachindex(scϕs)) do iϕ
            sϕ, cϕ = scϕs[iϕ]
            # To be consistent with our conventions should be '+'.
            t = (x * cϕ + y * sϕ + 1) * t₀
            t ∈ 1..nd ? interp(t, iϕ) : z
        end
        next!(p)
    end
    γ::T = min(rows, cols) / hypot(rows, cols) * rows / cols
    δt::T = π * γ / length(scϕs) / nd * λ^2
    tomog .*= δt
end

# function fbpa_fft(
#     sinog::AbstractMatrix{T},
#     xs::AbstractVector{U1},
#     ys::AbstractVector{U2},
#     ϕs::ClosedInterval = 0..2π;
#     background::Optional{U3} = nothing,
#     filter::Optional{F} = nothing,
#     interpolation::Optional{Interp} = nothing,
#     progress::Bool = true,
# ) where {
#     T <: Real,
#     U1 <: Real,
#     U2 <: Real,
#     U3 <: Real,
#     F <: AbstractCTFilter,
#     Interp <: AbstractInterp2DOrNone,
# }
#     nd, nϕ = size(sinog)
#     cols = length(xs)
#     rows = length(ys)
#     filtered = apply(maybe(Filters.RamLak(), filter)) do f
#         filter_freq = fft(sinog, 1) .* f(T, nd, nϕ)
#         ifft(filter_freq, 1) |> real
#     end
#     interpolation = maybe(interpolate, interpolation)
#     interp = interpolation(filtered)
#     t₀::T = (nd + 1) / 2
#     xys = Vector{NTuple{2,T}}(undef, rows * cols)
#     @inbounds @simd for k ∈ eachindex(xys)
#         xys[k] = xs[(k - 1) ÷ rows + 1], ys[(k - 1) % rows + 1]
#     end
#     scϕs = sincos.(linspace(T, ORI(ϕs), nϕ))
#     z::T = maybe(zero(T), background)
#     tomog = similar(sinog, rows, cols)
#     fill!(tomog, z)
#     p = _iradon_progress(length(tomog), progress)
#     Threads.@threads for k ∈ eachindex(xys)
#         x, y = xys[k]
#         @inbounds tomog[k] = sum(eachindex(scϕs)) do iϕ
#             sϕ, cϕ = scϕs[iϕ]
#             # To be consistent with our conventions should be '+'.
#             t = x * cϕ + y * sϕ + t₀
#             t ∈ 1..nd ? interp(t, iϕ) : z
#         end
#         next!(p)
#     end
#     CTTomogram(tomog)
# end


text:
"""
# module Scriptify

# Examples

```jldoctest
julia>
```
"""
module Scriptify

using Slurm.Commands
using Slurm.SystemModules
using Slurm.Shells

export scriptify

scriptify(s::Sbatch) = join(["#$(prefix(s)) $(string(dir, long = true, with_value = true))" for dir in s.directives], '\n')
scriptify(modules::AbstractVector{<: SystemModule}) = join(["module load $(string(m))" for m in modules], '\n')
scriptify(s::Shell) = string(s.path)

end


text:
using BernoulliFactory
using Random
import Random.GLOBAL_RNG
using Test

function makef(p::Float64, rng::RNG=GLOBAL_RNG) where RNG<:AbstractRNG
  @assert 0.0 <= p <= 1.0
  function f()
    return rand(rng) < p
  end
  return f
end

function testFunction(f::F1, m::Int64, v::Float64, fp::F2,
  args::Vararg{Any, N}) where {F1<:Function, F2<:Function, N}
  sumCoins::Int64 = 0
  for i in 1:m
    sumCoins += f(fp, args...)[1]
  end
  if v == 1.0 || v == 0.0
    @test sumCoins == m * v
  else
    @test abs((sumCoins/m - v))/sqrt(v*(1-v)/m) < 3
  end
end

ntrials = 100000

@testset "linear tests" begin
  @time include("linear_test.jl")
end

@testset "inverse tests" begin
  @time include("inverse_test.jl")
end

@testset "power tests" begin
  @time include("power_test.jl")
end

@testset "expminus tests" begin
  @time include("expminus_test.jl")
end

@testset "logistic tests" begin
  @time include("logistic_test.jl")
end

@testset "two coin tests" begin
  @time include("twocoin_test.jl")
end

@testset "signed estimate tests" begin
  @time include("signedEstimate_test.jl")
end


text:
export _zero, @fields

# update a field of a struct.
@inline @generated function field_update(main :: T, field::Val{Field}, value) where {T, Field}
    fields = fieldnames(T)
    Expr(:new, T, Any[field !== Field ? :(main.$field) : :value for field in fields]...)
end

# the default constructor of a struct
@inline @generated function default_constructor(::Type{T}, fields::Vararg{Any,N}) where {T,N}
    Expr(:new, T, Any[:(fields[$i]) for i=1:N]...)
end

"""
    _zero(T)
    _zero(x::T)

Create a `zero` of type `T` by recursively applying `zero` to its fields.
"""
@inline @generated function _zero(::Type{T}) where {T}
    Expr(:new, T, Any[:(_zero($field)) for field in T.types]...)
end

@inline @generated function _zero(x::T) where {T}
    Expr(:new, T, Any[:(_zero(x.$field)) for field in fieldnames(T)]...)
end

function lens_compile(ex, cache, value)
    @smatch ex begin
        :($a.$b.$c = $d) => begin
            updated =
                Expr(:let,
                    Expr(:block, :($cache = $cache.$b), :($value = $d)),
                    :($field_update($cache, $(Val(c)), $value)))
            lens_compile(:($a.$b = $updated), cache, value)
        end
        :($a.$b = $c) => begin
            Expr(:let,
                Expr(:block, :($cache = $a), :($value=$c)),
                :($field_update($cache, $(Val(b)), $value)))
        end
        _ => error("Malformed update notation $ex, expect the form like 'a.b = c'.")
    end
end

function with(ex)
    cache = gensym("cache")
    value = gensym("value")
    lens_compile(ex, cache, value)
end

"""
e.g. `@with x.y = val` will return a new object similar to `x`, with the `y` field changed to `val`.
"""
macro with(ex)
    with(ex) |> esc
end

@inline @generated function _zero(::Type{T}) where {T<:Tuple}
    Expr(:tuple, Any[:(_zero($field)) for field in T.types]...)
end
_zero(::Type{T}) where T<:Real = zero(T)
_zero(::Type{String}) = ""
_zero(::Type{Symbol}) = Symbol("")
_zero(::Type{Char}) = '\0'
_zero(::Type{T}) where {ET,N,T<:Array{ET,N}} = reshape(ET[], ntuple(x->0, N))
_zero(::Type{T}) where {A,B,T<:Dict{A,B}} = Dict{A,B}()

#_zero(x::T) where T = _zero(T) # not adding this line!
_zero(x::T) where T<:Real = zero(x)
_zero(::String) = ""
_zero(::Symbol) = Symbol("")
_zero(::Char) = '\0'
_zero(x::T) where T<:Array = zero(x)
function _zero(d::T) where {A,B,T<:Dict{A,B}}
    Dict{A,B}([x=>zero(y) for (x,y) in d])
end

@static if VERSION > v"1.6.100"
    @generated function chfield(x, ::Val{FIELD}, xval) where FIELD
        if ismutabletype(x)
            Expr(:block, :(x.$FIELD = xval), :x)
        else
            :(@with x.$FIELD = xval)
        end
    end
else
    @generated function chfield(x, ::Val{FIELD}, xval) where FIELD
        if x.mutable
            Expr(:block, :(x.$FIELD = xval), :x)
        else
            :(@with x.$FIELD = xval)
        end
    end
end
@generated function chfield(x, f::Function, xval)
    Expr(:block, _invcheck(:(f(x)), :xval), :x)
end

# convert field of an object to a tuple
@generated function type2tuple(x::T) where T
    Expr(:tuple, [:(x.$v) for v in fieldnames(T)]...)
end

macro fields(ex)
    esc(:($type2tuple($ex)))
end


text:
println("\n\ncheck_records.jl: start")

prepend!(LOAD_PATH, ["Project.toml"])
import LibPQ
using DataFrames, PrettyTables

DISASTER_PASSWORD = read("./etc/disaster-pass", String)

postgres = LibPQ.Connection("host=127.0.0.1 port=5432 dbname=disaster user=disaster password=$DISASTER_PASSWORD")

sql_queries = [
    """
    select 's_earthquake: record_count=' ||
        (select count(*) from public.s_earthquake) ||
        ', with_ogc_fid_count=' ||
        (select count(*) from public.s_earthquake where ogc_fid is not null);
    """,
    """
    select 's_storm: record_count=' ||
        (select count(*) from public.s_storm) ||
        ', with_ogc_fid_count=' ||
        (select count(*) from public.s_storm where ogc_fid is not null);
    """,
    """
    select 's_wildfire: record_count=' ||
        (select count(*) from public.s_wildfire) ||
        ', with_ogc_fid_count=' ||
        (select count(*) from public.s_wildfire where ogc_fid is not null);
    """,
    """
    select 'total: record_count=' || (
        (select count(*) from public.f_earthquake) +
        (select count(*) from public.f_storm) +
        (select count(*) from public.f_wildfire)
    );
    """,
    """
    select 'f_earthquake: record_count=' ||
        (select count(*) from public.f_earthquake);
    """,
    """
    select 'f_storm: record_count=' ||
        (select count(*) from public.f_storm);
    """,
    """
    select 'f_wildfire: record_count=' ||
        (select count(*) from public.f_wildfire);
    """,
    """
    select distinct(d_severity_id), count(1)
    from public.f_earthquake
    group by d_severity_id
    order by d_severity_id;
    """,
    """
    select distinct(d_severity_id), count(1)
    from public.f_wildfire
    group by d_severity_id
    order by d_severity_id;
    """,
    """
    select distinct(d_severity_id), count(1)
    from public.f_storm
    group by d_severity_id
    order by d_severity_id;
    """,
]

for sql_query in sql_queries
    result = LibPQ.execute(postgres, sql_query) |> DataFrame
    pretty_table(result)
end

LibPQ.close(postgres)

println("\ncheck_records: stop")