From 327aec327e0108df58906ca44ecb320fd257b94b Mon Sep 17 00:00:00 2001 From: Max Zhang Date: Mon, 27 Apr 2026 12:06:18 -0400 Subject: [PATCH 01/16] v0.21.1: ChebyshevTT.inner_product strict _dim_order check Raise ValueError with reorder() hint when self._dim_order != other._dim_order instead of silently computing a meaningless Frobenius product. Mirrors the v0.20.1 binary algebra behavior on _check_compatible_tt mismatches. Closes ultrareview Critical finding. --- src/pychebyshev/tensor_train.py | 8 ++++++++ tests/test_tensor_train.py | 35 +++++++++++++++++++++++++++++++++ 2 files changed, 43 insertions(+) diff --git a/src/pychebyshev/tensor_train.py b/src/pychebyshev/tensor_train.py index c6f0a59..c3b9ff2 100644 --- a/src/pychebyshev/tensor_train.py +++ b/src/pychebyshev/tensor_train.py @@ -1482,6 +1482,14 @@ def inner_product(self, other: "ChebyshevTT") -> float: "inner_product requires matching n_nodes; " f"got {self.n_nodes} vs {other.n_nodes}" ) + if list(self._dim_order) != list(other._dim_order): + raise ValueError( + f"inner_product requires matching _dim_order: " + f"{self._dim_order} vs {other._dim_order}. " + f"Call other = other.reorder(self._dim_order) " + f"(or self = self.reorder(other._dim_order)) to align " + f"before computing inner_product." + ) M = np.array([[1.0]]) # (r_self_0, r_other_0) = (1, 1) for k in range(self.num_dimensions): diff --git a/tests/test_tensor_train.py b/tests/test_tensor_train.py index 9dd6ce4..926e31d 100644 --- a/tests/test_tensor_train.py +++ b/tests/test_tensor_train.py @@ -587,6 +587,41 @@ def test_inner_product_raises_on_unbuilt_other(self): with pytest.raises(RuntimeError, match="Call build"): tt_a.inner_product(tt_b) + def test_inner_product_dim_order_mismatch_raises(self): + """Two TTs with different _dim_order must raise ValueError.""" + def f(x, _): return x[0] ** 2 + x[1] + x[2] ** 2 + tt_canonical = ChebyshevTT( + f, num_dimensions=3, domain=[(-1, 1)] * 3, n_nodes=[5, 5, 5] + ) + tt_canonical.build(verbose=False) + tt_reordered = tt_canonical.reorder([2, 0, 1]) + with pytest.raises(ValueError, match="dim_order"): + tt_canonical.inner_product(tt_reordered) + # Also the reverse direction + with pytest.raises(ValueError, match="dim_order"): + tt_reordered.inner_product(tt_canonical) + + def test_inner_product_after_explicit_reorder_alignment(self): + """After explicit reorder() to align dim_orders, inner_product + succeeds and matches the canonical-vs-canonical result.""" + def f(x, _): return x[0] ** 2 + x[1] + x[2] ** 2 + tt_a = ChebyshevTT( + f, num_dimensions=3, domain=[(-1, 1)] * 3, n_nodes=[5, 5, 5] + ) + tt_b = ChebyshevTT( + f, num_dimensions=3, domain=[(-1, 1)] * 3, n_nodes=[5, 5, 5] + ) + tt_a.build(verbose=False) + tt_b.build(verbose=False) + baseline = tt_a.inner_product(tt_b) + # Permute b, then realign before calling inner_product + tt_b_permuted = tt_b.reorder([2, 0, 1]) + tt_b_realigned = tt_b_permuted.reorder([0, 1, 2]) + aligned = tt_a.inner_product(tt_b_realigned) + assert abs(baseline - aligned) < 1e-9, ( + f"baseline {baseline} != aligned {aligned}" + ) + class TestALSInternals: """White-box tests for the internal ALS sweep primitive.""" From 7b72a06a7d45172eb10840de6e19f46fab4f0703 Mon Sep 17 00:00:00 2001 From: Max Zhang Date: Mon, 27 Apr 2026 12:08:39 -0400 Subject: [PATCH 02/16] =?UTF-8?q?v0.21.1:=20=5Fcheck=5Fcompatible=20?= =?UTF-8?q?=E2=80=94=20numerical=20domain=20comparison?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Replace exact == / != on .domain and .n_nodes with np.allclose / np.array_equal on coerced arrays. Tolerates tuple-of-tuples vs list-of-lists syntax that silently arises in scalar-mul (which normalizes via [list(b) for b in domain]). Closes issue #22. --- src/pychebyshev/_algebra.py | 7 ++++-- tests/test_algebra.py | 47 +++++++++++++++++++++++++++++++++++++ 2 files changed, 52 insertions(+), 2 deletions(-) diff --git a/src/pychebyshev/_algebra.py b/src/pychebyshev/_algebra.py index 4dd5141..7cd9c15 100644 --- a/src/pychebyshev/_algebra.py +++ b/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}" ) diff --git a/tests/test_algebra.py b/tests/test_algebra.py index ca23a35..49fbdd9 100644 --- a/tests/test_algebra.py +++ b/tests/test_algebra.py @@ -697,3 +697,50 @@ def test_portfolio_greeks(self, instruments): delta_strad = straddle.vectorized_eval(p, [1, 0]) weighted = 0.4 * delta_call + 0.3 * delta_put + 0.3 * delta_strad assert abs(delta_port - weighted) < 1e-10, f"Delta mismatch at {p}" + + +class TestMixedDomainSyntax: + """Issue #22: _check_compatible should treat tuple-of-tuples and + list-of-lists domain syntax as numerically equal.""" + + def test_chebyshev_algebra_tuple_vs_list_domain(self): + """ChebyshevApproximation: tuple domain + list domain combine OK.""" + def f(x, _): return x[0] + def g(x, _): return -x[0] + cheb_tuple = ChebyshevApproximation(f, 1, [(-1, 1)], [5]) + cheb_list = ChebyshevApproximation(g, 1, [[-1, 1]], [5]) + cheb_tuple.build(verbose=False) + cheb_list.build(verbose=False) + h = cheb_tuple + cheb_list + # f + g = x - x = 0 + assert abs(h.eval([0.5], [0])) < 1e-12 + + def test_slider_algebra_tuple_vs_list_domain(self): + """ChebyshevSlider: same parity through scalar mul + add.""" + def f(x, _): return x[0] + def g(x, _): return -x[0] + s_f = ChebyshevSlider( + f, num_dimensions=1, domain=[(-1, 1)], n_nodes=[5], + partition=[[0]], pivot_point=[0.0], + ) + s_g = ChebyshevSlider( + g, num_dimensions=1, domain=[(-1, 1)], n_nodes=[5], + partition=[[0]], pivot_point=[0.0], + ) + s_f.build(verbose=False) + s_g.build(verbose=False) + # 2.0 * s_g internally normalizes domain to list-of-lists; + # adding to s_f (still tuple-of-tuples) must succeed. + h = s_f + 2.0 * s_g + # f + 2*g = x - 2x = -x; at x=0.3 should be -0.3 + assert abs(h.eval([0.3], [0]) - (-0.3)) < 1e-12 + + def test_n_nodes_int_vs_list_difference_still_fails(self): + """Sanity: legitimate n_nodes mismatch still raises.""" + def f(x, _): return x[0] + cheb_5 = ChebyshevApproximation(f, 1, [[-1, 1]], [5]) + cheb_7 = ChebyshevApproximation(f, 1, [[-1, 1]], [7]) + cheb_5.build(verbose=False) + cheb_7.build(verbose=False) + with pytest.raises(ValueError, match="[Nn]ode"): + _ = cheb_5 + cheb_7 From ed9bd0adc75cdeb6f23da28cc53782962661b1ac Mon Sep 17 00:00:00 2001 From: Max Zhang Date: Mon, 27 Apr 2026 12:11:13 -0400 Subject: [PATCH 03/16] v0.21.1: ChebyshevTT.get_evaluation_points returns user-frame columns Permute columns by inverse _dim_order before returning so that eval(get_evaluation_points()[i]) round-trips for any TT regardless of storage permutation. Matches Approximation/Spline/Slider behavior. Closes ultrareview Important finding. Co-Authored-By: Claude Sonnet 4.6 --- src/pychebyshev/tensor_train.py | 10 ++++- tests/test_v0201_dim_threading.py | 69 +++++++++++++++++++++++++++++++ 2 files changed, 78 insertions(+), 1 deletion(-) diff --git a/src/pychebyshev/tensor_train.py b/src/pychebyshev/tensor_train.py index c3b9ff2..a773083 100644 --- a/src/pychebyshev/tensor_train.py +++ b/src/pychebyshev/tensor_train.py @@ -2722,6 +2722,8 @@ def get_evaluation_points(self) -> np.ndarray: ------- np.ndarray Shape ``(N, num_dimensions)`` where ``N = prod(n_nodes)``. + Columns are in user-frame order: column ``d`` corresponds to + user-frame dim ``d`` regardless of internal ``_dim_order``. """ per_dim = [] for d in range(self.num_dimensions): @@ -2729,7 +2731,13 @@ def get_evaluation_points(self) -> np.ndarray: a, b = self.domain[d] per_dim.append(np.sort(0.5 * (a + b) + 0.5 * (b - a) * nodes_std)) grids = np.meshgrid(*per_dim, indexing="ij") - return np.stack([g.ravel() for g in grids], axis=-1).astype(np.float64) + # grids columns are in storage frame. Permute back to user frame: + # user-frame dim u corresponds to storage position + # self._dim_order.index(u), so grids[self._dim_order.index(u)] is + # the array of values for user-frame dim u. + user_frame_grids = [grids[self._dim_order.index(u)] + for u in range(self.num_dimensions)] + return np.stack([g.ravel() for g in user_frame_grids], axis=-1).astype(np.float64) def clone(self) -> "ChebyshevTT": """Return an independent deep copy of this interpolant. diff --git a/tests/test_v0201_dim_threading.py b/tests/test_v0201_dim_threading.py index 48c3619..c844c8a 100644 --- a/tests/test_v0201_dim_threading.py +++ b/tests/test_v0201_dim_threading.py @@ -572,3 +572,72 @@ def test_save_load_pcb_not_supported_for_tt(self): restored = pickle.loads(blob) assert restored.dim_order == [2, 0, 1] assert abs(restored.eval([0.1, -0.2, 0.3]) - tt.eval([0.1, -0.2, 0.3])) < 1e-9 + + +class TestGetEvaluationPointsFrame: + """v0.21.1: get_evaluation_points must return columns in user-frame + order so that eval(get_evaluation_points()[i]) round-trips.""" + + def test_round_trip_canonical(self): + """Canonical (identity _dim_order) — round-trip works.""" + 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) + points = tt.get_evaluation_points() + # Sanity: shape and column count + assert points.shape == (5 * 5 * 5, 3) + # Round-trip a sample of points + for i in [0, 7, 31, 50, 124]: + pt = points[i] + expected = f(pt, None) + got = float(tt.eval(pt.tolist())) + assert abs(got - expected) < 1e-9, ( + f"point {i}={pt}: got {got}, expected {expected}" + ) + + def test_round_trip_after_reorder(self): + """Non-identity _dim_order via reorder([2, 0, 1]) — round-trip + must still work in user-frame coordinates.""" + 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_reordered = tt.reorder([2, 0, 1]) + # _dim_order is now non-identity + assert tt_reordered._dim_order != [0, 1, 2] + points = tt_reordered.get_evaluation_points() + assert points.shape == (5 * 5 * 5, 3) + # Round-trip in user-frame + for i in [0, 7, 31, 50, 124]: + pt = points[i] + expected = f(pt, None) + got = float(tt_reordered.eval(pt.tolist())) + assert abs(got - expected) < 1e-9, ( + f"point {i}={pt}: got {got}, expected {expected}" + ) + + def test_columns_match_user_frame_domain_after_reorder(self): + """After reorder, column d's range must match the user-frame + domain[d], not the storage-frame domain.""" + 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_reordered = tt.reorder([2, 0, 1]) + points = tt_reordered.get_evaluation_points() + # Column 0 must be in user-frame dim 0 range: [-1, 1] + assert points[:, 0].min() >= -1.0 - 1e-12 + assert points[:, 0].max() <= 1.0 + 1e-12 + # Column 1 must be in user-frame dim 1 range: [-2, 2] + assert points[:, 1].min() >= -2.0 - 1e-12 + assert points[:, 1].max() <= 2.0 + 1e-12 + # Column 2 must be in user-frame dim 2 range: [-3, 3] + assert points[:, 2].min() >= -3.0 - 1e-12 + assert points[:, 2].max() <= 3.0 + 1e-12 From 9eb2a9cd343e5690a6668a59434517a5b8a34e03 Mon Sep 17 00:00:00 2001 From: Max Zhang Date: Mon, 27 Apr 2026 12:16:05 -0400 Subject: [PATCH 04/16] v0.21.1: ChebyshevTT.roots/min/max validate against user-frame domain Add _user_frame_domain() private helper that returns self.domain permuted into user-frame order. The three v0.21 methods now pass this user-frame view to _calculus._validate_calculus_args instead of the raw storage-frame self.domain. Pre-fix: under non-identity _dim_order with non-uniform per-dim domains, validation referenced the wrong domain. Tests pass before the fix only because uniform domains masked the difference. Closes ultrareview Critical finding. Co-Authored-By: Claude Sonnet 4.6 --- src/pychebyshev/tensor_train.py | 18 +++++++-- tests/test_calculus_completion.py | 61 +++++++++++++++++++++++++++++++ 2 files changed, 76 insertions(+), 3 deletions(-) diff --git a/src/pychebyshev/tensor_train.py b/src/pychebyshev/tensor_train.py index a773083..cb9cfd1 100644 --- a/src/pychebyshev/tensor_train.py +++ b/src/pychebyshev/tensor_train.py @@ -1728,6 +1728,18 @@ def _to_1d_chebyshev(self, sliced_1d): n_nodes=[int(sliced_1d.n_nodes[0])], ) + def _user_frame_domain(self) -> list: + """Return ``self.domain`` reordered into user-frame indexing. + + For canonical ``_dim_order``, this returns ``self.domain`` unchanged + (in semantic terms — a fresh list either way). + For non-identity ``_dim_order``, ``self.domain[d]`` is the + storage-frame domain at storage position ``d``; user-frame dim ``u`` + lives at storage position ``self._dim_order.index(u)``. + """ + return [self.domain[self._dim_order.index(u)] + for u in range(self.num_dimensions)] + def roots(self, dim=None, fixed=None): """Find all roots of the TT-approximated function along *dim*. @@ -1764,7 +1776,7 @@ def roots(self, dim=None, fixed=None): from pychebyshev._calculus import _validate_calculus_args dim, slice_params = _validate_calculus_args( - self.num_dimensions, dim, fixed, self.domain + self.num_dimensions, dim, fixed, self._user_frame_domain() ) sliced = self.slice(slice_params) if slice_params else self @@ -1805,7 +1817,7 @@ def minimize(self, dim=None, fixed=None): from pychebyshev._calculus import _validate_calculus_args dim, slice_params = _validate_calculus_args( - self.num_dimensions, dim, fixed, self.domain + self.num_dimensions, dim, fixed, self._user_frame_domain() ) sliced = self.slice(slice_params) if slice_params else self @@ -1846,7 +1858,7 @@ def maximize(self, dim=None, fixed=None): from pychebyshev._calculus import _validate_calculus_args dim, slice_params = _validate_calculus_args( - self.num_dimensions, dim, fixed, self.domain + self.num_dimensions, dim, fixed, self._user_frame_domain() ) sliced = self.slice(slice_params) if slice_params else self diff --git a/tests/test_calculus_completion.py b/tests/test_calculus_completion.py index c2af0e1..e23f06e 100644 --- a/tests/test_calculus_completion.py +++ b/tests/test_calculus_completion.py @@ -1282,6 +1282,67 @@ def f(x, _): return x[0] ** 2 - 0.25 roots_after = tt_loaded.roots() np.testing.assert_array_almost_equal(roots_before, roots_after, decimal=10) + def test_roots_non_uniform_domain_after_reorder(self): + """Non-uniform domain + reorder: fixed= must validate against + user-frame domain, and roots must be in user-frame coordinates.""" + # f(x, y, z) = (x - 0.4) (1 + 0*y + 0*z) + # User-frame dim 0: domain [-1, 1], root at x=0.4 + # User-frame dim 1: domain [-2, 2] + # User-frame dim 2: domain [-3, 3] + 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_reordered = tt.reorder([2, 0, 1]) + # fixed[1] = 1.5 is valid in user-frame dim 1 ([-2, 2]). + # In storage-frame after reorder([2,0,1]), storage dim 1 is original dim 0 + # with domain [-1, 1]; 1.5 would be out of that storage-frame range. + # Pre-fix: validation against storage-frame raises misleading error. + roots = tt_reordered.roots(dim=0, fixed={1: 1.5, 2: 0.0}) + assert len(roots) == 1 + assert abs(roots[0] - 0.4) < 1e-7 + + def test_roots_non_uniform_domain_user_frame_out_of_range_raises(self): + """A fixed value outside the user-frame domain must raise ValueError.""" + 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_reordered = tt.reorder([2, 0, 1]) + # User-frame dim 1 has range [-2, 2]; passing 5.0 must fail + with pytest.raises(ValueError, match="domain"): + tt_reordered.roots(dim=0, fixed={1: 5.0, 2: 0.0}) + + def test_minimize_non_uniform_domain_after_reorder(self): + """Same non-uniform-domain pattern for minimize.""" + def f(x, _): return (x[0] - 0.3) ** 2 + 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_reordered = tt.reorder([2, 0, 1]) + val, loc = tt_reordered.minimize(dim=0, fixed={1: 1.5, 2: 0.0}) + assert abs(val - 0.0) < 1e-7 + assert abs(loc - 0.3) < 1e-7 + + def test_maximize_non_uniform_domain_after_reorder(self): + """Same non-uniform-domain pattern for maximize.""" + def f(x, _): return -(x[0] - 0.3) ** 2 + 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_reordered = tt.reorder([2, 0, 1]) + val, loc = tt_reordered.maximize(dim=0, fixed={1: 1.5, 2: 0.0}) + assert abs(val - 0.0) < 1e-7 + assert abs(loc - 0.3) < 1e-7 + # ============================================================================ # TestCrossClassCalculusConsistency From d3c8d7ee17479cd08b072a02739e029088503641 Mon Sep 17 00:00:00 2001 From: Max Zhang Date: Mon, 27 Apr 2026 12:18:40 -0400 Subject: [PATCH 05/16] v0.21.1: integrate error messages reference user-frame dim (issue #20) Add optional dim_labels= parameter to _normalize_bounds. ChebyshevTT.integrate passes user-frame dims_sorted as dim_labels so out-of-domain bounds errors reference the dim the user passed rather than the storage-frame index. Helper stays user-frame-naive otherwise. Closes #20. Co-Authored-By: Claude Sonnet 4.6 --- src/pychebyshev/_calculus.py | 13 ++++++++++--- src/pychebyshev/tensor_train.py | 5 ++++- tests/test_calculus_completion.py | 17 +++++++++++++++++ 3 files changed, 31 insertions(+), 4 deletions(-) diff --git a/src/pychebyshev/_calculus.py b/src/pychebyshev/_calculus.py index 1ae95b8..a872f6f 100644 --- a/src/pychebyshev/_calculus.py +++ b/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) diff --git a/src/pychebyshev/tensor_train.py b/src/pychebyshev/tensor_train.py index cb9cfd1..28b5837 100644 --- a/src/pychebyshev/tensor_train.py +++ b/src/pychebyshev/tensor_train.py @@ -1570,9 +1570,12 @@ def integrate( # Normalize bounds: validates against domain in *storage* frame using # storage positions (since `self.domain[storage_pos]` gives the # physical domain of the corresponding original dim after reorder). + # Pass dims_sorted as dim_labels so error messages reference the + # user-frame dim the caller passed, not the storage-frame index. bounds_storage_dims = [storage_for[d] for d in dims_sorted] normalized_bounds = _normalize_bounds( - bounds_storage_dims, bounds, self.domain + bounds_storage_dims, bounds, self.domain, + dim_labels=dims_sorted, ) # Compute scaled quadrature weights per *storage* position. diff --git a/tests/test_calculus_completion.py b/tests/test_calculus_completion.py index e23f06e..858894f 100644 --- a/tests/test_calculus_completion.py +++ b/tests/test_calculus_completion.py @@ -311,6 +311,23 @@ def f(x, _): loaded = ChebyshevTT.load(str(path)) assert loaded.eval([0.5]) == pytest.approx(result.eval([0.5]), abs=1e-12) + def test_out_of_domain_error_uses_user_frame_dim(self): + """Issue #20: error message must reference the user-frame dim + the user passed, not the storage-frame index.""" + 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_reordered = tt.reorder([2, 0, 1]) + # User passes dims=[1] (user-frame dim 1, range [-2, 2]). + # Pre-fix: error would reference dim 2 (storage position). + with pytest.raises(ValueError) as exc_info: + tt_reordered.integrate(dims=[1], bounds=[(5.0, 6.0)]) + msg = str(exc_info.value) + assert "dim 1" in msg, f"Error message references wrong dim: {msg}" + # ============================================================================ # T6: Slider full integration (returns scalar) From 1ddae32e85ade1dccea84628dc2ad81d7f75368c Mon Sep 17 00:00:00 2001 From: Max Zhang Date: Mon, 27 Apr 2026 12:22:32 -0400 Subject: [PATCH 06/16] =?UTF-8?q?v0.21.1:=20eval=5Fmulti=20structural=20fi?= =?UTF-8?q?x=20=E2=80=94=20eliminate=20=5Fdim=5Forder=20mutation=20(issue?= =?UTF-8?q?=20#19)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Extract _eval_storage_frame(point_storage, deriv_storage) private helper that takes already-permuted coords. eval() permutes user→storage once at entry and calls _eval_storage_frame; eval_multi() permutes once and calls the helper N times for each derivative_order. Drops the try/finally self._dim_order mutation entirely. Eliminates the race rather than papering over it with a lock. The FD machinery (_fd_single_dim, _fd_cross_deriv, _fd_nested) was updated to call _eval_storage_frame instead of self.eval, since the points it operates on are already in storage frame; routing through self.eval would re-permute and break correctness once the saved-state trick is removed. Closes #19. --- src/pychebyshev/tensor_train.py | 158 +++++++++++++++++++----------- tests/test_v0201_dim_threading.py | 42 ++++++++ 2 files changed, 142 insertions(+), 58 deletions(-) diff --git a/src/pychebyshev/tensor_train.py b/src/pychebyshev/tensor_train.py index 28b5837..f74332b 100644 --- a/src/pychebyshev/tensor_train.py +++ b/src/pychebyshev/tensor_train.py @@ -2157,23 +2157,59 @@ def eval(self, point: List[float]) -> float: # dim_order was set by with_auto_order(). Identity order is a no-op. canonical = list(range(self.num_dimensions)) if self._dim_order != canonical: - point = [point[self._dim_order[k]] for k in range(self.num_dimensions)] + point_storage = [ + point[self._dim_order[k]] for k in range(self.num_dimensions) + ] + else: + point_storage = list(point) - result = np.ones((1, 1)) - for d in range(self.num_dimensions): - # Map point[d] from [a, b] to [-1, 1] - a, b = self.domain[d] - scaled = 2.0 * (point[d] - a) / (b - a) - 1.0 - # Evaluate all Chebyshev polynomials T_0..T_{n-1} at scaled point - q = np.polynomial.chebyshev.chebval( - scaled, np.eye(self.n_nodes[d]) - ) # shape (n_d,) - # Contract polynomial values with coefficient core: - # v[i,k] = sum_j q[j] * core[i,j,k] → shape (r_{d-1}, r_d) - v = np.einsum("j,ijk->ik", q, self._coeff_cores[d]) - # Chain multiply: result = result @ v - result = result @ v - return float(result[0, 0]) + zero_deriv = [0] * self.num_dimensions + return self._eval_storage_frame(point_storage, zero_deriv) + + def _eval_storage_frame( + self, + point_storage: List[float], + derivative_order_storage: List[int], + ) -> float: + """Evaluate at a single point assuming storage-frame inputs. + + This is the structural workhorse used by both :meth:`eval` (after + a single user→storage permutation) and :meth:`eval_multi` (which + permutes once then calls this helper N times). Eliminates the + need for any temporary mutation of ``self._dim_order``. + + Parameters + ---------- + point_storage : list of float + Query coordinates in **storage frame** — i.e., already + permuted by ``self._dim_order`` if non-canonical. + derivative_order_storage : list of int + Derivative orders in **storage frame**. ``[0, 0, ..., 0]`` + triggers the core-contraction value path; otherwise the + finite-difference machinery is used. + + Returns + ------- + float + Interpolated value (or FD derivative) at the query point. + """ + if all(d == 0 for d in derivative_order_storage): + result = np.ones((1, 1)) + for d in range(self.num_dimensions): + # Map point_storage[d] from [a, b] to [-1, 1] + a, b = self.domain[d] + scaled = 2.0 * (point_storage[d] - a) / (b - a) - 1.0 + # Evaluate all Chebyshev polynomials T_0..T_{n-1} at scaled point + q = np.polynomial.chebyshev.chebval( + scaled, np.eye(self.n_nodes[d]) + ) # shape (n_d,) + # Contract polynomial values with coefficient core: + # v[i,k] = sum_j q[j] * core[i,j,k] → shape (r_{d-1}, r_d) + v = np.einsum("j,ijk->ik", q, self._coeff_cores[d]) + # Chain multiply: result = result @ v + result = result @ v + return float(result[0, 0]) + return self._fd_derivative(point_storage, derivative_order_storage) def eval_batch(self, points: np.ndarray) -> np.ndarray: """Evaluate at multiple points simultaneously. @@ -2256,41 +2292,29 @@ def eval_multi( """ self._check_built() - # If a non-identity ``_dim_order`` is set (e.g., from - # :meth:`with_auto_order` or :meth:`reorder`), permute both the - # input ``point`` and each ``deriv_order`` from user-frame into - # storage frame so the FD machinery (which references - # ``self.domain`` / ``self.n_nodes`` in storage frame) operates - # consistently. Then temporarily neutralize ``_dim_order`` so - # the inner ``self.eval`` calls do not re-permute the - # already-storage-frame points. + # Structural fix (issue #19): permute user-frame coords ONCE into + # storage frame, then call :meth:`_eval_storage_frame` for each + # derivative. The previous implementation temporarily mutated + # ``self._dim_order`` via try/finally to neutralize re-permutation + # inside the inner ``self.eval`` calls, which raced under + # concurrent invocations. No mutation occurs here. canonical = list(range(self.num_dimensions)) if self._dim_order != canonical: - point = [point[self._dim_order[k]] for k in range(self.num_dimensions)] - derivative_orders = [ + point_storage = [ + point[self._dim_order[k]] for k in range(self.num_dimensions) + ] + derivs_storage = [ [d[self._dim_order[k]] for k in range(self.num_dimensions)] for d in derivative_orders ] - saved_dim_order = self._dim_order - self._dim_order = canonical - try: - results = [] - for deriv_order in derivative_orders: - if all(d == 0 for d in deriv_order): - results.append(self.eval(point)) - else: - results.append(self._fd_derivative(point, deriv_order)) - finally: - self._dim_order = saved_dim_order - return results - - results = [] - for deriv_order in derivative_orders: - if all(d == 0 for d in deriv_order): - results.append(self.eval(point)) - else: - results.append(self._fd_derivative(point, deriv_order)) - return results + else: + point_storage = list(point) + derivs_storage = [list(d) for d in derivative_orders] + + return [ + self._eval_storage_frame(point_storage, d_storage) + for d_storage in derivs_storage + ] def _fd_derivative(self, point: List[float], deriv_order: List[int]) -> float: """Compute a single derivative via central finite differences. @@ -2343,34 +2367,48 @@ def _nudge_point(self, point: List[float], d: int, h: float) -> List[float]: return pt def _fd_single_dim(self, point: List[float], d: int, order: int) -> float: - """Central FD for a single dimension.""" + """Central FD for a single dimension. + + Operates entirely in storage frame: ``point`` is assumed already + permuted, ``d`` is a storage-frame dim index, and the inner + evaluations go through :meth:`_eval_storage_frame` to skip the + user→storage permutation. + """ h = self._fd_step(d) pt = self._nudge_point(point, d, h) + zero = [0] * self.num_dimensions if order == 1: pt_plus = list(pt) pt_minus = list(pt) pt_plus[d] += h pt_minus[d] -= h - return (self.eval(pt_plus) - self.eval(pt_minus)) / (2.0 * h) + return ( + self._eval_storage_frame(pt_plus, zero) + - self._eval_storage_frame(pt_minus, zero) + ) / (2.0 * h) elif order == 2: pt_plus = list(pt) pt_minus = list(pt) pt_plus[d] += h pt_minus[d] -= h - f_plus = self.eval(pt_plus) - f_center = self.eval(pt) - f_minus = self.eval(pt_minus) + f_plus = self._eval_storage_frame(pt_plus, zero) + f_center = self._eval_storage_frame(pt, zero) + f_minus = self._eval_storage_frame(pt_minus, zero) return (f_plus - 2.0 * f_center + f_minus) / (h * h) else: raise ValueError(f"Derivative order {order} not supported (use 1 or 2)") def _fd_cross_deriv(self, point: List[float], d1: int, d2: int) -> float: - """Central FD for mixed partial d^2f / dx_{d1} dx_{d2}.""" + """Central FD for mixed partial d^2f / dx_{d1} dx_{d2}. + + Operates entirely in storage frame. + """ h1 = self._fd_step(d1) h2 = self._fd_step(d2) pt = self._nudge_point(point, d1, h1) pt = self._nudge_point(pt, d2, h2) + zero = [0] * self.num_dimensions def make_pt(delta1: float, delta2: float) -> List[float]: p = list(pt) @@ -2378,19 +2416,23 @@ def make_pt(delta1: float, delta2: float) -> List[float]: p[d2] += delta2 return p - f_pp = self.eval(make_pt(+h1, +h2)) - f_pm = self.eval(make_pt(+h1, -h2)) - f_mp = self.eval(make_pt(-h1, +h2)) - f_mm = self.eval(make_pt(-h1, -h2)) + f_pp = self._eval_storage_frame(make_pt(+h1, +h2), zero) + f_pm = self._eval_storage_frame(make_pt(+h1, -h2), zero) + f_mp = self._eval_storage_frame(make_pt(-h1, +h2), zero) + f_mm = self._eval_storage_frame(make_pt(-h1, -h2), zero) return (f_pp - f_pm - f_mp + f_mm) / (4.0 * h1 * h2) def _fd_nested( self, point: List[float], active_dims: List[Tuple[int, int]] ) -> float: - """Nested FD for higher-order or multi-dimensional derivatives.""" + """Nested FD for higher-order or multi-dimensional derivatives. + + Operates entirely in storage frame. + """ # Apply FD one dimension at a time if len(active_dims) == 0: - return self.eval(point) + zero = [0] * self.num_dimensions + return self._eval_storage_frame(point, zero) d, order = active_dims[0] remaining = active_dims[1:] diff --git a/tests/test_v0201_dim_threading.py b/tests/test_v0201_dim_threading.py index c844c8a..c9cf26e 100644 --- a/tests/test_v0201_dim_threading.py +++ b/tests/test_v0201_dim_threading.py @@ -641,3 +641,45 @@ def f(x, _): return x[0] + x[1] + x[2] # Column 2 must be in user-frame dim 2 range: [-3, 3] assert points[:, 2].min() >= -3.0 - 1e-12 assert points[:, 2].max() <= 3.0 + 1e-12 + + +class TestEvalMultiNoDimOrderMutation: + """Issue #19: eval_multi must not mutate self._dim_order via try/finally + (race-prone under concurrent calls). The fix is structural: introduce + _eval_storage_frame and call it from eval_multi after permuting once.""" + + def test_eval_multi_source_does_not_assign_dim_order(self): + """Source-level check: the body of eval_multi contains no + 'self._dim_order = ' assignment.""" + import inspect + from pychebyshev.tensor_train import ChebyshevTT + + source = inspect.getsource(ChebyshevTT.eval_multi) + forbidden = [ + "self._dim_order = ", + "self._dim_order=", + ] + for pat in forbidden: + assert pat not in source, ( + f"eval_multi must not contain assignment '{pat}' " + f"after issue #19 fix" + ) + + def test_eval_multi_correctness_under_reorder(self): + """Behavioral check: eval_multi still returns correct values + for a reordered TT.""" + def f(x, _): return x[0] ** 2 + x[1] - x[2] + tt = ChebyshevTT( + f, num_dimensions=3, domain=[(-1, 1)] * 3, n_nodes=[7, 7, 7], + ) + tt.build(verbose=False) + tt_reordered = tt.reorder([2, 0, 1]) + # eval_multi for value + first-derivative-wrt-x0 (in user frame) + results = tt_reordered.eval_multi( + point=[0.3, 0.4, 0.5], + derivative_orders=[[0, 0, 0], [1, 0, 0]], + ) + # f(0.3, 0.4, 0.5) = 0.09 + 0.4 - 0.5 = -0.01 + # df/dx0(0.3, 0.4, 0.5) = 2*0.3 = 0.6 + assert abs(results[0] - (-0.01)) < 1e-9 + assert abs(results[1] - 0.6) < 1e-9 From 7b959b3977f490a2c383c6d9c4408454cf4ee54b Mon Sep 17 00:00:00 2001 From: Max Zhang Date: Mon, 27 Apr 2026 12:31:01 -0400 Subject: [PATCH 07/16] v0.21.1: vectorize _calculus._optimize_1d candidate evaluation MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Replace the Python list comprehension over candidates with a single vectorized barycentric evaluation (BLAS GEMV pattern). _optimize_1d is the canonical 1-D optimizer used by all four classes' minimize/maximize methods; v0.21 routed Slider/TT min/max through it, making the perf hoist worthwhile. Bit-identical results — non-regression tests pass. Co-Authored-By: Claude Sonnet 4.6 --- src/pychebyshev/_calculus.py | 23 ++++++++++++++++------- tests/test_calculus.py | 36 ++++++++++++++++++++++++++++++++++++ 2 files changed, 52 insertions(+), 7 deletions(-) diff --git a/src/pychebyshev/_calculus.py b/src/pychebyshev/_calculus.py index a872f6f..400cd20 100644 --- a/src/pychebyshev/_calculus.py +++ b/src/pychebyshev/_calculus.py @@ -266,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 @@ -278,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]) diff --git a/tests/test_calculus.py b/tests/test_calculus.py index 054ae9e..0850ae6 100644 --- a/tests/test_calculus.py +++ b/tests/test_calculus.py @@ -804,3 +804,39 @@ def f(x, _): return abs(x[0]) result = sp.integrate(bounds=(0.0, 0.5)) expected = 0.125 # int_0^0.5 x dx = 0.5^2/2 assert abs(result - expected) < 1e-10, f"Expected {expected}, got {result}" + + +# ====================================================================== +# TestOptimize1DVectorization +# ====================================================================== + +class TestOptimize1DVectorization: + """v0.21.1: _optimize_1d should produce bit-identical results before + and after vectorization. Non-regression tests.""" + + def test_optimize_1d_quadratic_min_unchanged(self): + """min of (x-0.3)^2 on [-1,1] is 0 at x=0.3.""" + def f(x, _): return (x[0] - 0.3) ** 2 + cheb = ChebyshevApproximation(f, 1, [(-1, 1)], [10]) + cheb.build(verbose=False) + val, loc = cheb.minimize() + assert abs(val - 0.0) < 1e-10 + assert abs(loc - 0.3) < 1e-10 + + def test_optimize_1d_quadratic_max_unchanged(self): + """max of -(x-0.3)^2 + 5 on [-1,1] is 5 at x=0.3.""" + def f(x, _): return -(x[0] - 0.3) ** 2 + 5.0 + cheb = ChebyshevApproximation(f, 1, [(-1, 1)], [10]) + cheb.build(verbose=False) + val, loc = cheb.maximize() + assert abs(val - 5.0) < 1e-10 + assert abs(loc - 0.3) < 1e-10 + + def test_optimize_1d_endpoint_min(self): + """Linear x on [0,1] — min at x=0 (endpoint, no critical points).""" + def f(x, _): return x[0] + cheb = ChebyshevApproximation(f, 1, [(0, 1)], [5]) + cheb.build(verbose=False) + val, loc = cheb.minimize() + assert abs(val - 0.0) < 1e-10 + assert abs(loc - 0.0) < 1e-10 From ef44ed45284b8acb86bd087ae9bdf292e2bc875b Mon Sep 17 00:00:00 2001 From: Max Zhang Date: Mon, 27 Apr 2026 12:38:00 -0400 Subject: [PATCH 08/16] v0.21.1: hoist derivative matrix application out of vectorized_eval_batch loop MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The differentiation matrix (D_T) matmul is point-independent — applying it once to self.tensor_values before the per-point loop, then doing only the barycentric reduction inside the loop, produces identical results with substantial speedup for batch derivative evaluations. Each D_T[d] is applied along axis d of the full tensor (via moveaxis to last, matmul, moveaxis back) to correctly match the per-dim interleaved semantics of vectorized_eval. ChebyshevSpline.eval_batch (which delegates per-piece) inherits the speedup. --- src/pychebyshev/barycentric.py | 33 ++++++++++++++++++++++++- tests/test_barycentric.py | 45 ++++++++++++++++++++++++++++++++++ 2 files changed, 77 insertions(+), 1 deletion(-) diff --git a/src/pychebyshev/barycentric.py b/src/pychebyshev/barycentric.py index 98b5e50..a9b1fcf 100644 --- a/src/pychebyshev/barycentric.py +++ b/src/pychebyshev/barycentric.py @@ -976,10 +976,41 @@ def vectorized_eval_batch( derivative_order = self._resolve_derivative_args( derivative_order, derivative_id ) + if self.tensor_values is None: + raise RuntimeError("Call build() first") + + # Hoist: apply all derivative-matrix matmuls once — they are + # point-independent (depend only on the grid, not on query coords). + # In vectorized_eval the D_T pass for dim d acts on the last axis + # of a tensor that has already had dims d+1..D-1 reduced away. + # Equivalently, on the full tensor it acts along axis d. We move + # axis d to the last position, apply the matmul, then move it back. + tensor_with_derivs = self.tensor_values + for d in range(self.num_dimensions - 1, -1, -1): + deriv = derivative_order[d] + if deriv > 0: + D_T = self.diff_matrices[d].T + for _ in range(deriv): + arr = np.moveaxis(tensor_with_derivs, d, -1) + arr = arr @ D_T + tensor_with_derivs = np.moveaxis(arr, -1, d) + + # Per-point: only the barycentric reduction (no derivative passes). + _matmul = self._matmul_last_axis N = points.shape[0] results = np.empty(N) for i in range(N): - results[i] = self.vectorized_eval(points[i], derivative_order=derivative_order) + current = tensor_with_derivs + for d in range(self.num_dimensions - 1, -1, -1): + x = points[i, d] + diff = x - self.nodes[d] + exact = np.where(np.abs(diff) < 1e-14)[0] + if len(exact) > 0: + current = current[..., exact[0]] + else: + w_over_diff = self.weights[d] / diff + current = _matmul(current, w_over_diff) / np.sum(w_over_diff) + results[i] = float(current) return results def vectorized_eval_multi( diff --git a/tests/test_barycentric.py b/tests/test_barycentric.py index 138e1fc..9a843c7 100644 --- a/tests/test_barycentric.py +++ b/tests/test_barycentric.py @@ -501,3 +501,48 @@ def test_per_dim_requires_build(self): ) with pytest.raises(RuntimeError, match="build"): cheb._error_estimate_per_dim() + + +class TestVectorizedEvalBatchDerivativeHoist: + """v0.21.1: vectorized_eval_batch should produce identical results + before and after the derivative-matrix hoist. Non-regression tests.""" + + def test_batch_eval_value_only(self): + """Batch eval with derivative_order=[0] returns same as per-point eval.""" + def f(x, _): return x[0] ** 2 + 0.5 * x[0] + cheb = ChebyshevApproximation(f, 1, [(-1, 1)], [10]) + cheb.build(verbose=False) + points = np.array([[-0.7], [-0.3], [0.0], [0.4], [0.8]]) + batch = cheb.vectorized_eval_batch(points, derivative_order=[0]) + per_point = np.array([ + cheb.vectorized_eval([float(p[0])], derivative_order=[0]) + for p in points + ]) + np.testing.assert_array_almost_equal(batch, per_point, decimal=12) + + def test_batch_eval_with_derivative(self): + """Batch eval with derivative_order=[1] returns same as per-point.""" + def f(x, _): return x[0] ** 2 + 0.5 * x[0] + cheb = ChebyshevApproximation(f, 1, [(-1, 1)], [10]) + cheb.build(verbose=False) + points = np.array([[-0.7], [-0.3], [0.0], [0.4], [0.8]]) + batch = cheb.vectorized_eval_batch(points, derivative_order=[1]) + per_point = np.array([ + cheb.vectorized_eval([float(p[0])], derivative_order=[1]) + for p in points + ]) + np.testing.assert_array_almost_equal(batch, per_point, decimal=12) + + def test_batch_eval_2d_with_derivative(self): + """2D function, mixed derivative — batch matches per-point.""" + def f(x, _): return x[0] ** 2 + x[1] ** 2 + cheb = ChebyshevApproximation(f, 2, [(-1, 1), (-1, 1)], [8, 8]) + cheb.build(verbose=False) + points = np.array([[-0.5, 0.3], [0.2, -0.4], [0.7, 0.1]]) + batch = cheb.vectorized_eval_batch(points, derivative_order=[1, 0]) + per_point = np.array([ + cheb.vectorized_eval([float(p[0]), float(p[1])], + derivative_order=[1, 0]) + for p in points + ]) + np.testing.assert_array_almost_equal(batch, per_point, decimal=12) From 6cdc5c646e470fd8ad91d3e58764ed756975e74e Mon Sep 17 00:00:00 2001 From: Max Zhang Date: Mon, 27 Apr 2026 12:43:47 -0400 Subject: [PATCH 09/16] v0.21.1: add comparison/demo script MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Demonstrates each fix in the TT _dim_order cluster with a non-uniform domain and explicit reorder() — the configuration that v0.20.1/v0.21.0 test coverage masked. Native TT sobol_indices deferred to v0.22; not included in this demo. --- compare_v0211_dim_cluster.py | 114 +++++++++++++++++++++++++++++++++++ 1 file changed, 114 insertions(+) create mode 100644 compare_v0211_dim_cluster.py diff --git a/compare_v0211_dim_cluster.py b/compare_v0211_dim_cluster.py new file mode 100644 index 0000000..dbc7bf4 --- /dev/null +++ b/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() From 609e75921339b815e4ec6f4f7c4aac521df2bfc0 Mon Sep 17 00:00:00 2001 From: Max Zhang Date: Mon, 27 Apr 2026 12:48:20 -0400 Subject: [PATCH 10/16] v0.21.1: docs + release housekeeping Add non-uniform-domain validation note to user-guide. Bump version to 0.21.1. CHANGELOG entry covering all 6 fixes + 2 perf wins, with sobol_indices deferred to v0.22. CLAUDE.md updated with new test count and v0.21.1 architecture line. --- CHANGELOG.md | 51 +++++++++++++++++++++++++++++++++++++ CLAUDE.md | 13 +++++++--- docs/user-guide/calculus.md | 2 ++ pyproject.toml | 2 +- src/pychebyshev/_version.py | 2 +- 5 files changed, 65 insertions(+), 5 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 02e2e35..9f005e0 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,57 @@ 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 + +### 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. +- `ChebyshevTT.sobol_indices()` deferred to v0.22 (native TT + contraction implementation exceeded v0.21.1 time budget). + ## [0.21.0] - 2026-04-27 ### Added diff --git a/CLAUDE.md b/CLAUDE.md index 439624c..3e5aa66 100644 --- a/CLAUDE.md +++ b/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. `ChebyshevTT.sobol_indices` deferred to v0.22. ### Benchmark Scripts (project root) @@ -129,12 +135,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. diff --git a/docs/user-guide/calculus.md b/docs/user-guide/calculus.md index a7a41c2..aa7b378 100644 --- a/docs/user-guide/calculus.md +++ b/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 diff --git a/pyproject.toml b/pyproject.toml index c51f9bf..2ac2e98 100644 --- a/pyproject.toml +++ b/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"} diff --git a/src/pychebyshev/_version.py b/src/pychebyshev/_version.py index 6a726d8..76f2458 100644 --- a/src/pychebyshev/_version.py +++ b/src/pychebyshev/_version.py @@ -1 +1 @@ -__version__ = "0.21.0" +__version__ = "0.21.1" From f35ce644ec818939e426b369dc8ba2c4bbb3fc9b Mon Sep 17 00:00:00 2001 From: Max Zhang Date: Mon, 27 Apr 2026 13:05:16 -0400 Subject: [PATCH 11/16] =?UTF-8?q?v0.21.1:=20address=20final=20review=20?= =?UTF-8?q?=E2=80=94=20extend=20np.allclose=20fix=20to=20TT=20internal=20s?= =?UTF-8?q?ites?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Final cumulative review (M1) noted that ChebyshevTT.inner_product and the binary-algebra _check_compatible_tt path also use exact == on self.domain, paralleling the issue #22 fix in _algebra._check_compatible. Apply the same np.allclose pattern to both TT sites for consistency with the v0.21.1 fix. (M2) Add compare_v0211_dim_cluster.py to the CLAUDE.md benchmark roster. --- CLAUDE.md | 1 + src/pychebyshev/tensor_train.py | 10 ++++++++-- 2 files changed, 9 insertions(+), 2 deletions(-) diff --git a/CLAUDE.md b/CLAUDE.md index 3e5aa66..da36da9 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -111,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/`) diff --git a/src/pychebyshev/tensor_train.py b/src/pychebyshev/tensor_train.py index f74332b..23268c6 100644 --- a/src/pychebyshev/tensor_train.py +++ b/src/pychebyshev/tensor_train.py @@ -1472,7 +1472,10 @@ def inner_product(self, other: "ChebyshevTT") -> float: f"other must be a ChebyshevTT, got {type(other).__name__}" ) other._check_built() - if self.domain != other.domain: + if not np.allclose( + np.asarray(self.domain, dtype=float), + np.asarray(other.domain, dtype=float), + ): raise ValueError( "inner_product requires matching domains; " f"got {self.domain} vs {other.domain}" @@ -3252,7 +3255,10 @@ def _check_compatible_tt(self, other) -> None: raise ValueError( f"n_nodes mismatch: {self.n_nodes} vs {other.n_nodes}" ) - if self.domain != other.domain: + if not np.allclose( + np.asarray(self.domain, dtype=float), + np.asarray(other.domain, dtype=float), + ): raise ValueError( f"domain mismatch: {self.domain} vs {other.domain}" ) From aa21e6f9e070b0f9983ee68ecbd798af287bc219 Mon Sep 17 00:00:00 2001 From: Max Zhang Date: Mon, 27 Apr 2026 13:59:38 -0400 Subject: [PATCH 12/16] v0.21.1: add ChebyshevTT.sobol_indices() (native TT contraction) Compute first-order + total-order Sobol indices by contracting through the TT coefficient cores in O(d * n * r^2) time. Cross-validated against ChebyshevApproximation.sobol_indices() on dense reconstructions. Closes the v0.21.1 sobol_indices parity item (originally deferred but attempted on user request). --- src/pychebyshev/_sensitivity.py | 107 +++++++++++++++++++++++++ src/pychebyshev/tensor_train.py | 47 +++++++++++ tests/test_tensor_train.py | 67 ++++++++++++++++ tests/test_v020_adaptive_refinement.py | 59 ++++++++++++++ 4 files changed, 280 insertions(+) diff --git a/src/pychebyshev/_sensitivity.py b/src/pychebyshev/_sensitivity.py index e94b1d0..dedb627 100644 --- a/src/pychebyshev/_sensitivity.py +++ b/src/pychebyshev/_sensitivity.py @@ -138,3 +138,110 @@ def _compute_sobol_from_coeffs( }, "variance": variance, } + + +def _compute_sobol_from_tt_cores(cores: list) -> dict: + """Compute first-order + total-order Sobol indices from TT coefficient cores. + + Mathematically equivalent to ``_compute_sobol_from_coeffs`` applied to the + dense coefficient tensor, but contracts through TT cores in coefficient + space with cost O(d * n * r^2) instead of O(n^d). + + Parameters + ---------- + cores : list of np.ndarray + TT cores in Chebyshev coefficient space. Each core has shape + (r_{k-1}, n_k, r_k) with cores[0] starting from r_0=1 and + cores[-1] ending at r_d=1. + + Returns + ------- + dict + Same shape as ``_compute_sobol_from_coeffs``: + ``{"first_order": {d: index}, "total_order": {d: index}, "variance": float}`` + where keys are storage-frame dim indices (0..len(cores)-1). + """ + d = len(cores) + pi = float(np.pi) + n_per_dim = [c.shape[1] for c in cores] + + # Per-dim Chebyshev inner-product weights: [pi, pi/2, pi/2, ..., pi/2] + w_full = [] + for n_k in n_per_dim: + w = np.full(n_k, pi / 2.0) + w[0] = pi + w_full.append(w) + + # ---- total_weighted_squared = sum over all alpha of + # coeffs[alpha]^2 * prod_k w_full[k][alpha_k] + M = np.array([[1.0]]) + for k in range(d): + A = cores[k] + Aw = A * w_full[k][None, :, None] + M = np.einsum("ij,ipa,jpb->ab", M, Aw, A) + total_weighted_squared = float(M[0, 0]) + + # ---- constant term c_0 (alpha = 0) -- contract along all-zero slices + v = np.array([1.0]) + for k in range(d): + v = v @ cores[k][:, 0, :] + c_0 = float(v[0]) + constant_weighted_squared = (c_0 ** 2) * (pi ** d) + + variance = total_weighted_squared - constant_weighted_squared + + if variance <= 0: + return { + "first_order": {j: 0.0 for j in range(d)}, + "total_order": {j: 0.0 for j in range(d)}, + "variance": float(max(variance, 0.0)), + } + + first_order_energy = {} + total_order_energy = {} + + for j in range(d): + # ---- first-order energy[j]: alpha_j >= 1 AND alpha_k = 0 for k != j + # left boundary: row vector formed by chaining cores[k][:, 0, :] for k < j + left = np.array([1.0]) + for k in range(j): + left = left @ cores[k][:, 0, :] + # left has shape (r_j,) + + # right boundary: column vector formed by chaining cores[k][:, 0, :] for k > j + right = np.array([1.0]) + for k in range(d - 1, j, -1): + right = cores[k][:, 0, :] @ right + # right has shape (r_{j+1},) + + G_j = cores[j] + sum_squared = 0.0 + for m in range(1, n_per_dim[j]): + coef_m = float(left @ G_j[:, m, :] @ right) + sum_squared += coef_m * coef_m + weight_first = (pi / 2.0) * (pi ** (d - 1)) + first_order_energy[j] = sum_squared * weight_first + + # ---- total-order energy[j]: alpha_j >= 1 (other dims unrestricted) + # = total_weighted_squared - sum_{alpha_j = 0} of + # coeffs[alpha]^2 * prod_k w_full[k][alpha_k] + # The "alpha_j = 0" sum: replace cores[j] with cores[j][:, 0:1, :] + # and use w_k = [pi] for dim j. + M_j = np.array([[1.0]]) + for k in range(d): + if k == j: + A_slice = cores[k][:, 0:1, :] + w_k = np.array([pi]) + else: + A_slice = cores[k] + w_k = w_full[k] + Aw_slice = A_slice * w_k[None, :, None] + M_j = np.einsum("ij,ipa,jpb->ab", M_j, Aw_slice, A_slice) + sum_alpha_j_zero_weighted = float(M_j[0, 0]) + total_order_energy[j] = total_weighted_squared - sum_alpha_j_zero_weighted + + return { + "first_order": {j: first_order_energy[j] / variance for j in range(d)}, + "total_order": {j: total_order_energy[j] / variance for j in range(d)}, + "variance": float(variance), + } diff --git a/src/pychebyshev/tensor_train.py b/src/pychebyshev/tensor_train.py index 23268c6..2bd21ae 100644 --- a/src/pychebyshev/tensor_train.py +++ b/src/pychebyshev/tensor_train.py @@ -2820,6 +2820,53 @@ def clone(self) -> "ChebyshevTT": import copy return copy.deepcopy(self) + def sobol_indices(self) -> dict: + """Compute first-order + total-order Sobol sensitivity indices. + + Returns + ------- + dict + ``{"first_order": {dim: index}, "total_order": {dim: index}, + "variance": float}`` -- keyed by user-frame dim indices. + + Raises + ------ + RuntimeError + If ``build()`` has not been called. + + Notes + ----- + Computed natively by contracting through the TT coefficient cores; + no dense materialization. Cost O(d * n * r^2) per dim. Mathematically + equivalent to ``ChebyshevApproximation.sobol_indices()`` on the + dense version of the same function. + + Under non-identity ``_dim_order`` (after ``with_auto_order`` / + ``reorder``), result keys are user-frame dim indices. + """ + if not self._built: + raise RuntimeError("Call build() first") + + from pychebyshev._sensitivity import _compute_sobol_from_tt_cores + + # Helper returns dict keyed by storage-frame indices (0..d-1) + storage = _compute_sobol_from_tt_cores(self._coeff_cores) + + # Translate keys to user-frame: storage position s holds + # original-dim self._dim_order[s]. + user_first = {} + user_total = {} + for s in range(self.num_dimensions): + user_d = self._dim_order[s] + user_first[user_d] = storage["first_order"][s] + user_total[user_d] = storage["total_order"][s] + + return { + "first_order": user_first, + "total_order": user_total, + "variance": storage["variance"], + } + @classmethod def from_values( cls, diff --git a/tests/test_tensor_train.py b/tests/test_tensor_train.py index 926e31d..cbe9a0a 100644 --- a/tests/test_tensor_train.py +++ b/tests/test_tensor_train.py @@ -986,3 +986,70 @@ def f(x, data=None): assert err_after <= err_before * 1.1 + 1e-14, ( f"completion increased error: before={err_before:.3e}, after={err_after:.3e}" ) + + +# ====================================================================== +# v0.21.1: Sobol indices via native TT contraction +# ====================================================================== + + +class TestSobolFromTTCores: + """v0.21.1: _compute_sobol_from_tt_cores must produce identical Sobol + indices to _compute_sobol_from_coeffs on the dense coefficient tensor.""" + + def test_separable_2d(self): + """Separable f(x, y) = x + y -- equal main-effect indices ~50/50.""" + from pychebyshev import ChebyshevApproximation + from pychebyshev._sensitivity import _compute_sobol_from_tt_cores + + def f(x, _): + return x[0] + x[1] + + cheb = ChebyshevApproximation(f, 2, [(-1, 1), (-1, 1)], [7, 7]) + cheb.build(verbose=False) + ref = cheb.sobol_indices() + tt = ChebyshevTT(f, num_dimensions=2, domain=[(-1, 1), (-1, 1)], n_nodes=[7, 7]) + tt.build(verbose=False) + tt_indices = _compute_sobol_from_tt_cores(tt._coeff_cores) + # Check first_order and total_order dicts agree + for key in ["first_order", "total_order"]: + assert set(tt_indices[key].keys()) == set(ref[key].keys()) + for d in tt_indices[key]: + assert abs(tt_indices[key][d] - ref[key][d]) < 1e-9, ( + f"{key}[{d}]: TT={tt_indices[key][d]}, Approx={ref[key][d]}" + ) + # Variance: should match too (same coefficients) + assert abs(tt_indices["variance"] - ref["variance"]) < 1e-9 * abs(ref["variance"]) + + def test_dominant_dim_3d(self): + """3D function dominated by dim 0 -- first_order[0] near 1.0.""" + from pychebyshev._sensitivity import _compute_sobol_from_tt_cores + + def f(x, _): + return x[0] + 0.01 * x[1] + 0.01 * x[2] + + tt = ChebyshevTT(f, num_dimensions=3, domain=[(-1, 1)] * 3, n_nodes=[8, 8, 8]) + tt.build(verbose=False) + idx = _compute_sobol_from_tt_cores(tt._coeff_cores) + # Dim 0 should dominate + assert idx["first_order"][0] > 0.99 + assert idx["first_order"][1] < 0.01 + assert idx["first_order"][2] < 0.01 + + def test_quadratic_3d_parity(self): + """3D quadratic -- match Approximation indices to ~1e-8.""" + from pychebyshev import ChebyshevApproximation + from pychebyshev._sensitivity import _compute_sobol_from_tt_cores + + def f(x, _): + return x[0] ** 2 + 0.5 * x[1] ** 2 + 0.25 * x[2] ** 2 + + cheb = ChebyshevApproximation(f, 3, [(-1, 1)] * 3, [9, 9, 9]) + cheb.build(verbose=False) + ref = cheb.sobol_indices() + tt = ChebyshevTT(f, num_dimensions=3, domain=[(-1, 1)] * 3, n_nodes=[9, 9, 9]) + tt.build(verbose=False) + idx = _compute_sobol_from_tt_cores(tt._coeff_cores) + for key in ["first_order", "total_order"]: + for d in idx[key]: + assert abs(idx[key][d] - ref[key][d]) < 1e-8 diff --git a/tests/test_v020_adaptive_refinement.py b/tests/test_v020_adaptive_refinement.py index 0e98d61..3526d35 100644 --- a/tests/test_v020_adaptive_refinement.py +++ b/tests/test_v020_adaptive_refinement.py @@ -515,3 +515,62 @@ def test_sobol_rejects_inf_coeffs(self): coeffs = np.array([1.0, float("inf"), 0.0]) with pytest.raises(ValueError, match="NaN or Inf"): _compute_sobol_from_coeffs(coeffs, num_dimensions=1) + + +# ============================================================================ +# v0.21.1: ChebyshevTT.sobol_indices parity (native TT contraction) +# ============================================================================ + +class TestTTSobolParity: + """v0.21.1: ChebyshevTT.sobol_indices must match + ChebyshevApproximation.sobol_indices on the same function.""" + + def test_separable_2d(self): + def f(x, _): + return x[0] + x[1] + + cheb = ChebyshevApproximation(f, 2, [(-1, 1), (-1, 1)], [7, 7]) + tt = ChebyshevTT(f, num_dimensions=2, domain=[(-1, 1), (-1, 1)], n_nodes=[7, 7]) + cheb.build(verbose=False) + tt.build(verbose=False) + cheb_idx = cheb.sobol_indices() + tt_idx = tt.sobol_indices() + for key in ["first_order", "total_order"]: + for d in cheb_idx[key]: + assert abs(cheb_idx[key][d] - tt_idx[key][d]) < 1e-9 + + def test_quadratic_3d(self): + def f(x, _): + return x[0] ** 2 + 0.5 * x[1] ** 2 + 0.25 * x[2] ** 2 + + cheb = ChebyshevApproximation(f, 3, [(-1, 1)] * 3, [10, 10, 10]) + tt = ChebyshevTT(f, num_dimensions=3, domain=[(-1, 1)] * 3, n_nodes=[10, 10, 10]) + cheb.build(verbose=False) + tt.build(verbose=False) + cheb_idx = cheb.sobol_indices() + tt_idx = tt.sobol_indices() + for key in ["first_order", "total_order"]: + for d in cheb_idx[key]: + assert abs(cheb_idx[key][d] - tt_idx[key][d]) < 1e-7 + + def test_user_frame_after_reorder(self): + """sobol_indices keys must be user-frame, even after reorder.""" + def f(x, _): + return x[0] + 0.01 * x[1] + 0.01 * x[2] + + tt = ChebyshevTT(f, num_dimensions=3, domain=[(-1, 1)] * 3, n_nodes=[8, 8, 8]) + tt.build(verbose=False) + tt_reordered = tt.reorder([2, 0, 1]) + idx = tt_reordered.sobol_indices() + # Dim 0 (user-frame) should still dominate + assert idx["first_order"][0] > 0.99 + assert idx["first_order"][1] < 0.01 + assert idx["first_order"][2] < 0.01 + + def test_before_build_raises(self): + def f(x, _): + return x[0] + + tt = ChebyshevTT(f, num_dimensions=1, domain=[(-1, 1)], n_nodes=[5]) + with pytest.raises(RuntimeError, match="build"): + tt.sobol_indices() From 88eba583765cde114d994cd484b04b340bb29b17 Mon Sep 17 00:00:00 2001 From: Max Zhang Date: Mon, 27 Apr 2026 14:00:06 -0400 Subject: [PATCH 13/16] =?UTF-8?q?v0.21.1:=20CHANGELOG/CLAUDE.md=20?= =?UTF-8?q?=E2=80=94=20sobol=5Findices=20shipped,=20not=20deferred?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- CHANGELOG.md | 10 ++++++++-- CLAUDE.md | 2 +- 2 files changed, 9 insertions(+), 3 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 9f005e0..6a7c837 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,14 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [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 @@ -53,8 +61,6 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 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. -- `ChebyshevTT.sobol_indices()` deferred to v0.22 (native TT - contraction implementation exceeded v0.21.1 time budget). ## [0.21.0] - 2026-04-27 diff --git a/CLAUDE.md b/CLAUDE.md index da36da9..2625000 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -89,7 +89,7 @@ The installable package. Public classes: `ChebyshevApproximation`, `ChebyshevSpl `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. `ChebyshevTT.sobol_indices` deferred to v0.22. + hoist. Adds `ChebyshevTT.sobol_indices()` parity (native TT contraction). ### Benchmark Scripts (project root) From 3139b5dbd8752b5336165948d6fc2a9661f634c3 Mon Sep 17 00:00:00 2001 From: Max Zhang Date: Mon, 27 Apr 2026 14:38:31 -0400 Subject: [PATCH 14/16] tests: add non-uniform-domain + reorder sobol parity test (Issue 1) Appends test_non_uniform_domain_after_reorder to TestTTSobolParity, covering the exact configuration that masked v0.21.0 _dim_order bugs: non-uniform per-dim domains [(-1,1), (-2,2), (-3,3)] combined with a non-identity _dim_order (reorder([2,0,1])). Verifies TT sobol keys stay in user-frame after reorder, matching ChebyshevApproximation reference. Co-Authored-By: Claude Sonnet 4.6 --- tests/test_v020_adaptive_refinement.py | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/tests/test_v020_adaptive_refinement.py b/tests/test_v020_adaptive_refinement.py index 3526d35..d89f85f 100644 --- a/tests/test_v020_adaptive_refinement.py +++ b/tests/test_v020_adaptive_refinement.py @@ -574,3 +574,26 @@ def f(x, _): tt = ChebyshevTT(f, num_dimensions=1, domain=[(-1, 1)], n_nodes=[5]) with pytest.raises(RuntimeError, match="build"): tt.sobol_indices() + + def test_non_uniform_domain_after_reorder(self): + """Cross-validate TT sobol against Approximation under non-uniform + per-dim domains AND non-identity _dim_order — the configuration that + v0.21.0 cluster bugs hid behind uniform-domain test coverage.""" + def f(x, _): return (x[0] - 0.2) ** 2 + 0.5 * (x[1] / 2.0) ** 2 + 0.25 * (x[2] / 3.0) ** 2 + cheb = ChebyshevApproximation(f, 3, [(-1, 1), (-2, 2), (-3, 3)], [10, 10, 10]) + tt = ChebyshevTT( + f, num_dimensions=3, + domain=[(-1, 1), (-2, 2), (-3, 3)], n_nodes=[10, 10, 10], + ) + cheb.build(verbose=False) + tt.build(verbose=False) + tt_reordered = tt.reorder([2, 0, 1]) + # Reference (canonical Approximation) + ref = cheb.sobol_indices() + # Reordered TT — keys must still be in user frame + idx = tt_reordered.sobol_indices() + for key in ["first_order", "total_order"]: + for d in ref[key]: + assert abs(ref[key][d] - idx[key][d]) < 1e-7, ( + f"{key}[{d}]: cheb={ref[key][d]}, tt_reordered={idx[key][d]}" + ) From 4b8a825cc6815350853fceb29107233a2d0b98ab Mon Sep 17 00:00:00 2001 From: Max Zhang Date: Mon, 27 Apr 2026 14:38:38 -0400 Subject: [PATCH 15/16] refactor: extract _apply_derivative_passes helper to barycentric.py (Issue 2) Adds private _apply_derivative_passes(tensor, derivative_order) that applies differentiation-matrix passes along each axis of the full coefficient tensor using np.moveaxis. Used by vectorized_eval_batch to replace the duplicated ~10-line derivative-hoist block, making the batch path self-documenting. vectorized_eval inlines its own per-axis logic (reduces axes as it goes, so the full-tensor moveaxis path doesn't apply there) and is left unchanged to avoid unneeded churn. Co-Authored-By: Claude Sonnet 4.6 --- src/pychebyshev/barycentric.py | 59 ++++++++++++++++++++++++++-------- 1 file changed, 46 insertions(+), 13 deletions(-) diff --git a/src/pychebyshev/barycentric.py b/src/pychebyshev/barycentric.py index a9b1fcf..d92ebe4 100644 --- a/src/pychebyshev/barycentric.py +++ b/src/pychebyshev/barycentric.py @@ -948,6 +948,47 @@ def vectorized_eval( return float(current) + def _apply_derivative_passes( + self, tensor: np.ndarray, derivative_order: List[int] | None + ) -> np.ndarray: + """Apply differentiation-matrix passes to the full coefficient tensor. + + For each dimension ``d`` with ``derivative_order[d] > 0``, applies the + transposed differentiation matrix ``D_T`` along axis ``d`` the required + number of times. Operates on the **full** tensor (no axes reduced), + using ``np.moveaxis`` to expose the target axis as the last axis for + efficient matmul. + + This helper is designed for use in :meth:`vectorized_eval_batch`, where + the derivative passes are hoisted outside the per-point loop. + :meth:`vectorized_eval` inlines the same logic per-dimension as it + reduces axes incrementally and does not call this helper. + + Parameters + ---------- + tensor : np.ndarray + Full coefficient tensor of shape matching ``self.tensor_values``. + derivative_order : list of int or None + Derivative order per dimension. ``None`` is treated as all-zeros. + + Returns + ------- + np.ndarray + Tensor with all derivative passes applied, same shape as input. + """ + if derivative_order is None: + return tensor + result = tensor + for d in range(self.num_dimensions - 1, -1, -1): + deriv = derivative_order[d] + if deriv > 0: + D_T = self.diff_matrices[d].T + for _ in range(deriv): + arr = np.moveaxis(result, d, -1) + arr = arr @ D_T + result = np.moveaxis(arr, -1, d) + return result + def vectorized_eval_batch( self, points: np.ndarray, @@ -981,19 +1022,11 @@ def vectorized_eval_batch( # Hoist: apply all derivative-matrix matmuls once — they are # point-independent (depend only on the grid, not on query coords). - # In vectorized_eval the D_T pass for dim d acts on the last axis - # of a tensor that has already had dims d+1..D-1 reduced away. - # Equivalently, on the full tensor it acts along axis d. We move - # axis d to the last position, apply the matmul, then move it back. - tensor_with_derivs = self.tensor_values - for d in range(self.num_dimensions - 1, -1, -1): - deriv = derivative_order[d] - if deriv > 0: - D_T = self.diff_matrices[d].T - for _ in range(deriv): - arr = np.moveaxis(tensor_with_derivs, d, -1) - arr = arr @ D_T - tensor_with_derivs = np.moveaxis(arr, -1, d) + # Delegates to _apply_derivative_passes which operates on the full + # tensor using moveaxis to target each axis in turn. + tensor_with_derivs = self._apply_derivative_passes( + self.tensor_values, derivative_order + ) # Per-point: only the barycentric reduction (no derivative passes). _matmul = self._matmul_last_axis From 0c5ba8f25e2c22cd83956e5277ae06820cc6d852 Mon Sep 17 00:00:00 2001 From: Max Zhang Date: Mon, 27 Apr 2026 14:38:45 -0400 Subject: [PATCH 16/16] perf: cache left/right partial contractions for O(d*n*r^2) total-order sobol (Issue 3) Replaces d independent full self-inner-product loops for total-order computation with precomputed L[k] (left partial, dims 0..k-1) and R[k] (right partial, dims k..d-1) matrices. For each dim j, combines L[j], the alpha_j=0 core slice, and R[j+1] via a single einsum instead of re-running all d cores. Reduces from O(d^2*n*r^2) to O(d*n*r^2). Existing TestSobolFromTTCores and TestTTSobolParity tests continue to pass at 1e-9 (2-D) and 1e-7 (3-D) tolerances. Co-Authored-By: Claude Sonnet 4.6 --- src/pychebyshev/_sensitivity.py | 53 +++++++++++++++++++++++---------- 1 file changed, 38 insertions(+), 15 deletions(-) diff --git a/src/pychebyshev/_sensitivity.py b/src/pychebyshev/_sensitivity.py index dedb627..9f73a5c 100644 --- a/src/pychebyshev/_sensitivity.py +++ b/src/pychebyshev/_sensitivity.py @@ -197,6 +197,37 @@ def _compute_sobol_from_tt_cores(cores: list) -> dict: "variance": float(max(variance, 0.0)), } + # ---- Precompute left and right partial self-inner-product matrices + # to reduce total-order computation from O(d^2 * n * r^2) to + # O(d * n * r^2). + # + # L[k] has shape (r_k, r_k): partial contraction of dims 0..k-1. + # R[k] has shape (r_k, r_k): partial contraction of dims k..d-1. + # + # L[0] = identity(1); L[k+1] = einsum("ij,ipa,jpb->ab", L[k], Aw_k, A_k) + # R[d] = identity(1); R[k] = einsum("ab,ipa,jpb->ij", R[k+1], Aw_k, A_k) + # + # For total-order[j]: the "alpha_j = 0" weighted sum decomposes as + # L[j] x cores[j][:,0,:] x R[j+1] (contracted with itself and scaled by pi). + # sum_alpha_j_zero = pi * einsum("ij,ia,jb,ab->", L[j], c_j0, c_j0, R[j+1]) + # where c_j0 = cores[j][:, 0, :]. + + # Build L[0..d] + L = [None] * (d + 1) + L[0] = np.array([[1.0]]) + for k in range(d): + A_k = cores[k] + Aw_k = A_k * w_full[k][None, :, None] + L[k + 1] = np.einsum("ij,ipa,jpb->ab", L[k], Aw_k, A_k) + + # Build R[0..d] + R = [None] * (d + 1) + R[d] = np.array([[1.0]]) + for k in range(d - 1, -1, -1): + A_k = cores[k] + Aw_k = A_k * w_full[k][None, :, None] + R[k] = np.einsum("ab,ipa,jpb->ij", R[k + 1], Aw_k, A_k) + first_order_energy = {} total_order_energy = {} @@ -223,21 +254,13 @@ def _compute_sobol_from_tt_cores(cores: list) -> dict: first_order_energy[j] = sum_squared * weight_first # ---- total-order energy[j]: alpha_j >= 1 (other dims unrestricted) - # = total_weighted_squared - sum_{alpha_j = 0} of - # coeffs[alpha]^2 * prod_k w_full[k][alpha_k] - # The "alpha_j = 0" sum: replace cores[j] with cores[j][:, 0:1, :] - # and use w_k = [pi] for dim j. - M_j = np.array([[1.0]]) - for k in range(d): - if k == j: - A_slice = cores[k][:, 0:1, :] - w_k = np.array([pi]) - else: - A_slice = cores[k] - w_k = w_full[k] - Aw_slice = A_slice * w_k[None, :, None] - M_j = np.einsum("ij,ipa,jpb->ab", M_j, Aw_slice, A_slice) - sum_alpha_j_zero_weighted = float(M_j[0, 0]) + # = total_weighted_squared - sum_{alpha_j = 0} weighted + # Using cached L[j] and R[j+1]: + # sum_alpha_j_zero = pi * einsum("ij,ia,jb,ab->", L[j], c_j0, c_j0, R[j+1]) + c_j0 = cores[j][:, 0, :] # shape (r_j, r_{j+1}) + sum_alpha_j_zero_weighted = pi * float( + np.einsum("ij,ia,jb,ab->", L[j], c_j0, c_j0, R[j + 1]) + ) total_order_energy[j] = total_weighted_squared - sum_alpha_j_zero_weighted return {