- expand gm to dirac interface, allow any covariance matrix

- expand tests
- minor fixes
- expand testing input parameters

Author: FEO\preus

src/deterministic_gaussian_sampling/approximation/base_approximation.py

@@ -1,8 +1,27 @@
 import ctypes
+from dataclasses import dataclass
 import numpy
 from typing import Optional
 from deterministic_gaussian_sampling.dll_handling import load_dll
 
+@dataclass
+class CovarianceData:
+    cov: numpy.ndarray
+    eigvals: numpy.ndarray
+    Q: numpy.ndarray
+    sqrt_eigvals: numpy.ndarray
+
+    def __init__(self, cov: numpy.ndarray):
+        self.cov = cov
+        self.eigvals, self.Q = numpy.linalg.eigh(cov)
+        self.eigvals = numpy.maximum(self.eigvals, 1e-14)
+        self.sqrt_eigvals = numpy.sqrt(self.eigvals)
+
+        self.cov = numpy.ascontiguousarray(self.cov, dtype=numpy.float64)
+        self.Q = numpy.ascontiguousarray(self.Q, dtype=numpy.float64)
+        self.eigvals = numpy.ascontiguousarray(self.eigvals, dtype=numpy.float64)
+        self.sqrt_eigvals = numpy.ascontiguousarray(self.sqrt_eigvals, dtype=numpy.float64)
+
 class BaseApproximation:
     cdll: Optional[ctypes.CDLL] = None
 
@@ -22,23 +41,61 @@ def _map_ctypes_numpy(self, ctype) -> type:
             raise TypeError("Unsupported ctype for mapping to numpy dtype")
 
     def _check_numpy_ndarray(
-        self, arr: Optional[numpy.ndarray], L: int, N: int
-    ) -> ctypes.Array:
+        self, arr: Optional[numpy.ndarray], L: int, N: Optional[int] = None
+    ) -> Optional[numpy.ndarray]:
         if arr is None:
             return None
         if not isinstance(arr, numpy.ndarray):
             raise TypeError("Input must be a numpy array")
-        if not numpy.issubdtype(arr.dtype, numpy.floating) and arr.dtype != float:
+        if not numpy.issubdtype(arr.dtype, numpy.floating):
             raise TypeError(
                 f"Input array must be of a floating type, but got {arr.dtype}."
             )
-        if arr.shape != (L, N):
-            row, cols = arr.shape
-            raise ValueError(
-                f"Input array must have size [{L}x{N}] but got [{row}x{cols}]"
-            )
-        return arr.ctypes.data_as(ctypes.POINTER(ctypes.c_double))
+        if L <= 0 or (N is not None and N <= 0):
+            raise ValueError("L and N must be positive integers.")
+        arrShape = arr.shape
+        if arr.ndim == 1:
+            if arrShape[0] != L:
+                raise ValueError(
+                    f"Input array must have size [{L}] but got [{arrShape[0]}]"
+                )
+            if N is not None and N != 1:
+                raise ValueError(
+                    f"Input array must have size [{L}x{N}] but got [{L}]"
+                )
+        elif arr.ndim == 2:
+            if arrShape != (L, N):
+                row, cols = arrShape
+                raise ValueError(
+                    f"Input array must have size [{L}x{N}] but got [{row}x{cols}]"
+                )
+        else:
+            raise ValueError("Input array must be 1D or 2D")
+
+        # ensure contiguous float64 memory
+        return numpy.ascontiguousarray(arr, dtype=numpy.float64)
+    
+    def _check_covariance_matrix(self, cov: numpy.ndarray, N: int, tol = 1e-6) -> CovarianceData:
+        if cov.shape != (N, N):
+            raise ValueError(f"Covariance matrix must be of shape [{N}x{N}] but got {cov.shape}")
+        if not numpy.allclose(cov, cov.T, atol=tol):
+            raise ValueError("Covariance matrix must be symmetric")
+        if numpy.any(numpy.linalg.eigvalsh(cov) <= -tol):
+            raise ValueError("Covariance matrix must be positive definite")
+        
+        covChecked = self._check_numpy_ndarray(cov, N, N)
+        return CovarianceData(covChecked)
+        
+    def _check_weights(self, weight: Optional[numpy.ndarray], size: int, tol=1e-6) -> Optional[numpy.ndarray]:
+        if weight is None:
+            return None
+        wChecked = self._check_numpy_ndarray(weight, size, None)
+        if numpy.any(wChecked < 0):
+            raise ValueError("Weights cannot be negative")
+        if numpy.abs(1 - numpy.sum(wChecked)) > tol:
+            raise ValueError(f"Sum of weights must be 1 within tolerance +/-{tol}, but got {numpy.sum(wChecked)}")
+        return wChecked
 
     @staticmethod
     def _register_cdll(cdll: ctypes.CDLL):
