chore: update the benchmark scripts and documentation

benches/_memory.py

@@ -94,6 +94,11 @@ def peak_rss_bytes() -> int:
 def reset_peak_rss() -> bool:
     """Reset ``VmHWM`` to the current RSS, starting a new measurement window.
 
+    This also resets ``getrusage(...).ru_maxrss`` (and so Julia's ``Sys.maxrss()``): both
+    report the same kernel field, ``mm->hiwater_rss``. Any process that opens a window
+    therefore loses ``ru_maxrss`` as a whole-run ceiling, and must take the maximum over
+    its windows instead.
+
     Returns:
         ``True`` if the reset took effect, ``False`` where ``/proc/self/clear_refs`` is
         unavailable (non-Linux, kernel < 4.0, or a restricted sandbox), in which case

benches/third_party/README.md

@@ -70,10 +70,14 @@ instead of squares with `--ny`.
 
 Two things to know when reading the output:
 
-- **The memory column is not one quantity.** Each backend reports what it exposes —
-  monoprop and cuPauliProp their own operator footprint, `PauliPropagation.jl`
-  `Base.summarysize`, and ppvm and Qiskit only a process-RSS proxy. `MEMORY_METRICS` in
-  `backends.py` records which is which, and every record also carries `peak_rss_MB`.
+- **The memory column is one quantity for every backend**: the peak resident set size over
+  each step, read from the kernel's `VmHWM` high-water mark and reset per step, so the
+  curves may be compared directly. Where a library also accounts for itself, that figure is
+  kept separately as `operator_memory_MB` / the `native_memory` series — those are *not*
+  commensurable across backends (one counts an operator, another an object graph, another a
+  device pool), so never plot them against each other. `HOST_MEMORY_METRIC` and
+  `OPERATOR_MEMORY_METRICS` in `backends.py` record which is which, and every record also
+  carries `peak_rss_MB` for the process lifetime.
 - **`PauliPropagation.jl` runs in its fastest documented configuration**, which is not its
   default: the `VectorPauliSum` container driven by `Performance.propagate!`, with
   coefficient truncation on. That combination needs the **dev branch (0.8.0)** — earlier

benches/third_party/bench_common.jl

@@ -0,0 +1,116 @@
+# Copyright 2026 Algorithmiq
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+# Shared across the third-party Julia benchmark scripts in this project (pauli_prop,
+# majorana_prop). Included via `include(joinpath(@__DIR__, "..", "bench_common.jl"))`.
+#
+# Line-for-line mirror of `benches/_memory.py`. The two must stay in step: a cross-language
+# memory comparison is only meaningful if both sides measure the same quantity the same way,
+# down to how the process is settled before the window opens.
+
+# Reset VmHWM to the current RSS, starting a new measurement window
+const _CLEAR_REFS_MM_HIWATER_RSS = "5\n"
+
+"""Return a `/proc/self` size field (kB -> bytes); 0 if unavailable."""
+function proc_field(path::AbstractString, key::AbstractString)::Int
+    try
+        for line in eachline(path)
+            if startswith(line, key)
+                return parse(Int, split(line)[2]) * 1024  # values are in kB
+            end
+        end
+    catch
+        return 0  # non-Linux or restricted /proc
+    end
+    return 0
+end
+
+"""Return this process's current resident set size (RSS) in bytes."""
+rss_bytes()::Int = proc_field("/proc/self/status", "VmRSS:")
+
+"""Return the kernel's high-water mark of this process's RSS, in bytes.
+
+VmHWM is maintained by the kernel on every RSS increase, so it is exact: unlike a polling
+sampler it cannot miss a transient that is allocated and freed between two observations.
+"""
+peak_rss_bytes()::Int = proc_field("/proc/self/status", "VmHWM:")
+
+"""Reset VmHWM to the current RSS, starting a new measurement window.
+
+Returns `false` where `/proc/self/clear_refs` is unavailable (non-Linux, kernel < 4.0, or a
+restricted sandbox), in which case VmHWM keeps counting from process start and callers must
+fall back.
+"""
+function reset_peak_rss()::Bool
+    try
+        write("/proc/self/clear_refs", _CLEAR_REFS_MM_HIWATER_RSS)
+    catch
+        return false
+    end
+    return true
+end
+
+"""Ask the C allocator to return unused heap pages to the OS."""
+function heap_trim()
+    try
+        ccall((:malloc_trim, "libc"), Cint, (Csize_t,), 0)
+    catch
+        nothing  # unsupported platform / allocator
+    end
+    return nothing
+end
+
+"""Return current RSS after collecting garbage and trimming the C heap."""
+function resting_rss_bytes()::Int
+    GC.gc()
+    heap_trim()
+    return rss_bytes()
+end
+
+"""Exact peak RSS over a measurement window, straight from the kernel.
+
+`start!` settles the process (`GC.gc()` + `malloc_trim`) and resets VmHWM to that floor, so
+the peak reported is this window's own and not an earlier window's retained garbage. The
+settling is what makes the number comparable against the Python side: without it the figure
+tracks the GC's willingness to return pages more than what the code needed.
+
+`exact` is `false` when the kernel would not reset the window (see `reset_peak_rss`); the
+peak then degrades to the RSS observed at `stop!`, which is a lower bound.
+"""
+mutable struct HighWaterMark
+    baseline_bytes::Int
+    peak_bytes::Int
+    exact::Bool
+end
+
+HighWaterMark() = HighWaterMark(0, 0, false)
+
+"""Open the window: settle, reset VmHWM, and record the floor."""
+function start!(hwm::HighWaterMark)
+    hwm.baseline_bytes = resting_rss_bytes()
+    hwm.exact = reset_peak_rss()
+    hwm.peak_bytes = hwm.baseline_bytes
+    return hwm
+end
+
+"""Close the window and latch the peak."""
+function stop!(hwm::HighWaterMark)
+    observed = hwm.exact ? peak_rss_bytes() : rss_bytes()
+    hwm.peak_bytes = max(hwm.baseline_bytes, observed)
+    return hwm
+end
+
+peak_mb(hwm::HighWaterMark) = hwm.peak_bytes / 1024^2
+baseline_mb(hwm::HighWaterMark) = hwm.baseline_bytes / 1024^2
+delta_mb(hwm::HighWaterMark) = (hwm.peak_bytes - hwm.baseline_bytes) / 1024^2

