14949967345

Committed 10 May 2025 11:08PM UTC coverage: 94.621% (-0.09%) from 94.706%

Build # 14949967345

Build Type

Pull #314

github

Committed by

joachimbrand

Commit Message

VMSize and get_vmsize

Pull Request Pull Request #314: Provide guidance for exact diagonalization algorithms

Run Details

76 of 84 new or added lines in 10 files covered. (90.48%)

1 existing line in 1 file now uncovered.

6983 of 7380 relevant lines covered (94.62%)

11371990.96 hits per line

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

94.44

/src/strategies_and_params/poststepstrategy.jl

"""
    PostStepStrategy

Subtypes of `PostStepStrategy` can be used to perform arbitrary computation on a single
state after an FCIQMC step is finished and report the results.

# Implemented strategies:

* [`ProjectedEnergy`](@ref)
* [`Projector`](@ref)
* [`SignCoherence`](@ref)
* [`WalkerLoneliness`](@ref)
* [`Timer`](@ref)

Note: A tuple of multiple strategies can be passed to [`ProjectorMonteCarloProblem`](@ref).
In that case, all reported column names must be distinct.

# Interface:

A subtype of this type must implement
[`post_step_action(::PostStepStrategy, ::SingleState, step::Int)`](@ref).
"""
abstract type PostStepStrategy end

"""
    post_step_action(::PostStepStrategy, ::SingleState, step) -> kvpairs

Compute statistics after FCIQMC step. Should return a tuple of `:key => value` pairs.
This function is only called every [`reporting_interval`](@ref) steps, as defined by the
`ReportingStrategy`.

See also [`PostStepStrategy`](@ref), [`ReportingStrategy`](@ref).
"""
post_step_action

# When strategies are a Tuple, apply all of them.
function post_step_action(::Tuple{}, _, _)
    return ()
end
function post_step_action((t,ts...)::Tuple, single_state, step)
    head = post_step_action(t, single_state, step)
    rest = post_step_action(ts, single_state, step)
    return (head..., rest...)
end

"""
    Projector(name=projector) <: PostStepStrategy

After each step, compute `dot(projector, dv)` and report it in the `DataFrame` under `name`.
`projector` can be an [`AbstractDVec`](@ref), or an [`AbstractProjector`](@ref).
"""
struct Projector{P} <: PostStepStrategy
    name::Symbol
    projector::P
end
function Projector(;kwarg...)
    length(kwarg) ≠ 1 && throw(
        ArgumentError("exactly one keyword argument must be passed to `Projector`")
    )
    return Projector(only(keys(kwarg)), freeze(only(values(kwarg))))
end

function post_step_action(p::Projector, single_state, _)
    return (p.name => dot(p.projector, single_state.v),)
end

"""
    ProjectedEnergy(hamiltonian, projector; hproj=:hproj, vproj=:vproj) <: PostStepStrategy

After every step, compute `hproj = dot(projector, hamiltonian, dv)` and `vproj =
dot(projector, dv)`, where `dv` is the instantaneous coefficient vector.  `projector` can be
an [`AbstractDVec`](@ref), or an [`AbstractProjector`](@ref).

Reports to columns `hproj` and `vproj`, which can be used to compute projective energy,
e.g. with [`projected_energy`](@ref). The keyword arguments `hproj` and
`vproj` can be used to change the names of these columns. This can be used to make the names
unique when computing projected energies with different projectors in the same run.

See also [`projected_energy`](@ref), [`ratio_of_means`](@ref), [`mixed_estimator`](@ref),
and [`PostStepStrategy`](@ref).
"""
struct ProjectedEnergy{H,P,Q} <: PostStepStrategy
    vproj_name::Symbol
    hproj_name::Symbol
    ham::H
    vproj::P
    hproj::Q
end

function ProjectedEnergy(
    hamiltonian, projector;
    vproj=:vproj, hproj=:hproj
)
    hproj_vec = compute_hproj(hamiltonian, projector)
    return ProjectedEnergy(vproj, hproj, hamiltonian, freeze(projector), hproj_vec)
end
compute_hproj(hamiltonian, projector::AbstractProjector) = nothing
# compute `dot` products with `AbstractProjector`s lazily
function compute_hproj(hamiltonian, projector)
    return compute_hproj(LOStructure(hamiltonian), hamiltonian, projector)
