From 996d0ea9ecb3d90970f8c6d12ebaa96f7615234e Mon Sep 17 00:00:00 2001 From: Hiroshi Shinaoka Date: Fri, 4 Sep 2026 14:42:19 +0200 Subject: [PATCH] Fix audit findings A-H at the ctypes boundary and de-sabotage the test suite Source fixes: * DLR `.u`/`.uhat` returned the *IR* basis functions, so `g_dlr @ dlr.u(tau)` silently produced wrong values. The C API accepts the DLR handle in `spir_basis_get_u`/`spir_basis_get_uhat`, so these properties now build the DLR function sets from the DLR handle (A). * Every remaining boundary crossing normalizes its input with `np.ascontiguousarray(x, dtype=)` and takes the pointer from the converted object: sampling evaluate/fit (both classes, real and complex paths), DLR `from_IR`/`to_IR` real paths, DLR poles, and the evaluation points in `poly.py`. A float32 or non-contiguous array was previously reinterpreted through a `c_double` pointer, reading 8 bytes per 4-byte element (B, C). * `FunctionSet.__call__`/`FunctionSetFT.__call__` had squeeze bugs; the returned shape is now `(n_funcs,) + np.shape(x)` with the function axis dropped only for a single function and the point axes only for scalar `x` (D). * `FiniteTempBasis.rescale` raised `NameError` on an undefined `new_lambda`; it now holds `lambda_ == beta * wmax` fixed, reuses the SVE, and validates `new_beta > 0` (E). * `TauConst` rejects fermionic statistics: its Fourier transform is `sqrt(beta) * (n == 0)` and fermionic reduced frequencies are odd, so the augmentation column vanished identically and the basis was rank-deficient (F). * Reduced Matsubara indices are validated for integrality and parity instead of being truncated by `int()`/`astype(np.int64)`, and the `i % n` index wrap-around in both `__getitem__` implementations is replaced by explicit negative-index resolution plus `IndexError` (G). Also fixed in passing: axis normalization (negative axes resolved before they reach C, out-of-range axes raise `IndexError`), non-finite input rejected before it reaches a factorization, null handle checks, status compared against `COMPUTATION_SUCCESS`, the bare `except: pass` in both `release()` methods, a stray `print` in `overlap`, the `AttributeError`s from `basis._slice_to_size`, `AugmentedTauFunction.xmax` returning `xmin`, the `*daug` argument mismatch in `AugmentedTauFunction.deriv`, and `__all__` naming nonexistent `TauPoles`/`MatsubaraPoles` (which made `from sparse_ir import *` raise). Test suite (H): * The 8 skips in `test_sve_advanced.py` were except-skip blocks hiding a test bug: `sve_result_new`/`basis_new` were called with the Python kernel wrapper instead of `kernel._ptr` (and with the wrong arity). Fixed and un-skipped; the suite now has no skips. * Removed the try/except-skip wrappers in `test_advanced_features.py`, `test_sampling_advanced.py` and `tests/conftest.py`. * Deleted the dead triple-quoted test blocks and the shadowed duplicate `test_broadcast_uv` in `test_poly.py` (tracked in #83), and added the missing `assert` on line 22. * New `tests/test_ffi_boundary.py`: dtype matrix (float32/float64/complex64/complex128 plus integer) and non-contiguous input for every boundary function touched, axis coverage, concrete exception types with `match=`, nonzero-norm assertions, DLR function reconstruction against the IR result, shape tests for the squeeze semantics, `rescale` round-trip, and an `__all__` smoke test. Test counts: 101 passed / 8 skipped -> 181 passed / 0 skipped. Version 2.1.3 -> 2.1.4. Co-Authored-By: Claude Fable 5 --- pyproject.toml | 2 +- src/sparse_ir/__init__.py | 4 +- src/sparse_ir/_util.py | 156 ++++++++- src/sparse_ir/augment.py | 47 ++- src/sparse_ir/basis.py | 18 +- src/sparse_ir/dlr.py | 110 +++++-- src/sparse_ir/poly.py | 201 +++++------- src/sparse_ir/sampling.py | 232 ++++++++----- tests/conftest.py | 19 +- tests/test_advanced_features.py | 239 ++++++-------- tests/test_augment.py | 57 ++-- tests/test_ffi_boundary.py | 564 ++++++++++++++++++++++++++++++++ tests/test_poly.py | 118 +------ tests/test_sampling_advanced.py | 20 +- tests/test_sve_advanced.py | 82 +++-- 15 files changed, 1279 insertions(+), 590 deletions(-) create mode 100644 tests/test_ffi_boundary.py diff --git a/pyproject.toml b/pyproject.toml index 373faac..761f26f 100644 --- a/pyproject.toml +++ b/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" diff --git a/src/sparse_ir/__init__.py b/src/sparse_ir/__init__.py index 247aad4..bbb48e3 100644 --- a/src/sparse_ir/__init__.py +++ b/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', ] \ No newline at end of file diff --git a/src/sparse_ir/_util.py b/src/sparse_ir/_util.py index def84dd..78b03c3 100644 --- a/src/sparse_ir/_util.py +++ b/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) diff --git a/src/sparse_ir/augment.py b/src/sparse_ir/augment.py index a8950af..94f174b 100644 --- a/src/sparse_ir/augment.py +++ b/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") \ No newline at end of file + raise ValueError(f"invalid statistics {statistics!r}, expected 'F' or 'B'") \ No newline at end of file diff --git a/src/sparse_ir/basis.py b/src/sparse_ir/basis.py index 4bedf8a..fb3bc4d 100644 --- a/src/sparse_ir/basis.py +++ 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): diff --git a/src/sparse_ir/dlr.py b/src/sparse_ir/dlr.py index 7c7acf1..c5d4d23 100644 --- a/src/sparse_ir/dlr.py +++ b/src/sparse_ir/dlr.py @@ -16,8 +16,17 @@ COMPUTATION_SUCCESS, get_default_blas_backend, c_double_complex, + basis_get_u, + basis_get_uhat, ) from pylibsparseir.constants import SPIR_ORDER_ROW_MAJOR +from . import _util +from .poly import ( + FunctionSet, + FunctionSetFT, + PiecewiseLegendrePolyVector, + PiecewiseLegendrePolyFTVector, +) class DiscreteLehmannRepresentation(AbstractBasis): """Discrete Lehmann representation (DLR), with poles being extrema of IR. @@ -46,18 +55,65 @@ def __init__(self, basis: AbstractBasis, poles=None): status = ctypes.c_int() if poles is None: poles = basis_get_default_omega_sampling_points(basis._ptr) + # Normalize first, then take the pointer from the *normalized* object. + # Taking it from the caller's array instead silently hands C the + # buffer of a non-contiguous or non-float64 array. + poles = _util.as_boundary_real(poles, "poles") + if poles.ndim != 1: + raise ValueError( + f"poles must be one-dimensional, got shape {poles.shape}") + if poles.size == 0: + raise ValueError("poles must not be empty") self._basis = basis - self._poles = np.ascontiguousarray(poles) + self._poles = poles + self._u = None + self._uhat = None self._backend = get_default_blas_backend() - self._ptr = _lib.spir_dlr_new_with_poles(basis._ptr, len(poles), poles.ctypes.data_as(ctypes.POINTER(ctypes.c_double)), status) + self._ptr = _lib.spir_dlr_new_with_poles( + basis._ptr, + poles.size, + poles.ctypes.data_as(ctypes.POINTER(ctypes.c_double)), + status, + ) if status.value != COMPUTATION_SUCCESS: raise RuntimeError(f"Failed to create DLR basis: {status.value}") + if not self._ptr: + raise RuntimeError("Failed to create DLR basis: null handle") @property - def u(self): return self._basis._u + def u(self): + r"""DLR basis functions on the imaginary-time axis. + + These are the *DLR* basis functions, i.e. ``u[i](tau)`` is the + imaginary-time kernel evaluated at the ``i``-th pole, so that:: + + gtau == g_dlr @ dlr.u(tau) + + holds for DLR coefficients ``g_dlr``. They are **not** the basis + functions of the underlying IR basis. + """ + if self._u is None: + beta = self._basis.beta + self._u = PiecewiseLegendrePolyVector( + FunctionSet(basis_get_u(self._ptr)), + -beta, beta, beta, default_overlap_range=(0, beta)) + return self._u @property - def uhat(self): return self._basis._uhat + def uhat(self): + r"""DLR basis functions on the reduced Matsubara frequency axis. + + ``uhat[i](n)`` is the Fourier transform of :py:attr:`u`, so that:: + + giv == g_dlr @ dlr.uhat(n) + + holds for DLR coefficients ``g_dlr``. They are **not** the Matsubara + basis functions of the underlying IR basis. + """ + if self._uhat is None: + self._uhat = PiecewiseLegendrePolyFTVector( + FunctionSetFT(basis_get_uhat(self._ptr))) + return self._uhat @property def statistics(self): @@ -115,19 +171,25 @@ def from_IR(self, gl: np.ndarray, axis=0) -> np.ndarray: array_like Expansion coefficients in DLR """ - gl = np.ascontiguousarray(gl) + gl = np.asarray(gl) + if gl.ndim == 0: + raise ValueError("IR coefficients must be at least one-dimensional") + axis = _util.normalize_axis(axis, gl.ndim) if gl.shape[axis] != self.basis.size: - raise ValueError(f"Input array has wrong size along dimension {axis}") + raise ValueError( + f"IR coefficients have length {gl.shape[axis]} along axis " + f"{axis}, expected {self.basis.size}") output_dims = list(gl.shape) output_dims[axis] = self.size - ndim = len(gl.shape) - input_dims = np.asarray(gl.shape, dtype=np.int32) + ndim = gl.ndim + input_dims = np.ascontiguousarray(gl.shape, dtype=np.int32) target_dim = axis order = SPIR_ORDER_ROW_MAJOR - if gl.dtype.kind == 'f': + if gl.dtype.kind != 'c': + gl = _util.as_boundary_real(gl, "IR coefficients") output = np.zeros(output_dims, dtype=np.float64) ret = _lib.spir_ir2dlr_dd( self._ptr, @@ -139,8 +201,8 @@ def from_IR(self, gl: np.ndarray, axis=0) -> np.ndarray: gl.ctypes.data_as(ctypes.POINTER(ctypes.c_double)), output.ctypes.data_as(ctypes.POINTER(ctypes.c_double)), ) - elif gl.dtype.kind == 'c': - gl = np.ascontiguousarray(gl, dtype=np.complex128) + else: + gl = _util.as_boundary_complex(gl, "IR coefficients") output_c = np.zeros(output_dims, dtype=c_double_complex) ret = _lib.spir_ir2dlr_zz( self._ptr, @@ -153,8 +215,6 @@ def from_IR(self, gl: np.ndarray, axis=0) -> np.ndarray: output_c.ctypes.data_as(ctypes.POINTER(c_double_complex)), ) output = output_c['real'] + 1j * output_c['imag'] - else: - raise ValueError(f"Unsupported dtype: {gl.dtype}") if ret != COMPUTATION_SUCCESS: raise RuntimeError(f"Failed to convert IR to DLR: {ret}") return output @@ -176,17 +236,23 @@ def to_IR(self, g_dlr: np.ndarray, axis=0) -> np.ndarray: array_like Expansion coefficients in IR """ - g_dlr = np.ascontiguousarray(g_dlr) + g_dlr = np.asarray(g_dlr) + if g_dlr.ndim == 0: + raise ValueError("DLR coefficients must be at least one-dimensional") + axis = _util.normalize_axis(axis, g_dlr.ndim) if g_dlr.shape[axis] != self.size: - raise ValueError(f"Input array has wrong size along dimension {axis}") - output_dims = np.asarray(g_dlr.shape, dtype=np.int32) + raise ValueError( + f"DLR coefficients have length {g_dlr.shape[axis]} along axis " + f"{axis}, expected {self.size}") + output_dims = list(g_dlr.shape) output_dims[axis] = self.basis.size - ndim = len(g_dlr.shape) - input_dims = np.asarray(g_dlr.shape, dtype=np.int32) + ndim = g_dlr.ndim + input_dims = np.ascontiguousarray(g_dlr.shape, dtype=np.int32) target_dim = axis order = SPIR_ORDER_ROW_MAJOR - if g_dlr.dtype.kind == 'f': + if g_dlr.dtype.kind != 'c': + g_dlr = _util.as_boundary_real(g_dlr, "DLR coefficients") output = np.zeros(output_dims, dtype=np.float64) ret = _lib.spir_dlr2ir_dd( self._ptr, @@ -198,8 +264,8 @@ def to_IR(self, g_dlr: np.ndarray, axis=0) -> np.ndarray: g_dlr.ctypes.data_as(ctypes.POINTER(ctypes.c_double)), output.ctypes.data_as(ctypes.POINTER(ctypes.c_double)), ) - elif g_dlr.dtype.kind == 'c': - g_dlr = np.ascontiguousarray(g_dlr, dtype=np.complex128) + else: + g_dlr = _util.as_boundary_complex(g_dlr, "DLR coefficients") output_c = np.zeros(output_dims, dtype=c_double_complex) ret = _lib.spir_dlr2ir_zz( self._ptr, @@ -212,8 +278,6 @@ def to_IR(self, g_dlr: np.ndarray, axis=0) -> np.ndarray: output_c.ctypes.data_as(ctypes.POINTER(c_double_complex)), ) output = output_c['real'] + 1j * output_c['imag'] - else: - raise ValueError(f"Unsupported dtype: {g_dlr.dtype}") if ret != COMPUTATION_SUCCESS: raise RuntimeError(f"Failed to convert DLR to IR: {ret}") return output diff --git a/src/sparse_ir/poly.py b/src/sparse_ir/poly.py index 7e8899f..ffc8254 100644 --- a/src/sparse_ir/poly.py +++ b/src/sparse_ir/poly.py @@ -15,10 +15,11 @@ import weakref import threading -from pylibsparseir.core import _lib, c_double_complex +from pylibsparseir.core import _lib, c_double_complex, COMPUTATION_SUCCESS from pylibsparseir.core import funcs_eval_single_float64, funcs_eval_single_complex128 from pylibsparseir.core import funcs_get_size, funcs_get_knots, SPIR_ORDER_COLUMN_MAJOR from ._gauss import kronrod_31_15 +from . import _util # Global registry to track pointer usage _pointer_registry = weakref.WeakSet() @@ -71,42 +72,41 @@ def __init__(self, funcs_ptr): def size(self): return self._size - """ - Size of returned array is (n_funcs, n_points). - """ def __call__(self, x): - """Evaluate basis functions at given points.""" + """Evaluate basis functions at the given point(s). + + The returned shape is ``(n_funcs,) + np.shape(x)``, with singleton + axes dropped only where they carry no information: + + - a set of ``n_funcs > 1`` functions at a scalar ``x`` gives shape + ``(n_funcs,)``; + - a single function (``n_funcs == 1``) gives exactly ``np.shape(x)``, + i.e. a scalar for scalar ``x``; + - otherwise the full ``(n_funcs,) + np.shape(x)`` is returned, so a + multi-dimensional ``x`` keeps its shape. + """ if self._released: raise RuntimeError("Function set has been released") - x = np.ascontiguousarray(x) + x = np.asarray(x) if x.ndim == 0: - o = funcs_eval_single_float64(self._ptr, x.item()) - if len(o) == 1: - return o[0] - else: - return o + o = np.asarray(funcs_eval_single_float64(self._ptr, float(x))) + return o[0] if self._size == 1 else o o = self.__call_batch(x) - - if x.size == 1 and self._size == 1: - return o.flat[0] - elif x.size == 1 and self._size > 1: - return o.ravel() - elif x.size > 1 and self._size == 1: - return o.ravel() - else: - return o + if self._size == 1: + # Single function: drop the leading function axis, keep x's shape. + return o.reshape(x.shape) + return o def __call_batch(self, x: np.ndarray): - # Use batch evaluation for arrays - x = np.ascontiguousarray(x) - original_shape = x.shape - x_flat = x.ravel() - n_points = len(x_flat) + original_shape = np.shape(x) n_funcs = self._size - # Prepare input array (double) - x_double = x_flat.astype(np.float64) + # Normalize with an explicit dtype and take the pointer from the + # normalized object: a float32/int array reinterpreted through a + # c_double pointer would read 8 bytes per 4-byte element. + x_double = _util.as_boundary_real(np.ravel(x), "evaluation points") + n_points = x_double.size # Prepare output array (double) output = np.zeros((n_funcs, n_points), dtype=np.float64) @@ -120,7 +120,7 @@ def __call_batch(self, x: np.ndarray): output.ctypes.data_as(POINTER(c_double)) ) - if status != 0: + if status != COMPUTATION_SUCCESS: raise RuntimeError(f"Batch evaluation failed with status {status}") # Reshape output to match input shape: (n_funcs, ...) + original_shape @@ -130,26 +130,17 @@ def __call_batch(self, x: np.ndarray): def __getitem__(self, index): - """Get a single basis function or slice of functions.""" + """Get a single basis function or slice of functions. + + Negative indices are resolved with Python semantics; an index outside + ``[-size, size)`` raises :class:`IndexError` rather than being wrapped + around with a modulo. + """ if self._released: raise RuntimeError("Function set has been released") sz = funcs_get_size(self._ptr) - - if isinstance(index, slice): - # Handle slice - start, stop, step = index.indices(sz) - indices = list(range(start, stop, step)) - else: - # Handle single index or list of indices - index = np.asarray(index) - if index.ndim == 0: - # Single index - indices = [int(index) % sz] - else: - # List/array of indices - indices = (index % sz).tolist() - - return funcs_get_slice(self._ptr, indices) + return funcs_get_slice(self._ptr, + _util.resolve_function_indices(index, sz)) def deriv(self, n=1): """Compute the n-th derivative of the basis functions. @@ -173,10 +164,7 @@ def deriv(self, n=1): def release(self): """Manually release the function set.""" if not self._released and self._ptr: - try: - _lib.spir_funcs_release(self._ptr) - except: - pass + _lib.spir_funcs_release(self._ptr) self._released = True self._ptr = None @@ -200,82 +188,65 @@ def size(self): return self._size def __call__(self, x): - """Evaluate basis functions at given points.""" + """Evaluate the basis functions at reduced Matsubara frequencies. + + The returned shape follows the same convention as + :py:meth:`FunctionSet.__call__`: ``(n_funcs,) + np.shape(x)``, with + the leading axis dropped if the set holds a single function and the + trailing axes dropped if ``x`` is a scalar. + + Raises: + ValueError: if any element of ``x`` is not an integer. Reduced + Matsubara frequencies are integers and are never truncated. + """ if self._released: raise RuntimeError("Function set has been released") - x = np.ascontiguousarray(x) - if x.ndim == 0: - o = funcs_eval_single_complex128(self._ptr, x.item()) - if len(o) == 1: - return o[0] - else: - return o - else: - # Use batch evaluation for arrays - original_shape = x.shape - x_flat = x.ravel() - n_points = len(x_flat) - n_funcs = self._size - - # Prepare input array - x_int64 = x_flat.astype(np.int64) - - # Prepare output array (complex128) - output = np.zeros((n_funcs, n_points), dtype=np.complex128) - - # Call batch evaluation function - status = _lib.spir_funcs_batch_eval_matsu( - self._ptr, - SPIR_ORDER_COLUMN_MAJOR, - n_points, - x_int64.ctypes.data_as(POINTER(c_int64)), - output.ctypes.data_as(POINTER(c_double_complex)) # FIX: Matsubara returns complex values - ) - - if status != 0: - raise RuntimeError(f"Batch evaluation failed with status {status}") - - # Reshape output to match input shape: (n_funcs, ...) + original_shape - output = output.reshape((n_funcs,) + original_shape) - - if x.size == 1 and self._size == 1: - return output.flat[0] - elif x.size == 1 and self._size > 1: - return output.ravel() - elif x.size > 1 and self._size == 1: - return output.ravel() - else: - return output + # Validate integrality *before* the int64 conversion: ``astype`` would + # silently turn 1.9 into 1. + x_checked = _util.check_reduced_matsubara(x) + original_shape = x_checked.shape + if x_checked.ndim == 0: + o = np.asarray( + funcs_eval_single_complex128(self._ptr, int(x_checked))) + return o[0] if self._size == 1 else o + + x_int64 = np.ascontiguousarray(np.ravel(x_checked), dtype=np.int64) + n_points = x_int64.size + n_funcs = self._size + output = np.zeros((n_funcs, n_points), dtype=np.complex128) + + status = _lib.spir_funcs_batch_eval_matsu( + self._ptr, + SPIR_ORDER_COLUMN_MAJOR, + n_points, + x_int64.ctypes.data_as(POINTER(c_int64)), + output.ctypes.data_as(POINTER(c_double_complex)) + ) + if status != COMPUTATION_SUCCESS: + raise RuntimeError(f"Batch evaluation failed with status {status}") + + output = output.reshape((n_funcs,) + original_shape) + if n_funcs == 1: + return output.reshape(original_shape) + return output def __getitem__(self, index): - """Get a single basis function or slice of functions.""" + """Get a single basis function or slice of functions. + + Negative indices are resolved with Python semantics; an index outside + ``[-size, size)`` raises :class:`IndexError` rather than being wrapped + around with a modulo. + """ if self._released: raise RuntimeError("Function set has been released") sz = funcs_get_size(self._ptr) - - if isinstance(index, slice): - # Handle slice - start, stop, step = index.indices(sz) - indices = list(range(start, stop, step)) - else: - # Handle single index or list of indices - index = np.asarray(index) - if index.ndim == 0: - # Single index - indices = [int(index) % sz] - else: - # List/array of indices - indices = (index % sz).tolist() - - return funcs_ft_get_slice(self._ptr, indices) + return funcs_ft_get_slice(self._ptr, + _util.resolve_function_indices(index, sz)) def release(self): """Manually release the function set.""" if not self._released and self._ptr: - try: - _lib.spir_funcs_release(self._ptr) - except: - pass + _lib.spir_funcs_release(self._ptr) self._released = True self._ptr = None @@ -365,8 +336,6 @@ def overlap(self, f, xmin: float = None, xmax: float = None, *, rtol=2.3e-16, re int_result = int_result.reshape(int_result.shape[1:]) int_error = int_error.reshape(int_error.shape[1:]) - print(type(int_result), type(int_error)) - if int_result.shape == (): int_result = int_result.item() int_error = int_error.item() diff --git a/src/sparse_ir/sampling.py b/src/sparse_ir/sampling.py index bb1b85a..ec97e9c 100644 --- a/src/sparse_ir/sampling.py +++ b/src/sparse_ir/sampling.py @@ -16,12 +16,47 @@ ) from pylibsparseir.constants import COMPUTATION_SUCCESS, SPIR_ORDER_ROW_MAJOR from . import augment +from . import _util + + +def _zeta(statistics): + """Reduced-frequency offset: 1 for fermions (odd n), 0 for bosons (even n).""" + if statistics == 'F': + return 1 + if statistics == 'B': + return 0 + raise ValueError(f"invalid statistics {statistics!r}, expected 'F' or 'B'") + + +def _prepare_input(a, axis, expected, what): + """Validate an input array's axis and length before it crosses to C. + + Returns ``(array, axis, ndim)`` with ``axis`` resolved to a non-negative + dimension index (the C API takes a non-negative target dimension). + """ + a = np.asarray(a) + if a.ndim == 0: + raise ValueError(f"{what} must be at least one-dimensional") + axis = _util.normalize_axis(axis, a.ndim) + if a.shape[axis] != expected: + raise ValueError( + f"{what} has length {a.shape[axis]} along axis {axis}, " + f"expected {expected}") + return a, axis, a.ndim + class TauSampling: """Sparse sampling in imaginary time. Allows the transformation between the IR basis and a set of sampling points in (scaled/unscaled) imaginary time. + + Note: + Real-valued input (any of ``bool``, integer, ``float32``, ``float64``) + is normalized to ``float64`` and complex input to ``complex128`` before + it crosses the C boundary; narrow types therefore agree with the + ``float64``/``complex128`` result to their own input precision rather + than producing garbage. """ def __init__(self, basis, sampling_points=None, use_positive_taus=True): @@ -43,18 +78,34 @@ 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( + points = basis.default_tau_sampling_points( use_positive_taus=use_positive_taus ) else: - self.sampling_points = np.asarray(sampling_points, dtype=np.float64) + points = sampling_points + points = _util.as_boundary_real(points, "sampling_points") + if points.ndim != 1: + raise ValueError( + f"sampling_points must be one-dimensional, got shape {points.shape}") + if points.size == 0: + raise ValueError("sampling_points must not be empty") + # np.sort returns a fresh C-contiguous array; the pointer below is + # taken from this object, not from the caller's array. + self.sampling_points = np.ascontiguousarray(np.sort(points), + dtype=np.float64) - self.sampling_points = np.sort(self.sampling_points) self._backend = get_default_blas_backend() if isinstance(basis, augment.AugmentedBasis): # Create sampling object # matrix: (n_points, n_funcs) - matrix = np.ascontiguousarray(basis.u(self.sampling_points).T) + matrix = np.asarray(basis.u(self.sampling_points).T) + if matrix.size and not np.all(np.isfinite(matrix)): + raise ValueError( + "tau sampling matrix is not finite: at least one " + "augmentation of this basis is undefined in imaginary " + "time (MatsubaraConst is NaN in tau), so tau sampling " + "cannot be constructed for it") + matrix = _util.as_boundary_real(matrix, "tau sampling matrix") status = c_int() sampling = _lib.spir_tau_sampling_new_with_matrix( SPIR_ORDER_ROW_MAJOR, @@ -67,6 +118,8 @@ def __init__(self, basis, sampling_points=None, use_positive_taus=True): ) if status.value != COMPUTATION_SUCCESS: raise RuntimeError(f"Failed to create tau sampling: {status.value}") + if not sampling: + raise RuntimeError("Failed to create tau sampling: null handle") self._ptr = sampling else: # Create sampling object @@ -91,89 +144,95 @@ def evaluate(self, al, axis=0): Returns: -------- ndarray - Values at sampling points + Values at sampling points. ``float64`` for real input, + ``complex128`` for complex input. """ - al = np.ascontiguousarray(al) + al, axis, ndim = _prepare_input(al, axis, self.basis.size, + "basis coefficients") output_dims = list(al.shape) - ndim = len(output_dims) - input_dims = np.asarray(al.shape, dtype=np.int32) output_dims[axis] = len(self.sampling_points) - if al.dtype.kind == "f": - output = np.zeros(output_dims, dtype=np.float64) + input_dims = np.ascontiguousarray(al.shape, dtype=np.int32) - status = _lib.spir_sampling_eval_dd( + if al.dtype.kind == "c": + al = _util.as_boundary_complex(al, "basis coefficients") + output = np.zeros(output_dims, dtype=c_double_complex) + status = _lib.spir_sampling_eval_zz( self._ptr, self._backend, SPIR_ORDER_ROW_MAJOR, ndim, input_dims.ctypes.data_as(POINTER(c_int)), axis, - al.ctypes.data_as(POINTER(c_double)), - output.ctypes.data_as(POINTER(c_double)) + al.ctypes.data_as(POINTER(c_double_complex)), + output.ctypes.data_as(POINTER(c_double_complex)) ) - elif al.dtype.kind == "c": - output = np.zeros(output_dims, dtype=c_double_complex) - - status = _lib.spir_sampling_eval_zz( + result = output['real'] + 1j * output['imag'] + else: + al = _util.as_boundary_real(al, "basis coefficients") + output = np.zeros(output_dims, dtype=np.float64) + status = _lib.spir_sampling_eval_dd( self._ptr, self._backend, SPIR_ORDER_ROW_MAJOR, ndim, input_dims.ctypes.data_as(POINTER(c_int)), axis, - al.ctypes.data_as(POINTER(c_double_complex)), - output.ctypes.data_as(POINTER(c_double_complex)) + al.ctypes.data_as(POINTER(c_double)), + output.ctypes.data_as(POINTER(c_double)) ) - output = output['real'] + 1j * output['imag'] - else: - raise ValueError(f"Unsupported dtype: {al.dtype}") + result = output if status != COMPUTATION_SUCCESS: raise RuntimeError(f"Failed to evaluate sampling: {status}") - return output + return result def fit(self, ax, axis=0): """ Fit basis coefficients from sampling point values. + + Returns ``float64`` for real input and ``complex128`` for complex + input. """ - ax = np.ascontiguousarray(ax) - ndim = len(ax.shape) - input_dims = np.asarray(ax.shape, dtype=np.int32) + ax, axis, ndim = _prepare_input(ax, axis, len(self.sampling_points), + "sampling point values") output_dims = list(ax.shape) output_dims[axis] = self.basis.size - if ax.dtype.kind == "f": - output = np.zeros(output_dims, dtype=np.float64) - status = _lib.spir_sampling_fit_dd( + input_dims = np.ascontiguousarray(ax.shape, dtype=np.int32) + + if ax.dtype.kind == "c": + ax = _util.as_boundary_complex(ax, "sampling point values") + output = np.zeros(output_dims, dtype=c_double_complex) + status = _lib.spir_sampling_fit_zz( self._ptr, self._backend, SPIR_ORDER_ROW_MAJOR, ndim, input_dims.ctypes.data_as(POINTER(c_int)), axis, - ax.ctypes.data_as(POINTER(c_double)), - output.ctypes.data_as(POINTER(c_double)) + ax.ctypes.data_as(POINTER(c_double_complex)), + output.ctypes.data_as(POINTER(c_double_complex)) ) - elif ax.dtype.kind == "c": - output = np.zeros(output_dims, dtype=c_double_complex) - status = _lib.spir_sampling_fit_zz( + result = output['real'] + 1j * output['imag'] + else: + ax = _util.as_boundary_real(ax, "sampling point values") + output = np.zeros(output_dims, dtype=np.float64) + status = _lib.spir_sampling_fit_dd( self._ptr, self._backend, SPIR_ORDER_ROW_MAJOR, ndim, input_dims.ctypes.data_as(POINTER(c_int)), axis, - ax.ctypes.data_as(POINTER(c_double_complex)), - output.ctypes.data_as(POINTER(c_double_complex)) + ax.ctypes.data_as(POINTER(c_double)), + output.ctypes.data_as(POINTER(c_double)) ) - output = output['real'] + 1j * output['imag'] - else: - raise ValueError(f"Unsupported dtype: {ax.dtype}") + result = output if status != COMPUTATION_SUCCESS: raise RuntimeError(f"Failed to fit sampling: {status}") - return output + return result @property def cond(self): @@ -202,6 +261,12 @@ class MatsubaraSampling: or equivalently, that they are purely real in imaginary time. In this case, sparse sampling is performed over non-negative frequencies only, cutting away half of the necessary sampling space. + + Note: + ``sampling_points`` are *reduced* Matsubara indices: odd integers for + a fermionic and even integers for a bosonic basis. A non-integral or + wrong-parity index raises :class:`ValueError`; it is never truncated + or adjusted. """ def __init__(self, basis, sampling_points=None, positive_only=False): @@ -218,25 +283,36 @@ def __init__(self, basis, sampling_points=None, positive_only=False): If True, use only positive frequencies """ self.basis = basis - self.positive_only = positive_only + self.positive_only = bool(positive_only) + zeta = _zeta(basis.statistics) if sampling_points is None: - self.sampling_points = basis.default_matsubara_sampling_points(positive_only=positive_only) + points = basis.default_matsubara_sampling_points( + positive_only=self.positive_only) else: - self.sampling_points = np.asarray(sampling_points, dtype=np.int64) + points = sampling_points + points = _util.as_boundary_matsubara(points, "sampling_points", + zeta=zeta) + if points.ndim != 1: + raise ValueError( + f"sampling_points must be one-dimensional, got shape {points.shape}") + if points.size == 0: + raise ValueError("sampling_points must not be empty") + self.sampling_points = points self._backend = get_default_blas_backend() if isinstance(basis, augment.AugmentedBasis): # Create sampling object - matrix = basis.uhat(self.sampling_points).T - matrix = np.ascontiguousarray(matrix, dtype=np.complex128) + matrix = _util.as_boundary_complex( + basis.uhat(self.sampling_points).T, + "Matsubara sampling matrix") status = c_int() sampling = _lib.spir_matsu_sampling_new_with_matrix( SPIR_ORDER_ROW_MAJOR, # order _statistics_to_c(basis.statistics), # statistics c_int(basis.size), # basis_size - c_bool(positive_only), # positive_only + 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 matrix.ctypes.data_as(POINTER(c_double_complex)), # matrix @@ -244,10 +320,14 @@ def __init__(self, basis, sampling_points=None, positive_only=False): ) if status.value != COMPUTATION_SUCCESS: raise RuntimeError(f"Failed to create matsubara sampling: {status.value}") + if not sampling: + raise RuntimeError( + "Failed to create matsubara sampling: null handle") self._ptr = sampling else: # Create sampling object - self._ptr = matsubara_sampling_new(basis._ptr, positive_only, self.sampling_points) + self._ptr = matsubara_sampling_new(basis._ptr, self.positive_only, + self.sampling_points) @property def wn(self): @@ -268,65 +348,61 @@ def evaluate(self, al, axis=0): Returns: -------- ndarray - Values at Matsubara frequencies (complex) + Values at Matsubara frequencies (always ``complex128``) """ - # For better numerical stability, we need to make the input array contiguous. - al = np.ascontiguousarray(al) + al, axis, ndim = _prepare_input(al, axis, self.basis.size, + "basis coefficients") output_dims = list(al.shape) - ndim = len(output_dims) - input_dims = np.asarray(al.shape, dtype=np.int32) output_dims[axis] = len(self.sampling_points) - output_cdouble_complex = np.zeros(output_dims, dtype=c_double_complex) - if al.dtype.kind == "f": - status = _lib.spir_sampling_eval_dz( + input_dims = np.ascontiguousarray(al.shape, dtype=np.int32) + output = np.zeros(output_dims, dtype=c_double_complex) + + if al.dtype.kind == "c": + al = _util.as_boundary_complex(al, "basis coefficients") + status = _lib.spir_sampling_eval_zz( self._ptr, self._backend, SPIR_ORDER_ROW_MAJOR, ndim, input_dims.ctypes.data_as(POINTER(c_int)), axis, - al.ctypes.data_as(POINTER(c_double)), - output_cdouble_complex.ctypes.data_as(POINTER(c_double_complex)) + al.ctypes.data_as(POINTER(c_double_complex)), + output.ctypes.data_as(POINTER(c_double_complex)) ) - output = output_cdouble_complex['real'] + 1j * output_cdouble_complex['imag'] - elif al.dtype.kind == "c": - status = _lib.spir_sampling_eval_zz( + else: + al = _util.as_boundary_real(al, "basis coefficients") + status = _lib.spir_sampling_eval_dz( self._ptr, self._backend, SPIR_ORDER_ROW_MAJOR, ndim, input_dims.ctypes.data_as(POINTER(c_int)), axis, - al.ctypes.data_as(POINTER(c_double_complex)), - output_cdouble_complex.ctypes.data_as(POINTER(c_double_complex)) + al.ctypes.data_as(POINTER(c_double)), + output.ctypes.data_as(POINTER(c_double_complex)) ) - output = output_cdouble_complex['real'] + 1j * output_cdouble_complex['imag'] - else: - raise ValueError(f"Unsupported dtype: {al.dtype}") if status != COMPUTATION_SUCCESS: raise RuntimeError(f"Failed to evaluate sampling: {status}") - return output + return output['real'] + 1j * output['imag'] def fit(self, ax, axis=0): """ Fit basis coefficients from Matsubara frequency values. + + Returns ``complex128``; the underlying C entry point + (``spir_sampling_fit_zz``) only exists in the complex flavour, so + real-valued input is widened to ``complex128`` here. Passing the raw + buffer of a real array through a complex pointer would read twice as + many bytes as were allocated. """ - ax = np.asarray(ax) - if ax.dtype.kind not in ("f", "c"): - raise ValueError(f"Unsupported dtype: {ax.dtype}") - # The underlying C entry point (spir_sampling_fit_zz) always expects - # complex128 data. Real-valued input must be normalized to - # complex128 here; otherwise the raw buffer of a float64 array would - # be reinterpreted as complex128 (reading twice as many bytes as - # were allocated), producing an out-of-bounds read and silent - # garbage output instead of a clear error. - ax = np.ascontiguousarray(ax, dtype=np.complex128) - ndim = len(ax.shape) - input_dims = np.asarray(ax.shape, dtype=np.int32) + ax, axis, ndim = _prepare_input(ax, axis, len(self.sampling_points), + "Matsubara frequency values") + ax = _util.as_boundary_complex(ax, "Matsubara frequency values") output_dims = list(ax.shape) output_dims[axis] = self.basis.size + input_dims = np.ascontiguousarray(ax.shape, dtype=np.int32) output = np.zeros(output_dims, dtype=c_double_complex) status = _lib.spir_sampling_fit_zz( diff --git a/tests/conftest.py b/tests/conftest.py index 354de7b..9cf7f61 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -40,9 +40,6 @@ def sve_reg_bose(): @pytest.fixture(scope="session") def test_bases(): """Precomputed test bases for common parameter sets.""" - print("Precomputing test bases ...") - bases = {} - test_params = [ ('F', 1.0, 10.0, 1e-6), # Small fermion ('F', 1.0, 42.0, 1e-8), # Medium fermion @@ -50,14 +47,14 @@ def test_bases(): ('F', 4.0, 20.0, 1e-6), # Different beta ] - for stat, beta, wmax, eps in test_params: - try: - basis = pylibsparseir.FiniteTempBasis(stat, beta, wmax, eps) - bases[(stat, beta, wmax)] = basis - except Exception as e: - print(f"Failed to create basis {(stat, beta, wmax)}: {e}") - - return bases + # A basis that cannot be constructed is a failure, not a missing + # precondition: let the exception propagate instead of silently handing + # tests an incomplete dict. + return { + (stat, beta, wmax): + pylibsparseir.FiniteTempBasis(stat, beta, wmax, eps) + for stat, beta, wmax, eps in test_params + } @pytest.fixture diff --git a/tests/test_advanced_features.py b/tests/test_advanced_features.py index 842dc91..776e86d 100644 --- a/tests/test_advanced_features.py +++ b/tests/test_advanced_features.py @@ -6,6 +6,11 @@ Tests for more sophisticated functionality including noise resilience, multi-dimensional operations, and edge cases. + +Note: these tests deliberately do *not* wrap their bodies in +``try/except: pytest.skip(...)``. A failure here is a defect, not an +unsupported precondition, and swallowing it makes the suite report success +while nothing was checked. """ import pytest @@ -22,42 +27,37 @@ def test_tau_noise_resilience(self, statistics, beta, wmax, rng): """Test that sampling is resilient to noise in tau domain.""" eps = 1e-6 - try: - basis = sparse_ir.FiniteTempBasis(statistics, beta, wmax, eps) - sampling = sparse_ir.TauSampling(basis) - except Exception as e: - pytest.skip(f"Failed to create basis/sampling: {e}") + basis = sparse_ir.FiniteTempBasis(statistics, beta, wmax, eps) + sampling = sparse_ir.TauSampling(basis) # Create synthetic IR coefficients (mimicking sparse-ir test) - try: - # Evaluate v functions at test frequencies - omega_test = np.array([-0.999 * wmax, -0.01 * wmax, 0.5 * wmax]) - v_vals = basis.v(omega_test) + # Evaluate v functions at test frequencies + omega_test = np.array([-0.999 * wmax, -0.01 * wmax, 0.5 * wmax]) + v_vals = basis.v(omega_test) - # Create IR coefficients as linear combination - rhol = v_vals @ np.array([0.8, -0.2, 0.5]) - Gl = basis.s * rhol - Gl_magn = np.linalg.norm(Gl) + # Create IR coefficients as linear combination + rhol = v_vals @ np.array([0.8, -0.2, 0.5]) + Gl = basis.s * rhol + Gl_magn = np.linalg.norm(Gl) + assert Gl_magn > 0 - # Evaluate to tau domain - Gtau = sampling.evaluate(Gl) + # Evaluate to tau domain + Gtau = sampling.evaluate(Gl) + assert np.linalg.norm(Gtau) > 0 - # Add noise - noise_level = 1e-5 - Gtau_noisy = Gtau + noise_level * np.linalg.norm(Gtau) * rng.randn(*Gtau.shape) + # Add noise + noise_level = 1e-5 + Gtau_noisy = Gtau + noise_level * np.linalg.norm(Gtau) * rng.randn(*Gtau.shape) - # Fit back - Gl_recovered = sampling.fit(Gtau_noisy) + # Fit back + Gl_recovered = sampling.fit(Gtau_noisy) - # Check that recovery is reasonable despite noise - recovery_error = np.linalg.norm(Gl - Gl_recovered) - expected_error = 12 * noise_level * Gl_magn # Following sparse-ir pattern + # Check that recovery is reasonable despite noise + recovery_error = np.linalg.norm(Gl - Gl_recovered) + expected_error = 12 * noise_level * Gl_magn # Following sparse-ir pattern - assert recovery_error <= expected_error, \ - f"Recovery error {recovery_error} exceeds expected {expected_error}" - - except Exception as e: - pytest.skip(f"Noise resilience test failed: {e}") + assert recovery_error <= expected_error, \ + f"Recovery error {recovery_error} exceeds expected {expected_error}" class TestEdgeCases: @@ -65,31 +65,24 @@ class TestEdgeCases: def test_small_lambda(self): """Test with very small Lambda values.""" - try: - basis = sparse_ir.FiniteTempBasis('F', 1.0, 1.0, 1e-4) # Lambda = 1 - assert basis.size >= 1 - assert len(basis.s) == basis.size - - # Should still be able to create sampling - sampling = sparse_ir.TauSampling(basis) - assert len(sampling.tau) == basis.size + basis = sparse_ir.FiniteTempBasis('F', 1.0, 1.0, 1e-4) # Lambda = 1 + assert basis.size >= 1 + assert len(basis.s) == basis.size - except Exception as e: - pytest.skip(f"Small lambda test failed: {e}") + # Should still be able to create sampling + sampling = sparse_ir.TauSampling(basis) + assert len(sampling.tau) == basis.size def test_large_lambda(self): """Test with large Lambda values.""" - try: - basis = sparse_ir.FiniteTempBasis('F', 10.0, 100.0, 1e-6) # Lambda = 1000 - assert basis.size > 10 # Should have reasonable number of basis functions - - # Test basic functionality - tau_points = np.linspace(0, basis.beta, 3) - u_vals = basis.u(tau_points) - assert u_vals.shape == (basis.size, 3) + basis = sparse_ir.FiniteTempBasis('F', 10.0, 100.0, 1e-6) # Lambda = 1000 + assert basis.size > 10 # Should have reasonable number of basis functions - except Exception as e: - pytest.skip(f"Large lambda test failed: {e}") + # Test basic functionality + tau_points = np.linspace(0, basis.beta, 3) + u_vals = basis.u(tau_points) + assert u_vals.shape == (basis.size, 3) + assert np.linalg.norm(u_vals) > 0 def test_high_precision(self): """Test with very high precision requirements.""" @@ -102,36 +95,29 @@ def test_high_precision(self): def test_boundary_tau_points(self): """Test evaluation at boundary tau points.""" - try: - basis = sparse_ir.FiniteTempBasis('F', 10.0, 8.0, 1e-6) + basis = sparse_ir.FiniteTempBasis('F', 10.0, 8.0, 1e-6) - # Test at boundaries - tau_boundary = np.array([0.0, basis.beta]) - u_vals = basis.u(tau_boundary) - assert u_vals.shape == (basis.size, 2) - assert np.all(np.isfinite(u_vals)) + # Test at boundaries + tau_boundary = np.array([0.0, basis.beta]) + u_vals = basis.u(tau_boundary) + assert u_vals.shape == (basis.size, 2) + assert np.all(np.isfinite(u_vals)) + assert np.linalg.norm(u_vals) > 0 - # Test very close to boundaries - eps_tau = 1e-10 - tau_near_boundary = np.array([eps_tau, basis.beta - eps_tau]) - u_vals_near = basis.u(tau_near_boundary) - assert np.all(np.isfinite(u_vals_near)) - - except Exception as e: - pytest.skip(f"Boundary tau test failed: {e}") + # Test very close to boundaries + eps_tau = 1e-10 + tau_near_boundary = np.array([eps_tau, basis.beta - eps_tau]) + u_vals_near = basis.u(tau_near_boundary) + assert np.all(np.isfinite(u_vals_near)) def test_zero_frequency(self): """Test evaluation at zero frequency.""" - try: - basis = sparse_ir.FiniteTempBasis('F', 10.0, 8.0, 1e-6) - - # Test v function at omega = 0 - v_zero = basis.v(np.array([0.0])) - assert v_zero.shape == (basis.size, 1) - assert np.all(np.isfinite(v_zero)) + basis = sparse_ir.FiniteTempBasis('F', 10.0, 8.0, 1e-6) - except Exception as e: - pytest.skip(f"Zero frequency test failed: {e}") + # Test v function at omega = 0 + v_zero = basis.v(np.array([0.0])) + assert v_zero.shape == (basis.size, 1) + assert np.all(np.isfinite(v_zero)) class TestConsistencyChecks: @@ -161,45 +147,40 @@ def test_reconstruction_accuracy(self): rel_error = abs(al_test[i] - al_recovered[i]) / max(abs(al_test[i]), 1e-14) assert rel_error < 1e-12, f"Relative error for coefficient {i}: {rel_error}" - def test_singular_value_ordering(self): + @pytest.mark.parametrize("statistics", ['F', 'B']) + @pytest.mark.parametrize("lambda_val", [10, 42]) + def test_singular_value_ordering(self, statistics, lambda_val): """Test that singular values are properly ordered.""" - for statistics in ['F', 'B']: - for lambda_val in [10, 42]: - try: - basis = sparse_ir.FiniteTempBasis(statistics, 1.0, lambda_val, 1e-6) + basis = sparse_ir.FiniteTempBasis(statistics, 1.0, lambda_val, 1e-6) - s = basis.s - # Check decreasing order - assert np.all(s[:-1] >= s[1:]), f"Singular values not decreasing for {statistics}, λ={lambda_val}" + s = basis.s + # Check decreasing order + assert np.all(s[:-1] >= s[1:]), \ + f"Singular values not decreasing for {statistics}, λ={lambda_val}" - # Check positivity - assert np.all(s > 0), f"Non-positive singular values for {statistics}, λ={lambda_val}" + # Check positivity + assert np.all(s > 0), \ + f"Non-positive singular values for {statistics}, λ={lambda_val}" - # Check normalization (first should be largest) - assert s[0] == np.max(s), f"First singular value not largest for {statistics}, λ={lambda_val}" - - except Exception as e: - pytest.skip(f"Singular value test failed for {statistics}, λ={lambda_val}: {e}") + # Check normalization (first should be largest) + assert s[0] == np.max(s), \ + f"First singular value not largest for {statistics}, λ={lambda_val}" def test_significance_properties(self): """Test properties of significance array.""" - try: - basis = sparse_ir.FiniteTempBasis('F', 1.0, 10.0, 1e-6) - - sig = basis.significance + basis = sparse_ir.FiniteTempBasis('F', 1.0, 10.0, 1e-6) - # First significance should be 1 - assert abs(sig[0] - 1.0) < 1e-14, f"First significance {sig[0]} != 1" + sig = basis.significance - # Should be decreasing - assert np.all(sig[:-1] >= sig[1:]), "Significance not decreasing" + # First significance should be 1 + assert abs(sig[0] - 1.0) < 1e-14, f"First significance {sig[0]} != 1" - # Should be between 0 and 1 - assert np.all(sig >= 0), "Negative significance values" - assert np.all(sig <= 1), "Significance values > 1" + # Should be decreasing + assert np.all(sig[:-1] >= sig[1:]), "Significance not decreasing" - except Exception as e: - pytest.skip(f"Significance test failed: {e}") + # Should be between 0 and 1 + assert np.all(sig >= 0), "Negative significance values" + assert np.all(sig <= 1), "Significance values > 1" class TestMemoryAndPerformance: @@ -207,44 +188,36 @@ class TestMemoryAndPerformance: def test_large_basis_creation(self): """Test creating large basis without memory issues.""" - try: - # Create a moderately large basis - basis = sparse_ir.FiniteTempBasis('F', 1.0, 100.0, 1e-8) - - # Should be able to access all properties - size = basis.size - assert size > 20 # Should be reasonably large + # Create a moderately large basis + basis = sparse_ir.FiniteTempBasis('F', 1.0, 100.0, 1e-8) - s = basis.s - assert len(s) == size + # Should be able to access all properties + size = basis.size + assert size > 20 # Should be reasonably large - # Should be able to evaluate functions - tau_test = np.array([0.1, 0.5, 0.9]) - u_vals = basis.u(tau_test) - assert u_vals.shape == (size, 3) + s = basis.s + assert len(s) == size - except MemoryError: - pytest.skip("Not enough memory for large basis test") - except Exception as e: - pytest.skip(f"Large basis test failed: {e}") + # Should be able to evaluate functions + tau_test = np.array([0.1, 0.5, 0.9]) + u_vals = basis.u(tau_test) + assert u_vals.shape == (size, 3) + assert np.linalg.norm(u_vals) > 0 - def test_repeated_operations(self): + def test_repeated_operations(self, rng): """Test that repeated operations don't leak memory or degrade performance.""" - try: - basis = sparse_ir.FiniteTempBasis('F', 1.0, 10.0, 1e-6) - sampling = sparse_ir.TauSampling(basis) - - # Test coefficients - al_test = np.random.random(basis.size) + basis = sparse_ir.FiniteTempBasis('F', 1.0, 10.0, 1e-6) + sampling = sparse_ir.TauSampling(basis) - # Perform many evaluate/fit cycles - for _ in range(10): - ax = sampling.evaluate(al_test) - al_recovered = sampling.fit(ax) + # Test coefficients + al_test = rng.rand(basis.size) - # Should maintain accuracy - error = np.max(np.abs(al_test - al_recovered)) - assert error < 1e-12, f"Accuracy degraded after repeated operations: {error}" + # Perform many evaluate/fit cycles + for _ in range(10): + ax = sampling.evaluate(al_test) + al_recovered = sampling.fit(ax) - except Exception as e: - pytest.skip(f"Repeated operations test failed: {e}") \ No newline at end of file + # Should maintain accuracy + error = np.max(np.abs(al_test - al_recovered)) + assert error < 1e-12, \ + f"Accuracy degraded after repeated operations: {error}" diff --git a/tests/test_augment.py b/tests/test_augment.py index 6ae1802..8cbf0d4 100644 --- a/tests/test_augment.py +++ b/tests/test_augment.py @@ -127,29 +127,27 @@ def test_normalize_tau_errors(): _util.normalize_tau('X', 0.0, beta) -@pytest.mark.parametrize("stat", ["F", "B"]) -def test_tau_const_periodicity(stat): - """Test TauConst with statistics-dependent periodicity""" +def test_tau_const_rejects_fermionic(): + """A fermionic TauConst would be identically zero in Matsubara.""" + with pytest.raises(ValueError, match="only allowed for a bosonic basis"): + augment.TauConst(10.0, 'F') + + basis = sparse_ir.FiniteTempBasis('F', 10.0, wmax=2.0, eps=1e-6) + with pytest.raises(ValueError, match="only allowed for a bosonic basis"): + augment.TauConst.create(basis) + with pytest.raises(ValueError, match="only allowed for a bosonic basis"): + augment.AugmentedBasis(basis, augment.TauConst) + + +def test_tau_const_periodicity(): + """Test TauConst periodicity (bosonic: periodic in beta)""" beta = 10.0 - tc = augment.TauConst(beta, stat) - - # Test at tau=0 - val0 = tc(0.0) - assert np.isclose(val0, 1.0 / np.sqrt(beta)) - - # Test at tau=5 - val_pos = tc(5.0) - assert np.isclose(val_pos, 1.0 / np.sqrt(beta)) - - # Test at tau=-5 (should apply periodicity) - val_neg = tc(-5.0) - - if stat == 'F': - # Fermionic: anti-periodic - assert np.isclose(val_neg, -1.0 / np.sqrt(beta)) - else: - # Bosonic: periodic - assert np.isclose(val_neg, 1.0 / np.sqrt(beta)) + tc = augment.TauConst(beta, 'B') + + assert np.isclose(tc(0.0), 1.0 / np.sqrt(beta)) + assert np.isclose(tc(5.0), 1.0 / np.sqrt(beta)) + # tau=-5 goes through the periodicity branch + assert np.isclose(tc(-5.0), 1.0 / np.sqrt(beta)) @pytest.mark.parametrize("stat", ["F", "B"]) @@ -198,19 +196,18 @@ def test_matsubara_const_range(): mc(beta + 1) -@pytest.mark.parametrize("stat", ["F", "B"]) -def test_tau_const_with_statistics(stat): +def test_tau_const_with_statistics(): """Test TauConst can be created with statistics parameter""" beta = 10.0 - basis = sparse_ir.FiniteTempBasis(stat, beta, wmax=2.0, eps=1e-6) - + basis = sparse_ir.FiniteTempBasis('B', beta, wmax=2.0, eps=1e-6) + # Test factory method tc = augment.TauConst.create(basis) - assert tc._statistics == stat - + assert tc._statistics == 'B' + # Test direct creation - tc2 = augment.TauConst(beta, stat) - assert tc2._statistics == stat + tc2 = augment.TauConst(beta, 'B') + assert tc2._statistics == 'B' # Test evaluation works val = tc(5.0) diff --git a/tests/test_ffi_boundary.py b/tests/test_ffi_boundary.py new file mode 100644 index 0000000..b5eca9e --- /dev/null +++ b/tests/test_ffi_boundary.py @@ -0,0 +1,564 @@ +# Copyright (C) 2020-2025 Satoshi Terasaki, Markus Wallerberger, Hiroshi Shinaoka, and others +# SPDX-License-Identifier: MIT +""" +Regression tests for the ctypes boundary and for the audit findings A-H. + +These cover, for every function that hands a numpy buffer to C: + +* the dtype matrix (``float32``/``float64``/``complex64``/``complex128`` plus + integer input) -- a narrow dtype must be widened, never reinterpreted; +* non-contiguous input -- the pointer must come from the normalized copy; +* concrete exception types with ``match=`` for rejected input; +* nonzero-norm assertions, so a silently zeroed result cannot pass. +""" + +import numpy as np +import pytest + +import sparse_ir +from sparse_ir import _util +from sparse_ir.dlr import DiscreteLehmannRepresentation + + +REAL_DTYPES = [np.float32, np.float64] +COMPLEX_DTYPES = [np.complex64, np.complex128] +ALL_DTYPES = REAL_DTYPES + COMPLEX_DTYPES + +# float32 carries ~7 decimal digits; a widened float32 input can only agree +# with the float64 reference to its own input precision. The coefficient +# vectors below span many orders of magnitude, so the tolerance is applied +# relative to the norm of the reference rather than element-wise. +DTYPE_TOL = { + np.float32: 1e-6, + np.complex64: 1e-6, + np.float64: 1e-14, + np.complex128: 1e-14, +} + + +def assert_close(got, ref, dtype): + """Compare against a float64 reference at the input dtype's precision.""" + scale = np.linalg.norm(ref) + assert scale > 0 + np.testing.assert_allclose(got, ref, rtol=0, atol=DTYPE_TOL[dtype] * scale) + + +@pytest.fixture(scope="module") +def basis(): + return sparse_ir.FiniteTempBasis('F', 2.0, 5.0, 1e-8) + + +@pytest.fixture(scope="module") +def basis_b(): + return sparse_ir.FiniteTempBasis('B', 2.0, 5.0, 1e-8) + + +@pytest.fixture(scope="module") +def gl(basis): + rng = np.random.RandomState(4711) + return basis.s * rng.randn(basis.size) + + +def noncontiguous(a): + """Return a non-contiguous view holding the same values as ``a``.""" + padded = np.zeros(tuple(2 * n for n in a.shape), dtype=a.dtype) + view = padded[(slice(None, None, 2),) * a.ndim] + view[...] = a + assert not view.flags['C_CONTIGUOUS'] + return view + + +# --------------------------------------------------------------------------- +# Finding B: dtype normalization at the boundary (tau sampling) +# --------------------------------------------------------------------------- + +@pytest.mark.parametrize("dtype", ALL_DTYPES) +def test_tau_evaluate_dtype_matrix(basis, gl, dtype): + smpl = sparse_ir.TauSampling(basis) + ref = smpl.evaluate(gl.astype(np.float64)) + assert np.linalg.norm(ref) > 0 + + got = smpl.evaluate(gl.astype(dtype)) + assert got.dtype == (np.complex128 if np.issubdtype(dtype, np.complexfloating) + else np.float64) + assert_close(got, ref, dtype) + + +@pytest.mark.parametrize("dtype", ALL_DTYPES) +def test_tau_fit_dtype_matrix(basis, gl, dtype): + smpl = sparse_ir.TauSampling(basis) + gtau = smpl.evaluate(gl) + ref = smpl.fit(gtau.astype(np.float64)) + assert np.linalg.norm(ref) > 0 + + got = smpl.fit(gtau.astype(dtype)) + assert_close(got, ref, dtype) + assert_close(np.real(got), gl, dtype) + + +@pytest.mark.parametrize("dtype", ALL_DTYPES) +def test_tau_noncontiguous(basis, gl, dtype): + """A strided input must be copied, and the pointer taken from the copy.""" + smpl = sparse_ir.TauSampling(basis) + ref = smpl.evaluate(gl) + got = smpl.evaluate(noncontiguous(gl.astype(dtype))) + assert_close(got, ref, dtype) + + +def test_tau_integer_input(basis): + """Integer coefficients are widened, not reinterpreted as doubles.""" + smpl = sparse_ir.TauSampling(basis) + al = np.zeros(basis.size, dtype=np.int64) + al[0] = 3 + ref = smpl.evaluate(al.astype(np.float64)) + np.testing.assert_allclose(smpl.evaluate(al), ref, rtol=1e-14, atol=0) + assert np.linalg.norm(ref) > 0 + + +@pytest.mark.parametrize("axis", [0, 1, -1, -2]) +def test_tau_axis_coverage(basis, gl, axis): + smpl = sparse_ir.TauSampling(basis) + stacked = np.stack([gl, 2 * gl], axis=1) # (size, 2) + if axis in (1, -1): + stacked = stacked.T # (2, size) + out = smpl.evaluate(stacked, axis=axis) + assert out.shape[axis] == len(smpl.tau) + assert np.linalg.norm(out) > 0 + back = smpl.fit(out, axis=axis) + np.testing.assert_allclose(back, stacked, rtol=0, atol=1e-11) + + +def test_tau_bad_axis(basis, gl): + smpl = sparse_ir.TauSampling(basis) + with pytest.raises(IndexError, match="axis 3 is out of bounds"): + smpl.evaluate(gl, axis=3) + with pytest.raises(IndexError, match="axis -2 is out of bounds"): + smpl.evaluate(gl, axis=-2) + + +def test_tau_bad_length(basis, gl): + smpl = sparse_ir.TauSampling(basis) + with pytest.raises(ValueError, match="expected"): + smpl.evaluate(gl[:-1]) + with pytest.raises(ValueError, match="at least one-dimensional"): + smpl.evaluate(1.0) + + +def test_tau_nonfinite_input(basis, gl): + smpl = sparse_ir.TauSampling(basis) + poisoned = gl.copy() + poisoned[2] = np.nan + with pytest.raises(ValueError, match="must be finite"): + smpl.evaluate(poisoned) + + +def test_tau_sampling_points_validation(basis): + with pytest.raises(ValueError, match="must not be empty"): + sparse_ir.TauSampling(basis, sampling_points=np.array([])) + with pytest.raises(ValueError, match="one-dimensional"): + sparse_ir.TauSampling(basis, sampling_points=np.zeros((2, 2))) + with pytest.raises(TypeError, match="must be real-valued"): + sparse_ir.TauSampling(basis, sampling_points=np.array([1 + 2j])) + + +# --------------------------------------------------------------------------- +# Finding B/G: dtype normalization and index validation (Matsubara sampling) +# --------------------------------------------------------------------------- + +@pytest.mark.parametrize("dtype", ALL_DTYPES) +def test_matsubara_evaluate_dtype_matrix(basis, gl, dtype): + smpl = sparse_ir.MatsubaraSampling(basis) + ref = smpl.evaluate(gl.astype(np.float64)) + assert np.linalg.norm(ref) > 0 + + got = smpl.evaluate(gl.astype(dtype)) + assert got.dtype == np.complex128 + assert_close(got, ref, dtype) + + +@pytest.mark.parametrize("dtype", COMPLEX_DTYPES) +def test_matsubara_fit_dtype_matrix(basis, gl, dtype): + smpl = sparse_ir.MatsubaraSampling(basis) + giv = smpl.evaluate(gl) + ref = smpl.fit(giv) + assert np.linalg.norm(ref) > 0 + + got = smpl.fit(giv.astype(dtype)) + assert got.dtype == np.complex128 + assert_close(got, ref, dtype) + + +@pytest.mark.parametrize("dtype", REAL_DTYPES + [np.int64]) +def test_matsubara_fit_real_input_dtypes(basis, dtype): + """Real-valued input to a complex fit is widened, not reinterpreted.""" + smpl = sparse_ir.MatsubaraSampling(basis) + giv = np.arange(1, len(smpl.wn) + 1).astype(dtype) + ref = smpl.fit(giv.astype(np.complex128)) + assert np.linalg.norm(ref) > 0 + got = smpl.fit(giv) + assert got.dtype == np.complex128 + assert_close(got, ref, np.float64 if dtype is np.int64 else dtype) + + +def test_matsubara_noncontiguous(basis, gl): + smpl = sparse_ir.MatsubaraSampling(basis) + giv = smpl.evaluate(gl) + ref = smpl.fit(giv) + got = smpl.fit(noncontiguous(giv)) + np.testing.assert_allclose(got, ref, rtol=1e-12, atol=0) + + +def test_matsubara_indices_are_not_truncated(basis): + """Finding G: ``int(1.9) == 1`` must not happen silently.""" + with pytest.raises(ValueError, match=r"1\.9"): + sparse_ir.MatsubaraSampling(basis, sampling_points=[1.9, 3.0]) + + +def test_matsubara_index_parity_is_checked(basis, basis_b): + with pytest.raises(ValueError, match="must be odd"): + sparse_ir.MatsubaraSampling(basis, sampling_points=[2, 4]) + with pytest.raises(ValueError, match="must be even"): + sparse_ir.MatsubaraSampling(basis_b, sampling_points=[1, 3]) + + +def test_matsubara_integral_float_indices_are_accepted(basis): + smpl = sparse_ir.MatsubaraSampling(basis, sampling_points=[1.0, -1.0, 3.0]) + assert smpl.wn.dtype == np.int64 + np.testing.assert_array_equal(np.sort(smpl.wn), [-1, 1, 3]) + + +def test_matsubara_complex_indices_rejected(basis): + with pytest.raises(TypeError, match="must be real"): + sparse_ir.MatsubaraSampling(basis, sampling_points=[1 + 0j]) + + +# --------------------------------------------------------------------------- +# Findings A, B, C: DLR +# --------------------------------------------------------------------------- + +def test_dlr_poles_pointer_provenance(basis, gl): + """Finding C: a non-contiguous/float32 pole array must be normalized.""" + poles = basis.default_omega_sampling_points() + ref = DiscreteLehmannRepresentation(basis, poles) + + for candidate in (noncontiguous(poles), poles.astype(np.float32), + list(poles)): + dlr = DiscreteLehmannRepresentation(basis, candidate) + np.testing.assert_allclose(dlr.sampling_points, poles, + rtol=1e-6, atol=1e-6) + np.testing.assert_allclose(dlr.to_IR(dlr.from_IR(gl)), gl, + rtol=0, atol=1e-6 * np.linalg.norm(gl)) + assert ref.size == len(poles) + + +def test_dlr_poles_validation(basis): + with pytest.raises(ValueError, match="must not be empty"): + DiscreteLehmannRepresentation(basis, np.array([])) + with pytest.raises(ValueError, match="one-dimensional"): + DiscreteLehmannRepresentation(basis, np.zeros((2, 2))) + with pytest.raises(ValueError, match="must be finite"): + DiscreteLehmannRepresentation(basis, np.array([1.0, np.nan])) + with pytest.raises(TypeError, match="must be real-valued"): + DiscreteLehmannRepresentation(basis, np.array([1 + 1j])) + + +@pytest.mark.parametrize("dtype", ALL_DTYPES) +def test_dlr_roundtrip_dtype_matrix(basis, gl, dtype): + dlr = DiscreteLehmannRepresentation(basis) + ref = dlr.from_IR(gl.astype(np.float64)) + assert np.linalg.norm(ref) > 0 + + got = dlr.from_IR(gl.astype(dtype)) + assert_close(got, ref, dtype) + + back = dlr.to_IR(got) + np.testing.assert_allclose(np.real(back), gl, + rtol=0, atol=1e-4 * np.linalg.norm(gl)) + + +@pytest.mark.parametrize("dtype", ALL_DTYPES) +def test_dlr_noncontiguous(basis, gl, dtype): + dlr = DiscreteLehmannRepresentation(basis) + ref = dlr.from_IR(gl) + got = dlr.from_IR(noncontiguous(gl.astype(dtype))) + assert_close(got, ref, dtype) + + +@pytest.mark.parametrize("axis", [0, 1, -1]) +def test_dlr_axis_coverage(basis, gl, axis): + dlr = DiscreteLehmannRepresentation(basis) + stacked = np.stack([gl, -gl], axis=0 if axis in (1, -1) else 1) + g_dlr = dlr.from_IR(stacked, axis=axis) + assert g_dlr.shape[axis] == dlr.size + back = dlr.to_IR(g_dlr, axis=axis) + np.testing.assert_allclose(back, stacked, rtol=0, + atol=1e-4 * np.linalg.norm(gl)) + + +def test_dlr_bad_axis_and_length(basis, gl): + dlr = DiscreteLehmannRepresentation(basis) + with pytest.raises(IndexError, match="out of bounds"): + dlr.from_IR(gl, axis=1) + with pytest.raises(ValueError, match="expected"): + dlr.from_IR(gl[:-1]) + with pytest.raises(ValueError, match="at least one-dimensional"): + dlr.to_IR(1.0) + + +def test_dlr_u_are_dlr_functions(basis, gl): + """Finding A: ``dlr.u`` must be the DLR functions, not the IR ones.""" + dlr = DiscreteLehmannRepresentation(basis) + g_dlr = dlr.from_IR(gl) + + assert dlr.u.shape == (dlr.size,) + tau = np.array([0.1, 0.7, 1.3]) + u_vals = dlr.u(tau) + assert u_vals.shape == (dlr.size, tau.size) + assert np.linalg.norm(u_vals) > 0 + + gtau_dlr = g_dlr @ u_vals + gtau_ir = gl @ basis.u(tau) + assert np.linalg.norm(gtau_ir) > 0 + np.testing.assert_allclose(gtau_dlr, gtau_ir, rtol=0, + atol=1e-6 * np.linalg.norm(gtau_ir)) + + +def test_dlr_uhat_are_dlr_functions(basis, gl): + """Finding A, Matsubara side.""" + dlr = DiscreteLehmannRepresentation(basis) + g_dlr = dlr.from_IR(gl) + + n = np.array([-3, 1, 5]) + uhat_vals = dlr.uhat(n) + assert uhat_vals.shape == (dlr.size, n.size) + assert np.linalg.norm(uhat_vals) > 0 + + giv_dlr = g_dlr @ uhat_vals + giv_ir = gl @ basis.uhat(n) + assert np.linalg.norm(giv_ir) > 0 + np.testing.assert_allclose(giv_dlr, giv_ir, rtol=0, + atol=1e-6 * np.linalg.norm(giv_ir)) + + +# --------------------------------------------------------------------------- +# Finding D: shape semantics of FunctionSet.__call__ +# --------------------------------------------------------------------------- + +def test_function_set_shapes(basis): + u = basis.u + size = basis.size + + # Whole set + assert np.shape(u(0.5)) == (size,) + assert u(np.array([0.5, 1.0])).shape == (size, 2) + assert u(np.array([0.5])).shape == (size, 1) + assert u(np.array([[0.1, 0.2], [0.3, 0.4]])).shape == (size, 2, 2) + + # Single function: the function axis disappears, x's shape survives + u0 = u[0] + assert np.isscalar(u0(0.5)) or np.shape(u0(0.5)) == () + assert u0(np.array([0.5])).shape == (1,) + assert u0(np.array([0.5, 1.0])).shape == (2,) + assert u0(np.array([[0.1, 0.2], [0.3, 0.4]])).shape == (2, 2) + + # A slice of several functions keeps both axes + assert u[0:3](np.array([0.5, 1.0])).shape == (3, 2) + assert u[0:3](0.5).shape == (3,) + + # v is defined on omega and follows the same rule + assert basis.v(np.array([0.0])).shape == (size, 1) + + +def test_function_set_ft_shapes(basis): + uhat = basis.uhat + size = basis.size + + assert np.shape(uhat(1)) == (size,) + assert uhat(np.array([1])).shape == (size, 1) + assert uhat(np.array([1, 3])).shape == (size, 2) + assert uhat(np.array([[1, 3], [5, 7]])).shape == (size, 2, 2) + + uhat0 = uhat[0] + assert np.shape(uhat0(1)) == () + assert uhat0(np.array([1])).shape == (1,) + assert uhat0(np.array([1, 3])).shape == (2,) + assert uhat0(np.array([[1, 3], [5, 7]])).shape == (2, 2) + + assert uhat[0:3](np.array([1, 3])).shape == (3, 2) + assert np.linalg.norm(uhat[0:3](np.array([1, 3]))) > 0 + + +def test_function_set_values_are_consistent_across_shapes(basis): + u = basis.u + x = np.array([0.25, 0.5, 1.0, 1.75]) + flat = u(x) + nested = u(x.reshape(2, 2)) + np.testing.assert_array_equal(nested, flat.reshape(basis.size, 2, 2)) + np.testing.assert_array_equal(u[1](x), [u[1](xi) for xi in x]) + + +@pytest.mark.parametrize("dtype", REAL_DTYPES + [np.int64]) +def test_function_set_evaluation_point_dtypes(basis, dtype): + """Finding B: evaluation points of any real dtype are widened.""" + u = basis.u + x = np.array([0.0, 1.0, 2.0]) + ref = u(x) + assert np.linalg.norm(ref) > 0 + np.testing.assert_allclose(u(x.astype(dtype)), ref, rtol=0, + atol=1e-6 * np.linalg.norm(ref)) + np.testing.assert_allclose(u(noncontiguous(x)), ref, rtol=1e-14, atol=0) + + +def test_function_set_rejects_complex_points(basis): + with pytest.raises(TypeError, match="must be real-valued"): + basis.u(np.array([0.5 + 0j])) + + +def test_function_set_ft_rejects_non_integer(basis): + """Finding G: ``astype(np.int64)`` must not truncate 1.9 to 1.""" + with pytest.raises(ValueError, match=r"1\.9"): + basis.uhat(np.array([1.9])) + with pytest.raises(ValueError, match="must be an integer"): + basis.uhat(np.array([1.0, 2.5])) + # integral floats are fine + np.testing.assert_allclose(basis.uhat(np.array([1.0, 3.0])), + basis.uhat(np.array([1, 3])), + rtol=1e-14, atol=0) + + +# --------------------------------------------------------------------------- +# Finding G: index wrap-around +# --------------------------------------------------------------------------- + +def test_function_set_index_out_of_range(basis): + size = basis.size + with pytest.raises(IndexError, match="out of range"): + basis.u[size] + with pytest.raises(IndexError, match="out of range"): + basis.u[-size - 1] + with pytest.raises(IndexError, match="out of range"): + basis.uhat[size + 5] + with pytest.raises(IndexError, match="out of range"): + basis.u[[0, size]] + + +def test_function_set_negative_index_resolves(basis): + last = basis.u[-1] + also_last = basis.u[basis.size - 1] + x = np.array([0.3, 1.1]) + np.testing.assert_array_equal(last(x), also_last(x)) + assert np.linalg.norm(last(x)) > 0 + + +def test_function_set_non_integer_index(basis): + with pytest.raises(ValueError, match="must be an integer"): + basis.u[1.5] + + +def test_resolve_function_indices_slices(): + assert _util.resolve_function_indices(slice(None), 4) == [0, 1, 2, 3] + assert _util.resolve_function_indices(slice(1, 3), 4) == [1, 2] + assert _util.resolve_function_indices(-1, 4) == [3] + assert _util.resolve_function_indices([0, -2], 4) == [0, 2] + + +# --------------------------------------------------------------------------- +# Finding E: rescale +# --------------------------------------------------------------------------- + +def test_rescale_keeps_lambda(basis): + rescaled = basis.rescale(4.0) + assert rescaled.beta == 4.0 + assert np.isclose(rescaled.lambda_, basis.lambda_) + assert np.isclose(rescaled.beta * rescaled.wmax, basis.lambda_) + assert rescaled.statistics == basis.statistics + assert rescaled.size == basis.size + np.testing.assert_allclose(rescaled.s / rescaled.s[0], + basis.s / basis.s[0], rtol=1e-10, atol=0) + + # The rescaled basis is usable + smpl = sparse_ir.TauSampling(rescaled) + al = np.ones(rescaled.size) + assert np.linalg.norm(smpl.evaluate(al)) > 0 + + +def test_rescale_rejects_nonpositive_beta(basis): + with pytest.raises(ValueError, match="must be positive"): + basis.rescale(0.0) + with pytest.raises(ValueError, match="must be positive"): + basis.rescale(-1.0) + + +# --------------------------------------------------------------------------- +# Finding F: fermionic TauConst +# --------------------------------------------------------------------------- + +def test_fermionic_tau_const_rejected(): + import sparse_ir.augment as aug + + fermionic = sparse_ir.FiniteTempBasis('F', 2.0, 5.0, 1e-6) + with pytest.raises(ValueError, match="only allowed for a bosonic basis"): + aug.TauConst(2.0, 'F') + with pytest.raises(ValueError, match="only allowed for a bosonic basis"): + aug.AugmentedBasis(fermionic, aug.TauConst) + + # The bosonic case still works and is full rank in Matsubara + bosonic = sparse_ir.FiniteTempBasis('B', 2.0, 5.0, 1e-6) + augmented = aug.AugmentedBasis(bosonic, aug.TauConst) + assert augmented.size == bosonic.size + 1 + n = bosonic.default_matsubara_sampling_points() + matrix = augmented.uhat(n) + assert np.linalg.norm(matrix[0]) > 0 + + +# --------------------------------------------------------------------------- +# Finding H / python.md: the public API is importable +# --------------------------------------------------------------------------- + +def test_all_names_exist(): + missing = [name for name in sparse_ir.__all__ + if not hasattr(sparse_ir, name)] + assert not missing, f"__all__ names missing from the module: {missing}" + + +def test_star_import_works(): + namespace = {} + exec("from sparse_ir import *", namespace) + for name in sparse_ir.__all__: + assert name in namespace + + +# --------------------------------------------------------------------------- +# Augmented bases: an augmentation undefined in tau cannot be tau-sampled +# --------------------------------------------------------------------------- + +def test_tau_sampling_of_vertex_basis_raises(): + """MatsubaraConst is NaN in imaginary time; sampling it must not silently + hand NaNs to the C factorization.""" + import sparse_ir.augment as aug + + bosonic = sparse_ir.FiniteTempBasis('B', 2.0, 5.0, 1e-6) + vertex = aug.AugmentedBasis(bosonic, aug.MatsubaraConst) + + with pytest.raises(ValueError, match="undefined in imaginary"): + sparse_ir.TauSampling(vertex) + + # The Matsubara side of the same basis is fine + smpl = sparse_ir.MatsubaraSampling(vertex) + al = np.ones(vertex.size) + assert np.linalg.norm(smpl.evaluate(al)) > 0 + + +def test_augmented_tau_sampling_roundtrip(): + """A bosonic TauConst/TauLinear augmented basis still round-trips in tau.""" + import sparse_ir.augment as aug + + bosonic = sparse_ir.FiniteTempBasis('B', 2.0, 5.0, 1e-6) + augmented = aug.AugmentedBasis(bosonic, aug.TauConst, aug.TauLinear) + smpl = sparse_ir.TauSampling(augmented) + + rng = np.random.RandomState(1234) + al = rng.randn(augmented.size) + ax = smpl.evaluate(al) + assert np.linalg.norm(ax) > 0 + np.testing.assert_allclose(smpl.fit(ax), al, rtol=0, atol=1e-8) diff --git a/tests/test_poly.py b/tests/test_poly.py index f97952c..8f63a84 100644 --- a/tests/test_poly.py +++ b/tests/test_poly.py @@ -1,13 +1,9 @@ # Copyright (C) 2020-2025 Satoshi Terasaki, Markus Wallerberger, Hiroshi Shinaoka, and others # SPDX-License-Identifier: MIT -from argparse import ArgumentError -from _pytest.mark import param import numpy as np import sparse_ir -from sparse_ir import poly -from scipy.integrate import quad import pytest @@ -19,7 +15,7 @@ def test_poly(): u1 = basis.u[1] assert np.allclose(u1(np.array([0.5, 0.3, 1.0, 2.0])), np.array([-0.43049722, -0.67225263, -0.18450157, -0.01225698])) - u1(1.0) == -0.18450156753665 + assert np.isclose(u1(1.0), -0.18450156753665) def test_poly_v(): beta = 2 @@ -54,18 +50,6 @@ def test_poly_v(): np.testing.assert_allclose(o, expected, atol=300*eps, rtol=0) -""" -TODO: support part() in C API -def test_shape(sve_logistic): - u, s, v = sve_logistic[42].part() - l = s.size - assert u.shape == (l,) - - assert u[3].shape == () - assert u[2:5].shape == (3,) -""" - - @pytest.mark.parametrize("lambda_, atol", [(1E+4, 5e-13)]) def test_overlap(sve_logistic, lambda_, atol): sve_result = sve_logistic[lambda_] @@ -87,26 +71,6 @@ def test_overlap(sve_logistic, lambda_, atol): -@pytest.mark.parametrize("fn", ["u", "v"]) -def test_broadcast_uv(sve_logistic, fn): - sve_result = sve_logistic[42] - basis = sparse_ir.FiniteTempBasis('F', 4.2, 10, sve_result=sve_result) - - f = getattr(basis, fn) - - u_overlap = basis.u[0:2].overlap(basis.u[0:2], 0.0, beta) - print(u_overlap) - - u_overlap = basis.u.overlap(basis.u, 0.0, beta) - assert u_overlap.shape == (basis.size, basis.size) - np.testing.assert_allclose(u_overlap, np.eye(basis.size), rtol=0.0, atol=atol) - - np.testing.assert_allclose(basis.u[0].overlap(basis.u[1], 0.0, beta), 0, rtol=0.0, atol=atol) - np.testing.assert_allclose(basis.u[0].overlap(basis.u[0], 0.0, beta), 1, rtol=0.0, atol=atol) - np.testing.assert_allclose(basis.u[-1].overlap(basis.u[-1], 0.0, beta), 1, rtol=0.0, atol=atol) - - - @pytest.mark.parametrize("fn", ["u", "v"]) def test_broadcast_uv(sve_logistic, fn): sve_result = sve_logistic[42] @@ -152,83 +116,3 @@ def assert_eq(x, y): np.testing.assert_allclose(x, y, rtol=0, atol=1e-15) # Tensorial assert_eq(f[l](np.reshape(x, (2, 2))), f[l](x).reshape(3, 2, 2)) - - -""" -TODO: support part() in C API -def test_violate(sve_logistic): - u, s, v = sve_logistic[42].part() - - with pytest.raises(ValueError): - u(1.5) - with pytest.raises(ValueError): - v(-3.0) - - -def test_eval(sve_logistic): - u, s, v = sve_logistic[42].part() - l = s.size - - # evaluate - np.testing.assert_array_equal( - u(0.4), [u[i](0.4) for i in range(l)]) - np.testing.assert_array_equal( - u([0.4, -0.2]), - [[u[i](x) for x in (0.4, -0.2)] for i in range(l)]) - - -def test_broadcast(sve_logistic): - u, s, v = sve_logistic[42].part() - - x = [0.3, 0.5] - l = [2, 7] - np.testing.assert_array_equal( - u.value(l, x), [u[ll](xx) for (ll, xx) in zip(l, x)]) - - -def test_matrix_hat(sve_logistic): - u, s, v = sve_logistic[42].part() - uhat = poly.PiecewiseLegendreFT(u, "odd") - - n = np.array([1, 3, 5, -1, -3, 5]) - result = uhat(n.reshape(3, 2)) - result_iter = uhat(n).reshape(-1, 3, 2) - assert result.shape == result_iter.shape - np.testing.assert_array_equal(result, result_iter) - - -@pytest.mark.parametrize("lambda_, atol", [(42, 1e-13), (1E+4, 1e-13)]) -def test_overlap(sve_logistic, lambda_, atol): - u, s, v = sve_logistic[lambda_].part() - - # Keep only even number of polynomials - u, s, v = u[:2*(s.size//2)], s[:2*(s.size//2)], v[:2*(s.size//2)] - - np.testing.assert_allclose(u[0].overlap(u[0]), 1, rtol=0, atol=atol) - - ref = (np.arange(s.size) == 0).astype(float) - np.testing.assert_allclose(u.overlap(u[0]), ref, rtol=0, atol=atol) - - -@pytest.mark.parametrize("lambda_, atol", [(42, 1e-13), (1E+4, 1e-13)]) -def test_overlap_break_points(sve_logistic, lambda_, atol): - u, s, v = sve_logistic[lambda_].part() - - D = 0.5 * v.xmax - rhow = lambda omega: np.where(abs(omega)<=D, 1, 0) - rhol = v.overlap(rhow, points=[-D, D]) - rhol_ref = [quad(v[l], -D, D)[0] for l in range(v.size)] - - np.testing.assert_allclose(rhol, rhol_ref, rtol=0, atol=1e-12*np.abs(rhol_ref).max()) - - -def test_eval_unique(sve_logistic): - u, s, v = sve_logistic[42].part() - uhat = poly.PiecewiseLegendreFT(u, "odd") - - # evaluate - res1 = uhat(np.array([1, 3, 3, 1])) - idx = np.array([0, 1, 1, 0]) - res2 = uhat(np.array([1,3]))[:,idx] - np.testing.assert_array_equal(res1, res2) -""" \ No newline at end of file diff --git a/tests/test_sampling_advanced.py b/tests/test_sampling_advanced.py index 6e085c5..1d955b9 100644 --- a/tests/test_sampling_advanced.py +++ b/tests/test_sampling_advanced.py @@ -100,20 +100,12 @@ def test_matsubara_creation_simple(self): """Test basic Matsubara sampling creation.""" basis = sparse_ir.FiniteTempBasis('F', 1.0, 10.0, 1e-6) - # This might fail due to current implementation issues - # But let's test what we can - try: - # Test with very simple custom points - simple_points = np.array([1, 3], dtype=np.int64) - smpl = sparse_ir.MatsubaraSampling(basis, sampling_points=simple_points) - - assert len(smpl.wn) == 2 - np.testing.assert_array_equal(smpl.wn, simple_points) - - except RuntimeError as e: - # If it fails, at least check the error is consistent - assert "Failed to create Matsubara sampling" in str(e) - pytest.skip(f"Matsubara sampling not yet fully implemented: {e}") + # Test with very simple custom points + simple_points = np.array([1, 3], dtype=np.int64) + smpl = sparse_ir.MatsubaraSampling(basis, sampling_points=simple_points) + + assert len(smpl.wn) == 2 + np.testing.assert_array_equal(smpl.wn, simple_points) class TestSamplingEdgeCases: diff --git a/tests/test_sve_advanced.py b/tests/test_sve_advanced.py index e39f191..6d306aa 100644 --- a/tests/test_sve_advanced.py +++ b/tests/test_sve_advanced.py @@ -24,22 +24,21 @@ def test_sve_accuracy_vs_epsilon(self, lambda_): epsilons = [1e-4, 1e-6, 1e-8, 1e-10] for eps in epsilons: - try: - sve = sve_result_new(kernel, eps) - size = sve_result_get_size(sve) - svals = sve_result_get_svals(sve) - - # Actual accuracy should be <= requested epsilon - actual_accuracy = svals[-1] / svals[0] - assert actual_accuracy <= eps, \ - f"Actual accuracy {actual_accuracy} > requested {eps}" - - # Should have at least 1 singular value - assert size >= 1 - assert len(svals) == size - - except Exception as e: - pytest.skip(f"SVE test failed for λ={lambda_}, ε={eps}: {e}") + # The C entry point takes the kernel *handle*, not the Python + # wrapper object. + sve = sve_result_new(kernel._ptr, eps) + size = sve_result_get_size(sve) + svals = sve_result_get_svals(sve) + + # Actual accuracy should be <= requested epsilon + actual_accuracy = svals[-1] / svals[0] + assert actual_accuracy <= eps, \ + f"Actual accuracy {actual_accuracy} > requested {eps}" + + # Should have at least 1 singular value + assert size >= 1 + assert len(svals) == size + assert np.all(svals > 0) @pytest.mark.parametrize("lambda_", [10, 42]) def test_sve_convergence(self, lambda_): @@ -50,12 +49,8 @@ def test_sve_convergence(self, lambda_): sizes = [] for eps in epsilons: - try: - sve = sve_result_new(kernel, eps) - size = sve_result_get_size(sve) - sizes.append(size) - except Exception as e: - pytest.skip(f"SVE convergence test failed: {e}") + sve = sve_result_new(kernel._ptr, eps) + sizes.append(sve_result_get_size(sve)) # Sizes should be non-decreasing as epsilon decreases for i in range(len(sizes) - 1): @@ -73,32 +68,31 @@ def test_basis_creation_from_sve(self, statistics): beta = 1.0 wmax = lambda_ / beta - try: - # Create kernel and SVE - if statistics == 'F': - kernel = LogisticKernel(lambda_) - else: - kernel = RegularizedBoseKernel(lambda_) + eps = 1e-6 - sve = sve_result_new(kernel, 1e-6) + # Create kernel and SVE + if statistics == 'F': + kernel = LogisticKernel(lambda_) + else: + kernel = RegularizedBoseKernel(lambda_) - # Create basis using C API directly - stats_int = 1 if statistics == 'F' else 0 - basis_c = basis_new(stats_int, beta, wmax, kernel, sve) + sve = sve_result_new(kernel._ptr, eps) - # Test basic properties - size = basis_get_size(basis_c) - assert size > 0 + # Create basis using C API directly. Both the kernel and the SVE are + # passed as handles, and ``basis_new`` takes epsilon and max_size. + stats_int = 1 if statistics == 'F' else 0 + basis_c = basis_new(stats_int, beta, wmax, eps, kernel._ptr, sve, -1) - svals = basis_get_svals(basis_c) - assert len(svals) == size - assert np.all(svals > 0) + # Test basic properties + size = basis_get_size(basis_c) + assert size > 0 - stats_recovered = basis_get_stats(basis_c) - assert stats_recovered == stats_int + svals = basis_get_svals(basis_c) + assert len(svals) == size + assert np.all(svals > 0) - except Exception as e: - pytest.skip(f"Basis from SVE test failed for {statistics}: {e}") + stats_recovered = basis_get_stats(basis_c) + assert stats_recovered == stats_int class TestSVEErrorHandling: @@ -110,8 +104,8 @@ def test_invalid_epsilon_values(self): # Test negative epsilon with pytest.raises(RuntimeError, match="Failed to create"): - sve_result_new(kernel, -1e-6) + sve_result_new(kernel._ptr, -1e-6) # Test zero epsilon with pytest.raises(RuntimeError, match="Failed to create"): - sve_result_new(kernel, 0.0) \ No newline at end of file + sve_result_new(kernel._ptr, 0.0) \ No newline at end of file