Fix audit findings A-H: ctypes boundary, DLR basis functions, and test-suite self-sabotage

pyproject.toml

@@ -1,6 +1,6 @@
 [project]
 name = "sparse-ir"
-version = "2.1.3"
+version = "2.1.4"
 description = "Python bindings for the libsparseir library, providing efficient sparse intermediate representation for many-body physics calculations"
 readme = "README.rst"
 requires-python = ">=3.10"

src/sparse_ir/__init__.py

@@ -30,12 +30,12 @@
     'AbstractBasis', 'FiniteTempBasis', 'finite_temp_bases',
     'TauSampling', 'MatsubaraSampling', 'FiniteTempBasisSet',
     'LogisticKernel', 'RegularizedBoseKernel',
-    'SVEResult', 'compute',
+    'SVEResult', 'compute', 'compute_sve',
 
     # Augmented functionality
     'AugmentedBasis', 'AugmentedTauFunction', 'AugmentedMatsubaraFunction',
     'AbstractAugmentation', 'TauConst', 'TauLinear', 'MatsubaraConst',
 
     # DLR functionality
-    'DiscreteLehmannRepresentation', 'TauPoles', 'MatsubaraPoles',
+    'DiscreteLehmannRepresentation',
 ]

src/sparse_ir/_util.py

@@ -37,6 +37,12 @@ def __call__(self, x):
             return res.reshape(x.shape + res.shape[1:])
 
 
+# Element-type kinds that can be widened to float64 without losing
+# information about what the caller meant: bool, signed/unsigned integer,
+# and floating point.  Complex is deliberately excluded.
+_REAL_KINDS = "biuf"
+
+
 def check_reduced_matsubara(n, zeta=None):
     """Checks that ``n`` is a reduced Matsubara frequency.
 
@@ -48,19 +54,157 @@ def check_reduced_matsubara(n, zeta=None):
     Note that this means that instead of a fermionic frequency (``zeta == 1``),
     we expect an odd integer, while for a bosonic frequency (``zeta == 0``),
     we expect an even one.  If ``zeta`` is omitted, any one is fine.
+
+    Raises:
+        TypeError: if ``n`` is complex.
+        ValueError: if ``n`` is not integral (naming the offending value) or
+            has the wrong parity.
     """
     n = np.asarray(n)
+    if n.dtype.kind == 'c':
+        raise TypeError(
+            f"reduced Matsubara frequency must be real, got dtype {n.dtype}")
     if not np.issubdtype(n.dtype, np.integer):
-        nfloat = n
-        n = nfloat.astype(int)
-        if not (n == nfloat).all():
-            raise ValueError("reduced frequency n must be integer")
+        if n.dtype.kind not in _REAL_KINDS:
+            raise TypeError(
+                f"reduced Matsubara frequency must be numeric, "
+                f"got dtype {n.dtype}")
+        nfloat = np.asarray(n, dtype=np.float64)
+        if not np.all(np.isfinite(nfloat)):
+            raise ValueError(
+                "reduced Matsubara frequency must be finite, got "
+                f"{nfloat[~np.isfinite(nfloat)][0]!r}")
+        n = np.rint(nfloat).astype(np.int64)
+        bad = n != nfloat
+        if bad.any():
+            offending = np.atleast_1d(nfloat)[np.atleast_1d(bad)][0]
+            raise ValueError(
+                "reduced Matsubara frequency must be an integer, got "
+                f"{offending!r} (no truncation is performed)")
     if zeta is not None:
-        if not (n & 1 == zeta).all():
-            raise ValueError("n have wrong parity")
+        parity = np.asarray(n) & 1
+        if not (parity == zeta).all():
+            expected = "odd" if zeta else "even"
+            offending = np.atleast_1d(n)[np.atleast_1d(parity != zeta)][0]
+            raise ValueError(
+                f"reduced Matsubara frequency must be {expected} for "
+                f"zeta={zeta}, got {offending!r}")
     return n
 
 
