feat: add use_positive_taus parameter and statistics-dependent periodicity

src/sparse_ir/_util.py

@@ -71,6 +71,72 @@ def check_range(x, xmin, xmax):
     return x
 
 
+def normalize_tau(statistics, tau, beta):
+    """Normalize τ to [0, β] with statistics-dependent periodicity.
+
+    Handles boundary conditions based on statistics:
+    - Fermions ('F'): Anti-periodic G(τ + β) = -G(τ)
+    - Bosons ('B'): Periodic G(τ + β) = G(τ)
+
+    This function maps τ values from the range [-β, β] to [0, β] with
+    appropriate sign factors, following the periodicity rules.
+
+    Arguments:
+        statistics (str):
+            'F' for Fermionic or 'B' for Bosonic statistics.
+        tau (array_like):
+            Imaginary time value(s) in range [-β, β].
+        beta (float):
+            Inverse temperature.
+
+    Returns:
+        tuple[ndarray, ndarray]:
+            (tau_normalized, sign) where:
+            - tau_normalized: τ values mapped to [0, β]
+            - sign: Sign factor (±1) for periodicity
+
+    Raises:
+        ValueError: If tau is outside [-β, β] or statistics is invalid.
+
+    Special cases:
+        - Negative zero (τ = -0.0) is treated as τ = β with appropriate sign
+        - For τ in [0, β]: returns (τ, +1)
+        - For τ in [-β, 0): returns (τ + β, sign) where sign depends on statistics
+
+    .. versionadded:: 1.2
+    """
+    tau = np.asarray(tau, dtype=np.float64)
+    beta = float(beta)
+
+    if statistics not in ('F', 'B'):
+        raise ValueError("statistics must be 'F' (Fermionic) or 'B' (Bosonic)")
+
+    if np.any(tau < -beta) or np.any(tau > beta):
+        raise ValueError(f"τ must be in [-β, β] = [{-beta}, {beta}]")
+
+    # Handle negative zero: τ = -0.0 → τ = β
+    is_neg_zero = (tau == 0.0) & np.signbit(tau)
+
+    tau_normalized = np.where(is_neg_zero, beta, tau)
+    sign = np.ones_like(tau, dtype=np.float64)
+
+    if statistics == 'F':
+        # Fermionic: anti-periodic
+        sign = np.where(is_neg_zero, -1.0, sign)
+    else:  # statistics == 'B'
+        # Bosonic: periodic
+        sign = np.where(is_neg_zero, 1.0, sign)
+
+    # Normalize negative tau to [0, β]
+    mask_neg = tau_normalized < 0
+    tau_normalized = np.where(mask_neg, tau_normalized + beta, tau_normalized)
+
+    if statistics == 'F':
+        sign = np.where(mask_neg, -sign, sign)
+
+    return tau_normalized, sign
+
+
 def check_svd_result(svd_result, matrix_shape=None):
     """Checks that argument is a valid SVD triple (u, s, vH)"""
     u, s, vH = map(np.asarray, svd_result)

src/sparse_ir/augment.py

@@ -111,14 +111,32 @@ def beta(self):
     def wmax(self):
         return self._basis.wmax
 
-    def default_tau_sampling_points(self, *, npoints=None):
+    def default_tau_sampling_points(self, *, npoints=None, use_positive_taus=True):
+        """Get default tau sampling points for augmented basis.
+
+        Arguments:
+            npoints (int):
+                Minimum number of sampling points to return. If None, uses self.size.
+            use_positive_taus (bool):
+                If True, fold points to [0, β] range and sort them (default: True).
+                If False, points are in symmetric range.
+
+                .. versionadded:: 1.2
+        """
         if npoints is None:
             npoints = self.size
 
         # Return the sampling points of the underlying basis, but since we
         # use the size of self, we add two further points.  One then has to
         # hope that these give good sampling points.
