12369488912

Committed 17 Dec 2024 08:51AM UTC coverage: 94.306% (-0.3%) from 94.567%

Build # 12369488912

Build Type

Pull #300

github

Committed by

mtsch

Commit Message

fix doctests

Pull Request Pull Request #300: Fix normalisation in `single_particle_density`

Run Details

3 of 3 new or added lines in 1 file covered. (100.0%)

97 existing lines in 22 files now uncovered.

6907 of 7324 relevant lines covered (94.31%)

13567973.06 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

92.36

/src/ExactDiagonalization/basis_set_representation.jl

"""
    build_sparse_matrix_from_LO(
        ham, address=starting_address(ham);
        cutoff, filter=nothing, nnzs, col_hint, sort=false, kwargs...
    ) -> sparse_matrix, basis
    build_sparse_matrix_from_LO(ham, addresses::AbstractVector; kwargs...)

Create a sparse matrix `sparse_matrix` of all reachable matrix elements of a linear operator `ham`
starting from `address`. Instead of a single address, a vector of `addresses` can be passed.
The vector `basis` contains the addresses of basis configurations.

Providing the number `nnzs` of expected calculated matrix elements and `col_hint` for the
estimated number of nonzero off-diagonal matrix elements in each matrix column may improve
performance.

Providing an energy cutoff will skip the columns and rows with diagonal elements greater
than `cutoff`. Alternatively, an arbitrary `filter` function can be used instead. These are
not enabled by default. To generate the matrix truncated to the subspace spanned by the
`addresses`, use `filter = Returns(false)`.

Setting `sort` to `true` will sort the `basis` and order the matrix rows and columns
accordingly. This is useful when the order of the columns matters, e.g. when comparing
matrices. Any additional keyword arguments are passed on to `Base.sortperm`.

See [`BasisSetRepresentation`](@ref).
"""
function build_sparse_matrix_from_LO(
    ham, addr_or_vec=starting_address(ham);
    cutoff=nothing,
    filter=isnothing(cutoff) ? nothing : (a -> diagonal_element(ham, a) ≤ cutoff),
    nnzs=0, col_hint=0, # sizehints are opt-in
    sort=false, kwargs...
)
    # Set up `adds` as queue of addresses. Also returned as the basis.
    adds = addr_or_vec isa Union{AbstractArray,Tuple} ? [addr_or_vec...] : [addr_or_vec]

    T = eltype(ham)
    dict = Dict(add => i for (i, add) in enumerate(adds)) # Map from addresses to indices
    col = Dict{Int,T}()       # Temporary column storage
    sizehint!(col, col_hint)

    is = Int[] # row indices
    js = Int[] # column indice
    vs = T[]   # non-zero values

    sizehint!(is, nnzs)
    sizehint!(js, nnzs)
    sizehint!(vs, nnzs)

    i = 0
    while i < length(adds)
        i += 1
        add = adds[i]
        push!(is, i)
        push!(js, i)
        push!(vs, diagonal_element(ham, add))

        empty!(col)
        for (off, v) in offdiagonals(ham, add)
            iszero(v) && continue
            j = get(dict, off, nothing)
            if isnothing(j)
                # Energy cutoff: remember skipped addresses, but avoid adding them to `adds`
                if !isnothing(filter) && !filter(off)
                    dict[off] = 0
                    j = 0
                else
                    push!(adds, off)
                    j = length(adds)
                    dict[off] = j
                end
            end
            if !iszero(j)
                col[j] = get(col, j, zero(T)) + v
            end
        end
        # Copy the column into the sparse matrix vectors.
        for (j, v) in col
            iszero(v) && continue
            push!(is, i)
            push!(js, j)
            push!(vs, v)
        end
    end

    matrix = sparse(js, is, vs, length(adds), length(adds))
    if sort
        perm = sortperm(adds; kwargs...)
        return permute!(matrix, perm, perm), permute!(adds, perm)
    else
        return matrix, adds
    end
end

"""
    build_basis(
        ham, address=starting_address(ham);
        cutoff, filter, sizelim, sort=false, kwargs...
    ) -> basis
    build_basis(ham, addresses::AbstractVector; kwargs...)

Get all basis element of a linear operator `ham` that are reachable (via
non-zero matrix elements) from the address `address`, returned as a vector.
Instead of a single address, a vector of `addresses` can be passed.
Does not return the matrix, for that purpose use [`BasisSetRepresentation`](@ref).

Providing an energy cutoff will skip addresses with diagonal elements greater
than `cutoff`. Alternatively, an arbitrary `filter` function can be used instead.
Addresses passed as arguments are not filtered.
A maximum basis size `sizelim` can be set which will throw an error if the expected
dimension of `ham` is larger than `sizelim`. This may be useful when memory may be a
concern. These options are disabled by default.

Setting `sort` to `true` will sort the basis. Any additional keyword arguments
are passed on to `Base.sort!`.
"""
function build_basis(
    ham, addr_or_vec=starting_address(ham);
    cutoff=nothing,
    filter=isnothing(cutoff) ? nothing : (a -> diagonal_element(ham, a) ≤ cutoff),
    sort=false,
    max_size=Inf, # retained for backwards compatibility; use sizelim instead
    sizelim=max_size,
    kwargs...
)
    check_address_type(ham, addr_or_vec)
    single_addr = addr_or_vec isa Union{AbstractArray,Tuple} ? addr_or_vec[1] : addr_or_vec
    if dimension(ham, single_addr) > sizelim
        throw(ArgumentError("dimension larger than sizelim"))
    end
    # Set up `adds` as queue of addresses. Also returned as the basis.
    adds = addr_or_vec isa Union{AbstractArray,Tuple} ? [addr_or_vec...] : [addr_or_vec]
    known_basis = Set(adds)     # known addresses

    i = 0
    while i < length(adds)
        i += 1
        add = adds[i]

        for (off, v) in offdiagonals(ham, add)
            (iszero(v) || off in known_basis) && continue   # check if valid
            push!(known_basis, off)
            !isnothing(filter) && !filter(off) && continue  # check filter
            push!(adds, off)
        end
    end

    if sort
        return sort!(adds, kwargs...)
    else
        return adds
    end