benches/third_party/majorana_prop/julia_hubbard1d_benchmark.jl

@@ -18,6 +18,8 @@ using BenchmarkTools
 using ArgParse
 using JSON
 
+include(joinpath(@__DIR__, "..", "bench_common.jl"))
+
 
 """CPU seconds (user + system) consumed by this process so far, summed over all threads.
 
@@ -30,17 +32,6 @@ function process_cpu_seconds()
 end
 
 
-"""Peak resident set size in MB, from VmHWM in /proc/self/status."""
-function peak_rss_mb()
-    for line in eachline("/proc/self/status")
-        if startswith(line, "VmHWM:")
-            return parse(Int, split(line)[2]) / 1024
-        end
-    end
-    return NaN
-end
-
-
 function experiment(N_spinful_sites, fock_state, circ_single, thetas_single, n_layers)
     site_index = N_spinful_sites ÷ 2
 
@@ -49,54 +40,91 @@ function experiment(N_spinful_sites, fock_state, circ_single, thetas_single, n_l
 
     obs = VectorMajoranaSum(MajoranaSum(N_spinful_sites, :nup, site_index))
 
-    res = zeros(n_layers + 1)
-    res[1] = overlapwithfock(obs, fock_state)
-
-
-    cpu_start = process_cpu_seconds()
-    stats = @timed for k = 1:n_layers
-        propagate!(circ_single, obs, thetas_single, min_abs_coeff=min_abs_coeff, max_unpaired=max_unpaired)
-        res[k+1] = overlapwithfock(obs, fock_state)
+    values = zeros(n_layers + 1)
+    term_counts = zeros(Int, n_layers + 1)
+    cumulative_runtimes = zeros(n_layers + 1)
+    memory_size = zeros(n_layers + 1)
+    native_memory_size = zeros(n_layers + 1)
+
+    sampler = HighWaterMark()
+
+    # Accumulated across the timed regions only: opening a memory window forces a full GC,
+    # so a counter spanning the whole loop would charge that settling cost to the workload
+    # and inflate both gc_seconds and busy_cores.
+    cpu_seconds = 0.0
+    gc_ns = 0
+
+    start!(sampler)
+    cpu_mark, gc_mark = process_cpu_seconds(), Base.gc_time_ns()
+    cumulative_runtimes[1] = @elapsed (values[1] = overlapwithfock(obs, fock_state))
+    cpu_seconds += process_cpu_seconds() - cpu_mark
+    gc_ns += Base.gc_time_ns() - gc_mark
+    stop!(sampler)
+    term_counts[1] = length(obs)
+    memory_size[1] = peak_mb(sampler)
+    native_memory_size[1] = Base.summarysize(obs) / 1024^2
+    for k = 1:n_layers
+        start!(sampler)
+        cpu_mark, gc_mark = process_cpu_seconds(), Base.gc_time_ns()
+        step_runtime = @elapsed propagate!(circ_single, obs, thetas_single, min_abs_coeff=min_abs_coeff, max_unpaired=max_unpaired)
+        cpu_seconds += process_cpu_seconds() - cpu_mark
+        gc_ns += Base.gc_time_ns() - gc_mark
+        stop!(sampler)
+        values[k+1] = overlapwithfock(obs, fock_state)
+        term_counts[k+1] = length(obs)
+        cumulative_runtimes[k+1] = cumulative_runtimes[k] + step_runtime
+        memory_size[k+1] = peak_mb(sampler)
+        native_memory_size[k+1] = Base.summarysize(obs) / 1024^2
     end
-    cpu_seconds = process_cpu_seconds() - cpu_start
-    memory_size = Base.summarysize(obs) / 1024^2
-    return (
-        res=res,
-        num_terms=length(obs),
-        runtime_seconds=stats.time,
-        gc_seconds=stats.gctime,
-        cpu_seconds=cpu_seconds,
-        memory_MB=memory_size,
+
+    total_runtime = cumulative_runtimes[end]
+    provenance = Dict(
+        "cpu_seconds" => cpu_seconds,
+        "gc_seconds" => gc_ns / 1e9,
+        "busy_cores" => total_runtime > 0 ? cpu_seconds / total_runtime : NaN,
         # Recorded, not hardcoded: only the Vector container dispatches into the
         # AcceleratedKernels path, so this is what proves the run was threaded at all.
-        container=string(nameof(typeof(obs))),
-    )
-
-
-end
-function save_result(output_path, N_spinful_sites, n_layers, result)
-    """Append one benchmark result as a JSON line, creating the parent directory if needed."""
-    record = Dict(
-        "n_spinful_sites" => N_spinful_sites,
-        "n_layers" => n_layers,
-        "num_terms" => result.num_terms,
-        "final_overlap" => result.res[end],
-        "runtime_seconds" => result.runtime_seconds,
-        "cpu_seconds" => result.cpu_seconds,
-        "busy_cores" => result.cpu_seconds / result.runtime_seconds,
-        "gc_seconds" => result.gc_seconds,
-        "num_threads" => Threads.nthreads(),
-        "container" => result.container,
-        "memory_MB" => result.memory_MB,
-        "peak_rss_MB" => peak_rss_mb(),
+        "container" => string(nameof(typeof(obs))),
         "library_version" => string(pkgversion(MajoranaPropagation)),
         "pauliprop_version" => string(pkgversion(PauliPropagation)),
         "host" => gethostname(),
     )
 
-    open(output_path, "a") do io
-        JSON.print(io, record)
-        println(io)
+    return values, term_counts, cumulative_runtimes, memory_size, native_memory_size, provenance
+
+
+end
+function save_result(output_path, source, N_spinful_sites, n_layers, term_counts, values, cumulative_runtimes, memory_size, native_memory_size, num_threads, provenance)
+    """Merge this run's per-step data into the shared results JSON file, keyed by source label."""
+    data = if isfile(output_path)
+        JSON.parsefile(output_path)
+    else
+        Dict(
+            "n_spinful_sites" => N_spinful_sites,
+            "n_layers" => n_layers,
+            "step_range" => collect(0:n_layers),
+            "num_threads" => Dict(),
+            "runtime_seconds" => Dict(),
+            "expectation_value" => Dict(),
+            "num_terms" => Dict(),
+            "memory_MB" => Dict(),
+            "native_memory_MB" => Dict(),
+        )
+    end
+    data["num_threads"] = get(data, "num_threads", Dict())
+    data["num_threads"][source] = num_threads
+    data["runtime_seconds"][source] = cumulative_runtimes
+    data["expectation_value"][source] = values
+    data["num_terms"][source] = term_counts
+    data["memory_MB"][source] = memory_size
+    data["native_memory_MB"] = get(data, "native_memory_MB", Dict())
+    data["native_memory_MB"][source] = native_memory_size
+    data["provenance"] = get(data, "provenance", Dict())
+    data["provenance"][source] = provenance
+
+    mkpath(dirname(output_path))
+    open(output_path, "w") do io
+        JSON.print(io, data, 4)
     end
 end
 
