Skip to content

Commit

Permalink
[PKG] Initial implementation of Hermite function (v0.0b1) (#1)
Browse files Browse the repository at this point in the history
* BLD:
- start `README`
- set up `.gitignore`
- made package buildable with dynamic version, readme, and requirements

* BUG:
- fixed wrong italics in `README`

* BLD:
- added full author name to `LICENSE`

* BLD:
- restructured requirements
- added `examples`  requirements

* DOC:
- tested enriched `README` with Hermite function pictures, equations, and links

* DOC:
- improved visuals of `README`?

* BLD:
- reverted changes in `README` to make it more readable ?

* BLD:
- fixed wrong image links in `README`

* BUG:
- fixed broken image link in `README`

* DOC:
- specified image widths in `README` ?
- fixed `\alpha` not rendering ?

* DOC:
- specified image sizes of equations ?
- fixed `\alpha` not rendering properly ?

* DOC:
- proved stability of Hermite functions and added new plot to `README`

* DOC:
- rewrote "come back later" to "check back later"

* BLD:
- updated requirements with stricter version of `numpy` and removal of `proplot`

* DEV:
- implemented a core function to evaluate the natural logarithm of Hermite polynomials (without interface with input validation)
- implemented a core function to evaluate arbitrary Hermite functions (without interface with input validation)

* TST:
- added test and test file generator for core function that computes the natural logarithm of Hermite polynomials

* DOC:
- added two examples of how to compute Hermite functions in a numerically stable fashion

* TST:
- added orthonormality and boundedness test for the core hermite function
- mentioned the test in the `README`
  • Loading branch information
MothNik authored Jul 6, 2024
1 parent e637e69 commit ab1c6be
Show file tree
Hide file tree
Showing 34 changed files with 1,365 additions and 1 deletion.
135 changes: 135 additions & 0 deletions .gitignore
Original file line number Diff line number Diff line change
@@ -0,0 +1,135 @@
# Byte-compiled / optimized / DLL files
__pycache__/
*.py[cod]
*$py.class

# C extensions
*.so

# Distribution / packaging
.Python
build/
develop-eggs/
dist/
downloads/
eggs/
.eggs/
lib/
lib64/
parts/
sdist/
var/
wheels/
pip-wheel-metadata/
share/python-wheels/
*.egg-info/
.installed.cfg
*.egg
MANIFEST

# PyInstaller
# Usually these files are written by a python script from a template
# before PyInstaller builds the exe, so as to inject date/other infos into it.
*.manifest
*.spec

# Installer logs
pip-log.txt
pip-delete-this-directory.txt

# Unit test / coverage reports
htmlcov/
.tox/
.nox/
.coverage
.coverage.*
.cache
nosetests.xml
coverage.xml
*.cover
*.py,cover
.hypothesis/
.pytest_cache/

# Translations
*.mo
*.pot

# Django stuff:
*.log
local_settings.py
db.sqlite3
db.sqlite3-journal

# Flask stuff:
instance/
.webassets-cache

# Scrapy stuff:
.scrapy

# Sphinx documentation
docs/_build/

# PyBuilder
target/

# Jupyter Notebook
.ipynb_checkpoints

# IPython
profile_default/
ipython_config.py

# pyenv
.python-version

# pipenv
# According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control.
# However, in case of collaboration, if having platform-specific dependencies or dependencies
# having no cross-platform support, pipenv may install dependencies that don't work, or not
# install all needed dependencies.
#Pipfile.lock

# PEP 582; used by e.g. github.com/David-OConnor/pyflow
__pypackages__/

# Celery stuff
celerybeat-schedule
celerybeat.pid

# SageMath parsed files
*.sage.py

# Environments
.env
.venv*
env/
venv/
ENV/
env.bak/
venv.bak/

# Spyder project settings
.spyderproject
.spyproject

# Rope project settings
.ropeproject

# mkdocs documentation
/site

# mypy
.mypy_cache/
.dmypy.json
dmypy.json

# Pyre type checker
.pyre/

# VScode
.vscode/

# Notebooks
*.ipynb
2 changes: 1 addition & 1 deletion LICENSE
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
MIT License

Copyright (c) 2024 Niklas Z
Copyright (c) 2024 Niklas Zell

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
Expand Down
60 changes: 60 additions & 0 deletions README.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1,60 @@
``robust_hermite_ft``
=====================

You want to compute the Fourier transform of a signal, but your signal can be corrupted
by outliers? If so, this package is for you even though you will have to say goodbye to
the *"fast"* in *Fast Fourier Transform* 🏃🙅‍♀️

🏗️🚧 👷👷‍♂️👷‍♀️🏗️🚧

Currently under construction. Please check back later.

〰️ Hermite functions
---------------------

Being the eigenfunctions of the Fourier transform, Hermite functions are an excellent
candidate as the basis functions for a Least Squares Regression approach to the Fourier
transform. However, their evaluation can be a bit tricky.

The module ``hermite_functions`` offers a numerically stable way to evaluate Hermite
functions or arbitrary order :math:`n` and argument - that can be scaled with a factor
:math:`{\alpha}`

.. image:: docs/hermite_functions/DilatedHermiteFunctions_DifferentScales.png
:width: 1000px
:align: center

The Hermite functions are defined as

.. image:: docs/hermite_functions/equations/DilatedHermiteFunctions.png
:width: 500px
: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

.. image:: docs/hermite_functions/equations/LogDilatedHermiteFunctions.png
:width: 902px
: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>`_.

This approach is tested against a symbolic evaluation with ``sympy`` that uses 100
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
170 in ``float64``-precision.

.. image:: docs/hermite_functions/DilatedHermiteFunctions_Stability.png
:width: 1000px
:align: center

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}}`.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
109 changes: 109 additions & 0 deletions examples/01_hermite_functions.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,109 @@
"""
This script generates plots of different Hermite functions with different scales.
"""

# === Imports ===

import os

import numpy as np
from matplotlib import pyplot as plt

from src.hermite_functions import _dilated_hermite_function_basis

# === Constants ===

# the x-values to evaluate the Hermite functions
X_FROM = -5.0
X_TO = 5.0
NUM_X = 1_001

# the scaling factors alpha to use
ALPHAS = [0.5, 1.0, 2.0]
# the orders of the Hermite functions to plot
ORDERS = 6
# the offset between the individual Hermite functions
OFFSET = -2.0

# the path where to store the plot and its resolution
PLOT_FILEPATH = "../docs//hermite_functions/01_DilatedHermiteFunctions.png"
DPI = 300

# === Main ===

if __name__ == "__main__":

fig, ax = plt.subplots(
ncols=len(ALPHAS),
sharex=True,
sharey=True,
figsize=(12, 8),
)

# the Hermite functions are evaluated and plotted for each scaling factor alpha
x_values = np.linspace(start=X_FROM, stop=X_TO, num=NUM_X)

colors = plt.cm.winter_r(np.linspace(0, 1, ORDERS + 1)) # type: ignore
for idx_alpha, alpha in enumerate(ALPHAS):
# the Hermite functions are computed and plotted
hermite_basis = _dilated_hermite_function_basis(
x=x_values,
n=ORDERS,
alpha=alpha,
)

# NOTE: x-axis are plotted for orientation
for idx_order in range(0, ORDERS + 1):
ax[idx_alpha].axhline(
y=idx_order * OFFSET,
color="black",
linewidth=0.5,
zorder=idx_order * 2,
)
ax[idx_alpha].plot(
x_values,
hermite_basis[::, idx_order] + idx_order * OFFSET,
color=colors[idx_order],
zorder=(idx_order + 1) * 2,
label=f"n = {idx_order}",
)

# the title, grid, x-labels, and ticks are set
ax[idx_alpha].set_title(
r"$\alpha$" + f"= {alpha:.1f}",
fontsize=16,
)
ax[idx_alpha].set_xlabel(
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)

# for the first plot, a y-label and a legend are added
if idx_alpha == 0:
ax[idx_alpha].set_ylabel(
r"$\psi_{n}^{\left(\alpha\right)}\left(x\right)$",
fontsize=16,
labelpad=10,
)
ax[idx_alpha].legend(
loc="upper left",
fontsize=14,
frameon=False,
)

plt.tight_layout()

# the plot is saved ...
plt.savefig(
os.path.join(os.path.dirname(__file__), PLOT_FILEPATH),
dpi=DPI,
bbox_inches="tight",
)

# ... and shown
plt.show()
Loading

0 comments on commit ab1c6be

Please sign in to comment.