diff --git a/.gitignore b/.gitignore
index 31221d3..3b49edf 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,4 +1,5 @@
 /target
+/.cache/
 Cargo.lock.orig
 __pycache__/
 .venv/
diff --git a/AGENTS.md b/AGENTS.md
index 06291b3..6979862 100644
--- a/AGENTS.md
+++ b/AGENTS.md
@@ -36,6 +36,26 @@ times of the memory heavy arms differ by factors between machines (README known
 issue 10). Never put a hostname in a result file, the machine label and the
 hardware fields in `run.yaml` carry the identity.
 
+## Expensive benchmark inputs are cached
+
+Every runner with nontrivial deterministic input construction must cache that
+input and reuse it across algorithm comparisons and reruns. Input generation,
+cache loading, format conversion and patch preparation stay outside timed
+regions.
+
+- Use one shared ignored cache directory, `.cache/inputs/` by default.
+- Include every input-defining parameter, seed and an explicit input schema
+  version in the cache key. Bump the schema when generation semantics change.
+- Write cache entries atomically and treat an existing entry as immutable.
+- Validate loaded shapes, ranks, dtypes and shared index identities before use.
+- Record the cache key, hit or miss, load time and build time in each result.
+- Provide explicit refresh and cache-directory environment knobs. Never silently
+  accept a corrupt or incompatible entry.
+- Cache inputs only, never algorithm outputs or timed intermediates.
+- Tests whose subject is input generation or serialization must force a rebuild;
+  ordinary unit tests remain isolated and must not depend on persistent cache
+  state.
+
 ## Upstream is pinned
 
 `Cargo.toml` pins every `tensor4all-rs` crate to one `rev`. Bump it deliberately,