@@ -105,28 +133,27 @@ function main(args)
     s = ArgParseSettings(description="Arguments for the 1D Hubbard model benchmark.")
 
     @add_arg_table! s begin
-        "--case", "-c"
-        help = "Case pair to run."
+        "--n-spins", "-n"
+        help = "Number of spinful sites."
+        arg_type = Int
+        default = 60
+        dest_name = "n_spins"
+        "--max-layers", "-l"
+        help = "Number of Trotter layers."
         arg_type = Int
-        default = 1
+        default = 20
+        dest_name = "max_layers"
         "--output", "-o"
-        help = "Path to the JSONL file results are appended to."
+        help = "Path to the shared JSON file results are merged into."
         arg_type = String
-        default = joinpath(@__DIR__, "julia_hubbard1d_benchmark_results.jsonl")
+        default = joinpath(@__DIR__, "results.json")
 
     end
 
     parsed_args = parse_args(s)
 
-    spin_layers_pairs = []
-    for i in [20, 40, 60]
-        for j in 10:2:18
-            push!(spin_layers_pairs, (i, j))
-        end
-    end
-
-    case_pair = parsed_args["case"]
-    N_spinful_sites, n_layers = spin_layers_pairs[case_pair]
+    N_spinful_sites = parsed_args["n_spins"]
+    n_layers = parsed_args["max_layers"]
 
     t = 1.
     U = 1.5