end
function compute_hproj(::AdjointUnknown, hamiltonian, projector)
    @warn "$(typeof(hamiltonian)) has an unknown adjoint. This will be slow."
    return nothing
end
function compute_hproj(::LOStructure, ham, projector)
    return freeze(ham' * projector)
end

function post_step_action(p::ProjectedEnergy{<:Any,<:Any,Nothing}, single_state, _)
    return (
        p.vproj_name => dot(p.vproj, single_state.v),
        p.hproj_name => dot(p.vproj, p.ham, single_state.v),
    )
end
function post_step_action(p::ProjectedEnergy, single_state, _)
    return (
        p.vproj_name => dot(p.vproj, single_state.v),
        p.hproj_name => dot(p.hproj, single_state.v),
    )
end

"""
    SignCoherence(reference[; name=:coherence]) <: PostStepStrategy

After each step, compute the proportion of configurations that have the same sign as they do
in the `reference_dvec`. Reports to a column named `name`, which defaults to `coherence`.
"""
struct SignCoherence{R} <: PostStepStrategy
    name::Symbol
    reference::R
end
SignCoherence(ref; name=:coherence) = SignCoherence(name, ref)

function post_step_action(sc::SignCoherence, single_state, _)
    vector = single_state.v
    return (sc.name => coherence(valtype(vector), sc.reference, vector),)
end

function coherence(::Type{<:Real}, reference, vector)
    accumulator, overlap = mapreduce(+, pairs(vector); init=MultiScalar(0.0, 0)) do ((k, v))
        ref = reference[k]
        MultiScalar(Float64(sign(ref) * sign(v)), Int(!iszero(ref)))
    end
    return iszero(overlap) ? 0.0 : accumulator / overlap
end
function coherence(::Type{<:Complex}, reference, vector)
    z = MultiScalar(0.0 + 0im, 0)
    accumulator, overlap = mapreduce(+, pairs(vector); init=z) do ((k, v))
        ref = sign(reference[k])
        MultiScalar(
            ComplexF64(sign(real(v)) * sign(ref) + im * sign(imag(v)) * sign(ref)),
            Int(!iszero(ref))
        )
    end
    return iszero(overlap) ? 0.0 : accumulator / overlap
end

"""
    WalkerLoneliness(threshold=1) <: PostStepStrategy

After each step, compute the proportion of configurations that are occupied by at most
`threshold` walkers. Reports to a column named `loneliness`.
"""
struct WalkerLoneliness{T} <: PostStepStrategy
    threshold::T
end
WalkerLoneliness() = WalkerLoneliness(1)

function post_step_action(wl::WalkerLoneliness, single_state, _)
    vector = single_state.v
    return (:loneliness => loneliness(valtype(vector), vector, wl.threshold),)
end

function loneliness(::Type{<:Real}, vector, threshold)
    num_lonely = mapreduce(+, values(vector), init=0) do v
        abs(v) ≤ threshold
    end
    return num_lonely / length(vector)
end

function loneliness(::Type{<:Complex}, vector, threshold)
    num_lonely = mapreduce(+, values(vector), init=0 + 0im) do v
        (abs(real(v)) ≤ threshold) + im*(abs(imag(v)) ≤ threshold)
    end
    return num_lonely / length(vector)
end

"""
    Timer() <: PostStepStrategy

Record current time after every step. See `Base.Libc.time` for information on what time
is recorded.
"""
struct Timer <: PostStepStrategy end

post_step_action(::Timer, _, _) = (:time => time(),)

"""
    SingleParticleDensity(; save_every=1, component) <: PostStepStrategy

[`PostStepStrategy`](@ref)  to  compute the diagonal [`single_particle_density`](@ref).
It records a `Tuple` with the same `eltype` as the vector.

Computing the density at every time step can be expensive. This cost can be reduced by
setting the `save_every` argument to a higher value. If the value is set, a vector of zeros
is recorded when the saving is skipped.

If the address type has multiple components, the `component` argument can be used to compute
the density on a per-component basis.

The density is not normalized, and must be divided by the vector `norm(⋅,2)` squared.

# See also

* [`single_particle_density`](@ref)
* [`DensityMatrixDiagonal`](@ref)
"""
struct SingleParticleDensity <: PostStepStrategy
    save_every::Int
    component::Int

    SingleParticleDensity(;save_every=1, component=0) = new(save_every, component)
end

function post_step_action(d::SingleParticleDensity, single_state, step)
    component = d.component
    if component == 0
        name = :single_particle_density
    else
        name = Symbol("single_particle_density_", component)
    end
    vector = single_state.v
    if step % d.save_every == 0
        return (name => single_particle_density(vector; component),)
    else
        V = valtype(vector)
        M = num_modes(keytype(vector))
        return (name => ntuple(_ -> 0.0, Val(M)),)
    end
end

"""
    single_particle_density(dvec; component)
    single_particle_density(add; component)

Compute the diagonal single particle density of vector `dvec` or address `add`. If the
`component` argument is given, only that component of the addresses is taken into
account. The result is always normalized so that `sum(result) ≈ num_particles(address)`.

# Examples

```jldoctest
julia> v = DVec(fs"|⋅↑⇅↓⋅⟩" => 1.0, fs"|↓↓⋅↑↑⟩" => 0.5)
DVec{FermiFS2C{2, 2, 5, 4, FermiFS{2, 5, BitString{5, 1, UInt8}}, FermiFS{2, 5, BitString{5, 1, UInt8}}},Float64} with 2 entries, style = IsDeterministic{Float64}()
  fs"|↓↓⋅↑↑⟩" => 0.5
  fs"|⋅↑⇅↓⋅⟩" => 1.0

julia> single_particle_density(v)
(0.2, 1.0, 1.6, 1.0, 0.2)

julia> single_particle_density(v; component=1)
(0.0, 0.8, 0.8, 0.2, 0.2)
```

# See also

* [`SingleParticleDensity`](@ref)
"""
function single_particle_density(dvec; component=0)
    K = keytype(dvec)
    V = float(valtype(dvec))
    M = num_modes(K)

    result = sum(pairs(dvec); init=MultiScalar(ntuple(_ -> zero(V), Val(M)))) do (k, v)
        MultiScalar(abs2(v) .* single_particle_density(k; component))
    end
    return result.tuple ./ sum(abs2, dvec)
end

function single_particle_density(add::SingleComponentFockAddress; component=0)
    return float.(Tuple(onr(add)))
end
function single_particle_density(add::Union{CompositeFS}; component=0)
    if component == 0
        return float.(Tuple(sum(onr(add))))
    else
        return float.(Tuple(onr(add)[component]))
    end
end

"""
    VMSize(; human_readable=true) <: PostStepStrategy

After each step, compute the maximum resident set size (VM size) of the current process in
bytes. Reports to a column named `vm_size`.
The `human_readable` keyword argument can be used to format the output in a human-readable
way. The default is `true`, which formats the output as a string using `Base.format_bytes`.

See also [`PostStepStrategy`](@ref).
"""
struct VMSize{HumanReadable} <: PostStepStrategy end

VMSize(human_readable::Bool) = VMSize{human_readable}()
function VMSize(; human_readable=true)
    return VMSize(human_readable)
end

post_step_action(::VMSize{false}, _, _) = (:vm_size => get_vmsize(),)
# `base.format_bytes` is not a public function, so might break in the future.
post_step_action(::VMSize{true}, _, _) = (:vm_size => Base.format_bytes(get_vmsize()),)

1	"""
2	PostStepStrategy
3
4	Subtypes of `PostStepStrategy` can be used to perform arbitrary computation on a single
5	state after an FCIQMC step is finished and report the results.
6
7	# Implemented strategies:
8
9	* [`ProjectedEnergy`](@ref)
10	* [`Projector`](@ref)
11	* [`SignCoherence`](@ref)
12	* [`WalkerLoneliness`](@ref)
13	* [`Timer`](@ref)
14
15	Note: A tuple of multiple strategies can be passed to [`ProjectorMonteCarloProblem`](@ref).
16	In that case, all reported column names must be distinct.
17
18	# Interface:
19
20	A subtype of this type must implement
21	[`post_step_action(::PostStepStrategy, ::SingleState, step::Int)`](@ref).
22	"""
23	abstract type PostStepStrategy end
24
25	"""
26	post_step_action(::PostStepStrategy, ::SingleState, step) -> kvpairs
27
28	Compute statistics after FCIQMC step. Should return a tuple of `:key => value` pairs.
29	This function is only called every [`reporting_interval`](@ref) steps, as defined by the
30	`ReportingStrategy`.
31
32	See also [`PostStepStrategy`](@ref), [`ReportingStrategy`](@ref).
33	"""
34	post_step_action
35
36	# When strategies are a Tuple, apply all of them.
37	function post_step_action(::Tuple{}, _, _)	318,625✔
38	return ()	318,625✔
39	end
40	function post_step_action((t,ts...)::Tuple, single_state, step)	63,622✔
41	head = post_step_action(t, single_state, step)	63,889✔
42	rest = post_step_action(ts, single_state, step)	63,889✔
43	return (head..., rest...)	63,622✔
44	end
45
46	"""
47	Projector(name=projector) <: PostStepStrategy
48
49	After each step, compute `dot(projector, dv)` and report it in the `DataFrame` under `name`.
50	`projector` can be an [`AbstractDVec`](@ref), or an [`AbstractProjector`](@ref).
51	"""
52	struct Projector{P} <: PostStepStrategy
53	name::Symbol	8✔
54	projector::P
55	end
56	function Projector(;kwarg...)	20✔
57	length(kwarg) ≠ 1 && throw(	10✔
58	ArgumentError("exactly one keyword argument must be passed to `Projector`")
59	)
60	return Projector(only(keys(kwarg)), freeze(only(values(kwarg))))	13✔
61	end
62
63	function post_step_action(p::Projector, single_state, _)	15,098✔
64	return (p.name => dot(p.projector, single_state.v),)	15,098✔
65	end
66
67	"""
68	ProjectedEnergy(hamiltonian, projector; hproj=:hproj, vproj=:vproj) <: PostStepStrategy
69
70	After every step, compute `hproj = dot(projector, hamiltonian, dv)` and `vproj =
71	dot(projector, dv)`, where `dv` is the instantaneous coefficient vector. `projector` can be
72	an [`AbstractDVec`](@ref), or an [`AbstractProjector`](@ref).
73
74	Reports to columns `hproj` and `vproj`, which can be used to compute projective energy,
75	e.g. with [`projected_energy`](@ref). The keyword arguments `hproj` and
76	`vproj` can be used to change the names of these columns. This can be used to make the names
77	unique when computing projected energies with different projectors in the same run.
78
79	See also [`projected_energy`](@ref), [`ratio_of_means`](@ref), [`mixed_estimator`](@ref),
80	and [`PostStepStrategy`](@ref).
81	"""
82	struct ProjectedEnergy{H,P,Q} <: PostStepStrategy
83	vproj_name::Symbol	9✔
84	hproj_name::Symbol
85	ham::H
86	vproj::P
87	hproj::Q
88	end
89
90	function ProjectedEnergy(	18✔
91	hamiltonian, projector;
92	vproj=:vproj, hproj=:hproj
93	)
94	hproj_vec = compute_hproj(hamiltonian, projector)	16✔
95	return ProjectedEnergy(vproj, hproj, hamiltonian, freeze(projector), hproj_vec)	10✔
96	end
97	compute_hproj(hamiltonian, projector::AbstractProjector) = nothing	1✔
98	# compute `dot` products with `AbstractProjector`s lazily
99	function compute_hproj(hamiltonian, projector)	8✔
100	return compute_hproj(LOStructure(hamiltonian), hamiltonian, projector)	15✔
101	end
102	function compute_hproj(::AdjointUnknown, hamiltonian, projector)	×
103	@warn "$(typeof(hamiltonian)) has an unknown adjoint. This will be slow."	×
104	return nothing	×
105	end
106	function compute_hproj(::LOStructure, ham, projector)	8✔
107	return freeze(ham' * projector)	15✔
108	end
109
110	function post_step_action(p::ProjectedEnergy{<:Any,<:Any,Nothing}, single_state, _)	100✔
111	return (	200✔
112	p.vproj_name => dot(p.vproj, single_state.v),
113	p.hproj_name => dot(p.vproj, p.ham, single_state.v),
114	)
115	end
116	function post_step_action(p::ProjectedEnergy, single_state, _)	27,424✔
117	return (	27,424✔
118	p.vproj_name => dot(p.vproj, single_state.v),
119	p.hproj_name => dot(p.hproj, single_state.v),
120	)
121	end
122
123	"""
124	SignCoherence(reference[; name=:coherence]) <: PostStepStrategy
125
126	After each step, compute the proportion of configurations that have the same sign as they do
127	in the `reference_dvec`. Reports to a column named `name`, which defaults to `coherence`.
128	"""
129	struct SignCoherence{R} <: PostStepStrategy
130	name::Symbol	4✔
131	reference::R
132	end
133	SignCoherence(ref; name=:coherence) = SignCoherence(name, ref)	8✔
134
135	function post_step_action(sc::SignCoherence, single_state, _)	10,400✔
136	vector = single_state.v	10,400✔
137	return (sc.name => coherence(valtype(vector), sc.reference, vector),)	10,567✔
138	end
139
140	function coherence(::Type{<:Real}, reference, vector)	10,200✔
141	accumulator, overlap = mapreduce(+, pairs(vector); init=MultiScalar(0.0, 0)) do ((k, v))	10,200✔
142	ref = reference[k]	2,150,289✔
143	MultiScalar(Float64(sign(ref) * sign(v)), Int(!iszero(ref)))	2,148,307✔
144	end
145	return iszero(overlap) ? 0.0 : accumulator / overlap	10,367✔
146	end
147	function coherence(::Type{<:Complex}, reference, vector)	200✔
148	z = MultiScalar(0.0 + 0im, 0)	200✔
149	accumulator, overlap = mapreduce(+, pairs(vector); init=z) do ((k, v))	200✔
150	ref = sign(reference[k])	6,337✔
151	MultiScalar(	4,158✔
152	ComplexF64(sign(real(v)) * sign(ref) + im * sign(imag(v)) * sign(ref)),
153	Int(!iszero(ref))
154	)
155	end
156	return iszero(overlap) ? 0.0 : accumulator / overlap	400✔
157	end
158
159	"""
160	WalkerLoneliness(threshold=1) <: PostStepStrategy
161
162	After each step, compute the proportion of configurations that are occupied by at most
163	`threshold` walkers. Reports to a column named `loneliness`.
164	"""
165	struct WalkerLoneliness{T} <: PostStepStrategy
166	threshold::T	3✔
167	end
168	WalkerLoneliness() = WalkerLoneliness(1)	3✔
169
170	function post_step_action(wl::WalkerLoneliness, single_state, _)	10,200✔
171	vector = single_state.v	10,200✔
172	return (:loneliness => loneliness(valtype(vector), vector, wl.threshold),)	10,200✔
173	end
174
175	function loneliness(::Type{<:Real}, vector, threshold)	10,100✔
176	num_lonely = mapreduce(+, values(vector), init=0) do v	10,100✔
177	abs(v) ≤ threshold	2,146,392✔
178	end
179	return num_lonely / length(vector)	10,100✔
180	end
181
182	function loneliness(::Type{<:Complex}, vector, threshold)	100✔
183	num_lonely = mapreduce(+, values(vector), init=0 + 0im) do v	100✔
184	(abs(real(v)) ≤ threshold) + im*(abs(imag(v)) ≤ threshold)	2,079✔
185	end
186	return num_lonely / length(vector)	100✔
187	end
188
189	"""
190	Timer() <: PostStepStrategy
191
192	Record current time after every step. See `Base.Libc.time` for information on what time
193	is recorded.
194	"""
195	struct Timer <: PostStepStrategy end	2✔
196
197	post_step_action(::Timer, _, _) = (:time => time(),)	200✔
198
199	"""
200	SingleParticleDensity(; save_every=1, component) <: PostStepStrategy
201
202	[`PostStepStrategy`](@ref) to compute the diagonal [`single_particle_density`](@ref).
203	It records a `Tuple` with the same `eltype` as the vector.
204
205	Computing the density at every time step can be expensive. This cost can be reduced by
206	setting the `save_every` argument to a higher value. If the value is set, a vector of zeros
207	is recorded when the saving is skipped.
208
209	If the address type has multiple components, the `component` argument can be used to compute
210	the density on a per-component basis.
211
212	The density is not normalized, and must be divided by the vector `norm(⋅,2)` squared.
213
214	# See also
215
216	* [`single_particle_density`](@ref)
217	* [`DensityMatrixDiagonal`](@ref)
218	"""
219	struct SingleParticleDensity <: PostStepStrategy
220	save_every::Int
221	component::Int
222
223	SingleParticleDensity(;save_every=1, component=0) = new(save_every, component)	2✔
224	end
225
226	function post_step_action(d::SingleParticleDensity, single_state, step)	100✔
227	component = d.component	100✔
228	if component == 0	100✔
229	name = :single_particle_density	100✔
230	else
231	name = Symbol("single_particle_density_", component)	×
232	end
233	vector = single_state.v	100✔
234	if step % d.save_every == 0	100✔
235	return (name => single_particle_density(vector; component),)	50✔
236	else
237	V = valtype(vector)	50✔
238	M = num_modes(keytype(vector))	50✔
239	return (name => ntuple(_ -> 0.0, Val(M)),)	50✔
240	end
241	end
242
243	"""
244	single_particle_density(dvec; component)
245	single_particle_density(add; component)
246
247	Compute the diagonal single particle density of vector `dvec` or address `add`. If the
248	`component` argument is given, only that component of the addresses is taken into
249	account. The result is always normalized so that `sum(result) ≈ num_particles(address)`.
250
251	# Examples
252
253	```jldoctest
254	julia> v = DVec(fs"\|⋅↑⇅↓⋅⟩" => 1.0, fs"\|↓↓⋅↑↑⟩" => 0.5)
255	DVec{FermiFS2C{2, 2, 5, 4, FermiFS{2, 5, BitString{5, 1, UInt8}}, FermiFS{2, 5, BitString{5, 1, UInt8}}},Float64} with 2 entries, style = IsDeterministic{Float64}()
256	fs"\|↓↓⋅↑↑⟩" => 0.5
257	fs"\|⋅↑⇅↓⋅⟩" => 1.0
258
259	julia> single_particle_density(v)
260	(0.2, 1.0, 1.6, 1.0, 0.2)
261
262	julia> single_particle_density(v; component=1)
263	(0.0, 0.8, 0.8, 0.2, 0.2)
264	```
265
266	# See also
267
268	* [`SingleParticleDensity`](@ref)
269	"""
270	function single_particle_density(dvec; component=0)	112✔
271	K = keytype(dvec)	56✔
272	V = float(valtype(dvec))	56✔
273	M = num_modes(K)	56✔
274
275	result = sum(pairs(dvec); init=MultiScalar(ntuple(_ -> zero(V), Val(M)))) do (k, v)	738✔
276	MultiScalar(abs2(v) .* single_particle_density(k; component))	936✔
277	end
278	return result.tuple ./ sum(abs2, dvec)	56✔
279	end
280
281	function single_particle_density(add::SingleComponentFockAddress; component=0)	10,201✔
282	return float.(Tuple(onr(add)))	9,274✔
283	end
284	function single_particle_density(add::Union{CompositeFS}; component=0)	24✔
285	if component == 0	10✔
286	return float.(Tuple(sum(onr(add))))	4✔
287	else
288	return float.(Tuple(onr(add)[component]))	6✔
289	end
290	end
291
292	"""
293	VMSize(; human_readable=true) <: PostStepStrategy
294
295	After each step, compute the maximum resident set size (VM size) of the current process in
296	bytes. Reports to a column named `vm_size`.
297	The `human_readable` keyword argument can be used to format the output in a human-readable
298	way. The default is `true`, which formats the output as a string using `Base.format_bytes`.
299
300	See also [`PostStepStrategy`](@ref).
301	"""
302	struct VMSize{HumanReadable} <: PostStepStrategy end	1✔
303
304	VMSize(human_readable::Bool) = VMSize{human_readable}()	1✔
305	function VMSize(; human_readable=true)	2✔
306	return VMSize(human_readable)	1✔
307	end
308
NEW 309	post_step_action(::VMSize{false}, _, _) = (:vm_size => get_vmsize(),)	×
310	# `base.format_bytes` is not a public function, so might break in the future.
311	post_step_action(::VMSize{true}, _, _) = (:vm_size => Base.format_bytes(get_vmsize()),)	100✔

RimuQMC / Rimu.jl / 14949967345

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