ENH: Implement faster Hermite functions in Cython (#5)

* tmp:
- started serial implementation of Hermite functions in Cython

* BLD:
- updated `requirements`
- updated `pyproject.toml` to handle Cython build

* ENH:
- switched NumPy and Numba Hermite basis function to faster implementation that avoids the expensive use of the logsumexp trick and rather uses a factor compensated recursion
- added function to compute a single Hermite function

* ENH:
- added Cythonized parallel implementation of the Hermite functions
- added a `setup.py` to specify the Cython build
- adapted the interface and "deprecated" the NumPy- and Numba-based implementations of the Hermite functions

* TST:
- rewrote generation of the symbolic reference data for the Hermite function tests to keep it more flexible
- created new test files for the Hermite functions

* TST:
- improved tests of NumPy- and Numba-based Hermite functions with new reference test files
- added test for single Hermite function computation
- made numerical tolerances of Hermite function tests way stricter
- extended orthonormality test from order 500 to 1_000

* DOC:
- updated examples to use Cython implementation of the Hermite functions
- transitioned performance testing of Hermite functions to Jupyter notebook and created a nice plot from there
- updated all example plots
- updated equations to make them more unified

* DOC:
- updated `README` with two new equations and removed the logsumexp trick link because it's no longer how the algorithm works

* TST:
- integrated single and multi-threaded Cython Hermite-functions into tests

* TST:
- added exception handling test for all Hermite function implementations

* DOC:
- updated number of significant digits in test description

* DOC:
- added reference to the publication the package is based on (right now)

* DOC:
- made DOI to reference in `README` a link ?

* DOC:
- fixed broken link to DOI of reference in `README` ?

.gitignore

@@ -129,7 +129,4 @@ dmypy.json
 .pyre/
 
 # VScode
-.vscode/
-
-# Notebooks
-*.ipynb
+.vscode/

README.rst

@@ -33,25 +33,32 @@ The Hermite functions are defined as
 with the Hermite polynomials
 
 .. image:: docs/hermite_functions/equations/DilatedHermitePolynomials.png
-    :width: 660px
+    :width: 681px
     :align: left
 
 By making use of logarithm tricks, the evaluation that might involve infinitely high
 polynomial values and at the same time infinitely small Gaussians - that are on top of
 that scaled by an infinitely high factorial - can be computed safely and yield accurate
 results.
 
-For doing so, the equation is rewritten in logarithmic form as
+For doing so, the relation between the dilated and the non-dilated Hermite functions
 
-.. image:: docs/hermite_functions/equations/LogDilatedHermiteFunctions.png
-    :width: 863px
+.. image:: docs/hermite_functions/equations/HermiteFunctions_UndilatedToDilated.png
+    :width: 321px
     :align: left
 
-where the evaluation of the natural logarithm of the Hermite polynomials is achieved by
-making use of the
-`logsumexp trick <https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.logsumexp.html>`_.
+and the recurrence relation for the Hermite functions
 
-This approach is tested against a symbolic evaluation with ``sympy`` that uses 100
+.. image:: docs/hermite_functions/equations/HermiteFunctions_RecurrenceRelation.png
+    :width: 699px
+    :align: left
+
+are used, but not directly. Instead, the latest evaluated Hermite function is kept at a
+value of either -1, 0, or +1 during the recursion and the logarithm of a correction
+factor is tracked and applied when the respective Hermite function is finally evaluated
+and stored. This approach is based on [1_].
+
+This approach is tested against a symbolic evaluation with ``sympy`` that uses 200
 digits of precision and it can be shown that even orders as high as 2,000 can still be
 computed even though neither the polynomial, the Gaussian nor the factorial can be
 evaluated for this anymore. The factorial for example would already have overflown for
@@ -64,3 +71,9 @@ orders of 170 in ``float64``-precision.
 As a sanity check, their orthogonality is part of the tests together with a test for
 the fact that the absolute values of the Hermite functions for real input cannot exceed
 the value :math:`\frac{\pi^{-\frac{1}{4}}}{\sqrt{\alpha}}`.
+
+References
+----------
+.. [1] Bunck B. F., A fast algorithm for evaluation of normalized Hermite
+    functions, BIT Numer Math (2009), 49, pp. 281–295, DOI:
+    `<https://doi.org/10.1007/s10543-009-0216-1>`_

examples/01_hermite_functions.py

@@ -17,7 +17,7 @@
 # the x-values to evaluate the Hermite functions
 X_FROM = -5.0
 X_TO = 5.0
-NUM_X = 1_001
+NUM_X = 10_001
 
 # the scaling factors alpha to use
 ALPHAS = [0.5, 1.0, 2.0]
@@ -48,21 +48,18 @@
     for idx_alpha, alpha in enumerate(ALPHAS):
         # the Hermite functions are computed and plotted
         hermite_basis = hermite_function_basis(
-            x=x_values,
-            n=ORDERS,
-            alpha=alpha,
-            jit=True,
+            x=x_values, n=ORDERS, alpha=alpha, workers=-1
         )
 
         # NOTE: x-axis are plotted for orientation
         for idx_order in range(0, ORDERS + 1):
-            ax[idx_alpha].axhline(
+            ax[idx_alpha].axhline(  # type: ignore
                 y=idx_order * OFFSET,
                 color="black",
                 linewidth=0.5,
                 zorder=idx_order * 2,
             )
-            ax[idx_alpha].plot(
+            ax[idx_alpha].plot(  # type: ignore
                 x_values,
                 hermite_basis[::, idx_order] + idx_order * OFFSET,
                 color=colors[idx_order],
@@ -71,27 +68,27 @@
             )
 
         # the title, grid, x-labels, and ticks are set
-        ax[idx_alpha].set_title(
+        ax[idx_alpha].set_title(  # type: ignore
             r"$\alpha$" + f"= {alpha:.1f}",
             fontsize=16,
         )
-        ax[idx_alpha].set_xlabel(
+        ax[idx_alpha].set_xlabel(  # type: ignore
             r"$x$",
             fontsize=16,
             labelpad=10,
         )
-        ax[idx_alpha].tick_params(axis="both", which="major", labelsize=14)
-        ax[idx_alpha].grid(which="major", axis="both")
-        ax[idx_alpha].set_xlim(X_FROM, X_TO)
+        ax[idx_alpha].tick_params(axis="both", which="major", labelsize=14)  # type: ignore
+        ax[idx_alpha].grid(which="major", axis="both")  # type: ignore
+        ax[idx_alpha].set_xlim(X_FROM, X_TO)  # type: ignore
 
         # for the first plot, a y-label and a legend are added
         if idx_alpha == 0:
-            ax[idx_alpha].set_ylabel(
+            ax[idx_alpha].set_ylabel(  # type: ignore
                 r"$\psi_{n}^{\left(\alpha\right)}\left(x\right)$",
                 fontsize=16,
                 labelpad=10,
             )
-            ax[idx_alpha].legend(
+            ax[idx_alpha].legend(  # type: ignore
                 loc="upper left",
                 fontsize=14,
                 frameon=False,

examples/02_hermite_function_stability.py

@@ -18,7 +18,7 @@
 # the x-values to evaluate the Hermite functions
 X_FROM = -65.0
 X_TO = 65.0
-NUM_X = 40_001
+NUM_X = 100_001
 
 # the scaling factor alpha to use
 ALPHA = 1.0
@@ -50,7 +50,7 @@
         x=x_values,
         n=max(ORDERS),
         alpha=ALPHA,
-        jit=True,
+        workers=-1,
     )
 
     # ... and the individual Hermite functions of interest are plotted

pyproject.toml

@@ -1,5 +1,10 @@
 [build-system]
-requires = ["setuptools>=65", "wheel"]
+requires = [
+    "setuptools>=65",
+    "wheel",
+    "numpy>=1.21.0,<2.0.0",
+    "Cython>=3.0.10,<3.1.0",
+]
 build-backend = "setuptools.build_meta"
 
 [project]
@@ -11,12 +16,32 @@ keywords = ["Fourier Transform", "Hermite Functions", "Robust", "Noise", "Outlie
 license = {file = "LICENSE"}
 dynamic = ["version", "readme", "dependencies", "optional-dependencies"]
 
+[tool.setuptools]
+include-package-data = true
+package-data = {"*" = ["AUTHORS.txt", "VERSION.txt"]}
+
 [tool.setuptools.dynamic]
 version = {file = "robust_hermite_ft/VERSION.txt"}
 readme = {file = ["README.rst"]}
 dependencies = {file = "requirements/base.txt"}
 optional-dependencies = {fast = {file = "requirements/fast.txt"}, dev = {file = "requirements/dev.txt"}, examples = {file = "requirements/examples.txt"}}
 
+[tool.isort]
+profile = "black"
+multi_line_output = 3
+
+[tool.black]
+target-version = [
+    "py38",
+    "py39",
+    "py310",
+    "py311",
+    "py312",
+]
+
+[tool.cython-lint]
+max-line-length = 120
+
 [tool.ruff]
 # Enable pycodestyle (`E`), Pyflakes (`F`) checks.
 select = ["E", "F"]

requirements/base.txt

@@ -1,2 +1,3 @@
 numpy>=1.21.0,<2.0.0
+psutil>=5.8.0
 scipy>=1.7.0

requirements/dev.txt

@@ -1,5 +1,7 @@
 black
 coverage
+cython>=3.0.10
+cython-lint>=0.16.0
 esbonio
 ipympl
 jupyter

robust_hermite_ft/__init__.py

@@ -10,7 +10,11 @@
 
 import os as _os
 
-from .hermite_functions import hermite_function_basis  # noqa: F401
+from .hermite_functions import (  # noqa: F401
+    hermite_function_basis,
+    single_hermite_function,
+    slow_hermite_function_basis,
+)
 
 # === Package Metadata ===
 
@@ -22,4 +26,3 @@
 
 with open(_VERSION_FILE_PATH, "r") as version_file:
     __version__ = version_file.read().strip()
-

robust_hermite_ft/hermite_functions/.gitignore

@@ -0,0 +1,2 @@
+_c_hermite.c
+_c_hermite.html

robust_hermite_ft/hermite_functions/__init__.py

@@ -21,4 +21,8 @@
 
 # === Imports ===
 
-from ._interface import hermite_function_basis  # noqa: F401
+from ._interface import (  # noqa: F401
+    hermite_function_basis,
+    single_hermite_function,
+    slow_hermite_function_basis,
+)