Skip to content

Improving Spectral #53

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
Jacob-Stevens-Haas merged 14 commits into from
Apr 2, 2025
Merged

Improving Spectral #53

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
14 commits
Select commit Hold shift + click to select a range
957caaf
spitballing
pavelkomarov Jan 16, 2025
859a87a
my interval was way wrong
pavelkomarov Jan 16, 2025
6766e52
updated spectral-derivatives such that some of the checks I had in th…
pavelkomarov Jan 18, 2025
ae65554
added some road signs
pavelkomarov Jan 18, 2025
2688f07
updated to spectral-derivatives 0.7
pavelkomarov Jan 22, 2025
cb59792
this test code encourages not-thorough testing
pavelkomarov Jan 22, 2025
0226b3b
fixed caching test
pavelkomarov Jan 25, 2025
b1154da
added tests for chebyshev case and removed caching check for spectral…
pavelkomarov Jan 27, 2025
cb184a0
added Chebyshev examples
pavelkomarov Jan 27, 2025
8bb4df7
removed print statement
pavelkomarov Jan 27, 2025
18926e3
bumped version number
pavelkomarov Feb 4, 2025
d11d1a7
added axis param to tests so they pass now
pavelkomarov Feb 6, 2025
22865d3
improved docstring to account for the fact filtering in chebyshev bas…
pavelkomarov Feb 28, 2025
825e6dd
added back a newline
pavelkomarov Apr 2, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
55 changes: 32 additions & 23 deletions derivative/dglobal.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,40 +7,50 @@
from scipy import interpolate, sparse
from scipy.special import legendre
from sklearn.linear_model import Lasso
from specderiv import cheb_deriv, fourier_deriv


@register("spectral")
class Spectral(Derivative):
def __init__(self, **kwargs):
def __init__(self, order=1, axis=0, basis='fourier', filter=None):
"""
Compute the numerical derivative by first computing the FFT. In Fourier space, derivatives are multiplication
by i*phase; compute the IFFT after.

Args:
**kwargs: Optional keyword arguments.

Compute the numerical derivative by spectral methods. In Fourier space, derivatives are multiplication
by i*phase; compute the inverse transform after. Use either Fourier modes of Chebyshev polynomials as
the basis.

Keyword Args:
filter: Optional. A function that takes in frequencies and outputs weights to scale the coefficient at
the input frequency in Fourier space. Input frequencies are the discrete fourier transform sample
frequencies associated with the domain variable. Look into python signal processing resources in
scipy.signal for common filters.

order (int): order of the derivative, defaults to 1st order
axis (int): the dimension of the data along which to differentiate, defaults to first dimension
basis (str): 'fourier' or 'chebyshev', the set of basis functions to use for differentiation
Note `basis='fourier'` assumes your function is periodic and sampled over a period of its domain,
[a, b), and `basis='chebyshev'` assumes your function is sampled at cosine-spaced points on the
domain [a, b].
filter: Optional. A function that takes in basis function indices and outputs weights, which scale
the corresponding modes in the basis-space interpolation before derivatives are taken, e.g.
`lambda k: k<10` will keep only the first ten modes. With the Fourier basis, k corresponds to
wavenumbers, so common filters from scipy.signal can be used. In the Chebyshev basis, modes do
not directly correspond to frequencies, so high frequency noise can not be separated quite as
cleanly, however it still may be helpful to dampen higher modes.
"""
# Filter function. Default: Identity filter
self.filter = kwargs.get('filter', np.vectorize(lambda f: 1))
self._x_hat = None
self._freq = None

def _dglobal(self, t, x):
self._x_hat = np.fft.fft(x)
self._freq = np.fft.fftfreq(t.size, d=(t[1] - t[0]))
self.order = order
self.axis = axis
if basis not in ['chebyshev', 'fourier']:
raise ValueError("Only chebyshev and fourier bases are allowed.")
self.basis = basis
self.filter = filter

@_memoize_arrays(1) # the memoization is 1 deep, as in only remembers the most recent args
def _global(self, t, x):
if self.basis == 'chebyshev':
return cheb_deriv(x, t, self.order, self.axis, self.filter)
Copy link
Contributor Author

@pavelkomarov pavelkomarov Feb 15, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It struck me just now that filtering by frequency on cosine-spaced points isn’t totally sound: high-frequency variations look lower frequency around the edges where sample density is higher.

Copy link
Contributor Author

