Vector Fitting Tools

CHANGELOG.rst

@@ -15,7 +15,7 @@ and this project adheres to `Semantic Versioning <https://semver.org/spec/v2.0.0
 - Convolve and DyConvolve can also take a 1D signal and will return a transformed 1D signal in that case
 - impulse function now 1D by default
 - Misc examples updated to show this feature 
-- Broke up benchmark into four seperatea performance tests
+- Broke up benchmark into four separate performance tests
 
 **Added**
 
@@ -44,7 +44,7 @@ and this project adheres to `Semantic Versioning <https://semver.org/spec/v2.0.0
 **Changed**
 
 - Expanded test suite to full code and branch coverage, including all edge cases and defensive branches.
-- Tests to do cover KLU wrapper, because it is not used in this version.
+- Tests do not cover KLU wrapper, because it is not used in this version.
 
 [0.3.4] - 2026-01-06
 --------------------

docs/conf.py

@@ -63,7 +63,7 @@
 master_doc = "index"
 
 project = "Sparse SGWT"
-copyright = "2025, Luke Lowery"
+copyright = "2025-2026, Luke Lowery"
 author = "Luke Lowery"
 version = importlib_metadata.version("sgwt")
 release = version

docs/dev/tests.rst

@@ -8,8 +8,12 @@ Test Modules
 
 - **test_cholconv.py**: Static and dynamic graph convolution (``Convolve``, ``DyConvolve``), VF kernels, topology updates.
 - **test_chebyshev.py**: Chebyshev polynomial approximation and ``ChebyConvolve``.
+- **test_chebymodel.py**: ``ChebyModel`` fitting logic and kernel construction.
 - **test_functions.py**: Analytical filter functions (low-pass, band-pass, high-pass).
+- **test_sgma.py**: Spectral Graph Modal Analysis (``sgma``, ``ModeTable``).
 - **test_util.py**: Resource loading, DLL utilities, ``ChebyKernel``, ``VFKernel``.
+- **test_vfit.py**: Vector Fitting rational approximation (``VFModel``).
+- **test_performance.py**: Performance benchmarks (excluded by default).
 
 Running Tests
 -------------

docs/library/json.rst

@@ -32,7 +32,7 @@ Below is a truncated example from the built-in ``MODIFIED_MORLET.json`` file.
 .. code-block:: json
 
     {
-        "description": "Modified Morlet Wavelet With Central Frequnecy of 2pi",
+        "description": "Modified Morlet Wavelet With Central Frequency of 2pi",
         "nfuncs": 1,
         "npoles": 14,
         "d": 0.0027966205394028575,

docs/usage/chebyshev_kernels.rst

@@ -32,17 +32,13 @@ Basic Usage
     def f(x):
         return np.array([my_filter(x)]).T
 
-    # Initialize convolution context to get spectrum bound
-    conv = sgwt.ChebConvolve(L)
-    ubnd = conv.spectrum_bound
-
     # Create the kernel
-    kernel = sgwt.ChebyKernel.from_function(f, order=50, spectrum_bound=ubnd)
+    kernel = sgwt.ChebyModel.kernel(L, f, order=50)
 
     # Apply convolution
-    with conv:
+    with sgwt.ChebyConvolve(L) as conv:
         result = conv.convolve(signal, kernel)
 
 .. seealso::
    :doc:`../theory/theory_kernel_fitting`
-      For a comparison between polynomial and rational approximations.
+      For a comparison between polynomial and rational approximations.

docs/usage/input_signals.rst

@@ -53,4 +53,4 @@ For testing and examples, the library provides the :func:`~sgwt.util.impulse` he
     X_impulse = impulse(L, n=600)
 
     print(X_impulse.shape)
-    # Output: (2000, 1)
+    # Output: (2000,)

examples/demo_cheby_coeffs.py

@@ -32,7 +32,7 @@ def f(x): return np.array([sgwt.functions.bandpass(x, scale=1.0, order=1)]).T
 
 for order in orders:
     # Fit kernel using the module's function (uses quadratic sampling by default)
-    kernel = sgwt.ChebyKernel.from_function(f, order, ubnd, min_lambda=1e-4)
+    kernel = sgwt.ChebyModel.kernel(L, f, order, min_lambda=1e-4)
     y_approx = kernel.evaluate(x_eval)
 
     ax1.plot(x_eval, y_approx, label=f'Order {order}')
@@ -66,4 +66,4 @@ def f(x): return np.array([sgwt.functions.bandpass(x, scale=1.0, order=1)]).T
 os.makedirs(static_images_dir, exist_ok=True)
 save_path = os.path.join(static_images_dir, 'demo_cheby_coeffs.png')
 plt.savefig(save_path, dpi=400, bbox_inches='tight')
-plt.show()
+plt.show()

examples/demo_cheby_convolve_1.py

@@ -25,7 +25,7 @@ def f(x):
         sgwt.functions.bandpass(x, scale=s, order=2) for s in SCALES], 
     axis=1)
 
