Fix bug in dimension

src/Hamiltonians/ExtendedHubbardReal1D.jl

@@ -33,6 +33,14 @@ function starting_address(h::ExtendedHubbardReal1D)
     return getfield(h, :add)
 end
 
+# `ExtendedHubbardReal1D` conserves particle number. Thus we can lower the bound on the dimension
+# for the non-conserving `OccupationNumberFS`.
+function dimension(::ExtendedHubbardReal1D, a::OccupationNumberFS)
+    m = num_modes(a)
+    n = num_particles(a)
+    return dimension(BoseFS{n,m})
+end
+
 LOStructure(::Type{<:ExtendedHubbardReal1D{<:Real}}) = IsHermitian()
 
 Base.getproperty(h::ExtendedHubbardReal1D, s::Symbol) = getproperty(h, Val(s))

src/Hamiltonians/HubbardMom1D.jl

@@ -70,6 +70,14 @@ function starting_address(h::HubbardMom1D)
     return h.add
 end
 
+# `HubbardMom1D` conserves particle number. Thus we can lower the bound on the dimension
+# for the non-conserving `OccupationNumberFS`.
+function dimension(::HubbardMom1D, a::OccupationNumberFS)
+    m = num_modes(a)
+    n = num_particles(a)
+    return dimension(BoseFS{n,m})
+end
+
 LOStructure(::Type{<:HubbardMom1D{<:Real}}) = IsHermitian()
 
 Base.getproperty(h::HubbardMom1D, s::Symbol) = getproperty(h, Val(s))

src/Hamiltonians/HubbardMom1DEP.jl

@@ -100,6 +100,15 @@ end
 function starting_address(h::HubbardMom1DEP)
     return h.add
 end
+
+# `HubbardMom1DEP` conserves particle number. Thus we can lower the bound on the dimension
+# for the non-conserving `OccupationNumberFS`.
+function dimension(::HubbardMom1DEP, a::OccupationNumberFS)
+    m = num_modes(a)
+    n = num_particles(a)
+    return dimension(BoseFS{n,m})
+end
+
 function get_offdiagonal(h::HubbardMom1DEP{<:Any,<:Any,F}, add::F, i) where {F}
     return offdiagonals(h, add)[i]
 end

src/Hamiltonians/HubbardReal1D.jl

@@ -36,6 +36,14 @@ function starting_address(h::HubbardReal1D)
     return getfield(h, :add)
 end
 
+# `HubbardReal1D` conserves particle number. Thus we can lower the bound on the dimension
+# for the non-conserving `OccupationNumberFS`.
+function dimension(::HubbardReal1D, a::OccupationNumberFS)
+    m = num_modes(a)
+    n = num_particles(a)
+    return dimension(BoseFS{n,m})
+end
+
 LOStructure(::Type{<:HubbardReal1D{<:Real}}) = IsHermitian()
 
 Base.getproperty(h::HubbardReal1D, s::Symbol) = getproperty(h, Val(s))

src/Hamiltonians/HubbardReal1DEP.jl

@@ -72,6 +72,14 @@ Base.getproperty(h::HubbardReal1DEP, ::Val{:ep}) = getfield(h, :ep)
 
 starting_address(h::HubbardReal1DEP) = h.add
 
+# `HubbardReal1DEP` conserves particle number. Thus we can lower the bound on the dimension
+# for the non-conserving `OccupationNumberFS`.
+function dimension(::HubbardReal1DEP, a::OccupationNumberFS)
+    m = num_modes(a)
+    n = num_particles(a)
+    return dimension(BoseFS{n,m})
+end
+
 function num_offdiagonals(::HubbardReal1DEP, address::SingleComponentFockAddress)
     return 2 * num_occupied_modes(address)
 end

src/Hamiltonians/HubbardRealSpace.jl

@@ -272,6 +272,16 @@ end
 Base.:(==)(H::HubbardRealSpace, G::HubbardRealSpace) = all(map(p -> getproperty(H, p) == getproperty(G, p), propertynames(H)))
 
 starting_address(h::HubbardRealSpace) = h.address
+
+# `HubbardRealSpace` conserves particle number. Thus we can lower the bound on the dimension
+# for the non-conserving `OccupationNumberFS`.
+# NOTE: We should think of a more general mechanism that works for composite addresses.
+function dimension(::HubbardRealSpace, a::OccupationNumberFS)
+    m = num_modes(a)
+    n = num_particles(a)
+    return dimension(BoseFS{n,m})
+end
+
 function diagonal_element(h::HubbardRealSpace, address)
     int = isnothing(h.u) ? 0.0 : local_interaction(address, h.u)
     pot = isnothing(h.v) ? 0.0 : external_potential(address, h.potential)

src/Hamiltonians/ParitySymmetry.jl

@@ -68,6 +68,8 @@ function starting_address(h::ParitySymmetry)
     return min(add, reverse(add))
 end
 
+dimension(h::ParitySymmetry, addr) = dimension(h.hamiltonian, addr) # upper bound
+
 get_offdiagonal(h::ParitySymmetry, add, i) = offdiagonals(h, add)[i]
 num_offdiagonals(h::ParitySymmetry, add) = num_offdiagonals(h.hamiltonian, add)
 

src/Hamiltonians/abstract.jl