-        return basis_get_default_tau_sampling_points_ext(self._basis._ptr, npoints)
+
+        points = basis_get_default_tau_sampling_points_ext(self._basis._ptr, npoints)
+
+        if use_positive_taus:
+            points = np.mod(points, self.beta)
+            points = np.sort(points)
+
+        return points
 
     def default_matsubara_sampling_points(self, *, positive_only=False):
         """Get default Matsubara sampling points for augmented basis.
@@ -247,20 +265,38 @@ def hat(self, n):
 
 
 class TauConst(AbstractAugmentation):
-    """Constant in imaginary time/discrete delta in frequency"""
+    """Constant in imaginary time with statistics-dependent periodicity.
+
+    Evaluates to a constant value in imaginary time with proper handling of
+    periodicity based on statistics:
+    - Fermions: Anti-periodic G(τ + β) = -G(τ)
+    - Bosons: Periodic G(τ + β) = G(τ)
+
+    .. versionchanged:: 1.2
+        Added statistics parameter and support for [-β, β] range.
+    """
     @classmethod
     def create(cls, basis):
-        _check_bosonic_statistics(basis.statistics)
-        return cls(basis.beta)
+        return cls(basis.beta, basis.statistics)
 
-    def __init__(self, beta):
+    def __init__(self, beta, statistics='B'):
+        """
+        Arguments:
+            beta (float):
+                Inverse temperature.
+            statistics (str):
+                'F' for Fermionic or 'B' for Bosonic (default: 'B' for backward compatibility).
+        """
         if beta <= 0:
             raise ValueError("temperature must be positive")
+        if statistics not in ('F', 'B'):
+            raise ValueError("statistics must be 'F' or 'B'")
         self._beta = beta
+        self._statistics = statistics
 
     def __call__(self, tau):
-        tau = _util.check_range(tau, -self._beta, self._beta)
-        return np.broadcast_to(1 / np.sqrt(self._beta), tau.shape)
+        tau_normalized, sign = _util.normalize_tau(self._statistics, tau, self._beta)
+        return sign / np.sqrt(self._beta)
 
     def deriv(self, n=1):
         if n == 0:
@@ -269,27 +305,46 @@ def deriv(self, n=1):
             return lambda tau: np.zeros_like(tau)
 
     def hat(self, n):
-        n = _util.check_reduced_matsubara(n, zeta=0)
+        zeta = 1 if self._statistics == 'F' else 0
+        n = _util.check_reduced_matsubara(n, zeta=zeta)
         return np.sqrt(self._beta) * (n == 0).astype(complex)
 
 
 class TauLinear(AbstractAugmentation):
-    """Linear function in imaginary time, antisymmetric around beta/2"""
+    """Linear function in imaginary time with statistics-dependent periodicity.
+
+    Evaluates to a linear function antisymmetric around β/2 with proper handling
+    of periodicity based on statistics:
+    - Fermions: Anti-periodic G(τ + β) = -G(τ)
+    - Bosons: Periodic G(τ + β) = G(τ)
+
+    .. versionchanged:: 1.2
+        Added statistics parameter and support for [-β, β] range.
+    """
     @classmethod
     def create(cls, basis):
-        _check_bosonic_statistics(basis.statistics)
-        return cls(basis.beta)
+        return cls(basis.beta, basis.statistics)
 
-    def __init__(self, beta):
+    def __init__(self, beta, statistics='B'):
+        """
+        Arguments:
+            beta (float):
+                Inverse temperature.
+            statistics (str):
+                'F' for Fermionic or 'B' for Bosonic (default: 'B' for backward compatibility).
+        """
         if beta <= 0:
             raise ValueError("temperature must be positive")
+        if statistics not in ('F', 'B'):
+            raise ValueError("statistics must be 'F' or 'B'")
         self._beta = beta
+        self._statistics = statistics
         self._norm = np.sqrt(3/beta)
 
     def __call__(self, tau):
-        tau = _util.check_range(tau, -self._beta, self._beta)
-        x = 2/self._beta * tau - 1
-        return self._norm * x
+        tau_normalized, sign = _util.normalize_tau(self._statistics, tau, self._beta)
+        x = 2/self._beta * tau_normalized - 1
+        return sign * self._norm * x
 
     def deriv(self, n=1):
         if n == 0:
@@ -301,22 +356,43 @@ def deriv(self, n=1):
             return lambda tau: np.zeros_like(tau)
 
     def hat(self, n):
-        n = _util.check_reduced_matsubara(n, zeta=0)
+        zeta = 1 if self._statistics == 'F' else 0
+        n = _util.check_reduced_matsubara(n, zeta=zeta)
         inv_w = np.pi/self._beta * n
         inv_w = np.reciprocal(inv_w, out=inv_w, where=n.astype(bool))
         return self._norm * 2/1j * inv_w
 
 
 class MatsubaraConst(AbstractAugmentation):
-    """Constant in Matsubara, undefined in imaginary time"""
+    """Constant in Matsubara, undefined in imaginary time.
+
+    This augmentation is constant in Matsubara frequency space and returns NaN
+    in imaginary time. The statistics parameter is accepted for type consistency
+    but does not affect the behavior.
+
+    .. versionchanged:: 1.2
+        Accepts tau in [-β, β] range (previously [0, β]).
+        Added statistics parameter for consistency.
+    """
     @classmethod
     def create(cls, basis):
-        return cls(basis.beta)
+        return cls(basis.beta, basis.statistics)
 
-    def __init__(self, beta):
+    def __init__(self, beta, statistics=None):
+        """
+        Arguments:
+            beta (float):
+                Inverse temperature.
+            statistics (str, optional):
+                'F' for Fermionic or 'B' for Bosonic. Accepted for type consistency
+                but behavior is identical for both.
+        """
         if beta <= 0:
             raise ValueError("temperature must be positive")
+        if statistics is not None and statistics not in ('F', 'B'):
+            raise ValueError("statistics must be 'F' or 'B'")
         self._beta = beta
+        self._statistics = statistics
 
     def __call__(self, tau):
         tau = _util.check_range(tau, -self._beta, self._beta)

src/sparse_ir/basis.py

@@ -176,9 +176,23 @@ def shape(self):
         """Shape of the basis function set"""
         return self.s.shape
 
-    def default_tau_sampling_points(self, npoints=None):
-        """Get default tau sampling points."""
-        return basis_get_default_tau_sampling_points(self._ptr)
+    def default_tau_sampling_points(self, npoints=None, use_positive_taus=True):
+        """Get default tau sampling points.
+
+        Arguments:
+            npoints (int):
+                Minimum number of sampling points to return (currently unused).
+            use_positive_taus (bool):
+                If True, fold points to [0, β] range and sort them (default: True).
+                If False, points are in symmetric range around β/2.
+
+                .. versionadded:: 1.2
+        """
+        points = basis_get_default_tau_sampling_points(self._ptr)
+        if use_positive_taus:
+            points = np.mod(points, self.beta)
+            points = np.sort(points)
+        return points
 
     def default_omega_sampling_points(self, npoints=None):
         """Return default sampling points in imaginary time.

src/sparse_ir/dlr.py

@@ -209,8 +209,8 @@ def to_IR(self, g_dlr: np.ndarray, axis=0) -> np.ndarray:
             raise RuntimeError(f"Failed to convert DLR to IR: {ret}")
         return output
 
-    def default_tau_sampling_points(self):
-        return self._basis.default_tau_sampling_points()
+    def default_tau_sampling_points(self, **kwargs):
+        return self._basis.default_tau_sampling_points(**kwargs)
 
     def default_matsubara_sampling_points(self, **kwargs):
         return self._basis.default_matsubara_sampling_points(**kwargs)

src/sparse_ir/sampling.py

@@ -43,9 +43,9 @@ def __init__(self, basis, sampling_points=None, use_positive_taus=True):
         self.basis = basis
 
         if sampling_points is None:
-            self.sampling_points = basis.default_tau_sampling_points()
-            if use_positive_taus:
-                self.sampling_points = np.mod(self.sampling_points, basis.beta)
+            self.sampling_points = basis.default_tau_sampling_points(
+                use_positive_taus=use_positive_taus
+            )
         else:
             self.sampling_points = np.asarray(sampling_points, dtype=np.float64)