Implement SingleComponentFockAddress interface for OccupationNumberFS

src/BitStringAddresses/fockaddress.jl

@@ -65,10 +65,10 @@ Find the `i`-th mode in address. Returns [`BoseFSIndex`](@ref) for [`BoseFS`](@r
 bounds.
 
 ```jldoctest
-julia> find_mode(BoseFS((1, 0, 2)), 2)
+julia> find_mode(BoseFS(1, 0, 2), 2)
 BoseFSIndex(occnum=0, mode=2, offset=2)
 
-julia> find_mode(FermiFS((1, 1, 1, 0)), (2,3))
+julia> find_mode(FermiFS(1, 1, 1, 0), (2,3))
 (FermiFSIndex(occnum=1, mode=2, offset=1), FermiFSIndex(occnum=1, mode=3, offset=2))
 ```
 
@@ -87,13 +87,13 @@ Returns [`BoseFSIndex`](@ref) for [`BoseFS`](@ref), and [`FermiFSIndex`](@ref) f
 # Example
 
 ```jldoctest
-julia> find_occupied_mode(FermiFS((1, 1, 1, 0)), 2)
+julia> find_occupied_mode(FermiFS(1, 1, 1, 0), 2)
 FermiFSIndex(occnum=1, mode=2, offset=1)
 
-julia> find_occupied_mode(BoseFS((1, 0, 2)), 1)
+julia> find_occupied_mode(BoseFS(1, 0, 2), 1)
 BoseFSIndex(occnum=1, mode=1, offset=0)
 
-julia> find_occupied_mode(BoseFS((1, 0, 2)), 1, 2)
+julia> find_occupied_mode(BoseFS(1, 0, 2), 1, 2)
 BoseFSIndex(occnum=2, mode=3, offset=3)
 ```
 
@@ -242,7 +242,8 @@ function OccupiedModeMap(addr::SingleComponentFockAddress{N,M}) where {N,M}
     modes = occupied_modes(addr)
     T = eltype(modes)
     # There are at most N occupied modes. This could be also @generated for cases where N ≫ M
-    indices = MVector{min(N,M),T}(undef)
+    L = ismissing(N) ? M : min(N, M)
+    indices = MVector{L,T}(undef)
     i = 0
     for index in modes
         i += 1
@@ -489,7 +490,7 @@ Struct used for indexing and performing [`excitation`](@ref)s on a [`BoseFS`](@r
  represented by `SortedParticleList`.
 
 """
-struct BoseFSIndex<:FieldVector{3,Int}
+Base.@kwdef struct BoseFSIndex<:FieldVector{3,Int}
     occnum::Int
     mode::Int
     offset::Int
@@ -604,7 +605,7 @@ Struct used for indexing and performing [`excitation`](@ref)s on a [`FermiFS`](@
   represented by a bitstring, and the position in the list when using `SortedParticleList`.
 
 """
-struct FermiFSIndex<:FieldVector{3,Int}
+Base.@kwdef struct FermiFSIndex<:FieldVector{3,Int}
     occnum::Int
     mode::Int
     offset::Int

src/BitStringAddresses/occupationnumberfs.jl

@@ -64,7 +64,6 @@ function print_address(io::IO, ofs::OccupationNumberFS{M,T}; compact=false) wher
 end
 
 onr(ofs::OccupationNumberFS) = ofs.onr
-num_occupied_modes(ofs::OccupationNumberFS) = count(!iszero, onr(ofs))
 num_particles(ofs::OccupationNumberFS) = Int(sum(onr(ofs)))
 # `num_modes` does not have to be defined here, because it is defined for the abstract type
 
@@ -113,12 +112,13 @@ See also: [`destroy`](@ref), [`excitation`](@ref).
 end
 
 """
-    excitation(addr::OccupationNumberFS, c::NTuple{<:Any,Int}, d::NTuple{<:Any,Int})
+    excitation(addr::OccupationNumberFS, c::NTuple, d::NTuple)
     → (nadd, α)
 Generate an excitation on an [`OccupationNumberFS`](@ref) by applying the creation and
 destruction operators specified by the tuples of mode numbers `c` and `d` to the Fock state
-`addr`. The modes are simply indexed by integers (starting at 1). The value of `α` is given
-by the square root of the product of mode occupations before destruction and after creation.
+`addr`. The modes are indexed by integers (starting at 1), or by indices of type
+`BoseFSIndex`. The value of `α` is given by the square root of the product of mode
+occupations before destruction and after creation.
 
 The number of particles may change by this type of excitation.
 
@@ -137,7 +137,11 @@ julia> num_particles(es)
 7
 ```
 """
-function excitation(fs::OccupationNumberFS{<:Any,T}, c::NTuple{<:Any,Int}, d::NTuple{<:Any,Int}) where {T}
+function excitation(
+    fs::OccupationNumberFS{<:Any,T},
+    c::NTuple{<:Any,Int},
+    d::NTuple{<:Any,Int}
+) where {T}
     accu = one(T)
     for i in d
         fs, val = destroy(fs, i)
@@ -149,6 +153,40 @@ function excitation(fs::OccupationNumberFS{<:Any,T}, c::NTuple{<:Any,Int}, d::NT
     end
     return fs, √accu
 end
+function excitation(
+    fs::OccupationNumberFS,
+    c::NTuple{N1,BoseFSIndex},
+    d::NTuple{N2,BoseFSIndex}
+) where {N1, N2}
+    creations = ntuple(i -> c[i].mode, Val(N1)) # convert BoseFSIndex to mode number
+    destructions = ntuple(i -> d[i].mode, Val(N2))
+    return excitation(fs, creations, destructions)
+end
+
+
+# `SingleComponentFockAddress` interface
+
+find_mode(ofs::OccupationNumberFS, n::Integer) = BoseFSIndex(ofs.onr[n], n, n)
+function find_mode(ofs::OccupationNumberFS, ns::NTuple{N,Integer}) where N
+    return ntuple(i -> find_mode(ofs, ns[i]), Val(N))
+end
+
+num_occupied_modes(ofs::OccupationNumberFS) = count(!iszero, onr(ofs))
+
+# for the lazy iterator `occupied_modes` we adapt the `BoseOccupiedModes` type
+function occupied_modes(ofs::OccupationNumberFS{M}) where {M}
+    return BoseOccupiedModes{missing, M, typeof(ofs)}(ofs)
+end
 
-# Do we need more methods for building Hamiltonians? (`find_occupied_mode`,
-# `OccupiedModMap`, `occupied_modes`?)
+function Base.length(bom::BoseOccupiedModes{<:Any,<:Any,<:OccupationNumberFS})
+    return num_occupied_modes(bom.storage)
+end
+
+function Base.iterate(bom::BoseOccupiedModes{<:Any,<:Any,<:OccupationNumberFS}, i=1)
+    s = onr(bom.storage) # is an SVector with the onr
+    while true
+        i > length(s) && return nothing
+        iszero(s[i]) || return BoseFSIndex(s[i], i, i), i + 1
+        i += 1
+    end
+end

test/BitStringAddresses.jl

@@ -525,6 +525,11 @@ end
         fs_after_excitation, sqrt_accu = excitation(ofs, c, d)
         @test fs_after_excitation.onr == SVector{3,UInt8}(2, 1, 3)
         @test sqrt_accu ≈ √4
+
+        # indexing with BoseFSIndex
+        i, j = find_mode(ofs, (1, 2))
+        @test j ==  BoseFSIndex(occnum=2, mode=2, offset=2)
+        @test excitation(ofs, (i,), (j,)) == (fs_after_excitation, sqrt_accu)
     end
 
     @testset "Test properties of OccupationNumberFS" begin
@@ -534,5 +539,7 @@ end
         @test onr(ofs) == SVector{3,UInt8}(1, 2, 3)
         lfs = OccupationNumberFS{6}([1 0 0; 1 1 0])
         @test onr(lfs, LadderBoundaries(2, 3)) == [1 0 0; 1 1 0]
+        @test num_occupied_modes(lfs) == length(occupied_modes(lfs)) == 3
+        @test OccupiedModeMap(lfs) == collect(occupied_modes(lfs))
     end
 end