end

"""
    BasisSetRepresentation(
        hamiltonian::AbstractHamiltonian, addr=starting_address(hamiltonian);
        sizelim=10^6, nnzs, cutoff, filter, sort=false, kwargs...
    )
    BasisSetRepresentation(hamiltonian::AbstractHamiltonian, addresses::AbstractVector; kwargs...)

Eagerly construct the basis set representation of the operator `hamiltonian` with all addresses
reachable from `addr`. Instead of a single address, a vector of `addresses` can be passed.

An `ArgumentError` is thrown if `dimension(hamiltonian) > sizelim` in order to prevent memory
overflow. Set `sizelim = Inf` in order to disable this behaviour.

Providing the number `nnzs` of expected calculated matrix elements and `col_hint` for the
estimated number of nonzero off-diagonal matrix elements in each matrix column may improve
performance.

Providing an energy cutoff will skip the columns and rows with diagonal elements greater
than `cutoff`. Alternatively, an arbitrary `filter` function can be used instead.
Addresses passed as arguments are not filtered. To generate the matrix truncated to the
subspace spanned by the `addresses`, use `filter = Returns(false)`.

Setting `sort` to `true` will sort the matrix rows and columns. This is useful when the
order of the columns matters, e.g. when comparing matrices. Any additional keyword arguments
are passed on to `Base.sortperm`.

## Fields
* `sparse_matrix`: sparse matrix representing `hamiltonian` in the basis `basis`
* `basis`: vector of addresses
* `hamiltonian`: the Hamiltonian `hamiltonian`

## Example
```jldoctest; filter = r"(\\d*)\\.(\\d{4})\\d+" => s"\\1.\\2***"
julia> hamiltonian = HubbardReal1D(BoseFS(1,0,0));

julia> bsr = BasisSetRepresentation(hamiltonian)
BasisSetRepresentation(HubbardReal1D(fs"|1 0 0⟩"; u=1.0, t=1.0)) with dimension 3 and 9 stored entries:3×3 SparseArrays.SparseMatrixCSC{Float64, Int64} with 9 stored entries:
  0.0  -1.0  -1.0
 -1.0   0.0  -1.0
 -1.0  -1.0   0.0

julia> BasisSetRepresentation(hamiltonian, bsr.basis[1:2]; filter = Returns(false)) # passing addresses and truncating
BasisSetRepresentation(HubbardReal1D(fs"|1 0 0⟩"; u=1.0, t=1.0)) with dimension 2 and 4 stored entries:2×2 SparseArrays.SparseMatrixCSC{Float64, Int64} with 4 stored entries:
  0.0  -1.0
 -1.0   0.0

julia> using LinearAlgebra; round.(eigvals(Matrix(bsr)); digits = 4) # eigenvalues
3-element Vector{Float64}:
 -2.0
  1.0
  1.0

julia> ev = eigvecs(Matrix(bsr))[:,1]; ev = ev .* sign(ev[1]) # ground state eigenvector
3-element Vector{Float64}:
 0.5773502691896257
 0.5773502691896255
 0.5773502691896257

julia> dv = DVec(zip(bsr.basis, ev)) # ground state as DVec
DVec{BoseFS{1, 3, BitString{3, 1, UInt8}},Float64} with 3 entries, style = IsDeterministic{Float64}()
  fs"|0 0 1⟩" => 0.57735
  fs"|0 1 0⟩" => 0.57735
  fs"|1 0 0⟩" => 0.57735
```
Has methods for [`dimension`](@ref), [`sparse`](@ref), [`Matrix`](@ref),
[`starting_address`](@ref).

Part of the [`AbstractHamiltonian`](@ref) interface. See also [`build_basis`](@ref).
"""
struct BasisSetRepresentation{A,SM,H}
    sparse_matrix::SM
    basis::Vector{A}
    hamiltonian::H
end

function BasisSetRepresentation(
    hamiltonian::AbstractHamiltonian, addr_or_vec=starting_address(hamiltonian);
    kwargs...
)
    # In the default case we pass `AdjointUnknown()` in order to skip the
    # symmetrisation of the sparse matrix
    return _bsr_ensure_symmetry(AdjointUnknown(), hamiltonian, addr_or_vec; kwargs...)