@pavelkomarov pavelkomarov Feb 20, 2025 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I've been pondering this more pavelkomarov/spectral-derivatives#14. Basically what I'm currently doing is keeping filter(k) amount of the $k^{th}$ Chebyshev polynomial, versus in the Fourier case keeping filter(k) amount of the $k^{th}$ wavenumber. This may be the best thing to do, most spiritually-aligned with "spectral". I so far haven't thought of a better alternative. But it's counter-intuitive: Because the Chebyshev polynomials have steeper slope at the edges of the domain, we can fit high-frequencies better there. So if the assumption is that noise is high frequency and needs to be separated from the data this way, we can't filter it as consistently using this basis.

Copy link
Contributor Author

@pavelkomarov pavelkomarov Mar 10, 2025 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I've made a notebook to empirically compare spectral filtering in the Fourier and Chebyshev bases, and pulled in a little of Floris' PyNumDiff to compare with Savitzky-Golay filtering + Gaussian smoothing. Fourier spectral is the gold standard, but even it can be fragile. Chebyshev can work okay (not great but really not worse than savgoldiff) for 1st and 2nd derivatives in the presence of a little noise but starts to blow up at the edges of the domain beyond that.

else: # self.basis == 'fourier'
return fourier_deriv(x, t, self.order, self.axis, self.filter)

def compute(self, t, x, i):
return next(self.compute_for(t, x, [i]))

def compute_for(self, t, x, indices):
self._dglobal(t, x)
res = np.fft.ifft(1j * 2 * np.pi * self._freq * self.filter(self._freq) * self._x_hat).real
res = self._global(t, x) # cached
for i in indices:
yield res[i]

Expand Down Expand Up @@ -212,7 +222,6 @@ def __init__(self, alpha=None):
"""
self.alpha = alpha


@_memoize_arrays(1)
def _global(self, t, z, alpha):
if alpha is None:
Expand Down
22 changes: 13 additions & 9 deletions derivative/differentiation.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@ def _gen_method(x, t, kind, axis, **kwargs):
return methods.get(kind)(**kwargs)


def dxdt(x, t, kind=None, axis=1, **kwargs):
def dxdt(x, t, kind=None, axis=0, **kwargs):
"""
Compute the derivative of x with respect to t along axis using the numerical derivative specified by "kind". This is
the functional interface of the Derivative class.
Expand All @@ -35,7 +35,7 @@ def dxdt(x, t, kind=None, axis=1, **kwargs):
x (:obj:`ndarray` of float): Ordered measurement values.
t (:obj:`ndarray` of float): Ordered measurement times.
kind (string): Derivative method name (see available kinds).
axis ({0,1}): Axis of x along which to differentiate. Default 1.
axis ({0,1}): Axis of x along which to differentiate. Default 0.
**kwargs: Keyword arguments for the derivative method "kind".

Available kinds
Expand All @@ -56,7 +56,7 @@ def dxdt(x, t, kind=None, axis=1, **kwargs):
return method.d(x, t, axis=axis)


def smooth_x(x, t, kind=None, axis=1, **kwargs):
def smooth_x(x, t, kind=None, axis=0, **kwargs):
"""
Compute the smoothed version of x given t along axis using the numerical
derivative specified by "kind". This is the functional interface of
Expand All @@ -71,7 +71,7 @@ def smooth_x(x, t, kind=None, axis=1, **kwargs):
x (:obj:`ndarray` of float): Ordered measurement values.
t (:obj:`ndarray` of float): Ordered measurement times.
kind (string): Derivative method name (see available kinds).
axis ({0,1}): Axis of x along which to differentiate. Default 1.
axis ({0,1}): Axis of x along which to differentiate. Default 0.
**kwargs: Keyword arguments for the derivative method "kind".

Available kinds
Expand Down Expand Up @@ -100,7 +100,7 @@ def compute(self, t, x, i):
"""
Compute the derivative of one-dimensional data x with respect to t at the index i of x, (dx/dt)[i].

Computation of a derivative should fail explicitely if the implementation is unable to compute a derivative at
Computation of a derivative should fail explicitly if the implementation is unable to compute a derivative at
the desired index. Used for global differentiation methods, for example.

This requires that x and t have equal lengths >= 2, and that the index i is a valid index.
Expand Down Expand Up @@ -174,7 +174,7 @@ def compute_x_for(self, t, x, indices):
for i in indices:
yield self.compute_x(t, x, i)