-        BaseApproximation.cdll = cdll
+        BaseApproximation.cdll = cdll

src/deterministic_gaussian_sampling/approximation/dirac_to_dirac.py

@@ -38,17 +38,21 @@ def approximate_double(
         cdll = self.__class__.cdll
         if cdll is None:
             raise OSError("C++-Library was not loaded. Unable to continue!!!")
+        yChecked = self._check_numpy_ndarray(y, M, N)
+        xChecked = self._check_numpy_ndarray(x, L, N)
+        wXChecked = self._check_weights(wX, L)
+        wYChecked = self._check_weights(wY, M)
         minimizer_result = ctypes_wrapper.GslMinimizerResultCTypes()
         success: ctypes.c_bool = cdll.dirac_to_dirac_approx_short_double_approximate(
             self.d2d_short_double,
-            self._check_numpy_ndarray(y, M, N),
+            yChecked.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
             ctypes.c_size_t(M),
             ctypes.c_size_t(L),
             ctypes.c_size_t(N),
             ctypes.c_size_t(100),
-            self._check_numpy_ndarray(x, L, N),
-            self._check_numpy_ndarray(wX, L, 1),
-            self._check_numpy_ndarray(wY, M, 1),
+            xChecked.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
+            None if wXChecked is None else wXChecked.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
+            None if wYChecked is None else wYChecked.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
             ctypes.byref(minimizer_result),
             (
                 None
@@ -69,24 +73,28 @@ def approximate_function_double(
         L: int,
         N: int,
         x: numpy.ndarray,
-        wX: Optional[numpy.ndarray] = None,
-        wY: Optional[numpy.ndarray] = None,
+        wX: python_variant.wXCallbackPythonType,
+        wXD: python_variant.wXDCallbackPythonType,
         options: Optional[python_variant.ApproximateOptionsPy] = None,
     ) -> python_variant.ApproximationResultPy:
         cdll = self.__class__.cdll
         if cdll is None:
             raise OSError("C++-Library was not loaded. Unable to continue!!!")
+        yChecked = self._check_numpy_ndarray(y, M, N)
+        xChecked = self._check_numpy_ndarray(x, L, N)
         minimizer_result = ctypes_wrapper.GslMinimizerResultCTypes()
+        wX = python_variant.wx_callback_python_wrapper(wX)
+        wXD = python_variant.wxd_callback_python_wrapper(wXD)
         success = cdll.dirac_to_dirac_approx_short_function_double_approximate(
             self.d2d_func_double,
-            self._check_numpy_ndarray(y, M, N),
+            yChecked.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
             ctypes.c_size_t(M),
             ctypes.c_size_t(L),
             ctypes.c_size_t(N),
             ctypes.c_size_t(100),
-            self._check_numpy_ndarray(x, L, N),
-            wX,  # TODO: replace with function wrapper
-            wY,  # TODO: replace with function wrapper
+            xChecked.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
+            wX,
+            wXD,
             ctypes.byref(minimizer_result),
             (
                 None
@@ -114,17 +122,21 @@ def approximate_thread_double(
         cdll = self.__class__.cdll
         if cdll is None:
             raise OSError("C++-Library was not loaded. Unable to continue!!!")
+        yChecked = self._check_numpy_ndarray(y, M, N)
+        xChecked = self._check_numpy_ndarray(x, L, N)
+        wXChecked = self._check_weights(wX, L)
+        wYChecked = self._check_weights(wY, M)
         minimizer_result = ctypes_wrapper.GslMinimizerResultCTypes()
         success = cdll.dirac_to_dirac_approx_short_thread_double_approximate(
             self.d2d_thread_double,
-            self._check_numpy_ndarray(y, M, N),
+            yChecked.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
             ctypes.c_size_t(M),
             ctypes.c_size_t(L),
             ctypes.c_size_t(N),
             ctypes.c_size_t(100),
-            self._check_numpy_ndarray(x, L, N),
-            self._check_numpy_ndarray(wX, L, 1),
-            self._check_numpy_ndarray(wY, M, 1),
+            xChecked.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
+            None if wXChecked is None else wXChecked.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
+            None if wYChecked is None else wYChecked.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
             ctypes.byref(minimizer_result),
             (
                 None

src/deterministic_gaussian_sampling/approximation/gaussian_to_dirac.py

@@ -28,7 +28,7 @@ def __del__(self):
 
     def approximate_double(
         self,
-        covDiag: numpy.ndarray,
+        cov: numpy.ndarray,
         L: int,
         N: int,
         x: numpy.ndarray,
@@ -38,25 +38,30 @@ def approximate_double(
         cdll = self.__class__.cdll
         if cdll is None:
             raise OSError("C++-Library was not loaded. Unable to continue!!!")
+        covData = self._check_covariance_matrix(cov, N)
+        xChecked = self._check_numpy_ndarray(x, L, N)
+        wXChecked = self._check_weights(wX, L)
         minimizer_result = ctypes_wrapper.GslMinimizerResultCTypes()
         success = cdll.gm_to_dirac_short_double_approximate(
             self.gm_to_dirac_double,
-            self._check_numpy_ndarray(covDiag, covDiag.shape[0], covDiag.shape[0]),
+            covData.sqrt_eigvals.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
             ctypes.c_size_t(L),
             ctypes.c_size_t(N),
             ctypes.c_size_t(100),
-            self._check_numpy_ndarray(x, L, N),
-            self._check_numpy_ndarray(wX, L, 1),
+            xChecked.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
+            None if wXChecked is None else wXChecked.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
             ctypes.byref(minimizer_result),
             (
                 None
                 if options is None
                 else ctypes_wrapper.ApproximateOptionsCTypes.from_py_type(options)
             ),
         )
-        return python_variant.ApproximationResultPy.from_ctypes(
+        result = python_variant.ApproximationResultPy.from_ctypes(
             success, minimizer_result, x, L, N
         )
+        result.x = result.x @ covData.Q.T
+        return result
 
     def approximate_snd_double(
         self,
@@ -69,14 +74,16 @@ def approximate_snd_double(
         cdll = self.__class__.cdll
         if cdll is None:
             raise OSError("C++-Library was not loaded. Unable to continue!!!")
+        xChecked = self._check_numpy_ndarray(x, L, N)
+        wXChecked = self._check_weights(wX, L)
         minimizer_result = ctypes_wrapper.GslMinimizerResultCTypes()
         success = cdll.gm_to_dirac_short_standard_normal_deviation_double_approximate(
             self.gm_to_dirac_snd_double,
             ctypes.c_size_t(L),
             ctypes.c_size_t(N),
             ctypes.c_size_t(100),
-            self._check_numpy_ndarray(x, L, N),
-            self._check_numpy_ndarray(wX, L, 1),
+            xChecked.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
+            None if wXChecked is None else wXChecked.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
             ctypes.byref(minimizer_result),
             (
                 None

src/deterministic_gaussian_sampling/type_wrapper/__init__.py

@@ -1,4 +1,4 @@
 from .ctypes_wrapper import GslMinimizerResultCTypes, ApproximateOptionsCTypes, wXCallbackCTypes
-from .python_variant import ApproximationResultPy, GslMinimizerResultPy, ApproximateOptionsPy
+from .python_variant import ApproximationResultPy, GslMinimizerResultPy, ApproximateOptionsPy, wx_callback_python_wrapper, wxd_callback_python_wrapper, wXCallbackPythonType, wXDCallbackPythonType
 
-__all__ = ["GslMinimizerResultCTypes", "ApproximateOptionsCTypes", "wXCallbackCTypes", "ApproximationResultPy", "GslMinimizerResultPy", "ApproximateOptionsPy"]
+__all__ = ["GslMinimizerResultCTypes", "ApproximateOptionsCTypes", "wXCallbackCTypes", "ApproximationResultPy", "GslMinimizerResultPy", "ApproximateOptionsPy", "wx_callback_python_wrapper", "wxd_callback_python_wrapper", "wXCallbackPythonType", "wXDCallbackPythonType"]

src/deterministic_gaussian_sampling/type_wrapper/ctypes_wrapper.py

@@ -26,8 +26,7 @@ class GslMinimizerResultCTypes(ctypes.Structure):
 
     @staticmethod
     def from_py_type(pyT: "deterministic_gaussian_sampling.type_wrapper.python_variant.GslMinimizerResultPy") -> "GslMinimizerResultCTypes":
-        # local import to avoid circular import at module import time
-        import deterministic_gaussian_sampling.type_wrapper.python_variant as python_variant  # noqa: F401
+        import deterministic_gaussian_sampling.type_wrapper.python_variant as python_variant
         return GslMinimizerResultCTypes(
             ctypes.c_double(pyT.initalStepSize),
             ctypes.c_double(pyT.stepTolerance),
@@ -47,7 +46,6 @@ def from_py_type(pyT: "deterministic_gaussian_sampling.type_wrapper.python_varia
         )
 
     def to_py_type(self):
-        # local import to avoid circular import at module import time
         import deterministic_gaussian_sampling.type_wrapper.python_variant as python_variant
         return python_variant.GslMinimizerResultPy(
             float(self.initalStepSize),
@@ -82,8 +80,7 @@ class ApproximateOptionsCTypes(ctypes.Structure):
 
     @staticmethod
     def from_py_type(pyT: "deterministic_gaussian_sampling.type_wrapper.python_variant.ApproximateOptionsPy") -> "ApproximateOptionsCTypes":
-        # local import
-        import deterministic_gaussian_sampling.type_wrapper.python_variant as python_variant  # noqa: F401
+        import deterministic_gaussian_sampling.type_wrapper.python_variant as python_variant
         return ApproximateOptionsCTypes(
             ctypes.c_double(pyT.xtolAbs),
             ctypes.c_double(pyT.xtolRel),

src/deterministic_gaussian_sampling/type_wrapper/python_variant.py

@@ -1,9 +1,12 @@
 # python_variant.py
 from __future__ import annotations
+from typing import Callable, Sequence, Union
 from dataclasses import dataclass
 import ctypes
 import numpy as np
 
+from deterministic_gaussian_sampling.type_wrapper import ctypes_wrapper
+
 @dataclass
 class ApproximationResultPy:
     success: bool
@@ -126,3 +129,50 @@ def __str__(self):
         return "\n".join(lines)
 
     __repr__ = __str__
+
+ArrayLike = Union[Sequence[float], np.ndarray]
+# Function type: (x, L, N) -> res
+wXCallbackPythonType = Callable[[ArrayLike, int, int], ArrayLike]
+wXDCallbackPythonType = wXCallbackPythonType
+
+def wx_callback_python_wrapper(func: wXCallbackPythonType) -> "ctypes_wrapper.wXCallbackCTypes":
+    """
+    Wrap a Python function of signature (x: ndarray[L,N], L, N) -> res
+    to a ctypes callback for C.
+    """
+    def c_callback(x_ptr, res_ptr, L, N):
+        L = int(L)
+        N = int(N)
+        size = L * N
+
+        # Convert C pointer to 1D array and reshape to (N,L) for row-major -> column-major
+        x_raw = np.ctypeslib.as_array(x_ptr, shape=(size,))
+        x = x_raw.reshape((L, N))
+
+        # Call user Python function
+        res_val = func(x, L, N)
+        res_val = np.asarray(res_val, dtype=np.float64)
+
+        # Prepare the result to write back to C in row-major
+        if res_val.ndim == 1:
+            # vector case (L,)
+            if res_val.size != L:
+                raise ValueError(f"Returned array has size {res_val.size}, expected {L}")
+            res_arr = np.ctypeslib.as_array(res_ptr, shape=(L,))
+            res_arr[:] = res_val
+        elif res_val.ndim == 2:
+            if res_val.shape != (L, N):
+                raise ValueError(f"Returned array has shape {res_val.shape}, expected ({L},{N})")
+            # Convert to row-major (flattened)
+            res_row_major = res_val.flatten()
+            res_arr = np.ctypeslib.as_array(res_ptr, shape=(size,))
+            res_arr[:] = res_row_major
+        else:
+            raise ValueError(f"Returned array has invalid ndim {res_val.ndim}")
+
+    cb = ctypes_wrapper.wXCallbackCTypes(c_callback)
+    cb._keepalive = func  # keep reference alive
+    return cb
+
+# Alias for your derivative callback if needed
+wxd_callback_python_wrapper = wx_callback_python_wrapper