end
# special cases are needed for symmetry wrappers

function BasisSetRepresentation(
    hamiltonian::ParitySymmetry, addr=starting_address(hamiltonian);
    kwargs...
)
    # For symmetry wrappers it is necessary to explicity symmetrise the matrix to
    # avoid the loss of matrix symmetry due to floating point rounding errors
    return _bsr_ensure_symmetry(LOStructure(hamiltonian), hamiltonian, addr; kwargs...)
end

function BasisSetRepresentation(
    hamiltonian::TimeReversalSymmetry, addr=starting_address(hamiltonian);
    kwargs...
)
    # For symmetry wrappers it is necessary to explicity symmetrise the matrix to
    # avoid the loss of matrix symmetry due to floating point rounding errors
    return _bsr_ensure_symmetry(LOStructure(hamiltonian), hamiltonian, addr; kwargs...)
end

# default, does not enforce symmetries
function _bsr_ensure_symmetry(
    ::LOStructure, hamiltonian::AbstractHamiltonian, addr_or_vec;
    sizelim=10^6, test_approx_symmetry=true, kwargs...
)
    single_addr = addr_or_vec isa Union{AbstractArray,Tuple} ? addr_or_vec[1] : addr_or_vec
    d = dimension(hamiltonian, single_addr)
    if d > sizelim
        throw(ArgumentError("Dimension = $d larger than sizelim = $sizelim. Set a larger `sizelim` if this is safe."))
    end
    check_address_type(hamiltonian, addr_or_vec)
    sparse_matrix, basis = build_sparse_matrix_from_LO(hamiltonian, addr_or_vec; kwargs...)
    return BasisSetRepresentation(sparse_matrix, basis, hamiltonian)
end

# build the BasisSetRepresentation while enforcing hermitian symmetry
function _bsr_ensure_symmetry(
    ::IsHermitian, hamiltonian::AbstractHamiltonian, addr_or_vec;
    sizelim=10^6, test_approx_symmetry=true, kwargs...
)
    single_addr = addr_or_vec isa Union{AbstractArray,Tuple} ? addr_or_vec[1] : addr_or_vec
    if dimension(hamiltonian, single_addr) > sizelim
        throw(ArgumentError("dimension larger than sizelim"))
    end
    check_address_type(hamiltonian, addr_or_vec)
    sparse_matrix, basis = build_sparse_matrix_from_LO(hamiltonian, addr_or_vec; kwargs...)
    fix_approx_hermitian!(sparse_matrix; test_approx_symmetry) # enforce hermitian symmetry after building
    return BasisSetRepresentation(sparse_matrix, basis, hamiltonian)
end

