Observables (#278)

* Allow Vector's of AbstractHamiltonian's in AllOverlaps * allows_address_type replaces allowed_address_type in the AbstractHamiltonians interface * valtype for AbstractHamiltonian and ProjectorMonteCarloProblem error message if using complex Hamiltonian with ProjectorMonteCarloProblem * fix docs * attempt to fix docs * AbstractOperator Define a new supertype for AbstractHamiltonian with less stringent requirements on the interface * use VectorInterface.scalartype for AbstractOperators (instead of Base.valtype) * move dot_from_right to Interfaces.jl as part of the AbstractOperator Interface * fix tests * make correlators AbstractOperator * fix docs * docs * fix error in test * test_operator_interface IsDeterministic default for complex PDVec * fix test * test operator interface for G2RealSpace empty constructor for AbstractDVec accepts a style argument * fix test * changes to call structure of StringCorrelator test more AbstractOperators * fix docs * demote ParticleNumberOperator to AbstractHamiltonian; it now will not accept passing an address * make SingleParticleExcitation and TwoParticleExcitation AbstractOperator * Momentum is AbstractOperator * Hamiltonian interface test for momentum
RimuQMC · Aug 19, 2024 · 7cb56f7 · 7cb56f7 · mtsch · Sep 13, 2024
1 parent ecda02c
commit 7cb56f7
Show file tree

Hide file tree

Showing 35 changed files with 796 additions and 437 deletions.
diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "Rimu"
 uuid = "c53c40cc-bd84-11e9-2cf4-a9fde2b9386e"
 authors = ["Joachim Brand <j.brand@massey.ac.nz>"]
-version = "0.12.1"
+version = "0.13.0"
 
 [deps]
 Arrow = "69666777-d1a9-59fb-9406-91d4454c9d45"

diff --git a/docs/src/hamiltonians.md b/docs/src/hamiltonians.md
@@ -69,10 +69,16 @@ Stoquastic
 ```
 
 ## Observables
-Observables are [`AbstractHamiltonian`](@ref)s that represent a physical
-observable. Their ground state expectation values can be sampled by passing
-them into [`AllOverlaps`](@ref).
+Observables are [`AbstractOperator`](@ref)s that represent a physical
+observable. Their expectation values can be sampled during a
+[`ProjectorMonteCarloProblem`](@ref) simulation by passing
+them into a suitable [`ReplicaStrategy`](@ref), e.g. 
+[`AllOverlaps`](@ref).  [`AbstractOperator`](@ref) is a supertype of 
+[`AbstractHamiltonian`](@ref) and has less stringent 
+requirements. Some observables are also [`AbstractHamiltonian`](@ref)s.
+
 ```@docs
+AbstractOperator
 ParticleNumberOperator
 G2RealCorrelator
 G2RealSpace
@@ -114,7 +120,10 @@ random_offdiagonal
 Hamiltonians.LOStructure
 dimension
 has_adjoint
-allowed_address_type
+allows_address_type
+Base.eltype
+VectorInterface.scalartype
+mul!
 ```
 
 This interface relies on unexported functionality, including

diff --git a/scripts/G2-example.jl b/scripts/G2-example.jl
@@ -47,8 +47,8 @@ G2list = ((G2RealCorrelator(d) for d in dvals)...,)
 # To obtain an unbiased result, at least two replicas should be used. One can also use more
 # than two to improve the statistics. This is particularly helpful when the walker number is
 # low.
-num_replicas = 3
-replica_strategy = AllOverlaps(num_replicas; operator = G2list)
+n_replicas = 3
+replica_strategy = AllOverlaps(n_replicas; operator=G2list)
 
 # Other FCIQMC parameters and strategies can be set in the same way as before.
 steps_equilibrate = 1_000
@@ -93,7 +93,7 @@ println(filter(contains("Op"), names(df)))
 # equilibration steps.
 
 # Now, we can calculate the correlation function for each value of ``d``.
-println("Two-body correlator from $num_replicas replicas:")
+println("Two-body correlator from $n_replicas replicas:")
 for d in dvals
     r = rayleigh_replica_estimator(df; op_name = "Op$(d+1)", skip=steps_equilibrate)
     println("   G2($d) = $(r.f) ± $(r.σ_f)")
@@ -105,8 +105,8 @@ end
 
 # Since we ran multiple independent replicas, we also have multiple estimates of the shift
 # energy.
-println("Shift energy from $num_replicas replicas:")
-for i in 1:num_replicas
+println("Shift energy from $n_replicas replicas:")
+for i in 1:n_replicas
     se = shift_estimator(df; shift="shift_$i", skip=steps_equilibrate)
     println("   Replica $i: $(se.mean) ± $(se.err)")
 end

diff --git a/src/BitStringAddresses/multicomponent.jl b/src/BitStringAddresses/multicomponent.jl
@@ -72,7 +72,7 @@ function CompositeFS(adds::Vararg{SingleComponentFockAddress})
     return CompositeFS{length(adds),N,M1,typeof(adds)}(adds)
 end
 
-num_components(::CompositeFS{C}) where {C} = C
+num_components(::Type{<:CompositeFS{C}}) where {C} = C
 Base.hash(c::CompositeFS, u::UInt) = hash(c.components, u)
 
 function print_address(io::IO, c::CompositeFS{C}; compact=false) where {C}

diff --git a/src/DictVectors/DictVectors.jl b/src/DictVectors/DictVectors.jl
@@ -21,7 +21,8 @@ import MPI
 using ..Interfaces: Interfaces, AbstractDVec, AdjointUnknown,
     CompressionStrategy, IsDiagonal, LOStructure,
     apply_column!, apply_operator!, compress!,
-    diagonal_element, offdiagonals, step_stats, AbstractHamiltonian
+    diagonal_element, offdiagonals, step_stats, AbstractHamiltonian, AbstractOperator,
+    dot_from_right
 using ..StochasticStyles: StochasticStyles, IsDeterministic
 
 import ..Interfaces: deposit!, storage, StochasticStyle, default_style, freeze, localpart,

diff --git a/src/DictVectors/abstractdvec.jl b/src/DictVectors/abstractdvec.jl
@@ -263,7 +263,7 @@ walkernumber_and_length(v) = walkernumber(v), length(v)
 ###
 ### Vector-operator functions
 ###
-function LinearAlgebra.mul!(w::AbstractDVec, h::AbstractHamiltonian, v::AbstractDVec)
+function LinearAlgebra.mul!(w::AbstractDVec, h::AbstractOperator, v::AbstractDVec)
     empty!(w)
     for (key, val) in pairs(v)
         w[key] += diagonal_element(h, key) * val
@@ -274,36 +274,34 @@ function LinearAlgebra.mul!(w::AbstractDVec, h::AbstractHamiltonian, v::Abstract
     return w
 end
 
-function Base.:*(h::AbstractHamiltonian, v::AbstractDVec)
-    return mul!(zerovector(v, promote_type(eltype(h), valtype(v))), h, v)
+function Base.:*(h::AbstractOperator, v::AbstractDVec)
+    T = promote_type(scalartype(h), scalartype(v))
+    if eltype(h) ≠ scalartype(h)
+        throw(ArgumentError("Operators with non-scalar eltype don't support "*
+                            "multiplication with `*`. Use `mul!` or `dot` instead."))
+    end
+    # first argument in mul! requires IsDeterministic style
+    w = empty(v, T; style=IsDeterministic{T}())
+    return mul!(w, h, v)
 end
 
-"""
-    dot(w, op::AbstractHamiltonian, v)
-
-Evaluate `w⋅op(v)` minimizing memory allocations.
-"""
-function LinearAlgebra.dot(w::AbstractDVec, op::AbstractHamiltonian, v::AbstractDVec)
+# docstring in Interfaces
+function LinearAlgebra.dot(w::AbstractDVec, op::AbstractOperator, v::AbstractDVec)
     return dot(LOStructure(op), w, op, v)
 end
-function LinearAlgebra.dot(::AdjointUnknown, w, op::AbstractHamiltonian, v)
+function LinearAlgebra.dot(::AdjointUnknown, w, op::AbstractOperator, v)
     return dot_from_right(w, op, v)
 end
-function LinearAlgebra.dot(::LOStructure, w, op::AbstractHamiltonian, v)
+function LinearAlgebra.dot(::LOStructure, w, op::AbstractOperator, v)
     if length(w) < length(v)
         return conj(dot_from_right(v, op', w)) # turn args around to execute faster
     else
         return dot_from_right(w, op, v) # original order
     end
 end
 
-"""
-    dot_from_right(w, op::AbstractHamiltonian, v)
-
-Internal function evaluates the 3-argument `dot()` function in order from right
-to left.
-"""
-function dot_from_right(w, op, v::AbstractDVec)
+# docstring in Interfaces
+function Interfaces.dot_from_right(w, op, v::AbstractDVec)
     T = typeof(zero(valtype(w)) * zero(eltype(op)) * zero(valtype(v)))
     result = zero(T)
     for (key, val) in pairs(v)

diff --git a/src/DictVectors/dvec.jl b/src/DictVectors/dvec.jl
@@ -41,7 +41,7 @@ DVec{Symbol,Float64} with 2 entries, style = IsDeterministic{Float64}()
   :b => 3.0
 ```
 """
-struct DVec{K,V,S<:StochasticStyle{V},D<:AbstractDict{K,V}} <: AbstractDVec{K,V}
+struct DVec{K,V,S<:StochasticStyle,D<:AbstractDict{K,V}} <: AbstractDVec{K,V}
     storage::D
     style::S
 end
@@ -50,16 +50,40 @@ end
 ### Constructors
 ###
 # Vararg
-function DVec(args...; kwargs...)
+function DVec(args...; style = nothing, kwargs...)
     storage = Dict(args...)
-    return DVec(storage; kwargs...)
+    if style === nothing
+        return DVec(storage; kwargs...)
+    else
+        # if style is given, check value type compatibility and convert if necessary
+        V = valtype(storage)
+        S = scalartype(style)
+        if V <: AbstractArray
+            if eltype(V) != S
+                throw(ArgumentError("Style $(style) is incompatible with value type $(V). "*
+                            "Use a style with scalartype(style) == $(eltype(V)) instead."))
+            end
+            # no conversion necessary, the types are compatible
+        else # V is a scalar type
+            # try to convert storage to the appropriate type
+            K = keytype(storage)
+            storage = Dict{K,S}(storage)
+        end
+        return DVec(storage; style, kwargs...)
+    end
 end
 
+# from Dict; check style compatibility
 function DVec(
-    dict::Dict{K}; style::StochasticStyle{V}=default_style(valtype(dict)), capacity=0
+    dict::Dict{K,V};
+    style::StochasticStyle=default_style(V),
+    capacity=0
 ) where {K,V}
-    storage = convert(Dict{K,V}, dict)
-    return DVec(storage, style)
+    if !(V == scalartype(style) || eltype(V) == scalartype(style))
+        throw(ArgumentError("Style $(style) is incompatible with value type $(V). "*
+                            "Use a style with scalartype(style) == $(eltype(V)) instead."))
+    end
+    return DVec{K,V,typeof(style),typeof(dict)}(dict, style)
 end
 # Empty constructor.
 function DVec{K,V}(; style::StochasticStyle=default_style(V), capacity=0) where {K,V}
@@ -71,17 +95,17 @@ function DVec(dv::AbstractDVec{K,V}; style=StochasticStyle(dv), capacity=0) wher
     return copyto!(dvec, dv)
 end
 
-function Base.empty(dvec::DVec{K,V}) where {K,V}
-    return DVec{K,V}(; style=StochasticStyle(dvec))
+function Base.empty(dvec::DVec{K,V}; style=dvec.style) where {K,V}
+    return DVec{K,V}(; style)
 end
-function Base.empty(dvec::DVec{K,V}, ::Type{V}) where {K,V}
-    return empty(dvec)
+function Base.empty(dvec::DVec{K,V}, ::Type{V}; style=dvec.style) where {K,V}
+    return empty(dvec; style)
 end
-function Base.empty(dvec::DVec{K,V}, ::Type{W}) where {K,V,W}
-    return DVec{K,W}()
+function Base.empty(::DVec{K}, ::Type{V}; style=default_style(V)) where {K,V}
+    return DVec{K,V}(; style)
 end
-function Base.empty(dvec::DVec, ::Type{K}, ::Type{V}) where {K,V}
-    return DVec{K,V}()
+function Base.empty(::DVec, ::Type{K}, ::Type{V}; style=default_style(V)) where {K,V}
+    return DVec{K,V}(; style)
 end
 
 ###

diff --git a/src/DictVectors/initiatordvec.jl b/src/DictVectors/initiatordvec.jl
@@ -120,17 +120,17 @@ function InitiatorDVec(
     return InitiatorDVec(copy(storage(dv)); style, initiator, capacity)
 end
 
-function Base.empty(dvec::InitiatorDVec{K,V}) where {K,V}
-    return InitiatorDVec{K,V}(; style=dvec.style, initiator=dvec.initiator)
+function Base.empty(dvec::InitiatorDVec{K,V}; style=dvec.style) where {K,V}
+    return InitiatorDVec{K,V}(; style, initiator=dvec.initiator)
 end
-function Base.empty(dvec::InitiatorDVec{K,V}, ::Type{V}) where {K,V}
-    return empty(dvec)
+function Base.empty(dvec::InitiatorDVec{K,V}, ::Type{V}; style=dvec.style) where {K,V}
+    return empty(dvec; style)
 end
-function Base.empty(dvec::InitiatorDVec{K}, ::Type{V}) where {K,V}
-    return InitiatorDVec{K,V}(; initiator=dvec.initiator)
+function Base.empty(dvec::InitiatorDVec{K}, ::Type{V}; style=default_style(V)) where {K,V}
+    return InitiatorDVec{K,V}(; style, initiator=dvec.initiator)
 end
-function Base.empty(dvec::InitiatorDVec, ::Type{K}, ::Type{V}) where {K,V}
-    return InitiatorDVec{K,V}(; initiator=dvec.initiator)
+function Base.empty(dvec::InitiatorDVec, ::Type{K}, ::Type{V}; style=default_style(V)) where {K,V}
+    return InitiatorDVec{K,V}(; style, initiator=dvec.initiator)
 end
 
 ###

diff --git a/src/DictVectors/pdvec.jl b/src/DictVectors/pdvec.jl
@@ -377,11 +377,15 @@ end
 function Base.empty(t::PDVec{K,V}, ::Type{V}; kwargs...) where {K,V}
     return empty(t; kwargs...)
 end
-function Base.empty(t::PDVec{K,<:Any,N}, ::Type{V}; kwargs...) where {K,V,N}
-    return PDVec{K,V,N}(; kwargs...)
+function Base.empty(
+    t::PDVec{K,<:Any,N}, ::Type{V}; style=default_style(V), kwargs...
+) where {K,V,N}
+    return PDVec{K,V,N}(; style, kwargs...)
 end
-function Base.empty(t::PDVec{<:Any,<:Any,N}, ::Type{K}, ::Type{V}; kwargs...) where {K,V,N}
-    return PDVec{K,V,N}(; kwargs...)
+function Base.empty(
+    t::PDVec{<:Any,<:Any,N}, ::Type{K}, ::Type{V}; style=default_style(V), kwargs...
+) where {K,V,N}
+    return PDVec{K,V,N}(; style, kwargs...)
 end
 function Base.empty!(t::PDVec)
     Folds.foreach(empty!, t.segments, )
@@ -744,53 +748,53 @@ end
 ### Operator linear algebra operations
 ###
 """
-    mul!(y::PDVec, A::AbstractHamiltonian, x::PDVec[, w::PDWorkingMemory])
+    mul!(y::PDVec, A::AbstractOperator, x::PDVec[, w::PDWorkingMemory])
 
 Perform `y = A * x` in-place. The working memory `w` is required to facilitate
 threaded/distributed operations. If not passed a new instance will be allocated. `y` and `x`
 may be the same vector.
 
-See [`PDVec`](@ref), [`PDWorkingMemory`](@ref).
+See [`PDVec`](@ref), [`PDWorkingMemory`](@ref), [`AbstractOperator`](@ref).
 """
 function LinearAlgebra.mul!(
-    y::PDVec, op::AbstractHamiltonian, x::PDVec,
-    w=PDWorkingMemory(y; style=IsDeterministic()),
+    y::PDVec, op::AbstractOperator, x::PDVec,
+    w=PDWorkingMemory(y; style=StochasticStyle(y)),
 )
-    if w.style ≢ IsDeterministic()
+    if !(w.style isa IsDeterministic)
         throw(ArgumentError(
             "Attempted to use `mul!` with non-deterministic working memory. " *
             "Use `apply_operator!` instead."
         ))
     end
-    _, _, wm, y = apply_operator!(w, y, x, op)
+    _, _, _, y = apply_operator!(w, y, x, op)
     return y
 end
 
 """
-    dot(y::PDVec, A::AbstractHamiltonian, x::PDVec[, w::PDWorkingMemory])
+    dot(y::PDVec, A::AbstractOperator, x::PDVec[, w::PDWorkingMemory])
 
 Perform `y ⋅ A ⋅ x`. The working memory `w` is required to facilitate threaded/distributed
 operations with non-diagonal `A`. If needed and not passed a new instance will be
 allocated. `A` can be replaced with a tuple of operators.
 
 See [`PDVec`](@ref), [`PDWorkingMemory`](@ref).
 """
-function LinearAlgebra.dot(t::PDVec, op::AbstractHamiltonian, u::PDVec, w)
+function LinearAlgebra.dot(t::PDVec, op::AbstractOperator, u::PDVec, w)
     return dot(LOStructure(op), t, op, u, w)
 end
-function LinearAlgebra.dot(t::PDVec, op::AbstractHamiltonian, u::PDVec)
+function LinearAlgebra.dot(t::PDVec, op::AbstractOperator, u::PDVec)
     return dot(LOStructure(op), t, op, u)
 end
 function LinearAlgebra.dot(
-    ::IsDiagonal, t::PDVec, op::AbstractHamiltonian, u::PDVec, w=nothing
+    ::IsDiagonal, t::PDVec, op::AbstractOperator, u::PDVec, w=nothing
 )
     T = typeof(zero(valtype(t)) * zero(valtype(u)) * zero(eltype(op)))
     return sum(pairs(u); init=zero(T)) do (k, v)
         T(conj(t[k]) * diagonal_element(op, k) * v)
     end
 end
 function LinearAlgebra.dot(
-    ::LOStructure, left::PDVec, op::AbstractHamiltonian, right::PDVec, w=nothing
+    ::LOStructure, left::PDVec, op::AbstractOperator, right::PDVec, w=nothing
 )
     # First two cases: only one vector is distrubuted. Avoid shuffling things around
     # by placing that one on the left to reduce the need for communication.
@@ -809,7 +813,7 @@ function LinearAlgebra.dot(
 end
 # Default variant: also called from other LOStructures.
 function LinearAlgebra.dot(
-    ::AdjointUnknown, t::PDVec, op::AbstractHamiltonian, source::PDVec, w=nothing
+    ::AdjointUnknown, t::PDVec, op::AbstractOperator, source::PDVec, w=nothing
 )
     if is_distributed(t)
         if isnothing(w)
@@ -823,7 +827,7 @@ function LinearAlgebra.dot(
     return dot_from_right(target, op, source)
 end
 
-function dot_from_right(target, op, source::PDVec)
+function Interfaces.dot_from_right(target, op, source::PDVec)
     T = typeof(zero(valtype(target)) * zero(valtype(source)) * zero(eltype(op)))
 
     result = sum(pairs(source); init=zero(T)) do (k, v)

diff --git a/src/DictVectors/projectors.jl b/src/DictVectors/projectors.jl
@@ -40,7 +40,7 @@ struct UniformProjector <: AbstractProjector end
 VectorInterface.inner(::UniformProjector, y::DVecOrVec) = sum(values(y))
 Base.getindex(::UniformProjector, add) = 1
 
-function VectorInterface.inner(::UniformProjector, op::AbstractHamiltonian, v::AbstractDVec)
+function LinearAlgebra.dot(::UniformProjector, op::AbstractOperator, v::AbstractDVec)
     return sum(pairs(v)) do (key, val)
         diag = diagonal_element(op, key) * val
         offdiag = sum(offdiagonals(op, key)) do (add, elem)

diff --git a/src/Hamiltonians/DensityMatrixDiagonal.jl b/src/Hamiltonians/DensityMatrixDiagonal.jl
@@ -21,6 +21,10 @@ struct DensityMatrixDiagonal{C} <: AbstractHamiltonian{Float64}
 end
 DensityMatrixDiagonal(mode; component=0) = DensityMatrixDiagonal{component}(mode)
 
+function allows_address_type(dmd::DensityMatrixDiagonal{C}, ::Type{A}) where {C, A}
+    return C ≤ num_components(A) && dmd.mode ≤ num_modes(A)
+end
+
 function diagonal_element(dmd::DensityMatrixDiagonal{1}, add::SingleComponentFockAddress)
     return float(find_mode(add, dmd.mode).occnum)
 end