def d(self, X, t, axis=1):
def d(self, X, t, axis=0):
"""
Compute the derivative of measurements X taken at times t.

Expand All @@ -184,7 +184,7 @@ def d(self, X, t, axis=1):
Args:
X (:obj:`ndarray` of float): Ordered measurements values. Multiple measurements allowed.
t (:obj:`ndarray` of float): Ordered measurement times.
axis ({0,1}). axis of X along which to differentiate. default 1.
axis ({0,1}). axis of X along which to differentiate. default 0.

Returns:
:obj:`ndarray` of float: Returns dX/dt along axis.
Expand All @@ -202,7 +202,7 @@ def d(self, X, t, axis=1):
return _restore_axes(dX, axis, flat)


def x(self, X, t, axis=1):
def x(self, X, t, axis=0):
"""
Compute the smoothed X values from measurements X taken at times t.

Expand All @@ -212,7 +212,7 @@ def x(self, X, t, axis=1):
Args:
X (:obj:`ndarray` of float): Ordered measurements values. Multiple measurements allowed.
t (:obj:`ndarray` of float): Ordered measurement times.
axis ({0,1}). axis of X along which to smooth. default 1.
axis ({0,1}). axis of X along which to smooth. default 0.

Returns:
:obj:`ndarray` of float: Returns dX/dt along axis.
Expand All @@ -228,6 +228,8 @@ def x(self, X, t, axis=1):


def _align_axes(X, t, axis) -> tuple[NDArray, tuple[int, ...]]:
"""Reshapes the data so the derivative always happens along axis 1.
"""
X = np.array(X)
orig_shape = X.shape
# By convention, differentiate axis 1
Expand All @@ -244,6 +246,8 @@ def _align_axes(X, t, axis) -> tuple[NDArray, tuple[int, ...]]:


def _restore_axes(dX: NDArray, axis: int, orig_shape: tuple[int, ...]) -> NDArray:
"""Undo the operation of _align_axes, so data can be returned in its original shape
"""
if len(orig_shape) == 1:
return dX.flatten()
else:
Expand Down
101 changes: 61 additions & 40 deletions docs/notebooks/Examples.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,20 +1,5 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"ExecuteTime": {
"end_time": "2020-05-25T22:56:27.818221Z",
"start_time": "2020-05-25T22:56:27.413152Z"
},
"tags": []
},
"outputs": [],
"source": [
"%matplotlib inline"
]
},
{
"cell_type": "code",
"execution_count": null,
Expand Down Expand Up @@ -65,17 +50,16 @@
"outputs": [],
"source": [
"def plot_example(diff_method, t, data_f, res_f, sigmas, y_label=None):\n",
" '''\n",
" Utility function for concise plotting of examples.\n",
" '''\n",
" '''Utility function for concise plotting of examples.'''\n",
" fig, axes = plt.subplots(1, len(sigmas), figsize=[len(sigmas)*4, 3])\n",
" \n",
" # Compute the derivative\n",
" res = diff_method.d(np.vstack([data_f(t, s) for s in sigmas]), t)\n",
" res = diff_method.d(np.vstack([data_f(t, s) for s in sigmas]), t, axis=1)\n",
" for i, s in enumerate(sigmas):\n",
" axes[i].plot(t, res[i])\n",
" axes[i].plot(t, res_f(t))\n",
" axes[i].set_title(\"Noise: $\\sigma$={}\".format(s))\n",
" axes[i].plot(t, res[i])\n",
" axes[i].set_title(r\"Noise: $\\sigma$={}\".format(s))\n",
" axes[i].set_ylim([-1.25,1.3])\n",
" if y_label:\n",
" axes[0].set_ylabel(y_label, fontsize=12)"
]
Expand Down Expand Up @@ -133,7 +117,7 @@
"\n",
"fig,ax = plt.subplots(1, figsize=[5,3])\n",
"kind = FiniteDifference(k=1)\n",
"ax.plot(t, kind.d(x,t))"
"ax.plot(t, kind.d(x,t));"
]
},
{
Expand All @@ -159,7 +143,7 @@
"from derivative import dxdt\n",
"\n",
"fig,ax = plt.subplots(1, figsize=[5,3])\n",
"ax.plot(t, dxdt(x, t, \"finite_difference\", k=1))"
"ax.plot(t, dxdt(x, t, \"finite_difference\", k=1));"
]
},
{
Expand Down Expand Up @@ -203,10 +187,10 @@
"sigmas = [0, 0.01, 0.1]\n",
"fig, ax = plt.subplots(1, len(sigmas), figsize=[len(sigmas)*4, 3])\n",
"\n",
"t = np.linspace(0, 2*np.pi, 50)\n",
"t = np.linspace(0, 2*np.pi, 50, endpoint=False)\n",
"for axs, s in zip(ax, sigmas): \n",
" axs.scatter(t, noisy_sin(t, s))\n",
" axs.set_title(\"Noise: $\\sigma$={}\".format(s))"
" axs.set_title(r\"Noise: $\\sigma$={}\".format(s))"
]
},
{
Expand Down Expand Up @@ -295,7 +279,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"### Spectral method\n",
"### Spectral method - Fourier basis\n",
"Add your own filter!"
]
},
Expand All @@ -312,12 +296,35 @@
"outputs": [],
"source": [
"no_filter = derivative.Spectral()\n",
"yes_filter = derivative.Spectral(filter=np.vectorize(lambda f: 1 if abs(f) < 0.5 else 0))\n",
"yes_filter = derivative.Spectral(filter=np.vectorize(lambda k: 1 if abs(k) < 3 else 0))\n",
"\n",
"plot_example(no_filter, t, noisy_sin, np.cos, sigmas, 'No filter')\n",
"plot_example(yes_filter, t, noisy_sin, np.cos, sigmas, 'Low-pass filter')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Spectral method - Chebyshev basis\n",
"\n",
"Now let's do with the Chebyshev basis, which requires cosine-spaced points on [a, b] rather than equispaced points on [a, b)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"t_cos = np.cos(np.pi * np.arange(50) / 49) * np.pi + np.pi # choose a = 0, b = 2*pi\n",
"no_filter = derivative.Spectral(basis='chebyshev')\n",
"yes_filter = derivative.Spectral(basis='chebyshev', filter=np.vectorize(lambda k: 1 if abs(k) < 6 else 0))\n",
"\n",
"plot_example(no_filter, t_cos, noisy_sin, np.cos, sigmas, 'No filter')\n",
"plot_example(yes_filter, t_cos, noisy_sin, np.cos, sigmas, 'Low-pass filter')"
]
},
{
"cell_type": "markdown",
"metadata": {},
Expand Down Expand Up @@ -394,7 +401,7 @@
"outputs": [],
"source": [
"def noisy_abs(t, sigma):\n",
" '''Sine with gaussian noise.'''\n",
" '''Abs with gaussian noise.'''\n",
" np.random.seed(17)\n",
" return np.abs(t) + np.random.normal(loc=0, scale=sigma, size=x.shape)\n",
"\n",
Expand All @@ -406,7 +413,7 @@
"t = np.linspace(-1, 1, 50)\n",
"for axs, s in zip(ax, sigmas): \n",
" axs.scatter(t, noisy_abs(t, s))\n",
" axs.set_title(\"Noise: $\\sigma$={}\".format(s))"
" axs.set_title(r\"Noise: $\\sigma$={}\".format(s))"
]
},
{
Expand Down Expand Up @@ -482,7 +489,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"### Spectral Method"
"### Spectral method - Fourier basis"
]
},
{
Expand All @@ -497,12 +504,33 @@
"outputs": [],
"source": [
"no_filter = derivative.Spectral()\n",
"yes_filter = derivative.Spectral(filter=np.vectorize(lambda f: 1 if abs(f) < 1 else 0))\n",
"yes_filter = derivative.Spectral(filter=np.vectorize(lambda k: 1 if abs(k) < 6 else 0))\n",
"\n",
"plot_example(no_filter, t, noisy_abs, d_abs, sigmas, 'No filter')\n",
"plot_example(yes_filter, t, noisy_abs, d_abs, sigmas, 'Low-pass filter')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Spectral method - Chebyshev basis"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"t_cos = np.cos(np.pi * np.arange(50)/49)\n",
"no_filter = derivative.Spectral(basis='chebyshev')\n",
"yes_filter = derivative.Spectral(basis='chebyshev', filter=np.vectorize(lambda k: 1 if abs(k) < 15 else 0))\n",
"\n",
"plot_example(no_filter, t_cos, noisy_abs, d_abs, sigmas, 'No filter')\n",
"plot_example(yes_filter, t_cos, noisy_abs, d_abs, sigmas, 'Low-pass filter')"
]
},
{
"cell_type": "markdown",
"metadata": {},
Expand Down Expand Up @@ -549,18 +577,11 @@
"kal = derivative.Kalman(alpha=1.)\n",
"plot_example(kal, t, noisy_abs, d_abs, sigmas, 'alpha: 1.')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
Expand All @@ -574,7 +595,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.10"
"version": "3.13.1"
},
"latex_envs": {
"LaTeX_envs_menu_present": true,
Expand Down
Loading