"""
    fix_approx_hermitian!(A; test_approx_symmetry=true, kwargs...)
Replaces the matrix `A` by `½(A + A')` in place. This will be successful and the result
is guaranteed to pass the `ishermitian` test only if the matrix is square and already
approximately hermitian.

By default logs an error message if the matrix `A` is found to not be approximately
hermitian. Set `test_approx_symmetry=false` to bypass testing.
Other keyword arguments are passed on to `isapprox`.
"""
function fix_approx_hermitian!(A; test_approx_symmetry=true, kwargs...)
    # Generic and inefficient version. Make sure to replace by efficient specialised code.
    if test_approx_symmetry
        passed = isapprox(A, A'; kwargs...)
        if !passed
            throw(ArgumentError("Matrix is not approximately hermitian."))
        end
    end
    @. A = 1 / 2 * (A + A')
    return A
end

# specialised code for sparse matrices
using SparseArrays: AbstractSparseMatrixCSC, getcolptr, rowvals, nonzeros

# special case for sparse matrices; avoids most allocations, testing is free
function fix_approx_hermitian!(A::AbstractSparseMatrixCSC; test_approx_symmetry=false, kwargs...)
    passed = isapprox_enforce_hermitian!(A; kwargs...)
    if test_approx_symmetry && !passed
        throw(ArgumentError("Matrix is not approximately hermitian."))
    end
    return A
end

"""
    isapprox_enforce_hermitian!(A::AbstractSparseMatrixCSC; kwargs...) -> Bool
Checks whether the matrix `A` is approximately hermitian by checking each pair of transposed
matrix elements with `isapprox`. Keyword arguments are passed on to `isapprox`.
Returns boolean `true` is the test is passed and `false` if not.

Furthermore, the matrix `A` is modified to become exactly equal to `½(A + A')` if the test
is passed.
"""
function isapprox_enforce_hermitian!(A::AbstractSparseMatrixCSC; kwargs...)
    # based on `ishermsym()` from `SparseArrays`; relies on `SparseArrays` internals
    # https://github.com/JuliaSparse/SparseArrays.jl/blob/1bae96dc8f9a8ca8b7879eef4cf71e186598e982/src/sparsematrix.jl#L3793
    m, n = size(A)
    if m != n
        return false
    end

    colptr = getcolptr(A)
    rowval = rowvals(A)
    nzval = nonzeros(A)
    tracker = copy(getcolptr(A))
    for col = 1:size(A, 2)
        # `tracker` is updated such that, for symmetric matrices,
        # the loop below starts from an element at or below the
        # diagonal element of column `col`"
        for p = tracker[col]:colptr[col+1]-1
            val = nzval[p]
            row = rowval[p]

            # Ignore stored zeros
            if iszero(val)
                continue
            end

            # If the matrix was symmetric we should have updated
            # the tracker to start at the diagonal or below. Here
            # we are above the diagonal so the matrix can't be symmetric.
            if row < col
                return false
            end

            # Diagonal element
            if row == col
                if isapprox(val, conj(val); kwargs...)
                    nzval[p] = real(val)
                else
                    return false
                end
                # if val != conj(val)
                #     return false
                # end
            else
                offset = tracker[row]

                # If the matrix is unsymmetric, there might not exist
                # a rowval[offset]
                if offset > length(rowval)
                    return false
                end

                row2 = rowval[offset]

                # row2 can be less than col if the tracker didn't
                # get updated due to stored zeros in previous elements.
                # We therefore "catch up" here while making sure that
                # the elements are actually zero.
                while row2 < col
                    if !iszero(nzval[offset])
                        return false
                    end
                    offset += 1
                    row2 = rowval[offset]
                    tracker[row] += 1
                end

                # Non zero A[i,j] exists but A[j,i] does not exist
                if row2 > col
                    return false
                end

                # A[i,j] and A[j,i] exists
                if row2 == col
                    if isapprox(val, conj(nzval[offset]); kwargs...)
                        val = 1 / 2 * (val + conj(nzval[offset]))
                        nzval[p] = val
                        nzval[offset] = conj(val)
                    else
                        return false
                    end
                    # if val != conj(nzval[offset])
                    #     return false
                    # end
                    tracker[row] += 1
                end
            end
        end
    end
    return true
end

function Base.show(io::IO, b::BasisSetRepresentation)
    print(io, "BasisSetRepresentation($(b.hamiltonian)) with dimension $(dimension(b))" *
        " and $(nnz(b.sparse_matrix)) stored entries:")
    show(io, MIME"text/plain"(), b.sparse_matrix)
end

Interfaces.starting_address(bsr::BasisSetRepresentation) = bsr.basis[1]

Hamiltonians.dimension(bsr::BasisSetRepresentation) = length(bsr.basis)

"""
    sparse(hamiltonian::AbstractHamiltonian, addr=starting_address(hamiltonian); kwargs...)
    sparse(bsr::BasisSetRepresentation)

Return a sparse matrix representation of `hamiltonian` or `bsr`. `kwargs` are passed to
[`BasisSetRepresentation`](@ref).

See [`BasisSetRepresentation`](@ref).
"""
function SparseArrays.sparse(hamiltonian::AbstractHamiltonian, args...; kwargs...)
    return sparse(BasisSetRepresentation(hamiltonian, args...; kwargs...))
end
SparseArrays.sparse(bsr::BasisSetRepresentation) = bsr.sparse_matrix

"""
    Matrix(
        hamiltonian::AbstractHamiltonian, addr=starting_address(hamiltonian);
        sizelim=10^4, kwargs...
    )
    Matrix(bsr::BasisSetRepresentation)

Return a dense matrix representation of `hamiltonian` or `bsr`. `kwargs` are passed to
[`BasisSetRepresentation`](@ref).

See [`BasisSetRepresentation`](@ref).
"""
function Base.Matrix(hamiltonian::AbstractHamiltonian, args...; sizelim=1e4, kwargs...)
    return Matrix(BasisSetRepresentation(hamiltonian, args...; sizelim, kwargs...))
end
Base.Matrix(bsr::BasisSetRepresentation) = Matrix(bsr.sparse_matrix)

1	"""
2	build_sparse_matrix_from_LO(
3	ham, address=starting_address(ham);
4	cutoff, filter=nothing, nnzs, col_hint, sort=false, kwargs...
5	) -> sparse_matrix, basis
6	build_sparse_matrix_from_LO(ham, addresses::AbstractVector; kwargs...)
7
8	Create a sparse matrix `sparse_matrix` of all reachable matrix elements of a linear operator `ham`
9	starting from `address`. Instead of a single address, a vector of `addresses` can be passed.
10	The vector `basis` contains the addresses of basis configurations.
11
12	Providing the number `nnzs` of expected calculated matrix elements and `col_hint` for the
13	estimated number of nonzero off-diagonal matrix elements in each matrix column may improve
14	performance.
15
16	Providing an energy cutoff will skip the columns and rows with diagonal elements greater
17	than `cutoff`. Alternatively, an arbitrary `filter` function can be used instead. These are
18	not enabled by default. To generate the matrix truncated to the subspace spanned by the
19	`addresses`, use `filter = Returns(false)`.
20
21	Setting `sort` to `true` will sort the `basis` and order the matrix rows and columns
22	accordingly. This is useful when the order of the columns matters, e.g. when comparing
23	matrices. Any additional keyword arguments are passed on to `Base.sortperm`.
24
25	See [`BasisSetRepresentation`](@ref).
26	"""
27	function build_sparse_matrix_from_LO(	287✔
28	ham, addr_or_vec=starting_address(ham);
29	cutoff=nothing,
30	filter=isnothing(cutoff) ? nothing : (a -> diagonal_element(ham, a) ≤ cutoff),	128✔
31	nnzs=0, col_hint=0, # sizehints are opt-in
32	sort=false, kwargs...
33	)
34	# Set up `adds` as queue of addresses. Also returned as the basis.
35	adds = addr_or_vec isa Union{AbstractArray,Tuple} ? [addr_or_vec...] : [addr_or_vec]	143✔
36
37	T = eltype(ham)	143✔
38	dict = Dict(add => i for (i, add) in enumerate(adds)) # Map from addresses to indices	143✔
39	col = Dict{Int,T}() # Temporary column storage	143✔
40	sizehint!(col, col_hint)	286✔
41
42	is = Int[] # row indices	143✔
43	js = Int[] # column indice	143✔
44	vs = T[] # non-zero values	143✔
45
46	sizehint!(is, nnzs)	143✔
47	sizehint!(js, nnzs)	143✔
48	sizehint!(vs, nnzs)	143✔
49
50	i = 0	143✔
51	while i < length(adds)	6,389✔
52	i += 1	6,246✔
53	add = adds[i]	6,246✔
54	push!(is, i)	6,246✔
55	push!(js, i)	6,246✔
56	push!(vs, diagonal_element(ham, add))	6,556✔
57
58	empty!(col)	460,320✔
59	for (off, v) in offdiagonals(ham, add)	11,543✔
60	iszero(v) && continue	1,076,683✔
61	j = get(dict, off, nothing)	340,949✔
62	if isnothing(j)	582,950✔
63	# Energy cutoff: remember skipped addresses, but avoid adding them to `adds`
64	if !isnothing(filter) && !filter(off)	6,123✔
65	dict[off] = 0	27✔
66	j = 0	27✔
67	else
68	push!(adds, off)	6,090✔
69	j = length(adds)	6,090✔
70	dict[off] = j	6,090✔
71	end
72	end
73	if !iszero(j)	294,665✔
74	col[j] = get(col, j, zero(T)) + v	294,611✔
75	end
76	end	2,146,440✔
77	# Copy the column into the sparse matrix vectors.
78	for (j, v) in col	12,489✔
79	iszero(v) && continue	109,217✔
80	push!(is, i)	109,189✔
81	push!(js, j)	109,189✔
82	push!(vs, v)	109,189✔
83	end	212,191✔
84	end	6,246✔
85
86	matrix = sparse(js, is, vs, length(adds), length(adds))	143✔
87	if sort	143✔
88	perm = sortperm(adds; kwargs...)	13✔
89	return permute!(matrix, perm, perm), permute!(adds, perm)	13✔
90	else
91	return matrix, adds	130✔
92	end
93	end
94
95	"""
96	build_basis(
97	ham, address=starting_address(ham);
98	cutoff, filter, sizelim, sort=false, kwargs...
99	) -> basis
100	build_basis(ham, addresses::AbstractVector; kwargs...)
101
102	Get all basis element of a linear operator `ham` that are reachable (via
103	non-zero matrix elements) from the address `address`, returned as a vector.
104	Instead of a single address, a vector of `addresses` can be passed.
105	Does not return the matrix, for that purpose use [`BasisSetRepresentation`](@ref).
106
107	Providing an energy cutoff will skip addresses with diagonal elements greater
108	than `cutoff`. Alternatively, an arbitrary `filter` function can be used instead.
109	Addresses passed as arguments are not filtered.
110	A maximum basis size `sizelim` can be set which will throw an error if the expected
111	dimension of `ham` is larger than `sizelim`. This may be useful when memory may be a
112	concern. These options are disabled by default.
113
114	Setting `sort` to `true` will sort the basis. Any additional keyword arguments
115	are passed on to `Base.sort!`.
116	"""
117	function build_basis(	305✔
118	ham, addr_or_vec=starting_address(ham);
119	cutoff=nothing,
120	filter=isnothing(cutoff) ? nothing : (a -> diagonal_element(ham, a) ≤ cutoff),	130✔
121	sort=false,
122	max_size=Inf, # retained for backwards compatibility; use sizelim instead
123	sizelim=max_size,
124	kwargs...
125	)
126	check_address_type(ham, addr_or_vec)	277✔
127	single_addr = addr_or_vec isa Union{AbstractArray,Tuple} ? addr_or_vec[1] : addr_or_vec	151✔
128	if dimension(ham, single_addr) > sizelim	301✔
129	throw(ArgumentError("dimension larger than sizelim"))	1✔
130	end
131	# Set up `adds` as queue of addresses. Also returned as the basis.
132	adds = addr_or_vec isa Union{AbstractArray,Tuple} ? [addr_or_vec...] : [addr_or_vec]	150✔
133	known_basis = Set(adds) # known addresses	150✔
134
135	i = 0	150✔
136	while i < length(adds)	1,151✔
137	i += 1	1,001✔
138	add = adds[i]	1,001✔
139
140	for (off, v) in offdiagonals(ham, add)	1,825✔
141	(iszero(v) \|\| off in known_basis) && continue # check if valid	21,111✔
142	push!(known_basis, off)	756✔
143	!isnothing(filter) && !filter(off) && continue # check filter	756✔
144	push!(adds, off)	726✔
145	end	14,174✔
146	end	1,001✔
147
148	if sort	150✔
149	return sort!(adds, kwargs...)	1✔
150	else
151	return adds	149✔
152	end
153	end
154
155	"""
156	BasisSetRepresentation(
157	hamiltonian::AbstractHamiltonian, addr=starting_address(hamiltonian);
158	sizelim=10^6, nnzs, cutoff, filter, sort=false, kwargs...
159	)
160	BasisSetRepresentation(hamiltonian::AbstractHamiltonian, addresses::AbstractVector; kwargs...)
161
162	Eagerly construct the basis set representation of the operator `hamiltonian` with all addresses
163	reachable from `addr`. Instead of a single address, a vector of `addresses` can be passed.
164
165	An `ArgumentError` is thrown if `dimension(hamiltonian) > sizelim` in order to prevent memory
166	overflow. Set `sizelim = Inf` in order to disable this behaviour.
167
168	Providing the number `nnzs` of expected calculated matrix elements and `col_hint` for the
169	estimated number of nonzero off-diagonal matrix elements in each matrix column may improve
170	performance.
171
172	Providing an energy cutoff will skip the columns and rows with diagonal elements greater
173	than `cutoff`. Alternatively, an arbitrary `filter` function can be used instead.
174	Addresses passed as arguments are not filtered. To generate the matrix truncated to the
175	subspace spanned by the `addresses`, use `filter = Returns(false)`.
176
177	Setting `sort` to `true` will sort the matrix rows and columns. This is useful when the
178	order of the columns matters, e.g. when comparing matrices. Any additional keyword arguments
179	are passed on to `Base.sortperm`.
180
181	## Fields
182	* `sparse_matrix`: sparse matrix representing `hamiltonian` in the basis `basis`
183	* `basis`: vector of addresses
184	* `hamiltonian`: the Hamiltonian `hamiltonian`
185
186	## Example
187	```jldoctest; filter = r"(\\d)\\.(\\d{4})\\d+" => s"\\1.\\2**"
188	julia> hamiltonian = HubbardReal1D(BoseFS(1,0,0));
189
190	julia> bsr = BasisSetRepresentation(hamiltonian)
191	BasisSetRepresentation(HubbardReal1D(fs"\|1 0 0⟩"; u=1.0, t=1.0)) with dimension 3 and 9 stored entries:3×3 SparseArrays.SparseMatrixCSC{Float64, Int64} with 9 stored entries:
192	0.0 -1.0 -1.0
193	-1.0 0.0 -1.0
194	-1.0 -1.0 0.0
195
196	julia> BasisSetRepresentation(hamiltonian, bsr.basis[1:2]; filter = Returns(false)) # passing addresses and truncating
197	BasisSetRepresentation(HubbardReal1D(fs"\|1 0 0⟩"; u=1.0, t=1.0)) with dimension 2 and 4 stored entries:2×2 SparseArrays.SparseMatrixCSC{Float64, Int64} with 4 stored entries:
198	0.0 -1.0
199	-1.0 0.0
200
201	julia> using LinearAlgebra; round.(eigvals(Matrix(bsr)); digits = 4) # eigenvalues
202	3-element Vector{Float64}:
203	-2.0
204	1.0
205	1.0
206
207	julia> ev = eigvecs(Matrix(bsr))[:,1]; ev = ev .* sign(ev[1]) # ground state eigenvector
208	3-element Vector{Float64}:
209	0.5773502691896257
210	0.5773502691896255
211	0.5773502691896257
212
213	julia> dv = DVec(zip(bsr.basis, ev)) # ground state as DVec
214	DVec{BoseFS{1, 3, BitString{3, 1, UInt8}},Float64} with 3 entries, style = IsDeterministic{Float64}()
215	fs"\|0 0 1⟩" => 0.57735
216	fs"\|0 1 0⟩" => 0.57735
217	fs"\|1 0 0⟩" => 0.57735
218	```
219	Has methods for [`dimension`](@ref), [`sparse`](@ref), [`Matrix`](@ref),
220	[`starting_address`](@ref).
221
222	Part of the [`AbstractHamiltonian`](@ref) interface. See also [`build_basis`](@ref).
223	"""
224	struct BasisSetRepresentation{A,SM,H}
225	sparse_matrix::SM	140✔
226	basis::Vector{A}
227	hamiltonian::H
228	end
229
230	function BasisSetRepresentation(	313✔
231	hamiltonian::AbstractHamiltonian, addr_or_vec=starting_address(hamiltonian);
232	kwargs...
233	)
234	# In the default case we pass `AdjointUnknown()` in order to skip the
235	# symmetrisation of the sparse matrix
236	return _bsr_ensure_symmetry(AdjointUnknown(), hamiltonian, addr_or_vec; kwargs...)	120✔
237	end
238	# special cases are needed for symmetry wrappers
239
240	function BasisSetRepresentation(	39✔
241	hamiltonian::ParitySymmetry, addr=starting_address(hamiltonian);
242	kwargs...
243	)
244	# For symmetry wrappers it is necessary to explicity symmetrise the matrix to
245	# avoid the loss of matrix symmetry due to floating point rounding errors
246	return _bsr_ensure_symmetry(LOStructure(hamiltonian), hamiltonian, addr; kwargs...)	13✔
247	end
248
249	function BasisSetRepresentation(	32✔
250	hamiltonian::TimeReversalSymmetry, addr=starting_address(hamiltonian);
251	kwargs...
252	)
253	# For symmetry wrappers it is necessary to explicity symmetrise the matrix to
254	# avoid the loss of matrix symmetry due to floating point rounding errors
255	return _bsr_ensure_symmetry(LOStructure(hamiltonian), hamiltonian, addr; kwargs...)	9✔
256	end
257
258	# default, does not enforce symmetries
259	function _bsr_ensure_symmetry(	240✔
260	::LOStructure, hamiltonian::AbstractHamiltonian, addr_or_vec;
261	sizelim=10^6, test_approx_symmetry=true, kwargs...
262	)
263	single_addr = addr_or_vec isa Union{AbstractArray,Tuple} ? addr_or_vec[1] : addr_or_vec	120✔
264	d = dimension(hamiltonian, single_addr)	120✔
265	if d > sizelim	161✔
266	throw(ArgumentError("Dimension = $d larger than sizelim = $sizelim. Set a larger `sizelim` if this is safe."))	1✔
267	end
268	check_address_type(hamiltonian, addr_or_vec)	132✔
269	sparse_matrix, basis = build_sparse_matrix_from_LO(hamiltonian, addr_or_vec; kwargs...)	118✔
270	return BasisSetRepresentation(sparse_matrix, basis, hamiltonian)	118✔
271	end
272
273	# build the BasisSetRepresentation while enforcing hermitian symmetry
274	function _bsr_ensure_symmetry(	44✔
275	::IsHermitian, hamiltonian::AbstractHamiltonian, addr_or_vec;
276	sizelim=10^6, test_approx_symmetry=true, kwargs...
277	)
278	single_addr = addr_or_vec isa Union{AbstractArray,Tuple} ? addr_or_vec[1] : addr_or_vec	22✔
279	if dimension(hamiltonian, single_addr) > sizelim	39✔
UNCOV 280	throw(ArgumentError("dimension larger than sizelim"))	×
281	end
282	check_address_type(hamiltonian, addr_or_vec)	30✔
283	sparse_matrix, basis = build_sparse_matrix_from_LO(hamiltonian, addr_or_vec; kwargs...)	22✔
284	fix_approx_hermitian!(sparse_matrix; test_approx_symmetry) # enforce hermitian symmetry after building	22✔
285	return BasisSetRepresentation(sparse_matrix, basis, hamiltonian)	22✔
286	end
287
288	"""
289	fix_approx_hermitian!(A; test_approx_symmetry=true, kwargs...)
290	Replaces the matrix `A` by `½(A + A')` in place. This will be successful and the result
291	is guaranteed to pass the `ishermitian` test only if the matrix is square and already
292	approximately hermitian.
293
294	By default logs an error message if the matrix `A` is found to not be approximately
295	hermitian. Set `test_approx_symmetry=false` to bypass testing.
296	Other keyword arguments are passed on to `isapprox`.
297	"""
298	function fix_approx_hermitian!(A; test_approx_symmetry=true, kwargs...)	4✔
299	# Generic and inefficient version. Make sure to replace by efficient specialised code.
300	if test_approx_symmetry	2✔
301	passed = isapprox(A, A'; kwargs...)	1✔
302	if !passed	1✔
303	throw(ArgumentError("Matrix is not approximately hermitian."))	1✔
304	end
305	end
306	@. A = 1 / 2 * (A + A')	2✔
307	return A	1✔
308	end
309
310	# specialised code for sparse matrices
311	using SparseArrays: AbstractSparseMatrixCSC, getcolptr, rowvals, nonzeros
312
313	# special case for sparse matrices; avoids most allocations, testing is free
314	function fix_approx_hermitian!(A::AbstractSparseMatrixCSC; test_approx_symmetry=false, kwargs...)	46✔
315	passed = isapprox_enforce_hermitian!(A; kwargs...)	23✔
316	if test_approx_symmetry && !passed	23✔
317	throw(ArgumentError("Matrix is not approximately hermitian."))	1✔
318	end
319	return A	22✔
320	end
321
322	"""
323	isapprox_enforce_hermitian!(A::AbstractSparseMatrixCSC; kwargs...) -> Bool
324	Checks whether the matrix `A` is approximately hermitian by checking each pair of transposed
325	matrix elements with `isapprox`. Keyword arguments are passed on to `isapprox`.
326	Returns boolean `true` is the test is passed and `false` if not.
327
328	Furthermore, the matrix `A` is modified to become exactly equal to `½(A + A')` if the test
329	is passed.
330	"""
331	function isapprox_enforce_hermitian!(A::AbstractSparseMatrixCSC; kwargs...)	46✔
332	# based on `ishermsym()` from `SparseArrays`; relies on `SparseArrays` internals
333	# https://github.com/JuliaSparse/SparseArrays.jl/blob/1bae96dc8f9a8ca8b7879eef4cf71e186598e982/src/sparsematrix.jl#L3793
334	m, n = size(A)	23✔
335	if m != n	23✔
UNCOV 336	return false	×
337	end
338
339	colptr = getcolptr(A)	23✔
340	rowval = rowvals(A)	23✔
341	nzval = nonzeros(A)	23✔
342	tracker = copy(getcolptr(A))	23✔
343	for col = 1:size(A, 2)	23✔
344	# `tracker` is updated such that, for symmetric matrices,
345	# the loop below starts from an element at or below the
346	# diagonal element of column `col`"
347	for p = tracker[col]:colptr[col+1]-1	154✔
348	val = nzval[p]	494✔
349	row = rowval[p]	494✔
350
351	# Ignore stored zeros
352	if iszero(val)	494✔
353	continue	4✔
354	end
355
356	# If the matrix was symmetric we should have updated
357	# the tracker to start at the diagonal or below. Here
358	# we are above the diagonal so the matrix can't be symmetric.
359	if row < col	490✔
UNCOV 360	return false	×
361	end
362
363	# Diagonal element
364	if row == col	490✔
365	if isapprox(val, conj(val); kwargs...)	150✔
366	nzval[p] = real(val)	150✔
367	else
UNCOV 368	return false	×
369	end
370	# if val != conj(val)
371	# return false
372	# end
373	else
374	offset = tracker[row]	340✔
375
376	# If the matrix is unsymmetric, there might not exist
377	# a rowval[offset]
378	if offset > length(rowval)	340✔
UNCOV 379	return false	×
380	end
381
382	row2 = rowval[offset]	340✔
383
384	# row2 can be less than col if the tracker didn't
385	# get updated due to stored zeros in previous elements.
386	# We therefore "catch up" here while making sure that
387	# the elements are actually zero.
388	while row2 < col	340✔
UNCOV 389	if !iszero(nzval[offset])	×
UNCOV 390	return false	×
391	end
UNCOV 392	offset += 1	×
UNCOV 393	row2 = rowval[offset]	×
UNCOV 394	tracker[row] += 1	×
UNCOV 395	end	×
396
397	# Non zero A[i,j] exists but A[j,i] does not exist
398	if row2 > col	340✔
UNCOV 399	return false	×
400	end
401
402	# A[i,j] and A[j,i] exists
403	if row2 == col	340✔
404	if isapprox(val, conj(nzval[offset]); kwargs...)	407✔
405	val = 1 / 2 * (val + conj(nzval[offset]))	339✔
406	nzval[p] = val	339✔
407	nzval[offset] = conj(val)	339✔
408	else
409	return false	1✔
410	end
411	# if val != conj(nzval[offset])
412	# return false
413	# end
414	tracker[row] += 1	339✔
415	end
416	end
417	end	833✔
418	end	284✔
419	return true	22✔
420	end
421
422	function Base.show(io::IO, b::BasisSetRepresentation)	2✔
423	print(io, "BasisSetRepresentation($(b.hamiltonian)) with dimension $(dimension(b))" *	2✔
424	" and $(nnz(b.sparse_matrix)) stored entries:")
425	show(io, MIME"text/plain"(), b.sparse_matrix)	2✔
426	end
427
428	Interfaces.starting_address(bsr::BasisSetRepresentation) = bsr.basis[1]	1✔
429
430	Hamiltonians.dimension(bsr::BasisSetRepresentation) = length(bsr.basis)	24✔
431
432	"""
433	sparse(hamiltonian::AbstractHamiltonian, addr=starting_address(hamiltonian); kwargs...)
434	sparse(bsr::BasisSetRepresentation)
435
436	Return a sparse matrix representation of `hamiltonian` or `bsr`. `kwargs` are passed to
437	[`BasisSetRepresentation`](@ref).
438
439	See [`BasisSetRepresentation`](@ref).
440	"""
441	function SparseArrays.sparse(hamiltonian::AbstractHamiltonian, args...; kwargs...)	48✔
442	return sparse(BasisSetRepresentation(hamiltonian, args...; kwargs...))	24✔
443	end
444	SparseArrays.sparse(bsr::BasisSetRepresentation) = bsr.sparse_matrix	26✔
445
446	"""
447	Matrix(
448	hamiltonian::AbstractHamiltonian, addr=starting_address(hamiltonian);
449	sizelim=10^4, kwargs...
450	)
451	Matrix(bsr::BasisSetRepresentation)
452
453	Return a dense matrix representation of `hamiltonian` or `bsr`. `kwargs` are passed to
454	[`BasisSetRepresentation`](@ref).
455
456	See [`BasisSetRepresentation`](@ref).
457	"""
458	function Base.Matrix(hamiltonian::AbstractHamiltonian, args...; sizelim=1e4, kwargs...)	126✔
459	return Matrix(BasisSetRepresentation(hamiltonian, args...; sizelim, kwargs...))	63✔
460	end
461	Base.Matrix(bsr::BasisSetRepresentation) = Matrix(bsr.sparse_matrix)	139✔

RimuQMC / Rimu.jl / 12369488912

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous