Keep the given order of Matsubara sampling points

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
         """

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
         """

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

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

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)