From eb4cf69b773d310dac91b09222f3f3e9d544c4fc Mon Sep 17 00:00:00 2001 From: Hiroshi Shinaoka Date: Fri, 25 Sep 2026 21:42:33 +0900 Subject: [PATCH 1/2] Keep the given order of Matsubara sampling points The C library orders Matsubara sampling points ascending. MatsubaraSampling kept `sampling_points` in the caller's order but returned evaluate() and expected fit() data in the sorted order, so custom points given in any other order silently mismatched the data (plain, positive-only, augmented and DLR bases alike). The wrapper now hands the C library sorted points and permutes the data between the two orders, as TauSampling already follows the given order. Default points are sorted already and take no extra work. Co-Authored-By: Claude Opus 5.5 (1M context) --- src/sparse_ir/sampling.py | 33 ++++++++++++++++++++++++------- tests/test_boundary.py | 41 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 67 insertions(+), 7 deletions(-) diff --git a/src/sparse_ir/sampling.py b/src/sparse_ir/sampling.py index fb513ba..82881db 100644 --- a/src/sparse_ir/sampling.py +++ b/src/sparse_ir/sampling.py @@ -104,7 +104,8 @@ def __init__(self, basis, sampling_points=None, use_positive_taus=True): If `use_positive_taus=True`, the sampling points are folded to the positive tau domain [0, β) [default]; they lie in (0, β). If `use_positive_taus=False`, the sampling points are - unfolded: they lie in (-β/2, β/2] and are symmetric about 0. + unfolded: they lie in (-β/2, β/2] and come in pairs ±τ, plus β/2 + when their number is odd. """ self.basis = basis @@ -308,7 +309,9 @@ class MatsubaraSampling: ``sampling_points`` are *reduced* Matsubara frequencies n, ν = nπ/β: odd integers for a fermionic and even integers for a bosonic basis. A non-integral or wrong-parity value raises :class:`ValueError`; it is - never truncated or adjusted. + never truncated or adjusted. The points may be given in any order; + :py:meth:`evaluate` and :py:meth:`fit` follow the order of + ``sampling_points``. """ def __init__(self, basis, sampling_points=None, positive_only=False): @@ -348,12 +351,23 @@ def __init__(self, basis, sampling_points=None, positive_only=False): "positive_only=True requires non-negative sampling points, " f"got {points[points < 0][0]!r}") self.sampling_points = points + # The C library orders Matsubara points ascending. Hand it sorted + # points and keep the permutation, so that evaluate() and fit() follow + # the order of ``sampling_points``. + order = np.argsort(points, kind="stable") + if np.array_equal(order, np.arange(points.size)): + self._order = self._inverse_order = None + c_points = points + else: + self._order, self._inverse_order = order, np.argsort(order) + c_points = np.ascontiguousarray(points[order]) + self._c_points = c_points self._backend = get_default_blas_backend() if isinstance(basis, augment.AugmentedBasis): # Create sampling object matrix = _util.as_boundary_complex( - basis.uhat(self.sampling_points).T, + basis.uhat(c_points).T, "Matsubara sampling matrix") status = c_int() @@ -362,8 +376,8 @@ def __init__(self, basis, sampling_points=None, positive_only=False): _statistics_to_c(basis.statistics), # statistics c_int(basis.size), # basis_size c_bool(self.positive_only), # positive_only - c_int(len(self.sampling_points)), # num_points - self.sampling_points.ctypes.data_as(POINTER(c_int64)), # points + c_int(len(c_points)), # num_points + c_points.ctypes.data_as(POINTER(c_int64)), # points matrix.ctypes.data_as(POINTER(c_double_complex)), # matrix byref(status) # status ) @@ -376,7 +390,7 @@ def __init__(self, basis, sampling_points=None, positive_only=False): else: # Create sampling object self._ptr = matsubara_sampling_new(basis._ptr, self.positive_only, - self.sampling_points) + c_points) @property def wn(self): @@ -438,7 +452,10 @@ def evaluate(self, al, axis=0): if status != COMPUTATION_SUCCESS: raise RuntimeError(f"Failed to evaluate sampling: {status}") - return output['real'] + 1j * output['imag'] + result = output['real'] + 1j * output['imag'] + if self._inverse_order is not None: + result = np.take(result, self._inverse_order, axis=axis) + return result def fit(self, ax, axis=0): """ @@ -453,6 +470,8 @@ def fit(self, ax, axis=0): """ ax, axis, ndim = _prepare_input(ax, axis, len(self.sampling_points), "Matsubara frequency values") + if self._order is not None: + ax = np.take(ax, self._order, axis=axis) ax = _util.as_boundary_complex(ax, "Matsubara frequency values") output_dims = list(ax.shape) output_dims[axis] = self.basis.size diff --git a/tests/test_boundary.py b/tests/test_boundary.py index 8a84d15..79c8a1a 100644 --- a/tests/test_boundary.py +++ b/tests/test_boundary.py @@ -191,6 +191,47 @@ def test_tau_sampling_keeps_the_given_order(bases): "fit follows the order of the given points") +def _check_matsubara_order(smpl, basis, points): + """evaluate and fit follow the order of the given points, on any axis.""" + np.testing.assert_array_equal(smpl.sampling_points, points) + gl = np.random.default_rng(7).normal(size=basis.size) + ref = gl @ basis.uhat(points) + atol_fit = 100 * smpl.cond * np.finfo(np.float64).eps * np.abs(gl).max() + assert_close(smpl.evaluate(gl), ref, 1e-13 * np.abs(ref).max(), + "evaluate follows the order of the given points") + assert_close(smpl.fit(ref), gl, atol_fit, + "fit follows the order of the given points") + gl2 = np.stack([gl, -2 * gl]) + ref2 = np.stack([ref, -2 * ref]) + assert_close(smpl.evaluate(gl2, axis=1), ref2, 2e-13 * np.abs(ref).max(), + "evaluate along axis 1") + assert_close(smpl.fit(ref2, axis=1), gl2, 2 * atol_fit, "fit along axis 1") + + +@pytest.mark.parametrize("positive_only", [False, True]) +@pytest.mark.parametrize("stat", ["F", "B"]) +def test_matsubara_sampling_keeps_the_given_order(bases, stat, positive_only): + """The C library sorts Matsubara points internally; the results must + still follow the order of the points the caller gave.""" + basis = bases[stat] + points = basis.default_matsubara_sampling_points(positive_only=positive_only) + points = np.random.default_rng(11).permutation(points) + smpl = sparse_ir.MatsubaraSampling(basis, points, positive_only=positive_only) + _check_matsubara_order(smpl, basis, points) + + +def test_matsubara_sampling_of_augmented_basis_and_dlr_keeps_the_given_order(bases): + from sparse_ir import augment + for basis in (augment.AugmentedBasis(bases["B"], augment.MatsubaraConst), + DLR(bases["F"])): + points = bases[basis.statistics].default_matsubara_sampling_points() + if isinstance(basis, augment.AugmentedBasis): + points = basis.default_matsubara_sampling_points() + points = np.random.default_rng(11).permutation(points) + smpl = sparse_ir.MatsubaraSampling(basis, points) + _check_matsubara_order(smpl, basis, points) + + def test_tau_sampling_does_not_alias_the_callers_array(bases): basis = bases["F"] points = np.linspace(0.1, 9.9, basis.size) From 95b553e14bc877c4cf4ec6c27cc11b172361414c Mon Sep 17 00:00:00 2001 From: Hiroshi Shinaoka Date: Fri, 25 Sep 2026 21:42:33 +0900 Subject: [PATCH 2/2] =?UTF-8?q?docs:=20unfolded=20=CF=84=20points=20come?= =?UTF-8?q?=20in=20pairs=20=C2=B1=CF=84,=20plus=20=CE=B2/2=20for=20an=20od?= =?UTF-8?q?d=20number?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The docstrings said the unfolded default τ points are symmetric about 0. For an odd number of points one of them is β/2, which has no partner, so reversing the unfolded points maps τ to -τ only for an even number; reversing the folded default points maps τ to β - τ in every case. Co-Authored-By: Claude Opus 5.5 (1M context) --- src/sparse_ir/augment.py | 2 +- src/sparse_ir/basis.py | 6 ++++-- src/sparse_ir/basis_set.py | 3 ++- 3 files changed, 7 insertions(+), 4 deletions(-) diff --git a/src/sparse_ir/augment.py b/src/sparse_ir/augment.py index d52a55d..399179e 100644 --- a/src/sparse_ir/augment.py +++ b/src/sparse_ir/augment.py @@ -141,7 +141,7 @@ def default_tau_sampling_points(self, *, npoints=None, use_positive_taus=True): If True (default), fold the points into [0, β) with ``np.mod`` and sort them; they lie in (0, β). If False, the points are unfolded: they lie in (-β/2, β/2] - and are symmetric about 0. + and come in pairs ±τ, plus β/2 when their number is odd. .. versionadded:: 1.2 """ diff --git a/src/sparse_ir/basis.py b/src/sparse_ir/basis.py index 9ac2bc1..e1bc9af 100644 --- a/src/sparse_ir/basis.py +++ b/src/sparse_ir/basis.py @@ -343,8 +343,10 @@ def default_tau_sampling_points(self, npoints=None, use_positive_taus=True): If True (default), fold the points into [0, β) with ``np.mod`` and sort them; they lie in (0, β). If False, the points are unfolded: they lie in (-β/2, β/2] - and are symmetric about 0 (for odd L, the point β/2 is its - own mirror image modulo β). + and come in pairs ±τ, plus β/2 when their number is odd. + Reversing the folded points maps τ to β - τ for any number + of points; reversing the unfolded ones maps τ to -τ only for + an even number. .. versionadded:: 1.2 """ diff --git a/src/sparse_ir/basis_set.py b/src/sparse_ir/basis_set.py index 5446a8e..6c82757 100644 --- a/src/sparse_ir/basis_set.py +++ b/src/sparse_ir/basis_set.py @@ -56,7 +56,8 @@ def __init__(self, beta, wmax, eps=None, sve_result=None, use_positive_taus=True If `use_positive_taus=True`, the sampling points are folded to the positive tau domain [0, β) [default]; they lie in (0, β). If `use_positive_taus=False`, the sampling points are - unfolded: they lie in (-β/2, β/2] and are symmetric about 0. + unfolded: they lie in (-β/2, β/2] and come in pairs ±τ, plus β/2 + when their number is odd. """ if sve_result is None: # Create bases by sve of the logistic kernel