v0.21.1: TT _dim_order cluster + open issues + perf

CHANGELOG.md

@@ -5,6 +5,63 @@ All notable changes to PyChebyshev will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [0.21.1] - 2026-04-27
+
+### Added
+
+- `ChebyshevTT.sobol_indices()` — first-order and total-order Sobol
+  sensitivity indices computed natively by contracting through the TT
+  coefficient cores. O(d · n · r²) per dim, no dense materialization.
+  Mirrors v0.20 `ChebyshevApproximation.sobol_indices` API; keys are
+  user-frame dim indices regardless of internal `_dim_order`.
+
+### Fixed
+
+- `ChebyshevTT.roots()`/`minimize()`/`maximize()` previously validated
+  `fixed` values against the storage-frame domain instead of the
+  user-frame physical domain. Under `with_auto_order()` or `reorder()`
+  with non-uniform per-dim domains, this could raise misleading errors
+  or accept invalid inputs. Now validates against user-frame domain.
+- `ChebyshevTT.inner_product()` previously returned a meaningless
+  Frobenius product when `self._dim_order != other._dim_order`, with
+  no error. Now raises `ValueError` with a `reorder()` alignment hint,
+  matching v0.20.1 binary algebra behavior.
+- `ChebyshevTT.get_evaluation_points()` previously returned columns in
+  storage order, breaking `eval(get_evaluation_points()[i])` for any
+  TT with non-identity `_dim_order`. Now returns columns in user-frame
+  order.
+- `ChebyshevTT.eval_multi()` previously mutated `self._dim_order` via
+  try/finally, racing under concurrent calls (issue #19). Now uses a
+  private `_eval_storage_frame` helper with no mutation.
+- `ChebyshevTT.integrate()` error messages on out-of-domain bounds
+  previously referenced the storage-frame dim index instead of the
+  user-frame index passed by the caller (issue #20). Now references
+  the user-frame index.
+- `_algebra._check_compatible` previously raised "Domain mismatch"
+  when comparing two interpolants with mixed `tuple` vs `list` domain
+  syntax even when bounds were numerically identical (issue #22). Now
+  uses `np.allclose` for comparison.
+
+### Performance
+
+- `vectorized_eval_batch` now hoists the differentiation matrix matmul
+  outside the per-point loop. Significant speedup for derivative
+  batch evaluations on large point sets.
+- `_calculus._optimize_1d` (used by all four classes' `minimize/maximize`)
+  now uses a single vectorized barycentric evaluation over critical
+  points + endpoints instead of a Python list comprehension.
+
+### Notes
+
+- Closes the v0.20+v0.20.1 `_dim_order` cluster on `ChebyshevTT`.
+  All TT methods that read `self.domain[d]` or `self.n_nodes[d]` now
+  consistently translate user-frame indices to storage-frame
+  internally.
+- No breaking API changes: all "Fixed" items change wrong behavior to
+  correct behavior. The `inner_product` strict-mode raises on mismatched
+  `_dim_order` (previously silently returned wrong numbers), which is
+  a behavior change in the failure path only.
+
 ## [0.21.0] - 2026-04-27
 
 ### Added

CLAUDE.md

@@ -12,7 +12,7 @@ PyChebyshev is a pip-installable Python library for multi-dimensional Chebyshev
 # Setup
 uv sync
 
-# Run tests (~1112 tests, ~115s due to 5D Black-Scholes builds)
+# Run tests (~1133 tests, ~115s due to 5D Black-Scholes builds)
 uv run pytest tests/ -v
 
 # Run a single test
@@ -84,6 +84,12 @@ The installable package. Public classes: `ChebyshevApproximation`, `ChebyshevSpl
 - v0.21 adds `ChebyshevSlider.roots()/minimize()/maximize()` and
   `ChebyshevTT.roots()/minimize()/maximize()`. After v0.21, all four
   classes support the full calculus surface (integrate + roots + min/max).
+- v0.21.1 closes the v0.20+v0.20.1 `_dim_order` cluster on `ChebyshevTT`:
+  `roots/minimize/maximize` validate against user-frame domain;
+  `inner_product` raises on mismatched `_dim_order`; `get_evaluation_points`
+  returns user-frame columns; `eval_multi` no longer mutates `_dim_order`.
+  Perf: vectorized `_optimize_1d` and `vectorized_eval_batch` derivative
+  hoist. Adds `ChebyshevTT.sobol_indices()` parity (native TT contraction).
 
 ### Benchmark Scripts (project root)
 
@@ -105,6 +111,7 @@ Not part of the library. Compare Chebyshev barycentric against alternative metho
 - `compare_calculus_completion.py` — PyChebyshev v0.17 Slider/TT integrate vs MoCaX 4.3.1 (no equivalent — beyond-MoCaX feature)
 - `compare_v018_tt_parity.py` — PyChebyshev v0.18 TT surface (extrude/slice/algebra/from_values/to_dense) vs MoCaX 4.3.1
 - `compare_v019_build_diagnostics.py` — PyChebyshev v0.19 build optimization (parallel eval, progress bars, visualization) — no MoCaX equivalent
+- `compare_v0211_dim_cluster.py` — PyChebyshev v0.21.1 TT `_dim_order` cluster fixes demo (no MoCaX equivalent — internal-correctness fixes)
 
 ### Tests (`tests/`)
 
@@ -129,12 +136,13 @@ Not part of the library. Compare Chebyshev barycentric against alternative metho
   `get_evaluation_points`, `get_num_evaluation_points`), `peek_format_version`,
   `is_dimensionality_allowed`, `defer_build` + `set_original_function_values`,
   `Domain`/`Ns`/`SpecialPoints` typed helpers.
-- `test_calculus_completion.py` — ~101 tests: `ChebyshevSlider.integrate/roots/minimize/maximize`,
+- `test_calculus_completion.py` — ~106 tests: `ChebyshevSlider.integrate/roots/minimize/maximize`,
   `ChebyshevTT.integrate/roots/minimize/maximize` (full and partial,
   user-frame dim/fixed transparent under `_dim_order`), cross-class
   consistency checks, bounds validation. v0.21 additions: 57 tests
   across 9 new test classes covering Slider/TT roots/min/max parity
-  with Approximation/Spline.
+  with Approximation/Spline. v0.21.1 additions: 5 tests covering
+  user-frame domain validation fix under non-identity `_dim_order`.
 - `test_v018_tt_parity.py` — ~52 tests: `ChebyshevTT.nodes()`, `from_values()`,
   `extrude()`, `slice()`, algebra (`+`, `-`, `*` scalar, in-place, `__neg__`),
   `to_dense()`; cross-feature and round-trip checks.

compare_v0211_dim_cluster.py

@@ -0,0 +1,114 @@
+"""v0.21.1 demo: TT _dim_order cluster fixes.
+
+Demonstrates each correctness fix with a non-uniform-domain TT under
+explicit reorder([2, 0, 1]) — the case that v0.20.1 / v0.21.0 tests
+masked with uniform domains.
+"""
+
+from __future__ import annotations
+
+import time
+
+import numpy as np
+
+from pychebyshev import ChebyshevApproximation, ChebyshevTT
+
+
+def _check(label: str, ok: bool, detail: str = "") -> None:
+    status = "OK " if ok else "FAIL"
+    print(f"  [{status}] {label}{(' — ' + detail) if detail else ''}")
+
+
+def demo_inner_product_strict() -> None:
+    print("\n=== inner_product strict _dim_order check (Item B) ===")
+    def f(x, _): return x[0] ** 2 + x[1]
+    tt = ChebyshevTT(f, num_dimensions=2, domain=[(-1, 1), (-1, 1)], n_nodes=[5, 5])
+    tt.build(verbose=False)
+    tt_p = tt.reorder([1, 0])
+    try:
+        tt.inner_product(tt_p)
+        _check("ValueError raised on dim_order mismatch", False)
+    except ValueError as e:
+        _check("ValueError raised on dim_order mismatch", True, str(e)[:80])
+
+
+def demo_get_evaluation_points_round_trip() -> None:
+    print("\n=== get_evaluation_points user-frame round-trip (Item C) ===")
+    def f(x, _): return 0.3 * x[0] + 0.7 * x[1] - 0.2 * x[2]
+    tt = ChebyshevTT(
+        f, num_dimensions=3,
+        domain=[(-1, 1), (-2, 2), (-3, 3)], n_nodes=[5, 5, 5],
+    )
+    tt.build(verbose=False)
+    tt_p = tt.reorder([2, 0, 1])
+    points = tt_p.get_evaluation_points()
+    max_err = 0.0
+    for i in range(0, len(points), 25):
+        pt = points[i]
+        expected = f(pt, None)
+        got = float(tt_p.eval(pt.tolist()))
+        max_err = max(max_err, abs(got - expected))
+    _check("round-trip eval == f for non-identity _dim_order",
+           max_err < 1e-9, f"max_err={max_err:.2e}")
+
+
+def demo_roots_user_frame_validation() -> None:
+    print("\n=== roots/min/max validate against user-frame domain (Item A) ===")
+    # User-frame dim 1 has domain [-2, 2]; storage-frame after reorder has different range
+    def f(x, _): return (x[0] - 0.4) * (1.0 + 0.0 * x[1] + 0.0 * x[2])
+    tt = ChebyshevTT(
+        f, num_dimensions=3,
+        domain=[(-1, 1), (-2, 2), (-3, 3)], n_nodes=[8, 8, 8],
+    )
+    tt.build(verbose=False)
+    tt_p = tt.reorder([2, 0, 1])
+    # fixed=1.5 is valid in user-frame dim 1, NOT in storage-frame after reorder
+    try:
+        roots = tt_p.roots(dim=0, fixed={1: 1.5, 2: 0.0})
+        _check("roots accepts user-frame-valid fixed value",
+               abs(float(roots[0]) - 0.4) < 1e-7,
+               f"root={float(roots[0]):.4f}")
+    except Exception as e:
+        _check("roots accepts user-frame-valid fixed value", False,
+               f"raised {type(e).__name__}: {e}")
+
+
+def demo_integrate_error_user_frame() -> None:
+    print("\n=== integrate error message uses user-frame dim (Item E / #20) ===")
+    def f(x, _): return x[0] + x[1] + x[2]
+    tt = ChebyshevTT(
+        f, num_dimensions=3,
+        domain=[(-1, 1), (-2, 2), (-3, 3)], n_nodes=[5, 5, 5],
+    )
+    tt.build(verbose=False)
+    tt_p = tt.reorder([2, 0, 1])
+    try:
+        tt_p.integrate(dims=[1], bounds=[(5.0, 6.0)])
+        _check("ValueError raised on out-of-domain bounds", False)
+    except ValueError as e:
+        msg = str(e)
+        _check("error references user-frame dim 1",
+               "dim 1" in msg, msg[:100])
+
+
+def demo_eval_multi_no_mutation() -> None:
+    print("\n=== eval_multi no longer mutates _dim_order (Item D / #19) ===")
+    import inspect
+    source = inspect.getsource(ChebyshevTT.eval_multi)
+    _check("eval_multi source contains no 'self._dim_order = ' assignment",
+           "self._dim_order = " not in source and "self._dim_order=" not in source)
+
+
+def main() -> None:
+    print("PyChebyshev v0.21.1 cluster fix demo")
+    t0 = time.time()
+    demo_inner_product_strict()
+    demo_get_evaluation_points_round_trip()
+    demo_roots_user_frame_validation()
+    demo_integrate_error_user_frame()
+    demo_eval_multi_no_mutation()
+    print(f"\nTotal: {time.time() - t0:.2f}s")
+
+
+if __name__ == "__main__":
+    main()

docs/user-guide/calculus.md

@@ -570,6 +570,8 @@ roots_b = tt_optimized.roots(dim=0, fixed={1: 0.0, 2: 0.0})
 np.testing.assert_array_almost_equal(roots_a, roots_b)
 ```
 
+> **Non-uniform domains:** v0.21.1 closed a latent bug where TT calculus methods validated `fixed` values against the storage-frame domain. With non-uniform per-dim domains and a non-identity `_dim_order` (after `with_auto_order` / `reorder`), this could either reject valid user-frame inputs or silently accept invalid ones. Since v0.21.1, validation always uses the user-frame physical domain.
+
 ### Slider example
 
 ```python

pyproject.toml

@@ -4,7 +4,7 @@ build-backend = "hatchling.build"
 
 [project]
 name = "pychebyshev"
-version = "0.21.0"
+version = "0.21.1"
 description = "Fast multi-dimensional Chebyshev tensor interpolation with analytical derivatives"
 readme = "README.md"
 license = {text = "MIT"}

src/pychebyshev/_algebra.py

@@ -38,12 +38,15 @@ def _check_compatible(a, b) -> None:
             f"Dimension mismatch: {a.num_dimensions} vs {b.num_dimensions}"
         )
 
-    if a.n_nodes != b.n_nodes:
+    if not np.array_equal(np.asarray(a.n_nodes, dtype=int), np.asarray(b.n_nodes, dtype=int)):
         raise ValueError(
             f"Node count mismatch: {a.n_nodes} vs {b.n_nodes}"
         )
 
-    if a.domain != b.domain:
+    if not np.allclose(
+        np.asarray(a.domain, dtype=float),
+        np.asarray(b.domain, dtype=float),
+    ):
         raise ValueError(
             f"Domain mismatch: {a.domain} vs {b.domain}"
         )

src/pychebyshev/_calculus.py

@@ -133,7 +133,7 @@ def _compute_sub_interval_weights(n: int, t_lo: float,
     return weights_desc[::-1].copy()
 
 
-def _normalize_bounds(dims, bounds, domain):
+def _normalize_bounds(dims, bounds, domain, dim_labels=None):
     """Normalize and validate the *bounds* parameter for ``integrate()``.
 
     Parameters
@@ -144,6 +144,12 @@ def _normalize_bounds(dims, bounds, domain):
         User-provided bounds specification.
     domain : list
         Per-dimension ``[lo, hi]`` bounds.
+    dim_labels : list of int or None
+        Optional override for dim labels used in error messages. When
+        ``None`` (default), error messages use the corresponding entry
+        from ``dims``. When provided, ``dim_labels[i]`` is used instead
+        of ``dims[i]`` so callers can present user-frame indices when
+        ``dims`` is in storage frame.
 
     Returns
     -------
@@ -173,14 +179,15 @@ def _normalize_bounds(dims, bounds, domain):
             result.append(None)
             continue
         lo, hi = bd
+        label = dim_labels[i] if dim_labels is not None else dims[i]
         if lo > hi:
-            raise ValueError(f"bounds lo={lo} > hi={hi} for dim {dims[i]}")
+            raise ValueError(f"bounds lo={lo} > hi={hi} for dim {label}")
         d = dims[i]
         dom_lo, dom_hi = domain[d]
         if lo < dom_lo - 1e-14 or hi > dom_hi + 1e-14:
             raise ValueError(
                 f"bounds ({lo}, {hi}) outside domain [{dom_lo}, {dom_hi}] "
-                f"for dim {d}"
+                f"for dim {label}"
             )
         lo = max(lo, dom_lo)
         hi = min(hi, dom_hi)
@@ -259,8 +266,6 @@ def _optimize_1d(values: np.ndarray, nodes: np.ndarray,
     -------
     (value, location) : (float, float)
     """
-    from pychebyshev.barycentric import barycentric_interpolate
-
     # Derivative values at nodes
     deriv_values = diff_matrix @ values
 
@@ -271,11 +276,22 @@ def _optimize_1d(values: np.ndarray, nodes: np.ndarray,
     a, b = domain
     candidates = np.concatenate([[a], critical, [b]])
 
-    # Evaluate original function at all candidates
-    vals = np.array([
-        barycentric_interpolate(float(x), nodes, values, bary_weights)
-        for x in candidates
-    ])
+    # Vectorized barycentric evaluation at all candidates simultaneously.
+    candidates_arr = np.asarray(candidates, dtype=float).reshape(-1)
+    diff = candidates_arr[:, None] - nodes[None, :]   # shape (M, n)
+    abs_diff = np.abs(diff)
+    exact_mask = abs_diff < 1e-14                      # (M, n)
+    has_exact = exact_mask.any(axis=1)                 # (M,)
+    # Replace zero diffs with 1.0 to avoid division by zero; overwritten below.
+    safe_diff = np.where(abs_diff < 1e-14, 1.0, diff)
+    w_over_diff = bary_weights[None, :] / safe_diff    # (M, n)
+    numer = (w_over_diff * values[None, :]).sum(axis=1)
+    denom = w_over_diff.sum(axis=1)
+    vals = numer / denom                                # (M,)
+    # For candidates that hit a node exactly, take the node value directly.
+    if has_exact.any():
+        exact_node_idx = exact_mask.argmax(axis=1)
+        vals = np.where(has_exact, values[exact_node_idx], vals)
 
     idx = np.argmin(vals) if mode == "min" else np.argmax(vals)
     return float(vals[idx]), float(candidates[idx])