From 114c8ad721ee760c579048d5b10ea83f45976b22 Mon Sep 17 00:00:00 2001 From: Hiroshi Shinaoka Date: Thu, 20 Aug 2026 10:17:17 +0900 Subject: [PATCH] Measure anisotropic MPO patch scaling to chi 900 --- Cargo.lock | 11 + Cargo.toml | 1 + README.md | 19 +- .../mpo_mpo_aniso_input_scaling.md | 51 ++ ...po_mpo_aniso_input-tci-n100000-chi900.json | 59 ++ ...po_mpo_aniso_input-tci-n125000-chi900.json | 59 ++ ...mpo_mpo_aniso_input-tci-n20480-chi458.json | 59 ++ .../mpo_mpo_aniso_input-tci-n4096-chi256.json | 59 ++ .../mpo_mpo_aniso_input-tci-n512-chi89.json | 59 ++ scripts/report_mpo_input_scaling.py | 113 ++++ src/bin/mpo_mpo_aniso_patched.rs | 336 ++++++++-- src/gaussian.rs | 633 +++++++++++++++++- src/patched_mpo.rs | 145 +++- 13 files changed, 1498 insertions(+), 106 deletions(-) create mode 100644 result/linux-epyc-7713p/mpo_mpo_aniso_input_scaling.md create mode 100644 result/linux-epyc-7713p/raw/mpo_mpo_aniso_input-tci-n100000-chi900.json create mode 100644 result/linux-epyc-7713p/raw/mpo_mpo_aniso_input-tci-n125000-chi900.json create mode 100644 result/linux-epyc-7713p/raw/mpo_mpo_aniso_input-tci-n20480-chi458.json create mode 100644 result/linux-epyc-7713p/raw/mpo_mpo_aniso_input-tci-n4096-chi256.json create mode 100644 result/linux-epyc-7713p/raw/mpo_mpo_aniso_input-tci-n512-chi89.json create mode 100644 scripts/report_mpo_input_scaling.py diff --git a/Cargo.lock b/Cargo.lock index de9bb9e..f0da780 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -1369,6 +1369,7 @@ dependencies = [ "tensor4all-aci", "tensor4all-core", "tensor4all-hdf5", + "tensor4all-interpolativeqtt", "tensor4all-itensorlike", "tensor4all-partitionedtreetn", "tensor4all-partitionedtt", @@ -1584,6 +1585,16 @@ dependencies = [ "thiserror 2.0.20", ] +[[package]] +name = "tensor4all-interpolativeqtt" +version = "0.2.0" +source = "git+https://github.com/tensor4all/tensor4all-rs?rev=9e9aedaebe0d3918b34dd399ff0981e337f3835b#9e9aedaebe0d3918b34dd399ff0981e337f3835b" +dependencies = [ + "tensor4all-simplett", + "tensor4all-tensorbackend", + "thiserror 2.0.20", +] + [[package]] name = "tensor4all-itensorlike" version = "0.2.0" diff --git a/Cargo.toml b/Cargo.toml index 5f3cce8..47b00dd 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -13,6 +13,7 @@ tensor4all-itensorlike = { git = "https://github.com/tensor4all/tensor4all-rs", tensor4all-partitionedtt = { git = "https://github.com/tensor4all/tensor4all-rs", rev = "9e9aedaebe0d3918b34dd399ff0981e337f3835b", package = "tensor4all-partitionedtt" } tensor4all-partitionedtreetn = { git = "https://github.com/tensor4all/tensor4all-rs", rev = "9e9aedaebe0d3918b34dd399ff0981e337f3835b", package = "tensor4all-partitionedtreetn" } tensor4all-quanticstci = { git = "https://github.com/tensor4all/tensor4all-rs", rev = "9e9aedaebe0d3918b34dd399ff0981e337f3835b", package = "tensor4all-quanticstci" } +tensor4all-interpolativeqtt = { git = "https://github.com/tensor4all/tensor4all-rs", rev = "9e9aedaebe0d3918b34dd399ff0981e337f3835b", package = "tensor4all-interpolativeqtt" } tensor4all-hdf5 = { git = "https://github.com/tensor4all/tensor4all-rs", rev = "9e9aedaebe0d3918b34dd399ff0981e337f3835b", package = "tensor4all-hdf5" } tensor4all-core = { git = "https://github.com/tensor4all/tensor4all-rs", rev = "9e9aedaebe0d3918b34dd399ff0981e337f3835b", package = "tensor4all-core" } num-complex = "0.4" diff --git a/README.md b/README.md index ce8079e..b73dd68 100644 --- a/README.md +++ b/README.md @@ -16,7 +16,7 @@ revision that produced it. The pinned `tensor4all-rs` revision lives in `Cargo.t | 3. `elementwise_gauss2d` | Elementwise product of two 2D quantics Gaussian mixtures at a fixed output budget, swept over bits per variable `R`, against the exact pointwise product | [description](#case-3-elementwise-product-of-2d-quantics-gaussian-mixtures) | [`src/bin/elementwise_gauss2d.rs`](src/bin/elementwise_gauss2d.rs) | [`result/mac-cpu/elementwise_gauss2d.md`](result/mac-cpu/elementwise_gauss2d.md) | [time](result/mac-cpu/elementwise_gauss2d-time.svg), [error](result/mac-cpu/elementwise_gauss2d-error.svg) | | 4. `elementwise_gauss2d_scaling` | Density-constant scaling study of case 3: how the quantics input rank `chi_in` grows with the number of Gaussians `N` when the box area grows proportionally to `N` | [description](#case-4-density-constant-scaling-of-the-quantics-rank) | [`src/bin/elementwise_gauss2d_scaling.rs`](src/bin/elementwise_gauss2d_scaling.rs) | [`result/mac-cpu/elementwise_gauss2d_scaling.md`](result/mac-cpu/elementwise_gauss2d_scaling.md) | [chi](result/mac-cpu/elementwise_gauss2d_scaling-chi.svg), [time](result/mac-cpu/elementwise_gauss2d_scaling-time.svg), [error](result/mac-cpu/elementwise_gauss2d_scaling-error.svg) | | 5. `elementwise_gauss2d_patched` | Patched (domain decomposed) elementwise product, controlled by a global relative tolerance instead of a fixed output budget, comparing patched and global representations at equal accuracy on size and time, over two instance families: anisotropic narrow spikes by default and case 4's smooth Gaussians on request | [description](#case-5-patched-elementwise-product-at-equal-accuracy) | [`src/bin/elementwise_gauss2d_patched.rs`](src/bin/elementwise_gauss2d_patched.rs) | [`result/mac-cpu/elementwise_gauss2d_patched.md`](result/mac-cpu/elementwise_gauss2d_patched.md) | [time](result/mac-cpu/elementwise_gauss2d_patched-time.svg), [error](result/mac-cpu/elementwise_gauss2d_patched-error.svg), [parameters](result/mac-cpu/elementwise_gauss2d_patched-params-aniso.svg) | -| 6. `mpo_mpo_aniso_patched` | Fit contraction of anisotropic Gaussian MPOs, comparing the global baseline with patched chain TreeTN | [description](#case-6-patched-mpo-mpo-fit-contraction) | [`src/bin/mpo_mpo_aniso_patched.rs`](src/bin/mpo_mpo_aniso_patched.rs) | [`result/linux-epyc-7713p/mpo_mpo_aniso_patched.md`](result/linux-epyc-7713p/mpo_mpo_aniso_patched.md) | [results table](result/linux-epyc-7713p/mpo_mpo_aniso_patched.md#results) | +| 6. `mpo_mpo_aniso_patched` | Fit contraction of anisotropic Gaussian MPOs, comparing the global baseline with patched chain TreeTN | [description](#case-6-patched-mpo-mpo-fit-contraction) | [`src/bin/mpo_mpo_aniso_patched.rs`](src/bin/mpo_mpo_aniso_patched.rs) | [contraction](result/linux-epyc-7713p/mpo_mpo_aniso_patched.md), [input scaling](result/linux-epyc-7713p/mpo_mpo_aniso_input_scaling.md) | [crossover](result/linux-epyc-7713p/mpo_mpo_aniso_patched.md#results), [patch scaling](result/linux-epyc-7713p/mpo_mpo_aniso_input_scaling.md#results) | ## Benchmark cases @@ -531,7 +531,13 @@ comparison that established numerical equivalence. The default `N = 512` instance targets the practical global-rank range discussed for case 5. The default per-patch cap is 128; cap 64 is intentionally excluded -from future case-6 measurements. The committed `N = 2048` crossover sweep is +from future case-6 measurements. The separate [input-scaling report](result/linux-epyc-7713p/mpo_mpo_aniso_input_scaling.md) +records the established TCI `1e-8` → global relative-L2/SVD `1e-6` → adaptive +cap-128 pipeline through compressed input χ 900, where each input has 120 +patches. For large TCI inputs the runner uses a spatially indexed evaluator with +a recorded rigorous absolute tail bound; sampled errors are checked against the +unmodified full mixture. Regenerate that report from its committed raw records +with `python3 scripts/report_mpo_input_scaling.py result/linux-epyc-7713p`. The committed `N = 2048` crossover sweep is linked in the [top-level case table](#cases-at-a-glance): patched fit is slower at input χ 192, faster at χ 224, and 1.60× faster at χ 256 on the recorded single-core profile. The global output cap is the measured pre-patching @@ -779,6 +785,15 @@ Environment knobs: | `BENCH_NGAUSS` | cases 2, 3 and 6 | `8` (cases 2 and 3), `512` (case 6) | number of Gaussians or anisotropic spikes per mixture | | `BENCH_INPUT_CACHE_DIR` | case 6 | `.cache/inputs` | directory for deterministic global input MPS cache entries | | `BENCH_INPUT_CACHE_REFRESH` | case 6 | unset | when set, rebuild and atomically replace the selected input cache entry | +| `BENCH_INPUT_GENERATOR` | case 6 | `tci` | input constructor: global fused `tci`, or per-Gaussian `multiscale` using principal-axis unsafe points | +| `BENCH_INPUT_TCI_RTOL` | case 6 | `1e-8` | relative TCI tolerance; also the interpolation tolerance supplied to the multiscale constructor | +| `BENCH_TCI_LOCAL_ABS_TOL` | case 6 | `1e-12` | rigorous pointwise absolute bound for tails omitted by the spatially indexed TCI callback | +| `BENCH_TCI_INITIAL_PIVOTS` | case 6 | `8` | deterministic mixture-center grid pivots used by TCI; random initial pivots are disabled | +| `BENCH_INPUT_POLY_DEGREE` | case 6 | `28` | polynomial degree for the multiscale single-Gaussian constructor | +| `BENCH_INPUT_ADD_RTOL` | case 6 | `1e-10` | per-bond SVD tolerance used while adding multiscale single-Gaussian QTTs | +| `BENCH_INPUT_SVD_RTOL` | case 6 | `1e-6` | global relative-L2 SVD tolerance applied before adaptive patching | +| `BENCH_INPUT_ONLY` | case 6 | unset | when set, write the input/patch record and skip contraction | +| `BENCH_INPUT_SANITY` | case 6 | `1e-4` | maximum sampled relative error accepted for each compressed input MPO | | `BENCH_BOX_L` | cases 2 and 3 | `6.0` | half-width `L` of the box `[-L, L]` | | `BENCH_NS` | cases 4 and 5 | `8,16,32,64` (case 4, and case 5 on the `smooth` family), `8,16,32,64,128,256,512` (case 5 on the default `aniso` family) | comma-separated Gaussian or spike counts `N` to sweep. Both cases derive `L` and `R` from `N`, so they ignore `BENCH_NGAUSS`, `BENCH_BOX_L` and `BENCH_RS` | | `BENCH_FAMILY` | case 5 | `aniso` | which instance family to sweep. `aniso` is `N` anisotropic narrow spikes at a fixed spacing-to-width ratio, the family the case is written for; `smooth` is case 4's isotropic Gaussians at constant density, unchanged, and it also switches `BENCH_NS` to case 4's four points. Recorded in the params of every record as `family`, and part of every record's filename, so one profile can hold both | diff --git a/result/linux-epyc-7713p/mpo_mpo_aniso_input_scaling.md b/result/linux-epyc-7713p/mpo_mpo_aniso_input_scaling.md new file mode 100644 index 0000000..bea4183 --- /dev/null +++ b/result/linux-epyc-7713p/mpo_mpo_aniso_input_scaling.md @@ -0,0 +1,51 @@ +# Anisotropic Gaussian MPO input and patch scaling + +This focused case-6 input study measures adaptive patching through compressed input bond dimension `chi_in ≈ 1000`. It uses the unchanged constant-density random anisotropic family (`sigma_minor = 0.05`, log-uniform aspect ratio in `[1, 8]`, independent angle, center, and positive weight per Gaussian). + +## Construction and accuracy + +The input pipeline is: + +1. fused two-variable TCI at relative tolerance `1e-8`; +2. global relative-L2 SVD truncation at `1e-6`; +3. adaptive input truncation with per-patch cap 128. + +TCI evaluates the same random mixture through a spatial index. A component is omitted only outside a radius where its individual Gaussian is below a common exponent threshold. Since all weights are positive, the sum of all omitted tails is bounded pointwise by `1e-12`. The committed sampled input errors include this bound and remain below `1e-5` after global SVD truncation. + +The alternative multiscale constructor builds one randomly rotated Gaussian directly. It marks points every half minor-axis standard deviation along the major-axis ridge as unsafe. The focused test covers `rho = 8`, a 45-degree rotation, and `R = 10`, and requires maximum sampled absolute error below `5e-8` at polynomial degree 28. Mixtures can be formed by adding these single-Gaussian QTTs with intermediate SVD compression. This path preserves the instance family but is not used for the large sweep because its per-Gaussian construction cost scales linearly with `N`; spatially indexed TCI reaches the requested rank regime efficiently. + +## Reproducible settings + +These input-only records were collected separately from the profile's earlier contraction sweep. The settings below, including core 0, are authoritative for these records rather than the contraction affinity in `run.yaml`. + +- Machine: AMD EPYC 7713P 64-Core Processor +- CPU affinity: core 0 (`taskset -c 0`) +- `RAYON_NUM_THREADS=1`, `OPENBLAS_NUM_THREADS=1`, `OMP_NUM_THREADS=1` +- `BENCH_INPUT_GENERATOR=tci` +- `BENCH_INPUT_TCI_RTOL=1e-8` +- `BENCH_TCI_LOCAL_ABS_TOL=1e-12` +- `BENCH_TCI_INITIAL_PIVOTS=8` (deterministic center-derived pivots) +- `BENCH_INPUT_SVD_RTOL=1e-6` +- `BENCH_PATCH_MAX_BOND=128` +- `BENCH_MAX_INPUT_CHI=1400` +- `BENCH_MAX_BOND=512, 512, 1280, 3072, 3584` for `N=512, 4096, 20480, 100000, 125000` +- one deterministic seed (`BENCH_SEED=0`) +- `BENCH_ERROR_SAMPLES=256` at every size + +Input generation, global SVD truncation, sampled validation, and adaptive patch preparation are all included in `preparation time`; no contraction is run (`BENCH_INPUT_ONLY=1`). + +## Results + +Each pair is `(left, right)` for the two independently drawn mixtures. + +| N | R | raw chi | compressed chi | global parameters | patches | max patch chi | patch parameters | sampled relative error | preparation time | +|---:|---:|---:|---:|---:|---:|---:|---:|---:|---:| +| 512 | 9 | 181 | 89 | (48,060, 47,608) | (1, 1) | (84, 84) | (45,056, 45,056) | (3.36e-06, 2.56e-06) | 5.08 s | +| 4,096 | 10 | 280 | 256 | (254,476, 255,736) | (6, 6) | (124, 124) | (352,704, 352,088) | (1.78e-06, 3.38e-06) | 46.36 s | +| 20,480 | 11 | 889 | 458 | (888,076, 886,408) | (30, 30) | (121, 121) | (1,775,032, 1,772,136) | (2.75e-06, 6.95e-06) | 263.99 s | +| 100,000 | 12 | 1,074 | 900 | (3,450,320, 3,454,660) | (120, 120) | (120, 121) | (9,866,732, 9,982,956) | (3.06e-06, 3.74e-06) | 1,568.20 s | +| 125,000 | 13 | 1,342 | 900 | (4,120,944, 4,113,404) | (120, 120) | (121, 123) | (11,553,076, 11,596,696) | (5.51e-06, 9.14e-06) | 2,539.63 s | + +The compressed-rank plateau reaches `chi_in = 900`, within 10% of the requested `chi_in ≈ 1000` regime. Patch count grows stepwise from 1 to 6, 30, and 120 while every patch remains below cap 128. The `N = 125,000`, `R = 13` point increases raw rank and parameter count but remains on the same compressed-rank and patch-count plateau, so increasing `N` alone does not produce a smooth patch-count curve. + +Raw JSON records are in [`raw/`](raw/) with filenames `mpo_mpo_aniso_input-tci-n*-chi*.json`. diff --git a/result/linux-epyc-7713p/raw/mpo_mpo_aniso_input-tci-n100000-chi900.json b/result/linux-epyc-7713p/raw/mpo_mpo_aniso_input-tci-n100000-chi900.json new file mode 100644 index 0000000..cb5b975 --- /dev/null +++ b/result/linux-epyc-7713p/raw/mpo_mpo_aniso_input-tci-n100000-chi900.json @@ -0,0 +1,59 @@ +{ + "box_l": 23.71708245126285, + "box_padding": 1.0, + "case": "mpo_mpo_aniso_input", + "input_add_rtol": 1e-10, + "input_build_secs": 1568.202631011, + "input_cache_hit": false, + "input_cache_key": "aniso-v5-tci-n100000-r12-e0-s4008000000000000-p3ff0000000000000-sig3fa999999999999a-rho4020000000000000-tol3e45798ee2308c3a-loc3d719799812dea11-piv8-cap3072-seed0", + "input_cache_load_secs": 0.0, + "input_chi": 900, + "input_compression_secs": 9.464503739, + "input_error_samples": 256, + "input_generator": "tci", + "input_params": [ + 3450320, + 3454660 + ], + "input_patch_counts": [ + 120, + 120 + ], + "input_patch_max_bonds": [ + 120, + 121 + ], + "input_patch_params": [ + 9866732, + 9982956 + ], + "input_poly_degree": 28, + "input_sampled_relative_errors": [ + 3.0555493578220235e-6, + 3.7404138963447223e-6 + ], + "input_sanity": 0.0001, + "input_svd_l2_rtol": 1e-6, + "input_tci_cap": 3072, + "input_tci_initial_pivots": 8, + "input_tci_local_abs_tol": 1e-12, + "input_tci_rtol": 1e-8, + "max_input_chi": 1400, + "n_gauss": 100000, + "np_over_chi_squared": [ + 4.2596543209876545, + 4.265012345679012 + ], + "patch_cap": 128, + "r": 12, + "r_extra": 0, + "raw_input_chi": 1074, + "raw_input_params": [ + 7305760, + 7308792 + ], + "rho_max": 8.0, + "seed": 0, + "sigma": 0.05, + "spacing": 3.0 +} \ No newline at end of file diff --git a/result/linux-epyc-7713p/raw/mpo_mpo_aniso_input-tci-n125000-chi900.json b/result/linux-epyc-7713p/raw/mpo_mpo_aniso_input-tci-n125000-chi900.json new file mode 100644 index 0000000..354302b --- /dev/null +++ b/result/linux-epyc-7713p/raw/mpo_mpo_aniso_input-tci-n125000-chi900.json @@ -0,0 +1,59 @@ +{ + "box_l": 26.516504294495537, + "box_padding": 1.0, + "case": "mpo_mpo_aniso_input", + "input_add_rtol": 1e-10, + "input_build_secs": 2539.628990155, + "input_cache_hit": false, + "input_cache_key": "aniso-v5-tci-n125000-r13-e0-s4008000000000000-p3ff0000000000000-sig3fa999999999999a-rho4020000000000000-tol3e45798ee2308c3a-loc3d719799812dea11-piv8-cap3584-seed0", + "input_cache_load_secs": 0.0, + "input_chi": 900, + "input_compression_secs": 13.723229085, + "input_error_samples": 256, + "input_generator": "tci", + "input_params": [ + 4120944, + 4113404 + ], + "input_patch_counts": [ + 120, + 120 + ], + "input_patch_max_bonds": [ + 121, + 123 + ], + "input_patch_params": [ + 11553076, + 11596696 + ], + "input_poly_degree": 28, + "input_sampled_relative_errors": [ + 5.5092801319520735e-6, + 9.140267928712483e-6 + ], + "input_sanity": 0.0001, + "input_svd_l2_rtol": 1e-6, + "input_tci_cap": 3584, + "input_tci_initial_pivots": 8, + "input_tci_local_abs_tol": 1e-12, + "input_tci_rtol": 1e-8, + "max_input_chi": 1400, + "n_gauss": 125000, + "np_over_chi_squared": [ + 5.087585185185185, + 5.078276543209877 + ], + "patch_cap": 128, + "r": 13, + "r_extra": 0, + "raw_input_chi": 1342, + "raw_input_params": [ + 9424632, + 9396972 + ], + "rho_max": 8.0, + "seed": 0, + "sigma": 0.05, + "spacing": 3.0 +} \ No newline at end of file diff --git a/result/linux-epyc-7713p/raw/mpo_mpo_aniso_input-tci-n20480-chi458.json b/result/linux-epyc-7713p/raw/mpo_mpo_aniso_input-tci-n20480-chi458.json new file mode 100644 index 0000000..7e164ee --- /dev/null +++ b/result/linux-epyc-7713p/raw/mpo_mpo_aniso_input-tci-n20480-chi458.json @@ -0,0 +1,59 @@ +{ + "box_l": 10.733126291998992, + "box_padding": 1.0, + "case": "mpo_mpo_aniso_input", + "input_add_rtol": 1e-10, + "input_build_secs": 263.992662118, + "input_cache_hit": false, + "input_cache_key": "aniso-v5-tci-n20480-r11-e0-s4008000000000000-p3ff0000000000000-sig3fa999999999999a-rho4020000000000000-tol3e45798ee2308c3a-loc3d719799812dea11-piv8-cap1280-seed0", + "input_cache_load_secs": 0.0, + "input_chi": 458, + "input_compression_secs": 1.802141383, + "input_error_samples": 256, + "input_generator": "tci", + "input_params": [ + 888076, + 886408 + ], + "input_patch_counts": [ + 30, + 30 + ], + "input_patch_max_bonds": [ + 121, + 121 + ], + "input_patch_params": [ + 1775032, + 1772136 + ], + "input_poly_degree": 28, + "input_sampled_relative_errors": [ + 2.7528033675655293e-6, + 6.954288054456006e-6 + ], + "input_sanity": 0.0001, + "input_svd_l2_rtol": 1e-6, + "input_tci_cap": 1280, + "input_tci_initial_pivots": 8, + "input_tci_local_abs_tol": 1e-12, + "input_tci_rtol": 1e-8, + "max_input_chi": 1400, + "n_gauss": 20480, + "np_over_chi_squared": [ + 4.2336911958200645, + 4.2257394023760035 + ], + "patch_cap": 128, + "r": 11, + "r_extra": 0, + "raw_input_chi": 889, + "raw_input_params": [ + 2315988, + 2322644 + ], + "rho_max": 8.0, + "seed": 0, + "sigma": 0.05, + "spacing": 3.0 +} \ No newline at end of file diff --git a/result/linux-epyc-7713p/raw/mpo_mpo_aniso_input-tci-n4096-chi256.json b/result/linux-epyc-7713p/raw/mpo_mpo_aniso_input-tci-n4096-chi256.json new file mode 100644 index 0000000..7c2b5ca --- /dev/null +++ b/result/linux-epyc-7713p/raw/mpo_mpo_aniso_input-tci-n4096-chi256.json @@ -0,0 +1,59 @@ +{ + "box_l": 4.800000000000001, + "box_padding": 1.0, + "case": "mpo_mpo_aniso_input", + "input_add_rtol": 1e-10, + "input_build_secs": 46.359608388, + "input_cache_hit": false, + "input_cache_key": "aniso-v5-tci-n4096-r10-e0-s4008000000000000-p3ff0000000000000-sig3fa999999999999a-rho4020000000000000-tol3e45798ee2308c3a-loc3d719799812dea11-piv8-cap512-seed0", + "input_cache_load_secs": 0.0, + "input_chi": 256, + "input_compression_secs": 0.288869464, + "input_error_samples": 256, + "input_generator": "tci", + "input_params": [ + 254476, + 255736 + ], + "input_patch_counts": [ + 6, + 6 + ], + "input_patch_max_bonds": [ + 124, + 124 + ], + "input_patch_params": [ + 352704, + 352088 + ], + "input_poly_degree": 28, + "input_sampled_relative_errors": [ + 1.781935608908162e-6, + 3.3819348731372854e-6 + ], + "input_sanity": 0.0001, + "input_svd_l2_rtol": 1e-6, + "input_tci_cap": 512, + "input_tci_initial_pivots": 8, + "input_tci_local_abs_tol": 1e-12, + "input_tci_rtol": 1e-8, + "max_input_chi": 1400, + "n_gauss": 4096, + "np_over_chi_squared": [ + 3.88299560546875, + 3.9022216796875 + ], + "patch_cap": 128, + "r": 10, + "r_extra": 0, + "raw_input_chi": 280, + "raw_input_params": [ + 510520, + 522008 + ], + "rho_max": 8.0, + "seed": 0, + "sigma": 0.05, + "spacing": 3.0 +} \ No newline at end of file diff --git a/result/linux-epyc-7713p/raw/mpo_mpo_aniso_input-tci-n512-chi89.json b/result/linux-epyc-7713p/raw/mpo_mpo_aniso_input-tci-n512-chi89.json new file mode 100644 index 0000000..2afb4c1 --- /dev/null +++ b/result/linux-epyc-7713p/raw/mpo_mpo_aniso_input-tci-n512-chi89.json @@ -0,0 +1,59 @@ +{ + "box_l": 1.6970562748477145, + "box_padding": 1.0, + "case": "mpo_mpo_aniso_input", + "input_add_rtol": 1e-10, + "input_build_secs": 5.078735227, + "input_cache_hit": false, + "input_cache_key": "aniso-v5-tci-n512-r9-e0-s4008000000000000-p3ff0000000000000-sig3fa999999999999a-rho4020000000000000-tol3e45798ee2308c3a-loc3d719799812dea11-piv8-cap512-seed0", + "input_cache_load_secs": 0.0, + "input_chi": 89, + "input_compression_secs": 0.082263936, + "input_error_samples": 256, + "input_generator": "tci", + "input_params": [ + 48060, + 47608 + ], + "input_patch_counts": [ + 1, + 1 + ], + "input_patch_max_bonds": [ + 84, + 84 + ], + "input_patch_params": [ + 45056, + 45056 + ], + "input_poly_degree": 28, + "input_sampled_relative_errors": [ + 3.3550290079883713e-6, + 2.5588681566275107e-6 + ], + "input_sanity": 0.0001, + "input_svd_l2_rtol": 1e-6, + "input_tci_cap": 512, + "input_tci_initial_pivots": 8, + "input_tci_local_abs_tol": 1e-12, + "input_tci_rtol": 1e-8, + "max_input_chi": 1400, + "n_gauss": 512, + "np_over_chi_squared": [ + 6.067415730337078, + 6.010352228254009 + ], + "patch_cap": 128, + "r": 9, + "r_extra": 0, + "raw_input_chi": 181, + "raw_input_params": [ + 133328, + 132608 + ], + "rho_max": 8.0, + "seed": 0, + "sigma": 0.05, + "spacing": 3.0 +} \ No newline at end of file diff --git a/scripts/report_mpo_input_scaling.py b/scripts/report_mpo_input_scaling.py new file mode 100644 index 0000000..a43893b --- /dev/null +++ b/scripts/report_mpo_input_scaling.py @@ -0,0 +1,113 @@ +#!/usr/bin/env python3 +"""Render the case-6 anisotropic MPO input scaling report from raw JSON records.""" + +import json +import sys +from pathlib import Path + + +def pair(values, integer=False): + if integer: + return f"({values[0]:,}, {values[1]:,})" + return f"({values[0]:.2e}, {values[1]:.2e})" + + +def main() -> None: + if len(sys.argv) != 2: + raise SystemExit("usage: report_mpo_input_scaling.py result/") + profile = Path(sys.argv[1]) + paths = sorted( + (profile / "raw").glob("mpo_mpo_aniso_input-tci-n*-chi*.json"), + key=lambda path: json.loads(path.read_text())["n_gauss"], + ) + if not paths: + raise SystemExit(f"no case-6 input records under {profile / 'raw'}") + records = [json.loads(path.read_text()) for path in paths] + expected_caps = {512: 512, 4096: 512, 20480: 1280, 100000: 3072, 125000: 3584} + assert {record["n_gauss"] for record in records} == set(expected_caps) + for record in records: + assert record["case"] == "mpo_mpo_aniso_input" + assert record["input_generator"] == "tci" + assert record["sigma"] == 0.05 + assert record["rho_max"] == 8.0 + assert record["spacing"] == 3.0 + assert record["box_padding"] == 1.0 + assert record["r_extra"] == 0 + assert record["seed"] == 0 + assert record["input_tci_cap"] == expected_caps[record["n_gauss"]] + assert record["input_tci_rtol"] == 1e-8 + assert record["input_tci_local_abs_tol"] == 1e-12 + assert record["input_tci_initial_pivots"] == 8 + assert record["input_svd_l2_rtol"] == 1e-6 + assert record["patch_cap"] == 128 + assert record["max_input_chi"] == 1400 + assert record["input_error_samples"] == 256 + assert not record["input_cache_hit"] + assert max(record["input_sampled_relative_errors"]) <= record["input_sanity"] + + lines = [ + "# Anisotropic Gaussian MPO input and patch scaling", + "", + "This focused case-6 input study measures adaptive patching through compressed input bond dimension `chi_in ≈ 1000`. It uses the unchanged constant-density random anisotropic family (`sigma_minor = 0.05`, log-uniform aspect ratio in `[1, 8]`, independent angle, center, and positive weight per Gaussian).", + "", + "## Construction and accuracy", + "", + "The input pipeline is:", + "", + "1. fused two-variable TCI at relative tolerance `1e-8`;", + "2. global relative-L2 SVD truncation at `1e-6`;", + "3. adaptive input truncation with per-patch cap 128.", + "", + "TCI evaluates the same random mixture through a spatial index. A component is omitted only outside a radius where its individual Gaussian is below a common exponent threshold. Since all weights are positive, the sum of all omitted tails is bounded pointwise by `1e-12`. The committed sampled input errors include this bound and remain below `1e-5` after global SVD truncation.", + "", + "The alternative multiscale constructor builds one randomly rotated Gaussian directly. It marks points every half minor-axis standard deviation along the major-axis ridge as unsafe. The focused test covers `rho = 8`, a 45-degree rotation, and `R = 10`, and requires maximum sampled absolute error below `5e-8` at polynomial degree 28. Mixtures can be formed by adding these single-Gaussian QTTs with intermediate SVD compression. This path preserves the instance family but is not used for the large sweep because its per-Gaussian construction cost scales linearly with `N`; spatially indexed TCI reaches the requested rank regime efficiently.", + "", + "## Reproducible settings", + "", + "These input-only records were collected separately from the profile's earlier contraction sweep. The settings below, including core 0, are authoritative for these records rather than the contraction affinity in `run.yaml`.", + "", + "- Machine: AMD EPYC 7713P 64-Core Processor", + "- CPU affinity: core 0 (`taskset -c 0`)", + "- `RAYON_NUM_THREADS=1`, `OPENBLAS_NUM_THREADS=1`, `OMP_NUM_THREADS=1`", + "- `BENCH_INPUT_GENERATOR=tci`", + "- `BENCH_INPUT_TCI_RTOL=1e-8`", + "- `BENCH_TCI_LOCAL_ABS_TOL=1e-12`", + "- `BENCH_TCI_INITIAL_PIVOTS=8` (deterministic center-derived pivots)", + "- `BENCH_INPUT_SVD_RTOL=1e-6`", + "- `BENCH_PATCH_MAX_BOND=128`", + "- `BENCH_MAX_INPUT_CHI=1400`", + "- `BENCH_MAX_BOND=512, 512, 1280, 3072, 3584` for `N=512, 4096, 20480, 100000, 125000`", + "- one deterministic seed (`BENCH_SEED=0`)", + "- `BENCH_ERROR_SAMPLES=256` at every size", + "", + "Input generation, global SVD truncation, sampled validation, and adaptive patch preparation are all included in `preparation time`; no contraction is run (`BENCH_INPUT_ONLY=1`).", + "", + "## Results", + "", + "Each pair is `(left, right)` for the two independently drawn mixtures.", + "", + "| N | R | raw chi | compressed chi | global parameters | patches | max patch chi | patch parameters | sampled relative error | preparation time |", + "|---:|---:|---:|---:|---:|---:|---:|---:|---:|---:|", + ] + for record in records: + lines.append( + f"| {record['n_gauss']:,} | {record['r']} | {record['raw_input_chi']:,} | " + f"{record['input_chi']:,} | {pair(record['input_params'], True)} | " + f"{pair(record['input_patch_counts'], True)} | " + f"{pair(record['input_patch_max_bonds'], True)} | " + f"{pair(record['input_patch_params'], True)} | " + f"{pair(record['input_sampled_relative_errors'])} | " + f"{record['input_build_secs']:,.2f} s |" + ) + lines += [ + "", + "The compressed-rank plateau reaches `chi_in = 900`, within 10% of the requested `chi_in ≈ 1000` regime. Patch count grows stepwise from 1 to 6, 30, and 120 while every patch remains below cap 128. The `N = 125,000`, `R = 13` point increases raw rank and parameter count but remains on the same compressed-rank and patch-count plateau, so increasing `N` alone does not produce a smooth patch-count curve.", + "", + "Raw JSON records are in [`raw/`](raw/) with filenames `mpo_mpo_aniso_input-tci-n*-chi*.json`.", + "", + ] + (profile / "mpo_mpo_aniso_input_scaling.md").write_text("\n".join(lines)) + + +if __name__ == "__main__": + main() diff --git a/src/bin/mpo_mpo_aniso_patched.rs b/src/bin/mpo_mpo_aniso_patched.rs index c19b8e3..15431a8 100644 --- a/src/bin/mpo_mpo_aniso_patched.rs +++ b/src/bin/mpo_mpo_aniso_patched.rs @@ -5,11 +5,18 @@ use std::path::PathBuf; use std::time::Instant; use tensor4all_hdf5::{append_mps, load_mps, save_mps}; +use tensor4all_itensorlike::TensorTrain; +use tensor4all_simplett::mpo::MPO; -use t4a_bench::gaussian::{to_quantics_mpo_field, AnisoMixture2D}; -use t4a_bench::harness::time_median; -use t4a_bench::mpo_contract::{max_rel_error_vs_aniso_grid, mpo_n_params}; -use t4a_bench::patched_mpo::{PatchedMpoOutput, PatchedMpoPair}; +use t4a_bench::gaussian::{ + grid_coord, to_quantics_mpo_field_with_pivots, AnisoMixture2D, Field2D, LocalizedAnisoField, +}; +use t4a_bench::harness::{index_to_bits, sample_grid_indices, time_median}; +use t4a_bench::mpo_contract::{max_rel_error_vs_aniso_grid, mpo_n_params, tensortrain_to_mpo}; +use t4a_bench::patched_mpo::{ + mpo_pair_to_tensortrains, tensortrain_n_params, truncate_tensortrain_l2, PatchedMpoOutput, + PatchedMpoPair, +}; use t4a_bench::record::{write_record, RunRecord, SCHEMA_VERSION}; fn env_or(key: &str, default: T) -> T { @@ -19,6 +26,88 @@ fn env_or(key: &str, default: T) -> T { .unwrap_or(default) } +fn tensortrain_pair_to_mpos( + left: &TensorTrain, + right: &TensorTrain, +) -> anyhow::Result<(MPO, MPO)> { + let left_sites = left.site_indices(); + let right_sites = right.site_indices(); + let mut x = Vec::with_capacity(left.len()); + let mut y = Vec::with_capacity(left.len()); + let mut z = Vec::with_capacity(left.len()); + for site in 0..left.len() { + let common = left_sites[site] + .iter() + .find(|index| right_sites[site].contains(index)) + .ok_or_else(|| anyhow::anyhow!("input site {site} has no shared index"))?; + let left_unique = left_sites[site] + .iter() + .find(|index| *index != common) + .ok_or_else(|| anyhow::anyhow!("left input site {site} has no unique index"))?; + let right_unique = right_sites[site] + .iter() + .find(|index| *index != common) + .ok_or_else(|| anyhow::anyhow!("right input site {site} has no unique index"))?; + x.push(left_unique.clone()); + y.push(common.clone()); + z.push(right_unique.clone()); + } + Ok(( + tensortrain_to_mpo(left, &x, &y)?, + tensortrain_to_mpo(right, &y, &z)?, + )) +} + +fn max_rel_input_error( + mpo: &MPO, + field: &impl Field2D, + r: usize, + box_l: f64, + n_samples: usize, + seed: u64, +) -> anyhow::Result { + let xs = sample_grid_indices(r, n_samples, seed); + let ys = sample_grid_indices(r, n_samples, seed.wrapping_add(1)); + let mut max_error = 0.0_f64; + let mut max_reference = 0.0_f64; + for (&ix, &iy) in xs.iter().zip(&ys) { + let x = grid_coord(ix, r, box_l); + let y = grid_coord(iy, r, box_l); + let indices = index_to_bits(ix, r) + .into_iter() + .zip(index_to_bits(iy, r)) + .flat_map(|(x_bit, y_bit)| [x_bit, y_bit]) + .collect::>(); + let reference = field.eval(x, y); + let value = mpo.evaluate(&indices)?; + max_error = max_error.max((value - reference).abs()); + max_reference = max_reference.max(reference.abs()); + } + Ok(max_error / max_reference.max(f64::MIN_POSITIVE)) +} + +fn center_grid_pivots( + mixture: &AnisoMixture2D, + count: usize, + r: usize, + box_l: f64, +) -> Vec> { + let grid_size = 1usize << r; + (0..count.min(mixture.centers.len())) + .map(|pivot| { + let component = pivot * mixture.centers.len() / count.min(mixture.centers.len()); + let (x, y) = mixture.centers[component]; + [x, y] + .into_iter() + .map(|coordinate| { + ((coordinate + box_l) * grid_size as f64 / (2.0 * box_l)).floor() as usize + }) + .map(|index| index.min(grid_size - 1)) + .collect() + }) + .collect() +} + fn main() -> anyhow::Result<()> { let n_gauss: usize = env_or("BENCH_NGAUSS", 512); let sigma: f64 = env_or("BENCH_ANISO_SIGMA", 0.05); @@ -27,6 +116,14 @@ fn main() -> anyhow::Result<()> { let box_padding: f64 = env_or("BENCH_BOX_PADDING", 1.0); let extra_bits: usize = env_or("BENCH_R_EXTRA", 0); let rtol: f64 = env_or("BENCH_RTOL", 1e-8); + let input_generator = std::env::var("BENCH_INPUT_GENERATOR").unwrap_or_else(|_| "tci".into()); + let input_tci_rtol: f64 = env_or("BENCH_INPUT_TCI_RTOL", 1e-8); + let input_poly_degree: usize = env_or("BENCH_INPUT_POLY_DEGREE", 28); + let input_add_rtol: f64 = env_or("BENCH_INPUT_ADD_RTOL", 1e-10); + let input_tci_local_abs_tol: f64 = env_or("BENCH_TCI_LOCAL_ABS_TOL", 1e-12); + let input_tci_initial_pivots: usize = env_or("BENCH_TCI_INITIAL_PIVOTS", 8); + let input_svd_l2_rtol: f64 = env_or("BENCH_INPUT_SVD_RTOL", 1e-6); + let input_only = std::env::var_os("BENCH_INPUT_ONLY").is_some(); let input_tci_cap: usize = env_or("BENCH_MAX_BOND", 384); let patch_cap: usize = env_or("BENCH_PATCH_MAX_BOND", 128); let max_input_chi: usize = env_or("BENCH_MAX_INPUT_CHI", 256); @@ -34,6 +131,7 @@ fn main() -> anyhow::Result<()> { let warmups: usize = env_or("BENCH_WARMUPS", 0); let seed: u64 = env_or("BENCH_SEED", 0); let sanity: f64 = env_or("BENCH_SANITY", 1e-4); + let input_sanity: f64 = env_or("BENCH_INPUT_SANITY", 1e-4); let n_error_samples: usize = env_or("BENCH_ERROR_SAMPLES", 128); let out_dir = PathBuf::from(std::env::var("OUT_DIR").unwrap_or_else(|_| "result/dev/raw".into())); @@ -45,6 +143,13 @@ fn main() -> anyhow::Result<()> { anyhow::ensure!(n_gauss > 0 && sigma > 0.0 && spacing > 0.0); anyhow::ensure!(box_padding >= 1.0 && rho_max >= 1.0); anyhow::ensure!(rtol.is_finite() && rtol >= 0.0); + anyhow::ensure!(input_generator == "tci" || input_generator == "multiscale"); + anyhow::ensure!(input_tci_rtol.is_finite() && input_tci_rtol >= 0.0); + anyhow::ensure!(input_poly_degree >= 2); + anyhow::ensure!(input_add_rtol.is_finite() && input_add_rtol >= 0.0); + anyhow::ensure!(input_tci_local_abs_tol.is_finite() && input_tci_local_abs_tol > 0.0); + anyhow::ensure!(input_tci_initial_pivots > 0); + anyhow::ensure!(input_svd_l2_rtol.is_finite() && input_svd_l2_rtol >= 0.0); anyhow::ensure!(input_tci_cap > 0 && patch_cap > 0 && runs > 0); // Same constant-density anisotropic family and resolution as case 5. @@ -55,43 +160,75 @@ fn main() -> anyhow::Result<()> { let g = AnisoMixture2D::random(n_gauss, inner_l, sigma, rho_max, seed.wrapping_add(2)); fs::create_dir_all(&cache_dir)?; - let cache_key = format!( - "aniso-mpo-v1-n{n_gauss}-r{r}-extra{extra_bits}-spacing{:016x}-padding{:016x}-inner{:016x}-box{:016x}-sigma{:016x}-rho{:016x}-rtol{:016x}-cap{input_tci_cap}-seed{seed}", - spacing.to_bits(), - box_padding.to_bits(), - inner_l.to_bits(), - box_l.to_bits(), - sigma.to_bits(), - rho_max.to_bits(), - rtol.to_bits(), - ); + let cache_key = if input_generator == "tci" { + format!( + "aniso-v5-tci-n{n_gauss}-r{r}-e{extra_bits}-s{:x}-p{:x}-sig{:x}-rho{:x}-tol{:x}-loc{:x}-piv{input_tci_initial_pivots}-cap{input_tci_cap}-seed{seed}", + spacing.to_bits(), + box_padding.to_bits(), + sigma.to_bits(), + rho_max.to_bits(), + input_tci_rtol.to_bits(), + input_tci_local_abs_tol.to_bits(), + ) + } else { + format!( + "aniso-v4-ms-n{n_gauss}-r{r}-e{extra_bits}-s{:x}-p{:x}-sig{:x}-rho{:x}-tol{:x}-d{input_poly_degree}-add{:x}-seed{seed}", + spacing.to_bits(), + box_padding.to_bits(), + sigma.to_bits(), + rho_max.to_bits(), + input_tci_rtol.to_bits(), + input_add_rtol.to_bits(), + ) + }; let cache_path = cache_dir.join(format!("{cache_key}.h5")); let dy = 2.0 * box_l / 2.0_f64.powi(r as i32); let build_started = Instant::now(); let cache_hit = cache_path.is_file() && !refresh_cache; - let (patched, input_cache_load_secs) = if cache_hit { + let (mut left_train, mut right_train, input_cache_load_secs) = if cache_hit { eprintln!("loading cached input {}", cache_path.display()); let load_started = Instant::now(); let path = cache_path.to_string_lossy(); let left = load_mps(&path, "left")?; let right = load_mps(&path, "right")?; - let load_secs = load_started.elapsed().as_secs_f64(); - ( - PatchedMpoPair::from_tensortrains(left, right, rtol, patch_cap)?, - load_secs, - ) + (left, right, load_started.elapsed().as_secs_f64()) } else { - eprintln!("building N={n_gauss} R={r} box_l={box_l:.6} input_tci_cap={input_tci_cap}"); - let (left, generated_dy) = to_quantics_mpo_field(&f, r, box_l, rtol, input_tci_cap)?; - let (right, _) = to_quantics_mpo_field(&g, r, box_l, rtol, input_tci_cap)?; - anyhow::ensure!((generated_dy - dy).abs() <= f64::EPSILON * dy.abs().max(1.0)); - let patched = PatchedMpoPair::new(&left, &right, rtol, patch_cap)?; + eprintln!( + "building generator={input_generator} N={n_gauss} R={r} box_l={box_l:.6} input_rtol={input_tci_rtol:.1e} input_tci_cap={input_tci_cap}" + ); + let (left, right) = if input_generator == "tci" { + let localized_f = LocalizedAnisoField::new(f.clone(), input_tci_local_abs_tol)?; + let localized_g = LocalizedAnisoField::new(g.clone(), input_tci_local_abs_tol)?; + let (left, generated_dy) = to_quantics_mpo_field_with_pivots( + &localized_f, + r, + box_l, + input_tci_rtol, + input_tci_cap, + center_grid_pivots(&f, input_tci_initial_pivots, r, box_l), + )?; + let (right, _) = to_quantics_mpo_field_with_pivots( + &localized_g, + r, + box_l, + input_tci_rtol, + input_tci_cap, + center_grid_pivots(&g, input_tci_initial_pivots, r, box_l), + )?; + anyhow::ensure!((generated_dy - dy).abs() <= f64::EPSILON * dy.abs().max(1.0)); + (left, right) + } else { + ( + f.to_multiscale_mpo(r, box_l, input_poly_degree, input_tci_rtol, input_add_rtol)?, + g.to_multiscale_mpo(r, box_l, input_poly_degree, input_tci_rtol, input_add_rtol)?, + ) + }; + let (left_train, right_train) = mpo_pair_to_tensortrains(&left, &right)?; let temporary = cache_dir.join(format!(".{cache_key}.{}.tmp", std::process::id())); let temporary_path = temporary.to_string_lossy(); - let (left_train, right_train) = patched.global_inputs(); let write_result = (|| -> anyhow::Result<()> { - save_mps(&temporary_path, "left", left_train)?; - append_mps(&temporary_path, "right", right_train)?; + save_mps(&temporary_path, "left", &left_train)?; + append_mps(&temporary_path, "right", &right_train)?; fs::rename(&temporary, &cache_path)?; Ok(()) })(); @@ -99,24 +236,120 @@ fn main() -> anyhow::Result<()> { let _ = fs::remove_file(&temporary); } write_result?; - (patched, 0.0) + (left_train, right_train, 0.0) }; - let input_build_secs = build_started.elapsed().as_secs_f64(); - let raw_input_chi = patched.input_max_bond(); anyhow::ensure!( - raw_input_chi <= max_input_chi, - "input chi {raw_input_chi} exceeds BENCH_MAX_INPUT_CHI={max_input_chi}" + left_train.len() == r && right_train.len() == r, + "cached input length does not match R={r}" + ); + let raw_input_chi = left_train.max_bond_dim().max(right_train.max_bond_dim()); + let raw_input_params = ( + tensortrain_n_params(&left_train), + tensortrain_n_params(&right_train), + ); + let compression_started = Instant::now(); + truncate_tensortrain_l2(&mut left_train, input_svd_l2_rtol)?; + truncate_tensortrain_l2(&mut right_train, input_svd_l2_rtol)?; + let input_compression_secs = compression_started.elapsed().as_secs_f64(); + let input_chi = left_train.max_bond_dim().max(right_train.max_bond_dim()); + anyhow::ensure!( + input_chi <= max_input_chi, + "compressed input chi {input_chi} exceeds BENCH_MAX_INPUT_CHI={max_input_chi}" + ); + let input_params = ( + tensortrain_n_params(&left_train), + tensortrain_n_params(&right_train), + ); + let (left_mpo, right_mpo) = tensortrain_pair_to_mpos(&left_train, &right_train)?; + let input_errors = ( + max_rel_input_error( + &left_mpo, + &f, + r, + box_l, + n_error_samples, + seed.wrapping_add(11), + )?, + max_rel_input_error( + &right_mpo, + &g, + r, + box_l, + n_error_samples, + seed.wrapping_add(12), + )?, ); - let input_chi = patched.input_max_bond(); - let input_patches = patched.input_patch_counts(); - let input_patch_bonds = patched.input_patch_max_bonds(); anyhow::ensure!( - input_chi == raw_input_chi, - "MPO bridge changed the input rank" + input_errors.0 <= input_sanity && input_errors.1 <= input_sanity, + "input sampled relative errors {input_errors:?} exceed {input_sanity:.3e}" ); + let patched = PatchedMpoPair::from_tensortrains_with_input_rtol( + left_train, + right_train, + input_svd_l2_rtol, + rtol, + patch_cap, + )?; + let input_build_secs = build_started.elapsed().as_secs_f64(); + let input_patches = patched.input_patch_counts(); + let input_patch_bonds = patched.input_patch_max_bonds(); + let input_patch_params = patched.input_patch_n_params(); eprintln!( - "prepared input_chi={input_chi} patches={input_patches:?} patch_chi={input_patch_bonds:?} build={input_build_secs:.3}s" + "prepared raw_chi={raw_input_chi} input_chi={input_chi} raw_params={raw_input_params:?} input_params={input_params:?} input_errors={input_errors:?} np_over_chi2=({:.3},{:.3}) patches={input_patches:?} patch_chi={input_patch_bonds:?} patch_params={input_patch_params:?} compress={input_compression_secs:.3}s build={input_build_secs:.3}s", + input_params.0 as f64 / (input_chi * input_chi) as f64, + input_params.1 as f64 / (input_chi * input_chi) as f64, ); + let input_record = serde_json::json!({ + "case": "mpo_mpo_aniso_input", + "n_gauss": n_gauss, + "r": r, + "box_l": box_l, + "sigma": sigma, + "rho_max": rho_max, + "spacing": spacing, + "box_padding": box_padding, + "r_extra": extra_bits, + "seed": seed, + "input_generator": input_generator, + "input_cache_key": cache_key, + "input_tci_rtol": input_tci_rtol, + "input_tci_cap": input_tci_cap, + "input_poly_degree": input_poly_degree, + "input_add_rtol": input_add_rtol, + "input_tci_local_abs_tol": input_tci_local_abs_tol, + "input_tci_initial_pivots": input_tci_initial_pivots, + "input_svd_l2_rtol": input_svd_l2_rtol, + "raw_input_chi": raw_input_chi, + "input_chi": input_chi, + "raw_input_params": raw_input_params, + "input_params": input_params, + "input_sampled_relative_errors": input_errors, + "input_sanity": input_sanity, + "input_error_samples": n_error_samples, + "max_input_chi": max_input_chi, + "np_over_chi_squared": [ + input_params.0 as f64 / (input_chi * input_chi) as f64, + input_params.1 as f64 / (input_chi * input_chi) as f64, + ], + "patch_cap": patch_cap, + "input_patch_counts": input_patches, + "input_patch_max_bonds": input_patch_bonds, + "input_patch_params": input_patch_params, + "input_cache_hit": cache_hit, + "input_cache_load_secs": input_cache_load_secs, + "input_compression_secs": input_compression_secs, + "input_build_secs": input_build_secs, + }); + fs::create_dir_all(&out_dir)?; + fs::write( + out_dir.join(format!( + "mpo_mpo_aniso_input-{input_generator}-n{n_gauss}-chi{input_chi}.json" + )), + serde_json::to_string_pretty(&input_record)?, + )?; + if input_only { + return Ok(()); + } let (global_train, global_timing) = time_median(warmups, runs, || { patched @@ -162,10 +395,17 @@ fn main() -> anyhow::Result<()> { "rho_max": rho_max, "spacing": spacing, "r_extra": extra_bits, + "input_tci_rtol": input_tci_rtol, "input_tci_cap": input_tci_cap, + "input_svd_l2_rtol": input_svd_l2_rtol, + "raw_input_chi": raw_input_chi, + "raw_input_params": [raw_input_params.0, raw_input_params.1], + "input_params": [input_params.0, input_params.1], + "input_compression_secs": input_compression_secs, "input_patch_cap": patch_cap, "input_patch_counts": [input_patches.0, input_patches.1], "input_patch_max_bonds": [input_patch_bonds.0, input_patch_bonds.1], + "input_patch_params": [input_patch_params.0, input_patch_params.1], "global_contract_max_bond": input_chi, "patched_contract_max_bond": patch_cap, "fit_nfullsweeps": 1, @@ -283,6 +523,24 @@ fn write_arm( mod tests { use super::*; use std::time::{SystemTime, UNIX_EPOCH}; + use t4a_bench::gaussian::to_quantics_mpo_field; + + #[test] + fn l2_compression_reduces_or_preserves_input_size() { + let field = AnisoMixture2D::random(4, 1.0, 0.2, 3.0, 9); + let other = AnisoMixture2D::random(4, 1.0, 0.2, 3.0, 10); + let (left, _) = to_quantics_mpo_field(&field, 5, 1.0, 1.0e-8, 64).unwrap(); + let (right, _) = to_quantics_mpo_field(&other, 5, 1.0, 1.0e-8, 64).unwrap(); + let (mut left, mut right) = mpo_pair_to_tensortrains(&left, &right).unwrap(); + let before = (tensortrain_n_params(&left), tensortrain_n_params(&right)); + + truncate_tensortrain_l2(&mut left, 1.0e-6).unwrap(); + truncate_tensortrain_l2(&mut right, 1.0e-6).unwrap(); + + assert!(tensortrain_n_params(&left) <= before.0); + assert!(tensortrain_n_params(&right) <= before.1); + PatchedMpoPair::from_tensortrains_with_input_rtol(left, right, 1.0e-6, 1.0e-8, 16).unwrap(); + } #[test] fn cached_inputs_round_trip_with_shared_index_identity() { diff --git a/src/gaussian.rs b/src/gaussian.rs index e7f78b9..1b8a3f3 100644 --- a/src/gaussian.rs +++ b/src/gaussian.rs @@ -1,13 +1,21 @@ //! 2D Gaussian mixtures: analytic y-integral and quantics MPO construction. +use std::collections::HashMap; + use rand::Rng; use rand::SeedableRng; use rand_chacha::ChaCha8Rng; +use tensor4all_interpolativeqtt::{ + direct_product_core_tensors, interpolate_multi_scale, interpolate_multi_scale_nd, + InterpolativeQttOptions, +}; use tensor4all_quanticstci::{ quanticscrossinterpolate, DiscretizedGrid, QtciOptions, UnfoldingScheme, }; use tensor4all_simplett::mpo::{tensor4_from_data, MPO}; -use tensor4all_simplett::{AbstractTensorTrain, SimpleTensorTrain, Tensor3Ops}; +use tensor4all_simplett::{ + AbstractTensorTrain, CompressionMethod, CompressionOptions, SimpleTensorTrain, Tensor3Ops, +}; /// A scalar function of two variables that a benchmark case can compress. /// @@ -22,6 +30,199 @@ pub trait Field2D: Clone + 'static { fn eval(&self, x: f64, y: f64) -> f64; } +/// Construct one axis-aligned Gaussian as a quantics MPO. +/// +/// Each one-variable factor is built by multiscale interpolative QTT with its +/// center marked as the interval that must remain refined. The factors are then +/// combined core-wise, without another interpolation or approximation. +/// +/// This construction assumes `r` is large enough to resolve both Gaussian +/// widths on `[-box_l, box_l)`. It does not attempt to detect an under-resolved +/// grid. +pub fn axis_aligned_gaussian_mpo( + r: usize, + box_l: f64, + center: (f64, f64), + alpha: (f64, f64), + weight: f64, + polynomial_degree: usize, + tolerance: f64, +) -> anyhow::Result> { + anyhow::ensure!(r >= 2, "an interpolative QTT needs at least two bits"); + anyhow::ensure!(box_l.is_finite() && box_l > 0.0, "invalid box size"); + anyhow::ensure!( + center.0.abs() < box_l && center.1.abs() < box_l, + "Gaussian center is outside the box" + ); + anyhow::ensure!( + alpha.0.is_finite() && alpha.0 > 0.0 && alpha.1.is_finite() && alpha.1 > 0.0, + "invalid Gaussian width" + ); + anyhow::ensure!(weight.is_finite(), "invalid Gaussian weight"); + anyhow::ensure!( + tolerance.is_finite() && tolerance >= 0.0, + "invalid interpolation tolerance" + ); + + let options = InterpolativeQttOptions::default().with_tolerance(tolerance); + let x = interpolate_multi_scale( + |value| weight * (-alpha.0 * (value - center.0).powi(2)).exp(), + -box_l, + box_l, + r, + polynomial_degree, + &[center.0], + &options, + )?; + let y = interpolate_multi_scale( + |value| (-alpha.1 * (value - center.1).powi(2)).exp(), + -box_l, + box_l, + r, + polynomial_degree, + &[center.1], + &options, + )?; + + let tensors = x + .site_tensors() + .iter() + .zip(y.site_tensors()) + .map(|(x_core, y_core)| { + let fused = direct_product_core_tensors(&[x_core.clone(), y_core.clone()])?; + let mut data = + Vec::with_capacity(fused.left_dim() * fused.site_dim() * fused.right_dim()); + for right in 0..fused.right_dim() { + for y_bit in 0..2 { + for x_bit in 0..2 { + for left in 0..fused.left_dim() { + data.push(*fused.get3(left, x_bit + 2 * y_bit, right)); + } + } + } + } + tensor4_from_data(data, fused.left_dim(), 2, 2, fused.right_dim()) + .map_err(anyhow::Error::from) + }) + .collect::>>()?; + MPO::new(tensors).map_err(anyhow::Error::from) +} + +fn rotated_gaussian_qtt( + r: usize, + box_l: f64, + center: (f64, f64), + quad: (f64, f64, f64), + weight: f64, + polynomial_degree: usize, + tolerance: f64, +) -> anyhow::Result> { + anyhow::ensure!(r >= 2, "an interpolative QTT needs at least two bits"); + anyhow::ensure!(box_l.is_finite() && box_l > 0.0, "invalid box size"); + anyhow::ensure!( + center.0.abs() < box_l && center.1.abs() < box_l, + "Gaussian center is outside the box" + ); + let (a, b, c) = quad; + let discriminant = ((a - c).powi(2) + 4.0 * b * b).sqrt(); + let lambda_major = 0.5 * (a + c - discriminant); + let lambda_minor = 0.5 * (a + c + discriminant); + anyhow::ensure!( + lambda_major.is_finite() && lambda_major > 0.0 && lambda_minor.is_finite(), + "Gaussian quadratic form is not positive definite" + ); + anyhow::ensure!(weight.is_finite(), "invalid Gaussian weight"); + anyhow::ensure!( + tolerance.is_finite() && tolerance >= 0.0, + "invalid interpolation tolerance" + ); + + let mut major = if b.abs() > f64::EPSILON * a.abs().max(c.abs()) { + (b, lambda_major - a) + } else if a <= c { + (1.0, 0.0) + } else { + (0.0, 1.0) + }; + let norm = major.0.hypot(major.1); + major.0 /= norm; + major.1 /= norm; + let sigma_major = (2.0 * lambda_major).sqrt().recip(); + let sigma_minor = (2.0 * lambda_minor).sqrt().recip(); + let unsafe_spacing = 0.5 * sigma_minor; + let n = (4.0 * sigma_major / unsafe_spacing).ceil() as isize; + let cusp_locations = (-n..=n) + .filter_map(|i| { + let distance = i as f64 * unsafe_spacing; + let point = [center.0 + distance * major.0, center.1 + distance * major.1]; + (point[0].abs() < box_l && point[1].abs() < box_l).then(|| point.to_vec()) + }) + .collect::>(); + let options = InterpolativeQttOptions::default().with_tolerance(tolerance); + interpolate_multi_scale_nd( + |xy| { + let (dx, dy) = (xy[0] - center.0, xy[1] - center.1); + weight * (-(a * dx * dx + 2.0 * b * dx * dy + c * dy * dy)).exp() + }, + &[-box_l, -box_l], + &[box_l, box_l], + r, + polynomial_degree, + &cusp_locations, + &options, + ) + .map_err(anyhow::Error::from) +} + +fn fused_qtt_to_mpo(qtt: &SimpleTensorTrain) -> anyhow::Result> { + let tensors = qtt + .site_tensors() + .iter() + .map(|core| { + anyhow::ensure!(core.site_dim() == 4, "expected fused site dimension 4"); + let mut data = Vec::with_capacity(core.left_dim() * 4 * core.right_dim()); + for right in 0..core.right_dim() { + for y_bit in 0..2 { + for x_bit in 0..2 { + for left in 0..core.left_dim() { + data.push(*core.get3(left, x_bit + 2 * y_bit, right)); + } + } + } + } + tensor4_from_data(data, core.left_dim(), 2, 2, core.right_dim()) + .map_err(anyhow::Error::from) + }) + .collect::>>()?; + MPO::new(tensors).map_err(anyhow::Error::from) +} + +/// Construct one rotated anisotropic Gaussian as a multiscale quantics MPO. +/// +/// Unsafe interpolation points are spaced by half the minor-axis standard +/// deviation along the major axis through the center. This keeps every box +/// intersecting the narrow ridge refined without making the entire 2D box unsafe. +/// `r` must be large enough to resolve the minor width. +pub fn rotated_gaussian_mpo( + r: usize, + box_l: f64, + center: (f64, f64), + quad: (f64, f64, f64), + weight: f64, + polynomial_degree: usize, + tolerance: f64, +) -> anyhow::Result> { + fused_qtt_to_mpo(&rotated_gaussian_qtt( + r, + box_l, + center, + quad, + weight, + polynomial_degree, + tolerance, + )?) +} + /// A sum of isotropic 2D Gaussians `w_i exp(-a_i ((x-cx_i)^2 + (y-cy_i)^2))`. #[derive(Clone, Debug)] pub struct GaussianMixture2D { @@ -105,6 +306,58 @@ pub struct AnisoMixture2D { } impl AnisoMixture2D { + /// Construct this mixture by adding individually interpolated Gaussian QTTs. + /// + /// The sum is SVD-compressed after every addition with the supplied relative + /// per-bond tolerance; the benchmark records this separately from its final + /// global relative-L2 truncation. + pub fn to_multiscale_mpo( + &self, + r: usize, + box_l: f64, + polynomial_degree: usize, + interpolation_tolerance: f64, + addition_tolerance: f64, + ) -> anyhow::Result> { + anyhow::ensure!( + self.weights.len() == self.quad.len() && self.weights.len() == self.centers.len(), + "anisotropic mixture arrays have different lengths" + ); + anyhow::ensure!(!self.weights.is_empty(), "anisotropic mixture is empty"); + anyhow::ensure!( + addition_tolerance.is_finite() && addition_tolerance >= 0.0, + "invalid addition tolerance" + ); + let compression = CompressionOptions { + method: CompressionMethod::SVD, + tolerance: addition_tolerance, + max_bond_dim: None, + normalize_error: true, + }; + let mut sum = rotated_gaussian_qtt( + r, + box_l, + self.centers[0], + self.quad[0], + self.weights[0], + polynomial_degree, + interpolation_tolerance, + )?; + for i in 1..self.weights.len() { + sum = sum.add(&rotated_gaussian_qtt( + r, + box_l, + self.centers[i], + self.quad[i], + self.weights[i], + polynomial_degree, + interpolation_tolerance, + )?)?; + sum.compress(&compression)?; + } + fused_qtt_to_mpo(&sum) + } + /// Draw `n` anisotropic spikes on the box `[-box_l, box_l)^2`. /// /// Weights are `U[0.5, 1.5]`, the aspect ratio `rho` is log-uniform in @@ -182,6 +435,124 @@ impl Field2D for AnisoMixture2D { } } +/// Spatially indexed evaluator with a rigorous global absolute tail bound. +/// +/// Components whose center is farther than their major-axis cutoff radius are +/// omitted. The common exponent threshold is chosen so the sum of all omitted +/// positive Gaussian tails is at most `absolute_tolerance` at every point. +#[derive(Clone, Debug)] +pub struct LocalizedAnisoField { + mixture: AnisoMixture2D, + absolute_tolerance: f64, + bin_width: f64, + cutoff_squared: Vec, + bins: HashMap<(i64, i64), Vec>, +} + +impl LocalizedAnisoField { + /// Build a spatial index without changing the underlying random mixture. + pub fn new(mixture: AnisoMixture2D, absolute_tolerance: f64) -> anyhow::Result { + anyhow::ensure!( + absolute_tolerance.is_finite() && absolute_tolerance > 0.0, + "localized evaluator tolerance must be positive and finite" + ); + anyhow::ensure!( + mixture.weights.len() == mixture.quad.len() + && mixture.weights.len() == mixture.centers.len() + && !mixture.weights.is_empty(), + "invalid anisotropic mixture" + ); + anyhow::ensure!( + mixture + .weights + .iter() + .all(|weight| weight.is_finite() && *weight >= 0.0), + "localized evaluator requires finite non-negative weights" + ); + anyhow::ensure!( + mixture + .centers + .iter() + .all(|(x, y)| x.is_finite() && y.is_finite()), + "localized evaluator requires finite centers" + ); + let total_weight = mixture.weights.iter().sum::(); + anyhow::ensure!( + total_weight.is_finite() && total_weight > absolute_tolerance, + "localized evaluator tolerance must be smaller than total weight" + ); + let exponent_cutoff = (total_weight / absolute_tolerance).ln(); + let cutoff_squared = mixture + .quad + .iter() + .map(|&(a, b, c)| { + let discriminant = ((a - c).powi(2) + 4.0 * b * b).sqrt(); + let lambda_min = 0.5 * (a + c - discriminant); + anyhow::ensure!(lambda_min.is_finite() && lambda_min > 0.0); + Ok(exponent_cutoff / lambda_min) + }) + .collect::>>()?; + let bin_width = cutoff_squared + .iter() + .copied() + .fold(0.0_f64, f64::max) + .sqrt(); + let mut bins = HashMap::<(i64, i64), Vec>::new(); + for (i, &(x, y)) in mixture.centers.iter().enumerate() { + bins.entry(Self::bin_key(x, y, bin_width)) + .or_default() + .push(i); + } + Ok(Self { + mixture, + absolute_tolerance, + bin_width, + cutoff_squared, + bins, + }) + } + + fn bin_key(x: f64, y: f64, bin_width: f64) -> (i64, i64) { + ( + (x / bin_width).floor() as i64, + (y / bin_width).floor() as i64, + ) + } + + /// Guaranteed pointwise absolute error bound relative to exact evaluation. + pub fn absolute_tolerance(&self) -> f64 { + self.absolute_tolerance + } + + /// Evaluate using only components whose tail can exceed the global budget. + pub fn eval(&self, x: f64, y: f64) -> f64 { + let (bx, by) = Self::bin_key(x, y, self.bin_width); + let mut value = 0.0; + for dx_bin in -1..=1 { + for dy_bin in -1..=1 { + if let Some(indices) = self.bins.get(&(bx + dx_bin, by + dy_bin)) { + for &i in indices { + let (cx, cy) = self.mixture.centers[i]; + let (dx, dy) = (x - cx, y - cy); + if dx * dx + dy * dy <= self.cutoff_squared[i] { + let (a, b, c) = self.mixture.quad[i]; + value += self.mixture.weights[i] + * (-(a * dx * dx + 2.0 * b * dx * dy + c * dy * dy)).exp(); + } + } + } + } + } + value + } +} + +impl Field2D for LocalizedAnisoField { + fn eval(&self, x: f64, y: f64) -> f64 { + LocalizedAnisoField::eval(self, x, y) + } +} + /// Closed form of `int f(x,y) g(y,z) dy` over the whole real line. pub fn analytic_contraction(f: &GaussianMixture2D, g: &GaussianMixture2D, x: f64, z: f64) -> f64 { let mut s = 0.0; @@ -307,6 +678,17 @@ pub fn to_quantics_fused_tt_field( box_l: f64, tol: f64, max_bond: usize, +) -> anyhow::Result<(SimpleTensorTrain, f64)> { + to_quantics_fused_tt_field_with_pivots(mix, r, box_l, tol, max_bond, None) +} + +fn to_quantics_fused_tt_field_with_pivots( + mix: &M, + r: usize, + box_l: f64, + tol: f64, + max_bond: usize, + initial_pivots: Option>>, ) -> anyhow::Result<(SimpleTensorTrain, f64)> { let grid = DiscretizedGrid::builder(&[r, r]) .with_lower_bound(&[-box_l, -box_l]) @@ -314,12 +696,19 @@ pub fn to_quantics_fused_tt_field( .with_unfolding_scheme(UnfoldingScheme::Fused) .build()?; let m = mix.clone(); - let opts = QtciOptions::default() + let mut opts = QtciOptions::default() .with_tolerance(tol) .with_max_bond_dim(max_bond) .with_unfoldingscheme(UnfoldingScheme::Fused); - let (qtci, _ranks, _errs) = - quanticscrossinterpolate(&grid, move |xy: &[f64]| m.eval(xy[0], xy[1]), None, opts)?; + if initial_pivots.is_some() { + opts = opts.with_nrandominitpivot(0); + } + let (qtci, _ranks, _errs) = quanticscrossinterpolate( + &grid, + move |xy: &[f64]| m.eval(xy[0], xy[1]), + initial_pivots, + opts, + )?; let step = 2.0 * box_l / (1u64 << r) as f64; Ok((qtci.tensor_train(), step)) } @@ -345,35 +734,31 @@ pub fn to_quantics_mpo_field( max_bond: usize, ) -> anyhow::Result<(MPO, f64)> { let (tt, step) = to_quantics_fused_tt_field(mix, r, box_l, tol, max_bond)?; - let mut cores4 = Vec::with_capacity(tt.len()); - for c in tt.site_tensors() { - let (l, s, rd) = (c.left_dim(), c.site_dim(), c.right_dim()); - anyhow::ensure!(s == 4, "expected fused site dim 4, got {s}"); - // Fused local index is `s = s1 + 2*s2` with s1 the bit of variable 1 (x) - // as the least significant digit. Verified both against the quanticsgrids - // source (add_fused_indices pushes variables in reverse order, so the last - // variable lands at pos_in_site 0 and gets place value 2^(len-1) = 2) and - // empirically by `quantics_mpo_evaluates_to_function_values`. - // Task 8's Julia mirror must use the same order. - let mut data = vec![0.0f64; l * 4 * rd]; - for rr in 0..rd { - for s2 in 0..2 { - for s1 in 0..2 { - for ll in 0..l { - let fused = s1 + 2 * s2; - data[ll + l * (s1 + 2 * (s2 + 2 * rr))] = *c.get3(ll, fused, rr); - } - } - } - } - cores4.push(tensor4_from_data(data, l, 2, 2, rd)?); - } - Ok((MPO::new(cores4)?, step)) + Ok((fused_qtt_to_mpo(&tt)?, step)) +} + +/// Deterministic [`to_quantics_mpo_field`] variant using explicit zero-based grid pivots. +/// +/// Random initial pivots are disabled when this entry point is used, so identical +/// inputs and pivots produce cacheable input tensors. +pub fn to_quantics_mpo_field_with_pivots( + mix: &M, + r: usize, + box_l: f64, + tol: f64, + max_bond: usize, + initial_pivots: Vec>, +) -> anyhow::Result<(MPO, f64)> { + anyhow::ensure!(!initial_pivots.is_empty(), "initial pivot list is empty"); + let (tt, step) = + to_quantics_fused_tt_field_with_pivots(mix, r, box_l, tol, max_bond, Some(initial_pivots))?; + Ok((fused_qtt_to_mpo(&tt)?, step)) } #[cfg(test)] mod tests { use super::*; + use crate::harness::index_to_bits; fn quadrature(f: impl Fn(f64) -> f64, lo: f64, hi: f64, intervals: usize) -> f64 { let step = (hi - lo) / intervals as f64; @@ -385,6 +770,198 @@ mod tests { .sum() } + #[test] + fn multiscale_axis_aligned_gaussian_mpo_matches_grid_values() { + let (r, box_l) = (8, 1.0); + let center = (0.17, -0.23); + let alpha = (8.0, 20.0); + let weight = 1.3; + let mpo = axis_aligned_gaussian_mpo(r, box_l, center, alpha, weight, 12, 1e-12).unwrap(); + let step = 2.0 * box_l / (1usize << r) as f64; + for (ix, iy) in [(0, 0), (23, 197), (128, 128), (150, 98), (255, 255)] { + let x = -box_l + ix as f64 * step; + let y = -box_l + iy as f64 * step; + let indices = index_to_bits(ix, r) + .into_iter() + .zip(index_to_bits(iy, r)) + .flat_map(|(x_bit, y_bit)| [x_bit, y_bit]) + .collect::>(); + let got = mpo.evaluate(&indices).unwrap(); + let expected = weight + * (-alpha.0 * (x - center.0).powi(2) - alpha.1 * (y - center.1).powi(2)).exp(); + assert!( + (got - expected).abs() < 2e-8, + "({x}, {y}): {got} != {expected}" + ); + } + } + + #[test] + fn multiscale_gaussian_mpo_stays_accurate_as_r_increases() { + let box_l = 1.0; + let center = (0.173, -0.227); + let alpha = (32.0, 200.0); + let weight = 1.3; + let sigma_min = 1.0 / (2.0_f64 * alpha.1).sqrt(); + + for r in [8, 10, 12] { + let step = 2.0 * box_l / (1usize << r) as f64; + assert!( + step <= sigma_min / 4.0, + "R={r} does not resolve the narrow width" + ); + let mpo = + axis_aligned_gaussian_mpo(r, box_l, center, alpha, weight, 24, 1e-12).unwrap(); + let grid_size = 1usize << r; + let mut max_error = 0.0_f64; + for sample in 0..128 { + let ix = (73 * sample + 11) % grid_size; + let iy = (151 * sample + 29) % grid_size; + let x = -box_l + ix as f64 * step; + let y = -box_l + iy as f64 * step; + let indices = index_to_bits(ix as u64, r) + .into_iter() + .zip(index_to_bits(iy as u64, r)) + .flat_map(|(x_bit, y_bit)| [x_bit, y_bit]) + .collect::>(); + let got = mpo.evaluate(&indices).unwrap(); + let expected = weight + * (-alpha.0 * (x - center.0).powi(2) - alpha.1 * (y - center.1).powi(2)).exp(); + max_error = max_error.max((got - expected).abs()); + } + eprintln!("R={r} rank={} max_error={max_error:.3e}", mpo.rank()); + assert!( + max_error < 5e-8, + "R={r}, rank={}, max error={max_error:.3e}", + mpo.rank() + ); + assert!( + mpo.rank() <= 400, + "R={r} unexpectedly increased rank to {}", + mpo.rank() + ); + } + } + + #[test] + fn multiscale_rotated_gaussian_mpo_resolves_a_narrow_ridge() { + let (r, box_l) = (10, 1.0); + let center = (0.07, -0.09); + let (sigma_minor, rho, angle) = (0.05_f64, 8.0_f64, std::f64::consts::FRAC_PI_4); + let lambda_minor = 1.0 / (2.0 * sigma_minor.powi(2)); + let lambda_major = lambda_minor / rho.powi(2); + let (sin, cos) = angle.sin_cos(); + let quad = ( + lambda_major * cos * cos + lambda_minor * sin * sin, + (lambda_major - lambda_minor) * sin * cos, + lambda_major * sin * sin + lambda_minor * cos * cos, + ); + let mpo = rotated_gaussian_mpo(r, box_l, center, quad, 1.3, 28, 1e-12).unwrap(); + let step = 2.0 * box_l / (1usize << r) as f64; + let mut max_error = 0.0_f64; + { + let mut check_point = |ix: usize, iy: usize| { + let x = -box_l + ix as f64 * step; + let y = -box_l + iy as f64 * step; + let indices = index_to_bits(ix as u64, r) + .into_iter() + .zip(index_to_bits(iy as u64, r)) + .flat_map(|(x_bit, y_bit)| [x_bit, y_bit]) + .collect::>(); + let got = mpo.evaluate(&indices).unwrap(); + let (dx, dy) = (x - center.0, y - center.1); + let expected = + 1.3 * (-(quad.0 * dx * dx + 2.0 * quad.1 * dx * dy + quad.2 * dy * dy)).exp(); + max_error = max_error.max((got - expected).abs()); + }; + for sample in 0..256 { + check_point( + (73 * sample + 11) % (1usize << r), + (151 * sample + 29) % (1usize << r), + ); + } + for major_step in -6..=6 { + for minor_step in -4..=4 { + let u = major_step as f64 * 0.5 * sigma_minor * rho; + let v = minor_step as f64 * 0.5 * sigma_minor; + let x = center.0 + cos * u - sin * v; + let y = center.1 + sin * u + cos * v; + let ix = ((x + box_l) / step).round() as usize; + let iy = ((y + box_l) / step).round() as usize; + check_point(ix.min((1usize << r) - 1), iy.min((1usize << r) - 1)); + } + } + } + eprintln!("rotated rank={} max_error={max_error:.3e}", mpo.rank()); + assert!( + max_error < 5e-8, + "rank={}, max error={max_error:.3e}", + mpo.rank() + ); + } + + #[test] + fn multiscale_anisotropic_mixture_matches_grid_values() { + let (r, box_l) = (8, 1.0); + let mixture = AnisoMixture2D::random(2, 0.7, 0.08, 3.0, 17); + let mpo = mixture + .to_multiscale_mpo(r, box_l, 24, 1e-12, 1e-10) + .unwrap(); + let step = 2.0 * box_l / (1usize << r) as f64; + let mut max_error = 0.0_f64; + for sample in 0..128 { + let ix = (73 * sample + 11) % (1usize << r); + let iy = (151 * sample + 29) % (1usize << r); + let x = -box_l + ix as f64 * step; + let y = -box_l + iy as f64 * step; + let indices = index_to_bits(ix as u64, r) + .into_iter() + .zip(index_to_bits(iy as u64, r)) + .flat_map(|(x_bit, y_bit)| [x_bit, y_bit]) + .collect::>(); + max_error = max_error.max((mpo.evaluate(&indices).unwrap() - mixture.eval(x, y)).abs()); + } + assert!( + max_error < 1e-7, + "rank={}, max error={max_error:.3e}", + mpo.rank() + ); + } + + #[test] + fn localized_anisotropic_evaluator_respects_absolute_tail_bound() { + let mixture = AnisoMixture2D::random(512, 4.0, 0.05, 8.0, 23); + let localized = LocalizedAnisoField::new(mixture.clone(), 1e-10).unwrap(); + assert_eq!(localized.absolute_tolerance(), 1e-10); + for sample in 0..256 { + let x = -5.0 + 10.0 * ((73 * sample + 11) % 257) as f64 / 256.0; + let y = -5.0 + 10.0 * ((151 * sample + 29) % 257) as f64 / 256.0; + let error = (localized.eval(x, y) - mixture.eval(x, y)).abs(); + assert!(error <= 1e-10, "({x}, {y}): error={error:.3e}"); + } + } + + #[test] + fn explicit_tci_pivots_make_input_construction_deterministic() { + let mixture = AnisoMixture2D::random(4, 1.0, 0.1, 3.0, 29); + let pivots = vec![vec![16, 16], vec![8, 24], vec![24, 8]]; + let (first, _) = + to_quantics_mpo_field_with_pivots(&mixture, 5, 1.0, 1e-8, 64, pivots.clone()).unwrap(); + let (second, _) = + to_quantics_mpo_field_with_pivots(&mixture, 5, 1.0, 1e-8, 64, pivots).unwrap(); + assert_eq!(first.rank(), second.rank()); + for (left, right) in first.site_tensors().iter().zip(second.site_tensors()) { + assert_eq!(left.to_col_major_vec(), right.to_col_major_vec()); + } + } + + #[test] + fn localized_anisotropic_evaluator_rejects_negative_weights() { + let mut mixture = AnisoMixture2D::random(1, 1.0, 0.05, 8.0, 23); + mixture.weights[0] = -1.0; + assert!(LocalizedAnisoField::new(mixture, 1e-10).is_err()); + } + #[test] fn analytic_contraction_matches_quadrature() { let f = GaussianMixture2D::random(3, 4.0, (0.5, 2.0), 1); diff --git a/src/patched_mpo.rs b/src/patched_mpo.rs index fc9cb99..a6f8fa9 100644 --- a/src/patched_mpo.rs +++ b/src/patched_mpo.rs @@ -1,7 +1,7 @@ //! Patched MPO-MPO contraction on the legacy TT and chain TreeTN wrappers. use tensor4all_core::{DynIndex, SvdTruncationPolicy}; -use tensor4all_itensorlike::{ContractOptions, TensorTrain}; +use tensor4all_itensorlike::{ContractOptions, TensorTrain, TruncateOptions}; use tensor4all_partitionedtreetn as tree; use tensor4all_partitionedtt as tt; use tensor4all_simplett::mpo::MPO; @@ -42,6 +42,62 @@ pub struct PatchedMpoOutput { pub max_patch_bond: usize, } +/// Convert a compatible MPO pair to tensor trains with one shared site index. +pub fn mpo_pair_to_tensortrains( + left: &MPO, + right: &MPO, +) -> anyhow::Result<(TensorTrain, TensorTrain)> { + anyhow::ensure!(!left.is_empty(), "patched MPO input is empty"); + anyhow::ensure!(left.len() == right.len(), "MPO lengths differ"); + + let n = left.len(); + let x: Vec<_> = (0..n) + .map(|site| DynIndex::new_dyn(left.site_dim(site).0)) + .collect(); + let y: Vec<_> = (0..n) + .map(|site| DynIndex::new_dyn(left.site_dim(site).1)) + .collect(); + let z: Vec<_> = (0..n) + .map(|site| DynIndex::new_dyn(right.site_dim(site).1)) + .collect(); + for site in 0..n { + anyhow::ensure!( + left.site_dim(site).1 == right.site_dim(site).0, + "contracted site dimensions differ at site {site}" + ); + } + + Ok(( + mpo_to_tensortrain(left, &x, &y)?, + mpo_to_tensortrain(right, &y, &z)?, + )) +} + +/// Truncate a tensor train to a relative L2 error tolerance with no rank cap. +pub fn truncate_tensortrain_l2(tt: &mut TensorTrain, rtol: f64) -> anyhow::Result<()> { + anyhow::ensure!(rtol.is_finite() && rtol >= 0.0, "invalid L2 tolerance"); + let policy = SvdTruncationPolicy::new(rtol * rtol) + .with_relative() + .with_squared_values() + .with_discarded_tail_sum(); + tt.truncate(&TruncateOptions::svd().with_svd_policy(policy))?; + Ok(()) +} + +/// Number of stored scalar parameters in a tensor train. +pub fn tensortrain_n_params(tt: &TensorTrain) -> usize { + tt.tensors() + .into_iter() + .map(|tensor| { + tensor + .indices() + .iter() + .map(|index| index.dim) + .product::() + }) + .sum() +} + impl PatchedMpoPair { /// Build both patched representations from exactly the same MPO cores and indices. pub fn new( @@ -50,29 +106,7 @@ impl PatchedMpoPair { rtol: f64, patch_max_bond: usize, ) -> anyhow::Result { - anyhow::ensure!(!left.is_empty(), "patched MPO input is empty"); - anyhow::ensure!(left.len() == right.len(), "MPO lengths differ"); - anyhow::ensure!(patch_max_bond > 0, "patch max bond must be positive"); - - let n = left.len(); - let x: Vec<_> = (0..n) - .map(|site| DynIndex::new_dyn(left.site_dim(site).0)) - .collect(); - let y: Vec<_> = (0..n) - .map(|site| DynIndex::new_dyn(left.site_dim(site).1)) - .collect(); - let z: Vec<_> = (0..n) - .map(|site| DynIndex::new_dyn(right.site_dim(site).1)) - .collect(); - for site in 0..n { - anyhow::ensure!( - left.site_dim(site).1 == right.site_dim(site).0, - "contracted site dimensions differ at site {site}" - ); - } - - let left_train = mpo_to_tensortrain(left, &x, &y)?; - let right_train = mpo_to_tensortrain(right, &y, &z)?; + let (left_train, right_train) = mpo_pair_to_tensortrains(left, right)?; Self::from_tensortrains(left_train, right_train, rtol, patch_max_bond) } @@ -82,6 +116,17 @@ impl PatchedMpoPair { right_train: TensorTrain, rtol: f64, patch_max_bond: usize, + ) -> anyhow::Result { + Self::from_tensortrains_with_input_rtol(left_train, right_train, rtol, rtol, patch_max_bond) + } + + /// Build patched inputs with a tolerance independent of contraction truncation. + pub fn from_tensortrains_with_input_rtol( + left_train: TensorTrain, + right_train: TensorTrain, + input_rtol: f64, + contract_rtol: f64, + patch_max_bond: usize, ) -> anyhow::Result { anyhow::ensure!(!left_train.is_empty(), "cached MPO input is empty"); anyhow::ensure!( @@ -128,50 +173,62 @@ impl PatchedMpoPair { let input_max_bond = left_train.max_bond_dim().max(right_train.max_bond_dim()); let left_order = y.iter().chain(&x).cloned().collect::>(); let right_order = y.iter().chain(&z).cloned().collect::>(); - let tt_options = tt::PatchingOptions { - rtol, + let tt_input_options = tt::PatchingOptions { + rtol: input_rtol, max_bond_dim: Some(patch_max_bond), patch_order: left_order.clone(), split_strategy: tt::PatchSplitStrategy::Sequential, }; - let tt_right_options = tt::PatchingOptions { + let tt_right_input_options = tt::PatchingOptions { patch_order: right_order.clone(), - ..tt_options.clone() + ..tt_input_options.clone() }; let tt_left = tt::add_with_patching( vec![tt::SubDomainTT::from_tt(left_train.clone())], - &tt_options, + &tt_input_options, )?; let tt_right = tt::add_with_patching( vec![tt::SubDomainTT::from_tt(right_train.clone())], - &tt_right_options, + &tt_right_input_options, )?; let global_left = left_train.clone(); let global_right = right_train.clone(); + let tree_input_options = tree::PatchingOptions { + rtol: input_rtol, + max_bond_dim: Some(patch_max_bond), + patch_order: left_order.clone(), + split_strategy: tree::PatchSplitStrategy::Sequential, + }; + let tree_right_input_options = tree::PatchingOptions { + patch_order: right_order.clone(), + ..tree_input_options.clone() + }; + let tt_options = tt::PatchingOptions { + rtol: contract_rtol, + max_bond_dim: Some(patch_max_bond), + patch_order: left_order.clone(), + split_strategy: tt::PatchSplitStrategy::Sequential, + }; let tree_options = tree::PatchingOptions { - rtol, + rtol: contract_rtol, max_bond_dim: Some(patch_max_bond), patch_order: left_order, split_strategy: tree::PatchSplitStrategy::Sequential, }; - let tree_right_options = tree::PatchingOptions { - patch_order: right_order, - ..tree_options.clone() - }; let center = n - 1; let tree_left = tree::add_with_patching( vec![tree::SubDomainTreeTN::from_treetn( left_train.into_treetn(), )?], ¢er, - &tree_options, + &tree_input_options, )?; let tree_right = tree::add_with_patching( vec![tree::SubDomainTreeTN::from_treetn( right_train.into_treetn(), )?], ¢er, - &tree_right_options, + &tree_right_input_options, )?; anyhow::ensure!( (tree_left.len(), tree_right.len()) == (tt_left.len(), tt_right.len()), @@ -224,6 +281,20 @@ impl PatchedMpoPair { ) } + /// Total stored parameters in the left and right patched inputs. + pub fn input_patch_n_params(&self) -> (usize, usize) { + ( + self.tt_left + .values() + .map(|subdomain| tensortrain_n_params(subdomain.data())) + .sum(), + self.tt_right + .values() + .map(|subdomain| tensortrain_n_params(subdomain.data())) + .sum(), + ) + } + /// Run one full-sweep global fit on the prepared unpatched inputs. pub fn contract_fit_global( &self,