@@ -22,12 +22,13 @@ BitStringAddresses.num_modes(h::AbstractHamiltonian) = num_modes(starting_addres
 """
     dimension(h::AbstractHamiltonian, addr=starting_address(h))
     dimension(addr::AbstractFockAddress)
+    dimension(::Type{<:AbstractFockAddress})
 
 Return the estimated dimension of Hilbert space. May return a `BigInt` number.
 
-When called on an address, the dimension of the linear space spanned by the address type is
-returned. When called on an `AbstractHamiltonian`, an upper bound on the dimension of
-the matrix representing the Hamiltonian is returned.
+When called on an address or address type, the dimension of the linear space spanned by the
+address type is returned. When called on an `AbstractHamiltonian`, an upper bound on the
+dimension of the matrix representing the Hamiltonian is returned.
 
 # Examples
 
@@ -44,32 +45,36 @@ julia> dimension(HubbardReal1D(near_uniform(BoseFS{200,100})))
 julia> dimension(HubbardReal1D(near_uniform(BoseFS{200,100})))|>Float64
 1.3860838210861882e81
 ```
-# Interface
-
-When extending `AbstractHamiltonian`, define a method for the two-argument form
-`dimension(h::MyNewHamiltonian, addr)`.
 
 See also [`BasisSetRep`](@ref).
+# Extended Help
+
+When extending `AbstractHamiltonian`, define a method for the two-argument form
+`dimension(h::MyNewHamiltonian, addr)`. When extending `AbstractFockAddress`, define a
+method for `dimension(::Type{MyNewFockAddress})`.
 """
 dimension(h::AbstractHamiltonian) = dimension(h, starting_address(h))
 dimension(::AbstractHamiltonian, addr) = dimension(addr)
+dimension(addr::AbstractFockAddress) = dimension(typeof(addr))
 # dimension(_) = Inf # fallback
 
-function dimension(::BoseFS{N,M}) where {N,M}
+function dimension(::Type{<:BoseFS{N,M}}) where {N,M}
     return binomial(BigInt(N + M - 1), BigInt(N))
 end
-function dimension(::OccupationNumberFS{M,T}) where {M,T}
+function dimension(::Type{<:OccupationNumberFS{M,T}}) where {M,T}
     n = typemax(T)
-    return binomial(BigInt(n + M - 1), BigInt(n))
+    return BigInt(n + 1)^BigInt(M)
 end
-function dimension(::FermiFS{N,M}) where {N,M}
+function dimension(::Type{<:FermiFS{N,M}}) where {N,M}
     return binomial(BigInt(M), BigInt(N))
 end
-function dimension(b::BoseFS2C)
-    return dimension(b.bsa) * dimension(b.bsb)
+function dimension(::Type{<:BoseFS2C{NA,NB,M}}) where {NA,NB,M}
+    return dimension(BoseFS{NA,M}) * dimension(BoseFS{NB,M})
 end
-function dimension(c::CompositeFS)
-    return prod(x -> dimension(x), c.components)
+function dimension(::Type{<:CompositeFS{<:Any,<:Any,<:Any,T}}) where {T}
+    return prod(dimension, T.parameters)
+    # This relies on an implementation detail of the Tuple type and may break in future
+    # julia versions.
 end
 
 # for backward compatibility

test/BitStringAddresses.jl

@@ -567,6 +567,7 @@ end
         @test sparse(ExtendedHubbardReal1D(ofs)) == sparse(ExtendedHubbardReal1D(bfs))
         oham = HubbardReal1D(OccupationNumberFS(0, 2, 1))
         bham = HubbardReal1D(BoseFS(0, 2, 1))
-        @test sparse(ParitySymmetry(oham; odd=true)) == sparse(ParitySymmetry(bham; odd=true))
+        @test sparse(ParitySymmetry(oham; odd=true)) ==
+            sparse(ParitySymmetry(bham; odd=true))
     end
 end

test/Hamiltonians.jl

@@ -60,8 +60,6 @@ function test_hamiltonian_interface(H)
         end
         @testset "dimension" begin
             @test dimension(H) ≥ dimension(H, starting_address(H))
-            @test dimension(Float64, H) isa Float64
-            @test dimension(Int, H) == dimension(H)
         end
         @testset "allowed_address_type" begin
             @test addr isa allowed_address_type(H)
@@ -880,7 +878,7 @@ end
     m = 6
     n1 = 4
     n2 = m
-    
+
     # unital refers to n̄=1
     non_unital_localised_state = BoseFS((n1,0,0,0,0,0))
     non_unital_uniform_state = near_uniform(non_unital_localised_state)
@@ -893,7 +891,7 @@ end
     S2 = StringCorrelator(2)
 
     @test num_offdiagonals(S0, localised_state) == 0
-    
+
     # non unital localised state
     @test @inferred diagonal_element(S0, non_unital_localised_state) ≈ 20/9
     @test @inferred diagonal_element(S1, non_unital_localised_state) ≈ (-4/9)*exp(im * -2pi/3)
@@ -918,7 +916,7 @@ end
     d = 5
     output = @capture_out print(StringCorrelator(d))
     @test output == "StringCorrelator($d)"
-    
+
 end
 
 @testset "Momentum" begin
@@ -1588,3 +1586,10 @@ end
         @test isempty(fock_to_cart(null_addr, S))
     end
 end
+
+@testset "dimension and multi-component addresses" begin
+    addresses = [CompositeFS(FermiFS((1,0,1)), FermiFS((0,1,0))), BoseFS((1,0,1)),
+                FermiFS2C((1,0,1), (0,1,0)), BoseFS2C((1,0,1), (0,1,0))
+    ]
+    [@test dimension(addr) == dimension(typeof(addr)) for addr in addresses]
+end