diff --git a/Cargo.lock b/Cargo.lock
index b368d8c..de9bb9e 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -287,13 +287,10 @@ dependencies = [
  "generativity",
  "libm",
  "nano-gemm",
- "npyz",
  "num-complex",
  "num-traits",
  "private-gemm-x86",
  "pulp",
- "rand 0.9.5",
- "rand_distr",
  "rayon",
  "reborrow",
  "spindle",
@@ -612,26 +609,6 @@ version = "1.0.18"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "8f42a60cbdf9a97f5d2305f08a87dc4e09308d1276d28c869c684d7777685682"
 
-[[package]]
-name = "lapack"
-version = "0.20.0"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "508f7f48955e5e4ec3feb771a4d981c6e59ff2400d750cf41733a3350fdb2dc4"
-dependencies = [
- "lapack-sys",
- "libc",
- "num-complex",
-]
-
-[[package]]
-name = "lapack-sys"
-version = "0.15.0"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "314b879030845b68571809a6978e52d3b67ac5fba07e77b1b317b484092e2fb5"
-dependencies = [
- "libc",
-]
-
 [[package]]
 name = "lazy_static"
 version = "1.5.0"
@@ -807,17 +784,6 @@ dependencies = [
  "rawpointer",
 ]
 
-[[package]]
-name = "npyz"
-version = "0.8.4"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "9f0e759e014e630f90af745101b614f761306ddc541681e546649068e25ec1b9"
-dependencies = [
- "byteorder",
- "num-bigint",
- "py_literal",
-]
-
 [[package]]
 name = "nu-ansi-term"
 version = "0.50.3"
@@ -827,16 +793,6 @@ dependencies = [
  "windows-sys 0.61.2",
 ]
 
-[[package]]
-name = "num-bigint"
-version = "0.4.8"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "c89e69e7e0f03bea5ef08013795c25018e101932225a656383bd384495ecc367"
-dependencies = [
- "num-integer",
- "num-traits",
-]
-
 [[package]]
 name = "num-complex"
 version = "0.4.6"
@@ -845,7 +801,6 @@ checksum = "73f88a1307638156682bada9d7604135552957b7818057dcef22705b4d509495"
 dependencies = [
  "bytemuck",
  "num-traits",
- "rand 0.8.7",
 ]
 
 [[package]]
@@ -884,7 +839,7 @@ source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "cd9b9ad7716ee4ed2a6a5277f4fd1a65f683eda2c05f0f45c8fc6482795a1edc"
 dependencies = [
  "priority-queue",
- "rand 0.9.5",
+ "rand",
  "rayon",
  "serde",
  "serde_json",
@@ -932,48 +887,6 @@ version = "0.2.3"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "2ee67f1008b1ba2321834326597b8e186293b049a023cdef258527550b9935b4"
 
-[[package]]
-name = "pest"
-version = "2.8.8"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "7df728be843c7070fab6ab7c328c4e9e9d78e23bf749c0669c86ee7ebfa050a2"
-dependencies = [
- "memchr",
- "ucd-trie",
-]
-
-[[package]]
-name = "pest_derive"
-version = "2.8.8"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "9e2dd6fc3b26b3462ee188aac870f5a41d398f1cd5e2408d16531bd71c9591fd"
-dependencies = [
- "pest",
- "pest_generator",
-]
-
-[[package]]
-name = "pest_generator"
-version = "2.8.8"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "6a7a9205cfb6f596a9e8b689c0a15f9ceb7a1aafae7aaf788150ac65b29975b6"
-dependencies = [
- "pest",
- "pest_meta",
- "proc-macro2",
- "quote",
- "syn 2.0.119",
-]
-
-[[package]]
-name = "pest_meta"
-version = "2.8.8"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "85abd351c0de1e8384fc791a0737111a350394937e92b956b743dac12429f57c"
-dependencies = [
- "pest",
-]
-
 [[package]]
 name = "petgraph"
 version = "0.7.1"
@@ -1088,19 +1001,6 @@ version = "0.1.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "1d8f70e07b9c3962945a74e59ca1c511bba65b6419468acc217c457d93f3c740"
 
-[[package]]
-name = "py_literal"
-version = "0.4.0"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "102df7a3d46db9d3891f178dcc826dc270a6746277a9ae6436f8d29fd490a8e1"
-dependencies = [
- "num-bigint",
- "num-complex",
- "num-traits",
- "pest",
- "pest_derive",
-]
-
 [[package]]
 name = "qd"
 version = "0.8.0"
@@ -1136,15 +1036,6 @@ version = "5.3.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "69cdb34c158ceb288df11e18b4bd39de994f6657d83847bdffdbd7f346754b0f"
 
-[[package]]
-name = "rand"
-version = "0.8.7"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "22f6172bdec972074665ed81ed53b71da00bfc44b65a753cfde883ec4c702a1a"
-dependencies = [
- "rand_core 0.6.4",
-]
-
 [[package]]
 name = "rand"
 version = "0.9.5"
@@ -1152,7 +1043,7 @@ source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "b9ef1d0d795eb7d84685bca4f72f3649f064e6641543d3a8c415898726a57b41"
 dependencies = [
  "rand_chacha",
- "rand_core 0.9.5",
+ "rand_core",
 ]
 
 [[package]]
@@ -1162,15 +1053,9 @@ source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "d3022b5f1df60f26e1ffddd6c66e8aa15de382ae63b3a0c1bfc0e4d3e3f325cb"
 dependencies = [
  "ppv-lite86",
- "rand_core 0.9.5",
+ "rand_core",
 ]
 
-[[package]]
-name = "rand_core"
-version = "0.6.4"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "ec0be4795e2f6a28069bec0b5ff3e2ac9bafc99e6a9a7dc3547996c5c816922c"
-
 [[package]]
 name = "rand_core"
 version = "0.9.5"
@@ -1187,7 +1072,7 @@ source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "6a8615d50dcf34fa31f7ab52692afec947c4dd0ab803cc87cb3b0b4570ff7463"
 dependencies = [
  "num-traits",
- "rand 0.9.5",
+ "rand",
 ]
 
 [[package]]
@@ -1381,28 +1266,32 @@ dependencies = [
 
 [[package]]
 name = "strided-kernel"
-version = "0.1.0"
-source = "git+https://github.com/tensor4all/strided-rs?rev=ea37986f4d9cb99cfc1f62a1d4aea561cb3a9722#ea37986f4d9cb99cfc1f62a1d4aea561cb3a9722"
+version = "0.4.0"
+source = "git+https://github.com/tensor4all/strided-rs.git?rev=b29e7601ba090aa5eafc65b9bde5d9450282e0d8#b29e7601ba090aa5eafc65b9bde5d9450282e0d8"
 dependencies = [
  "num-complex",
  "num-traits",
  "pulp",
+ "rayon",
+ "smallvec",
  "strided-perm",
  "strided-view",
 ]
 
 [[package]]
 name = "strided-perm"
-version = "0.1.0"
-source = "git+https://github.com/tensor4all/strided-rs?rev=ea37986f4d9cb99cfc1f62a1d4aea561cb3a9722#ea37986f4d9cb99cfc1f62a1d4aea561cb3a9722"
+version = "0.4.0"
+source = "git+https://github.com/tensor4all/strided-rs.git?rev=b29e7601ba090aa5eafc65b9bde5d9450282e0d8#b29e7601ba090aa5eafc65b9bde5d9450282e0d8"
 dependencies = [
+ "rayon",
  "strided-view",
+ "thiserror 1.0.69",
 ]
 
 [[package]]
 name = "strided-traits"
-version = "0.1.0"
-source = "git+https://github.com/tensor4all/strided-rs?rev=ea37986f4d9cb99cfc1f62a1d4aea561cb3a9722#ea37986f4d9cb99cfc1f62a1d4aea561cb3a9722"
+version = "0.4.0"
+source = "git+https://github.com/tensor4all/strided-rs.git?rev=b29e7601ba090aa5eafc65b9bde5d9450282e0d8#b29e7601ba090aa5eafc65b9bde5d9450282e0d8"
 dependencies = [
  "num-complex",
  "num-traits",
@@ -1410,8 +1299,8 @@ dependencies = [
 
 [[package]]
 name = "strided-view"
-version = "0.1.0"
-source = "git+https://github.com/tensor4all/strided-rs?rev=ea37986f4d9cb99cfc1f62a1d4aea561cb3a9722#ea37986f4d9cb99cfc1f62a1d4aea561cb3a9722"
+version = "0.4.0"
+source = "git+https://github.com/tensor4all/strided-rs.git?rev=b29e7601ba090aa5eafc65b9bde5d9450282e0d8#b29e7601ba090aa5eafc65b9bde5d9450282e0d8"
 dependencies = [
  "num-complex",
  "num-traits",
@@ -1471,8 +1360,9 @@ name = "t4a-bench"
 version = "0.1.0"
 dependencies = [
  "anyhow",
+ "libm",
  "num-complex",
- "rand 0.9.5",
+ "rand",
  "rand_chacha",
  "serde",
  "serde_json",
@@ -1480,6 +1370,7 @@ dependencies = [
  "tensor4all-core",
  "tensor4all-hdf5",
  "tensor4all-itensorlike",
+ "tensor4all-partitionedtreetn",
  "tensor4all-partitionedtt",
  "tensor4all-quanticstci",
  "tensor4all-simplett",
@@ -1488,8 +1379,8 @@ dependencies = [
 
 [[package]]
 name = "tenferro-ad"
-version = "0.1.0"
-source = "git+https://github.com/tensor4all/tenferro-rs.git?rev=f68b33d9e06cc0af8d2b9c40f7c52d8972cd1266#f68b33d9e06cc0af8d2b9c40f7c52d8972cd1266"
+version = "0.3.0"
+source = "git+https://github.com/tensor4all/tenferro-rs.git?rev=b5a106be3133979d78832a0ca3f4d6b57613b3d7#b5a106be3133979d78832a0ca3f4d6b57613b3d7"
 dependencies = [
  "computegraph",
  "lru",
@@ -1500,36 +1391,41 @@ dependencies = [
  "tenferro-runtime",
  "tenferro-tensor",
  "thiserror 2.0.20",
- "tidu",
 ]
 
 [[package]]
 name = "tenferro-core-ops"
-version = "0.1.0"
-source = "git+https://github.com/tensor4all/tenferro-rs.git?rev=f68b33d9e06cc0af8d2b9c40f7c52d8972cd1266#f68b33d9e06cc0af8d2b9c40f7c52d8972cd1266"
+version = "0.3.0"
+source = "git+https://github.com/tensor4all/tenferro-rs.git?rev=b5a106be3133979d78832a0ca3f4d6b57613b3d7#b5a106be3133979d78832a0ca3f4d6b57613b3d7"
 
 [[package]]
 name = "tenferro-cpu"
-version = "0.1.0"
-source = "git+https://github.com/tensor4all/tenferro-rs.git?rev=f68b33d9e06cc0af8d2b9c40f7c52d8972cd1266#f68b33d9e06cc0af8d2b9c40f7c52d8972cd1266"
+version = "0.3.0"
+source = "git+https://github.com/tensor4all/tenferro-rs.git?rev=b5a106be3133979d78832a0ca3f4d6b57613b3d7#b5a106be3133979d78832a0ca3f4d6b57613b3d7"
 dependencies = [
  "faer",
+ "lru",
  "num-complex",
  "num-traits",
+ "rayon",
  "smallvec",
  "strided-kernel",
+ "tenferro-core-ops",
+ "tenferro-internal-cpu-kernels",
+ "tenferro-runtime",
  "tenferro-tensor",
  "thiserror 2.0.20",
 ]
 
 [[package]]
 name = "tenferro-einsum"
-version = "0.1.0"
-source = "git+https://github.com/tensor4all/tenferro-rs.git?rev=f68b33d9e06cc0af8d2b9c40f7c52d8972cd1266#f68b33d9e06cc0af8d2b9c40f7c52d8972cd1266"
+version = "0.3.0"
+source = "git+https://github.com/tensor4all/tenferro-rs.git?rev=b5a106be3133979d78832a0ca3f4d6b57613b3d7#b5a106be3133979d78832a0ca3f4d6b57613b3d7"
 dependencies = [
  "computegraph",
  "lru",
  "omeco",
+ "smallvec",
  "tenferro-ad",
  "tenferro-cpu",
  "tenferro-internal-extension-macros",
@@ -1537,13 +1433,24 @@ dependencies = [
  "tenferro-runtime",
  "tenferro-tensor",
  "thiserror 2.0.20",
- "tidu",
+]
+
+[[package]]
+name = "tenferro-internal-cpu-kernels"
+version = "0.3.0"
+source = "git+https://github.com/tensor4all/tenferro-rs.git?rev=b5a106be3133979d78832a0ca3f4d6b57613b3d7#b5a106be3133979d78832a0ca3f4d6b57613b3d7"
+dependencies = [
+ "num-complex",
+ "num-traits",
+ "strided-kernel",
+ "tenferro-tensor",
+ "thiserror 2.0.20",
 ]
 
 [[package]]
 name = "tenferro-internal-extension-macros"
-version = "0.1.0"
-source = "git+https://github.com/tensor4all/tenferro-rs.git?rev=f68b33d9e06cc0af8d2b9c40f7c52d8972cd1266#f68b33d9e06cc0af8d2b9c40f7c52d8972cd1266"
+version = "0.3.0"
+source = "git+https://github.com/tensor4all/tenferro-rs.git?rev=b5a106be3133979d78832a0ca3f4d6b57613b3d7#b5a106be3133979d78832a0ca3f4d6b57613b3d7"
 dependencies = [
  "proc-macro2",
  "quote",
@@ -1552,8 +1459,8 @@ dependencies = [
 
 [[package]]
 name = "tenferro-internal-ops"
-version = "0.1.0"
-source = "git+https://github.com/tensor4all/tenferro-rs.git?rev=f68b33d9e06cc0af8d2b9c40f7c52d8972cd1266#f68b33d9e06cc0af8d2b9c40f7c52d8972cd1266"
+version = "0.3.0"
+source = "git+https://github.com/tensor4all/tenferro-rs.git?rev=b5a106be3133979d78832a0ca3f4d6b57613b3d7#b5a106be3133979d78832a0ca3f4d6b57613b3d7"
 dependencies = [
  "computegraph",
  "num-complex",
@@ -1561,17 +1468,15 @@ dependencies = [
  "tenferro-internal-extension-macros",
  "tenferro-tensor",
  "thiserror 2.0.20",
- "tidu",
 ]
 
 [[package]]
 name = "tenferro-linalg"
-version = "0.1.0"
-source = "git+https://github.com/tensor4all/tenferro-rs.git?rev=f68b33d9e06cc0af8d2b9c40f7c52d8972cd1266#f68b33d9e06cc0af8d2b9c40f7c52d8972cd1266"
+version = "0.3.0"
+source = "git+https://github.com/tensor4all/tenferro-rs.git?rev=b5a106be3133979d78832a0ca3f4d6b57613b3d7#b5a106be3133979d78832a0ca3f4d6b57613b3d7"
 dependencies = [
  "computegraph",
  "faer",
- "lapack",
  "num-complex",
  "num-traits",
  "tenferro-ad",
@@ -1580,18 +1485,19 @@ dependencies = [
  "tenferro-internal-ops",
  "tenferro-runtime",
  "tenferro-tensor",
- "tidu",
+ "thiserror 2.0.20",
 ]
 
 [[package]]
 name = "tenferro-runtime"
-version = "0.1.0"
-source = "git+https://github.com/tensor4all/tenferro-rs.git?rev=f68b33d9e06cc0af8d2b9c40f7c52d8972cd1266#f68b33d9e06cc0af8d2b9c40f7c52d8972cd1266"
+version = "0.3.0"
+source = "git+https://github.com/tensor4all/tenferro-rs.git?rev=b5a106be3133979d78832a0ca3f4d6b57613b3d7#b5a106be3133979d78832a0ca3f4d6b57613b3d7"
 dependencies = [
  "computegraph",
  "lru",
  "num-complex",
  "num-traits",
+ "smallvec",
  "tenferro-core-ops",
  "tenferro-internal-ops",
  "tenferro-tensor",
@@ -1600,12 +1506,13 @@ dependencies = [
 
 [[package]]
 name = "tenferro-tensor"
-version = "0.1.0"
-source = "git+https://github.com/tensor4all/tenferro-rs.git?rev=f68b33d9e06cc0af8d2b9c40f7c52d8972cd1266#f68b33d9e06cc0af8d2b9c40f7c52d8972cd1266"
+version = "0.3.0"
+source = "git+https://github.com/tensor4all/tenferro-rs.git?rev=b5a106be3133979d78832a0ca3f4d6b57613b3d7#b5a106be3133979d78832a0ca3f4d6b57613b3d7"
 dependencies = [
  "num-complex",
  "num-traits",
  "smallvec",
+ "strided-kernel",
  "tenferro-core-ops",
  "tenferro-tensor-core",
  "thiserror 2.0.20",
@@ -1613,8 +1520,8 @@ dependencies = [
 
 [[package]]
 name = "tenferro-tensor-core"
-version = "0.1.0"
-source = "git+https://github.com/tensor4all/tenferro-rs.git?rev=f68b33d9e06cc0af8d2b9c40f7c52d8972cd1266#f68b33d9e06cc0af8d2b9c40f7c52d8972cd1266"
+version = "0.3.0"
+source = "git+https://github.com/tensor4all/tenferro-rs.git?rev=b5a106be3133979d78832a0ca3f4d6b57613b3d7#b5a106be3133979d78832a0ca3f4d6b57613b3d7"
 dependencies = [
  "num-complex",
  "smallvec",
@@ -1624,14 +1531,14 @@ dependencies = [
 [[package]]
 name = "tensor4all-aci"
 version = "0.2.0"
-source = "git+https://github.com/tensor4all/tensor4all-rs?rev=c9ecb7f0d4eec4bbf6db51dc7648a4b2d84b4f34#c9ecb7f0d4eec4bbf6db51dc7648a4b2d84b4f34"
+source = "git+https://github.com/tensor4all/tensor4all-rs?rev=9e9aedaebe0d3918b34dd399ff0981e337f3835b#9e9aedaebe0d3918b34dd399ff0981e337f3835b"
 dependencies = [
  "num-complex",
- "rand 0.9.5",
+ "rand",
  "rand_chacha",
  "rand_distr",
+ "tensor4all-core",
  "tensor4all-simplett",
- "tensor4all-tcicore",
  "tensor4all-tensorbackend",
  "thiserror 2.0.20",
 ]
@@ -1639,15 +1546,17 @@ dependencies = [
 [[package]]
 name = "tensor4all-core"
 version = "0.2.0"
-source = "git+https://github.com/tensor4all/tensor4all-rs?rev=c9ecb7f0d4eec4bbf6db51dc7648a4b2d84b4f34#c9ecb7f0d4eec4bbf6db51dc7648a4b2d84b4f34"
+source = "git+https://github.com/tensor4all/tensor4all-rs?rev=9e9aedaebe0d3918b34dd399ff0981e337f3835b#9e9aedaebe0d3918b34dd399ff0981e337f3835b"
 dependencies = [
  "anyhow",
+ "bnum",
  "dyn-clone",
  "num-complex",
  "num-traits",
  "omeco",
+ "paste",
  "petgraph",
- "rand 0.9.5",
+ "rand",
  "rand_distr",
  "smallvec",
  "tenferro-ad",
@@ -1655,7 +1564,6 @@ dependencies = [
  "tenferro-linalg",
  "tenferro-runtime",
  "tenferro-tensor",
- "tensor4all-tcicore",
  "tensor4all-tensorbackend",
  "thiserror 2.0.20",
 ]
@@ -1663,7 +1571,7 @@ dependencies = [
 [[package]]
 name = "tensor4all-hdf5"
 version = "0.2.0"
-source = "git+https://github.com/tensor4all/tensor4all-rs?rev=c9ecb7f0d4eec4bbf6db51dc7648a4b2d84b4f34#c9ecb7f0d4eec4bbf6db51dc7648a4b2d84b4f34"
+source = "git+https://github.com/tensor4all/tensor4all-rs?rev=9e9aedaebe0d3918b34dd399ff0981e337f3835b#9e9aedaebe0d3918b34dd399ff0981e337f3835b"
 dependencies = [
  "anyhow",
  "hdf5-metno",
@@ -1679,85 +1587,78 @@ dependencies = [
 [[package]]
 name = "tensor4all-itensorlike"
 version = "0.2.0"
-source = "git+https://github.com/tensor4all/tensor4all-rs?rev=c9ecb7f0d4eec4bbf6db51dc7648a4b2d84b4f34#c9ecb7f0d4eec4bbf6db51dc7648a4b2d84b4f34"
+source = "git+https://github.com/tensor4all/tensor4all-rs?rev=9e9aedaebe0d3918b34dd399ff0981e337f3835b#9e9aedaebe0d3918b34dd399ff0981e337f3835b"
 dependencies = [
  "anyhow",
  "num-complex",
  "petgraph",
- "rand 0.9.5",
+ "rand",
  "tensor4all-core",
  "tensor4all-treetn",
  "thiserror 2.0.20",
 ]
 
 [[package]]
-name = "tensor4all-partitionedtt"
+name = "tensor4all-partitionedtreetn"
 version = "0.2.0"
-source = "git+https://github.com/tensor4all/tensor4all-rs?rev=c9ecb7f0d4eec4bbf6db51dc7648a4b2d84b4f34#c9ecb7f0d4eec4bbf6db51dc7648a4b2d84b4f34"
+source = "git+https://github.com/tensor4all/tensor4all-rs?rev=9e9aedaebe0d3918b34dd399ff0981e337f3835b#9e9aedaebe0d3918b34dd399ff0981e337f3835b"
 dependencies = [
  "anyhow",
- "num-complex",
- "num-traits",
- "rand 0.9.5",
  "tensor4all-core",
- "tensor4all-itensorlike",
- "tensor4all-simplett",
- "tensor4all-tcicore",
  "tensor4all-tensorbackend",
- "tensor4all-tensorci",
  "tensor4all-treetn",
  "thiserror 2.0.20",
 ]
 
 [[package]]
-name = "tensor4all-quanticstci"
+name = "tensor4all-partitionedtt"
 version = "0.2.0"
-source = "git+https://github.com/tensor4all/tensor4all-rs?rev=c9ecb7f0d4eec4bbf6db51dc7648a4b2d84b4f34#c9ecb7f0d4eec4bbf6db51dc7648a4b2d84b4f34"
+source = "git+https://github.com/tensor4all/tensor4all-rs?rev=9e9aedaebe0d3918b34dd399ff0981e337f3835b#9e9aedaebe0d3918b34dd399ff0981e337f3835b"
 dependencies = [
  "anyhow",
  "num-complex",
  "num-traits",
- "quanticsgrids",
- "rand 0.9.5",
+ "rand",
  "tensor4all-core",
+ "tensor4all-itensorlike",
  "tensor4all-simplett",
- "tensor4all-tcicore",
  "tensor4all-tensorbackend",
- "tensor4all-treetci",
+ "tensor4all-tensorci",
  "tensor4all-treetn",
  "thiserror 2.0.20",
 ]
 
 [[package]]
-name = "tensor4all-simplett"
+name = "tensor4all-quanticstci"
 version = "0.2.0"
-source = "git+https://github.com/tensor4all/tensor4all-rs?rev=c9ecb7f0d4eec4bbf6db51dc7648a4b2d84b4f34#c9ecb7f0d4eec4bbf6db51dc7648a4b2d84b4f34"
+source = "git+https://github.com/tensor4all/tensor4all-rs?rev=9e9aedaebe0d3918b34dd399ff0981e337f3835b#9e9aedaebe0d3918b34dd399ff0981e337f3835b"
 dependencies = [
  "anyhow",
- "bnum",
  "num-complex",
  "num-traits",
- "rand 0.9.5",
- "tenferro-einsum",
- "tenferro-tensor",
+ "quanticsgrids",
+ "rand",
  "tensor4all-core",
- "tensor4all-tcicore",
+ "tensor4all-simplett",
  "tensor4all-tensorbackend",
+ "tensor4all-treetci",
+ "tensor4all-treetn",
  "thiserror 2.0.20",
 ]
 
 [[package]]
-name = "tensor4all-tcicore"
+name = "tensor4all-simplett"
 version = "0.2.0"
-source = "git+https://github.com/tensor4all/tensor4all-rs?rev=c9ecb7f0d4eec4bbf6db51dc7648a4b2d84b4f34#c9ecb7f0d4eec4bbf6db51dc7648a4b2d84b4f34"
+source = "git+https://github.com/tensor4all/tensor4all-rs?rev=9e9aedaebe0d3918b34dd399ff0981e337f3835b#9e9aedaebe0d3918b34dd399ff0981e337f3835b"
 dependencies = [
  "anyhow",
  "bnum",
  "num-complex",
  "num-traits",
- "paste",
+ "rand",
  "tenferro-einsum",
  "tenferro-tensor",
+ "tensor4all-core",
  "tensor4all-tensorbackend",
  "thiserror 2.0.20",
 ]
@@ -1765,13 +1666,13 @@ dependencies = [
 [[package]]
 name = "tensor4all-tensorbackend"
 version = "0.2.0"
-source = "git+https://github.com/tensor4all/tensor4all-rs?rev=c9ecb7f0d4eec4bbf6db51dc7648a4b2d84b4f34#c9ecb7f0d4eec4bbf6db51dc7648a4b2d84b4f34"
+source = "git+https://github.com/tensor4all/tensor4all-rs?rev=9e9aedaebe0d3918b34dd399ff0981e337f3835b#9e9aedaebe0d3918b34dd399ff0981e337f3835b"
 dependencies = [
  "anyhow",
  "num-complex",
  "num-traits",
  "omeco",
- "rand 0.9.5",
+ "rand",
  "rand_distr",
  "tenferro-ad",
  "tenferro-cpu",
@@ -1785,14 +1686,14 @@ dependencies = [
 [[package]]
 name = "tensor4all-tensorci"
 version = "0.2.0"
-source = "git+https://github.com/tensor4all/tensor4all-rs?rev=c9ecb7f0d4eec4bbf6db51dc7648a4b2d84b4f34#c9ecb7f0d4eec4bbf6db51dc7648a4b2d84b4f34"
+source = "git+https://github.com/tensor4all/tensor4all-rs?rev=9e9aedaebe0d3918b34dd399ff0981e337f3835b#9e9aedaebe0d3918b34dd399ff0981e337f3835b"
 dependencies = [
  "anyhow",
  "num-complex",
  "num-traits",
- "rand 0.9.5",
+ "rand",
+ "tensor4all-core",
  "tensor4all-simplett",
- "tensor4all-tcicore",
  "tensor4all-tensorbackend",
  "thiserror 2.0.20",
 ]
@@ -1800,16 +1701,15 @@ dependencies = [
 [[package]]
 name = "tensor4all-treetci"
 version = "0.2.0"
-source = "git+https://github.com/tensor4all/tensor4all-rs?rev=c9ecb7f0d4eec4bbf6db51dc7648a4b2d84b4f34#c9ecb7f0d4eec4bbf6db51dc7648a4b2d84b4f34"
+source = "git+https://github.com/tensor4all/tensor4all-rs?rev=9e9aedaebe0d3918b34dd399ff0981e337f3835b#9e9aedaebe0d3918b34dd399ff0981e337f3835b"
 dependencies = [
  "anyhow",
  "num-complex",
  "petgraph",
- "rand 0.9.5",
+ "rand",
  "tenferro-linalg",
  "tenferro-tensor",
  "tensor4all-core",
- "tensor4all-tcicore",
  "tensor4all-tensorbackend",
  "tensor4all-treetn",
  "thiserror 2.0.20",
@@ -1818,13 +1718,13 @@ dependencies = [
 [[package]]
 name = "tensor4all-treetn"
 version = "0.2.0"
-source = "git+https://github.com/tensor4all/tensor4all-rs?rev=c9ecb7f0d4eec4bbf6db51dc7648a4b2d84b4f34#c9ecb7f0d4eec4bbf6db51dc7648a4b2d84b4f34"
+source = "git+https://github.com/tensor4all/tensor4all-rs?rev=9e9aedaebe0d3918b34dd399ff0981e337f3835b#9e9aedaebe0d3918b34dd399ff0981e337f3835b"
 dependencies = [
  "anyhow",
  "dyn-clone",
  "num-complex",
  "petgraph",
- "rand 0.9.5",
+ "rand",
  "tensor4all-core",
  "tensor4all-simplett",
  "tensor4all-tensorbackend",
@@ -1880,14 +1780,6 @@ dependencies = [
  "cfg-if",
 ]
 
-[[package]]
-name = "tidu"
-version = "0.1.0"
-source = "git+https://github.com/tensor4all/tidu-rs.git?rev=9934d8b775fcfbbad8e2daf3f9af54477da63e35#9934d8b775fcfbbad8e2daf3f9af54477da63e35"
-dependencies = [
- "computegraph",
-]
-
 [[package]]
 name = "toml_datetime"
 version = "1.1.1+spec-1.1.0"
@@ -1967,12 +1859,6 @@ dependencies = [
  "tracing-log",
 ]
 
-[[package]]
-name = "ucd-trie"
-version = "0.1.7"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "2896d95c02a80c6d6a5d6e953d479f5ddf2dfdb6a244441010e373ac0fb88971"
-
 [[package]]
 name = "unicode-ident"
 version = "1.0.24"
diff --git a/Cargo.toml b/Cargo.toml
index 16c6147..5f3cce8 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -6,20 +6,22 @@ license = "MIT"
 publish = false
 
 [dependencies]
-tensor4all-simplett = { git = "https://github.com/tensor4all/tensor4all-rs", rev = "c9ecb7f0d4eec4bbf6db51dc7648a4b2d84b4f34", package = "tensor4all-simplett" }
-tensor4all-aci = { git = "https://github.com/tensor4all/tensor4all-rs", rev = "c9ecb7f0d4eec4bbf6db51dc7648a4b2d84b4f34", package = "tensor4all-aci" }
-tensor4all-treetn = { git = "https://github.com/tensor4all/tensor4all-rs", rev = "c9ecb7f0d4eec4bbf6db51dc7648a4b2d84b4f34", package = "tensor4all-treetn" }
-tensor4all-itensorlike = { git = "https://github.com/tensor4all/tensor4all-rs", rev = "c9ecb7f0d4eec4bbf6db51dc7648a4b2d84b4f34", package = "tensor4all-itensorlike" }
-tensor4all-partitionedtt = { git = "https://github.com/tensor4all/tensor4all-rs", rev = "c9ecb7f0d4eec4bbf6db51dc7648a4b2d84b4f34", package = "tensor4all-partitionedtt" }
-tensor4all-quanticstci = { git ="https://github.com/tensor4all/tensor4all-rs", rev = "c9ecb7f0d4eec4bbf6db51dc7648a4b2d84b4f34", package = "tensor4all-quanticstci" }
-tensor4all-hdf5 = { git = "https://github.com/tensor4all/tensor4all-rs", rev = "c9ecb7f0d4eec4bbf6db51dc7648a4b2d84b4f34", package = "tensor4all-hdf5" }
-tensor4all-core = { git = "https://github.com/tensor4all/tensor4all-rs", rev = "c9ecb7f0d4eec4bbf6db51dc7648a4b2d84b4f34", package = "tensor4all-core" }
+tensor4all-simplett = { git = "https://github.com/tensor4all/tensor4all-rs", rev = "9e9aedaebe0d3918b34dd399ff0981e337f3835b", package = "tensor4all-simplett" }
+tensor4all-aci = { git = "https://github.com/tensor4all/tensor4all-rs", rev = "9e9aedaebe0d3918b34dd399ff0981e337f3835b", package = "tensor4all-aci" }
+tensor4all-treetn = { git = "https://github.com/tensor4all/tensor4all-rs", rev = "9e9aedaebe0d3918b34dd399ff0981e337f3835b", package = "tensor4all-treetn" }
+tensor4all-itensorlike = { git = "https://github.com/tensor4all/tensor4all-rs", rev = "9e9aedaebe0d3918b34dd399ff0981e337f3835b", package = "tensor4all-itensorlike" }
+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-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"
 rand = "0.9"
 rand_chacha = "0.9"
 serde = { version = "1", features = ["derive"] }
 serde_json = "1"
 anyhow = "1"
+libm = "0.2"
 
 [profile.release]
 debug = true
diff --git a/README.md b/README.md
index ccdd791..ce8079e 100644
--- a/README.md
+++ b/README.md
@@ -15,7 +15,8 @@ revision that produced it. The pinned `tensor4all-rs` revision lives in `Cargo.t
 | 2. `mpo_mpo_quantics` | Contraction of two 2D quantics Gaussian-mixture MPOs over their shared variable, swept over bits per variable `R`, against the closed-form Gaussian integral | [description](#case-2-mpo-mpo-contraction-of-2d-quantics-gaussian-mixtures) | [`src/bin/mpo_mpo_quantics.rs`](src/bin/mpo_mpo_quantics.rs) | [`result/mac-cpu/mpo_mpo_quantics.md`](result/mac-cpu/mpo_mpo_quantics.md) | [time](result/mac-cpu/mpo_mpo_quantics-time.svg), [error](result/mac-cpu/mpo_mpo_quantics-error.svg) |
 | 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) | not yet recorded in any profile | not yet recorded in any profile |
+| 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) |
 
 ## Benchmark cases
 
@@ -510,6 +511,55 @@ it stays one `BENCH_NS` away.
 Runner: [`src/bin/elementwise_gauss2d_patched.rs`](src/bin/elementwise_gauss2d_patched.rs),
 sweep over `N`.
 
+### Case 6: patched MPO-MPO fit contraction
+
+This case reuses case 5's constant-density anisotropic Gaussian family for
+`F(x,y)` and `G(y,z)`. Each function is built once as a quantics MPO. The same
+cores and the same shared `y` index objects are copied into
+`tensor4all-partitionedtreetn`, using a linear chain and its last site as the
+explicit center. Both inputs are split with `add_with_patching`, sequentially
+over the shared `y` bits first, and the patched output is computed with
+`contract_adaptive` using one full fit sweep. A global `fit_treetn` arm
+contracts the same unpatched MPOs with the same output cap and sweep count.
+Only the two contraction calls are timed. Input interpolation, HDF5 cache
+loading, bridging, patching and output format conversion are outside the timed
+region. Deterministic global input trains are cached under `.cache/inputs/` and
+reused across reruns; patched inputs are always rebuilt from the cached global
+trains. Set `BENCH_INPUT_CACHE_REFRESH=1` to rebuild an entry or
+`BENCH_INPUT_CACHE_DIR` to move the cache. A focused test retains the legacy TT
+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
+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
+`chi_in`; every local and merged patched contraction is capped at
+`BENCH_PATCH_MAX_BOND`. The global arm shows the effect of patching. Naive and
+zip-up arms are excluded:
+naive scales too steeply at this rank, while the purpose here is the patched fit
+path.
+
+Accuracy is checked against the finite quantics `y` grid, not against an
+unbounded integral. Each anisotropic Gaussian pair is integrated in closed form
+on the box and corrected to the left-endpoint grid sum with the first two
+Euler-Maclaurin endpoint terms. A unit test compares this reference with the
+explicit grid sum.
+
+Run single-threaded on one chosen CPU core, for example core 7:
+
+```bash
+RAYON_NUM_THREADS=1 OMP_NUM_THREADS=1 OPENBLAS_NUM_THREADS=1 \
+MKL_NUM_THREADS=1 VECLIB_MAXIMUM_THREADS=1 \
+OUT_DIR=/tmp/mpo-mpo-aniso-patched \
+taskset -c 7 cargo run --release --bin mpo_mpo_aniso_patched
+```
+
+Runner: [`src/bin/mpo_mpo_aniso_patched.rs`](src/bin/mpo_mpo_aniso_patched.rs),
+one default point at `N = 512`. This exploratory case is standalone and is not
+included in `scripts/run_all.sh`.
+
 ## Latest results
 
 One profile per physical machine, so numbers from different hardware never overwrite
@@ -726,18 +776,24 @@ Environment knobs:
 | `BENCH_KS` | case 1 | `4,8,16,32,64,128` | comma-separated Fourier mode counts `K` to sweep |
 | `BENCH_R` | case 1 | `20` | number of quantics bits |
 | `BENCH_RS` | cases 2 and 3 | `6,8,10,12,14` | comma-separated bits per variable `R` to sweep |
-| `BENCH_NGAUSS` | cases 2 and 3 | `8` | number of Gaussians per mixture |
+| `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_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 |
-| `BENCH_ANISO_SIGMA` | case 5, `aniso` family | `0.05` | minor width of every spike, in box units. This is the length the grid has to resolve: `R` is chosen so a grid step is at most a quarter of it |
-| `BENCH_ANISO_RHO_MAX` | case 5, `aniso` family | `8.0` | upper end of the aspect ratio, drawn log-uniform in `[1, BENCH_ANISO_RHO_MAX]` per spike. Setting it to `1` is legal and is the isotropic control: circular spikes of one common shape, everything else unchanged |
-| `BENCH_ANISO_SPACING` | case 5, `aniso` family | `3.0` | mean spacing between spikes, in minor widths, held fixed as `N` grows. It fixes both the box, `L = BENCH_ANISO_SPACING * sigma * sqrt(N) / 2`, and the overlap of the two mixtures, which is what keeps the product non degenerate |
+| `BENCH_ANISO_SIGMA` | case 5 `aniso` family and case 6 | `0.05` | minor width of every spike, in box units. This is the length the grid has to resolve: `R` is chosen so a grid step is at most a quarter of it |
+| `BENCH_ANISO_RHO_MAX` | case 5 `aniso` family and case 6 | `8.0` | upper end of the aspect ratio, drawn log-uniform in `[1, BENCH_ANISO_RHO_MAX]` per spike. Setting it to `1` is legal and is the isotropic control: circular spikes of one common shape, everything else unchanged |
+| `BENCH_ANISO_SPACING` | case 5 `aniso` family and case 6 | `3.0` | mean spacing between spikes, in minor widths, held fixed as `N` grows. It fixes both the box, `L = BENCH_ANISO_SPACING * sigma * sqrt(N) / 2`, and the overlap of the two mixtures |
+| `BENCH_BOX_PADDING` | case 6 | `1.0` | multiplier applied to the case-5 box after drawing the spikes. Keep at 1 for the identical family; larger values add an empty integration margin |
+| `BENCH_R_EXTRA` | case 6 | `0` | extra quantics bits beyond case 5's quarter-minor-width resolution |
+| `BENCH_MAX_INPUT_CHI` | case 6 | `256` | safety ceiling checked before patch construction or fit; the runner stops instead of starting a contraction above this rank |
+| `BENCH_ERROR_SAMPLES` | case 6 | `128` | sampled `(x,z)` points used for the grid-reference error |
 | `BENCH_L0` | cases 4 and 5 (`smooth` family) | `6.0` | reference box half-width, the `L` at `N` = `BENCH_N0` |
 | `BENCH_N0` | cases 4 and 5 (`smooth` family) | `8` | reference Gaussian count, the `N` at which `L` = `BENCH_L0` and `R` = `BENCH_R0` |
 | `BENCH_R0` | cases 4 and 5 (`smooth` family) | `10` | reference bits per variable, the `R` at `N` = `BENCH_N0`. Lower it for a cheap probe of the whole sweep |
-| `BENCH_RTOL` | case 5 | `1e-8` | the one accuracy knob of case 5: the absolute-per-patch input TCI tolerance of the patched arms, the tolerance of the global baseline inputs, the absolute pivot budget of both `aci` arms, and the global output budget of both products. Every arm is compared at this one value |
-| `BENCH_PATCH_MAX_BOND` | case 5 | `64` | per-patch rank cap of the patched input construction. A subdomain that does not fit under it is split again, so this is what decides how deep the patch tree goes: a smaller value means more, smaller patches |
+| `BENCH_RTOL` | cases 5 and 6 | `1e-8` | case 5's one accuracy knob: the absolute-per-patch input TCI tolerance, global baseline input tolerance, absolute pivot budget of both `aci` arms, and global output budget of both products. Case 6 uses it for global patching, output truncation, and the fit SVD policy |
+| `BENCH_PATCH_MAX_BOND` | cases 5 and 6 | `64` (case 5), `128` (case 6) | per-patch rank cap of the patched input construction. A subdomain that does not fit under it is split again; case 6 no longer measures cap 64 |
 | `BENCH_PATCH_INPUT` | case 5 | `norm` | which construction produces the patched inputs. `norm` builds one global train per input, exactly as case 4 does, and splits it with `partitionedtt::add_with_patching`, whose decisions come from Frobenius norms alone. `tci` runs `partitionedtt::adaptiveinterpolate` instead, which never forms a global train and splits a patch whose own TCI does not converge under the cap; it is the construction the case is eventually written for, it works at the current pin (known issue 11 records the upstream fix), and it is not the default because for the same cap it splits far harder and returns a much larger representation. Recorded as `input_path` |
 | `BENCH_PATCH_SPLIT` | case 5 | `gain` | how the `norm` path picks the site to split. `gain` is the upstream `ExactParameterGain`: it forms and budget-truncates the children of every candidate site and keeps the cheapest. `sequential` takes the first unprojected site of the patch order instead, so the splitting runs strictly coarse to fine and a patch is a single quadrant. Ignored by the `tci` path, whose splitting is sequential by construction. Recorded as `split_strategy` |
 | `BENCH_PATCH_MAX_ITER` | case 5, `tci` path | `20` | half-sweep limit of each patch's own TCI run on the `tci` path, the upstream default. It exists as a knob because a run that stops at its limit is not converged and its patch is split whatever its rank was, which would make the patch tree a measurement of the limit; measured, it does not bind, and raising it to 200 changed neither the patch counts nor the patch ranks at six times the build cost. Ignored by the default `norm` path, which runs no TCI, and recorded as `patch_max_iter` either way |
@@ -746,14 +802,14 @@ Environment knobs:
 | `BENCH_BASELINES` | case 5 | `1` | whether to also measure the two global arms at the same `rtol`. Each has its own `N` ceiling whatever this says, since they are two cost classes: on the `smooth` family both stop at `N` = 64, and on `aniso` the interpolating `aci` arm runs everywhere while the uncapped global fit stops at `N` = 1024, the last point measured to stay under two minutes for it |
 | `BENCH_ALPHA_LO` | cases 2, 3, 4 and 5 (`smooth` family) | `0.5` | lower bound of the Gaussian width parameter |
 | `BENCH_ALPHA_HI` | cases 2, 3, 4 and 5 (`smooth` family) | `8.0` | upper bound of the Gaussian width parameter |
-| `BENCH_SANITY` | cases 2, 3, 4 and 5 | `1e-2` | relative error gate. Cases 2 and 5 apply it to every algorithm; cases 3 and 4 apply it to all but `zipup_treetn`, which is gated at a hardcoded `5.0` |
+| `BENCH_SANITY` | cases 2, 3, 4, 5 and 6 | `1e-2` (cases 2 to 5), `1e-4` (case 6) | relative error gate. Cases 2 and 5 apply it to every algorithm; cases 3 and 4 apply it to all but `zipup_treetn`, which is gated at a hardcoded `5.0`; case 6 applies it to both global and patched grid-reference errors |
 | `BENCH_TOL` | cases 1, 2, 3 and 4 | `1e-8` | instance tolerance. In case 1 it is the working tolerance of the whole case, both the input compression and the product. In cases 2, 3 and 4 it is scoped to the input TCI construction only, where it fixes `chi_in` and therefore the instance; the arms take `BENCH_CONTRACT_TOL` instead. It is what each record's top-level `tolerance` field and the Julia-check metadata report, since both describe the inputs. Case 5 has no separate instance tolerance: `BENCH_RTOL` is its one knob and is what its records report |
 | `BENCH_CONTRACT_TOL` | cases 2, 3 and 4 | `1e-15` | truncation tolerance handed to every contraction or product arm, recorded as `contract_tol` in the params of each record. At the default it never fires, so the rank cap `chi_in` is the only binding truncation control and the fixed-budget cases are `chi_out`-driven by construction. Raise it if you want a tolerance-driven variant of those cases |
-| `BENCH_MAX_BOND` | all | `4096` (case 1), `512` (cases 2, 3, 4 and 5) | bond dimension cap. In cases 2, 3 and 4 it caps only the input TCI construction, since the arms themselves run at the fixed output budget `chi_in`. Case 4 fails rather than reports if an instance reaches it, since `chi_in` is what that case measures. In case 5 it caps the global input trains, which are the global baselines' inputs and, on the default `norm` path, also the trains the patches are split out of; the patched counterpart is `BENCH_PATCH_MAX_BOND`. Reaching it is likewise a failure, since a cap-limited train is no longer tolerance-driven |
-| `BENCH_RUNS` | all | `5` (case 1), `3` (cases 2, 3 and 4), `1` (case 5) | timed repetitions, the median is reported. Case 5 defaults to a single pass because it is the most expensive case per point while its run to run spread is a few percent, smaller than the spread between two constructions of the same instance |
-| `BENCH_WARMUPS` | all | `1` (case 1), `0` (cases 2, 3, 4 and 5) | untimed warmup repetitions |
+| `BENCH_MAX_BOND` | all | `4096` (case 1), `512` (cases 2 to 5), `384` (case 6) | bond cap for input construction. In cases 2 to 4 it caps only input TCI; case 4 fails if it binds. In case 5 it caps both global baseline inputs and the global trains split by the default `norm` path, and binding is likewise a failure because the train is no longer tolerance-driven. Case 6 records cap-limited crossover points explicitly and separately stops when measured rank exceeds `BENCH_MAX_INPUT_CHI` |
+| `BENCH_RUNS` | all | `5` (case 1), `3` (cases 2, 3 and 4), `1` (cases 5 and 6) | timed repetitions, the median is reported |
+| `BENCH_WARMUPS` | all | `1` (case 1), `0` (cases 2 to 6) | untimed warmup repetitions |
 | `BENCH_SEED` | all | `0` | base seed for instance generation |
-| `BENCH_ALGOS` | all | `naive,zipup,fit,aci` (case 1), `naive,zipup_simplett,zipup_treetn,fit_treetn` (case 2), `naive,zipup_treetn,fit_treetn,aci` (case 3), `zipup_treetn,fit_treetn,aci` (case 4), `patched_fit_treetn,patched_naive,patched_aci` (case 5, whose fourth arm `patched_zipup_treetn` is excluded from the defaults on cost and the global baselines are controlled by `BENCH_BASELINES` instead) | comma-separated algorithms to run |
+| `BENCH_ALGOS` | cases 1 to 5 | `naive,zipup,fit,aci` (case 1), `naive,zipup_simplett,zipup_treetn,fit_treetn` (case 2), `naive,zipup_treetn,fit_treetn,aci` (case 3), `zipup_treetn,fit_treetn,aci` (case 4), `patched_fit_treetn,patched_naive,patched_aci` (case 5) | comma-separated algorithms to run. Case 6 always runs one global and one patched fit arm |
 | `OUT_DIR` | all | `result/dev/raw` | directory for the `RunRecord` JSON files |
 | `EXPORT_HDF5` | cases 1, 2 and 3 | unset | directory for ITensors-compatible HDF5 instance dumps, plus their JSON metadata. Set it to enable the Julia checks. An empty value counts as unset. Cases 2 and 3 use the same file names, so give them separate directories. Cases 4 and 5 export nothing and ignore it |
 
diff --git a/result/linux-epyc-7713p/mpo_mpo_aniso_patched.md b/result/linux-epyc-7713p/mpo_mpo_aniso_patched.md
new file mode 100644
index 0000000..f2d38fb
--- /dev/null
+++ b/result/linux-epyc-7713p/mpo_mpo_aniso_patched.md
@@ -0,0 +1,54 @@
+# Patched MPO-MPO fit contraction
+
+Profile metadata: [`run.yaml`](run.yaml). Raw records: [`raw/`](raw/).
+
+## Conditions
+
+- `N = 2048`, `R = 10`, relative squared-value discarded-tail cutoff `1e-8`
+- patch bond cap `128`; cap `64` is intentionally excluded
+- one pinned CPU core; Rayon, OpenMP, OpenBLAS, MKL and vecLib limited to one thread
+- two measured runs, no warmup
+- identical cached global MPO inputs for the global and patched arms
+- input generation, HDF5 loading, patch construction and output conversion excluded from contraction timing
+
+## Results
+
+The reported median is the harness's upper median for two runs.
+
+| Input χ | Input patches (left/right) | Global runs [s] | Global median [s] | Patched runs [s] | Patched median [s] | Global / patched | Output patches | Max output-patch χ | Patched error |
+| ---: | ---: | --- | ---: | --- | ---: | ---: | ---: | ---: | ---: |
+| 128 | 1 / 1 | 3.576, 3.515 | 3.576 | 3.321, 3.292 | 3.321 | 1.08× | 1 | 78 | 3.455e-4 |
+| 192 | 8 / 8 | 11.791, 11.682 | 11.791 | 13.353, 13.244 | 13.353 | 0.88× | 1 | 68 | 5.547e-4 |
+| 224 | 8 / 8 | 15.348, 15.468 | 15.468 | 13.074, 13.103 | 13.103 | 1.18× | 1 | 68 | 4.917e-4 |
+| 256 | 8 / 8 | 19.606, 19.558 | 19.606 | 12.274, 12.096 | 12.274 | **1.60×** | 1 | 68 | 6.281e-4 |
+
+The meaningful crossover is between input χ `192` and `224`: patched takes
+1.13× the global time at χ `192`, while global takes 1.18× the patched time at
+χ `224` and 1.60× at χ `256`. The χ `128` point has one input patch per side
+and therefore does not measure a partitioning advantage.
+
+The χ `128`, `192`, and `224` inputs hit their respective TCI construction caps;
+they are cap-limited views of the same instance, while the χ `256` input was
+built with cap `384` and converged below it.
+
+## Input MPO size
+
+Pure `f64` core storage, excluding index and HDF5 metadata:
+
+| Input χ | Left values / MiB | Right values / MiB | Pair MiB | HDF5 MiB |
+| ---: | ---: | ---: | ---: | ---: |
+| 128 | 166,440 / 1.270 | 167,552 / 1.278 | 2.548 | 2.930 |
+| 192 | 282,964 / 2.159 | 281,472 / 2.147 | 4.306 | 4.689 |
+| 224 | 314,912 / 2.403 | 314,912 / 2.403 | 4.805 | 5.187 |
+| 256 | 346,176 / 2.641 | 347,976 / 2.655 | 5.296 | 5.678 |
+
+## Accuracy status
+
+All points exceed the runner's unchanged `1e-4` sanity gate under this
+ITensor-style cutoff, so each process exits nonzero after writing its records.
+The records are retained because the comparison and crossover are valid, but
+these runs do **not** claim that the accuracy gate passed. No tolerance or gate
+was relaxed.
+
+The recorded rerun uses the merged tensor4all-rs revision `9e9aeda`, which is
+also the revision pinned by this repository.
diff --git a/result/linux-epyc-7713p/raw/mpo_mpo_aniso_patched-global_fit_treetn-n2048-chi128.json b/result/linux-epyc-7713p/raw/mpo_mpo_aniso_patched-global_fit_treetn-n2048-chi128.json
new file mode 100644
index 0000000..f664c9a
--- /dev/null
+++ b/result/linux-epyc-7713p/raw/mpo_mpo_aniso_patched-global_fit_treetn-n2048-chi128.json
@@ -0,0 +1,63 @@
+{
+  "schema_version": 1,
+  "case": "mpo_mpo_aniso_patched",
+  "algorithm": "global_fit_treetn",
+  "params": {
+    "box_l": 3.394112549695429,
+    "box_padding": 1.0,
+    "contraction_svd_policy": "relative_squared_discarded_tail_sum",
+    "error_metric": "max_rel_vs_quantics_grid_reference",
+    "error_seed": 99,
+    "fit_nfullsweeps": 1,
+    "global_contract_max_bond": 128,
+    "inner_l": 3.394112549695429,
+    "input_cache_hit": true,
+    "input_cache_key": "aniso-mpo-v1-n2048-r10-extra0-spacing4008000000000000-padding3ff0000000000000-inner400b27247aff1491-box400b27247aff1491-sigma3fa999999999999a-rho4020000000000000-rtol3e45798ee2308c3a-cap128-seed0",
+    "input_cache_load_secs": 0.028256491,
+    "input_patch_cap": 128,
+    "input_patch_counts": [
+      1,
+      1
+    ],
+    "input_patch_max_bonds": [
+      128,
+      128
+    ],
+    "input_tci_cap": 128,
+    "n_error_samples": 128,
+    "n_gauss": 2048,
+    "patched_contract_max_bond": 128,
+    "r": 10,
+    "r_extra": 0,
+    "rho_max": 8.0,
+    "runs": 2,
+    "sigma": 0.05,
+    "spacing": 3.0,
+    "split_strategy": "sequential_y_first",
+    "warmups": 0
+  },
+  "seed": 0,
+  "tolerance": 1e-8,
+  "wall_time_median_secs": 3.576332009,
+  "wall_times_secs": [
+    3.576332009,
+    3.514710307
+  ],
+  "max_error": 0.0003454729475465829,
+  "input_max_bond_dim": 128,
+  "output_max_bond_dim": 78,
+  "output_bond_dims": [
+    4,
+    16,
+    64,
+    78,
+    30,
+    15,
+    9,
+    6,
+    3
+  ],
+  "n_params": 36336,
+  "rtol": 1e-8,
+  "input_build_secs": 1.185248102
+}
\ No newline at end of file
diff --git a/result/linux-epyc-7713p/raw/mpo_mpo_aniso_patched-global_fit_treetn-n2048-chi192.json b/result/linux-epyc-7713p/raw/mpo_mpo_aniso_patched-global_fit_treetn-n2048-chi192.json
new file mode 100644
index 0000000..03ffe39
--- /dev/null
+++ b/result/linux-epyc-7713p/raw/mpo_mpo_aniso_patched-global_fit_treetn-n2048-chi192.json
@@ -0,0 +1,63 @@
+{
+  "schema_version": 1,
+  "case": "mpo_mpo_aniso_patched",
+  "algorithm": "global_fit_treetn",
+  "params": {
+    "box_l": 3.394112549695429,
+    "box_padding": 1.0,
+    "contraction_svd_policy": "relative_squared_discarded_tail_sum",
+    "error_metric": "max_rel_vs_quantics_grid_reference",
+    "error_seed": 99,
+    "fit_nfullsweeps": 1,
+    "global_contract_max_bond": 192,
+    "inner_l": 3.394112549695429,
+    "input_cache_hit": true,
+    "input_cache_key": "aniso-mpo-v1-n2048-r10-extra0-spacing4008000000000000-padding3ff0000000000000-inner400b27247aff1491-box400b27247aff1491-sigma3fa999999999999a-rho4020000000000000-rtol3e45798ee2308c3a-cap192-seed0",
+    "input_cache_load_secs": 0.031952882,
+    "input_patch_cap": 128,
+    "input_patch_counts": [
+      8,
+      8
+    ],
+    "input_patch_max_bonds": [
+      123,
+      123
+    ],
+    "input_tci_cap": 192,
+    "n_error_samples": 128,
+    "n_gauss": 2048,
+    "patched_contract_max_bond": 128,
+    "r": 10,
+    "r_extra": 0,
+    "rho_max": 8.0,
+    "runs": 2,
+    "sigma": 0.05,
+    "spacing": 3.0,
+    "split_strategy": "sequential_y_first",
+    "warmups": 0
+  },
+  "seed": 0,
+  "tolerance": 1e-8,
+  "wall_time_median_secs": 11.79123273,
+  "wall_times_secs": [
+    11.79123273,
+    11.682349005
+  ],
+  "max_error": 0.0004383987773122374,
+  "input_max_bond_dim": 192,
+  "output_max_bond_dim": 78,
+  "output_bond_dims": [
+    4,
+    16,
+    64,
+    78,
+    31,
+    15,
+    9,
+    6,
+    3
+  ],
+  "n_params": 36708,
+  "rtol": 1e-8,
+  "input_build_secs": 6.63499393
+}
\ No newline at end of file
diff --git a/result/linux-epyc-7713p/raw/mpo_mpo_aniso_patched-global_fit_treetn-n2048-chi224.json b/result/linux-epyc-7713p/raw/mpo_mpo_aniso_patched-global_fit_treetn-n2048-chi224.json
new file mode 100644
index 0000000..db6a2a8
--- /dev/null
+++ b/result/linux-epyc-7713p/raw/mpo_mpo_aniso_patched-global_fit_treetn-n2048-chi224.json
@@ -0,0 +1,63 @@
+{
+  "schema_version": 1,
+  "case": "mpo_mpo_aniso_patched",
+  "algorithm": "global_fit_treetn",
+  "params": {
+    "box_l": 3.394112549695429,
+    "box_padding": 1.0,
+    "contraction_svd_policy": "relative_squared_discarded_tail_sum",
+    "error_metric": "max_rel_vs_quantics_grid_reference",
+    "error_seed": 99,
+    "fit_nfullsweeps": 1,
+    "global_contract_max_bond": 224,
+    "inner_l": 3.394112549695429,
+    "input_cache_hit": true,
+    "input_cache_key": "aniso-mpo-v1-n2048-r10-extra0-spacing4008000000000000-padding3ff0000000000000-inner400b27247aff1491-box400b27247aff1491-sigma3fa999999999999a-rho4020000000000000-rtol3e45798ee2308c3a-cap224-seed0",
+    "input_cache_load_secs": 0.032833929,
+    "input_patch_cap": 128,
+    "input_patch_counts": [
+      8,
+      8
+    ],
+    "input_patch_max_bonds": [
+      121,
+      121
+    ],
+    "input_tci_cap": 224,
+    "n_error_samples": 128,
+    "n_gauss": 2048,
+    "patched_contract_max_bond": 128,
+    "r": 10,
+    "r_extra": 0,
+    "rho_max": 8.0,
+    "runs": 2,
+    "sigma": 0.05,
+    "spacing": 3.0,
+    "split_strategy": "sequential_y_first",
+    "warmups": 0
+  },
+  "seed": 0,
+  "tolerance": 1e-8,
+  "wall_time_median_secs": 15.468127432,
+  "wall_times_secs": [
+    15.347890068,
+    15.468127432
+  ],
+  "max_error": 0.0004367723269018427,
+  "input_max_bond_dim": 224,
+  "output_max_bond_dim": 78,
+  "output_bond_dims": [
+    4,
+    16,
+    64,
+    78,
+    31,
+    15,
+    9,
+    6,
+    3
+  ],
+  "n_params": 36708,
+  "rtol": 1e-8,
+  "input_build_secs": 6.446811312
+}
\ No newline at end of file
diff --git a/result/linux-epyc-7713p/raw/mpo_mpo_aniso_patched-global_fit_treetn-n2048-chi256.json b/result/linux-epyc-7713p/raw/mpo_mpo_aniso_patched-global_fit_treetn-n2048-chi256.json
new file mode 100644
index 0000000..7cb7476
--- /dev/null
+++ b/result/linux-epyc-7713p/raw/mpo_mpo_aniso_patched-global_fit_treetn-n2048-chi256.json
@@ -0,0 +1,63 @@
+{
+  "schema_version": 1,
+  "case": "mpo_mpo_aniso_patched",
+  "algorithm": "global_fit_treetn",
+  "params": {
+    "box_l": 3.394112549695429,
+    "box_padding": 1.0,
+    "contraction_svd_policy": "relative_squared_discarded_tail_sum",
+    "error_metric": "max_rel_vs_quantics_grid_reference",
+    "error_seed": 99,
+    "fit_nfullsweeps": 1,
+    "global_contract_max_bond": 256,
+    "inner_l": 3.394112549695429,
+    "input_cache_hit": true,
+    "input_cache_key": "aniso-mpo-v1-n2048-r10-extra0-spacing4008000000000000-padding3ff0000000000000-inner400b27247aff1491-box400b27247aff1491-sigma3fa999999999999a-rho4020000000000000-rtol3e45798ee2308c3a-cap384-seed0",
+    "input_cache_load_secs": 0.034797543,
+    "input_patch_cap": 128,
+    "input_patch_counts": [
+      8,
+      8
+    ],
+    "input_patch_max_bonds": [
+      119,
+      120
+    ],
+    "input_tci_cap": 384,
+    "n_error_samples": 128,
+    "n_gauss": 2048,
+    "patched_contract_max_bond": 128,
+    "r": 10,
+    "r_extra": 0,
+    "rho_max": 8.0,
+    "runs": 2,
+    "sigma": 0.05,
+    "spacing": 3.0,
+    "split_strategy": "sequential_y_first",
+    "warmups": 0
+  },
+  "seed": 0,
+  "tolerance": 1e-8,
+  "wall_time_median_secs": 19.606482696,
+  "wall_times_secs": [
+    19.606482696,
+    19.557511476
+  ],
+  "max_error": 0.0004394615152723373,
+  "input_max_bond_dim": 256,
+  "output_max_bond_dim": 78,
+  "output_bond_dims": [
+    4,
+    16,
+    64,
+    78,
+    31,
+    15,
+    9,
+    6,
+    3
+  ],
+  "n_params": 36708,
+  "rtol": 1e-8,
+  "input_build_secs": 6.272382668
+}
\ No newline at end of file
diff --git a/result/linux-epyc-7713p/raw/mpo_mpo_aniso_patched-patched_fit_treetn-n2048-chi128.json b/result/linux-epyc-7713p/raw/mpo_mpo_aniso_patched-patched_fit_treetn-n2048-chi128.json
new file mode 100644
index 0000000..a347dcb
--- /dev/null
+++ b/result/linux-epyc-7713p/raw/mpo_mpo_aniso_patched-patched_fit_treetn-n2048-chi128.json
@@ -0,0 +1,64 @@
+{
+  "schema_version": 1,
+  "case": "mpo_mpo_aniso_patched",
+  "algorithm": "patched_fit_treetn",
+  "params": {
+    "box_l": 3.394112549695429,
+    "box_padding": 1.0,
+    "contraction_svd_policy": "relative_squared_discarded_tail_sum",
+    "error_metric": "max_rel_vs_quantics_grid_reference",
+    "error_seed": 99,
+    "fit_nfullsweeps": 1,
+    "global_contract_max_bond": 128,
+    "inner_l": 3.394112549695429,
+    "input_cache_hit": true,
+    "input_cache_key": "aniso-mpo-v1-n2048-r10-extra0-spacing4008000000000000-padding3ff0000000000000-inner400b27247aff1491-box400b27247aff1491-sigma3fa999999999999a-rho4020000000000000-rtol3e45798ee2308c3a-cap128-seed0",
+    "input_cache_load_secs": 0.028256491,
+    "input_patch_cap": 128,
+    "input_patch_counts": [
+      1,
+      1
+    ],
+    "input_patch_max_bonds": [
+      128,
+      128
+    ],
+    "input_tci_cap": 128,
+    "n_error_samples": 128,
+    "n_gauss": 2048,
+    "patched_contract_max_bond": 128,
+    "r": 10,
+    "r_extra": 0,
+    "rho_max": 8.0,
+    "runs": 2,
+    "sigma": 0.05,
+    "spacing": 3.0,
+    "split_strategy": "sequential_y_first",
+    "warmups": 0
+  },
+  "seed": 0,
+  "tolerance": 1e-8,
+  "wall_time_median_secs": 3.321234724,
+  "wall_times_secs": [
+    3.321234724,
+    3.292352661
+  ],
+  "max_error": 0.0003454900427644558,
+  "input_max_bond_dim": 128,
+  "output_max_bond_dim": 78,
+  "output_bond_dims": [
+    4,
+    16,
+    64,
+    78,
+    30,
+    15,
+    9,
+    6,
+    3
+  ],
+  "n_patches": 1,
+  "max_patch_bond": 78,
+  "rtol": 1e-8,
+  "input_build_secs": 1.185248102
+}
\ No newline at end of file
diff --git a/result/linux-epyc-7713p/raw/mpo_mpo_aniso_patched-patched_fit_treetn-n2048-chi192.json b/result/linux-epyc-7713p/raw/mpo_mpo_aniso_patched-patched_fit_treetn-n2048-chi192.json
new file mode 100644
index 0000000..22de149
--- /dev/null
+++ b/result/linux-epyc-7713p/raw/mpo_mpo_aniso_patched-patched_fit_treetn-n2048-chi192.json
@@ -0,0 +1,64 @@
+{
+  "schema_version": 1,
+  "case": "mpo_mpo_aniso_patched",
+  "algorithm": "patched_fit_treetn",
+  "params": {
+    "box_l": 3.394112549695429,
+    "box_padding": 1.0,
+    "contraction_svd_policy": "relative_squared_discarded_tail_sum",
+    "error_metric": "max_rel_vs_quantics_grid_reference",
+    "error_seed": 99,
+    "fit_nfullsweeps": 1,
+    "global_contract_max_bond": 192,
+    "inner_l": 3.394112549695429,
+    "input_cache_hit": true,
+    "input_cache_key": "aniso-mpo-v1-n2048-r10-extra0-spacing4008000000000000-padding3ff0000000000000-inner400b27247aff1491-box400b27247aff1491-sigma3fa999999999999a-rho4020000000000000-rtol3e45798ee2308c3a-cap192-seed0",
+    "input_cache_load_secs": 0.031952882,
+    "input_patch_cap": 128,
+    "input_patch_counts": [
+      8,
+      8
+    ],
+    "input_patch_max_bonds": [
+      123,
+      123
+    ],
+    "input_tci_cap": 192,
+    "n_error_samples": 128,
+    "n_gauss": 2048,
+    "patched_contract_max_bond": 128,
+    "r": 10,
+    "r_extra": 0,
+    "rho_max": 8.0,
+    "runs": 2,
+    "sigma": 0.05,
+    "spacing": 3.0,
+    "split_strategy": "sequential_y_first",
+    "warmups": 0
+  },
+  "seed": 0,
+  "tolerance": 1e-8,
+  "wall_time_median_secs": 13.352732666,
+  "wall_times_secs": [
+    13.352732666,
+    13.243957682
+  ],
+  "max_error": 0.0005547409992032523,
+  "input_max_bond_dim": 192,
+  "output_max_bond_dim": 68,
+  "output_bond_dims": [
+    4,
+    16,
+    64,
+    68,
+    28,
+    15,
+    8,
+    6,
+    3
+  ],
+  "n_patches": 1,
+  "max_patch_bond": 68,
+  "rtol": 1e-8,
+  "input_build_secs": 6.63499393
+}
\ No newline at end of file
diff --git a/result/linux-epyc-7713p/raw/mpo_mpo_aniso_patched-patched_fit_treetn-n2048-chi224.json b/result/linux-epyc-7713p/raw/mpo_mpo_aniso_patched-patched_fit_treetn-n2048-chi224.json
new file mode 100644
index 0000000..7f2e047
--- /dev/null
+++ b/result/linux-epyc-7713p/raw/mpo_mpo_aniso_patched-patched_fit_treetn-n2048-chi224.json
@@ -0,0 +1,64 @@
+{
+  "schema_version": 1,
+  "case": "mpo_mpo_aniso_patched",
+  "algorithm": "patched_fit_treetn",
+  "params": {
+    "box_l": 3.394112549695429,
+    "box_padding": 1.0,
+    "contraction_svd_policy": "relative_squared_discarded_tail_sum",
+    "error_metric": "max_rel_vs_quantics_grid_reference",
+    "error_seed": 99,
+    "fit_nfullsweeps": 1,
+    "global_contract_max_bond": 224,
+    "inner_l": 3.394112549695429,
+    "input_cache_hit": true,
+    "input_cache_key": "aniso-mpo-v1-n2048-r10-extra0-spacing4008000000000000-padding3ff0000000000000-inner400b27247aff1491-box400b27247aff1491-sigma3fa999999999999a-rho4020000000000000-rtol3e45798ee2308c3a-cap224-seed0",
+    "input_cache_load_secs": 0.032833929,
+    "input_patch_cap": 128,
+    "input_patch_counts": [
+      8,
+      8
+    ],
+    "input_patch_max_bonds": [
+      121,
+      121
+    ],
+    "input_tci_cap": 224,
+    "n_error_samples": 128,
+    "n_gauss": 2048,
+    "patched_contract_max_bond": 128,
+    "r": 10,
+    "r_extra": 0,
+    "rho_max": 8.0,
+    "runs": 2,
+    "sigma": 0.05,
+    "spacing": 3.0,
+    "split_strategy": "sequential_y_first",
+    "warmups": 0
+  },
+  "seed": 0,
+  "tolerance": 1e-8,
+  "wall_time_median_secs": 13.102732971,
+  "wall_times_secs": [
+    13.073508836,
+    13.102732971
+  ],
+  "max_error": 0.0004916712233304541,
+  "input_max_bond_dim": 224,
+  "output_max_bond_dim": 68,
+  "output_bond_dims": [
+    4,
+    16,
+    64,
+    68,
+    28,
+    15,
+    8,
+    6,
+    3
+  ],
+  "n_patches": 1,
+  "max_patch_bond": 68,
+  "rtol": 1e-8,
+  "input_build_secs": 6.446811312
+}
\ No newline at end of file
diff --git a/result/linux-epyc-7713p/raw/mpo_mpo_aniso_patched-patched_fit_treetn-n2048-chi256.json b/result/linux-epyc-7713p/raw/mpo_mpo_aniso_patched-patched_fit_treetn-n2048-chi256.json
new file mode 100644
index 0000000..4303bd0
--- /dev/null
+++ b/result/linux-epyc-7713p/raw/mpo_mpo_aniso_patched-patched_fit_treetn-n2048-chi256.json
@@ -0,0 +1,64 @@
+{
+  "schema_version": 1,
+  "case": "mpo_mpo_aniso_patched",
+  "algorithm": "patched_fit_treetn",
+  "params": {
+    "box_l": 3.394112549695429,
+    "box_padding": 1.0,
+    "contraction_svd_policy": "relative_squared_discarded_tail_sum",
+    "error_metric": "max_rel_vs_quantics_grid_reference",
+    "error_seed": 99,
+    "fit_nfullsweeps": 1,
+    "global_contract_max_bond": 256,
+    "inner_l": 3.394112549695429,
+    "input_cache_hit": true,
+    "input_cache_key": "aniso-mpo-v1-n2048-r10-extra0-spacing4008000000000000-padding3ff0000000000000-inner400b27247aff1491-box400b27247aff1491-sigma3fa999999999999a-rho4020000000000000-rtol3e45798ee2308c3a-cap384-seed0",
+    "input_cache_load_secs": 0.034797543,
+    "input_patch_cap": 128,
+    "input_patch_counts": [
+      8,
+      8
+    ],
+    "input_patch_max_bonds": [
+      119,
+      120
+    ],
+    "input_tci_cap": 384,
+    "n_error_samples": 128,
+    "n_gauss": 2048,
+    "patched_contract_max_bond": 128,
+    "r": 10,
+    "r_extra": 0,
+    "rho_max": 8.0,
+    "runs": 2,
+    "sigma": 0.05,
+    "spacing": 3.0,
+    "split_strategy": "sequential_y_first",
+    "warmups": 0
+  },
+  "seed": 0,
+  "tolerance": 1e-8,
+  "wall_time_median_secs": 12.273990265,
+  "wall_times_secs": [
+    12.273990265,
+    12.095820478
+  ],
+  "max_error": 0.0006281323596457093,
+  "input_max_bond_dim": 256,
+  "output_max_bond_dim": 68,
+  "output_bond_dims": [
+    4,
+    16,
+    64,
+    68,
+    28,
+    15,
+    8,
+    6,
+    3
+  ],
+  "n_patches": 1,
+  "max_patch_bond": 68,
+  "rtol": 1e-8,
+  "input_build_secs": 6.272382668
+}
\ No newline at end of file
diff --git a/result/linux-epyc-7713p/run.yaml b/result/linux-epyc-7713p/run.yaml
new file mode 100644
index 0000000..4295b75
--- /dev/null
+++ b/result/linux-epyc-7713p/run.yaml
@@ -0,0 +1,13 @@
+profile: linux-epyc-7713p
+date: 2026-08-19T10:30:32Z
+machine: linux-epyc-7713p
+os: Linux x86_64
+chip: AMD EPYC 7713P 64-Core Processor
+memory_gb: 995
+repo_rev: add289e57c9b090a10d048892746d8e396255682
+working_tree: verify/partitioned-mpo-correctness
+tensor4all_rs_rev: 9e9aedaebe0d3918b34dd399ff0981e337f3835b
+threads: 1
+cpu_affinity: 62
+runs: 2
+warmups: 0
diff --git a/scripts/report.py b/scripts/report.py
index 55d2096..5c60339 100644
--- a/scripts/report.py
+++ b/scripts/report.py
@@ -30,6 +30,11 @@
     "elementwise_gauss2d_patched": "N",
 }
 
+# Case 6 has an intentionally hand-maintained crossover table because its x
+# points are selected by input TCI caps and every committed run records a failed
+# accuracy gate. Do not replace that context with the generic scaling report.
+HAND_MAINTAINED_CASES = {"mpo_mpo_aniso_patched"}
+
 ERROR_LABEL = {
     # case name -> what `max_error` in the records actually measures
     "elementwise_fourier": "max abs error",
@@ -444,6 +449,12 @@ def main():
     if not cases:
         sys.exit(f"no records under {profile_dir}/raw")
     for case, algos in cases.items():
+        if case in HAND_MAINTAINED_CASES:
+            report = profile_dir / f"{case}.md"
+            if not report.is_file():
+                sys.exit(f"missing hand-maintained report {report}")
+            print(f"kept hand-maintained {report}")
+            continue
         render_case(case, algos, profile_dir)
 
 
diff --git a/src/bin/mpo_mpo_aniso_patched.rs b/src/bin/mpo_mpo_aniso_patched.rs
new file mode 100644
index 0000000..c19b8e3
--- /dev/null
+++ b/src/bin/mpo_mpo_aniso_patched.rs
@@ -0,0 +1,321 @@
+//! Global and patched chain-TreeTN fit contraction of anisotropic Gaussian MPOs.
+
+use std::fs;
+use std::path::PathBuf;
+use std::time::Instant;
+
+use tensor4all_hdf5::{append_mps, load_mps, save_mps};
+
+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::record::{write_record, RunRecord, SCHEMA_VERSION};
+
+fn env_or<T: std::str::FromStr>(key: &str, default: T) -> T {
+    std::env::var(key)
+        .ok()
+        .and_then(|value| value.parse().ok())
+        .unwrap_or(default)
+}
+
+fn main() -> anyhow::Result<()> {
+    let n_gauss: usize = env_or("BENCH_NGAUSS", 512);
+    let sigma: f64 = env_or("BENCH_ANISO_SIGMA", 0.05);
+    let rho_max: f64 = env_or("BENCH_ANISO_RHO_MAX", 8.0);
+    let spacing: f64 = env_or("BENCH_ANISO_SPACING", 3.0);
+    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_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);
+    let runs: usize = env_or("BENCH_RUNS", 1);
+    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 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()));
+    let cache_dir = PathBuf::from(
+        std::env::var("BENCH_INPUT_CACHE_DIR").unwrap_or_else(|_| ".cache/inputs".into()),
+    );
+    let refresh_cache = std::env::var_os("BENCH_INPUT_CACHE_REFRESH").is_some();
+
+    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_tci_cap > 0 && patch_cap > 0 && runs > 0);
+
+    // Same constant-density anisotropic family and resolution as case 5.
+    let inner_l = spacing * sigma * (n_gauss as f64).sqrt() / 2.0;
+    let box_l = box_padding * inner_l;
+    let r = (8.0 * box_l / sigma).log2().ceil().max(1.0) as usize + extra_bits;
+    let f = AnisoMixture2D::random(n_gauss, inner_l, sigma, rho_max, seed.wrapping_add(1));
+    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_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 {
+        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,
+        )
+    } 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)?;
+        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)?;
+            fs::rename(&temporary, &cache_path)?;
+            Ok(())
+        })();
+        if write_result.is_err() {
+            let _ = fs::remove_file(&temporary);
+        }
+        write_result?;
+        (patched, 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}"
+    );
+    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"
+    );
+    eprintln!(
+        "prepared input_chi={input_chi} patches={input_patches:?} patch_chi={input_patch_bonds:?} build={input_build_secs:.3}s"
+    );
+
+    let (global_train, global_timing) = time_median(warmups, runs, || {
+        patched
+            .contract_fit_global(rtol, input_chi)
+            .expect("global fit contraction failed")
+    });
+    let (tree_partitioned, tree_timing) = time_median(warmups, runs, || {
+        patched
+            .contract_fit_treetn_partitioned(rtol, patch_cap)
+            .expect("PartitionedTreeTN fit contraction failed")
+    });
+    let global_output = patched.finish_global_output(&global_train)?;
+    let tree_output = patched.finish_treetn_output(tree_partitioned)?;
+
+    let error_seed = seed.wrapping_add(99);
+    let global_error = max_rel_error_vs_aniso_grid(
+        &global_output,
+        dy,
+        &f,
+        &g,
+        r,
+        box_l,
+        n_error_samples,
+        error_seed,
+    );
+    let tree_error = max_rel_error_vs_aniso_grid(
+        &tree_output.mpo,
+        dy,
+        &f,
+        &g,
+        r,
+        box_l,
+        n_error_samples,
+        error_seed,
+    );
+    let common_params = serde_json::json!({
+        "n_gauss": n_gauss,
+        "r": r,
+        "inner_l": inner_l,
+        "box_l": box_l,
+        "box_padding": box_padding,
+        "sigma": sigma,
+        "rho_max": rho_max,
+        "spacing": spacing,
+        "r_extra": extra_bits,
+        "input_tci_cap": input_tci_cap,
+        "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],
+        "global_contract_max_bond": input_chi,
+        "patched_contract_max_bond": patch_cap,
+        "fit_nfullsweeps": 1,
+        "contraction_svd_policy": "relative_squared_discarded_tail_sum",
+        "split_strategy": "sequential_y_first",
+        "runs": runs,
+        "warmups": warmups,
+        "n_error_samples": n_error_samples,
+        "error_seed": error_seed,
+        "input_cache_key": cache_key,
+        "input_cache_hit": cache_hit,
+        "input_cache_load_secs": input_cache_load_secs,
+        "error_metric": "max_rel_vs_quantics_grid_reference",
+    });
+    let global_record = RunRecord {
+        schema_version: SCHEMA_VERSION,
+        case: "mpo_mpo_aniso_patched".into(),
+        algorithm: "global_fit_treetn".into(),
+        params: common_params.clone(),
+        seed,
+        tolerance: rtol,
+        wall_time_median_secs: global_timing.median_secs,
+        wall_times_secs: global_timing.runs_secs.clone(),
+        max_error: global_error,
+        input_max_bond_dim: input_chi,
+        output_max_bond_dim: global_output.rank(),
+        output_bond_dims: global_output.link_dims(),
+        n_params: Some(mpo_n_params(&global_output)),
+        n_patches: None,
+        max_patch_bond: None,
+        rtol: Some(rtol),
+        input_build_secs: Some(input_build_secs),
+    };
+    write_record(
+        &out_dir,
+        &format!("mpo_mpo_aniso_patched-global_fit_treetn-n{n_gauss}-chi{input_chi}"),
+        &global_record,
+    )?;
+    write_arm(
+        &out_dir,
+        "patched_fit_treetn",
+        &tree_output,
+        &tree_timing,
+        tree_error,
+        input_chi,
+        input_build_secs,
+        seed,
+        rtol,
+        &common_params,
+    )?;
+
+    eprintln!(
+        "global_fit_treetn: t={:.6}s error={global_error:.3e} output_chi={}",
+        global_timing.median_secs,
+        global_output.rank()
+    );
+    eprintln!(
+        "patched_fit_treetn: t={:.6}s error={tree_error:.3e} output_patches={} output_chi={}",
+        tree_timing.median_secs, tree_output.n_patches, tree_output.max_patch_bond
+    );
+    anyhow::ensure!(
+        global_error <= sanity,
+        "global fit relative error {global_error:.3e} > {sanity:.3e}"
+    );
+    anyhow::ensure!(
+        tree_error <= sanity,
+        "TreeTN relative error {tree_error:.3e} > {sanity:.3e}"
+    );
+    Ok(())
+}
+
+#[allow(clippy::too_many_arguments)]
+fn write_arm(
+    out_dir: &std::path::Path,
+    algorithm: &str,
+    output: &PatchedMpoOutput,
+    timing: &t4a_bench::harness::Timing,
+    max_error: f64,
+    input_chi: usize,
+    input_build_secs: f64,
+    seed: u64,
+    rtol: f64,
+    params: &serde_json::Value,
+) -> anyhow::Result<()> {
+    let record = RunRecord {
+        schema_version: SCHEMA_VERSION,
+        case: "mpo_mpo_aniso_patched".into(),
+        algorithm: algorithm.into(),
+        params: params.clone(),
+        seed,
+        tolerance: rtol,
+        wall_time_median_secs: timing.median_secs,
+        wall_times_secs: timing.runs_secs.clone(),
+        max_error,
+        input_max_bond_dim: input_chi,
+        output_max_bond_dim: output.mpo.rank(),
+        output_bond_dims: output.mpo.link_dims(),
+        n_params: None,
+        n_patches: Some(output.n_patches),
+        max_patch_bond: Some(output.max_patch_bond),
+        rtol: Some(rtol),
+        input_build_secs: Some(input_build_secs),
+    };
+    write_record(
+        out_dir,
+        &format!(
+            "mpo_mpo_aniso_patched-{algorithm}-n{}-chi{input_chi}",
+            params["n_gauss"]
+        ),
+        &record,
+    )
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use std::time::{SystemTime, UNIX_EPOCH};
+
+    #[test]
+    fn cached_inputs_round_trip_with_shared_index_identity() {
+        let field = AnisoMixture2D::random(2, 1.0, 0.2, 3.0, 7);
+        let other = AnisoMixture2D::random(2, 1.0, 0.2, 3.0, 8);
+        let (left, _) = to_quantics_mpo_field(&field, 4, 1.0, 1.0e-8, 32).unwrap();
+        let (right, _) = to_quantics_mpo_field(&other, 4, 1.0, 1.0e-8, 32).unwrap();
+        let original = PatchedMpoPair::new(&left, &right, 1.0e-8, 16).unwrap();
+        let nonce = SystemTime::now()
+            .duration_since(UNIX_EPOCH)
+            .unwrap()
+            .as_nanos();
+        let path =
+            std::env::temp_dir().join(format!("t4a-mpo-cache-{}-{nonce}.h5", std::process::id()));
+        let path_str = path.to_string_lossy();
+        let (left_train, right_train) = original.global_inputs();
+        save_mps(&path_str, "left", left_train).unwrap();
+        append_mps(&path_str, "right", right_train).unwrap();
+
+        let cached = PatchedMpoPair::from_tensortrains(
+            load_mps(&path_str, "left").unwrap(),
+            load_mps(&path_str, "right").unwrap(),
+            1.0e-8,
+            16,
+        )
+        .unwrap();
+        fs::remove_file(path).unwrap();
+
+        assert_eq!(cached.input_max_bond(), original.input_max_bond());
+        assert_eq!(cached.input_patch_counts(), original.input_patch_counts());
+        assert_eq!(
+            cached.input_patch_max_bonds(),
+            original.input_patch_max_bonds()
+        );
+    }
+}
diff --git a/src/gaussian.rs b/src/gaussian.rs
index 2237dba..e7f78b9 100644
--- a/src/gaussian.rs
+++ b/src/gaussian.rs
@@ -160,6 +160,20 @@ impl AnisoMixture2D {
             })
             .sum()
     }
+
+    /// First derivatives with respect to the first and second coordinates.
+    pub fn gradient(&self, x: f64, y: f64) -> (f64, f64) {
+        (0..self.weights.len()).fold((0.0, 0.0), |(gx, gy), i| {
+            let (cx, cy) = self.centers[i];
+            let (a, b, c) = self.quad[i];
+            let (dx, dy) = (x - cx, y - cy);
+            let value = self.weights[i] * (-(a * dx * dx + 2.0 * b * dx * dy + c * dy * dy)).exp();
+            (
+                gx - 2.0 * (a * dx + b * dy) * value,
+                gy - 2.0 * (b * dx + c * dy) * value,
+            )
+        })
+    }
 }
 
 impl Field2D for AnisoMixture2D {
@@ -188,6 +202,83 @@ pub fn analytic_contraction(f: &GaussianMixture2D, g: &GaussianMixture2D, x: f64
     s
 }
 
+/// Closed form of `int f(x,y) g(y,z) dy` for anisotropic mixtures.
+pub fn analytic_contraction_aniso(f: &AnisoMixture2D, g: &AnisoMixture2D, x: f64, z: f64) -> f64 {
+    analytic_contraction_aniso_with(f, g, x, z, |quadratic, linear, constant| {
+        (std::f64::consts::PI / quadratic).sqrt()
+            * (-(constant - linear * linear / quadratic)).exp()
+    })
+}
+
+/// Closed form of `int f(x,y) g(y,z) dy` over `[-box_l, box_l]`.
+pub fn analytic_contraction_aniso_box(
+    f: &AnisoMixture2D,
+    g: &AnisoMixture2D,
+    x: f64,
+    z: f64,
+    box_l: f64,
+) -> f64 {
+    analytic_contraction_aniso_with(f, g, x, z, |quadratic, linear, constant| {
+        let root = quadratic.sqrt();
+        let center = -linear / quadratic;
+        let bounded = libm::erf(root * (box_l - center)) - libm::erf(root * (-box_l - center));
+        0.5 * (std::f64::consts::PI / quadratic).sqrt()
+            * (-(constant - linear * linear / quadratic)).exp()
+            * bounded
+    })
+}
+
+/// Euler-Maclaurin reference for the left-endpoint quantics grid sum times its step.
+pub fn discrete_contraction_aniso_reference(
+    f: &AnisoMixture2D,
+    g: &AnisoMixture2D,
+    x: f64,
+    z: f64,
+    r: usize,
+    box_l: f64,
+) -> f64 {
+    let step = 2.0 * box_l / (1u64 << r) as f64;
+    let boundary = |y: f64| {
+        let fv = f.eval(x, y);
+        let gv = g.eval(y, z);
+        let fy = f.gradient(x, y).1;
+        let gy = g.gradient(y, z).0;
+        (fv * gv, fy * gv + fv * gy)
+    };
+    let (left, left_derivative) = boundary(-box_l);
+    let (right, right_derivative) = boundary(box_l);
+    analytic_contraction_aniso_box(f, g, x, z, box_l)
+        + 0.5 * step * (left - right)
+        + step * step * (right_derivative - left_derivative) / 12.0
+}
+
+fn analytic_contraction_aniso_with(
+    f: &AnisoMixture2D,
+    g: &AnisoMixture2D,
+    x: f64,
+    z: f64,
+    integral: impl Fn(f64, f64, f64) -> f64,
+) -> f64 {
+    let mut sum = 0.0;
+    for i in 0..f.weights.len() {
+        let (fcx, fcy) = f.centers[i];
+        let (fa, fb, fc) = f.quad[i];
+        let dx = x - fcx;
+        for j in 0..g.weights.len() {
+            let (gcy, gcz) = g.centers[j];
+            let (ga, gb, gc) = g.quad[j];
+            let dz = z - gcz;
+            let quadratic = fc + ga;
+            let linear = fb * dx - fc * fcy + gb * dz - ga * gcy;
+            let constant = fa * dx * dx - 2.0 * fb * dx * fcy + fc * fcy * fcy + ga * gcy * gcy
+                - 2.0 * gb * gcy * dz
+                + gc * dz * dz;
+            sum += f.weights[i] * g.weights[j] * integral(quadratic, linear, constant);
+        }
+    }
+    sum
+}
+
 /// Coordinate of grid index `i`: `-L + i * 2L/2^R` (the grid is `[-L, L)`).
 pub fn grid_coord(i: u64, r: usize, box_l: f64) -> f64 {
     let step = 2.0 * box_l / (1u64 << r) as f64;
@@ -242,7 +333,18 @@ pub fn to_quantics_mpo(
     tol: f64,
     max_bond: usize,
 ) -> anyhow::Result<(MPO<f64>, f64)> {
-    let (tt, step) = to_quantics_fused_tt(mix, r, box_l, tol, max_bond)?;
+    to_quantics_mpo_field(mix, r, box_l, tol, max_bond)
+}
+
+/// [`to_quantics_mpo`] for any [`Field2D`].
+pub fn to_quantics_mpo_field<M: Field2D>(
+    mix: &M,
+    r: usize,
+    box_l: f64,
+    tol: f64,
+    max_bond: usize,
+) -> anyhow::Result<(MPO<f64>, 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());
@@ -273,26 +375,60 @@ pub fn to_quantics_mpo(
 mod tests {
     use super::*;
 
+    fn quadrature(f: impl Fn(f64) -> f64, lo: f64, hi: f64, intervals: usize) -> f64 {
+        let step = (hi - lo) / intervals as f64;
+        (0..=intervals)
+            .map(|i| {
+                let weight = if i == 0 || i == intervals { 0.5 } else { 1.0 };
+                weight * f(lo + i as f64 * step) * step
+            })
+            .sum()
+    }
+
     #[test]
     fn analytic_contraction_matches_quadrature() {
         let f = GaussianMixture2D::random(3, 4.0, (0.5, 2.0), 1);
         let g = GaussianMixture2D::random(2, 4.0, (0.5, 2.0), 2);
         let (x, z) = (0.3, -0.7);
-        // trapezoid quadrature over y in [-8, 8]
-        let n = 20_000;
-        let (lo, hi) = (-8.0, 8.0);
-        let h = (hi - lo) / n as f64;
-        let mut s = 0.0;
-        for i in 0..=n {
-            let y = lo + i as f64 * h;
-            let w = if i == 0 || i == n { 0.5 } else { 1.0 };
-            s += w * f.eval(x, y) * g.eval(y, z) * h;
-        }
-        let a = analytic_contraction(&f, &g, x, z);
-        assert!(
-            (s - a).abs() < 1e-8 * a.abs().max(1.0),
-            "quad {s} vs analytic {a}"
-        );
+        let numeric = quadrature(|y| f.eval(x, y) * g.eval(y, z), -8.0, 8.0, 20_000);
+        let analytic = analytic_contraction(&f, &g, x, z);
+        assert!((numeric - analytic).abs() < 1e-8 * analytic.abs().max(1.0));
+    }
+
+    #[test]
+    fn anisotropic_analytic_contraction_matches_quadrature() {
+        let f = AnisoMixture2D::random(3, 4.0, 0.3, 4.0, 1);
+        let g = AnisoMixture2D::random(2, 4.0, 0.3, 4.0, 2);
+        let (x, z) = (0.3, -0.7);
+        let numeric = quadrature(|y| f.eval(x, y) * g.eval(y, z), -12.0, 12.0, 40_000);
+        let analytic = analytic_contraction_aniso(&f, &g, x, z);
+        assert!((numeric - analytic).abs() < 1e-8 * analytic.abs().max(1.0));
+    }
+
+    #[test]
+    fn bounded_anisotropic_contraction_matches_quadrature() {
+        let f = AnisoMixture2D::random(3, 1.0, 0.3, 4.0, 1);
+        let g = AnisoMixture2D::random(2, 1.0, 0.3, 4.0, 2);
+        let (x, z, box_l) = (0.3, -0.7, 1.0);
+        let numeric = quadrature(|y| f.eval(x, y) * g.eval(y, z), -box_l, box_l, 40_000);
+        let analytic = analytic_contraction_aniso_box(&f, &g, x, z, box_l);
+        assert!((numeric - analytic).abs() < 1e-8 * analytic.abs().max(1.0));
+    }
+
+    #[test]
+    fn anisotropic_grid_reference_matches_the_explicit_sum() {
+        let f = AnisoMixture2D::random(3, 1.0, 0.3, 4.0, 1);
+        let g = AnisoMixture2D::random(2, 1.0, 0.3, 4.0, 2);
+        let (x, z, r, box_l) = (0.3, -0.7, 10, 1.0);
+        let step = 2.0 * box_l / (1u64 << r) as f64;
+        let direct = (0..1u64 << r)
+            .map(|i| {
+                let y = -box_l + i as f64 * step;
+                f.eval(x, y) * g.eval(y, z) * step
+            })
+            .sum::<f64>();
+        let reference = discrete_contraction_aniso_reference(&f, &g, x, z, r, box_l);
+        assert!((direct - reference).abs() < 1e-10 * direct.abs().max(1.0));
     }
 
     /// The rotated quadratic form, against a value computed by hand.
diff --git a/src/lib.rs b/src/lib.rs
index a3c3c9d..fe8401a 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -5,5 +5,6 @@ pub mod harness;
 pub mod hdf5_export;
 pub mod mpo_contract;
 pub mod patched;
+pub mod patched_mpo;
 pub mod record;
 pub mod scalar;
diff --git a/src/mpo_contract.rs b/src/mpo_contract.rs
index 467baff..6470121 100644
--- a/src/mpo_contract.rs
+++ b/src/mpo_contract.rs
@@ -10,7 +10,10 @@ use tensor4all_simplett::mpo::{
     contract_naive, contract_zipup, tensor4_from_data, ContractionOptions, Tensor4Ops, MPO,
 };
 
-use crate::gaussian::{analytic_contraction, grid_coord, GaussianMixture2D};
+use crate::gaussian::{
+    analytic_contraction, analytic_contraction_aniso, analytic_contraction_aniso_box,
+    discrete_contraction_aniso_reference, grid_coord, AnisoMixture2D, GaussianMixture2D,
+};
 use crate::harness::{index_to_bits, sample_grid_indices};
 
 /// Number of full variational sweeps used by the `fit_treetn` arm.
@@ -179,7 +182,7 @@ fn contract_via_bridge(
 /// adjacent sites to share exactly one index and identifies links purely by
 /// shared index id, and a leading or trailing extent-1 axis does not change
 /// the column-major layout.
-fn mpo_to_tensortrain(
+pub fn mpo_to_tensortrain(
     m: &MPO<f64>,
     site1: &[DynIndex],
     site2: &[DynIndex],
@@ -220,7 +223,7 @@ fn mpo_to_tensortrain(
 /// The contraction preserves the external site index objects, so each result
 /// core is permuted into `(left, site1, site2, right)` and read out as
 /// column-major data, which is what `tensor4_from_data` expects.
-fn tensortrain_to_mpo(
+pub fn tensortrain_to_mpo(
     tt: &TensorTrain,
     site1: &[DynIndex],
     site2: &[DynIndex],
@@ -271,6 +274,71 @@ pub fn max_rel_error_vs_analytic(
     box_l: f64,
     n_samples: usize,
     seed: u64,
+) -> f64 {
+    max_rel_error_vs_reference(h, dy, r, box_l, n_samples, seed, |x, z| {
+        analytic_contraction(f, g, x, z)
+    })
+}
+
+/// Relative sampled error against the anisotropic closed-form contraction.
+#[allow(clippy::too_many_arguments)]
+pub fn max_rel_error_vs_aniso_analytic(
+    h: &MPO<f64>,
+    dy: f64,
+    f: &AnisoMixture2D,
+    g: &AnisoMixture2D,
+    r: usize,
+    box_l: f64,
+    n_samples: usize,
+    seed: u64,
+) -> f64 {
+    max_rel_error_vs_reference(h, dy, r, box_l, n_samples, seed, |x, z| {
+        analytic_contraction_aniso(f, g, x, z)
+    })
+}
+
+/// Relative sampled error against the finite-box anisotropic contraction.
+#[allow(clippy::too_many_arguments)]
+pub fn max_rel_error_vs_aniso_box(
+    h: &MPO<f64>,
+    dy: f64,
+    f: &AnisoMixture2D,
+    g: &AnisoMixture2D,
+    r: usize,
+    box_l: f64,
+    n_samples: usize,
+    seed: u64,
+) -> f64 {
+    max_rel_error_vs_reference(h, dy, r, box_l, n_samples, seed, |x, z| {
+        analytic_contraction_aniso_box(f, g, x, z, box_l)
+    })
+}
+
+/// Relative sampled error against the quantics grid-sum reference.
+#[allow(clippy::too_many_arguments)]
+pub fn max_rel_error_vs_aniso_grid(
+    h: &MPO<f64>,
+    dy: f64,
+    f: &AnisoMixture2D,
+    g: &AnisoMixture2D,
+    r: usize,
+    box_l: f64,
+    n_samples: usize,
+    seed: u64,
+) -> f64 {
+    max_rel_error_vs_reference(h, dy, r, box_l, n_samples, seed, |x, z| {
+        discrete_contraction_aniso_reference(f, g, x, z, r, box_l)
+    })
+}
+
+fn max_rel_error_vs_reference(
+    h: &MPO<f64>,
+    dy: f64,
+    r: usize,
+    box_l: f64,
+    n_samples: usize,
+    seed: u64,
+    reference: impl Fn(f64, f64) -> f64,
 ) -> f64 {
     let xs = sample_grid_indices(r, n_samples, seed);
     let zs = sample_grid_indices(r, n_samples, seed.wrapping_add(1));
@@ -287,13 +355,62 @@ pub fn max_rel_error_vs_analytic(
             idx.push(zb[n]);
         }
         let got = h.evaluate(&idx).unwrap() * dy;
-        let want = analytic_contraction(f, g, x, z);
+        let want = reference(x, z);
         max_abs = max_abs.max((got - want).abs());
         max_ref = max_ref.max(want.abs());
     }
     max_abs / max_ref.max(f64::MIN_POSITIVE)
 }
 
+/// Maximum absolute difference between two MPOs at sampled `(x, z)` grid points.
+pub fn max_sampled_mpo_diff(
+    left: &MPO<f64>,
+    right: &MPO<f64>,
+    r: usize,
+    n_samples: usize,
+    seed: u64,
+) -> f64 {
+    sampled_mpo_diff_and_scale(left, right, r, n_samples, seed).0
+}
+
+/// Sampled maximum difference divided by the largest sampled magnitude.
+pub fn max_sampled_mpo_relative_diff(
+    left: &MPO<f64>,
+    right: &MPO<f64>,
+    r: usize,
+    n_samples: usize,
+    seed: u64,
+) -> f64 {
+    let (difference, scale) = sampled_mpo_diff_and_scale(left, right, r, n_samples, seed);
+    difference / scale.max(f64::MIN_POSITIVE)
+}
+
+fn sampled_mpo_diff_and_scale(
+    left: &MPO<f64>,
+    right: &MPO<f64>,
+    r: usize,
+    n_samples: usize,
+    seed: u64,
+) -> (f64, f64) {
+    let xs = sample_grid_indices(r, n_samples, seed);
+    let zs = sample_grid_indices(r, n_samples, seed.wrapping_add(1));
+    xs.iter()
+        .zip(&zs)
+        .fold((0.0_f64, 0.0_f64), |(difference, scale), (&ix, &iz)| {
+            let xb = index_to_bits(ix, r);
+            let zb = index_to_bits(iz, r);
+            let indices = (0..r)
+                .flat_map(|site| [xb[site], zb[site]])
+                .collect::<Vec<_>>();
+            let left_value = left.evaluate(&indices).unwrap();
+            let right_value = right.evaluate(&indices).unwrap();
+            (
+                difference.max((left_value - right_value).abs()),
+                scale.max(left_value.abs()).max(right_value.abs()),
+            )
+        })
+}
+
 #[cfg(test)]
 mod tests {
     use super::*;
@@ -420,4 +537,31 @@ mod tests {
             "zipup arms disagree: max diff {diff:e} on scale {scale:e}"
         );
     }
+
+    #[test]
+    fn partitioned_tt_and_chain_treetn_match_each_other() {
+        use crate::patched_mpo::PatchedMpoPair;
+
+        let (r, l, patch_cap) = (6, 6.0, 36);
+        let f = GaussianMixture2D::random(3, l, (0.5, 2.0), 40);
+        let g = GaussianMixture2D::random(3, l, (0.5, 2.0), 41);
+        let (left, _dy) = to_quantics_mpo(&f, r, l, 1e-8, 128).unwrap();
+        let (right, _) = to_quantics_mpo(&g, r, l, 1e-8, 128).unwrap();
+        let patched = PatchedMpoPair::new(&left, &right, 1e-8, patch_cap).unwrap();
+        let input_chi = patched.input_max_bond();
+        assert!(input_chi > patch_cap);
+        assert!(patched.input_patch_counts().0 > 1);
+
+        let tt = patched.contract_fit_tt(1e-10, input_chi).unwrap();
+        let tree = patched.contract_fit_treetn(1e-10, input_chi).unwrap();
+        let cross_diff = max_sampled_mpo_diff(&tt.mpo, &tree.mpo, r, 64, 42);
+
+        assert_eq!(tt.n_patches, tree.n_patches);
+        assert_eq!(tt.max_patch_bond, tree.max_patch_bond);
+        assert!(tt.max_patch_bond <= input_chi);
+        assert!(
+            cross_diff < 1e-8,
+            "partitioned outputs differ by {cross_diff:e}"
+        );
+    }
 }
diff --git a/src/patched.rs b/src/patched.rs
index 720e1b4..987ae23 100644
--- a/src/patched.rs
+++ b/src/patched.rs
@@ -511,7 +511,7 @@ pub fn patched_elementwise_with_stats(
                     pair_t0.elapsed().as_secs_f64()
                 );
             }
-            patches.push(SubDomainTT::new(tt, merged));
+            patches.push(SubDomainTT::new(tt, merged)?);
         }
     }
     stats.n_pairs = patches.len();
@@ -827,7 +827,7 @@ pub fn eval_patched(
             .cloned()
             .zip(fused.iter().copied())
             .collect::<Vec<_>>(),
-    );
+    )?;
     match restrict_to_free(subdomain.data(), sites, &full, &[])? {
         Restricted::Scalar(value) => Ok(value),
         Restricted::Train(_) => anyhow::bail!("a fully fixed restriction kept free sites"),
@@ -1205,7 +1205,7 @@ mod tests {
 
         // One subdomain over the whole domain, so a full projector turns the
         // restriction into an evaluation of the bridged train itself.
-        let whole = PartitionedTT::from_subdomain(SubDomainTT::from_tt(bridged));
+        let whole = PartitionedTT::from_subdomain(SubDomainTT::from_tt(bridged)).unwrap();
         for &(ix, iy) in &[(0u64, 0u64), (5, 27), (31, 1), (16, 16)] {
             let xb = index_to_bits(ix, r);
             let yb = index_to_bits(iy, r);
diff --git a/src/patched_mpo.rs b/src/patched_mpo.rs
new file mode 100644
index 0000000..fc9cb99
--- /dev/null
+++ b/src/patched_mpo.rs
@@ -0,0 +1,344 @@
+//! Patched MPO-MPO contraction on the legacy TT and chain TreeTN wrappers.
+
+use tensor4all_core::{DynIndex, SvdTruncationPolicy};
+use tensor4all_itensorlike::{ContractOptions, TensorTrain};
+use tensor4all_partitionedtreetn as tree;
+use tensor4all_partitionedtt as tt;
+use tensor4all_simplett::mpo::MPO;
+use tensor4all_treetn::contraction::ContractionOptions as TreeContractOptions;
+
+use crate::mpo_contract::{mpo_to_tensortrain, tensortrain_to_mpo};
+
+fn itensor_cutoff_policy(cutoff: f64) -> SvdTruncationPolicy {
+    SvdTruncationPolicy::new(cutoff)
+        .with_relative()
+        .with_squared_values()
+        .with_discarded_tail_sum()
+}
+
+/// Prepared copies of one MPO pair for a fair TT versus chain TreeTN comparison.
+pub struct PatchedMpoPair {
+    x: Vec<DynIndex>,
+    z: Vec<DynIndex>,
+    global_left: TensorTrain,
+    global_right: TensorTrain,
+    tt_left: tt::PartitionedTT,
+    tt_right: tt::PartitionedTT,
+    tree_left: tree::PartitionedTreeTN<usize>,
+    tree_right: tree::PartitionedTreeTN<usize>,
+    tt_options: tt::PatchingOptions,
+    tree_options: tree::PatchingOptions,
+    center: usize,
+    input_max_bond: usize,
+}
+
+/// One contracted patched representation, summed back to an MPO for validation.
+pub struct PatchedMpoOutput {
+    /// Summed output MPO.
+    pub mpo: MPO<f64>,
+    /// Number of output patches before summing.
+    pub n_patches: usize,
+    /// Largest bond dimension among output patches.
+    pub max_patch_bond: usize,
+}
+
+impl PatchedMpoPair {
+    /// Build both patched representations from exactly the same MPO cores and indices.
+    pub fn new(
+        left: &MPO<f64>,
+        right: &MPO<f64>,
+        rtol: f64,
+        patch_max_bond: usize,
+    ) -> anyhow::Result<Self> {
+        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)?;
+        Self::from_tensortrains(left_train, right_train, rtol, patch_max_bond)
+    }
+
+    /// Build patched representations from cached global tensor trains.
+    pub fn from_tensortrains(
+        left_train: TensorTrain,
+        right_train: TensorTrain,
+        rtol: f64,
+        patch_max_bond: usize,
+    ) -> anyhow::Result<Self> {
+        anyhow::ensure!(!left_train.is_empty(), "cached MPO input is empty");
+        anyhow::ensure!(
+            left_train.len() == right_train.len(),
+            "cached MPO lengths differ"
+        );
+        anyhow::ensure!(patch_max_bond > 0, "patch max bond must be positive");
+
+        let left_sites = left_train.site_indices();
+        let right_sites = right_train.site_indices();
+        let mut x = Vec::with_capacity(left_train.len());
+        let mut y = Vec::with_capacity(left_train.len());
+        let mut z = Vec::with_capacity(left_train.len());
+        for site in 0..left_train.len() {
+            let common: Vec<_> = left_sites[site]
+                .iter()
+                .filter(|index| right_sites[site].contains(index))
+                .cloned()
+                .collect();
+            anyhow::ensure!(
+                common.len() == 1,
+                "cached site {site} does not have one shared index"
+            );
+            let left_unique: Vec<_> = left_sites[site]
+                .iter()
+                .filter(|index| *index != &common[0])
+                .cloned()
+                .collect();
+            let right_unique: Vec<_> = right_sites[site]
+                .iter()
+                .filter(|index| *index != &common[0])
+                .cloned()
+                .collect();
+            anyhow::ensure!(
+                left_unique.len() == 1 && right_unique.len() == 1,
+                "cached site {site} does not have MPO-like site indices"
+            );
+            x.push(left_unique[0].clone());
+            y.push(common[0].clone());
+            z.push(right_unique[0].clone());
+        }
+
+        let n = left_train.len();
+        let input_max_bond = left_train.max_bond_dim().max(right_train.max_bond_dim());
+        let left_order = y.iter().chain(&x).cloned().collect::<Vec<_>>();
+        let right_order = y.iter().chain(&z).cloned().collect::<Vec<_>>();
+        let tt_options = tt::PatchingOptions {
+            rtol,
+            max_bond_dim: Some(patch_max_bond),
+            patch_order: left_order.clone(),
+            split_strategy: tt::PatchSplitStrategy::Sequential,
+        };
+        let tt_right_options = tt::PatchingOptions {
+            patch_order: right_order.clone(),
+            ..tt_options.clone()
+        };
+        let tt_left = tt::add_with_patching(
+            vec![tt::SubDomainTT::from_tt(left_train.clone())],
+            &tt_options,
+        )?;
+        let tt_right = tt::add_with_patching(
+            vec![tt::SubDomainTT::from_tt(right_train.clone())],
+            &tt_right_options,
+        )?;
+        let global_left = left_train.clone();
+        let global_right = right_train.clone();
+        let tree_options = tree::PatchingOptions {
+            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(),
+            )?],
+            &center,
+            &tree_options,
+        )?;
+        let tree_right = tree::add_with_patching(
+            vec![tree::SubDomainTreeTN::from_treetn(
+                right_train.into_treetn(),
+            )?],
+            &center,
+            &tree_right_options,
+        )?;
+        anyhow::ensure!(
+            (tree_left.len(), tree_right.len()) == (tt_left.len(), tt_right.len()),
+            "TT and TreeTN patch counts differ"
+        );
+        Ok(Self {
+            x,
+            z,
+            global_left,
+            global_right,
+            tt_left,
+            tt_right,
+            tree_left,
+            tree_right,
+            tt_options,
+            tree_options,
+            center,
+            input_max_bond,
+        })
+    }
+
+    /// Borrow the global tensor trains for input caching.
+    pub fn global_inputs(&self) -> (&TensorTrain, &TensorTrain) {
+        (&self.global_left, &self.global_right)
+    }
+
+    /// Maximum bond dimension before patching.
+    pub fn input_max_bond(&self) -> usize {
+        self.input_max_bond
+    }
+
+    /// Left and right input patch counts.
+    pub fn input_patch_counts(&self) -> (usize, usize) {
+        (self.tt_left.len(), self.tt_right.len())
+    }
+
+    /// Largest left and right input patch bond dimensions.
+    pub fn input_patch_max_bonds(&self) -> (usize, usize) {
+        (
+            self.tt_left
+                .values()
+                .map(|patch| patch.max_bond_dim())
+                .max()
+                .unwrap_or(0),
+            self.tt_right
+                .values()
+                .map(|patch| patch.max_bond_dim())
+                .max()
+                .unwrap_or(0),
+        )
+    }
+
+    /// Run one full-sweep global fit on the prepared unpatched inputs.
+    pub fn contract_fit_global(
+        &self,
+        rtol: f64,
+        output_max_bond: usize,
+    ) -> anyhow::Result<TensorTrain> {
+        let options = ContractOptions::fit()
+            .with_max_bond_dim(output_max_bond)
+            .with_svd_policy(itensor_cutoff_policy(rtol))
+            .with_nsweeps(1);
+        self.global_left
+            .contract(&self.global_right, &options)
+            .map_err(Into::into)
+    }
+
+    /// Convert a prepared global result to an MPO outside the timed region.
+    pub fn finish_global_output(&self, output: &TensorTrain) -> anyhow::Result<MPO<f64>> {
+        tensortrain_to_mpo(output, &self.x, &self.z)
+    }
+
+    /// Run one full-sweep fit through `tensor4all-partitionedtt`.
+    pub fn contract_fit_tt_partitioned(
+        &self,
+        rtol: f64,
+        output_max_bond: usize,
+    ) -> anyhow::Result<tt::PartitionedTT> {
+        let output_options = tt::PatchingOptions {
+            rtol,
+            max_bond_dim: Some(output_max_bond),
+            ..self.tt_options.clone()
+        };
+        let options = ContractOptions::fit()
+            .with_max_bond_dim(output_max_bond)
+            .with_svd_policy(itensor_cutoff_policy(rtol))
+            .with_nsweeps(1);
+        tt::contract_adaptive(&self.tt_left, &self.tt_right, &options, &output_options)
+            .map_err(Into::into)
+    }
+
+    /// Sum and convert a legacy partitioned result outside the timed region.
+    pub fn finish_tt_output(&self, output: tt::PartitionedTT) -> anyhow::Result<PatchedMpoOutput> {
+        let max_patch_bond = output
+            .values()
+            .map(|patch| patch.max_bond_dim())
+            .max()
+            .unwrap_or(0);
+        let mpo = tensortrain_to_mpo(&output.to_tensor_train()?, &self.x, &self.z)?;
+        Ok(PatchedMpoOutput {
+            mpo,
+            n_patches: output.len(),
+            max_patch_bond,
+        })
+    }
+
+    /// Run and convert one legacy partitioned fit contraction.
+    pub fn contract_fit_tt(
+        &self,
+        rtol: f64,
+        output_max_bond: usize,
+    ) -> anyhow::Result<PatchedMpoOutput> {
+        let output = self.contract_fit_tt_partitioned(rtol, output_max_bond)?;
+        self.finish_tt_output(output)
+    }
+
+    /// Run one full-sweep fit through `tensor4all-partitionedtreetn` on a chain.
+    pub fn contract_fit_treetn_partitioned(
+        &self,
+        rtol: f64,
+        output_max_bond: usize,
+    ) -> anyhow::Result<tree::PartitionedTreeTN<usize>> {
+        let output_options = tree::PatchingOptions {
+            rtol,
+            max_bond_dim: Some(output_max_bond),
+            ..self.tree_options.clone()
+        };
+        let options = TreeContractOptions::fit()
+            .with_max_bond_dim(output_max_bond)
+            .with_svd_policy(itensor_cutoff_policy(rtol))
+            .with_nfullsweeps(1);
+        tree::contract_adaptive(
+            &self.tree_left,
+            &self.tree_right,
+            &self.center,
+            &options,
+            &output_options,
+        )
+        .map_err(Into::into)
+    }
+
+    /// Sum and convert a chain TreeTN partitioned result outside the timed region.
+    pub fn finish_treetn_output(
+        &self,
+        output: tree::PartitionedTreeTN<usize>,
+    ) -> anyhow::Result<PatchedMpoOutput> {
+        let max_patch_bond = output
+            .values()
+            .map(|patch| patch.max_bond_dim())
+            .max()
+            .unwrap_or(0);
+        let train = TensorTrain::from_treetn(output.to_treetn()?)?;
+        let mpo = tensortrain_to_mpo(&train, &self.x, &self.z)?;
+        Ok(PatchedMpoOutput {
+            mpo,
+            n_patches: output.len(),
+            max_patch_bond,
+        })
+    }
+
+    /// Run and convert one chain TreeTN partitioned fit contraction.
+    pub fn contract_fit_treetn(
+        &self,
+        rtol: f64,
+        output_max_bond: usize,
+    ) -> anyhow::Result<PatchedMpoOutput> {
+        let output = self.contract_fit_treetn_partitioned(rtol, output_max_bond)?;
+        self.finish_treetn_output(output)
+    }
+}