+def _check_finite(arr, name):
+    if arr.size and not np.all(np.isfinite(arr)):
+        pos = tuple(int(i) for i in np.argwhere(~np.isfinite(arr))[0])
+        raise ValueError(
+            f"{name} must be finite, but contains {arr[pos]!r} at index "
+            f"{pos[0] if arr.ndim == 1 else pos}")
+    return arr
+
+
+def as_boundary_real(a, name="array", check_finite=True):
+    """Normalize ``a`` into a C-contiguous ``float64`` array for the C boundary.
+
+    The returned object is the one whose pointer must be handed to C: a
+    pointer taken from the *original* array would be a defect if a copy was
+    made here (see ``rules/ffi-boundary.md``, Pointer Provenance).
+
+    Raises:
+        TypeError: if ``a`` is complex or of a non-numeric element type.
+        ValueError: if ``a`` contains a non-finite value.
+    """
+    arr = np.asarray(a)
+    if arr.dtype.kind == 'c':
+        raise TypeError(
+            f"{name} must be real-valued, got dtype {arr.dtype}; "
+            "the C entry point takes a double pointer")
+    if arr.dtype.kind not in _REAL_KINDS:
+        raise TypeError(f"{name} has unsupported dtype {arr.dtype}")
+    out = np.ascontiguousarray(arr, dtype=np.float64)
+    if check_finite:
+        _check_finite(out, name)
+    return out
+
+
+def as_boundary_complex(a, name="array", check_finite=True):
+    """Normalize ``a`` into a C-contiguous ``complex128`` array.
+
+    ``complex64`` is *not* ``complex128``: passing its buffer through a
+    ``c_double_complex`` pointer would read twice as many bytes per element
+    as were allocated, so the conversion here is explicit and the pointer
+    must be taken from the returned object.
+    """
+    arr = np.asarray(a)
+    if arr.dtype.kind not in _REAL_KINDS + "c":
+        raise TypeError(f"{name} has unsupported dtype {arr.dtype}")
+    out = np.ascontiguousarray(arr, dtype=np.complex128)
+    if check_finite:
+        _check_finite(out, name)
+    return out
+
+
+def as_boundary_matsubara(n, name="Matsubara indices", zeta=None):
+    """Normalize reduced Matsubara indices into a C-contiguous ``int64`` array.
+
+    Validates integrality (and, if ``zeta`` is given, parity) *before* the
+    conversion, so a non-integral index raises instead of being truncated.
+    """
+    checked = check_reduced_matsubara(n, zeta=zeta)
+    return np.ascontiguousarray(checked, dtype=np.int64)
+
+
+def normalize_axis(axis, ndim):
+    """Resolve a possibly negative ``axis`` against ``ndim`` and range-check it.
+
+    The C API takes a non-negative target dimension; a negative Python axis
+    must be resolved here rather than handed through.
+    """
+    axis = int(axis)
+    resolved = axis + ndim if axis < 0 else axis
+    if not 0 <= resolved < ndim:
+        raise IndexError(
+            f"axis {axis} is out of bounds for an array of dimension {ndim} "
+            f"(valid: {-ndim} .. {ndim - 1})")
+    return resolved
+
+
+def resolve_function_indices(index, size):
+    """Resolve a basis-function index, list of indices, or slice.
+
+    Negative indices are resolved explicitly (Python semantics); an index
+    outside ``[-size, size)`` raises :class:`IndexError` naming the requested
+    index and the valid range.  No modulo wrap-around is performed.
+    """
+    if isinstance(index, slice):
+        return list(range(*index.indices(size)))
+
+    idx = np.asarray(index)
+    if idx.dtype.kind == 'c':
+        raise TypeError(
+            f"basis-function index must be an integer, got dtype {idx.dtype}")
+    if idx.dtype.kind not in _REAL_KINDS:
+        raise TypeError(
+            f"basis-function index must be an integer, got dtype {idx.dtype}")
+    if idx.dtype.kind == 'f':
+        rounded = np.rint(idx)
+        if not np.array_equal(rounded, idx):
+            offending = np.atleast_1d(idx)[np.atleast_1d(rounded != idx)][0]
+            raise ValueError(
+                f"basis-function index must be an integer, got {offending!r} "
+                "(no truncation is performed)")
+        idx = rounded.astype(np.int64)
+
+    flat = np.atleast_1d(idx).ravel()
+    resolved = []
+    for i in flat.tolist():
+        j = i + size if i < 0 else i
+        if not 0 <= j < size:
+            raise IndexError(
+                f"basis-function index {i} is out of range for a function set "
+                f"of size {size} (valid: {-size} .. {size - 1})")
+        resolved.append(int(j))
+    return resolved
+
+
 def check_range(x, xmin, xmax):
     """Checks each element is in range [xmin, xmax]"""
     x = np.asarray(x)

src/sparse_ir/augment.py

@@ -77,7 +77,7 @@ def statistics(self):
         return self._basis.statistics
 
     def __getitem__(self, index):
