16921074123

Committed 12 Aug 2025 09:15PM UTC coverage: 93.978% (-0.2%) from 94.128%

Build # 16921074123

Build Type

push

github

Committed by

web-flow

Commit Message

New interface implementation for `HubbardRealSpace` (#334)

# Changes
- `HubbardRealSpace` now uses the new interface which is more efficient.
- `find_mode` now accepts an optional `occupied_mode_map`. Passing it
increases performance significantly for `BoseFS`.
- `each_mode` iterator for iterating through all modes in an address.

Run Details

136 of 144 new or added lines in 6 files covered. (94.44%)

23 existing lines in 2 files now uncovered.

7163 of 7622 relevant lines covered (93.98%)

15399400.14 hits per line

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

96.71

/src/BitStringAddresses/bosefs.jl

"""
    BoseFS{N,M,S} <: SingleComponentFockAddress

Address type that represents a Fock state of `N` spinless bosons in `M` modes by wrapping a
[`BitString`](@ref), or a [`SortedParticleList`](@ref). Which is wrapped is chosen
automatically based on the properties of the address.

# Constructors

* `BoseFS{[N,M]}(val::Integer...)`: Create `BoseFS{N,M}` from occupation numbers. This is
  type-stable if the number of modes `M` and the number of particles `N` are provided.
  Otherwise, `M` and `N` are inferred from the arguments.

* `BoseFS{[N,M]}(onr)`: Create `BoseFS{N,M}` from occupation number representation, see
  [`onr`](@ref). This is efficient if `N` and `M` are provided, and `onr` is a
  statically-sized collection, such as a `Tuple` or `SVector`.

* `BoseFS{[N,M]}([M, ]pairs...)`: Provide the number of modes `M` and `mode =>
  occupation_number` pairs. If `M` is provided as a type parameter, it should not be
  provided as the first argument.  Useful for creating sparse addresses. `pairs` can be
  multiple arguments or an iterator of pairs.

* `BoseFS{N,M,S}(bs::S)`: Unsafe constructor. Does not check whether the number of
  particles in `bs` is equal to `N`.

* [`@fs_str`](@ref): Addresses are sometimes printed in a compact manner. This
  representation can also be used as a constructor. See the examples below.

# Examples

```jldoctest
julia> BoseFS{6,5}(0, 1, 2, 3, 0)
BoseFS{6,5}(0, 1, 2, 3, 0)

julia> BoseFS(abs(i - 3) ≤ 1 ? i - 1 : 0 for i in 1:5)
BoseFS{6,5}(0, 1, 2, 3, 0)

julia> BoseFS(5, 2 => 1, 3 => 2, 4 => 3)
BoseFS{6,5}(0, 1, 2, 3, 0)

julia> BoseFS{6,5}(i => i - 1 for i in 2:4)
BoseFS{6,5}(0, 1, 2, 3, 0)

julia> fs"|0 1 2 3 0⟩"
BoseFS{6,5}(0, 1, 2, 3, 0)

julia> fs"|b 5: 2 3 3 4 4 4⟩"
BoseFS{6,5}(0, 1, 2, 3, 0)
```

See also: [`SingleComponentFockAddress`](@ref), [`OccupationNumberFS`](@ref),
[`FermiFS`](@ref), [`CompositeFS`](@ref), [`FermiFS2C`](@ref).
"""
struct BoseFS{N,M,S} <: SingleComponentFockAddress{N,M}
    bs::S
end

@inline function BoseFS{N,M,S}(onr::Union{SVector{M},MVector{M},NTuple{M}}) where {N,M,S}
    @boundscheck begin
        sum(onr) == N || throw(ArgumentError(
            "invalid ONR: $N particles expected, $(sum(onr)) given"
        ))
        if S <: BitString
            B = num_bits(S)
            M + N - 1 == B || throw(ArgumentError(
                "invalid ONR: $B-bit BitString does not fit $N particles in $M modes"
            ))
        elseif S <: SortedParticleList
            N == num_particles(S) && M == num_modes(S) || throw(ArgumentError(
                "invalid ONR: $S does not fit $N particles in $M modes"
            ))
        end
    end
    return BoseFS{N,M,S}(from_bose_onr(S, onr))
end
function BoseFS{N,M}(onr::Union{AbstractArray{<:Integer},NTuple{M,<:Integer}}) where {N,M}
    @boundscheck begin
        sum(onr) == N || throw(ArgumentError(
            "invalid ONR: $N particles expected, $(sum(onr)) given"
        ))
        length(onr) == M || throw(ArgumentError(
            "invalid ONR: $M modes expected, $(length(onr)) given"
        ))
    end
    spl_type = select_int_type(M)

    # Pick smaller address type, but prefer sparse.
    # Alway pick dense if it fits into one chunk.

    # Compute the size of container in words
    sparse_sizeof = ceil(Int, N * sizeof(spl_type) / 8)
    dense_sizeof = ceil(Int, (N + M - 1) / 64)
    if dense_sizeof == 1 || dense_sizeof < sparse_sizeof
        S = typeof(BitString{M + N - 1}(0))
    else
        S = SortedParticleList{N,M,spl_type}
    end
    return BoseFS{N,M,S}(from_bose_onr(S, onr))
end
function BoseFS(onr) # single argument constructor
    onr = Tuple(onr)
    M = length(onr)
    N = sum(onr)
    return BoseFS{N,M}(onr)
end
BoseFS(vals::Integer...) = BoseFS(vals) # specify occupation numbers
BoseFS(val::Integer) = BoseFS((val,)) # single mode address
BoseFS{N,M}(vals::Integer...) where {N,M} = BoseFS{N,M}(vals)

# Sparse constructors
BoseFS(M::Integer, pairs::Pair...) = BoseFS(M, pairs)
BoseFS(M::Integer, pairs) = BoseFS(sparse_to_onr(M, pairs))
BoseFS{N,M}(pairs::Pair...) where {N,M} = BoseFS{N,M}(pairs)
BoseFS{N,M}(pairs) where {N,M} = BoseFS{N,M}(sparse_to_onr(M, pairs))
BoseFS(pairs::Pair...) = throw(ArgumentError("number of modes must be provided"))


function print_address(io::IO, b::BoseFS{N,M}; compact=false) where {N,M}
    if compact && b.bs isa SortedParticleList
        print(io, "|b ", M, ": ", join(Int.(b.bs.storage), ' '), "⟩")
    elseif compact
        print(io, "|", join(onr(b), ' '), "⟩")
    elseif b.bs isa SortedParticleList
        print(io, "BoseFS{$N,$M}(", onr_sparse_string(onr(b)), ")")
    else
        print(io, "BoseFS{$N,$M}", tuple(onr(b)...))
    end
end

Base.bitstring(b::BoseFS) = bitstring(b.bs) # TODO rename?

Base.isless(a::BoseFS, b::BoseFS) = isless(a.bs, b.bs)
Base.hash(bba::BoseFS,  h::UInt) = hash(bba.bs, h)
Base.:(==)(a::BoseFS, b::BoseFS) = a.bs == b.bs

"""
    near_uniform_onr(N, M) -> onr::SVector{M,Int}

Create occupation number representation `onr` distributing `N` particles in `M`
modes in a close-to-uniform fashion with each mode filled with at least
`N ÷ M` particles and at most with `N ÷ M + 1` particles.
"""
function near_uniform_onr(n::Number, m::Number)
    return near_uniform_onr(Val(n),Val(m))
end
function near_uniform_onr(::Val{N}, ::Val{M}) where {N, M}
    fillingfactor, extras = divrem(N, M)
    # startonr = fill(fillingfactor,M)
    startonr = fillingfactor * @MVector ones(Int,M)
    startonr[1:extras] .+= 1
    return SVector{M}(startonr)
end

"""
    near_uniform(BoseFS{N,M}) -> BoseFS{N,M}

Create bosonic Fock state with near uniform occupation number of `M` modes with
a total of `N` particles.

# Examples
```jldoctest
julia> near_uniform(BoseFS{7,5})
BoseFS{7,5}(2, 2, 1, 1, 1)

julia> near_uniform(FermiFS{3,5})
FermiFS{3,5}(1, 1, 1, 0, 0)
```
"""
function near_uniform(::Type{<:BoseFS{N,M}}) where {N,M}
    return BoseFS{N,M}(near_uniform_onr(Val(N),Val(M)))
end
near_uniform(b::AbstractFockAddress) = near_uniform(typeof(b))

onr(b::BoseFS{<:Any,M}) where {M} = to_bose_onr(b.bs, Val(M))
const occupation_number_representation = onr # resides here because `onr` has to be defined

function Base.reverse(b::BoseFS)
    return typeof(b)(reverse(b.bs))
end

# For vacuum state
function num_occupied_modes(b::BoseFS{0})
    return 0
end
function num_occupied_modes(b::BoseFS)
    return bose_num_occupied_modes(b.bs)
end
function occupied_modes(b::BoseFS{N,M,S}) where {N,M,S}
    return BoseOccupiedModes{N,M,S}(b.bs)
end

function find_mode(b::BoseFS, index, occ=occupied_modes(b))
    last_occnum = last_mode = last_offset = 0
    for (occnum, mode, offset) in occ
        dist = index - mode
        if dist == 0
            return BoseFSIndex(occnum, index, offset)
        elseif dist < 0
            return BoseFSIndex(0, index, offset + dist)
        end
        last_occnum = occnum
        last_mode = mode
        last_offset = offset
    end
    offset = last_offset + last_occnum + index - last_mode
    return BoseFSIndex(0, index, offset)
end
# Multiple in a single pass
function find_mode(b::BoseFS, indices::NTuple{N}, occ=occupied_modes(b)) where {N}
    # Idea: find permutation, then use the permutation to find indices in order even though
    # they are not sorted.
    perm = sortperm(SVector(indices))
    # perm_i is the index in permutation and goes from 1:N.
    perm_i = 1
    # curr_i points to indices and result
    curr_i = perm[1]
    # index is the current index we are looking for.
    index = indices[curr_i]

    result = ntuple(_ -> BoseFSIndex(0, 0, 0), Val(N))
    last_occnum = last_mode = last_offset = 0
    @inbounds for (occnum, mode, offset) in occ
        dist = index - mode
        # While loop handles duplicate entries in indices.
        while dist ≤ 0
            if dist == 0
                @set! result[curr_i] = BoseFSIndex(occnum, mode, offset)
            else
                @set! result[curr_i] = BoseFSIndex(0, index, offset + dist)
            end
            perm_i += 1
            perm_i > N && return result
            curr_i = perm[perm_i]
            index = indices[curr_i]
            dist = index - mode
        end
        last_occnum = occnum
        last_mode = mode
        last_offset = offset
    end
    # Now we have to find all indices that appear after the last occupied site.
    # While true because we break out of the loop early anyway.
    @inbounds while true
        offset = last_offset + last_occnum + index - last_mode
        @set! result[curr_i] = BoseFSIndex(0, index, offset)
        perm_i += 1
        perm_i > N && return result
        curr_i = perm[perm_i]
        index = indices[curr_i]
    end
    return result # not reached
end

# Specialised version of each_mode for iterating modes in BoseFS with BitString storage.
# This is necessary because find_mode in a BitString-backed BoseFS is inefficient.
struct BoseBitStringEachMode{M,A<:BoseFS{<:Any,M,<:BitString}}
    address::A
end
Base.eltype(::BoseBitStringEachMode) = BoseFSIndex
Base.length(::BoseBitStringEachMode{M}) where {M} = M

function Base.iterate(em::BoseBitStringEachMode{M}, state=(0, 1, em.address.bs)) where {M}
    offset, mode, bitstring = state
    if mode > M
        return nothing
    else
        bosons = Int32(trailing_ones(bitstring))
        bitstring >>>= (bosons + 1) % UInt

        return BoseFSIndex(bosons, mode, offset), (offset + bosons + 1, mode + 1, bitstring)
    end
end

function each_mode(addr::BoseFS{<:Any,<:Any,<:BitString})
    return BoseBitStringEachMode(addr)
end

# find_occupied_mode provided by generic implementation

function excitation(b::B, creations, destructions) where {B<:BoseFS}
    new_bs, val = bose_excitation(b.bs, creations, destructions)
    return B(new_bs), val
end

"""
    new_address, value = hopnextneighbour(add, chosen, boundary_condition)

Compute the new address of a hopping event for the Hubbard model. Returns the new
address and the square root of product of occupation numbers of the involved modes
multiplied by a term consistent with boundary condition as the `value`.
The following boundary conditions are supported:

* `:periodic`: hopping over the boundary gives does not change the `value`.
* `:twisted`: hopping over the boundary flips the sign of the `value`.
* `:hard_wall`: hopping over the boundary gives a `value` of zero.
* `θ <: Number`: hopping over the boundary gives a `value` multiplied by ``\\exp(iθ)`` or ``\\exp(−iθ)`` depending on the direction of hopping.

The off-diagonals are indexed as follows:

* `(chosen + 1) ÷ 2` selects the hopping site.
* Even `chosen` indicates a hop to the left.
* Odd `chosen` indicates a hop to the right.

# Example

```jldoctest
julia> using Rimu.Hamiltonians: hopnextneighbour

julia> hopnextneighbour(BoseFS(1, 0, 1), 3)
(BoseFS{2,3}(2, 0, 0), 1.4142135623730951)

julia> hopnextneighbour(BoseFS(1, 0, 1), 4)
(BoseFS{2,3}(1, 1, 0), 1.0)

julia> hopnextneighbour(BoseFS(1, 0, 1), 3, :twisted)
(BoseFS{2,3}(2, 0, 0), -1.4142135623730951)

julia> hopnextneighbour(BoseFS(1, 0, 1), 3, :hard_wall)
(BoseFS{2,3}(2, 0, 0), 0.0)

julia> hopnextneighbour(BoseFS(1, 0, 1), 3, π/4)
(BoseFS{2,3}(2, 0, 0), 1.0000000000000002 + 1.0im)
```
"""
function hopnextneighbour(b::BoseFS{N,M,A}, chosen) where {N,M,A<:BitString}
    address = b.bs
    T = chunk_type(address)
    site = (chosen + 1) >>> 0x1
    if isodd(chosen) # Hopping to the right
        next = 0
        curr = 0
        offset = 0
        sc = 0
        reached_end = false
        for (i, (num, sn, bit)) in enumerate(occupied_modes(b))
            next = num * (sn == sc + 1) # only set next to > 0 if sites are neighbours
            reached_end = i == site + 1
            reached_end && break
            curr = num
            offset = bit + num
            sc = sn
        end
        if sc == M
            new_address = (address << 0x1) | A(T(1))
            prod = curr * (trailing_ones(address) + 1) # mul occupation num of first obital
        else
            next *= reached_end
            new_address = address ⊻ A(T(3)) << ((offset - 1) % T)
            prod = curr * (next + 1)
        end
    else # Hopping to the left
        if site == 1 && isodd(address)
            # For leftmost site, we shift the whole address circularly by one bit.
            new_address = (address >>> 0x1) | A(T(1)) << ((N + M - 2) % T)
            prod = trailing_ones(address) * leading_ones(new_address)
        else
            prev = 0
            curr = 0
            offset = 0
            sp = 0
            for (i, (num, sc, bit)) in enumerate(occupied_modes(b))
                prev = curr * (sc == sp + 1) # only set prev to > 0 if sites are neighbours
                curr = num
                offset = bit
                i == site && break
                sp = sc
            end
            new_address = address ⊻ A(T(3)) << ((offset - 1) % T)
            prod = curr * (prev + 1)
        end
    end
    return BoseFS{N,M,A}(new_address), √prod
end

function hopnextneighbour(b::SingleComponentFockAddress, i)
    src = find_occupied_mode(b, (i + 1) >>> 0x1)
    dst = find_mode(b, mod1(src.mode + ifelse(isodd(i), 1, -1), num_modes(b)))

    new_b, val = excitation(b, (dst,), (src,))
    return new_b, val
end

function hopnextneighbour(
    b::SingleComponentFockAddress, i, boundary_condition::Symbol)
    src = find_occupied_mode(b, (i + 1) >>> 0x1)
    dir = ifelse(isodd(i), 1, -1)
    dst = find_mode(b, mod1(src.mode + dir, num_modes(b)))
    new_b, val = excitation(b, (dst,), (src,))
    on_boundary = src.mode == 1 && dir == -1 || src.mode == num_modes(b) && dir == 1
    if boundary_condition == :twisted && on_boundary
        return new_b, -val
    elseif boundary_condition == :hard_wall && on_boundary
        return new_b, 0.0
    else
        return new_b, val
    end
end

function hopnextneighbour(b::SingleComponentFockAddress, i, boundary_condition::Number)
    src = find_occupied_mode(b, (i + 1) >>> 0x1)
    dir = ifelse(isodd(i), 1, -1)
    dst = find_mode(b, mod1(src.mode + dir, num_modes(b)))
    new_b, val = excitation(b, (dst,), (src,))
    if (src.mode == 1 && dir == -1)
        return new_b, val*exp(-im*boundary_condition)
    elseif (src.mode == num_modes(b) && dir == 1)
        return new_b, val*exp(im*boundary_condition)
    else
        return new_b, complex(val)
    end
end

"""
    bose_hubbard_interaction(address)

Return ``Σ_i n_i (n_i-1)`` for computing the Bose-Hubbard on-site interaction (without the
``U`` prefactor.)

# Example

```jldoctest
julia> Hamiltonians.bose_hubbard_interaction(BoseFS{4,4}((2,1,1,0)))
2
julia> Hamiltonians.bose_hubbard_interaction(BoseFS{4,4}((3,0,1,0)))
6
```
"""
function bose_hubbard_interaction(b::BoseFS{<:Any,<:Any,A}) where {A<:BitString}
    return bose_hubbard_interaction(Val(num_chunks(A)), b)
end
function bose_hubbard_interaction(b::SingleComponentFockAddress)
    return bose_hubbard_interaction(nothing, b)
end

@inline function bose_hubbard_interaction(_, b::SingleComponentFockAddress)
    result = 0
    for (n, _, _) in occupied_modes(b)
        result += n * (n - 1)
    end
    return result
end

@inline function bose_hubbard_interaction(::Val{1}, b::BoseFS{<:Any,<:Any,<:BitString})
    # currently this ammounts to counting occupation numbers of modes
    chunk = chunks(b.bs)[1]
    matrixelementint = 0
    while !iszero(chunk)
        chunk >>>= (trailing_zeros(chunk) % UInt) # proceed to next occupied mode
        bosonnumber = trailing_ones(chunk) # count how many bosons inside
        # surpsingly it is faster to not check whether this is nonzero and do the
        # following operations anyway
        chunk >>>= (bosonnumber % UInt) # remove the counted mode
        matrixelementint += bosonnumber * (bosonnumber - 1)
    end
    return matrixelementint
end

1	"""
2	BoseFS{N,M,S} <: SingleComponentFockAddress
3
4	Address type that represents a Fock state of `N` spinless bosons in `M` modes by wrapping a
5	[`BitString`](@ref), or a [`SortedParticleList`](@ref). Which is wrapped is chosen
6	automatically based on the properties of the address.
7
8	# Constructors
9
10	* `BoseFS{[N,M]}(val::Integer...)`: Create `BoseFS{N,M}` from occupation numbers. This is
11	type-stable if the number of modes `M` and the number of particles `N` are provided.
12	Otherwise, `M` and `N` are inferred from the arguments.
13
14	* `BoseFS{[N,M]}(onr)`: Create `BoseFS{N,M}` from occupation number representation, see
15	[`onr`](@ref). This is efficient if `N` and `M` are provided, and `onr` is a
16	statically-sized collection, such as a `Tuple` or `SVector`.
17
18	* `BoseFS{[N,M]}([M, ]pairs...)`: Provide the number of modes `M` and `mode =>
19	occupation_number` pairs. If `M` is provided as a type parameter, it should not be
20	provided as the first argument. Useful for creating sparse addresses. `pairs` can be
21	multiple arguments or an iterator of pairs.
22
23	* `BoseFS{N,M,S}(bs::S)`: Unsafe constructor. Does not check whether the number of
24	particles in `bs` is equal to `N`.
25
26	* [`@fs_str`](@ref): Addresses are sometimes printed in a compact manner. This
27	representation can also be used as a constructor. See the examples below.
28
29	# Examples
30
31	```jldoctest
32	julia> BoseFS{6,5}(0, 1, 2, 3, 0)
33	BoseFS{6,5}(0, 1, 2, 3, 0)
34
35	julia> BoseFS(abs(i - 3) ≤ 1 ? i - 1 : 0 for i in 1:5)
36	BoseFS{6,5}(0, 1, 2, 3, 0)
37
38	julia> BoseFS(5, 2 => 1, 3 => 2, 4 => 3)
39	BoseFS{6,5}(0, 1, 2, 3, 0)
40
41	julia> BoseFS{6,5}(i => i - 1 for i in 2:4)
42	BoseFS{6,5}(0, 1, 2, 3, 0)
43
44	julia> fs"\|0 1 2 3 0⟩"
45	BoseFS{6,5}(0, 1, 2, 3, 0)
46
47	julia> fs"\|b 5: 2 3 3 4 4 4⟩"
48	BoseFS{6,5}(0, 1, 2, 3, 0)
49	```
50
51	See also: [`SingleComponentFockAddress`](@ref), [`OccupationNumberFS`](@ref),
52	[`FermiFS`](@ref), [`CompositeFS`](@ref), [`FermiFS2C`](@ref).
53	"""
54	struct BoseFS{N,M,S} <: SingleComponentFockAddress{N,M}
55	bs::S	463,159,354✔
56	end
57
58	@inline function BoseFS{N,M,S}(onr::Union{SVector{M},MVector{M},NTuple{M}}) where {N,M,S}	11,523✔
59	@boundscheck begin	11,523✔
60	sum(onr) == N \|\| throw(ArgumentError(	11,523✔
61	"invalid ONR: $N particles expected, $(sum(onr)) given"
62	))
63	if S <: BitString	11,523✔
64	B = num_bits(S)	11,523✔
65	M + N - 1 == B \|\| throw(ArgumentError(	11,523✔
66	"invalid ONR: $B-bit BitString does not fit $N particles in $M modes"
67	))
68	elseif S <: SortedParticleList	×
69	N == num_particles(S) && M == num_modes(S) \|\| throw(ArgumentError(	×
70	"invalid ONR: $S does not fit $N particles in $M modes"
71	))
72	end
73	end
74	return BoseFS{N,M,S}(from_bose_onr(S, onr))	91,791✔
75	end
76	function BoseFS{N,M}(onr::Union{AbstractArray{<:Integer},NTuple{M,<:Integer}}) where {N,M}	244,160✔
77	@boundscheck begin	244,160✔
78	sum(onr) == N \|\| throw(ArgumentError(	244,161✔
79	"invalid ONR: $N particles expected, $(sum(onr)) given"
80	))
81	length(onr) == M \|\| throw(ArgumentError(	244,160✔
82	"invalid ONR: $M modes expected, $(length(onr)) given"
83	))
84	end
85	spl_type = select_int_type(M)	244,158✔
86
87	# Pick smaller address type, but prefer sparse.
88	# Alway pick dense if it fits into one chunk.
89
90	# Compute the size of container in words
91	sparse_sizeof = ceil(Int, N * sizeof(spl_type) / 8)	244,158✔
92	dense_sizeof = ceil(Int, (N + M - 1) / 64)	244,158✔
93	if dense_sizeof == 1 \|\| dense_sizeof < sparse_sizeof	244,158✔
94	S = typeof(BitString{M + N - 1}(0))	239,389✔
95	else
96	S = SortedParticleList{N,M,spl_type}	4,769✔
97	end
98	return BoseFS{N,M,S}(from_bose_onr(S, onr))	244,930✔
99	end
100	function BoseFS(onr) # single argument constructor	1,747✔
101	onr = Tuple(onr)	1,747✔
102	M = length(onr)	1,747✔
103	N = sum(onr)	5,925✔
104	return BoseFS{N,M}(onr)	2,494✔
105	end
106	BoseFS(vals::Integer...) = BoseFS(vals) # specify occupation numbers	260✔
107	BoseFS(val::Integer) = BoseFS((val,)) # single mode address	3✔
108	BoseFS{N,M}(vals::Integer...) where {N,M} = BoseFS{N,M}(vals)	71✔
109
110	# Sparse constructors
111	BoseFS(M::Integer, pairs::Pair...) = BoseFS(M, pairs)	36✔
112	BoseFS(M::Integer, pairs) = BoseFS(sparse_to_onr(M, pairs))	564✔
113	BoseFS{N,M}(pairs::Pair...) where {N,M} = BoseFS{N,M}(pairs)	11✔
114	BoseFS{N,M}(pairs) where {N,M} = BoseFS{N,M}(sparse_to_onr(M, pairs))	12✔
115	BoseFS(pairs::Pair...) = throw(ArgumentError("number of modes must be provided"))	1✔
116
117
118	function print_address(io::IO, b::BoseFS{N,M}; compact=false) where {N,M}	1,578✔
119	if compact && b.bs isa SortedParticleList	789✔
120	print(io, "\|b ", M, ": ", join(Int.(b.bs.storage), ' '), "⟩")	10✔
121	elseif compact	779✔
122	print(io, "\|", join(onr(b), ' '), "⟩")	500✔
123	elseif b.bs isa SortedParticleList	279✔
124	print(io, "BoseFS{$N,$M}(", onr_sparse_string(onr(b)), ")")	40✔
125	else
126	print(io, "BoseFS{$N,$M}", tuple(onr(b)...))	239✔
127	end
128	end
129
130	Base.bitstring(b::BoseFS) = bitstring(b.bs) # TODO rename?	×
131
132	Base.isless(a::BoseFS, b::BoseFS) = isless(a.bs, b.bs)	232,861✔
133	Base.hash(bba::BoseFS, h::UInt) = hash(bba.bs, h)	855,564,186✔
134	Base.:(==)(a::BoseFS, b::BoseFS) = a.bs == b.bs	34,246,794✔
135
136	"""
137	near_uniform_onr(N, M) -> onr::SVector{M,Int}
138
139	Create occupation number representation `onr` distributing `N` particles in `M`
140	modes in a close-to-uniform fashion with each mode filled with at least
141	`N ÷ M` particles and at most with `N ÷ M + 1` particles.
142	"""
143	function near_uniform_onr(n::Number, m::Number)	×
144	return near_uniform_onr(Val(n),Val(m))	×
145	end
146	function near_uniform_onr(::Val{N}, ::Val{M}) where {N, M}	28✔
147	fillingfactor, extras = divrem(N, M)	28✔
148	# startonr = fill(fillingfactor,M)
149	startonr = fillingfactor * @MVector ones(Int,M)	28✔
150	startonr[1:extras] .+= 1	28✔
151	return SVector{M}(startonr)	28✔
152	end
153
154	"""
155	near_uniform(BoseFS{N,M}) -> BoseFS{N,M}
156
157	Create bosonic Fock state with near uniform occupation number of `M` modes with
158	a total of `N` particles.
159
160	# Examples
161	```jldoctest
162	julia> near_uniform(BoseFS{7,5})
163	BoseFS{7,5}(2, 2, 1, 1, 1)
164
165	julia> near_uniform(FermiFS{3,5})
166	FermiFS{3,5}(1, 1, 1, 0, 0)
167	```
168	"""
169	function near_uniform(::Type{<:BoseFS{N,M}}) where {N,M}	28✔
170	return BoseFS{N,M}(near_uniform_onr(Val(N),Val(M)))	28✔
171	end
172	near_uniform(b::AbstractFockAddress) = near_uniform(typeof(b))	1✔
173
174	onr(b::BoseFS{<:Any,M}) where {M} = to_bose_onr(b.bs, Val(M))	29,121,708✔
175	const occupation_number_representation = onr # resides here because `onr` has to be defined
176
177	function Base.reverse(b::BoseFS)	3,082✔
178	return typeof(b)(reverse(b.bs))	3,082✔
179	end
180
181	# For vacuum state
182	function num_occupied_modes(b::BoseFS{0})	1✔
183	return 0	1✔
184	end
185	function num_occupied_modes(b::BoseFS)	284,643,146✔
186	return bose_num_occupied_modes(b.bs)	553,365,820✔
187	end
188	function occupied_modes(b::BoseFS{N,M,S}) where {N,M,S}	803,019,087✔
189	return BoseOccupiedModes{N,M,S}(b.bs)	803,182,134✔
190	end
191
192	function find_mode(b::BoseFS, index, occ=occupied_modes(b))	442,292,496✔
193	last_occnum = last_mode = last_offset = 0	732,475,749✔
194	for (occnum, mode, offset) in occ	441,712,665✔
195	dist = index - mode	417,490,060✔
196	if dist == 0	417,388,812✔
197	return BoseFSIndex(occnum, index, offset)	76,675,373✔
198	elseif dist < 0	340,946,366✔
199	return BoseFSIndex(0, index, offset + dist)	85,823,646✔
200	end
201	last_occnum = occnum	255,424,194✔
202	last_mode = mode	255,506,229✔
203	last_offset = offset	255,498,315✔
204	end	451,377,764✔
205	offset = last_offset + last_occnum + index - last_mode	59,408,493✔
206	return BoseFSIndex(0, index, offset)	59,408,012✔
207	end
208	# Multiple in a single pass
209	function find_mode(b::BoseFS, indices::NTuple{N}, occ=occupied_modes(b)) where {N}	480,359,117✔
210	# Idea: find permutation, then use the permutation to find indices in order even though
211	# they are not sorted.
212	perm = sortperm(SVector(indices))	480,316,169✔
213	# perm_i is the index in permutation and goes from 1:N.
214	perm_i = 1	243,300,654✔
215	# curr_i points to indices and result
216	curr_i = perm[1]	243,287,236✔
217	# index is the current index we are looking for.
218	index = indices[curr_i]	243,333,954✔
219
220	result = ntuple(_ -> BoseFSIndex(0, 0, 0), Val(N))	726,176,232✔
221	last_occnum = last_mode = last_offset = 0	243,881,000✔
222	@inbounds for (occnum, mode, offset) in occ	486,363,970✔
223	dist = index - mode	897,258,658✔
224	# While loop handles duplicate entries in indices.
225	while dist ≤ 0	1,112,704,851✔
226	if dist == 0	425,173,987✔
227	@set! result[curr_i] = BoseFSIndex(occnum, mode, offset)	217,475,060✔
228	else
229	@set! result[curr_i] = BoseFSIndex(0, index, offset + dist)	211,835,667✔
230	end
231	perm_i += 1	425,110,477✔
232	perm_i > N && return result	424,882,522✔
233	curr_i = perm[perm_i]	235,324,312✔
234	index = indices[curr_i]	235,394,679✔
235	dist = index - mode	235,662,144✔
236	end	235,257,102✔
237	last_occnum = occnum	711,176,299✔
238	last_mode = mode	710,765,635✔
239	last_offset = offset	713,190,587✔
240	end	323,042,366✔
241	# Now we have to find all indices that appear after the last occupied site.
242	# While true because we break out of the loop early anyway.
243	@inbounds while true	50,379,474✔
244	offset = last_offset + last_occnum + index - last_mode	59,533,302✔
245	@set! result[curr_i] = BoseFSIndex(0, index, offset)	59,535,633✔
246	perm_i += 1	59,537,010✔
247	perm_i > N && return result	59,541,102✔
248	curr_i = perm[perm_i]	9,168,396✔
249	index = indices[curr_i]	9,168,455✔
250	end	9,168,557✔
251	return result # not reached	×
252	end
253
254	# Specialised version of each_mode for iterating modes in BoseFS with BitString storage.
255	# This is necessary because find_mode in a BitString-backed BoseFS is inefficient.
256	struct BoseBitStringEachMode{M,A<:BoseFS{<:Any,M,<:BitString}}
257	address::A	100✔
258	end
NEW 259	Base.eltype(::BoseBitStringEachMode) = BoseFSIndex	×
260	Base.length(::BoseBitStringEachMode{M}) where {M} = M	100✔
261
262	function Base.iterate(em::BoseBitStringEachMode{M}, state=(0, 1, em.address.bs)) where {M}	6,200✔
263	offset, mode, bitstring = state	6,200✔
264	if mode > M	6,100✔
265	return nothing	100✔
266	else
267	bosons = Int32(trailing_ones(bitstring))	6,000✔
268	bitstring >>>= (bosons + 1) % UInt	6,000✔
269
270	return BoseFSIndex(bosons, mode, offset), (offset + bosons + 1, mode + 1, bitstring)	6,000✔
271	end
272	end
273
274	function each_mode(addr::BoseFS{<:Any,<:Any,<:BitString})	100✔
275	return BoseBitStringEachMode(addr)	100✔
276	end
277
278	# find_occupied_mode provided by generic implementation
279
280	function excitation(b::B, creations, destructions) where {B<:BoseFS}	461,201,019✔
281	new_bs, val = bose_excitation(b.bs, creations, destructions)	677,899,496✔
282	return B(new_bs), val	462,799,896✔
283	end
284
285	"""
286	new_address, value = hopnextneighbour(add, chosen, boundary_condition)
287
288	Compute the new address of a hopping event for the Hubbard model. Returns the new
289	address and the square root of product of occupation numbers of the involved modes
290	multiplied by a term consistent with boundary condition as the `value`.
291	The following boundary conditions are supported:
292
293	* `:periodic`: hopping over the boundary gives does not change the `value`.
294	* `:twisted`: hopping over the boundary flips the sign of the `value`.
295	* `:hard_wall`: hopping over the boundary gives a `value` of zero.
296	* `θ <: Number`: hopping over the boundary gives a `value` multiplied by ``\\exp(iθ)`` or ``\\exp(−iθ)`` depending on the direction of hopping.
297
298	The off-diagonals are indexed as follows:
299
300	* `(chosen + 1) ÷ 2` selects the hopping site.
301	* Even `chosen` indicates a hop to the left.
302	* Odd `chosen` indicates a hop to the right.
303
304	# Example
305
306	```jldoctest
307	julia> using Rimu.Hamiltonians: hopnextneighbour
308
309	julia> hopnextneighbour(BoseFS(1, 0, 1), 3)
310	(BoseFS{2,3}(2, 0, 0), 1.4142135623730951)
311
312	julia> hopnextneighbour(BoseFS(1, 0, 1), 4)
313	(BoseFS{2,3}(1, 1, 0), 1.0)
314
315	julia> hopnextneighbour(BoseFS(1, 0, 1), 3, :twisted)
316	(BoseFS{2,3}(2, 0, 0), -1.4142135623730951)
317
318	julia> hopnextneighbour(BoseFS(1, 0, 1), 3, :hard_wall)
319	(BoseFS{2,3}(2, 0, 0), 0.0)
320
321	julia> hopnextneighbour(BoseFS(1, 0, 1), 3, π/4)
322	(BoseFS{2,3}(2, 0, 0), 1.0000000000000002 + 1.0im)
323	```
324	"""
325	function hopnextneighbour(b::BoseFS{N,M,A}, chosen) where {N,M,A<:BitString}	3,752✔
326	address = b.bs	3,752✔
327	T = chunk_type(address)	3,752✔
328	site = (chosen + 1) >>> 0x1	3,752✔
329	if isodd(chosen) # Hopping to the right	3,752✔
330	next = 0	1,876✔
331	curr = 0	1,876✔
332	offset = 0	1,876✔
333	sc = 0	1,876✔
334	reached_end = false	1,876✔
335	for (i, (num, sn, bit)) in enumerate(occupied_modes(b))	2,009✔
336	next = num * (sn == sc + 1) # only set next to > 0 if sites are neighbours	89,001✔
337	reached_end = i == site + 1	89,001✔
338	reached_end && break	89,001✔
339	curr = num	87,176✔
340	offset = bit + num	87,176✔
341	sc = sn	87,176✔
342	end	87,176✔
343	if sc == M	1,876✔
344	new_address = (address << 0x1) \| A(T(1))	35✔
345	prod = curr * (trailing_ones(address) + 1) # mul occupation num of first obital	35✔
346	else
347	next *= reached_end	1,841✔
348	new_address = address ⊻ A(T(3)) << ((offset - 1) % T)	1,844✔
349	prod = curr * (next + 1)	1,841✔
350	end
351	else # Hopping to the left
352	if site == 1 && isodd(address)	1,876✔
353	# For leftmost site, we shift the whole address circularly by one bit.
354	new_address = (address >>> 0x1) \| A(T(1)) << ((N + M - 2) % T)	35✔
355	prod = trailing_ones(address) * leading_ones(new_address)	35✔
356	else
357	prev = 0	1,841✔
358	curr = 0	1,841✔
359	offset = 0	1,841✔
360	sp = 0	1,841✔
361	for (i, (num, sc, bit)) in enumerate(occupied_modes(b))	1,967✔
362	prev = curr * (sc == sp + 1) # only set prev to > 0 if sites are neighbours	87,141✔
363	curr = num	87,141✔
364	offset = bit	87,141✔
365	i == site && break	87,141✔
366	sp = sc	85,300✔
367	end	85,300✔
368	new_address = address ⊻ A(T(3)) << ((offset - 1) % T)	1,841✔
369	prod = curr * (prev + 1)	1,841✔
370	end
371	end
372	return BoseFS{N,M,A}(new_address), √prod	3,752✔
373	end
374
375	function hopnextneighbour(b::SingleComponentFockAddress, i)	386✔
376	src = find_occupied_mode(b, (i + 1) >>> 0x1)	386✔
377	dst = find_mode(b, mod1(src.mode + ifelse(isodd(i), 1, -1), num_modes(b)))	386✔
378
379	new_b, val = excitation(b, (dst,), (src,))	772✔
380	return new_b, val	386✔
381	end
382
383	function hopnextneighbour(	215,977,184✔
384	b::SingleComponentFockAddress, i, boundary_condition::Symbol)
385	src = find_occupied_mode(b, (i + 1) >>> 0x1)	342,164,426✔
386	dir = ifelse(isodd(i), 1, -1)	216,065,292✔
387	dst = find_mode(b, mod1(src.mode + dir, num_modes(b)))	501,062,060✔
388	new_b, val = excitation(b, (dst,), (src,))	431,501,024✔
389	on_boundary = src.mode == 1 && dir == -1 \|\| src.mode == num_modes(b) && dir == 1	398,662,719✔
390	if boundary_condition == :twisted && on_boundary	216,067,497✔
391	return new_b, -val	13✔
392	elseif boundary_condition == :hard_wall && on_boundary	216,044,117✔
393	return new_b, 0.0	9✔
394	else
395	return new_b, val	216,037,023✔
396	end
397	end
398
399	function hopnextneighbour(b::SingleComponentFockAddress, i, boundary_condition::Number)	2,261✔
400	src = find_occupied_mode(b, (i + 1) >>> 0x1)	4,074✔
401	dir = ifelse(isodd(i), 1, -1)	2,261✔
402	dst = find_mode(b, mod1(src.mode + dir, num_modes(b)))	4,661✔
403	new_b, val = excitation(b, (dst,), (src,))	4,522✔
404	if (src.mode == 1 && dir == -1)	2,261✔
405	return new_b, valexp(-imboundary_condition)	192✔
406	elseif (src.mode == num_modes(b) && dir == 1)	2,069✔
407	return new_b, valexp(imboundary_condition)	191✔
408	else
409	return new_b, complex(val)	1,878✔
410	end
411	end
412
413	"""
414	bose_hubbard_interaction(address)
415
416	Return ``Σ_i n_i (n_i-1)`` for computing the Bose-Hubbard on-site interaction (without the
417	``U`` prefactor.)
418
419	# Example
420
421	```jldoctest
422	julia> Hamiltonians.bose_hubbard_interaction(BoseFS{4,4}((2,1,1,0)))
423	2
424	julia> Hamiltonians.bose_hubbard_interaction(BoseFS{4,4}((3,0,1,0)))
425	6
426	```
427	"""
428	function bose_hubbard_interaction(b::BoseFS{<:Any,<:Any,A}) where {A<:BitString}	123,165,993✔
429	return bose_hubbard_interaction(Val(num_chunks(A)), b)	123,138,534✔
430	end
431	function bose_hubbard_interaction(b::SingleComponentFockAddress)	50✔
432	return bose_hubbard_interaction(nothing, b)	66✔
433	end
434
435	@inline function bose_hubbard_interaction(_, b::SingleComponentFockAddress)	83✔
436	result = 0	83✔
437	for (n, _, _) in occupied_modes(b)	131✔
438	result += n * (n - 1)	2,034✔
439	end	4,009✔
440	return result	83✔
441	end
442
443	@inline function bose_hubbard_interaction(::Val{1}, b::BoseFS{<:Any,<:Any,<:BitString})	123,136,383✔
444	# currently this ammounts to counting occupation numbers of modes
445	chunk = chunks(b.bs)[1]	123,071,821✔
446	matrixelementint = 0	123,071,043✔
447	while !iszero(chunk)	442,406,073✔
448	chunk >>>= (trailing_zeros(chunk) % UInt) # proceed to next occupied mode	320,421,973✔
449	bosonnumber = trailing_ones(chunk) # count how many bosons inside	320,289,547✔
450	# surpsingly it is faster to not check whether this is nonzero and do the
451	# following operations anyway
452	chunk >>>= (bosonnumber % UInt) # remove the counted mode	320,588,839✔
453	matrixelementint += bosonnumber * (bosonnumber - 1)	320,653,568✔
454	end	320,835,653✔
455	return matrixelementint	122,950,538✔
456	end

RimuQMC / Rimu.jl / 16921074123

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