@@ -161,12 +188,11 @@ function main(args)
 
     println("Number of threads: $(Threads.nthreads())")
 
-    result = experiment(N_spinful_sites, fock_state, circ_single, thetas_single, n_layers)
-    busy_cores = result.cpu_seconds / result.runtime_seconds
-    println("$N_spinful_sites n_spin $n_layers layers $(result.num_terms) num_terms $(result.res[end]) final overlap $(result.runtime_seconds) seconds")
-    println("container $(result.container)  cpu $(round(result.cpu_seconds, digits=1)) s  busy_cores $(round(busy_cores, digits=2))  gc $(round(result.gc_seconds, digits=1)) s")
+    values, term_counts, cumulative_runtimes, memory_size, native_memory_size, provenance = experiment(N_spinful_sites, fock_state, circ_single, thetas_single, n_layers)
+    println("$N_spinful_sites n_spin $n_layers layers $(term_counts[end]) num_terms $(values[end]) final overlap $(cumulative_runtimes[end]) seconds")
+    println("container $(provenance["container"])  cpu $(round(provenance["cpu_seconds"], digits=1)) s  busy_cores $(round(provenance["busy_cores"], digits=2))  gc $(round(provenance["gc_seconds"], digits=1)) s")
 
-    save_result(parsed_args["output"], N_spinful_sites, n_layers, result)
+    save_result(parsed_args["output"], "MajoranaPropagation.jl", N_spinful_sites, n_layers, term_counts, values, cumulative_runtimes, memory_size, native_memory_size, Threads.nthreads(), provenance)
 
 end