-        stop = basis._slice_to_size(index)
+        stop = _slice_to_size(index, self.size)
         if stop <= self._naug:
             raise ValueError("Cannot truncate to only augmentation")
         return AugmentedBasis(self._basis[:stop - self._naug],
@@ -202,7 +202,7 @@ def __call__(self, x):
     def __getitem__(self, l):
         # TODO make this more general
         if isinstance(l, slice):
-            stop = basis._slice_to_size(l)
+            stop = _slice_to_size(l, self.size)
             if stop <= self._naug:
                 raise NotImplementedError("Don't truncate to only augmentation")
             return _AugmentedFunction(self._fbasis[:stop-self._naug], self._faug)
@@ -221,13 +221,13 @@ def xmin(self):
 
     @property
     def xmax(self):
-        return self._fbasis.xmin
+        return self._fbasis.xmax
 
     def deriv(self, n=1):
         """Get polynomial for the n'th derivative"""
         dbasis = self._fbasis.deriv(n)
         daug = [faug_l.deriv(n) for faug_l in self._faug]
-        return AugmentedTauFunction(dbasis, *daug)
+        return AugmentedTauFunction(dbasis, daug)
 
 
 class AugmentedMatsubaraFunction(_AugmentedFunction):
@@ -291,6 +291,12 @@ def __init__(self, beta, statistics='B'):
             raise ValueError("temperature must be positive")
         if statistics not in ('F', 'B'):
             raise ValueError("statistics must be 'F' or 'B'")
+        # A fermionic TauConst is not merely ill-conditioned, it is useless:
+        # its Fourier transform is sqrt(beta) * (n == 0), and fermionic reduced
+        # frequencies are odd, so the augmentation column vanishes identically
+        # and the augmented basis is rank-deficient.  Refuse instead of
+        # silently returning a singular fit.
+        _check_bosonic_statistics(statistics)
         self._beta = beta
         self._statistics = statistics
 
@@ -406,6 +412,31 @@ def hat(self, n):
         return np.broadcast_to(1.0, n.shape)
 
 
+def _slice_to_size(index, size):
+    """Return the number of basis functions selected by ``index``.
+
+    Only ``basis[:stop]``-style truncation is supported, mirroring
+    :py:meth:`FiniteTempBasis.__getitem__`.
+    """
+    if not isinstance(index, slice):
+        raise TypeError(
+            f"only slice truncation is supported, got {index!r}")
+    if index.start not in (None, 0):
+        raise ValueError(
+            f"basis truncation must start at 0, got {index.start!r}")
+    if index.step not in (None, 1):
+        raise ValueError(
+            f"basis truncation must have unit step, got {index.step!r}")
+    if index.stop is None:
+        return size
+    stop = int(index.stop)
+    if not 0 < stop <= size:
+        raise IndexError(
+            f"truncation to {stop} functions is out of range for a basis of "
+            f"size {size}")
+    return stop
+
+
 def _augmentation_factory(basis, *augs):
     for aug in augs:
         if isinstance(aug, AbstractAugmentation):
@@ -418,6 +449,10 @@ def _check_bosonic_statistics(statistics):
     if statistics == 'B':
         return
     elif statistics == 'F':
-        raise ValueError("term only allowed for bosonic basis")
+        raise ValueError(
+            "TauConst augmentation is only allowed for a bosonic basis: for "
+            "fermionic statistics its Fourier transform vanishes at every "
+            "(odd) reduced Matsubara frequency, which makes the augmented "
+            "basis rank-deficient")
     else:
-        raise ValueError("invalid statistics")
+        raise ValueError(f"invalid statistics {statistics!r}, expected 'F' or 'B'")

src/sparse_ir/basis.py

@@ -303,16 +303,16 @@ def rescale(self, new_beta):
         temperature.  Note that this implies a different UV cutoff ``wmax``,
         since ``lambda_ == beta * wmax`` stays constant.
         """
-        # Calculate new beta and wmax that give the desired lambda
-        # We keep the ratio beta/wmax constant
-        ratio = self.beta / self.wmax
-        new_wmax = np.sqrt(new_lambda / ratio)
-        new_beta = new_lambda / new_wmax
-
-        # Get epsilon from the current basis accuracy
-        eps = self.accuracy
+        new_beta = float(new_beta)
+        if not new_beta > 0:
+            raise ValueError(
+                f"inverse temperature must be positive, got {new_beta!r}")
 
-        return FiniteTempBasis(self.statistics, new_beta, new_wmax, eps)
+        # lambda_ == beta * wmax is held fixed, so the SVE (which depends only
+        # on lambda_ and eps) can be reused as is.
+        new_wmax = self._lambda / new_beta
+        return FiniteTempBasis(self.statistics, new_beta, new_wmax, self._eps,
+                               kernel=self._kernel, sve_result=self._sve)
 
 
 def finite_temp_bases(beta, wmax, eps=None, sve_result=None):