chore: update the benchmark scripts and documentation

benches/third_party/bench_common.jl

@@ -0,0 +1,60 @@
+# 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"))`.
+
+if Threads.nthreads() == 1
+    @warn "Julia is running with only 1 thread: the peak-memory sampler runs as a background " *
+          "task, but a single-threaded Julia can only switch to it between steps, not while the " *
+          "propagation call itself is running, so the memory figures will undercount any " *
+          "transient spike freed before that call returns. Re-run with `julia --threads=auto` " *
+          "(or set JULIA_NUM_THREADS) for accurate peak-memory measurements."
+end
+
+# Linux-only: reads the kernel's live RSS for this process directly, in kB.
+function current_rss_bytes()::Int
+    for line in eachline("/proc/self/status")
+        if startswith(line, "VmRSS:")
+            return parse(Int, split(line)[2]) * 1024
+        end
+    end
+    return 0
+end
+
+# Tracks this process's peak RSS within a resettable window: the kernel's own high-water mark
+# (Sys.maxrss()) is monotonic for the whole process and can't be reset per step, so this instead
+# polls the *current* RSS from a background task and keeps the max seen since the last reset!.
+mutable struct RssPeakSampler
+    peak_bytes::Threads.Atomic{Int}
+    running::Threads.Atomic{Bool}
+end
+
+function start_sampler(interval_s::Float64=1e-3)
+    sampler = RssPeakSampler(Threads.Atomic{Int}(current_rss_bytes()), Threads.Atomic{Bool}(true))
+    Threads.@spawn begin
+        while sampler.running[]
+            rss = current_rss_bytes()
+            if rss > sampler.peak_bytes[]
+                sampler.peak_bytes[] = rss
+            end
+            sleep(interval_s)
+        end
+    end
+    return sampler
+end
+
+reset!(sampler::RssPeakSampler) = (sampler.peak_bytes[] = current_rss_bytes())
+peak_mb(sampler::RssPeakSampler) = sampler.peak_bytes[] / 1024^2
+stop!(sampler::RssPeakSampler) = (sampler.running[] = false)

benches/third_party/bench_common.py

@@ -0,0 +1,65 @@
+# 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 benchmark scripts in this project (pauli_prop, majorana_prop)."""
+
+from __future__ import annotations
+
+import threading
+
+import psutil
+
+
+class RssPeakSampler:
+    """Tracks this process's peak resident set size within a resettable window.
+
+    The kernel-tracked high-water mark (``resource.getrusage().ru_maxrss``) is monotonic for the
+    whole process lifetime and can never be reset, so it cannot isolate a single step's peak from
+    the steps around it. This instead polls the *current* RSS from a background thread at a high
+    fixed frequency and keeps the max seen since the last :meth:`reset`, so each step gets its own
+    peak, independent of what earlier or later steps did.
+
+    Caveat: this can only sample while the calling thread doesn't hold the GIL — a C extension
+    call that never releases it (e.g. monoprop's ``propagate()``) blocks this thread out for its
+    whole duration, so the sampler only catches whatever RSS growth is already visible by the time
+    control returns to Python.
+    """
+
+    def __init__(self, interval_s: float = 1e-3) -> None:
+        self._process = psutil.Process()
+        self._interval_s = interval_s
+        self._peak_bytes = 0
+        self._stop_event = threading.Event()
+        self._thread = threading.Thread(target=self._poll_loop, daemon=True)
+
+    def _poll_loop(self) -> None:
+        while not self._stop_event.wait(self._interval_s):
+            rss = self._process.memory_info().rss
+            if rss > self._peak_bytes:
+                self._peak_bytes = rss
+
+    def __enter__(self) -> RssPeakSampler:
+        self._thread.start()
+        return self
+
+    def __exit__(self, *exc_info: object) -> None:
+        self._stop_event.set()
+        self._thread.join()
+
+    def reset(self) -> None:
+        """Start a new measurement window, floored at the current RSS."""
+        self._peak_bytes = self._process.memory_info().rss
+
+    def peak_mb(self) -> float:
+        return self._peak_bytes / 1024**2

benches/third_party/majorana_prop/julia_hubbard1d_benchmark.jl

@@ -17,6 +17,8 @@ using BenchmarkTools
 using ArgParse
 using JSON
 
+include(joinpath(@__DIR__, "..", "bench_common.jl"))
+
 
 function experiment(N_spinful_sites, fock_state, circ_single, thetas_single, n_layers)
     site_index = N_spinful_sites ÷ 2
@@ -26,33 +28,64 @@ 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)
+    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 = start_sampler()
+
+    reset!(sampler)
+    cumulative_runtimes[1] = @elapsed (values[1] = overlapwithfock(obs, fock_state))
+    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
+        reset!(sampler)
+        step_runtime = @elapsed propagate!(circ_single, obs, thetas_single, min_abs_coeff=min_abs_coeff, max_unpaired=max_unpaired)
+        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
 
+    stop!(sampler)
 
-    loop_elapsed = @elapsed 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)
-    end
-    memory_size = Base.summarysize(obs) / 1024^2
-    return res, length(obs), loop_elapsed, memory_size
+    return values, term_counts, cumulative_runtimes, memory_size, native_memory_size
 
 
 end
-function save_result(output_path, N_spinful_sites, n_layers, obs_length, final_res, loop_elapsed, memory_size)
-    """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" => obs_length,
-        "final_overlap" => final_res,
-        "runtime_seconds" => loop_elapsed,
-        "memory_MB" => memory_size,
-    )
-
-    open(output_path, "a") do io
-        JSON.print(io, record)
-        println(io)
+function save_result(output_path, source, N_spinful_sites, n_layers, term_counts, values, cumulative_runtimes, memory_size, native_memory_size, num_threads)
+    """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
+
+    mkpath(dirname(output_path))
+    open(output_path, "w") do io
+        JSON.print(io, data, 4)
     end
 end
 
@@ -61,28 +94,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 = 1
+        default = 60
+        dest_name = "n_spins"
+        "--max-layers", "-l"
+        help = "Number of Trotter layers."
+        arg_type = Int
+        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
@@ -117,11 +149,10 @@ function main(args)
 
     println("Number of threads: $(Threads.nthreads())")
 
-    res, obs_length, loop_elapsed, memory_size = experiment(N_spinful_sites, fock_state, circ_single, thetas_single, n_layers)
-    final_res = res[end]
-    println("$N_spinful_sites n_spin $n_layers layers $obs_length num_terms $final_res final overlap $loop_elapsed seconds")
+    values, term_counts, cumulative_runtimes, memory_size, native_memory_size = 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")
 
-    save_result(parsed_args["output"], N_spinful_sites, n_layers, obs_length, final_res, loop_elapsed, memory_size)
+    save_result(parsed_args["output"], "MajoranaPropagation.jl", N_spinful_sites, n_layers, term_counts, values, cumulative_runtimes, memory_size, native_memory_size, Threads.nthreads())
 
 end