-kernel = sgwt.ChebyKernel.from_function_on_graph(L, f, ORDER, min_lambda=XMIN)
+kernel = sgwt.ChebyModel.kernel(L, f, ORDER, min_lambda=XMIN)
 
 X = impulse(L, n=600) 
 

examples/demo_cheby_time.py

@@ -31,7 +31,7 @@
 def f(x): return np.stack([sgwt.functions.bandpass(x, scale=s, order=1) for s in SCALES], axis=1)
 
 lbnd = 1e-3
-kernel = sgwt.ChebyKernel.from_function_on_graph(L, f, ORDER, min_lambda=lbnd)
+kernel = sgwt.ChebyModel.kernel(L, f, ORDER, min_lambda=lbnd)
 
 print(f"Benchmarking on {L.shape[0]} nodes...")
 
@@ -114,4 +114,4 @@ def f(x): return np.stack([sgwt.functions.bandpass(x, scale=s, order=1) for s in
 os.makedirs(static_images_dir, exist_ok=True)
 save_path = os.path.join(static_images_dir, 'demo_cheby_time.png')
 plt.savefig(save_path, dpi=400, bbox_inches='tight')
-plt.show()
+plt.show()

pyproject.toml

@@ -75,4 +75,5 @@ omit = [
 exclude_lines = [
     "pragma: no cover",
     "if __name__ == .__main__.:",
+    "if TYPE_CHECKING:",
 ]

sgwt/__init__.py

@@ -16,6 +16,9 @@
 # Spectral Graph Modal Analysis
 from .sgma import SGMA, ModeTable
 
+# Fitting models
+from .fitting import VFModel, ChebyModel
+
 # Analytical function generators
 from . import functions
 
@@ -28,8 +31,66 @@
     estimate_spectral_bound,
     get_cholmod_dll,
     get_klu_dll,
+    load_ply_laplacian,
+    load_ply_xyz,
+    list_graphs,
 )
 
+# Static type declarations for lazy-loaded resources (IDE intellisense only)
+from typing import TYPE_CHECKING
+if TYPE_CHECKING:
+    from typing import Dict, Any
+    import numpy as np
+    from scipy.sparse import csc_matrix
+
+    # Kernels
+    GAUSSIAN_WAV: Dict[str, Any]
+    MEXICAN_HAT: Dict[str, Any]
+    MODIFIED_MORLET: Dict[str, Any]
+    SHANNON: Dict[str, Any]
+
+    # Delay Laplacians
+    DELAY_EAST: csc_matrix
+    DELAY_EASTWEST: csc_matrix
+    DELAY_HAWAII: csc_matrix
+    DELAY_NEISO: csc_matrix
+    DELAY_TEXAS: csc_matrix
+    DELAY_USA: csc_matrix
+    DELAY_WECC: csc_matrix
+
+    # Impedance Laplacians
+    IMPEDANCE_EASTWEST: csc_matrix
+    IMPEDANCE_HAWAII: csc_matrix
+    IMPEDANCE_TEXAS: csc_matrix
+    IMPEDANCE_USA: csc_matrix
+    IMPEDANCE_WECC: csc_matrix
+
+    # Length Laplacians
+    LENGTH_EASTWEST: csc_matrix
+    LENGTH_HAWAII: csc_matrix
+    LENGTH_NEISO: csc_matrix
+    LENGTH_TEXAS: csc_matrix
+    LENGTH_USA: csc_matrix
+    LENGTH_WECC: csc_matrix
+
+    # Coordinate Signals
+    COORD_EASTWEST: np.ndarray
+    COORD_HAWAII: np.ndarray
+    COORD_TEXAS: np.ndarray
+    COORD_USA: np.ndarray
+
+    # Mesh Laplacians
+    MESH_BUNNY: csc_matrix
+    MESH_HORSE: csc_matrix
+    MESH_LBRAIN: csc_matrix
+    MESH_RBRAIN: csc_matrix
+
+    # Mesh Coordinates
+    BUNNY_XYZ: np.ndarray
+    HORSE_XYZ: np.ndarray
+    LBRAIN_XYZ: np.ndarray
+    RBRAIN_XYZ: np.ndarray
+
 # Delegate lazy-loading of data resources to the util module
 _LAZY_RESOURCES = list(util._ensure_registry().keys())
 
@@ -45,5 +106,7 @@ def __dir__():  # pragma: no cover
 
 __all__ = [
     "Convolve", "ChebyConvolve", "DyConvolve", "SGMA", "ModeTable", "functions",
-    "VFKernel", "ChebyKernel", "impulse", "get_klu_dll", "get_cholmod_dll", "estimate_spectral_bound"
+    "VFKernel", "ChebyKernel", "VFModel", "ChebyModel",
+    "impulse", "get_klu_dll", "get_cholmod_dll", "estimate_spectral_bound",
+    "load_ply_laplacian", "load_ply_xyz", "list_graphs"
 ] + _LAZY_RESOURCES

sgwt/chebyconv.py

@@ -25,37 +25,58 @@ class ChebyConvolve:
     """
 
     def __init__(self, L: csc_matrix) -> None:
+        """Initialize the Chebyshev convolution context.
+
+        Parameters
+        ----------
+        L : csc_matrix
+            Sparse Graph Laplacian. The spectral bound is estimated
+            automatically via :func:`~sgwt.util.estimate_spectral_bound`.
+        """
         self.n_vertices = L.shape[0]
         self.spectrum_bound = estimate_spectral_bound(L)
         self.chol = CholWrapper(L)
         self._cached_M_ptr = None
         self._cached_spectrum_bound = None
 
     def __enter__(self) -> "ChebyConvolve":
+        """Start CHOLMOD and perform symbolic factorization.
+
+        Returns
+        -------
+        ChebyConvolve
+            This context manager instance.
+        """
         self.chol.start()
         self.chol.sym_factor()
         return self
 
     def __exit__(self, exc_type, exc_val, exc_tb):
+        """Free CHOLMOD resources and cached recurrence matrix."""
         if self._cached_M_ptr is not None:
             self.chol.free_sparse(self._cached_M_ptr)
         self.chol.free_factor(self.chol.fact_ptr)
         self.chol.finish()
 
-    def _get_recurrence_matrix(self, spectrum_bound: float):
-        """Get cached recurrence matrix M = (2/λ_max)L - I."""
+    def _get_recurrence_matrix(self, spectrum_bound: float, min_lambda: float = 0.0):
+        """Get cached recurrence matrix M = alpha*L + beta*I for domain [min_lambda, spectrum_bound]."""
+        cache_key = (spectrum_bound, min_lambda)
         if self._cached_M_ptr is not None:
-            if self._cached_spectrum_bound == spectrum_bound:
+            if self._cached_spectrum_bound == cache_key:
                 return self._cached_M_ptr
             self.chol.free_sparse(self._cached_M_ptr)
-
+
+        range_ = spectrum_bound - min_lambda
+        alpha = 2.0 / range_
+        beta = -(spectrum_bound + min_lambda) / range_
+
         EYE = self.chol.speye(self.n_vertices, self.n_vertices)
         try:
             self._cached_M_ptr = self.chol.add(
                 byref(self.chol.A), EYE,
-                alpha=2.0 / spectrum_bound, beta=-1.0
+                alpha=alpha, beta=beta
             )
-            self._cached_spectrum_bound = spectrum_bound
+            self._cached_spectrum_bound = cache_key
             return self._cached_M_ptr
         finally:
             self.chol.free_sparse(EYE)
@@ -122,7 +143,7 @@ def _convolve_real(self, B: np.ndarray, C: ChebyKernel) -> np.ndarray:
             if n_order == 1:
                 return W
 
-            M = self._get_recurrence_matrix(C.spectrum_bound)
+            M = self._get_recurrence_matrix(C.spectrum_bound, C.min_lambda)
             T_km1 = self.chol.allocate_dense(n_vertex, n_signals)
 
             # T_1(L)B = M @ B
@@ -184,20 +205,20 @@ def convolve_multi(self, B: np.ndarray, kernels: list) -> list:
 
         n_vertex, n_signals = B.shape
 
-        # Group by spectrum_bound
+        # Group by (spectrum_bound, min_lambda)
         from collections import defaultdict
         groups = defaultdict(list)
         for idx, k in enumerate(kernels):
-            groups[k.spectrum_bound].append((idx, k))
-        
+            groups[(k.spectrum_bound, k.min_lambda)].append((idx, k))
+
         results = [None] * len(kernels)
-        
-        for bound, group in groups.items():
+
+        for (bound, min_lam), group in groups.items():
             max_order = max(k.C.shape[0] for _, k in group)
-            
+
             # Stack all coefficients for this group for vectorized application
             # Compute terms and apply via tensor contraction
-            T_terms = self._compute_terms(B, max_order, bound)
+            T_terms = self._compute_terms(B, max_order, bound, min_lam)
             T_stack = np.stack(T_terms, axis=0)  # (max_order, n_vertex, n_signals)
 
             for idx, kernel in group:
@@ -208,20 +229,20 @@ def convolve_multi(self, B: np.ndarray, kernels: list) -> list:
 
         return results
 
-    def _compute_terms(self, B: np.ndarray, max_order: int, bound: float) -> list:
+    def _compute_terms(self, B: np.ndarray, max_order: int, bound: float, min_lambda: float = 0.0) -> list:
         """Compute Chebyshev terms T_k(L)B for k = 0..max_order-1."""
         n_vertex, n_signals = B.shape
         terms = []
-        
+
         B_chol = byref(self.chol.numpy_to_chol_dense(B))
         T_km2, T_km1 = None, None
-        
+
         try:
             T_km2 = self.chol.copy_dense(B_chol)
             terms.append(self.chol.chol_dense_to_numpy(T_km2).copy())
-            
+
             if max_order > 1:
-                M = self._get_recurrence_matrix(bound)
+                M = self._get_recurrence_matrix(bound, min_lambda)
                 T_km1 = self.chol.allocate_dense(n_vertex, n_signals)
 
                 self.chol.sdmult(M, T_km2, T_km1, alpha=1.0, beta=0.0)