diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index a2bffe2..85af682 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -34,6 +34,11 @@ jobs: BENCH_INPUT_CACHE_DIR=/tmp/smoke-cache OUT_DIR=/tmp/smoke/raw EXPORT_HDF5=/tmp/smoke-export cargo run --release --locked --bin mpo_mpo_aniso_patched + - name: smoke vertex Bethe-Salpeter contraction (downloads one 3 MB input) + run: > + BENCH_VERTEX_TAUS=1e-4 BENCH_RUNS=1 BENCH_WARMUPS=0 + BENCH_INPUT_CACHE_DIR=/tmp/smoke-cache OUT_DIR=/tmp/smoke/raw + cargo run --release --locked --bin vertex_bse - uses: julia-actions/setup-julia@v2 with: version: "1" diff --git a/AGENTS.md b/AGENTS.md index 812ff5b..0d8dfdd 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -4,10 +4,12 @@ Read `README.md` before changing case definitions, timing boundaries, input gene ## Maintained scope -The repository has exactly three benchmark cases: Fourier elementwise, Gaussian elementwise global versus patched, and Gaussian MPO-MPO contraction global versus patched. Do not add another case when an arm or metric in one of these cases answers the question. +The repository has exactly four benchmark cases: Fourier elementwise, Gaussian elementwise global versus patched, Gaussian MPO-MPO contraction global versus patched, and the vertex Bethe-Salpeter contraction. Do not add another case when an arm or metric in one of these cases answers the question. Gaussian inputs use independent two-dimensional interpolative decomposition per term, deterministic balanced pairwise add-and-compress, final relative-L2/SVD tolerance `1e-6`, and fixed patch cap 128. Do not add TCI generation, whole-mixture interpolation, isotropic controls, independent `R` sweeps, or a runtime patch-cap knob without an explicit design change. +Vertex inputs are the fixed files of the `vertex-bse-inputs-v1` release. Do not regenerate or replace them; new inputs need a new release tag. + ACI residual tolerance is not an L2 tolerance. Keep its internal metric explicit and compare it with other arms only through the common achieved external error. ## Timing @@ -18,7 +20,7 @@ Run release builds with one pinned CPU core and single-threaded Rayon and BLAS. A machine may carry one additional scaling profile named `-t`, produced with `BENCH_THREADS=n`, which pins a contiguous core range and lets the -same three cases run on `n` Rayon and BLAS threads. Its `run.yaml` records the +same four cases run on `n` Rayon and BLAS threads. Its `run.yaml` records the thread count and the core range. Treat it as scaling evidence only: the single-core profile stays the reference, and the two profiles are never merged. diff --git a/Cargo.lock b/Cargo.lock index 2c351fa..d92d693 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -1369,6 +1369,7 @@ name = "t4a-bench" version = "0.1.0" dependencies = [ "anyhow", + "hdf5-metno", "libm", "num-complex", "rand", diff --git a/Cargo.toml b/Cargo.toml index 7a01119..65b2e16 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -22,6 +22,7 @@ serde = { version = "1", features = ["derive"] } serde_json = "1" anyhow = "1" libm = "0.2" +hdf5-metno = { version = "0.12", default-features = false } [profile.release] debug = 0 # CARGO_PROFILE_RELEASE_DEBUG=true restores full debug info diff --git a/README.md b/README.md index 7fb143e..9ea4173 100644 --- a/README.md +++ b/README.md @@ -23,6 +23,7 @@ has the layout. | `elementwise_fourier` | Elementwise product of two random one-dimensional Fourier QTTs | naive, zip-up, fit, ACI | measured input χ | | `gaussian_elementwise` | Elementwise product of randomly rotated anisotropic Gaussian QTTs | global fit, patched fit, global ACI, patched ACI | compressed input χ | | `gaussian_mpo_contraction` | MPO-MPO contraction of the same Gaussian input family | global fit, patched fit | compressed input χ | +| `vertex_bse_contraction` | Bethe-Salpeter contraction of impurity-model vertices (fixed inputs, downloaded) | zip-up, fit | input tolerance | There is no independent `R` sweep. `R` is the smallest value that gives at least eight grid cells across the fixed minor-axis standard deviation. Gaussian count `N` and `R` are construction metadata, not analysis axes. @@ -34,6 +35,19 @@ The terms are combined by deterministic balanced pairwise reduction. Every pair The expensive raw input pair is cached atomically in `.cache/inputs/`. The cache key includes the schema version, pinned tensor4all-rs revision, `N`, `R`, width, aspect range, spacing, polynomial degree, interpolation tolerance, addition tolerance, and seed. Cases 2 and 3 share the same cache entry. Set `BENCH_INPUT_CACHE_REFRESH=1` to rebuild it or `BENCH_INPUT_CACHE_DIR` to move the cache. +## Vertex input + +Case 4 contracts two-particle vertices of the single-impurity Anderson model (parquet approximation, magnetic channel, +β = 1000, 1024 Matsubara frequencies per axis) as a Bethe-Salpeter equation does, +Φ(ν, ν', ω) = Σ_ν₁ Γ(ν, ν₁, ω) χ0(ν₁, ω) F(ν₁, ν', ω): `hadamard` weights Γ by the bubble χ0, and `partial_contract` +joins the result with F, timed together and truncating at τ/10 per bond. The inputs are fixed files, one per input +tolerance τ (`1e-4`, `1e-6`, `1e-8`), from the +[`vertex-bse-inputs-v1`](https://github.com/tensor4all/tensor4all-benchmark/releases/tag/vertex-bse-inputs-v1) +release, whose notes describe them; they are downloaded once into `.cache/inputs/vertex_bse/`. The error is the +relative L2 error against the exact Φ at 4096 points, and a run fails above 10 τ. On one thread τ = `1e-4` takes +17 s (zip-up) and 33 s (fit) with a 1.5 GB peak, and τ = `1e-6` 9 and 14 minutes with an 18 GB peak; the default +sweep runs these two, `1e-8` is opt-in. + ## Allocation counters Every record carries the allocation calls and the requested bytes of the timed @@ -95,7 +109,7 @@ Requirements are Rust, HDF5, and a BLAS/LAPACK implementation. A complete single BENCH_CPU_CORE=0 scripts/run_all.sh linux-epyc-7713p ``` -The script pins Rayon and common BLAS implementations and the tenferro thread count to one thread, runs the three maintained binaries in release mode, writes raw records under `result//raw/`, records machine metadata in `run.yaml`, and generates `report.md`. +The script pins Rayon and common BLAS implementations and the tenferro thread count to one thread, runs the four maintained binaries in release mode, writes raw records under `result//raw/`, records machine metadata in `run.yaml`, and generates `report.md`. Set `BENCH_THREADS=n` for the optional scaling sweep. It uses `n` threads on the contiguous core range starting at `BENCH_CPU_CORE`, and belongs in its own profile named `-t`. The single-core profile stays the source of truth. @@ -128,6 +142,8 @@ Case 2 additionally accepts `BENCH_ACI_TOL`, default `1e-8`. This is an ACI resi Case 1 accepts `BENCH_KS`, `BENCH_R`, `BENCH_TOL`, `BENCH_MAX_BOND`, `BENCH_RUNS`, `BENCH_WARMUPS`, `BENCH_SEED`, and `BENCH_ALGOS`. Its report is sorted by measured input χ even though mode count generates the instances. +Case 4 accepts `BENCH_VERTEX_TAUS` (default `1e-4,1e-6`; available `1e-4`, `1e-6`, `1e-8`), `BENCH_ALGOS` (default `zipup,fit`), `BENCH_RUNS`, `BENCH_WARMUPS`, `BENCH_INPUT_CACHE_DIR`, and `OUT_DIR`. + ## Validation ```bash diff --git a/result/README.md b/result/README.md index d1c4f4b..77bbe07 100644 --- a/result/README.md +++ b/result/README.md @@ -1,7 +1,7 @@ # result/ One directory per (hardware profile, measured tensor4all-rs revision), each -holding one complete sweep of the three maintained cases: +holding one complete sweep of the maintained cases: ```text result/INDEX.md generated table of contents over every sweep diff --git a/scripts/report.py b/scripts/report.py index 8a23d64..97c3282 100755 --- a/scripts/report.py +++ b/scripts/report.py @@ -1,5 +1,5 @@ #!/usr/bin/env python3 -"""Render and validate the three maintained benchmark cases.""" +"""Render and validate the four maintained benchmark cases.""" import json import sys @@ -10,6 +10,7 @@ "elementwise_fourier": {"naive", "zipup", "fit", "aci"}, "gaussian_elementwise": {"global_fit", "patched_fit", "global_aci", "patched_aci"}, "gaussian_mpo_contraction": {"global_fit", "patched_fit"}, + "vertex_bse_contraction": {"zipup", "fit"}, } @@ -47,7 +48,8 @@ def main(profile: Path) -> None: assert {record["algorithm"] for record in grouped[case]} == arms by_instance = defaultdict(list) for record in grouped[case]: - key = record["params"].get("k_max", record["params"].get("n_gauss")) + params = record["params"] + key = params.get("k_max", params.get("n_gauss", params.get("input_tolerance"))) by_instance[key].append(record["algorithm"]) for key, actual in by_instance.items(): assert Counter(actual) == Counter(arms), (case, key, actual) @@ -57,6 +59,10 @@ def main(profile: Path) -> None: assert params["input_l2_rtol"] == 1e-6 assert params["external_error_metric"] == "sampled_relative_l2" assert record["max_error"] <= 1e-4 + for record in grouped["vertex_bse_contraction"]: + params = record["params"] + assert params["external_error_metric"] == "reference_relative_l2" + assert record["max_error"] <= 10 * params["input_tolerance"] run = metadata(profile / "run.yaml") lines = [ "# Tensor4all benchmark results", @@ -96,6 +102,18 @@ def main(profile: Path) -> None: lines.append( f"| {record['input_max_bond_dim']} | {raw_chi} | {params['n_gauss']} | {params['r']} | {record['algorithm']} | {record['wall_time_median_secs']:.6f} | {record['max_error']:.3e} | {input_patches} | {input_patch_chi} | {record.get('n_patches', 1)} | {record.get('n_params', '')} | {alloc_calls(record)} | {alloc_mib(record)} | {speedup:.3f} |" ) + lines += [ + "", + "## Case 4: vertex Bethe-Salpeter contraction", + "", + "| input tolerance | input χ (Γ / F) | arm | time (s) | relative L2 vs reference | best attainable | output χ | parameters | allocations | allocated (MiB) |", + "|---:|---:|---|---:|---:|---:|---:|---:|---:|---:|", + ] + for record in sorted(grouped["vertex_bse_contraction"], key=lambda item: (-item["params"]["input_tolerance"], item["algorithm"])): + params = record["params"] + lines.append( + f"| {params['input_tolerance']:.0e} | {params['gamma_chi']} / {params['f_chi']} | {record['algorithm']} | {record['wall_time_median_secs']:.6f} | {record['max_error']:.3e} | {params['best_attainable_error']:.3e} | {record['output_max_bond_dim']} | {record.get('n_params', '')} | {alloc_calls(record)} | {alloc_mib(record)} |" + ) (profile / "report.md").write_text("\n".join(lines) + "\n") diff --git a/scripts/run_all.sh b/scripts/run_all.sh index 7dcf46f..4db3a9e 100755 --- a/scripts/run_all.sh +++ b/scripts/run_all.sh @@ -74,5 +74,6 @@ EOF run cargo run --release --locked --bin elementwise_fourier run cargo run --release --locked --bin elementwise_gauss2d_patched run cargo run --release --locked --bin mpo_mpo_aniso_patched +run cargo run --release --locked --bin vertex_bse python3 scripts/report.py "$out" python3 scripts/gen_index.py diff --git a/src/bin/vertex_bse.rs b/src/bin/vertex_bse.rs new file mode 100644 index 0000000..2a5ae09 --- /dev/null +++ b/src/bin/vertex_bse.rs @@ -0,0 +1,229 @@ +//! Case 4: the Bethe-Salpeter contraction Φ(ν,ν',ω) = Σ_ν₁ Γ(ν,ν₁,ω) χ0(ν₁,ω) F(ν₁,ν',ω) of impurity-model vertices, +//! zip-up and fit. Inputs: the fixed files of the `vertex-bse-inputs-v1` release (layout in its notes), downloaded once +//! into the input cache. + +use std::path::{Path, PathBuf}; +use std::time::Instant; + +use anyhow::{anyhow, ensure}; +use t4a_bench::harness::time_median; +use t4a_bench::record::{write_record, RunRecord, SCHEMA_VERSION}; +use tensor4all_core::SvdTruncationPolicy; +use tensor4all_core::{ColMajorArrayRef, DynIndex, IdxTensor, IndexLike, LinearizationOrder}; +use tensor4all_treetn::contraction::ContractionOptions; +use tensor4all_treetn::{hadamard, partial_contract, PartialContractionSpec, TreeTN}; + +type Tree = TreeTN; + +const CASE: &str = "vertex_bse_contraction"; +const URL: &str = + "https://github.com/tensor4all/tensor4all-benchmark/releases/download/vertex-bse-inputs-v1"; +const FILES: [(f64, &str); 3] = [ + (1e-4, "vertex_bse_tau1e-04.h5"), + (1e-6, "vertex_bse_tau1e-06.h5"), + (1e-8, "vertex_bse_tau1e-08.h5"), +]; +/// Quantics bits per frequency. +const R: usize = 10; +const MAX_BOND: usize = 4096; + +fn env_or(name: &str, default: T) -> T { + std::env::var(name) + .ok() + .and_then(|v| v.parse().ok()) + .unwrap_or(default) +} + +/// The cached input file, downloaded first if missing; returns the path and whether it was cached. +fn fetch(name: &str, cache: &Path) -> anyhow::Result<(PathBuf, bool)> { + let path = cache.join("vertex_bse").join(name); + if path.exists() { + return Ok((path, true)); + } + std::fs::create_dir_all(path.parent().unwrap())?; + let partial = path.with_extension("partial"); + let status = std::process::Command::new("curl") + .args(["-fsSL", "--retry", "3", "-o"]) + .arg(&partial) + .arg(format!("{URL}/{name}")) + .status()?; + ensure!(status.success(), "downloading {name} failed"); + std::fs::rename(&partial, &path)?; + Ok((path, false)) +} + +/// One train of the file with its fused sites split into legs of dimension 2 (the first leg fastest). +fn load(path: &str, name: &str, legs: &[&[DynIndex]]) -> anyhow::Result { + let fused = tensor4all_hdf5::load_mps(path, name)?.into_treetn(); + ensure!(fused.node_count() == R, "{name}: expected {R} sites"); + let mut tensors = Vec::with_capacity(R); + for site in 0..R { + let index = fused + .site_space(&site) + .and_then(|s| s.iter().next()) + .cloned(); + let index = index.ok_or_else(|| anyhow!("{name}: site {site} has no index"))?; + ensure!( + index.dim() == 1 << legs.len(), + "{name}: wrong site dimension" + ); + let tensor = fused.tensor(fused.node_index(&site).unwrap()).unwrap(); + ensure!(tensor.is_f64(), "{name}: not real double precision"); + let parts: Vec = legs.iter().map(|series| series[site].clone()).collect(); + tensors.push(tensor.unfuse_index(&index, &parts, LinearizationOrder::ColumnMajor)?); + } + Ok(TreeTN::from_tensors(tensors, (0..R).collect())?) +} + +fn series() -> Vec { + (0..R).map(|_| DynIndex::new_dyn(2)).collect() +} + +fn max_bond(tree: &Tree) -> usize { + tree.link_dims().into_iter().max().unwrap_or(1) +} + +fn main() -> anyhow::Result<()> { + let taus = std::env::var("BENCH_VERTEX_TAUS").unwrap_or_else(|_| "1e-4,1e-6".into()); + let algos = std::env::var("BENCH_ALGOS").unwrap_or_else(|_| "zipup,fit".into()); + let out_dir = PathBuf::from(env_or("OUT_DIR", "result/dev/raw".to_string())); + let cache = PathBuf::from(env_or("BENCH_INPUT_CACHE_DIR", ".cache/inputs".to_string())); + let (runs, warmups) = ( + env_or("BENCH_RUNS", 1usize), + env_or("BENCH_WARMUPS", 0usize), + ); + for tau in taus.split(',').map(|v| v.trim().parse::()) { + let tau = tau?; + let (_, name) = FILES + .iter() + .find(|(t, _)| (t / tau - 1.0).abs() < 1e-9) + .ok_or_else(|| anyhow!("no input file for tau {tau:e}"))?; + let (path, cache_hit) = fetch(name, &cache)?; + let started = Instant::now(); + let file = path.to_str().unwrap(); + // Γ(ν, ν₁, ω) and F(ν₁, ν', ω), fused site index b_ω + 2 b_ν + 4 b_ν'; χ0(ν₁, ω): b_ω + 2 b_ν₁ + let (g_omega, g_nu, g_nu1) = (series(), series(), series()); + let (f_omega, f_nu1, f_nup) = (series(), series(), series()); + let (b_omega, b_nu1) = (series(), series()); + let gamma = load(file, "Gamma", &[&g_omega, &g_nu, &g_nu1])?; + let full = load(file, "F", &[&f_omega, &f_nu1, &f_nup])?; + let bubble = load(file, "bubble", &[&b_omega, &b_nu1])?; + let h5 = hdf5_metno::File::open(&path)?; + let best: f64 = h5.attr("contraction_error_of_inputs")?.read_scalar()?; + let reference = |n: &str| h5.dataset(&format!("reference/{n}"))?.read_raw::(); + let points = [reference("omega")?, reference("nu")?, reference("nup")?]; + let exact: Vec = h5.dataset("reference/value")?.read_raw()?; + ensure!( + points.iter().all(|p| p.len() == exact.len()), + "inconsistent reference" + ); + let load_secs = started.elapsed().as_secs_f64(); + + let weight_pairs: Vec<_> = (0..R) + .flat_map(|s| { + [ + (g_omega[s].clone(), b_omega[s].clone()), + (g_nu1[s].clone(), b_nu1[s].clone()), + ] + }) + .collect(); + let join = PartialContractionSpec { + contract_pairs: (0..R) + .map(|s| (g_nu1[s].clone(), f_nu1[s].clone())) + .collect(), + diagonal_pairs: (0..R) + .map(|s| (g_omega[s].clone(), f_omega[s].clone())) + .collect(), + output_order: None, + }; + let rtol = 0.1 * tau; + let policy = SvdTruncationPolicy::new(rtol * rtol / (R - 1) as f64) + .with_squared_values() + .with_discarded_tail_sum(); + for algo in algos.split(',').map(str::trim) { + let options = match algo { + "zipup" => ContractionOptions::zipup(), + "fit" => ContractionOptions::fit() + .with_nfullsweeps(1) + .with_convergence_tol(rtol), + _ => anyhow::bail!("unknown arm {algo}"), + } + .with_svd_policy(policy) + .with_max_bond_dim(MAX_BOND); + let (phi, timing) = time_median(warmups, runs, || -> anyhow::Result { + let weighted = hadamard(&gamma, &bubble, &weight_pairs, &0, options.clone()) + .map_err(|e| anyhow!("{e}"))?; + Ok(partial_contract( + &weighted, + &full, + &join, + &0, + options.clone(), + )?) + }); + let phi = phi?; + + // relative L2 error against the exact values; the result has legs ω, ν of Γ and ν' of F + let indices: Vec = (0..R) + .flat_map(|s| [g_omega[s].clone(), g_nu[s].clone(), f_nup[s].clone()]) + .collect(); + let values: Vec = (0..exact.len()) + .flat_map(|p| (0..R).flat_map(move |s| (0..3).map(move |k| (p, s, k)))) + .map(|(p, s, k)| ((points[k][p] as usize) >> (R - 1 - s)) & 1) + .collect(); + let shape = [indices.len(), exact.len()]; + let got = phi.evaluate(&indices, ColMajorArrayRef::new(&values, &shape)?)?; + let diff: f64 = got + .iter() + .zip(&exact) + .map(|(g, e)| (g.real() - e).powi(2)) + .sum(); + let error = (diff / exact.iter().map(|e| e * e).sum::()).sqrt(); + ensure!( + error <= 10.0 * tau, + "{algo} at tau {tau:e}: error {error:.3e}" + ); + + let n_params = (0..R) + .map(|s| { + phi.tensor(phi.node_index(&s).unwrap()) + .unwrap() + .dims() + .iter() + .product::() + }) + .sum(); + let record = RunRecord { + schema_version: SCHEMA_VERSION, + case: CASE.into(), + algorithm: algo.into(), + params: serde_json::json!({ + "input_tolerance": tau, "input_file": name, "cache_hit": cache_hit, + "gamma_chi": max_bond(&gamma), "f_chi": max_bond(&full), "bubble_chi": max_bond(&bubble), + "best_attainable_error": best, "external_error_metric": "reference_relative_l2", + }), + seed: 0, + tolerance: rtol, + wall_time_median_secs: timing.median_secs, + wall_times_secs: timing.runs_secs, + max_error: error, + input_max_bond_dim: max_bond(&gamma).max(max_bond(&full)), + output_max_bond_dim: max_bond(&phi), + output_bond_dims: phi.link_dims(), + n_params: Some(n_params), + n_patches: None, + max_patch_bond: None, + rtol: Some(rtol), + input_build_secs: Some(load_secs), + n_allocations: Some(timing.alloc_calls_median), + allocated_bytes: Some(timing.alloc_bytes_median), + }; + write_record(&out_dir, &format!("{CASE}-{algo}-tau{tau:e}"), &record)?; + eprintln!( + "tau {tau:e} {algo}: {:.2} s, error {error:.3e} (best {best:.3e})", + timing.median_secs + ); + } + } + Ok(()) +}