From 9ff0d24cc424b156ac2e5c42f4623963b81278d4 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sun, 14 Jul 2024 15:16:46 +0300 Subject: [PATCH 01/26] feat: add default side to derivatives --- src/pycaputo/derivatives.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/pycaputo/derivatives.py b/src/pycaputo/derivatives.py index 40a2539..91f17d5 100644 --- a/src/pycaputo/derivatives.py +++ b/src/pycaputo/derivatives.py @@ -51,7 +51,7 @@ class RiemannLiouvilleDerivative(FractionalOperator): alpha: float """Order of the Riemann-Liouville derivative.""" - side: Side + side: Side = Side.Left """Side on which to compute the derivative.""" @property @@ -83,7 +83,7 @@ class CaputoDerivative(FractionalOperator): alpha: float """Order of the Caputo derivative.""" - side: Side + side: Side = Side.Left """Side on which to compute the derivative.""" @property @@ -114,7 +114,7 @@ class GrunwaldLetnikovDerivative(FractionalOperator): alpha: float """Order of the Grünwald-Letnikov derivative.""" - side: Side + side: Side = Side.Left """Side on which to compute the derivative.""" @property @@ -143,7 +143,7 @@ class HadamardDerivative(FractionalOperator): alpha: float """Order of the Hadamard derivative.""" - side: Side + side: Side = Side.Left """Side on which to compute the derivative.""" @property @@ -172,7 +172,7 @@ class CaputoHadamardDerivative(FractionalOperator): alpha: float """Order of the Caputo-Hadamard derivative.""" - side: Side + side: Side = Side.Left """Side on which to compute the derivative.""" @property From 6aef8fdf450355312ab3aa17acd5c411481a5faa Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sun, 14 Jul 2024 15:18:03 +0300 Subject: [PATCH 02/26] feat: define new diff interface --- src/pycaputo/differentiation/base.py | 54 +++++++++++++++++++++++++--- 1 file changed, 50 insertions(+), 4 deletions(-) diff --git a/src/pycaputo/differentiation/base.py b/src/pycaputo/differentiation/base.py index 370d7ad..399189b 100644 --- a/src/pycaputo/differentiation/base.py +++ b/src/pycaputo/differentiation/base.py @@ -19,7 +19,10 @@ class DerivativeMethod(ABC): @property def name(self) -> str: """An (non-unique) identifier for the differentiation method.""" - return type(self).__name__.replace("Method", "") + cls = self.__class__ + name = cls.__name__ + module = cls.__module__.split(".")[-1] + return f"diff_{module}_{name}".lower() @property @abstractmethod @@ -27,13 +30,56 @@ def d(self) -> FractionalOperator: """The fractional operator that is being discretized by the method.""" +@singledispatch +def quadrature_weights(m: DerivativeMethod, p: Array) -> Array: + r"""Evaluate the quadrature weights for the method *m* with points *p*. + + In general, a fractional operator is an integral operator that we can write + as + + .. math:: + + D^*[f](x_n) \approx \sum_{k = 0}^n w_{n, k} f_k + + where :math:`w_{n, k}` is a row of the weight matrix that is computed by this + function. In some cases of interest, e.g. uniform, it is sufficient to + compute the :math:`w_{N, k}` row and perform the convolution by FFT. + """ + raise NotImplementedError( + f"Cannot evaluate quadrature weights for method '{type(m).__name__}'" + ) + + +@singledispatch +def diffs(m: DerivativeMethod, f: ArrayOrScalarFunction, p: Points) -> Array: + r"""Evaluate the fractional derivative of *f* using method *m* at the point + *p.a* with the underlying grid. + + This function evaluates :math:`D^*[f](b)`, where :math:`f: [a, b] \to \mathbb{R}`, + i.e. at the end of the interval. This is in contract to :func:`diff`, which + evaluates the fractional derivative at all the points in the grid *p*. + + .. note:: + + The function *f* can be provided as a callable of as an array of values + evaluated on the grid *p*. However, not all methods support both variants + (e.g. if they require evaluating the function at additional points). + + :arg m: method used to evaluate the derivative. + :arg f: a simple function for which to evaluate the derivative. + :arg p: an array of points at which to evaluate the derivative. + """ + raise NotImplementedError( + f"Cannot evaluate fractional derivative with method '{type(m).__name__}'" + ) + + @singledispatch def diff(m: DerivativeMethod, f: ArrayOrScalarFunction, p: Points) -> Array: """Evaluate the fractional derivative of *f* at *p* using method *m*. - This is a low-level function that should be implemented by methods. For - a higher-level function see :func:`pycaputo.diff`. Note that this function - evaluates the derivative at all points in *p*, not just at *p[-1]*. + This function uses :func:`diffs` to evaluate the fractional derivative at + all the points in *p*. In most cases it is a trivial wrapper. :arg m: method used to evaluate the derivative. :arg f: a simple function for which to evaluate the derivative. From 7e96547964078a57bda2e18f66f8eee355758ffe Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sun, 14 Jul 2024 18:42:21 +0300 Subject: [PATCH 03/26] feat(diff): port Caputo L1 to new diff interface --- src/pycaputo/differentiation/__init__.py | 15 +++-- src/pycaputo/differentiation/base.py | 9 ++- src/pycaputo/differentiation/caputo.py | 76 ++++++++++++++++-------- 3 files changed, 68 insertions(+), 32 deletions(-) diff --git a/src/pycaputo/differentiation/__init__.py b/src/pycaputo/differentiation/__init__.py index 4991780..1ecfbc2 100644 --- a/src/pycaputo/differentiation/__init__.py +++ b/src/pycaputo/differentiation/__init__.py @@ -3,8 +3,13 @@ from __future__ import annotations -from pycaputo.derivatives import FractionalOperator, Side -from pycaputo.differentiation.base import DerivativeMethod, diff +from pycaputo.derivatives import FractionalOperator +from pycaputo.differentiation.base import ( + DerivativeMethod, + diff, + diffs, + quadrature_weights, +) from pycaputo.grid import Points @@ -22,7 +27,7 @@ def guess_method_for_order( :arg p: a set of points on which to evaluate the fractional operator. :arg d: a fractional operator to discretize. If only a float is given, the - common Caputo derivative is used. + common (left) Caputo derivative is used. """ from pycaputo import grid from pycaputo.derivatives import CaputoDerivative, RiemannLiouvilleDerivative @@ -31,7 +36,7 @@ def guess_method_for_order( m: DerivativeMethod | None = None if not isinstance(d, FractionalOperator): - d = CaputoDerivative(alpha=d, side=Side.Left) + d = CaputoDerivative(alpha=d) if isinstance(d, CaputoDerivative): if isinstance(p, grid.JacobiGaussLobattoPoints): @@ -61,5 +66,7 @@ def guess_method_for_order( __all__ = ( "DerivativeMethod", "diff", + "diffs", "guess_method_for_order", + "quadrature_weights", ) diff --git a/src/pycaputo/differentiation/base.py b/src/pycaputo/differentiation/base.py index 399189b..bcd0fc7 100644 --- a/src/pycaputo/differentiation/base.py +++ b/src/pycaputo/differentiation/base.py @@ -9,7 +9,7 @@ from pycaputo.derivatives import FractionalOperator from pycaputo.grid import Points -from pycaputo.typing import Array, ArrayOrScalarFunction +from pycaputo.typing import Array, ArrayOrScalarFunction, Scalar @dataclass(frozen=True) @@ -31,7 +31,7 @@ def d(self) -> FractionalOperator: @singledispatch -def quadrature_weights(m: DerivativeMethod, p: Array) -> Array: +def quadrature_weights(m: DerivativeMethod, p: Points, n: int) -> Array: r"""Evaluate the quadrature weights for the method *m* with points *p*. In general, a fractional operator is an integral operator that we can write @@ -44,6 +44,9 @@ def quadrature_weights(m: DerivativeMethod, p: Array) -> Array: where :math:`w_{n, k}` is a row of the weight matrix that is computed by this function. In some cases of interest, e.g. uniform, it is sufficient to compute the :math:`w_{N, k}` row and perform the convolution by FFT. + + :arg p: a grid on which to compute the quadrature weights at the point + ``p.x[n]``. """ raise NotImplementedError( f"Cannot evaluate quadrature weights for method '{type(m).__name__}'" @@ -51,7 +54,7 @@ def quadrature_weights(m: DerivativeMethod, p: Array) -> Array: @singledispatch -def diffs(m: DerivativeMethod, f: ArrayOrScalarFunction, p: Points) -> Array: +def diffs(m: DerivativeMethod, f: ArrayOrScalarFunction, p: Points, n: int) -> Scalar: r"""Evaluate the fractional derivative of *f* using method *m* at the point *p.a* with the underlying grid. diff --git a/src/pycaputo/differentiation/caputo.py b/src/pycaputo/differentiation/caputo.py index c08355c..0d80f2a 100644 --- a/src/pycaputo/differentiation/caputo.py +++ b/src/pycaputo/differentiation/caputo.py @@ -10,12 +10,17 @@ import numpy as np -from pycaputo.derivatives import CaputoDerivative, Side +from pycaputo.derivatives import CaputoDerivative from pycaputo.grid import Points, UniformMidpoints, UniformPoints from pycaputo.logging import get_logger -from pycaputo.typing import Array, ArrayOrScalarFunction, DifferentiableScalarFunction +from pycaputo.typing import ( + Array, + ArrayOrScalarFunction, + DifferentiableScalarFunction, + Scalar, +) -from .base import DerivativeMethod, diff +from .base import DerivativeMethod, diff, diffs, quadrature_weights logger = get_logger(__name__) @@ -35,7 +40,7 @@ def __post_init__(self) -> None: @property def d(self) -> CaputoDerivative: - return CaputoDerivative(self.alpha, side=Side.Left) + return CaputoDerivative(self.alpha) # {{{ L1 @@ -60,37 +65,58 @@ def __post_init__(self) -> None: ) -def _weights_l1(m: L1, p: Points) -> Iterator[Array]: - x, dx = p.x, p.dx +@quadrature_weights.register(L1) +def _quadrature_weights_caputo_l1(m: L1, p: Points, n: int) -> Array: + if n < 0: + n = p.size + n + + if not 0 <= n <= p.size: + raise IndexError(f"Index 'n' out of range: 0 <= {n} < {p.size}") + alpha = m.alpha - w0 = 1 / math.gamma(2 - alpha) + x, dx = p.x[:n], p.dx[: n - 1] + xn = x[-1] + + w = np.empty_like(x) + w[n - 1] = 0.0 + w[: n - 1] = ( + ((xn - x[:-1]) ** (1 - alpha) - (xn - x[1:]) ** (1 - alpha)) + / dx + / math.gamma(2 - alpha) + ) + w[1:n] = w[: n - 1] - w[1:n] + w[0] = -w[0] - # NOTE: weights given in [Li2020] Equation 4.20 - if isinstance(p, UniformPoints): - w0 = w0 / p.dx[0] ** alpha - k = np.arange(x.size - 1) + return w - for n in range(1, x.size): - w = (n - k[:n]) ** (1 - alpha) - (n - k[:n] - 1) ** (1 - alpha) - yield w0 * w - else: - for n in range(1, x.size): - w = ( - (x[n] - x[:n]) ** (1 - alpha) - (x[n] - x[1 : n + 1]) ** (1 - alpha) - ) / dx[:n] - yield w0 * w +@diffs.register(L1) +def _diffs_caputo_l1(m: L1, f: ArrayOrScalarFunction, p: Points, n: int) -> Scalar: + if n < 0: + n = p.size + n + + if not 0 <= n <= p.size: + raise IndexError(f"Index 'n' out of range: 0 <= {n} < {p.size}") + + if n == 1: + return np.array([np.nan]) + + w = quadrature_weights(m, p, n + 1) + fx = f(p.x[: n + 1]) if callable(f) else f + + return np.sum(w * fx) # type: ignore[no-any-return] @diff.register(L1) -def _diff_l1_method(m: L1, f: ArrayOrScalarFunction, p: Points) -> Array: - dfx = np.diff(f(p.x) if callable(f) else f) +def _diff_caputo_l1(m: L1, f: ArrayOrScalarFunction, p: Points) -> Array: + fx = f(p.x) if callable(f) else f - df = np.empty(p.x.shape, dtype=dfx.dtype) + df = np.empty_like(fx) df[0] = np.nan - for n, w in enumerate(_weights_l1(m, p)): - df[n + 1] = np.sum(w * dfx[: n + 1]) + for n in range(1, df.size): + w = quadrature_weights(m, p, n + 1) + df[n] = np.sum(w * fx[: n + 1]) return df From c4dfa64d38b1899b497b364193c4ea5af20c8ee1 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sun, 14 Jul 2024 19:33:11 +0300 Subject: [PATCH 04/26] feat: add example for constructing matrix --- examples/caputo-derivative-l1-spectrum.py | 85 +++++++++++++++++++++++ 1 file changed, 85 insertions(+) create mode 100644 examples/caputo-derivative-l1-spectrum.py diff --git a/examples/caputo-derivative-l1-spectrum.py b/examples/caputo-derivative-l1-spectrum.py new file mode 100644 index 0000000..dec8d22 --- /dev/null +++ b/examples/caputo-derivative-l1-spectrum.py @@ -0,0 +1,85 @@ +# SPDX-FileCopyrightText: 2023 Alexandru Fikl +# SPDX-License-Identifier: MIT + +r"""Computes the matrix operator for the Caputo derivative and look at its +spectrum for various values of :math:`\alpha`. + +The spectrum of the matrix is quite boring, since the matrix itself is lower +triangular, so all the eigenvalues are on the diagonal. From the method, we know +that + +.. math:: + + W_{n, n} = \frac{1}{\Gamma(2 - \alpha)} \frac{1}{\Delta x^\alpha}. + +i.e. the diagonal is constant. This mostly serves as an example of how the +matrices can be constructed. +""" + +from __future__ import annotations + +import numpy as np + +from pycaputo.differentiation import quadrature_weights +from pycaputo.differentiation.caputo import L1 +from pycaputo.grid import make_uniform_points +from pycaputo.logging import get_logger +from pycaputo.utils import Array + +logger = get_logger("caputo_derivative_l1_spectrum") + + +# {{{ matrix + + +def get_l1_matrix(alpha: float, *, n: int, a: float, b: float) -> Array: + W = np.zeros((n, n)) + + meth = L1(alpha=alpha) + points = make_uniform_points(n, a=0, b=1) + + for i in range(1, n): + W[i, : i + 1] = quadrature_weights(meth, points, i + 1) + W[0, 0] = W[1, 1] + + return W + + +# }}} + + +# {{{ conditioning + +n = 256 +alphas = [0.1, 0.25, 0.5, 0.75, 0.9, 0.95] + +for alpha in alphas: + W = get_l1_matrix(alpha, n=n, a=0.0, b=1.0) + kappa = np.linalg.cond(W) + logger.info("alpha = %.2f kappa = %.8e", alpha, kappa) + +# }}} + +# {{{ plot + +try: + import matplotlib # noqa: F401 +except ImportError as exc: + raise SystemExit(0) from exc + +from pycaputo.utils import figure, set_recommended_matplotlib + +set_recommended_matplotlib() + +with figure("caputo-derivative-l1-spectrum") as fig: + ax = fig.gca() + + for alpha in alphas: + W = get_l1_matrix(alpha, n=n, a=0.0, b=1.0) + sigma = np.linalg.eigvals(W) + ax.plot(sigma.real, sigma.imag, "o") + + ax.set_xlabel(r"$\Re$") + ax.set_ylabel(r"$\Im$") + +# }}} From bef7801153e8f539b7e586d0fb4817ae428dbe71 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sun, 14 Jul 2024 20:49:53 +0300 Subject: [PATCH 05/26] feat(grid): add helper to create a midpoint grid --- src/pycaputo/grid.py | 24 ++++++++++++++++++------ 1 file changed, 18 insertions(+), 6 deletions(-) diff --git a/src/pycaputo/grid.py b/src/pycaputo/grid.py index 7994e01..d6ffd8c 100644 --- a/src/pycaputo/grid.py +++ b/src/pycaputo/grid.py @@ -148,23 +148,35 @@ def make_uniform_points(n: int, a: float = 0.0, b: float = 1.0) -> UniformPoints @dataclass(frozen=True) -class UniformMidpoints(Points): - """A set of points on :math:`[a, b]` formed from midpoints of a uniform grid. +class MidPoints(Points): + """A set of points on :math:`[a, b]` formed from midpoints of another grid. This set of points consists of :math:`x_0 = a` and :math:`x_j` are the - midpoints of the uniform grid :class:`UniformPoints` + midpoints of another grid. """ -def make_uniform_midpoints(n: int, a: float = 0.0, b: float = 1.0) -> UniformMidpoints: - """Construct a :class:`UniformMidpoints` on :math:`[a, b]`. +def make_uniform_midpoints(n: int, a: float = 0.0, b: float = 1.0) -> MidPoints: + """Construct a :class:`MidPoints` on :math:`[a, b]`. :arg n: number of points in :math:`[a, b]`. """ x = np.linspace(a, b, n) x[1:] = (x[1:] + x[:-1]) / 2 - return UniformMidpoints(a=a, b=b, x=x) + return MidPoints(a=a, b=b, x=x) + + +def make_midpoints_from(p: Points) -> MidPoints: + """Construct a grid of midpoints from a given set of points *p*. + + The new grid will contain the point :math:`x_0` and the midpoints of *p*. + """ + x = np.empty_like(p.x) + x[0] = p.x[0] + x[1:] = p.xm + + return MidPoints(a=p.a, b=p.b, x=x) # }}} From c263236919febdf415f48e059a588e1ff7a793af Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sun, 14 Jul 2024 20:50:43 +0300 Subject: [PATCH 06/26] feat(diff): adapt ModifiedL1 to new interface --- src/pycaputo/differentiation/__init__.py | 2 +- src/pycaputo/differentiation/caputo.py | 83 ++++++++++++++---------- 2 files changed, 48 insertions(+), 37 deletions(-) diff --git a/src/pycaputo/differentiation/__init__.py b/src/pycaputo/differentiation/__init__.py index 1ecfbc2..c2e6ee3 100644 --- a/src/pycaputo/differentiation/__init__.py +++ b/src/pycaputo/differentiation/__init__.py @@ -42,7 +42,7 @@ def guess_method_for_order( if isinstance(p, grid.JacobiGaussLobattoPoints): m = caputo.SpectralJacobi(d.alpha) elif 0 < d.alpha < 1: - if isinstance(p, grid.UniformMidpoints): + if isinstance(p, grid.MidPoints): m = caputo.ModifiedL1(d.alpha) else: m = caputo.L1(d.alpha) diff --git a/src/pycaputo/differentiation/caputo.py b/src/pycaputo/differentiation/caputo.py index 0d80f2a..4bdfb62 100644 --- a/src/pycaputo/differentiation/caputo.py +++ b/src/pycaputo/differentiation/caputo.py @@ -6,12 +6,11 @@ import math from abc import abstractmethod from dataclasses import dataclass -from typing import Iterator import numpy as np from pycaputo.derivatives import CaputoDerivative -from pycaputo.grid import Points, UniformMidpoints, UniformPoints +from pycaputo.grid import MidPoints, Points, UniformPoints, make_midpoints_from from pycaputo.logging import get_logger from pycaputo.typing import ( Array, @@ -109,11 +108,22 @@ def _diffs_caputo_l1(m: L1, f: ArrayOrScalarFunction, p: Points, n: int) -> Scal @diff.register(L1) def _diff_caputo_l1(m: L1, f: ArrayOrScalarFunction, p: Points) -> Array: - fx = f(p.x) if callable(f) else f + if callable(f): + fx = f(p.x) + elif isinstance(f, np.ndarray): + if f.shape[0] != p.size: + raise ValueError( + f"Shape of 'f' does match points: got {f.shape} expected{p.shape}" + ) + fx = f + else: + raise TypeError(f"Unknown type for 'f': {type(f)}") df = np.empty_like(fx) df[0] = np.nan + # FIXME: in the uniform case, we can also do an FFT, but we need different + # weights for that, so we leave it like this for now for n in range(1, df.size): w = quadrature_weights(m, p, n + 1) df[n] = np.sum(w * fx[: n + 1]) @@ -132,8 +142,9 @@ class ModifiedL1(CaputoMethod): r"""Implements the modified L1 method for the Caputo fractional derivative of order :math:`\alpha \in (0, 1)`. - This method is defined in Section 4.1.1 (III) from [Li2020]_ for quasi-uniform - grids. These grids can be constructed by :class:`~pycaputo.grid.UniformMidpoints`. + This method is defined in Section 4.1.1 (III) from [Li2020]_. Note that this + method evaluates the fractional derivative at the midpoints + :math:`(x_i + x_{i - 1}) / 2` for any given input grid. """ if __debug__: @@ -146,48 +157,48 @@ def __post_init__(self) -> None: ) -def _weights_modified_l1(m: ModifiedL1, p: Points) -> Iterator[Array]: - if not isinstance(p, UniformMidpoints): - raise NotImplementedError( - f"'{type(m).__name__}' does not implement 'weights' for" - f" '{type(p).__name__}' grids" - ) +@quadrature_weights.register(ModifiedL1) +def _quadrature_weights_caputo_modified_l1(m: ModifiedL1, p: Points, n: int) -> Array: + if isinstance(p, MidPoints): + return quadrature_weights.dispatch(L1)(m, p, n) + else: + p = make_midpoints_from(p) + w = quadrature_weights.dispatch(L1)(m, p, n) - # NOTE: weights from [Li2020] Equation 4.51 - # FIXME: this does not use the formula from the book; any benefit to it? - alpha = m.alpha - wc = 1 / math.gamma(2 - alpha) / p.dx[-1] ** alpha - k = np.arange(p.x.size) + # FIXME: not clear if this does anything? Might be the same as just + # doing L1 on the original grid.. + wp = np.empty((n + 1,), dtype=w.dtype) + wp[0] = w[0] / 2.0 + wp[-1] = w[-1] / 2.0 + wp[1:-1] = (w[1:] + w[:-1]) / 2.0 + + return wp - # NOTE: first interval has a size of h / 2 and is weighted differently - w0 = 2 * ((k[:-1] + 0.5) ** (1 - alpha) - k[:-1] ** (1 - alpha)) - for n in range(1, p.x.size): - w = (n - k[:n]) ** (1 - alpha) - (n - k[:n] - 1) ** (1 - alpha) - w[0] = w0[n - 1] +@diffs.register(ModifiedL1) +def _diffs_caputo_modified_l1( + m: ModifiedL1, f: ArrayOrScalarFunction, p: Points, n: int +) -> Scalar: + if not isinstance(p, MidPoints): + p = make_midpoints_from(p) - yield wc * w + if not callable(f): + f[1:] = (f[:-1] + f[1]) / 2.0 + + return diffs.dispatch(L1)(m, f, p, n) @diff.register(ModifiedL1) -def _diff_modified_l1_method( +def _diff_caputo_modified_l1( m: ModifiedL1, f: ArrayOrScalarFunction, p: Points ) -> Array: - if not isinstance(p, UniformMidpoints): - raise NotImplementedError( - f"'{type(m).__name__}' does not implement 'diff' for '{type(p).__name__}'" - " grids" - ) + if not isinstance(p, MidPoints): + p = make_midpoints_from(p) - dfx = np.diff(f(p.x) if callable(f) else f) + if not callable(f): + f[1:] = (f[:-1] + f[1]) / 2.0 - df = np.empty(p.x.size) - df[0] = np.nan - - for n, w in enumerate(_weights_modified_l1(m, p)): - df[n + 1] = np.sum(w * dfx[: n + 1]) - - return df + return diff.dispatch(L1)(m, f, p) # }}} From b0396b7b5f77b835246e81021e7431fedaa9f809 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Thu, 18 Jul 2024 15:32:04 +0300 Subject: [PATCH 07/26] feat(diff): adapt L2 to new interface --- src/pycaputo/differentiation/caputo.py | 141 +++++++++++++++++++++---- tests/test_diff_caputo.py | 6 ++ 2 files changed, 124 insertions(+), 23 deletions(-) diff --git a/src/pycaputo/differentiation/caputo.py b/src/pycaputo/differentiation/caputo.py index 4bdfb62..f8f5d9f 100644 --- a/src/pycaputo/differentiation/caputo.py +++ b/src/pycaputo/differentiation/caputo.py @@ -64,25 +64,37 @@ def __post_init__(self) -> None: ) -@quadrature_weights.register(L1) -def _quadrature_weights_caputo_l1(m: L1, p: Points, n: int) -> Array: - if n < 0: - n = p.size + n +def _caputo_piecewise_constant_integral(p: Points, n: int, alpha: float) -> Array: + r"""Computes the integrals - if not 0 <= n <= p.size: - raise IndexError(f"Index 'n' out of range: 0 <= {n} < {p.size}") + .. math:: + + a_{nk} = \int_{x_k}^{x_{k + 1}} \frac{1}{(x_n - s)^\alpha} ds + + for :math:`0 \le k \le n - 1`. + """ - alpha = m.alpha x, dx = p.x[:n], p.dx[: n - 1] xn = x[-1] - w = np.empty_like(x) - w[n - 1] = 0.0 - w[: n - 1] = ( + return ( ((xn - x[:-1]) ** (1 - alpha) - (xn - x[1:]) ** (1 - alpha)) / dx / math.gamma(2 - alpha) ) + + +@quadrature_weights.register(L1) +def _quadrature_weights_caputo_l1(m: L1, p: Points, n: int) -> Array: + if n < 0: + n = p.size + n + + if not 0 <= n <= p.size: + raise IndexError(f"Index 'n' out of range: 0 <= {n} < {p.size}") + + w = np.empty(n, dtype=p.x.dtype) + w[n - 1] = 0.0 + w[: n - 1] = _caputo_piecewise_constant_integral(p, n, m.alpha) w[1:n] = w[: n - 1] - w[1:n] w[0] = -w[0] @@ -97,27 +109,22 @@ def _diffs_caputo_l1(m: L1, f: ArrayOrScalarFunction, p: Points, n: int) -> Scal if not 0 <= n <= p.size: raise IndexError(f"Index 'n' out of range: 0 <= {n} < {p.size}") - if n == 1: + if n == 0: return np.array([np.nan]) w = quadrature_weights(m, p, n + 1) - fx = f(p.x[: n + 1]) if callable(f) else f + fx = f(p.x[: n + 1]) if callable(f) else f[:n + 1] return np.sum(w * fx) # type: ignore[no-any-return] @diff.register(L1) def _diff_caputo_l1(m: L1, f: ArrayOrScalarFunction, p: Points) -> Array: - if callable(f): - fx = f(p.x) - elif isinstance(f, np.ndarray): - if f.shape[0] != p.size: - raise ValueError( - f"Shape of 'f' does match points: got {f.shape} expected{p.shape}" - ) - fx = f - else: - raise TypeError(f"Unknown type for 'f': {type(f)}") + fx = f(p.x) if callable(f) else f + if fx.shape[0] != p.size: + raise ValueError( + f"Shape of 'f' does match points: got {f.shape} expected {p.shape}" + ) df = np.empty_like(fx) df[0] = np.nan @@ -145,6 +152,10 @@ class ModifiedL1(CaputoMethod): This method is defined in Section 4.1.1 (III) from [Li2020]_. Note that this method evaluates the fractional derivative at the midpoints :math:`(x_i + x_{i - 1}) / 2` for any given input grid. + + As noted in [Li2020]_, this method has the same order of convergence as the + standard L1 method. However the weights can be used to construct a + Crank-Nicolson type method. """ if __debug__: @@ -225,6 +236,90 @@ def __post_init__(self) -> None: ) +def _caputo_d2_coefficients(x: Array) -> Array: + x0, x1, x2, x3 = x[:4] + a0 = 2.0 * (3.0 * x0 - x1 - x2 - x3) / ((x0 - x1) * (x0 - x2) * (x0 - x3)) + b0 = 2.0 * (x3 - 2.0 * x0 + x2) / ((x0 - x1) * (x1 - x2) * (x1 - x3)) + c0 = 2.0 * (x3 - 2.0 * x0 + x1) / ((x0 - x2) * (x2 - x1) * (x2 - x3)) + d0 = 2.0 * (x2 - 2.0 * x0 + x1) / ((x0 - x3) * (x3 - x1) * (x3 - x2)) + + return np.array([a0, b0, c0, d0]) + + +@quadrature_weights.register(L2) +def _quadrature_weights_caputo_l2(m: L2, p: Points, n: int) -> Array: + if n < 0: + n = p.size + n + + if not 0 <= n <= p.size: + raise IndexError(f"Index 'n' out of range: 0 <= {n} < {p.size}") + + dx = p.dx[:n - 1] + dxm = (dx[1:] + dx[:-1]) / 2.0 + + # get finite difference coefficients for center stencil + a = 1.0 / (dx[:-2] * dxm) + b = 1.0 / (dx[1:-1] * dxm) + + # get finite difference coefficients for boundary stencils + a_l, b_l, c_l, d_l = _caputo_d2_coefficients(p.x) + d_r, c_r, b_r, a_r = _caputo_d2_coefficients(p.x[::-1]) + + wi = _caputo_piecewise_constant_integral(p, n, m.alpha) + + w = np.empty_like(p.x[:n]) + # center stencil + w[2:n - 1] = a * wi[2:] - (a + b) * w[1:-1] + b * w[:-2] + # add left boundary stencil + w[0] = a_l + a[1] * wi[1] + w[1] = b_l + a[2] * wi[2] - (a[1] + b[1]) * wi[1] + w[2] += c_l + w[3] += d_l + # add right-boundary stencil + if n == p.size: + w[-4] += d_r + w[-3] += c_r + w[-2] = b_r + a[-2] * wi[-2] - (a[-1] + b[-1]) * wi[-1] + w[-1] = a_r + a[-1] * wi[-1] + + return w + + +@diffs.register(L2) +def _diffs_caputo_l2(m: L2, f: ArrayOrScalarFunction, p: Points, n: int) -> Array: + if n < 0: + n = p.size + n + + if not 0 <= n <= p.size: + raise IndexError(f"Index 'n' out of range: 0 <= {n} < {p.size}") + + if n == 0: + return np.array([np.nan]) + + w = quadrature_weights(m, p, n + 1) + fx = f(p.x[: n + 1]) if callable(f) else f[:n + 1] + + return np.sum(w * fx) # type: ignore[no-any-return] + + +@diff.register(L2) +def _diff_caputo_l2(m: L2, f: ArrayOrScalarFunction, p: Points) -> Array: + fx = f(p.x) if callable(f) else f + if fx.shape[0] != p.size: + raise ValueError( + f"Shape of 'f' does match points: got {f.shape} expected {p.shape}" + ) + + df = np.empty_like(fx) + df[0] = np.nan + + # FIXME: in the uniform case, we can also do an FFT, but we need different + # weights for that, so we leave it like this for now + for n in range(1, df.size): + w = quadrature_weights(m, p, n + 1) + df[n] = np.sum(w * fx[: n + 1]) + + def _weights_l2(alpha: float, i: int | Array, k: int | Array) -> Array: return np.array((i - k) ** (2 - alpha) - (i - k - 1) ** (2 - alpha)) @@ -377,7 +472,7 @@ def _diff_jacobi(m: SpectralJacobi, f: ArrayOrScalarFunction, p: Points) -> Arra # }}} -# {{{ +# {{{ DiffusiveCaputoMethod @dataclass(frozen=True) diff --git a/tests/test_diff_caputo.py b/tests/test_diff_caputo.py index 744fb61..f6983d3 100644 --- a/tests/test_diff_caputo.py +++ b/tests/test_diff_caputo.py @@ -66,6 +66,7 @@ def df_test(x: Array, *, alpha: float, mu: float = 3.5) -> Array: ("L2", "uniform"), ("L2C", "uniform"), ("ModifiedL1", "midpoints"), + ("ModifiedL1", "stretch"), ], ) @pytest.mark.parametrize("alpha", [0.1, 0.25, 0.5, 0.75, 0.9]) @@ -115,6 +116,11 @@ def test_caputo_lmethods( for n in [16, 32, 64, 128, 256, 512, 768, 1024]: p = make_points_from_name(grid_type, n, a=0.0, b=0.5) + if name == "ModifiedL1" and grid_type != "midpoints": + from pycaputo.grid import make_midpoints_from + + p = make_midpoints_from(p) + df_num = diff(meth, f_test, p) df_ref = df_test(p.x, alpha=alpha) From 8b214d20af69a4d1f0c6371a280d637d38bffec9 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Wed, 7 Aug 2024 17:11:58 +0300 Subject: [PATCH 08/26] docs: improve docs for new diff functionality --- docs/guide_differentiation.rst | 35 ++++++++-- examples/guide-differentiation.py | 23 ++++++- src/pycaputo/__init__.py | 18 +++--- src/pycaputo/differentiation/base.py | 90 +++++++++++++++++++++++--- src/pycaputo/differentiation/caputo.py | 16 +++-- 5 files changed, 148 insertions(+), 34 deletions(-) diff --git a/docs/guide_differentiation.rst b/docs/guide_differentiation.rst index 70ce97e..ba55917 100644 --- a/docs/guide_differentiation.rst +++ b/docs/guide_differentiation.rst @@ -1,8 +1,8 @@ A New Differentiation Method ============================ -Computing the derivative of a function can be done in two ways. The recommended -method (with included magic) is calling :func:`pycaputo.diff` as +Computing the fractional derivative of a function can be done in two ways. The +recommended method (with included magic) is calling :func:`pycaputo.diff` as .. code:: python @@ -10,8 +10,8 @@ method (with included magic) is calling :func:`pycaputo.diff` as which will automatically select an appropriate method to use given the point set ``p`` and the order ``alpha`` (see also -:func:`pycaputo.differentiation.guess_method_for_order`). To manually call a -specific method, use :func:`pycaputo.differentiation.diff` instead as +:func:`~pycaputo.differentiation.guess_method_for_order`). To manually call a +specific method, use the lower level :func:`pycaputo.differentiation.diff` instead as .. code:: python @@ -40,8 +40,31 @@ abstract methods of the base class. For example, :start-after: [class-definition-start] :end-before: [class-definition-end] -Then, we can implement the :func:`~pycaputo.differentiation.diff` method by -registering it with the :func:`~functools.singledispatch` mechanism as +Then, we have defined a set of functions useful for working with a discrete +fractional operator approximations. As mentioned above, the main functionality +is provided by the :func:`~pycaputo.differentiation.diff` method. However, we +also have + +* :func:`~pycaputo.differentiation.quadrature_weights`: a function that can provide + the underlying quadrature weights for a given method. Note that not all methods + can be efficiently expressed as a weighted sum with quadrature weights, so they + should not implement this method. However, for those that can, this offers a + pleasant way of reusing and analyzing the weights. + +* :func:`~pycaputo.differentiation.diffs`: a method very similar to + :func:`~pycaputo.differentiation.diff`, but that only computes the fractional + operator at a given point on the grid. Note that not all methods can efficiently + evaluate the fractional operator at a single point (e.g. global spectral methods), + so they should not implement this functionality. + +* :func:`~pycaputo.differentiation.diff`: for many methods, the ``diff`` function + is a simple wrapper around :func:`~pycaputo.differentiation.diffs`, but this + doesn't always make sense. For example, if we want to compute the Caputo + fractional derivative on a uniform grid, this can be evaluated more efficiently + using the FFT. + +These methods are all registered with the :func:`~functools.singledispatch` mechanism. +An example setup is provided below. .. literalinclude:: ../examples/guide-differentiation.py :language: python diff --git a/examples/guide-differentiation.py b/examples/guide-differentiation.py index da86e0e..e53608d 100644 --- a/examples/guide-differentiation.py +++ b/examples/guide-differentiation.py @@ -42,7 +42,28 @@ def d(self) -> RiemannLiouvilleDerivative: # {{{ # [register-start] -from pycaputo.differentiation import diff +from pycaputo.differentiation import diff, diffs, quadrature_weights + + +@quadrature_weights.register(RiemannLiouvilleDerivativeMethod) +def _quadrature_weights_rl( + m: RiemannLiouvilleDerivativeMethod, + p: Points, + n: int, +) -> Array: + # ... add an actual implementation here ... + return np.zeros(n, dtype=p.x.dtype) + + +@diffs.register(RiemannLiouvilleDerivativeMethod) +def _diffs_rl( + m: RiemannLiouvilleDerivativeMethod, + f: ArrayOrScalarFunction, + p: Points, +) -> Array: + fx = f(p.x) if callable(f) else f + # ... add an actual implementation here ... + return np.zeros_like(fx) @diff.register(RiemannLiouvilleDerivativeMethod) diff --git a/src/pycaputo/__init__.py b/src/pycaputo/__init__.py index b43e5e3..0943086 100644 --- a/src/pycaputo/__init__.py +++ b/src/pycaputo/__init__.py @@ -29,8 +29,7 @@ def diff( :arg d: a fractional operator. If this is just a number, the standard Caputo derivative will be used. """ - from pycaputo.differentiation import diff as _diff - from pycaputo.differentiation import guess_method_for_order + import pycaputo.differentiation as fracd if d is None: raise ValueError("'d' is required if 'method is not given") @@ -38,9 +37,9 @@ def diff( if not isinstance(d, FractionalOperator): d = CaputoDerivative(d, side=Side.Left) - m = guess_method_for_order(p, d) + m = fracd.guess_method_for_order(p, d) - return _diff(m, f, p) + return fracd.diff(m, f, p) def quad( @@ -56,8 +55,7 @@ def quad( :arg d: a fractional operator. If this is just a number, the standard Riemann-Liouville integral will be used. """ - from pycaputo.quadrature import guess_method_for_order - from pycaputo.quadrature import quad as _quad + import pycaputo.quadrature as fracq if d is None: raise ValueError("'d' is required if 'method' is not given") @@ -65,9 +63,9 @@ def quad( if not isinstance(d, FractionalOperator): d = RiemannLiouvilleDerivative(d, side=Side.Left) - m = guess_method_for_order(p, d) + m = fracq.guess_method_for_order(p, d) - return _quad(m, f, p) + return fracq.quad(m, f, p) def grad( @@ -118,12 +116,12 @@ def f_i(y: Array) -> Array: def make_component_p(i: tuple[int, ...]) -> Points: return p.translate(a[i], x[i]) - from pycaputo.differentiation import diff as _diff + import pycaputo.differentiation as fracd result = np.empty_like(x) for i in np.ndindex(x.shape): # FIXME: this should just compute the gradient at -1 - result[i] = _diff(m, make_component_f(i), make_component_p(i))[-1] + result[i] = fracd.diff(m, make_component_f(i), make_component_p(i))[-1] return result diff --git a/src/pycaputo/differentiation/base.py b/src/pycaputo/differentiation/base.py index bcd0fc7..5256d8b 100644 --- a/src/pycaputo/differentiation/base.py +++ b/src/pycaputo/differentiation/base.py @@ -7,6 +7,8 @@ from dataclasses import dataclass from functools import singledispatch +import numpy as np + from pycaputo.derivatives import FractionalOperator from pycaputo.grid import Points from pycaputo.typing import Array, ArrayOrScalarFunction, Scalar @@ -45,22 +47,29 @@ def quadrature_weights(m: DerivativeMethod, p: Points, n: int) -> Array: function. In some cases of interest, e.g. uniform, it is sufficient to compute the :math:`w_{N, k}` row and perform the convolution by FFT. - :arg p: a grid on which to compute the quadrature weights at the point - ``p.x[n]``. + .. note:: + + Not all methods can compute the weights for a single evaluation in + this fashion. This function is mainly implemented for methods where the + *nth* row can be computed efficiently. + + :arg p: a grid on which to compute the quadrature weights for the point ``p.x[n]``. + :arg n: the index of the point within *p*. """ raise NotImplementedError( - f"Cannot evaluate quadrature weights for method '{type(m).__name__}'" + f"Cannot evaluate quadrature weights for method '{type(m).__name__}' " + "('quadrature_weights' is not implemented)" ) @singledispatch def diffs(m: DerivativeMethod, f: ArrayOrScalarFunction, p: Points, n: int) -> Scalar: r"""Evaluate the fractional derivative of *f* using method *m* at the point - *p.a* with the underlying grid. + *p.x[n]* with the underlying grid. - This function evaluates :math:`D^*[f](b)`, where :math:`f: [a, b] \to \mathbb{R}`, - i.e. at the end of the interval. This is in contract to :func:`diff`, which - evaluates the fractional derivative at all the points in the grid *p*. + This function evaluates :math:`D^*[f](x_n)`, where :math:`f: [a, b] \to \mathbb{R}`. + This is in contract to :func:`diff`, which evaluates the fractional derivative + at all the points in the grid *p*. .. note:: @@ -72,10 +81,28 @@ def diffs(m: DerivativeMethod, f: ArrayOrScalarFunction, p: Points, n: int) -> S :arg f: a simple function for which to evaluate the derivative. :arg p: an array of points at which to evaluate the derivative. """ - raise NotImplementedError( - f"Cannot evaluate fractional derivative with method '{type(m).__name__}'" + try: + result = diff(m, f, p) + except NotImplementedError: + raise NotImplementedError( + "Cannot evaluate pointwise fractional derivative with method " + f"'{type(m).__name__}' ('diffs' is not implemented)" + ) from None + + from warnings import warn + + warn( + f"'diffs' is not implemented for '{type(m).__name__}' (aka '{m.name}'). " + "Falling back to 'diff', which is likely significantly slower! Use " + "'diff' directly if this is acceptable.", + stacklevel=2, ) + if not 0 <= n <= p.size: + raise IndexError(f"Index 'n' out of range: 0 <= {n} < {p.size}") + + return np.array(result[n]) + @singledispatch def diff(m: DerivativeMethod, f: ArrayOrScalarFunction, p: Points) -> Array: @@ -84,10 +111,53 @@ def diff(m: DerivativeMethod, f: ArrayOrScalarFunction, p: Points) -> Array: This function uses :func:`diffs` to evaluate the fractional derivative at all the points in *p*. In most cases it is a trivial wrapper. + .. note:: + + All methods are expected to implement this method, while + :func:`quadrature_weights` and :func:`diffs` are optional and not + supported by all methods. + :arg m: method used to evaluate the derivative. :arg f: a simple function for which to evaluate the derivative. :arg p: an array of points at which to evaluate the derivative. """ raise NotImplementedError( - f"Cannot evaluate fractional derivative with method '{type(m).__name__}'" + f"Cannot evaluate fractional derivative with method '{type(m).__name__}' " + "('diff' is not implemented)" ) + + +def quadrature_matrix(m: DerivativeMethod, p: Points, w00: float = 0.0) -> Array: + """Evaluate the quadrature matrix for a given method *m* at points *p*. + + This requires that the method implements the :func:`quadrature_weights`. As + fractional operator are usually convolutional, this results in a lower + triangular matrix. + + The value of the derivative at ``p.a`` is not defined. The value of the + matrix at ``W[0, 0]`` can be controlled using *w00*, which is set to 0 by + default. This default choice results in a singular operator and may not + always be desired. + + Using this matrix operator, we can evaluate the derivative equivalently as + + .. code:: python + + mat = quadrature_matrix(m, p) + + df_mat = mat @ f + df_fun = diff(m, f, p) + assert np.allclose(df_mat, df_fun) + + :returns: a two-dimensional array of shape ``(n, n)``, where *n* is the + number of points in *p*. + """ + + n = p.size + W = np.zeros((n, n)) + W[0, 0] = w00 + + for i in range(1, n): + W[i, : i + 1] = quadrature_weights(m, p, i + 1) + + return W diff --git a/src/pycaputo/differentiation/caputo.py b/src/pycaputo/differentiation/caputo.py index f8f5d9f..34c409a 100644 --- a/src/pycaputo/differentiation/caputo.py +++ b/src/pycaputo/differentiation/caputo.py @@ -77,7 +77,7 @@ def _caputo_piecewise_constant_integral(p: Points, n: int, alpha: float) -> Arra x, dx = p.x[:n], p.dx[: n - 1] xn = x[-1] - return ( + return np.array( ((xn - x[:-1]) ** (1 - alpha) - (xn - x[1:]) ** (1 - alpha)) / dx / math.gamma(2 - alpha) @@ -113,7 +113,7 @@ def _diffs_caputo_l1(m: L1, f: ArrayOrScalarFunction, p: Points, n: int) -> Scal return np.array([np.nan]) w = quadrature_weights(m, p, n + 1) - fx = f(p.x[: n + 1]) if callable(f) else f[:n + 1] + fx = f(p.x[: n + 1]) if callable(f) else f[: n + 1] return np.sum(w * fx) # type: ignore[no-any-return] @@ -123,7 +123,7 @@ def _diff_caputo_l1(m: L1, f: ArrayOrScalarFunction, p: Points) -> Array: fx = f(p.x) if callable(f) else f if fx.shape[0] != p.size: raise ValueError( - f"Shape of 'f' does match points: got {f.shape} expected {p.shape}" + f"Shape of 'f' does match points: got {fx.shape} expected {p.shape}" ) df = np.empty_like(fx) @@ -254,7 +254,7 @@ def _quadrature_weights_caputo_l2(m: L2, p: Points, n: int) -> Array: if not 0 <= n <= p.size: raise IndexError(f"Index 'n' out of range: 0 <= {n} < {p.size}") - dx = p.dx[:n - 1] + dx = p.dx[: n - 1] dxm = (dx[1:] + dx[:-1]) / 2.0 # get finite difference coefficients for center stencil @@ -269,7 +269,7 @@ def _quadrature_weights_caputo_l2(m: L2, p: Points, n: int) -> Array: w = np.empty_like(p.x[:n]) # center stencil - w[2:n - 1] = a * wi[2:] - (a + b) * w[1:-1] + b * w[:-2] + w[2 : n - 1] = a * wi[2:] - (a + b) * w[1:-1] + b * w[:-2] # add left boundary stencil w[0] = a_l + a[1] * wi[1] w[1] = b_l + a[2] * wi[2] - (a[1] + b[1]) * wi[1] @@ -297,7 +297,7 @@ def _diffs_caputo_l2(m: L2, f: ArrayOrScalarFunction, p: Points, n: int) -> Arra return np.array([np.nan]) w = quadrature_weights(m, p, n + 1) - fx = f(p.x[: n + 1]) if callable(f) else f[:n + 1] + fx = f(p.x[: n + 1]) if callable(f) else f[: n + 1] return np.sum(w * fx) # type: ignore[no-any-return] @@ -307,7 +307,7 @@ def _diff_caputo_l2(m: L2, f: ArrayOrScalarFunction, p: Points) -> Array: fx = f(p.x) if callable(f) else f if fx.shape[0] != p.size: raise ValueError( - f"Shape of 'f' does match points: got {f.shape} expected {p.shape}" + f"Shape of 'f' does match points: got {fx.shape} expected {p.shape}" ) df = np.empty_like(fx) @@ -319,6 +319,8 @@ def _diff_caputo_l2(m: L2, f: ArrayOrScalarFunction, p: Points) -> Array: w = quadrature_weights(m, p, n + 1) df[n] = np.sum(w * fx[: n + 1]) + return df + def _weights_l2(alpha: float, i: int | Array, k: int | Array) -> Array: return np.array((i - k) ** (2 - alpha) - (i - k - 1) ** (2 - alpha)) From 29cf7ff989b80f2e846c71f12fdffaf245f610ae Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Wed, 7 Aug 2024 20:26:55 +0300 Subject: [PATCH 09/26] feat: add fallbacks for diff functions --- docs/guide_differentiation.rst | 7 +- examples/guide-differentiation.py | 24 ++++- src/pycaputo/differentiation/__init__.py | 113 ++++++++++++++++++++++ src/pycaputo/differentiation/base.py | 114 ++++++++--------------- 4 files changed, 179 insertions(+), 79 deletions(-) diff --git a/docs/guide_differentiation.rst b/docs/guide_differentiation.rst index ba55917..9c0fc90 100644 --- a/docs/guide_differentiation.rst +++ b/docs/guide_differentiation.rst @@ -51,6 +51,9 @@ also have should not implement this method. However, for those that can, this offers a pleasant way of reusing and analyzing the weights. +* :func:`~pycaputo.differentiation.differentiation_matrix`: a function that + gathers all the quadrature weights into a matrix operator. + * :func:`~pycaputo.differentiation.diffs`: a method very similar to :func:`~pycaputo.differentiation.diff`, but that only computes the fractional operator at a given point on the grid. Note that not all methods can efficiently @@ -58,11 +61,13 @@ also have so they should not implement this functionality. * :func:`~pycaputo.differentiation.diff`: for many methods, the ``diff`` function - is a simple wrapper around :func:`~pycaputo.differentiation.diffs`, but this + can be a simple wrapper around :func:`~pycaputo.differentiation.diffs`, but this doesn't always make sense. For example, if we want to compute the Caputo fractional derivative on a uniform grid, this can be evaluated more efficiently using the FFT. +In general, we require that the :func:`~pycaputo.differentiation.diff` function be +implemented for all available methods, while the remaining functions are optional. These methods are all registered with the :func:`~functools.singledispatch` mechanism. An example setup is provided below. diff --git a/examples/guide-differentiation.py b/examples/guide-differentiation.py index e53608d..c613254 100644 --- a/examples/guide-differentiation.py +++ b/examples/guide-differentiation.py @@ -14,7 +14,7 @@ from pycaputo.derivatives import RiemannLiouvilleDerivative, Side from pycaputo.grid import Points -from pycaputo.typing import Array, ArrayOrScalarFunction +from pycaputo.typing import Array, ArrayOrScalarFunction, Scalar # {{{ @@ -42,7 +42,12 @@ def d(self) -> RiemannLiouvilleDerivative: # {{{ # [register-start] -from pycaputo.differentiation import diff, diffs, quadrature_weights +from pycaputo.differentiation import ( + diff, + differentiation_matrix, + diffs, + quadrature_weights, +) @quadrature_weights.register(RiemannLiouvilleDerivativeMethod) @@ -55,15 +60,24 @@ def _quadrature_weights_rl( return np.zeros(n, dtype=p.x.dtype) +@differentiation_matrix.register(RiemannLiouvilleDerivativeMethod) +def _differentiation_matrix_rl( + m: RiemannLiouvilleDerivativeMethod, + p: Points, +) -> Array: + # ... add an actual implementation here ... + return np.zeros((p.size, p.size), dtype=p.x.dtype) + + @diffs.register(RiemannLiouvilleDerivativeMethod) def _diffs_rl( m: RiemannLiouvilleDerivativeMethod, f: ArrayOrScalarFunction, p: Points, -) -> Array: - fx = f(p.x) if callable(f) else f + n: int, +) -> Scalar: # ... add an actual implementation here ... - return np.zeros_like(fx) + return np.array(0.0) @diff.register(RiemannLiouvilleDerivativeMethod) diff --git a/src/pycaputo/differentiation/__init__.py b/src/pycaputo/differentiation/__init__.py index c2e6ee3..fc013c0 100644 --- a/src/pycaputo/differentiation/__init__.py +++ b/src/pycaputo/differentiation/__init__.py @@ -3,14 +3,18 @@ from __future__ import annotations +import numpy as np + from pycaputo.derivatives import FractionalOperator from pycaputo.differentiation.base import ( DerivativeMethod, diff, + differentiation_matrix, diffs, quadrature_weights, ) from pycaputo.grid import Points +from pycaputo.utils import Array, ArrayOrScalarFunction, Scalar def guess_method_for_order( @@ -63,10 +67,119 @@ def guess_method_for_order( return m +def quadrature_weights_fallback(m: DerivativeMethod, p: Points, n: int) -> Array: + """Evaluate the quadrature weights for the method *m* at the points *p.x[n]*. + + This function attempts to construct the quadrature weights using + :func:`~pycaputo.differentiation.quadrature_weights`. If not available, it + falls back to constructing the full differentiation matrix using + :func:`~pycaputo.differentiation.differentiation_matrix` and selecting + the *nth* row. + + .. warning:: + + Falling back to :func:`~pycaputo.differentiation.differentiation_matrix` + will be significantly slower, as it required computing the full matrix. + Use this function with care. + """ + + try: + return quadrature_weights(m, p, n) + except NotImplementedError: + pass + + if not 0 <= n <= p.size: + raise IndexError(f"Index 'n' out of range: 0 <= {n} < {p.size}") + + W = differentiation_matrix(m, p) + return W[n, : n + 1] + + +def differentiation_matrix_fallback(m: DerivativeMethod, p: Points) -> Array: + """Evaluate the differentiation matrix for the method *m* at points *p*. + + This function attempts to constructs the differentiation matrix using + :func:`~pycaputo.differentiation.differentiation_matrix`. If not available, + it falls back to :func:`~pycaputo.differentiation.quadrature_weights`. + """ + try: + return differentiation_matrix(m, p) + except NotImplementedError: + pass + + n = p.size + W = np.zeros((n, n)) + W[0, :] = np.nan + + for i in range(2, n): + W[i, : i + 1] = quadrature_weights(m, p, i + 1) + + return W + + +def diffs_fallback( + m: DerivativeMethod, f: ArrayOrScalarFunction, p: Points, n: int +) -> Scalar: + """Evaluate the fractional derivative of *f* using method *m* at the point + *p.x[n]* using the underlying grid. + + This function attempts to evaluate the fractional derivative using + :func:`~pycaputo.differentiation.diffs`. If not available, it falls back + to evaluating the derivative at all points using + :func:`~pycaputo.differentiation.diff` and selecting the desired point. + + .. warning:: + + Falling back to the ``diff`` function will be significantly slower, + since all the points on the grid *p* are evaluated. Use this function + with care. + """ + try: + return diffs(m, f, p, n) + except NotImplementedError: + pass + + if not 0 <= n <= p.size: + raise IndexError(f"Index 'n' out of range: 0 <= {n} < {p.size}") + + result = diff(m, f, p) + return np.array(result[n]) + + +def diff_fallback(m: DerivativeMethod, f: ArrayOrScalarFunction, p: Points) -> Array: + """Evaluate the fractional derivative of *f* at *p* using method *m*. + + This function attempts to evaluate the fractional derivative using + :func:`~pycaputo.differentiation.diff`. If not available, it falls back + to evaluating the derivative point by points using + :func:`~pycaputo.differentiation.diffs`. + """ + try: + return diff(m, f, p) + except NotImplementedError: + pass + + n = 1 + df1 = diffs(m, f, p, n + 1) + df = np.empty((p.size, *df1.shape), dtype=df1.dtype) + df[0] = np.nan + df[1] = df1 + + for n in range(2, df.size): + df[n] = diffs(m, f, p, n + 1) + + return df + + __all__ = ( "DerivativeMethod", "diff", + "diff_fallback", + "differentiation_matrix", + "differentiation_matrix_fallback", "diffs", + "diffs_fallback", "guess_method_for_order", "quadrature_weights", + "quadrature_weights_fallback", ) diff --git a/src/pycaputo/differentiation/base.py b/src/pycaputo/differentiation/base.py index 5256d8b..e2edf37 100644 --- a/src/pycaputo/differentiation/base.py +++ b/src/pycaputo/differentiation/base.py @@ -7,8 +7,6 @@ from dataclasses import dataclass from functools import singledispatch -import numpy as np - from pycaputo.derivatives import FractionalOperator from pycaputo.grid import Points from pycaputo.typing import Array, ArrayOrScalarFunction, Scalar @@ -50,18 +48,45 @@ def quadrature_weights(m: DerivativeMethod, p: Points, n: int) -> Array: .. note:: Not all methods can compute the weights for a single evaluation in - this fashion. This function is mainly implemented for methods where the - *nth* row can be computed efficiently. + this fashion. This function should be implemented for methods where the + *nth* row of :math:`w` can be computed efficiently. Other methods may + implement :func:`differentiation_matrix` instead. + + The weights can also be constructed using :func:`differentiation_matrix`, if + implemented. For a safe function with this fallback, see + :func:`~pycaputo.differentiation.quadrature_weights_fallback`. :arg p: a grid on which to compute the quadrature weights for the point ``p.x[n]``. :arg n: the index of the point within *p*. """ + raise NotImplementedError( f"Cannot evaluate quadrature weights for method '{type(m).__name__}' " "('quadrature_weights' is not implemented)" ) +@singledispatch +def differentiation_matrix(m: DerivativeMethod, p: Points) -> Array: + """Evaluate the differentiation matrix for the method *m* at points *p*. + + Without additional knowledge, the value of the derivative at *p.a* is not + known. Therefore, the first row of the matrix should not be used in computation. + + This matrix can also be constructed using :func:`quadrature_weights`, if + implemented. For a safe function with this fallback, see + :func:`~pycaputo.differentiation.differentiation_matrix_fallback`. + + :returns: a two-dimensional array of shape ``(n, n)``, where *n* is the + number of points in *p*. + """ + + raise NotImplementedError( + f"Cannot evaluate quadrature weights for method '{type(m).__name__}' " + "('differentiation_matrix' is not implemented)" + ) + + @singledispatch def diffs(m: DerivativeMethod, f: ArrayOrScalarFunction, p: Points, n: int) -> Scalar: r"""Evaluate the fractional derivative of *f* using method *m* at the point @@ -71,53 +96,32 @@ def diffs(m: DerivativeMethod, f: ArrayOrScalarFunction, p: Points, n: int) -> S This is in contract to :func:`diff`, which evaluates the fractional derivative at all the points in the grid *p*. - .. note:: + The function *f* can be provided as a callable of as an array of values + evaluated on the grid *p*. However, not all methods support both variants + (e.g. if they require evaluating the function at additional points). - The function *f* can be provided as a callable of as an array of values - evaluated on the grid *p*. However, not all methods support both variants - (e.g. if they require evaluating the function at additional points). + The function values can also be obtained from :func:`diff` directly, if + implemented. For a safe function with this fallback, see + :func:`~pycaputo.differentiation.diffs_fallback`. :arg m: method used to evaluate the derivative. :arg f: a simple function for which to evaluate the derivative. :arg p: an array of points at which to evaluate the derivative. """ - try: - result = diff(m, f, p) - except NotImplementedError: - raise NotImplementedError( - "Cannot evaluate pointwise fractional derivative with method " - f"'{type(m).__name__}' ('diffs' is not implemented)" - ) from None - - from warnings import warn - - warn( - f"'diffs' is not implemented for '{type(m).__name__}' (aka '{m.name}'). " - "Falling back to 'diff', which is likely significantly slower! Use " - "'diff' directly if this is acceptable.", - stacklevel=2, + raise NotImplementedError( + "Cannot evaluate pointwise fractional derivative with method " + f"'{type(m).__name__}' ('diffs' is not implemented)" ) - if not 0 <= n <= p.size: - raise IndexError(f"Index 'n' out of range: 0 <= {n} < {p.size}") - - return np.array(result[n]) - @singledispatch def diff(m: DerivativeMethod, f: ArrayOrScalarFunction, p: Points) -> Array: """Evaluate the fractional derivative of *f* at *p* using method *m*. - This function uses :func:`diffs` to evaluate the fractional derivative at - all the points in *p*. In most cases it is a trivial wrapper. + The function values can also be obtained from :func:`diffs` directly, if + implemented. For a safe function with this fallback, see + :func:`~pycaputo.differentiation.diff_fallback`. - .. note:: - - All methods are expected to implement this method, while - :func:`quadrature_weights` and :func:`diffs` are optional and not - supported by all methods. - - :arg m: method used to evaluate the derivative. :arg f: a simple function for which to evaluate the derivative. :arg p: an array of points at which to evaluate the derivative. """ @@ -125,39 +129,3 @@ def diff(m: DerivativeMethod, f: ArrayOrScalarFunction, p: Points) -> Array: f"Cannot evaluate fractional derivative with method '{type(m).__name__}' " "('diff' is not implemented)" ) - - -def quadrature_matrix(m: DerivativeMethod, p: Points, w00: float = 0.0) -> Array: - """Evaluate the quadrature matrix for a given method *m* at points *p*. - - This requires that the method implements the :func:`quadrature_weights`. As - fractional operator are usually convolutional, this results in a lower - triangular matrix. - - The value of the derivative at ``p.a`` is not defined. The value of the - matrix at ``W[0, 0]`` can be controlled using *w00*, which is set to 0 by - default. This default choice results in a singular operator and may not - always be desired. - - Using this matrix operator, we can evaluate the derivative equivalently as - - .. code:: python - - mat = quadrature_matrix(m, p) - - df_mat = mat @ f - df_fun = diff(m, f, p) - assert np.allclose(df_mat, df_fun) - - :returns: a two-dimensional array of shape ``(n, n)``, where *n* is the - number of points in *p*. - """ - - n = p.size - W = np.zeros((n, n)) - W[0, 0] = w00 - - for i in range(1, n): - W[i, : i + 1] = quadrature_weights(m, p, i + 1) - - return W From 8b775c5dd54ebd0f7d33913af8b5b8ad92ff7398 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Wed, 7 Aug 2024 21:00:11 +0300 Subject: [PATCH 10/26] feat: more work on the L2 method --- src/pycaputo/differentiation/__init__.py | 2 +- src/pycaputo/differentiation/caputo.py | 36 ++++++++++++++++-------- src/pycaputo/fode/caputo.py | 2 ++ 3 files changed, 27 insertions(+), 13 deletions(-) diff --git a/src/pycaputo/differentiation/__init__.py b/src/pycaputo/differentiation/__init__.py index fc013c0..83ce94c 100644 --- a/src/pycaputo/differentiation/__init__.py +++ b/src/pycaputo/differentiation/__init__.py @@ -160,7 +160,7 @@ def diff_fallback(m: DerivativeMethod, f: ArrayOrScalarFunction, p: Points) -> A pass n = 1 - df1 = diffs(m, f, p, n + 1) + df1 = np.array(diffs(m, f, p, n + 1)) df = np.empty((p.size, *df1.shape), dtype=df1.dtype) df[0] = np.nan df[1] = df1 diff --git a/src/pycaputo/differentiation/caputo.py b/src/pycaputo/differentiation/caputo.py index 34c409a..a1c1e12 100644 --- a/src/pycaputo/differentiation/caputo.py +++ b/src/pycaputo/differentiation/caputo.py @@ -69,7 +69,9 @@ def _caputo_piecewise_constant_integral(p: Points, n: int, alpha: float) -> Arra .. math:: - a_{nk} = \int_{x_k}^{x_{k + 1}} \frac{1}{(x_n - s)^\alpha} ds + a_{nk} = + \frac{1}{\Gamma(1 - \alpha) (x_{k + 1} - x_k)} + \int_{x_k}^{x_{k + 1}} \frac{1}{(x_n - s)^\alpha} ds for :math:`0 \le k \le n - 1`. """ @@ -151,7 +153,8 @@ class ModifiedL1(CaputoMethod): This method is defined in Section 4.1.1 (III) from [Li2020]_. Note that this method evaluates the fractional derivative at the midpoints - :math:`(x_i + x_{i - 1}) / 2` for any given input grid. + :math:`(x_i + x_{i - 1}) / 2` for any given input grid except + :class:`~pycaputo.grids.MidPoints`. As noted in [Li2020]_, this method has the same order of convergence as the standard L1 method. However the weights can be used to construct a @@ -223,7 +226,13 @@ class L2(CaputoMethod): r"""Implements the L2 method for the Caputo fractional derivative of order :math:`\alpha \in (1, 2)`. - This method is defined in Section 4.1.2 from [Li2020]_ for uniform grids. + This method is defined in Section 4.1.2 from [Li2020]_ for general + non-uniform grids. + + .. note:: + + Unlike the method from [Li2020_], we do not assume knowledge of points + outside of the domain. Instead a biased stencil is used at the boundary. """ if __debug__: @@ -232,18 +241,21 @@ def __post_init__(self) -> None: super().__post_init__() if not 1 < self.alpha < 2: raise ValueError( - f"'{type(self).__name__}' only supports 0 < alpha < 1: {self.alpha}" + f"'{type(self).__name__}' only supports 1 < alpha < 2: {self.alpha}" ) -def _caputo_d2_coefficients(x: Array) -> Array: +def _caputo_d2_coefficients(x: Array, s: Scalar) -> tuple[Scalar, ...]: + """Get coefficients for a fourth-order approximation of the second derivative + at the boundary. + """ x0, x1, x2, x3 = x[:4] - a0 = 2.0 * (3.0 * x0 - x1 - x2 - x3) / ((x0 - x1) * (x0 - x2) * (x0 - x3)) - b0 = 2.0 * (x3 - 2.0 * x0 + x2) / ((x0 - x1) * (x1 - x2) * (x1 - x3)) - c0 = 2.0 * (x3 - 2.0 * x0 + x1) / ((x0 - x2) * (x2 - x1) * (x2 - x3)) - d0 = 2.0 * (x2 - 2.0 * x0 + x1) / ((x0 - x3) * (x3 - x1) * (x3 - x2)) + a0 = 2.0 * (3.0 * s - x1 - x2 - x3) / ((x0 - x1) * (x0 - x2) * (x0 - x3)) + b0 = 2.0 * (x0 + x2 + x3 - 3.0 * s) / ((x0 - x1) * (x1 - x2) * (x1 - x3)) + c0 = 2.0 * (x0 + x1 + x3 - 3.0 * s) / ((x0 - x2) * (x2 - x1) * (x2 - x3)) + d0 = 2.0 * (x0 + x1 + x2 - 3.0 * s) / ((x0 - x3) * (x3 - x1) * (x3 - x2)) - return np.array([a0, b0, c0, d0]) + return a0, b0, c0, d0 @quadrature_weights.register(L2) @@ -262,8 +274,8 @@ def _quadrature_weights_caputo_l2(m: L2, p: Points, n: int) -> Array: b = 1.0 / (dx[1:-1] * dxm) # get finite difference coefficients for boundary stencils - a_l, b_l, c_l, d_l = _caputo_d2_coefficients(p.x) - d_r, c_r, b_r, a_r = _caputo_d2_coefficients(p.x[::-1]) + a_l, b_l, c_l, d_l = _caputo_d2_coefficients(p.x, p.x[0]) + d_r, c_r, b_r, a_r = _caputo_d2_coefficients(p.x[::-1], p.x[-1]) wi = _caputo_piecewise_constant_integral(p, n, m.alpha) diff --git a/src/pycaputo/fode/caputo.py b/src/pycaputo/fode/caputo.py index 2b611d3..90130f9 100644 --- a/src/pycaputo/fode/caputo.py +++ b/src/pycaputo/fode/caputo.py @@ -39,6 +39,8 @@ def _truncation_error( ) -> Array: from pycaputo.controller import JannelliIntegralController + # FIXME: this should not be our job: either let the controller or the method + # figure out how to compute the truncation error estimates. assert t > tprev if isinstance(c, JannelliIntegralController): trunc = np.array( From 07d1a71776499b11372aae09e7f48a50953fb44b Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Wed, 7 Aug 2024 21:53:38 +0300 Subject: [PATCH 11/26] fix: typing imports --- examples/caputo-derivative-l1-spectrum.py | 2 +- src/pycaputo/differentiation/__init__.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/caputo-derivative-l1-spectrum.py b/examples/caputo-derivative-l1-spectrum.py index dec8d22..4661ed1 100644 --- a/examples/caputo-derivative-l1-spectrum.py +++ b/examples/caputo-derivative-l1-spectrum.py @@ -24,7 +24,7 @@ from pycaputo.differentiation.caputo import L1 from pycaputo.grid import make_uniform_points from pycaputo.logging import get_logger -from pycaputo.utils import Array +from pycaputo.typing import Array logger = get_logger("caputo_derivative_l1_spectrum") diff --git a/src/pycaputo/differentiation/__init__.py b/src/pycaputo/differentiation/__init__.py index 83ce94c..65ccf09 100644 --- a/src/pycaputo/differentiation/__init__.py +++ b/src/pycaputo/differentiation/__init__.py @@ -14,7 +14,7 @@ quadrature_weights, ) from pycaputo.grid import Points -from pycaputo.utils import Array, ArrayOrScalarFunction, Scalar +from pycaputo.typing import Array, ArrayOrScalarFunction, Scalar def guess_method_for_order( From 7bd1e5feffe2ceccba91638866da55065e74880a Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Thu, 8 Aug 2024 09:35:52 +0300 Subject: [PATCH 12/26] feat(diff): move fallbacks --- src/pycaputo/differentiation/__init__.py | 57 ++-------------------- src/pycaputo/differentiation/base.py | 62 ++++++++++++++++++------ 2 files changed, 52 insertions(+), 67 deletions(-) diff --git a/src/pycaputo/differentiation/__init__.py b/src/pycaputo/differentiation/__init__.py index 65ccf09..f809c17 100644 --- a/src/pycaputo/differentiation/__init__.py +++ b/src/pycaputo/differentiation/__init__.py @@ -8,6 +8,7 @@ from pycaputo.derivatives import FractionalOperator from pycaputo.differentiation.base import ( DerivativeMethod, + FunctionCallableError, diff, differentiation_matrix, diffs, @@ -95,28 +96,6 @@ def quadrature_weights_fallback(m: DerivativeMethod, p: Points, n: int) -> Array return W[n, : n + 1] -def differentiation_matrix_fallback(m: DerivativeMethod, p: Points) -> Array: - """Evaluate the differentiation matrix for the method *m* at points *p*. - - This function attempts to constructs the differentiation matrix using - :func:`~pycaputo.differentiation.differentiation_matrix`. If not available, - it falls back to :func:`~pycaputo.differentiation.quadrature_weights`. - """ - try: - return differentiation_matrix(m, p) - except NotImplementedError: - pass - - n = p.size - W = np.zeros((n, n)) - W[0, :] = np.nan - - for i in range(2, n): - W[i, : i + 1] = quadrature_weights(m, p, i + 1) - - return W - - def diffs_fallback( m: DerivativeMethod, f: ArrayOrScalarFunction, p: Points, n: int ) -> Scalar: @@ -130,9 +109,9 @@ def diffs_fallback( .. warning:: - Falling back to the ``diff`` function will be significantly slower, - since all the points on the grid *p* are evaluated. Use this function - with care. + Falling back to the :func:`~pycaputo.differentiation.diff`` function + will be significantly slower, since all the points on the grid *p* are + evaluated. Use this function with care. """ try: return diffs(m, f, p, n) @@ -146,37 +125,11 @@ def diffs_fallback( return np.array(result[n]) -def diff_fallback(m: DerivativeMethod, f: ArrayOrScalarFunction, p: Points) -> Array: - """Evaluate the fractional derivative of *f* at *p* using method *m*. - - This function attempts to evaluate the fractional derivative using - :func:`~pycaputo.differentiation.diff`. If not available, it falls back - to evaluating the derivative point by points using - :func:`~pycaputo.differentiation.diffs`. - """ - try: - return diff(m, f, p) - except NotImplementedError: - pass - - n = 1 - df1 = np.array(diffs(m, f, p, n + 1)) - df = np.empty((p.size, *df1.shape), dtype=df1.dtype) - df[0] = np.nan - df[1] = df1 - - for n in range(2, df.size): - df[n] = diffs(m, f, p, n + 1) - - return df - - __all__ = ( "DerivativeMethod", + "FunctionCallableError", "diff", - "diff_fallback", "differentiation_matrix", - "differentiation_matrix_fallback", "diffs", "diffs_fallback", "guess_method_for_order", diff --git a/src/pycaputo/differentiation/base.py b/src/pycaputo/differentiation/base.py index e2edf37..bef904d 100644 --- a/src/pycaputo/differentiation/base.py +++ b/src/pycaputo/differentiation/base.py @@ -7,11 +7,20 @@ from dataclasses import dataclass from functools import singledispatch +import numpy as np + from pycaputo.derivatives import FractionalOperator from pycaputo.grid import Points from pycaputo.typing import Array, ArrayOrScalarFunction, Scalar +class FunctionCallableError(TypeError): + """An error raised by :func:`diffs` or :func:`diff` when the argument provided + for evaluation must be a callable + (satisfy :class:`~pycaputo.typing.ScalarFunction`). + """ + + @dataclass(frozen=True) class DerivativeMethod(ABC): """A generic method used to evaluate a fractional derivative at a set of points.""" @@ -72,19 +81,30 @@ def differentiation_matrix(m: DerivativeMethod, p: Points) -> Array: Without additional knowledge, the value of the derivative at *p.a* is not known. Therefore, the first row of the matrix should not be used in computation. - - This matrix can also be constructed using :func:`quadrature_weights`, if - implemented. For a safe function with this fallback, see - :func:`~pycaputo.differentiation.differentiation_matrix_fallback`. + By default, this function is implemented in terms of :func:`quadrature_weights` :returns: a two-dimensional array of shape ``(n, n)``, where *n* is the number of points in *p*. """ + i = 1 - raise NotImplementedError( - f"Cannot evaluate quadrature weights for method '{type(m).__name__}' " - "('differentiation_matrix' is not implemented)" - ) + try: + w1 = quadrature_weights(m, p, i + 1) + except NotImplementedError: + raise NotImplementedError( + f"Cannot evaluate differentiation matrix for method '{type(m).__name__}' " + "('differentiation_matrix' is not implemented)" + ) from None + + n = p.size + W = np.zeros((n, n), dtype=w1.dtype) + W[0, :1] = np.nan + W[0, :2] = w1 + + for i in range(2, W.shape[0]): + W[i, : i + 1] = quadrature_weights(m, p, i + 1) + + return W @singledispatch @@ -118,14 +138,26 @@ def diffs(m: DerivativeMethod, f: ArrayOrScalarFunction, p: Points, n: int) -> S def diff(m: DerivativeMethod, f: ArrayOrScalarFunction, p: Points) -> Array: """Evaluate the fractional derivative of *f* at *p* using method *m*. - The function values can also be obtained from :func:`diffs` directly, if - implemented. For a safe function with this fallback, see - :func:`~pycaputo.differentiation.diff_fallback`. + By default, this function is implemented in terms of :func:`diffs`. :arg f: a simple function for which to evaluate the derivative. :arg p: an array of points at which to evaluate the derivative. """ - raise NotImplementedError( - f"Cannot evaluate fractional derivative with method '{type(m).__name__}' " - "('diff' is not implemented)" - ) + n = 1 + + try: + df1 = np.array(diffs(m, f, p, n + 1)) + except NotImplementedError: + raise NotImplementedError( + f"Cannot evaluate fractional derivative with method '{type(m).__name__}' " + "('diff' is not implemented)" + ) from None + + df = np.empty((p.size, *df1.shape), dtype=df1.dtype) + df[0] = np.nan + df[1] = df1 + + for n in range(2, df.size): + df[n] = diffs(m, f, p, n + 1) + + return df From 62bf5f04e43a8c566b69fb5cae4558b7fdde7e24 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Thu, 8 Aug 2024 09:54:20 +0300 Subject: [PATCH 13/26] docs: update diff docs --- docs/guide_differentiation.rst | 27 +++++++++++++------------- src/pycaputo/differentiation/base.py | 9 +++++++-- src/pycaputo/differentiation/caputo.py | 6 +++--- 3 files changed, 24 insertions(+), 18 deletions(-) diff --git a/docs/guide_differentiation.rst b/docs/guide_differentiation.rst index 9c0fc90..210ebcf 100644 --- a/docs/guide_differentiation.rst +++ b/docs/guide_differentiation.rst @@ -52,19 +52,20 @@ also have pleasant way of reusing and analyzing the weights. * :func:`~pycaputo.differentiation.differentiation_matrix`: a function that - gathers all the quadrature weights into a matrix operator. - -* :func:`~pycaputo.differentiation.diffs`: a method very similar to - :func:`~pycaputo.differentiation.diff`, but that only computes the fractional - operator at a given point on the grid. Note that not all methods can efficiently - evaluate the fractional operator at a single point (e.g. global spectral methods), - so they should not implement this functionality. - -* :func:`~pycaputo.differentiation.diff`: for many methods, the ``diff`` function - can be a simple wrapper around :func:`~pycaputo.differentiation.diffs`, but this - doesn't always make sense. For example, if we want to compute the Caputo - fractional derivative on a uniform grid, this can be evaluated more efficiently - using the FFT. + gathers all the quadrature weights into a matrix operator. By default, this + function is implemented in terms of ``quadrature_weights``, but can be + overwritten by more efficient implementations. + +* :func:`~pycaputo.differentiation.diffs`: a function that computes the + fractional operator at a given point on a grid. Note that not all methods can + efficiently evaluate the fractional operator at a single point (e.g. global + spectral methods), so they should not implement this functionality. + +* :func:`~pycaputo.differentiation.diff`: a function that computes the fractional + operator at a set of points. By default this is implemented in terms of the + ``diffs`` function, but this doesn't always make sense. For example, if we want + to compute the Caputo fractional derivative on a uniform grid, this can be + evaluated more efficiently using the FFT. In general, we require that the :func:`~pycaputo.differentiation.diff` function be implemented for all available methods, while the remaining functions are optional. diff --git a/src/pycaputo/differentiation/base.py b/src/pycaputo/differentiation/base.py index bef904d..e5d099b 100644 --- a/src/pycaputo/differentiation/base.py +++ b/src/pycaputo/differentiation/base.py @@ -67,6 +67,7 @@ def quadrature_weights(m: DerivativeMethod, p: Points, n: int) -> Array: :arg p: a grid on which to compute the quadrature weights for the point ``p.x[n]``. :arg n: the index of the point within *p*. + :returns: an array of quadrature weights of shape ``(n,)``. """ raise NotImplementedError( @@ -127,6 +128,8 @@ def diffs(m: DerivativeMethod, f: ArrayOrScalarFunction, p: Points, n: int) -> S :arg m: method used to evaluate the derivative. :arg f: a simple function for which to evaluate the derivative. :arg p: an array of points at which to evaluate the derivative. + :returns: a scalar representing the fractional derivative approximation at + the *nth* point in *p*. """ raise NotImplementedError( "Cannot evaluate pointwise fractional derivative with method " @@ -140,8 +143,10 @@ def diff(m: DerivativeMethod, f: ArrayOrScalarFunction, p: Points) -> Array: By default, this function is implemented in terms of :func:`diffs`. - :arg f: a simple function for which to evaluate the derivative. - :arg p: an array of points at which to evaluate the derivative. + :arg f: an array or function for which to evaluate the derivative. + :arg p: a set points at which to evaluate the derivative. + :returns: an array of shape ``(n,)``, where *n* is the number of points in *p*, + that contains the evaluation of the fractional derivative at each point. """ n = 1 diff --git a/src/pycaputo/differentiation/caputo.py b/src/pycaputo/differentiation/caputo.py index a1c1e12..bcf807f 100644 --- a/src/pycaputo/differentiation/caputo.py +++ b/src/pycaputo/differentiation/caputo.py @@ -88,12 +88,12 @@ def _caputo_piecewise_constant_integral(p: Points, n: int, alpha: float) -> Arra @quadrature_weights.register(L1) def _quadrature_weights_caputo_l1(m: L1, p: Points, n: int) -> Array: - if n < 0: - n = p.size + n - if not 0 <= n <= p.size: raise IndexError(f"Index 'n' out of range: 0 <= {n} < {p.size}") + if n == 0: + return np.array([]) + w = np.empty(n, dtype=p.x.dtype) w[n - 1] = 0.0 w[: n - 1] = _caputo_piecewise_constant_integral(p, n, m.alpha) From 5bccc7734d6d752798bacb76be11b418046f7580 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Thu, 8 Aug 2024 10:12:05 +0300 Subject: [PATCH 14/26] feat: add L1D method that takes a callable only --- src/pycaputo/differentiation/caputo.py | 97 +++++++++++++++++++++++--- 1 file changed, 89 insertions(+), 8 deletions(-) diff --git a/src/pycaputo/differentiation/caputo.py b/src/pycaputo/differentiation/caputo.py index bcf807f..fbe681d 100644 --- a/src/pycaputo/differentiation/caputo.py +++ b/src/pycaputo/differentiation/caputo.py @@ -19,7 +19,13 @@ Scalar, ) -from .base import DerivativeMethod, diff, diffs, quadrature_weights +from .base import ( + DerivativeMethod, + FunctionCallableError, + diff, + diffs, + quadrature_weights, +) logger = get_logger(__name__) @@ -70,18 +76,17 @@ def _caputo_piecewise_constant_integral(p: Points, n: int, alpha: float) -> Arra .. math:: a_{nk} = - \frac{1}{\Gamma(1 - \alpha) (x_{k + 1} - x_k)} + \frac{1}{\Gamma(1 - \alpha)} \int_{x_k}^{x_{k + 1}} \frac{1}{(x_n - s)^\alpha} ds for :math:`0 \le k \le n - 1`. """ - x, dx = p.x[:n], p.dx[: n - 1] + x = p.x[:n] xn = x[-1] return np.array( ((xn - x[:-1]) ** (1 - alpha) - (xn - x[1:]) ** (1 - alpha)) - / dx / math.gamma(2 - alpha) ) @@ -94,9 +99,11 @@ def _quadrature_weights_caputo_l1(m: L1, p: Points, n: int) -> Array: if n == 0: return np.array([]) + a = _caputo_piecewise_constant_integral(p, n, m.alpha) / p.dx[: n - 1] + w = np.empty(n, dtype=p.x.dtype) w[n - 1] = 0.0 - w[: n - 1] = _caputo_piecewise_constant_integral(p, n, m.alpha) + w[: n - 1] = a w[1:n] = w[: n - 1] - w[1:n] w[0] = -w[0] @@ -105,9 +112,6 @@ def _quadrature_weights_caputo_l1(m: L1, p: Points, n: int) -> Array: @diffs.register(L1) def _diffs_caputo_l1(m: L1, f: ArrayOrScalarFunction, p: Points, n: int) -> Scalar: - if n < 0: - n = p.size + n - if not 0 <= n <= p.size: raise IndexError(f"Index 'n' out of range: 0 <= {n} < {p.size}") @@ -143,6 +147,83 @@ def _diff_caputo_l1(m: L1, f: ArrayOrScalarFunction, p: Points) -> Array: # }}} +# {{{ L1D + + +@dataclass(frozen=True) +class L1D(CaputoMethod): + r"""Implements the L1 method for the Caputo fractional derivative + of order :math:`\alpha \in (0, 1)`. + + This method is defined in Section 4.1.1 (II) from [Li2020]_ for general + non-uniform grids. This method is equivalent to :class:`L1`, but evaluation + requires explicit knowledge of the first derivative. + """ + + if __debug__: + + def __post_init__(self) -> None: + super().__post_init__() + if not 0 < self.alpha < 1: + raise ValueError( + f"'{type(self).__name__}' only supports 0 < alpha < 1: {self.alpha}" + ) + + +@quadrature_weights.register(L1D) +def _quadrature_weights_caputo_l1d(m: L1D, p: Points, n: int) -> Array: + if not 0 <= n <= p.size: + raise IndexError(f"Index 'n' out of range: 0 <= {n} < {p.size}") + + if n == 0: + return np.array([]) + + return _caputo_piecewise_constant_integral(p, n, m.alpha) + + +@diffs.register(L1D) +def _diffs_caputo_l1d(m: L1D, f: ArrayOrScalarFunction, p: Points, n: int) -> Scalar: + if not 0 <= n <= p.size: + raise IndexError(f"Index 'n' out of range: 0 <= {n} < {p.size}") + + if n == 0: + return np.array([np.nan]) + + if not isinstance(f, DifferentiableScalarFunction): + raise FunctionCallableError( + f"The '{type(m).__name__}' method requires evaluating the first " + f"derivative. 'f' is not callable: {type(f)}" + ) + + w = quadrature_weights(m, p, n + 1) + dfx = f((p.x[1 : n + 1] + p.x[:n]) / 2, d=1) + + return np.sum(w * dfx) # type: ignore[no-any-return] + + +@diff.register(L1D) +def _diff_caputo_l1d(m: L1D, f: ArrayOrScalarFunction, p: Points) -> Array: + if not isinstance(f, DifferentiableScalarFunction): + raise FunctionCallableError( + f"The '{type(m).__name__}' method requires evaluating the first " + f"derivative. 'f' is not callable: {type(f)}" + ) + + dfx = f((p.x[1:] + p.x[:-1]) / 2, d=1) + + df = np.empty((p.size, *dfx.shape[1:]), dtype=dfx.dtype) + df[0] = np.nan + + for n in range(1, df.size): + w = quadrature_weights(m, p, n + 1) + df[n] = np.sum(w * dfx[: n + 1]) + + return df + + +# }}} + + # {{{ ModifiedL1 From e778dfae280b9d434b1c922639d4fc89a818fc54 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Thu, 8 Aug 2024 21:15:48 +0300 Subject: [PATCH 15/26] feat: finish up L2 reimplementation --- src/pycaputo/differentiation/base.py | 16 +-- src/pycaputo/differentiation/caputo.py | 161 ++++++++++++------------- src/pycaputo/grid.py | 5 + 3 files changed, 85 insertions(+), 97 deletions(-) diff --git a/src/pycaputo/differentiation/base.py b/src/pycaputo/differentiation/base.py index e5d099b..811eebd 100644 --- a/src/pycaputo/differentiation/base.py +++ b/src/pycaputo/differentiation/base.py @@ -87,24 +87,20 @@ def differentiation_matrix(m: DerivativeMethod, p: Points) -> Array: :returns: a two-dimensional array of shape ``(n, n)``, where *n* is the number of points in *p*. """ - i = 1 + n = p.size + W = np.zeros((n, n), dtype=p.dtype) + W[0, 0] = np.nan try: - w1 = quadrature_weights(m, p, i + 1) + for i in range(1, W.shape[0]): + w = quadrature_weights(m, p, i) + W[i, : w.size] = w except NotImplementedError: raise NotImplementedError( f"Cannot evaluate differentiation matrix for method '{type(m).__name__}' " "('differentiation_matrix' is not implemented)" ) from None - n = p.size - W = np.zeros((n, n), dtype=w1.dtype) - W[0, :1] = np.nan - W[0, :2] = w1 - - for i in range(2, W.shape[0]): - W[i, : i + 1] = quadrature_weights(m, p, i + 1) - return W diff --git a/src/pycaputo/differentiation/caputo.py b/src/pycaputo/differentiation/caputo.py index fbe681d..a7fdb82 100644 --- a/src/pycaputo/differentiation/caputo.py +++ b/src/pycaputo/differentiation/caputo.py @@ -97,11 +97,11 @@ def _quadrature_weights_caputo_l1(m: L1, p: Points, n: int) -> Array: raise IndexError(f"Index 'n' out of range: 0 <= {n} < {p.size}") if n == 0: - return np.array([]) + return np.array([], dtype=p.dtype) a = _caputo_piecewise_constant_integral(p, n, m.alpha) / p.dx[: n - 1] - w = np.empty(n, dtype=p.x.dtype) + w = np.empty(n, dtype=p.dtype) w[n - 1] = 0.0 w[: n - 1] = a w[1:n] = w[: n - 1] - w[1:n] @@ -178,7 +178,7 @@ def _quadrature_weights_caputo_l1d(m: L1D, p: Points, n: int) -> Array: if n == 0: return np.array([]) - return _caputo_piecewise_constant_integral(p, n, m.alpha) + return _caputo_piecewise_constant_integral(p, n + 1, m.alpha) @diffs.register(L1D) @@ -195,7 +195,7 @@ def _diffs_caputo_l1d(m: L1D, f: ArrayOrScalarFunction, p: Points, n: int) -> Sc f"derivative. 'f' is not callable: {type(f)}" ) - w = quadrature_weights(m, p, n + 1) + w = quadrature_weights(m, p, n) dfx = f((p.x[1 : n + 1] + p.x[:n]) / 2, d=1) return np.sum(w * dfx) # type: ignore[no-any-return] @@ -215,8 +215,8 @@ def _diff_caputo_l1d(m: L1D, f: ArrayOrScalarFunction, p: Points) -> Array: df[0] = np.nan for n in range(1, df.size): - w = quadrature_weights(m, p, n + 1) - df[n] = np.sum(w * dfx[: n + 1]) + w = quadrature_weights(m, p, n) + df[n] = np.sum(w * dfx[:n]) return df @@ -326,54 +326,77 @@ def __post_init__(self) -> None: ) -def _caputo_d2_coefficients(x: Array, s: Scalar) -> tuple[Scalar, ...]: - """Get coefficients for a fourth-order approximation of the second derivative +def _caputo_d2_boundary_coefficients(x: Array, s: Scalar) -> Array: + r"""Get coefficients for a fourth-order approximation of the second derivative at the boundary. + + .. math:: + + \frac{\partial f}{\partial x}(s) = + c_0 f(x_0) + c_1 f(x_1) + c_2 f(x_2) + c_3 f(x_3) + """ x0, x1, x2, x3 = x[:4] - a0 = 2.0 * (3.0 * s - x1 - x2 - x3) / ((x0 - x1) * (x0 - x2) * (x0 - x3)) - b0 = 2.0 * (x0 + x2 + x3 - 3.0 * s) / ((x0 - x1) * (x1 - x2) * (x1 - x3)) - c0 = 2.0 * (x0 + x1 + x3 - 3.0 * s) / ((x0 - x2) * (x2 - x1) * (x2 - x3)) - d0 = 2.0 * (x0 + x1 + x2 - 3.0 * s) / ((x0 - x3) * (x3 - x1) * (x3 - x2)) + c0 = 2.0 * (3.0 * s - x1 - x2 - x3) / ((x0 - x1) * (x0 - x2) * (x0 - x3)) + c1 = 2.0 * (x0 + x2 + x3 - 3.0 * s) / ((x0 - x1) * (x1 - x2) * (x1 - x3)) + c2 = 2.0 * (x0 + x1 + x3 - 3.0 * s) / ((x0 - x2) * (x2 - x1) * (x2 - x3)) + c3 = 2.0 * (x0 + x1 + x2 - 3.0 * s) / ((x0 - x3) * (x3 - x1) * (x3 - x2)) + + return np.array([c0, c1, c2, c3]) + + +def _caputo_d2_interior_coefficients(p: Points, n: int) -> tuple[Array, Array, Array]: + r"""Get coefficients for the second-order approximation of the second + derivative at the interior. + + .. math:: + + \frac{\partial f}{\partial x}(s) = + d_- f_{k - 1} + d_m f_k + d_+ f_{k + 1} + + which will simplify to the standard :math:`(1, -2, 1)` in the uniform case. + """ + + dx = p.dx[:n] + dxm = p.dxm[: n - 1] + + d_l = 1.0 / (dx[:-1] * dxm) + d_r = 1.0 / (dx[1:] * dxm) - return a0, b0, c0, d0 + return d_l, -(d_l + d_r), d_r @quadrature_weights.register(L2) def _quadrature_weights_caputo_l2(m: L2, p: Points, n: int) -> Array: - if n < 0: - n = p.size + n - if not 0 <= n <= p.size: raise IndexError(f"Index 'n' out of range: 0 <= {n} < {p.size}") - dx = p.dx[: n - 1] - dxm = (dx[1:] + dx[:-1]) / 2.0 - - # get finite difference coefficients for center stencil - a = 1.0 / (dx[:-2] * dxm) - b = 1.0 / (dx[1:-1] * dxm) - - # get finite difference coefficients for boundary stencils - a_l, b_l, c_l, d_l = _caputo_d2_coefficients(p.x, p.x[0]) - d_r, c_r, b_r, a_r = _caputo_d2_coefficients(p.x[::-1], p.x[-1]) - - wi = _caputo_piecewise_constant_integral(p, n, m.alpha) - - w = np.empty_like(p.x[:n]) - # center stencil - w[2 : n - 1] = a * wi[2:] - (a + b) * w[1:-1] + b * w[:-2] - # add left boundary stencil - w[0] = a_l + a[1] * wi[1] - w[1] = b_l + a[2] * wi[2] - (a[1] + b[1]) * wi[1] - w[2] += c_l - w[3] += d_l - # add right-boundary stencil - if n == p.size: - w[-4] += d_r - w[-3] += c_r - w[-2] = b_r + a[-2] * wi[-2] - (a[-1] + b[-1]) * wi[-1] - w[-1] = a_r + a[-1] * wi[-1] + if n == 0: + return np.array([], dtype=p.dtype) + + # weights have size at least 4 for the initial biased stencil. + w = np.zeros(max(4, n + 1), dtype=p.dtype) + + # coefficients + a = _caputo_piecewise_constant_integral(p, n + 1, m.alpha - 1) + c_l = _caputo_d2_boundary_coefficients(p.x, p.x[0]) + d_l, d_m, d_r = _caputo_d2_interior_coefficients(p, n) + + # add boundary stencil + w[:4] += a[0] * c_l + + # add interior stencils + if n > 1: + w[0] += a[1] * d_l[0] + w[1] += a[1] * d_m[0] + w[n] += a[-1] * d_r[-1] + + if n > 2: + w[1] += a[2] * d_l[1] + w[n - 1] += a[-1] * d_m[-1] + a[-2] * d_r[-2] + + if n > 3: + w[2 : n - 1] += a[1:-2] * d_r[:-2] + a[2:-1] * d_m[1:-1] + a[3:] * d_l[2:] return w @@ -389,8 +412,8 @@ def _diffs_caputo_l2(m: L2, f: ArrayOrScalarFunction, p: Points, n: int) -> Arra if n == 0: return np.array([np.nan]) - w = quadrature_weights(m, p, n + 1) - fx = f(p.x[: n + 1]) if callable(f) else f[: n + 1] + w = quadrature_weights(m, p, n) + fx = f(p.x[: w.size]) if callable(f) else f[: w.size] return np.sum(w * fx) # type: ignore[no-any-return] @@ -409,48 +432,8 @@ def _diff_caputo_l2(m: L2, f: ArrayOrScalarFunction, p: Points) -> Array: # FIXME: in the uniform case, we can also do an FFT, but we need different # weights for that, so we leave it like this for now for n in range(1, df.size): - w = quadrature_weights(m, p, n + 1) - df[n] = np.sum(w * fx[: n + 1]) - - return df - - -def _weights_l2(alpha: float, i: int | Array, k: int | Array) -> Array: - return np.array((i - k) ** (2 - alpha) - (i - k - 1) ** (2 - alpha)) - - -@diff.register(L2) -def _diff_l2_method(m: L2, f: ArrayOrScalarFunction, p: Points) -> Array: - # precompute variables - x = p.x - fx = f(x) if callable(f) else f - - alpha = m.alpha - w0 = 1 / math.gamma(3 - alpha) - - # NOTE: [Li2020] Section 4.2 - # NOTE: the method is not written as in [Li2020] and has several tweaks: - # * terms are written as `sum(w * f'')` instead of `sum(w * f)`, which - # makes it easier to express w - # * boundary terms are approximated with a biased stencil. - - df = np.empty_like(x) - df[0] = np.nan - - if isinstance(p, UniformPoints): - w0 = w0 / p.dx[0] ** alpha - k = np.arange(fx.size) - - ddf = np.zeros(fx.size - 1, dtype=fx.dtype) - ddf[:-1] = fx[2:] - 2 * fx[1:-1] + fx[:-2] - ddf[-1] = 2 * fx[-1] - 5 * fx[-2] + 4 * fx[-3] - fx[-4] - - for n in range(1, df.size): - df[n] = w0 * np.sum(_weights_l2(alpha, n, k[:n]) * ddf[:n]) - else: - raise NotImplementedError( - f"'{type(m).__name__}' not implemented for '{type(p).__name__}' grids" - ) + w = quadrature_weights(m, p, n) + df[n] = np.sum(w * fx[: w.size]) return df @@ -479,6 +462,10 @@ def __post_init__(self) -> None: ) +def _weights_l2(alpha: float, i: int | Array, k: int | Array) -> Array: + return np.array((i - k) ** (2 - alpha) - (i - k - 1) ** (2 - alpha)) + + @diff.register(L2C) def _diff_l2c_method(m: L2C, f: ArrayOrScalarFunction, p: Points) -> Array: # precompute variables diff --git a/src/pycaputo/grid.py b/src/pycaputo/grid.py index d6ffd8c..fe395ec 100644 --- a/src/pycaputo/grid.py +++ b/src/pycaputo/grid.py @@ -56,6 +56,11 @@ def xm(self) -> Array: """Array of midpoints.""" return np.array(self.x[1:] + self.x[:-1]) / 2 + @cached_property + def dxm(self) -> Array: + """Distance between midpoints.""" + return np.diff(self.xm) + def translate(self, a: float, b: float) -> Points: """Linearly translate the set of points to the new interval :math:`[a, b]`.""" # translate from [a', b'] to [0, 1] to [a, b] From ea94b7f70c465b59e57f441b581b433be71a174f Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Fri, 9 Aug 2024 19:58:32 +0300 Subject: [PATCH 16/26] feat: finish up the L2 method --- src/pycaputo/differentiation/caputo.py | 43 ++++++++++++++++++++------ tests/test_diff_caputo.py | 6 +++- 2 files changed, 39 insertions(+), 10 deletions(-) diff --git a/src/pycaputo/differentiation/caputo.py b/src/pycaputo/differentiation/caputo.py index a7fdb82..f5bcc41 100644 --- a/src/pycaputo/differentiation/caputo.py +++ b/src/pycaputo/differentiation/caputo.py @@ -99,11 +99,14 @@ def _quadrature_weights_caputo_l1(m: L1, p: Points, n: int) -> Array: if n == 0: return np.array([], dtype=p.dtype) - a = _caputo_piecewise_constant_integral(p, n, m.alpha) / p.dx[: n - 1] + a = _caputo_piecewise_constant_integral(p, n, m.alpha) + # NOTE: the first step of the discretization is just + # sum a_{ik} f'_k + # and we need to re-arrange the sum with the approximation of the derivative w = np.empty(n, dtype=p.dtype) w[n - 1] = 0.0 - w[: n - 1] = a + w[: n - 1] = a / p.dx[: n - 1] w[1:n] = w[: n - 1] - w[1:n] w[0] = -w[0] @@ -156,8 +159,11 @@ class L1D(CaputoMethod): of order :math:`\alpha \in (0, 1)`. This method is defined in Section 4.1.1 (II) from [Li2020]_ for general - non-uniform grids. This method is equivalent to :class:`L1`, but evaluation - requires explicit knowledge of the first derivative. + non-uniform grids. + + This method is equivalent to :class:`L1`, but evaluation requires explicit + knowledge of the first derivative, i.e. the *f* function must be a callable + satisfying :class:`~pycaputo.tuping.DifferentiableScalarFunction`. """ if __debug__: @@ -189,27 +195,46 @@ def _diffs_caputo_l1d(m: L1D, f: ArrayOrScalarFunction, p: Points, n: int) -> Sc if n == 0: return np.array([np.nan]) - if not isinstance(f, DifferentiableScalarFunction): + if not callable(f): raise FunctionCallableError( f"The '{type(m).__name__}' method requires evaluating the first " f"derivative. 'f' is not callable: {type(f)}" ) - w = quadrature_weights(m, p, n) - dfx = f((p.x[1 : n + 1] + p.x[:n]) / 2, d=1) + # FIXME: isinstance(f, DifferentiableScalarFunction) does not work? + assert isinstance(f, DifferentiableScalarFunction) + + try: + w = quadrature_weights(m, p, n) + dfx = f((p.x[1 : n + 1] + p.x[:n]) / 2, d=1) + except TypeError as exc: + raise FunctionCallableError( + f"The '{type(m).__name__}' method requires evaluating the first " + f"derivative. 'f' is not a 'DifferentiableScalarFunction'." + ) from exc return np.sum(w * dfx) # type: ignore[no-any-return] @diff.register(L1D) def _diff_caputo_l1d(m: L1D, f: ArrayOrScalarFunction, p: Points) -> Array: - if not isinstance(f, DifferentiableScalarFunction): + # FIXME: isinstance(f, DifferentiableScalarFunction) does not work? + if not callable(f): raise FunctionCallableError( f"The '{type(m).__name__}' method requires evaluating the first " f"derivative. 'f' is not callable: {type(f)}" ) - dfx = f((p.x[1:] + p.x[:-1]) / 2, d=1) + # FIXME: isinstance(f, DifferentiableScalarFunction) does not work? + assert isinstance(f, DifferentiableScalarFunction) + + try: + dfx = f((p.x[1:] + p.x[:-1]) / 2, d=1) + except TypeError as exc: + raise FunctionCallableError( + f"The '{type(m).__name__}' method requires evaluating the first " + f"derivative. 'f' is not a 'DifferentiableScalarFunction'." + ) from exc df = np.empty((p.size, *dfx.shape[1:]), dtype=dfx.dtype) df[0] = np.nan diff --git a/tests/test_diff_caputo.py b/tests/test_diff_caputo.py index f6983d3..8801bfd 100644 --- a/tests/test_diff_caputo.py +++ b/tests/test_diff_caputo.py @@ -63,6 +63,7 @@ def df_test(x: Array, *, alpha: float, mu: float = 3.5) -> Array: [ ("L1", "stretch"), ("L1", "uniform"), + ("L1D", "uniform"), ("L2", "uniform"), ("L2C", "uniform"), ("ModifiedL1", "midpoints"), @@ -93,6 +94,9 @@ def test_caputo_lmethods( if name == "L1": meth = caputo.L1(alpha=alpha) order = 2.0 - alpha + elif name == "L1D": + meth = caputo.L1D(alpha=alpha) + order = 2.0 - alpha elif name == "ModifiedL1": meth = caputo.ModifiedL1(alpha=alpha) order = 2 - alpha @@ -474,7 +478,7 @@ def test_caputo_vs_differint( else: # NOTE: we use slightly different boundary handling for the L2 methods # so these get larger errors compared to differint - assert error_vs_di < 1.0e-2 + assert error_vs_di < 1.1e-2 # }}} From 2909845954ad507a4e31e18bdff79327175dc2b8 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sun, 11 Aug 2024 20:22:15 +0300 Subject: [PATCH 17/26] feat(caputo): add some methods for testing --- src/pycaputo/differentiation/caputo.py | 312 ++++++++++++++++--------- tests/test_diff_caputo.py | 15 +- 2 files changed, 217 insertions(+), 110 deletions(-) diff --git a/src/pycaputo/differentiation/caputo.py b/src/pycaputo/differentiation/caputo.py index f5bcc41..2420380 100644 --- a/src/pycaputo/differentiation/caputo.py +++ b/src/pycaputo/differentiation/caputo.py @@ -150,105 +150,6 @@ def _diff_caputo_l1(m: L1, f: ArrayOrScalarFunction, p: Points) -> Array: # }}} -# {{{ L1D - - -@dataclass(frozen=True) -class L1D(CaputoMethod): - r"""Implements the L1 method for the Caputo fractional derivative - of order :math:`\alpha \in (0, 1)`. - - This method is defined in Section 4.1.1 (II) from [Li2020]_ for general - non-uniform grids. - - This method is equivalent to :class:`L1`, but evaluation requires explicit - knowledge of the first derivative, i.e. the *f* function must be a callable - satisfying :class:`~pycaputo.tuping.DifferentiableScalarFunction`. - """ - - if __debug__: - - def __post_init__(self) -> None: - super().__post_init__() - if not 0 < self.alpha < 1: - raise ValueError( - f"'{type(self).__name__}' only supports 0 < alpha < 1: {self.alpha}" - ) - - -@quadrature_weights.register(L1D) -def _quadrature_weights_caputo_l1d(m: L1D, p: Points, n: int) -> Array: - if not 0 <= n <= p.size: - raise IndexError(f"Index 'n' out of range: 0 <= {n} < {p.size}") - - if n == 0: - return np.array([]) - - return _caputo_piecewise_constant_integral(p, n + 1, m.alpha) - - -@diffs.register(L1D) -def _diffs_caputo_l1d(m: L1D, f: ArrayOrScalarFunction, p: Points, n: int) -> Scalar: - if not 0 <= n <= p.size: - raise IndexError(f"Index 'n' out of range: 0 <= {n} < {p.size}") - - if n == 0: - return np.array([np.nan]) - - if not callable(f): - raise FunctionCallableError( - f"The '{type(m).__name__}' method requires evaluating the first " - f"derivative. 'f' is not callable: {type(f)}" - ) - - # FIXME: isinstance(f, DifferentiableScalarFunction) does not work? - assert isinstance(f, DifferentiableScalarFunction) - - try: - w = quadrature_weights(m, p, n) - dfx = f((p.x[1 : n + 1] + p.x[:n]) / 2, d=1) - except TypeError as exc: - raise FunctionCallableError( - f"The '{type(m).__name__}' method requires evaluating the first " - f"derivative. 'f' is not a 'DifferentiableScalarFunction'." - ) from exc - - return np.sum(w * dfx) # type: ignore[no-any-return] - - -@diff.register(L1D) -def _diff_caputo_l1d(m: L1D, f: ArrayOrScalarFunction, p: Points) -> Array: - # FIXME: isinstance(f, DifferentiableScalarFunction) does not work? - if not callable(f): - raise FunctionCallableError( - f"The '{type(m).__name__}' method requires evaluating the first " - f"derivative. 'f' is not callable: {type(f)}" - ) - - # FIXME: isinstance(f, DifferentiableScalarFunction) does not work? - assert isinstance(f, DifferentiableScalarFunction) - - try: - dfx = f((p.x[1:] + p.x[:-1]) / 2, d=1) - except TypeError as exc: - raise FunctionCallableError( - f"The '{type(m).__name__}' method requires evaluating the first " - f"derivative. 'f' is not a 'DifferentiableScalarFunction'." - ) from exc - - df = np.empty((p.size, *dfx.shape[1:]), dtype=dfx.dtype) - df[0] = np.nan - - for n in range(1, df.size): - w = quadrature_weights(m, p, n) - df[n] = np.sum(w * dfx[:n]) - - return df - - -# }}} - - # {{{ ModifiedL1 @@ -332,8 +233,8 @@ class L2(CaputoMethod): r"""Implements the L2 method for the Caputo fractional derivative of order :math:`\alpha \in (1, 2)`. - This method is defined in Section 4.1.2 from [Li2020]_ for general - non-uniform grids. + This method is defined in Section 4.1.2 from [Li2020]_ for uniform grids. + The variant implemented here supports arbitrary non-uniform grids. .. note:: @@ -401,16 +302,14 @@ def _quadrature_weights_caputo_l2(m: L2, p: Points, n: int) -> Array: # weights have size at least 4 for the initial biased stencil. w = np.zeros(max(4, n + 1), dtype=p.dtype) - - # coefficients a = _caputo_piecewise_constant_integral(p, n + 1, m.alpha - 1) - c_l = _caputo_d2_boundary_coefficients(p.x, p.x[0]) - d_l, d_m, d_r = _caputo_d2_interior_coefficients(p, n) # add boundary stencil + c_l = _caputo_d2_boundary_coefficients(p.x, p.x[0]) w[:4] += a[0] * c_l # add interior stencils + d_l, d_m, d_r = _caputo_d2_interior_coefficients(p, n) if n > 1: w[0] += a[1] * d_l[0] w[1] += a[1] * d_m[0] @@ -428,9 +327,6 @@ def _quadrature_weights_caputo_l2(m: L2, p: Points, n: int) -> Array: @diffs.register(L2) def _diffs_caputo_l2(m: L2, f: ArrayOrScalarFunction, p: Points, n: int) -> Array: - if n < 0: - n = p.size + n - if not 0 <= n <= p.size: raise IndexError(f"Index 'n' out of range: 0 <= {n} < {p.size}") @@ -526,6 +422,206 @@ def _diff_l2c_method(m: L2C, f: ArrayOrScalarFunction, p: Points) -> Array: # }}} +# {{{ L2F + + +@dataclass(frozen=True) +class L2F(CaputoMethod): + r"""Implements the L2 method for the Caputo fractional derivative + of order :math:`\alpha \in (1, 2)`. + + This is similar to :class:`L2`, but it assumes that the function *f* can + be evaluated at :math:`f(x_{-1})` outside of the domain :math:`[a, b]`. + For symmetry, we set :math:`x_{-1} = a - \Delta x_0`. This allows using + the same centered stencil for all the points in the domain. + + This method mainly exists for comparison with the literature, as [Li2020]_ + also uses the value outside of the interval. + """ + + +@quadrature_weights.register(L2F) +def _quadrature_weights_caputo_l2f(m: L2F, p: Points, n: int) -> Array: + if not 0 <= n <= p.size: + raise IndexError(f"Index 'n' out of range: 0 <= {n} < {p.size}") + + if n == 0: + return np.array([], dtype=p.dtype) + + # weights have at least size 3 for the initial stencil + w = np.zeros(max(3, n + 2), dtype=p.dtype) + a = _caputo_piecewise_constant_integral(p, n + 1, m.alpha - 1) + + # add boundary stencil + dx2 = p.dx[0] ** 2 + c_l, c_m, c_r = 1.0 / dx2, -2.0 / dx2, 1.0 / dx2 + w[0] = a[0] * c_l + w[1] = a[0] * c_m + w[2] = a[0] * c_r + + # add interior stencils + d_l, d_m, d_r = _caputo_d2_interior_coefficients(p, n) + if n > 1: + w[1] += a[1] * d_l[0] + w[2] += a[1] * d_m[0] + w[n + 1] += a[-1] * d_r[-1] + + if n > 2: + w[2] += a[2] * d_l[1] + w[n] += a[-1] * d_m[-1] + a[-2] * d_r[-2] + + if n > 3: + w[3:n] += a[1:-2] * d_r[:-2] + a[2:-1] * d_m[1:-1] + a[3:] * d_l[2:] + + return w + + +@diffs.register(L2F) +def _diffs_caputo_l2f(m: L2F, f: ArrayOrScalarFunction, p: Points, n: int) -> Array: + if not callable(f): + raise FunctionCallableError( + f"The '{type(m).__name__}' method requires evaluating the function " + f"values. 'f' is not callable: {type(f)}" + ) + + if not 0 <= n <= p.size: + raise IndexError(f"Index 'n' out of range: 0 <= {n} < {p.size}") + + if n == 0: + return np.array([np.nan]) + + w = quadrature_weights(m, p, n) + x = np.empty(w.size, dtype=p.dtype) + x[0] = p.a - p.dx[0] + x[1:] = p.x + fx = f(x) + + return np.sum(w * fx) # type: ignore[no-any-return] + + +@diff.register(L2F) +def _diff_caputo_l2f(m: L2F, f: ArrayOrScalarFunction, p: Points) -> Array: + if not callable(f): + raise FunctionCallableError( + f"The '{type(m).__name__}' method requires evaluating the function " + f"values. 'f' is not callable: {type(f)}" + ) + + x = np.empty(p.size + 1, dtype=p.dtype) + x[0] = p.a - p.dx[0] + x[1:] = p.x + fx = f(x) + + df = np.empty(p.shape, dtype=fx.dtype) + df[0] = np.nan + + for n in range(1, df.size): + w = quadrature_weights(m, p, n) + df[n] = np.sum(w * fx[: w.size]) + + return df + + +# }}} + + +# {{{ LXD + + +@dataclass(frozen=True) +class LXD(CaputoMethod): + r"""Implements the LX method for the Caputo fractional derivative + of order :math:`\alpha \in (0, \infty)`. + + This method is equivalent to :class:`L1` (or :class:`L2` method), but + evaluation requires explicit knowledge of the required derivative, i.e. the + *f* function must be a callable satisfying + :class:`~pycaputo.tuping.DifferentiableScalarFunction`. With explicit + knowledge of the derivative, we can evaluate any fractional order. + + The derivatives are always evaluated at the midpoint of each interval. + """ + + if __debug__: + + def __post_init__(self) -> None: + super().__post_init__() + + +@quadrature_weights.register(LXD) +def _quadrature_weights_caputo_lxd(m: LXD, p: Points, n: int) -> Array: + if not 0 <= n <= p.size: + raise IndexError(f"Index 'n' out of range: 0 <= {n} < {p.size}") + + if n == 0: + return np.array([]) + + d = m.d + return _caputo_piecewise_constant_integral(p, n + 1, d.alpha - d.n + 1) + + +@diffs.register(LXD) +def _diffs_caputo_lxd(m: LXD, f: ArrayOrScalarFunction, p: Points, n: int) -> Scalar: + if not 0 <= n <= p.size: + raise IndexError(f"Index 'n' out of range: 0 <= {n} < {p.size}") + + if n == 0: + return np.array([np.nan]) + + if not callable(f): + raise FunctionCallableError( + f"The '{type(m).__name__}' method requires evaluating the first " + f"derivative. 'f' is not callable: {type(f)}" + ) + + # FIXME: isinstance(f, DifferentiableScalarFunction) does not work? + assert isinstance(f, DifferentiableScalarFunction) + d = m.d + + try: + w = quadrature_weights(m, p, n) + dfx = f((p.x[1 : n + 1] + p.x[:n]) / 2, d=d.n) + except TypeError as exc: + raise FunctionCallableError( + f"The '{type(m).__name__}' method requires evaluating the first " + f"derivative. 'f' is not a 'DifferentiableScalarFunction'." + ) from exc + + return np.sum(w * dfx) # type: ignore[no-any-return] + + +@diff.register(LXD) +def _diff_caputo_lxd(m: LXD, f: ArrayOrScalarFunction, p: Points) -> Array: + if not callable(f): + raise FunctionCallableError( + f"The '{type(m).__name__}' method requires evaluating the first " + f"derivative. 'f' is not callable: {type(f)}" + ) + + # FIXME: isinstance(f, DifferentiableScalarFunction) does not work? + assert isinstance(f, DifferentiableScalarFunction) + d = m.d + + try: + dfx = f((p.x[1:] + p.x[:-1]) / 2, d=d.n) + except TypeError as exc: + raise FunctionCallableError( + f"The '{type(m).__name__}' method requires evaluating the first " + f"derivative. 'f' is not a 'DifferentiableScalarFunction'." + ) from exc + + df = np.empty((p.size, *dfx.shape[1:]), dtype=dfx.dtype) + df[0] = np.nan + + for n in range(1, df.size): + w = quadrature_weights(m, p, n) + df[n] = np.sum(w * dfx[:n]) + + return df + + +# }}} + # {{{ SpectralJacobi diff --git a/tests/test_diff_caputo.py b/tests/test_diff_caputo.py index 8801bfd..638f775 100644 --- a/tests/test_diff_caputo.py +++ b/tests/test_diff_caputo.py @@ -29,6 +29,8 @@ def f_test(x: Array, d: int = 0, *, mu: float = 3.5) -> Array: return (0.5 - x) ** mu elif d == 1: return -mu * (0.5 - x) ** (mu - 1) + elif d == 2: + return mu * (mu - 1) * (0.5 - x) ** (mu - 2) else: raise NotImplementedError @@ -66,6 +68,8 @@ def df_test(x: Array, *, alpha: float, mu: float = 3.5) -> Array: ("L1D", "uniform"), ("L2", "uniform"), ("L2C", "uniform"), + ("L2D", "uniform"), + ("L2F", "uniform"), ("ModifiedL1", "midpoints"), ("ModifiedL1", "stretch"), ], @@ -85,7 +89,7 @@ def test_caputo_lmethods( from pycaputo.grid import make_points_from_name - if name in {"L2", "L2C"}: + if name in {"L2", "L2C", "L2D", "L2F"}: alpha += 1 from pycaputo.utils import EOCRecorder, savefig, stringify_eoc @@ -95,7 +99,7 @@ def test_caputo_lmethods( meth = caputo.L1(alpha=alpha) order = 2.0 - alpha elif name == "L1D": - meth = caputo.L1D(alpha=alpha) + meth = caputo.LXD(alpha=alpha) order = 2.0 - alpha elif name == "ModifiedL1": meth = caputo.ModifiedL1(alpha=alpha) @@ -107,6 +111,12 @@ def test_caputo_lmethods( elif name == "L2C": meth = caputo.L2C(alpha=alpha) order = 3.0 - alpha + elif name == "L2F": + meth = caputo.L2F(alpha=alpha) + order = 1.0 + elif name == "L2D": + meth = caputo.LXD(alpha=alpha) + order = 3.0 - alpha else: raise ValueError(f"Unsupported method: '{name}'") @@ -127,6 +137,7 @@ def test_caputo_lmethods( df_num = diff(meth, f_test, p) df_ref = df_test(p.x, alpha=alpha) + print(df_num[:10] - df_ref[:10]) h = np.max(p.dx) e = la.norm(df_num[1:] - df_ref[1:]) / la.norm(df_ref[1:]) From 338bcfedc7dc3511d332730c4b6236d0930f0b8f Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Mon, 12 Aug 2024 10:17:06 +0300 Subject: [PATCH 18/26] docs: fix some typos --- src/pycaputo/differentiation/__init__.py | 2 +- src/pycaputo/differentiation/caputo.py | 16 ++++++++-------- src/pycaputo/lagrange.py | 5 +++-- 3 files changed, 12 insertions(+), 11 deletions(-) diff --git a/src/pycaputo/differentiation/__init__.py b/src/pycaputo/differentiation/__init__.py index f809c17..f9aff7d 100644 --- a/src/pycaputo/differentiation/__init__.py +++ b/src/pycaputo/differentiation/__init__.py @@ -109,7 +109,7 @@ def diffs_fallback( .. warning:: - Falling back to the :func:`~pycaputo.differentiation.diff`` function + Falling back to the :func:`~pycaputo.differentiation.diff` function will be significantly slower, since all the points on the grid *p* are evaluated. Use this function with care. """ diff --git a/src/pycaputo/differentiation/caputo.py b/src/pycaputo/differentiation/caputo.py index 2420380..b865968 100644 --- a/src/pycaputo/differentiation/caputo.py +++ b/src/pycaputo/differentiation/caputo.py @@ -161,7 +161,7 @@ class ModifiedL1(CaputoMethod): This method is defined in Section 4.1.1 (III) from [Li2020]_. Note that this method evaluates the fractional derivative at the midpoints :math:`(x_i + x_{i - 1}) / 2` for any given input grid except - :class:`~pycaputo.grids.MidPoints`. + :class:`~pycaputo.grid.MidPoints`. As noted in [Li2020]_, this method has the same order of convergence as the standard L1 method. However the weights can be used to construct a @@ -186,8 +186,6 @@ def _quadrature_weights_caputo_modified_l1(m: ModifiedL1, p: Points, n: int) -> p = make_midpoints_from(p) w = quadrature_weights.dispatch(L1)(m, p, n) - # FIXME: not clear if this does anything? Might be the same as just - # doing L1 on the original grid.. wp = np.empty((n + 1,), dtype=w.dtype) wp[0] = w[0] / 2.0 wp[-1] = w[-1] / 2.0 @@ -238,7 +236,7 @@ class L2(CaputoMethod): .. note:: - Unlike the method from [Li2020_], we do not assume knowledge of points + Unlike the method from [Li2020]_, we do not assume knowledge of points outside of the domain. Instead a biased stencil is used at the boundary. """ @@ -430,13 +428,14 @@ class L2F(CaputoMethod): r"""Implements the L2 method for the Caputo fractional derivative of order :math:`\alpha \in (1, 2)`. - This is similar to :class:`L2`, but it assumes that the function *f* can - be evaluated at :math:`f(x_{-1})` outside of the domain :math:`[a, b]`. + This is similar to :class:`L2`, but it assumes that *f* is a callable that can + be evaluated at :math:`f(x_{-1})`, i.e. outside of the domain :math:`[a, b]`. For symmetry, we set :math:`x_{-1} = a - \Delta x_0`. This allows using the same centered stencil for all the points in the domain. This method mainly exists for comparison with the literature, as [Li2020]_ - also uses the value outside of the interval. + also uses the value outside of the interval. Note that this is not always + possible, e.g. for :math:`\sqrt{x}` on :math:`[0, 1]`. """ @@ -536,7 +535,7 @@ class LXD(CaputoMethod): This method is equivalent to :class:`L1` (or :class:`L2` method), but evaluation requires explicit knowledge of the required derivative, i.e. the *f* function must be a callable satisfying - :class:`~pycaputo.tuping.DifferentiableScalarFunction`. With explicit + :class:`~pycaputo.typing.DifferentiableScalarFunction`. With explicit knowledge of the derivative, we can evaluate any fractional order. The derivatives are always evaluated at the midpoint of each interval. @@ -622,6 +621,7 @@ def _diff_caputo_lxd(m: LXD, f: ArrayOrScalarFunction, p: Points) -> Array: # }}} + # {{{ SpectralJacobi diff --git a/src/pycaputo/lagrange.py b/src/pycaputo/lagrange.py index 430880c..ddec7d3 100644 --- a/src/pycaputo/lagrange.py +++ b/src/pycaputo/lagrange.py @@ -76,7 +76,8 @@ def lagrange_riemann_liouville_integral( .. math:: - I^\alpha[\phi_{kj}](x_n) = \frac{1}{\Gamma(\alpha)} \int_{0}^{x_n} + I^\alpha[\phi_{kj}](x_n) = + \frac{1}{\Gamma(\alpha)} \int_{0}^{x_n} (x_n - s)^{\alpha - 1} \phi_{kj}(s) \,\mathrm{d}s. As the polynomials are zero except in the domain of definition, this simplifies @@ -133,7 +134,7 @@ def lagrange_caputo_derivative( .. math:: - D^{q, \alpha}_{nk} \triangleq + L^{q, \alpha}_{nk} \triangleq \frac{1}{\Gamma(m - \alpha)} \int_{x_k}^{x_{k + 1}} (x_n - s)^{m - \alpha - 1} \frac{\mathrm{d}^m}{\mathrm{d} s^m} From 5543d8d7fd3bbbb7615b9d6de49b72bbea3a1e5d Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Mon, 12 Aug 2024 10:23:14 +0300 Subject: [PATCH 19/26] docs: add changelog --- docs/changelog.rst | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/docs/changelog.rst b/docs/changelog.rst index 6b2a5b5..bf5d599 100644 --- a/docs/changelog.rst +++ b/docs/changelog.rst @@ -1,6 +1,29 @@ Changelog ========= +pycaputo (TBD) +-------------- + +Features +^^^^^^^^ + +* Add a :mod:`pycaputo.typing` module containing some helpful typing definitions + (mostly previously in :mod:`pycaputo.utils`). +* Reworked :func:`~functools.singledispatch` functions for + :class:`~pycaputo.differentiation.DerivativeMethod`. We now have + :func:`~pycaputo.differentiation.quadrature_weights`, + :func:`~pycaputo.differentiation.differentiation_matrix`, + :func:`~pycaputo.differentiation.diffs`, and + :func:`~pycaputo.differentiation.diff`. Not all of these need to be implemented + and most methods are not ported yet. Currently only the + :class:`~pycaputo.differentiation.caputo.L1` and + :class:`~pycaputo.differentiation.caputo.L2` methods implement these new functions. +* Introduce some more specific methods for the Caputo derivative. The + :class:`~pycaputo.differentiation.caputo.L2F` uses the L2 method with function + evaluations outside of the interval of definition. The + :class:`~pycaputo.differentiation.caputo.LXD` allows evaluating arbitrary + Caputo derivatives when the integer derivatives are known. + pycaputo 0.7.0 (July 13, 2024) ------------------------------ From 1329bc8dda972d60fec78cceb2e74d7c136606bd Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Mon, 12 Aug 2024 10:43:39 +0300 Subject: [PATCH 20/26] feat(diff): rename SpectralJacobi to Jacobi --- docs/changelog.rst | 6 +++ docs/misc_others.rst | 4 +- examples/caputo-derivative-quadratic.py | 4 +- examples/caputo-gradient-spectral.py | 2 +- src/pycaputo/differentiation/__init__.py | 2 +- src/pycaputo/differentiation/caputo.py | 30 +++++++------ src/pycaputo/quadrature/riemann_liouville.py | 16 ++++--- src/pycaputo/typing.py | 47 ++++++++++++-------- tests/test_diff_caputo.py | 2 +- 9 files changed, 69 insertions(+), 44 deletions(-) diff --git a/docs/changelog.rst b/docs/changelog.rst index bf5d599..5461b56 100644 --- a/docs/changelog.rst +++ b/docs/changelog.rst @@ -24,6 +24,12 @@ Features :class:`~pycaputo.differentiation.caputo.LXD` allows evaluating arbitrary Caputo derivatives when the integer derivatives are known. +Changes +^^^^^^^ + +* Renamed `pycaputo.differentiation.caputo.SpectralJacobi` to + :class:`~pycaputo.differentiation.caputo.Jacobi`. + pycaputo 0.7.0 (July 13, 2024) ------------------------------ diff --git a/docs/misc_others.rst b/docs/misc_others.rst index 1d70b66..da6f6e7 100644 --- a/docs/misc_others.rst +++ b/docs/misc_others.rst @@ -14,14 +14,14 @@ Logging Typing ------ -.. automodule:: pycaputo.typing - .. class:: pycaputo.utils.P A generic invarint :class:`typing.ParamSpec`. alias of ParamSpec('P') +.. automodule:: pycaputo.typing + Utils ----- diff --git a/examples/caputo-derivative-quadratic.py b/examples/caputo-derivative-quadratic.py index 21160cd..9959640 100644 --- a/examples/caputo-derivative-quadratic.py +++ b/examples/caputo-derivative-quadratic.py @@ -13,7 +13,7 @@ import numpy as np from pycaputo.differentiation import diff -from pycaputo.differentiation.caputo import SpectralJacobi +from pycaputo.differentiation.caputo import Jacobi from pycaputo.grid import make_jacobi_gauss_lobatto_points from pycaputo.typing import Array @@ -58,7 +58,7 @@ def f(x: Array) -> Array: p = make_jacobi_gauss_lobatto_points(32, a=0.0, b=2.0) for alpha in alphas: - method = SpectralJacobi(alpha=alpha) + method = Jacobi(alpha=alpha) df_num = diff(method, f, p) ax.plot(p.x, df_num, color="k", alpha=0.2) diff --git a/examples/caputo-gradient-spectral.py b/examples/caputo-gradient-spectral.py index 0cc9a04..b005f13 100644 --- a/examples/caputo-gradient-spectral.py +++ b/examples/caputo-gradient-spectral.py @@ -53,7 +53,7 @@ def df(x: Array, alpha: float) -> Array: from pycaputo.differentiation import caputo alpha = 0.9 -method = caputo.SpectralJacobi(alpha=0.9) +method = caputo.Jacobi(alpha=0.9) # make a set of points at which to evaluate the gradient rng = np.random.default_rng(42) diff --git a/src/pycaputo/differentiation/__init__.py b/src/pycaputo/differentiation/__init__.py index f9aff7d..38ffed9 100644 --- a/src/pycaputo/differentiation/__init__.py +++ b/src/pycaputo/differentiation/__init__.py @@ -45,7 +45,7 @@ def guess_method_for_order( if isinstance(d, CaputoDerivative): if isinstance(p, grid.JacobiGaussLobattoPoints): - m = caputo.SpectralJacobi(d.alpha) + m = caputo.Jacobi(d.alpha) elif 0 < d.alpha < 1: if isinstance(p, grid.MidPoints): m = caputo.ModifiedL1(d.alpha) diff --git a/src/pycaputo/differentiation/caputo.py b/src/pycaputo/differentiation/caputo.py index b865968..3c30adc 100644 --- a/src/pycaputo/differentiation/caputo.py +++ b/src/pycaputo/differentiation/caputo.py @@ -581,9 +581,9 @@ def _diffs_caputo_lxd(m: LXD, f: ArrayOrScalarFunction, p: Points, n: int) -> Sc w = quadrature_weights(m, p, n) dfx = f((p.x[1 : n + 1] + p.x[:n]) / 2, d=d.n) except TypeError as exc: - raise FunctionCallableError( - f"The '{type(m).__name__}' method requires evaluating the first " - f"derivative. 'f' is not a 'DifferentiableScalarFunction'." + raise TypeError( + f"{type(m).__name__!r} requires a 'DifferentiableScalarFunction': " + f"f is a {type(f).__name__!r}" ) from exc return np.sum(w * dfx) # type: ignore[no-any-return] @@ -604,9 +604,9 @@ def _diff_caputo_lxd(m: LXD, f: ArrayOrScalarFunction, p: Points) -> Array: try: dfx = f((p.x[1:] + p.x[:-1]) / 2, d=d.n) except TypeError as exc: - raise FunctionCallableError( - f"The '{type(m).__name__}' method requires evaluating the first " - f"derivative. 'f' is not a 'DifferentiableScalarFunction'." + raise TypeError( + f"{type(m).__name__!r} requires a 'DifferentiableScalarFunction': " + f"f is a {type(f).__name__!r}" ) from exc df = np.empty((p.size, *dfx.shape[1:]), dtype=dfx.dtype) @@ -622,11 +622,11 @@ def _diff_caputo_lxd(m: LXD, f: ArrayOrScalarFunction, p: Points) -> Array: # }}} -# {{{ SpectralJacobi +# {{{ Jacobi @dataclass(frozen=True) -class SpectralJacobi(CaputoMethod): +class Jacobi(CaputoMethod): r"""Caputo derivative approximation using spectral methods based on Jacobi polynomials. @@ -650,8 +650,8 @@ class SpectralJacobi(CaputoMethod): """ -@diff.register(SpectralJacobi) -def _diff_jacobi(m: SpectralJacobi, f: ArrayOrScalarFunction, p: Points) -> Array: +@diff.register(Jacobi) +def _diff_jacobi(m: Jacobi, f: ArrayOrScalarFunction, p: Points) -> Array: from pycaputo.grid import JacobiGaussLobattoPoints if not isinstance(p, JacobiGaussLobattoPoints): @@ -817,14 +817,18 @@ def _diff_caputo_yuan_agrawal( f: ArrayOrScalarFunction, p: Points, ) -> Array: - if not isinstance(f, DifferentiableScalarFunction): + # FIXME: isinstance(f, DifferentiableScalarFunction) does not work? + assert isinstance(f, DifferentiableScalarFunction) + + try: + dtype = np.array(f(p.x[0], d=0)).dtype + except TypeError: raise TypeError( f"{type(m).__name__!r} requires a 'DifferentiableScalarFunction': " f"f is a {type(f).__name__!r}" - ) + ) from None x = p.x - dtype = np.array(f(p.x[0], d=0)).dtype # solve ODE at quadrature nodes omega, w = m.nodes_and_weights() diff --git a/src/pycaputo/quadrature/riemann_liouville.py b/src/pycaputo/quadrature/riemann_liouville.py index cbae88a..5b0feb2 100644 --- a/src/pycaputo/quadrature/riemann_liouville.py +++ b/src/pycaputo/quadrature/riemann_liouville.py @@ -271,14 +271,20 @@ def _quad_rl_cubic_hermite( ) -> Array: from pycaputo.grid import UniformPoints - if not isinstance(f, DifferentiableScalarFunction): - raise TypeError(f"Input 'f' needs to be a callable: {type(f).__name__}") - if not isinstance(p, UniformPoints): raise TypeError(f"Only uniform points are supported: {type(p).__name__}") - fx = f(p.x, d=0) - fp = f(p.x, d=1) + # FIXME: isinstance(f, DifferentiableScalarFunction) does not work? + assert isinstance(f, DifferentiableScalarFunction) + + try: + fx = f(p.x, d=0) + fp = f(p.x, d=1) + except TypeError as exc: + raise TypeError( + f"{type(m).__name__!r} requires a 'DifferentiableScalarFunction': " + f"f is a {type(f).__name__!r}" + ) from exc alpha = -m.alpha h = p.dx[0] diff --git a/src/pycaputo/typing.py b/src/pycaputo/typing.py index f02cf6e..8e7be01 100644 --- a/src/pycaputo/typing.py +++ b/src/pycaputo/typing.py @@ -26,17 +26,22 @@ # {{{ TypeVars +# NOTE: sphinx doesn't seem to render this correctly at the moment, so it's +# written explicitly in `misc_others.rst` P = ParamSpec("P") T = TypeVar("T") """A generic invariant :class:`typing.TypeVar`.""" R = TypeVar("R") """A generic invariant :class:`typing.TypeVar`.""" + PathLike = Union[pathlib.Path, str] """A union of types supported as paths.""" + # }}} + # {{{ numpy @@ -52,6 +57,7 @@ # }}} + # {{{ dataclass @@ -68,6 +74,7 @@ class DataclassInstance(Protocol): # {{{ callable protocols +@runtime_checkable class ScalarFunction(Protocol): """A generic callable that can be evaluated at :math:`x`. @@ -80,6 +87,26 @@ def __call__(self, x: Array, /) -> Array: """ +@runtime_checkable +class DifferentiableScalarFunction(Protocol): + """A :class:`ScalarFunction` that can also compute its integer order derivatives. + + .. automethod:: __call__ + """ + + def __call__(self, x: Array, /, d: int = 0) -> Array: + """Evaluate the function or any of its derivatives. + + :arg x: a scalar or array at which to evaluate the function. + :arg d: order of the derivative. + """ + + +ArrayOrScalarFunction = Union[Array, ScalarFunction, DifferentiableScalarFunction] +"""A union of scalar functions.""" + + +@runtime_checkable class StateFunction(Protocol): r"""A generic callable for right-hand side functions :math:`\mathbf{f}(t, \mathbf{y})`. @@ -94,6 +121,7 @@ def __call__(self, t: float, y: Array, /) -> Array: """ +@runtime_checkable class ScalarStateFunction(Protocol): """A generic callable similar to :class:`StateFunction` that returns a scalar. @@ -108,27 +136,8 @@ def __call__(self, t: float, y: Array, /) -> float: """ -@runtime_checkable -class DifferentiableScalarFunction(Protocol): - """A :class:`ScalarFunction` that can also compute its integer order derivatives. - - By default no derivatives are implemented, so subclasses can handle any such - cases. - """ - - def __call__(self, x: Array, /, d: int = 0) -> Array: - """Evaluate the function or any of its derivatives. - - :arg x: a scalar or array at which to evaluate the function. - :arg d: order of the derivative. - """ - - StateFunctionT = TypeVar("StateFunctionT", bound=StateFunction) """An invariant :class:`~typing.TypeVar` bound to :class:`StateFunction`.""" -ArrayOrScalarFunction = Union[Array, ScalarFunction, DifferentiableScalarFunction] -"""A union of scalar functions.""" - # }}} diff --git a/tests/test_diff_caputo.py b/tests/test_diff_caputo.py index 638f775..182f110 100644 --- a/tests/test_diff_caputo.py +++ b/tests/test_diff_caputo.py @@ -198,7 +198,7 @@ def test_caputo_spectral( from pycaputo.grid import make_jacobi_gauss_lobatto_points from pycaputo.utils import EOCRecorder, savefig - meth = caputo.SpectralJacobi(alpha=alpha) + meth = caputo.Jacobi(alpha=alpha) eoc = EOCRecorder() if visualize: From a41ea6ef8430f11ce0efacdb3c2b074356070fa9 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Mon, 12 Aug 2024 11:35:51 +0300 Subject: [PATCH 21/26] feat(diff): expand RL methods to match Caputo --- .../differentiation/riemann_liouville.py | 117 ++++++++++++++---- src/pycaputo/finite_difference.py | 2 +- 2 files changed, 93 insertions(+), 26 deletions(-) diff --git a/src/pycaputo/differentiation/riemann_liouville.py b/src/pycaputo/differentiation/riemann_liouville.py index c0d339c..afe6e7c 100644 --- a/src/pycaputo/differentiation/riemann_liouville.py +++ b/src/pycaputo/differentiation/riemann_liouville.py @@ -7,13 +7,15 @@ from dataclasses import dataclass from functools import cached_property +import numpy as np + from pycaputo.derivatives import RiemannLiouvilleDerivative, Side from pycaputo.grid import Points from pycaputo.logging import get_logger -from pycaputo.typing import Array, ArrayOrScalarFunction +from pycaputo.typing import Array, ArrayOrScalarFunction, Scalar from . import caputo -from .base import DerivativeMethod, diff +from .base import DerivativeMethod, diff, diffs, quadrature_weights logger = get_logger(__name__) @@ -96,6 +98,93 @@ def base(self) -> DerivativeMethod: return caputo.L2C(self.alpha) +def _rl_d1_boundary_coefficients(x: Array) -> Array: + # NOTE: this is a 3rd order approximation of the derivative at f(a) + # that works on non-uniform grids too (derived by Mathematica) + + x0, x1, x2, x3 = x[:4] + c0 = 1.0 / (x0 - x1) + 1.0 / (x0 - x2) + 1.0 / (x0 - x3) + c1 = -(x0 - x2) * (x0 - x3) / ((x0 - x1) * (x1 - x2) * (x1 - x3)) + c2 = -(x0 - x1) * (x0 - x3) / ((x0 - x2) * (x2 - x1) * (x2 - x3)) + c3 = -(x0 - x1) * (x0 - x2) / ((x0 - x3) * (x3 - x1) * (x3 - x2)) + + return np.array([c0, c1, c2, c3]) + + +def _rl_correction_weights( + m: RiemannLiouvilleFromCaputoMethod, + x: Array, + n: int, +) -> Array: + xn = x[n] + alpha = m.alpha + w0 = 1.0 / (xn - x[0]) ** alpha / math.gamma(1 - alpha) + + if 0.0 < alpha < 1.0: + w = np.array([w0], dtype=x.dtype) + elif 1.0 < alpha < 2.0: + w = ( + _rl_d1_boundary_coefficients(x) + / (xn - x[0]) ** (alpha - 1) + / math.gamma(2 - alpha) + ) + w[0] += w0 + else: + raise NotImplementedError + + return w + + +def _rl_correction(m: RiemannLiouvilleFromCaputoMethod, x: Array, fx: Array) -> Array: + alpha = m.alpha + df = fx[0] / (x[1:] - x[0]) ** alpha / math.gamma(1 - alpha) + + if 0.0 < alpha <= 1.0: + # NOTE: handled above + pass + elif 1.0 < alpha <= 2.0: + dfx = _rl_d1_boundary_coefficients(x) @ fx[:4] + df += dfx / (x[1:] - x[0]) ** (alpha - 1) / math.gamma(2 - alpha) + else: + raise NotImplementedError(f"Unsupported derivative order: {alpha}") + + return np.array(df) + + +@quadrature_weights.register(L1) +@quadrature_weights.register(L2) +@quadrature_weights.register(L2C) +def _quadrature_weights_rl_from_caputo( + m: RiemannLiouvilleFromCaputoMethod, + p: Points, + n: int, +) -> Array: + w = quadrature_weights(m.base, p, n) + + wc = _rl_correction_weights(m, p.x, n) + w[: wc.size] += wc + + return w + + +@diffs.register(L1) +@diffs.register(L2) +@diffs.register(L2C) +def _diffs_rl_from_caputo( + m: RiemannLiouvilleFromCaputoMethod, + f: ArrayOrScalarFunction, + p: Points, + n: int, +) -> Scalar: + x = p.x + fx = f(x[: n + 1]) if callable(f) else f[: n + 1] + + df = diffs(m.base, fx, p, n) + + wc = _rl_correction_weights(m, p.x, n) + return df + wc @ fx[: wc.size] # type: ignore[no-any-return] + + @diff.register(L1) @diff.register(L2) @diff.register(L2C) @@ -104,33 +193,11 @@ def _diff_rl_from_caputo( f: ArrayOrScalarFunction, p: Points, ) -> Array: - alpha = m.alpha x = p.x fx = f(x) if callable(f) else f df = diff(m.base, fx, p) - df[1:] += fx[0] / (x[1:] - x[0]) ** alpha / math.gamma(1 - alpha) - - if 0.0 < alpha <= 1.0: - # NOTE: handled above - pass - elif 1.0 < alpha <= 2.0: - # NOTE: this is a 3rd order approximation of the derivative at f(a) - # that works on non-uniform grids too (derived by Mathematica) - dfx = ( - # fmt: off - (fx[0] - fx[1]) / (x[0] - x[1]) - + (fx[0] - fx[2]) / (x[0] - x[2]) - + (fx[0] - fx[3]) / (x[0] - x[3]) - - (fx[1] - fx[2]) / (x[1] - x[2]) - - (fx[1] - fx[3]) * (x[0] - x[2]) / ((x[1] - x[2]) * (x[1] - x[3])) - + (fx[2] - fx[3]) * (x[0] - x[1]) / ((x[1] - x[2]) * (x[2] - x[3])) - # fmt: on - ) - - df[1:] += dfx / (x[1:] - x[0]) ** (alpha - 1) / math.gamma(2 - alpha) - else: - raise NotImplementedError(f"Unsupported derivative order: {alpha}") + df[1:] += _rl_correction(m, p.x, fx) return df diff --git a/src/pycaputo/finite_difference.py b/src/pycaputo/finite_difference.py index a05ec4c..d54727b 100644 --- a/src/pycaputo/finite_difference.py +++ b/src/pycaputo/finite_difference.py @@ -1,4 +1,4 @@ -# SPDX-FileCopyrightText: 2023 Alexandru Fikl +# q SPDX-FileCopyrightText: 2023 Alexandru Fikl # SPDX-License-Identifier: MIT from __future__ import annotations From a5a5c5119f2279d49d84021d3e8d003f2f3e0abb Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Mon, 12 Aug 2024 19:33:45 +0300 Subject: [PATCH 22/26] test: check diff interface consistency --- src/pycaputo/differentiation/base.py | 4 +- src/pycaputo/differentiation/caputo.py | 23 ++- tests/test_diff_caputo.py | 261 +++++++++++++++++++------ 3 files changed, 211 insertions(+), 77 deletions(-) diff --git a/src/pycaputo/differentiation/base.py b/src/pycaputo/differentiation/base.py index 811eebd..930b264 100644 --- a/src/pycaputo/differentiation/base.py +++ b/src/pycaputo/differentiation/base.py @@ -89,7 +89,9 @@ def differentiation_matrix(m: DerivativeMethod, p: Points) -> Array: """ n = p.size W = np.zeros((n, n), dtype=p.dtype) - W[0, 0] = np.nan + + # NOTE: make it very obvious that the first row is garbage + W[0, :] = np.nan try: for i in range(1, W.shape[0]): diff --git a/src/pycaputo/differentiation/caputo.py b/src/pycaputo/differentiation/caputo.py index 3c30adc..3bb8216 100644 --- a/src/pycaputo/differentiation/caputo.py +++ b/src/pycaputo/differentiation/caputo.py @@ -93,36 +93,35 @@ def _caputo_piecewise_constant_integral(p: Points, n: int, alpha: float) -> Arra @quadrature_weights.register(L1) def _quadrature_weights_caputo_l1(m: L1, p: Points, n: int) -> Array: - if not 0 <= n <= p.size: + if not 0 <= n < p.size: raise IndexError(f"Index 'n' out of range: 0 <= {n} < {p.size}") if n == 0: return np.array([], dtype=p.dtype) - a = _caputo_piecewise_constant_integral(p, n, m.alpha) + a = _caputo_piecewise_constant_integral(p, n + 1, m.alpha) / p.dx[:n] # NOTE: the first step of the discretization is just # sum a_{ik} f'_k # and we need to re-arrange the sum with the approximation of the derivative - w = np.empty(n, dtype=p.dtype) - w[n - 1] = 0.0 - w[: n - 1] = a / p.dx[: n - 1] - w[1:n] = w[: n - 1] - w[1:n] - w[0] = -w[0] + w = np.empty(n + 1, dtype=p.dtype) + w[1:n] = a[:-1] - a[1:] + w[0] = -a[0] + w[n] = a[-1] return w @diffs.register(L1) def _diffs_caputo_l1(m: L1, f: ArrayOrScalarFunction, p: Points, n: int) -> Scalar: - if not 0 <= n <= p.size: + if not 0 <= n < p.size: raise IndexError(f"Index 'n' out of range: 0 <= {n} < {p.size}") if n == 0: return np.array([np.nan]) - w = quadrature_weights(m, p, n + 1) - fx = f(p.x[: n + 1]) if callable(f) else f[: n + 1] + w = quadrature_weights(m, p, n) + fx = f(p.x[: w.size]) if callable(f) else f[: w.size] return np.sum(w * fx) # type: ignore[no-any-return] @@ -141,8 +140,8 @@ def _diff_caputo_l1(m: L1, f: ArrayOrScalarFunction, p: Points) -> Array: # FIXME: in the uniform case, we can also do an FFT, but we need different # weights for that, so we leave it like this for now for n in range(1, df.size): - w = quadrature_weights(m, p, n + 1) - df[n] = np.sum(w * fx[: n + 1]) + w = quadrature_weights(m, p, n) + df[n] = np.sum(w * fx[: w.size]) return df diff --git a/tests/test_diff_caputo.py b/tests/test_diff_caputo.py index 182f110..8d997f9 100644 --- a/tests/test_diff_caputo.py +++ b/tests/test_diff_caputo.py @@ -21,7 +21,7 @@ set_recommended_matplotlib() -# {{{ test_caputo_lmethods +# {{{ utils def f_test(x: Array, d: int = 0, *, mu: float = 3.5) -> Array: @@ -60,40 +60,36 @@ def df_test(x: Array, *, alpha: float, mu: float = 3.5) -> Array: raise ValueError(f"Unsupported order: {alpha}") -@pytest.mark.parametrize( - ("name", "grid_type"), - [ - ("L1", "stretch"), - ("L1", "uniform"), - ("L1D", "uniform"), - ("L2", "uniform"), - ("L2C", "uniform"), - ("L2D", "uniform"), - ("L2F", "uniform"), - ("ModifiedL1", "midpoints"), - ("ModifiedL1", "stretch"), - ], -) -@pytest.mark.parametrize("alpha", [0.1, 0.25, 0.5, 0.75, 0.9]) -def test_caputo_lmethods( +def make_test_points_from_name( name: str, - grid_type: str, - alpha: float, + n: int, *, - visualize: bool = False, -) -> None: - r""" - Test convergence of the LXXX methods for the Caputo derivative. - The convergence is checked in the :math:`\ell^2` norm using :func:`f_test`. - """ + j_alpha: float = 0.0, + j_beta: float = 0.0, +) -> Points: + from pycaputo.grid import make_jacobi_gauss_lobatto_points, make_points_from_name - from pycaputo.grid import make_points_from_name + p: Points + if name == "jacobi": + p = make_jacobi_gauss_lobatto_points( + n, + a=0.0, + b=0.5, + alpha=j_alpha, + beta=j_beta, + ) + else: + p = make_points_from_name(name, n, a=0.0, b=0.5) - if name in {"L2", "L2C", "L2D", "L2F"}: - alpha += 1 + return p - from pycaputo.utils import EOCRecorder, savefig, stringify_eoc +def make_method_from_name( + name: str, + alpha: float, + *, + quad_order: int | None = None, +) -> tuple[caputo.CaputoMethod, float]: meth: caputo.CaputoMethod if name == "L1": meth = caputo.L1(alpha=alpha) @@ -117,9 +113,71 @@ def test_caputo_lmethods( elif name == "L2D": meth = caputo.LXD(alpha=alpha) order = 3.0 - alpha + elif name == "Jacobi": + meth = caputo.Jacobi(alpha=alpha) + order = 5.0 + elif name == "YuanAgrawal": + assert quad_order is not None + meth = caputo.YuanAgrawal( + alpha, + quad_order=quad_order, + method="Radau", + ) + beta = 2.0 * alpha - 1.0 + order = 1.0 - beta + elif name == "Diethelm": + assert quad_order is not None + meth = caputo.Diethelm(alpha, quad_order=quad_order, method="Radau") + order = 0 + elif name == "BirkSong": + assert quad_order is not None + meth = caputo.BirkSong(alpha, quad_order=quad_order, method="Radau") + order = 0 else: raise ValueError(f"Unsupported method: '{name}'") + return meth, order + + +# }}} + + +# {{{ test_caputo_lmethods + + +@pytest.mark.parametrize( + ("name", "grid_type"), + [ + ("L1", "stretch"), + ("L1", "uniform"), + ("L1D", "uniform"), + ("L2", "uniform"), + ("L2C", "uniform"), + ("L2D", "uniform"), + ("L2F", "uniform"), + ("ModifiedL1", "midpoints"), + ("ModifiedL1", "stretch"), + ], +) +@pytest.mark.parametrize("alpha", [0.1, 0.25, 0.5, 0.75, 0.9]) +def test_caputo_lmethods( + name: str, + grid_type: str, + alpha: float, + *, + visualize: bool = False, +) -> None: + r""" + Test convergence of the LXXX methods for the Caputo derivative. + The convergence is checked in the :math:`\ell^2` norm using :func:`f_test`. + """ + + if name in {"L2", "L2C", "L2D", "L2F"}: + alpha += 1 + + from pycaputo.utils import EOCRecorder, savefig, stringify_eoc + + meth, order = make_method_from_name(name, alpha) eoc = EOCRecorder(order=order) if visualize: @@ -129,15 +187,16 @@ def test_caputo_lmethods( ax = fig.gca() for n in [16, 32, 64, 128, 256, 512, 768, 1024]: - p = make_points_from_name(grid_type, n, a=0.0, b=0.5) + p = make_test_points_from_name(grid_type, n) if name == "ModifiedL1" and grid_type != "midpoints": from pycaputo.grid import make_midpoints_from - p = make_midpoints_from(p) + p_ref: Points = make_midpoints_from(p) + else: + p_ref = p df_num = diff(meth, f_test, p) - df_ref = df_test(p.x, alpha=alpha) - print(df_num[:10] - df_ref[:10]) + df_ref = df_test(p_ref.x, alpha=alpha) h = np.max(p.dx) e = la.norm(df_num[1:] - df_ref[1:]) / la.norm(df_ref[1:]) @@ -145,13 +204,13 @@ def test_caputo_lmethods( logger.info("n %4d h %.5e e %.12e", n, h, e) if visualize: - ax.plot(p.x[1:], df_num[1:]) - # ax.semilogy(p.x, abs(df_num - df_ref)) + ax.plot(p_ref.x[1:], df_num[1:]) + # ax.semilogy(p_ref.x, abs(df_num - df_ref)) logger.info("\n%s", stringify_eoc(eoc)) if visualize: - ax.plot(p.x[1:], df_ref[1:], "k--") + ax.plot(p_ref.x[1:], df_ref[1:], "k--") ax.set_xlabel("$x$") ax.set_ylabel(rf"$D^{{{alpha}}}_C f$") # ax.set_ylim([1.0e-16, 1]) @@ -195,10 +254,9 @@ def test_caputo_spectral( it requires a slightly different setup. """ - from pycaputo.grid import make_jacobi_gauss_lobatto_points from pycaputo.utils import EOCRecorder, savefig - meth = caputo.Jacobi(alpha=alpha) + meth, _ = make_method_from_name("Jacobi", alpha) eoc = EOCRecorder() if visualize: @@ -208,13 +266,7 @@ def test_caputo_spectral( ax = fig.gca() for n in [8, 12, 16, 24, 32]: - p = make_jacobi_gauss_lobatto_points( - n, - a=0.0, - b=0.5, - alpha=j_alpha, - beta=j_beta, - ) + p = make_test_points_from_name("jacobi", n, j_alpha=j_alpha, j_beta=j_beta) df_num = diff(meth, f_test, p) df_ref = df_test(p.x, alpha=alpha) @@ -274,7 +326,6 @@ def test_caputo_diffusive( from pycaputo.grid import make_points_from_name from pycaputo.utils import EOCRecorder, savefig - meth: caputo.DiffusiveCaputoMethod resolutions = [8, 16, 24, 32, 48, 64] n = 128 @@ -289,23 +340,7 @@ def test_caputo_diffusive( ax = fig.gca() for quad_order in resolutions: - if name == "YuanAgrawal": - meth = caputo.YuanAgrawal( - alpha, - quad_order=quad_order, - method="Radau", - ) - beta = 2.0 * alpha - 1.0 - order = 1.0 - beta - elif name == "Diethelm": - meth = caputo.Diethelm(alpha, quad_order=quad_order, method="Radau") - order = None - elif name == "BirkSong": - meth = caputo.BirkSong(alpha, quad_order=quad_order, method="Radau") - order = None - else: - raise ValueError(f"Unsupported method: '{name}'") - + meth, order = make_method_from_name(name, alpha, quad_order=quad_order) df_num = diff(meth, f_test, p) h = 1.0 / quad_order @@ -329,7 +364,7 @@ def test_caputo_diffusive( filename = f"test_caputo_{meth.name}_{alpha}".replace(".", "_") savefig(fig, dirname / filename.lower()) - if order is not None: + if order > 0: assert order - 0.25 < eoc.estimated_order < order + 1.0 else: assert eoc.estimated_order > 1.9 @@ -494,6 +529,104 @@ def test_caputo_vs_differint( # }}} + +# {{{ test_caputo_consistency + + +@pytest.mark.parametrize( + ("name", "grid_type"), + [ + ("L1", "stretch"), + ("L2", "stretch"), + ], +) +@pytest.mark.parametrize("alpha", [0.3, 0.9]) +def test_caputo_consistency( + name: str, + grid_type: str, + alpha: float, + *, + visualize: bool = False, +) -> None: + """ + Tests that the `quadrature_weights`, `differentiation_matrix`, `diffs` and + `diff` methods actually give the same results. + """ + from pycaputo.differentiation import ( + differentiation_matrix, + diffs, + quadrature_weights, + ) + + if name in {"L2", "L2F", "L2D"}: + alpha += 1 + + meth, _ = make_method_from_name(name, alpha) + p = make_test_points_from_name(grid_type, 256) + + # check quadrature_weights vs differentiation_matrix + W_weights = np.zeros((p.size, p.size), dtype=p.dtype) + for n in range(1, p.size): + w = quadrature_weights(meth, p, n) + W_weights[n, : w.size] = w + + W_mat = differentiation_matrix(meth, p) + + error = la.norm(W_weights[1:, :] - W_mat[1:, :]) + logger.info("quadrature_weights vs differentiation_matrix") + logger.info(" error %.12e", error) + assert error < 1.0e-15 + + # check quadrature_weights vs diffs + logger.info("quadrature_weights vs diffs") + fx = f_test(p.x) + for n in [0, 1, p.size // 2 + 3, p.size - 1]: + df_from_diffs = diffs(meth, fx, p, n) + + if n == 0: + w = quadrature_weights(meth, p, n) + assert w.shape == (0,) + assert df_from_diffs.shape == (1,) + assert np.isnan(df_from_diffs) + else: + w = quadrature_weights(meth, p, n) + df_from_weights = w @ fx[: w.size] + logger.info( + " value %.12e %.12e (ref %.12e) w shape %s", + df_from_weights, + df_from_diffs, + df_test(p.x[n], alpha=alpha), + w.shape, + ) + + error = la.norm(df_from_weights - df_from_diffs) + logger.info(" n %d error %.12e", n, error) + assert error < 6.0e-14 + + # check diff vs differentiation_matrix + df = diff(meth, fx, p) + df_mat = W_mat @ fx + + error = la.norm(df[1:] - df_mat[1:]) / la.norm(df_mat[1:]) + logger.info("diff vs differentiation_matrix") + logger.info(" error %.12e", error) + assert error < 1.0e-12 + + # check diff vs diffs + logger.info("diff vs diffs") + df_ref = df + assert isinstance(df_ref, np.ndarray) + for n in [1, p.size // 2 - 5, p.size - 1]: + df_from_diffs = diffs(meth, fx, p, n) + + error = la.norm(df_from_diffs - df_ref[n]) / la.norm(df_ref[n]) + logger.info(" error %.12e", error) + assert error < 1.0e-15 + + +# }}} + + if __name__ == "__main__": import sys From a37460cbb62f7b9b5b60219dba546efe265b0912 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Mon, 12 Aug 2024 20:08:16 +0300 Subject: [PATCH 23/26] feat: speed up caputo derivative evaluation --- src/pycaputo/differentiation/caputo.py | 207 ++++++++++++++++--------- tests/test_diff_caputo.py | 6 +- 2 files changed, 138 insertions(+), 75 deletions(-) diff --git a/src/pycaputo/differentiation/caputo.py b/src/pycaputo/differentiation/caputo.py index 3bb8216..d13a9a7 100644 --- a/src/pycaputo/differentiation/caputo.py +++ b/src/pycaputo/differentiation/caputo.py @@ -9,7 +9,7 @@ import numpy as np -from pycaputo.derivatives import CaputoDerivative +from pycaputo.derivatives import CaputoDerivative, Side from pycaputo.grid import MidPoints, Points, UniformPoints, make_midpoints_from from pycaputo.logging import get_logger from pycaputo.typing import ( @@ -70,7 +70,7 @@ def __post_init__(self) -> None: ) -def _caputo_piecewise_constant_integral(p: Points, n: int, alpha: float) -> Array: +def _caputo_piecewise_constant_integral(x: Array, alpha: float) -> Array: r"""Computes the integrals .. math:: @@ -82,24 +82,17 @@ def _caputo_piecewise_constant_integral(p: Points, n: int, alpha: float) -> Arra for :math:`0 \le k \le n - 1`. """ - x = p.x[:n] xn = x[-1] - return np.array( ((xn - x[:-1]) ** (1 - alpha) - (xn - x[1:]) ** (1 - alpha)) / math.gamma(2 - alpha) ) -@quadrature_weights.register(L1) -def _quadrature_weights_caputo_l1(m: L1, p: Points, n: int) -> Array: - if not 0 <= n < p.size: - raise IndexError(f"Index 'n' out of range: 0 <= {n} < {p.size}") - - if n == 0: - return np.array([], dtype=p.dtype) - - a = _caputo_piecewise_constant_integral(p, n + 1, m.alpha) / p.dx[:n] +def _caputo_l1_weights(p: Points, n: int, alpha: float) -> Array: + x = p.x[: n + 1] + dx = p.dx[:n] + a = _caputo_piecewise_constant_integral(x, alpha) / dx # NOTE: the first step of the discretization is just # sum a_{ik} f'_k @@ -112,6 +105,17 @@ def _quadrature_weights_caputo_l1(m: L1, p: Points, n: int) -> Array: return w +@quadrature_weights.register(L1) +def _quadrature_weights_caputo_l1(m: L1, p: Points, n: int) -> Array: + if not 0 <= n < p.size: + raise IndexError(f"Index 'n' out of range: 0 <= {n} < {p.size}") + + if n == 0: + return np.array([], dtype=p.dtype) + + return _caputo_l1_weights(p, n, m.alpha) + + @diffs.register(L1) def _diffs_caputo_l1(m: L1, f: ArrayOrScalarFunction, p: Points, n: int) -> Scalar: if not 0 <= n < p.size: @@ -120,7 +124,7 @@ def _diffs_caputo_l1(m: L1, f: ArrayOrScalarFunction, p: Points, n: int) -> Scal if n == 0: return np.array([np.nan]) - w = quadrature_weights(m, p, n) + w = _caputo_l1_weights(p, n, m.alpha) fx = f(p.x[: w.size]) if callable(f) else f[: w.size] return np.sum(w * fx) # type: ignore[no-any-return] @@ -134,13 +138,15 @@ def _diff_caputo_l1(m: L1, f: ArrayOrScalarFunction, p: Points) -> Array: f"Shape of 'f' does match points: got {fx.shape} expected {p.shape}" ) - df = np.empty_like(fx) - df[0] = np.nan + alpha = m.alpha # FIXME: in the uniform case, we can also do an FFT, but we need different # weights for that, so we leave it like this for now + df = np.empty_like(fx) + df[0] = np.nan + for n in range(1, df.size): - w = quadrature_weights(m, p, n) + w = _caputo_l1_weights(p, n, alpha) df[n] = np.sum(w * fx[: w.size]) return df @@ -268,45 +274,24 @@ def _caputo_d2_boundary_coefficients(x: Array, s: Scalar) -> Array: return np.array([c0, c1, c2, c3]) -def _caputo_d2_interior_coefficients(p: Points, n: int) -> tuple[Array, Array, Array]: - r"""Get coefficients for the second-order approximation of the second - derivative at the interior. - - .. math:: - - \frac{\partial f}{\partial x}(s) = - d_- f_{k - 1} + d_m f_k + d_+ f_{k + 1} - - which will simplify to the standard :math:`(1, -2, 1)` in the uniform case. - """ - +def _caputo_l2_weights(p: Points, n: int, alpha: float) -> Array: + x = p.x[: n + 1] dx = p.dx[:n] dxm = p.dxm[: n - 1] - d_l = 1.0 / (dx[:-1] * dxm) - d_r = 1.0 / (dx[1:] * dxm) - - return d_l, -(d_l + d_r), d_r - - -@quadrature_weights.register(L2) -def _quadrature_weights_caputo_l2(m: L2, p: Points, n: int) -> Array: - if not 0 <= n <= p.size: - raise IndexError(f"Index 'n' out of range: 0 <= {n} < {p.size}") - - if n == 0: - return np.array([], dtype=p.dtype) - - # weights have size at least 4 for the initial biased stencil. - w = np.zeros(max(4, n + 1), dtype=p.dtype) - a = _caputo_piecewise_constant_integral(p, n + 1, m.alpha - 1) + # weights have size at least 4 for the initial biased stencil + w = np.zeros(max(4, n + 1), dtype=x.dtype) + a = _caputo_piecewise_constant_integral(x, alpha - 1) # add boundary stencil - c_l = _caputo_d2_boundary_coefficients(p.x, p.x[0]) + c_l = _caputo_d2_boundary_coefficients(p.x[:4], x[0]) w[:4] += a[0] * c_l # add interior stencils - d_l, d_m, d_r = _caputo_d2_interior_coefficients(p, n) + d_l = 1.0 / (dx[:-1] * dxm) + d_r = 1.0 / (dx[1:] * dxm) + d_m = -(d_l + d_r) + if n > 1: w[0] += a[1] * d_l[0] w[1] += a[1] * d_m[0] @@ -322,6 +307,17 @@ def _quadrature_weights_caputo_l2(m: L2, p: Points, n: int) -> Array: return w +@quadrature_weights.register(L2) +def _quadrature_weights_caputo_l2(m: L2, p: Points, n: int) -> Array: + if not 0 <= n <= p.size: + raise IndexError(f"Index 'n' out of range: 0 <= {n} < {p.size}") + + if n == 0: + return np.array([], dtype=p.dtype) + + return _caputo_l2_weights(p, n, m.alpha) + + @diffs.register(L2) def _diffs_caputo_l2(m: L2, f: ArrayOrScalarFunction, p: Points, n: int) -> Array: if not 0 <= n <= p.size: @@ -344,14 +340,33 @@ def _diff_caputo_l2(m: L2, f: ArrayOrScalarFunction, p: Points) -> Array: f"Shape of 'f' does match points: got {fx.shape} expected {p.shape}" ) - df = np.empty_like(fx) - df[0] = np.nan + # variables + x = p.x + dx = p.dx + dxm = p.dxm + alpha = m.alpha + + # estimate second-order derivative + # NOTE: this is faster than using the quadrature weights + ddfx = np.empty(p.size - 1, dtype=fx.dtype) + + # interior + cm = 1.0 / (dx[:-1] * dxm) + cp = 1.0 / (dx[1:] * dxm) + ddfx[1:] = cm * fx[:-2] - (cm + cp) * fx[1:-1] + cp * fx[2:] + + # boundary + c = _caputo_d2_boundary_coefficients(x[:4], x[0]) + ddfx[0] = c @ fx[: c.size] # FIXME: in the uniform case, we can also do an FFT, but we need different # weights for that, so we leave it like this for now + df = np.empty_like(fx) + df[0] = np.nan + for n in range(1, df.size): - w = quadrature_weights(m, p, n) - df[n] = np.sum(w * fx[: w.size]) + a = _caputo_piecewise_constant_integral(x[: n + 1], alpha - 1) + df[n] = np.sum(a * ddfx[: a.size]) return df @@ -437,28 +452,33 @@ class L2F(CaputoMethod): possible, e.g. for :math:`\sqrt{x}` on :math:`[0, 1]`. """ + side: Side -@quadrature_weights.register(L2F) -def _quadrature_weights_caputo_l2f(m: L2F, p: Points, n: int) -> Array: - if not 0 <= n <= p.size: - raise IndexError(f"Index 'n' out of range: 0 <= {n} < {p.size}") - - if n == 0: - return np.array([], dtype=p.dtype) +def _caputo_l2f_weights(p: Points, n: int, alpha: float) -> Array: # weights have at least size 3 for the initial stencil w = np.zeros(max(3, n + 2), dtype=p.dtype) - a = _caputo_piecewise_constant_integral(p, n + 1, m.alpha - 1) + + # variables + x = p.x[: n + 1] + dx = p.dx[:n] + dxm = p.dxm[: n - 1] + + a = _caputo_piecewise_constant_integral(x, alpha - 1) # add boundary stencil dx2 = p.dx[0] ** 2 c_l, c_m, c_r = 1.0 / dx2, -2.0 / dx2, 1.0 / dx2 + w[0] = a[0] * c_l w[1] = a[0] * c_m w[2] = a[0] * c_r # add interior stencils - d_l, d_m, d_r = _caputo_d2_interior_coefficients(p, n) + d_l = 1.0 / (dx[:-1] * dxm) + d_r = 1.0 / (dx[1:] * dxm) + d_m = -(d_l + d_r) + if n > 1: w[1] += a[1] * d_l[0] w[2] += a[1] * d_m[0] @@ -474,8 +494,25 @@ def _quadrature_weights_caputo_l2f(m: L2F, p: Points, n: int) -> Array: return w +@quadrature_weights.register(L2F) +def _quadrature_weights_caputo_l2f(m: L2F, p: Points, n: int) -> Array: + if m.side == Side.Right: + raise NotImplementedError("Right boundary approximation not implemented") + + if not 0 <= n <= p.size: + raise IndexError(f"Index 'n' out of range: 0 <= {n} < {p.size}") + + if n == 0: + return np.array([], dtype=p.dtype) + + return _caputo_l2f_weights(p, n, m.alpha) + + @diffs.register(L2F) def _diffs_caputo_l2f(m: L2F, f: ArrayOrScalarFunction, p: Points, n: int) -> Array: + if m.side == Side.Right: + raise NotImplementedError("Right boundary approximation not implemented") + if not callable(f): raise FunctionCallableError( f"The '{type(m).__name__}' method requires evaluating the function " @@ -488,7 +525,7 @@ def _diffs_caputo_l2f(m: L2F, f: ArrayOrScalarFunction, p: Points, n: int) -> Ar if n == 0: return np.array([np.nan]) - w = quadrature_weights(m, p, n) + w = _caputo_l2f_weights(p, n, m.alpha) x = np.empty(w.size, dtype=p.dtype) x[0] = p.a - p.dx[0] x[1:] = p.x @@ -505,17 +542,38 @@ def _diff_caputo_l2f(m: L2F, f: ArrayOrScalarFunction, p: Points) -> Array: f"values. 'f' is not callable: {type(f)}" ) - x = np.empty(p.size + 1, dtype=p.dtype) - x[0] = p.a - p.dx[0] - x[1:] = p.x - fx = f(x) + # variables + x = p.x + dx = p.dx + dxm = p.dxm + alpha = m.alpha + + fx = f(p.x) + + # estimate second-order derivative + # NOTE: this is faster than using the quadrature weights + ddfx = np.empty(p.size - 1, dtype=fx.dtype) + + cm = 1.0 / (dx[:-1] * dxm) + cp = 1.0 / (dx[1:] * dxm) - df = np.empty(p.shape, dtype=fx.dtype) + if m.side == Side.Left: + fa = f(p.a - dx[0]) + ddfx[1:] = cm * fx[:-2] - (cm + cp) * fx[1:-1] + cp * fx[2:] + ddfx[0] = (fx[1] - 2.0 * fx[0] + fa) / dx[0] ** 2 + else: + fb = f(p.b + dx[-1]) + ddfx[:-1] = cm * fx[:-2] - (cm + cp) * fx[1:-1] + cp * fx[2:] + ddfx[-1] = (fx[-2] - 2.0 * fx[-1] + fb) / dx[-1] ** 2 + + # FIXME: in the uniform case, we can also do an FFT, but we need different + # weights for that, so we leave it like this for now + df = np.empty_like(fx) df[0] = np.nan for n in range(1, df.size): - w = quadrature_weights(m, p, n) - df[n] = np.sum(w * fx[: w.size]) + a = _caputo_piecewise_constant_integral(x[: n + 1], alpha - 1) + df[n] = np.sum(a * ddfx[: a.size]) return df @@ -555,7 +613,8 @@ def _quadrature_weights_caputo_lxd(m: LXD, p: Points, n: int) -> Array: return np.array([]) d = m.d - return _caputo_piecewise_constant_integral(p, n + 1, d.alpha - d.n + 1) + x = p.x[: n + 1] + return _caputo_piecewise_constant_integral(x, d.alpha - d.n + 1) @diffs.register(LXD) @@ -598,10 +657,12 @@ def _diff_caputo_lxd(m: LXD, f: ArrayOrScalarFunction, p: Points) -> Array: # FIXME: isinstance(f, DifferentiableScalarFunction) does not work? assert isinstance(f, DifferentiableScalarFunction) - d = m.d + alpha_hat = m.alpha - m.d.n + 1 + x = p.x + xm = p.xm try: - dfx = f((p.x[1:] + p.x[:-1]) / 2, d=d.n) + dfx = f(xm, d=m.d.n) except TypeError as exc: raise TypeError( f"{type(m).__name__!r} requires a 'DifferentiableScalarFunction': " @@ -612,7 +673,7 @@ def _diff_caputo_lxd(m: LXD, f: ArrayOrScalarFunction, p: Points) -> Array: df[0] = np.nan for n in range(1, df.size): - w = quadrature_weights(m, p, n) + w = _caputo_piecewise_constant_integral(x[: n + 1], alpha_hat) df[n] = np.sum(w * dfx[:n]) return df diff --git a/tests/test_diff_caputo.py b/tests/test_diff_caputo.py index 8d997f9..0594110 100644 --- a/tests/test_diff_caputo.py +++ b/tests/test_diff_caputo.py @@ -108,7 +108,9 @@ def make_method_from_name( meth = caputo.L2C(alpha=alpha) order = 3.0 - alpha elif name == "L2F": - meth = caputo.L2F(alpha=alpha) + from pycaputo.derivatives import Side + + meth = caputo.L2F(alpha=alpha, side=Side.Left) order = 1.0 elif name == "L2D": meth = caputo.LXD(alpha=alpha) @@ -621,7 +623,7 @@ def test_caputo_consistency( error = la.norm(df_from_diffs - df_ref[n]) / la.norm(df_ref[n]) logger.info(" error %.12e", error) - assert error < 1.0e-15 + assert error < 6.0e-12 # }}} From dc8899e7816ce7aa328d8a3325b8eb9a49ff00d3 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Mon, 12 Aug 2024 21:03:12 +0300 Subject: [PATCH 24/26] test: fix differint comparison for L2 method --- pyproject.toml | 1 + src/pycaputo/differentiation/caputo.py | 9 ++++----- tests/test_diff_caputo.py | 26 ++++++++++++++++---------- 3 files changed, 21 insertions(+), 15 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 8c02aba..900b4be 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -47,6 +47,7 @@ dependencies = [ ] optional-dependencies.dev = [ "differint", + "differint @ git+https://github.com/alexfikl/differint.git@fix-l2-index#egg=differint", "doc8", "matplotlib", "mypy", diff --git a/src/pycaputo/differentiation/caputo.py b/src/pycaputo/differentiation/caputo.py index d13a9a7..392375c 100644 --- a/src/pycaputo/differentiation/caputo.py +++ b/src/pycaputo/differentiation/caputo.py @@ -547,23 +547,22 @@ def _diff_caputo_l2f(m: L2F, f: ArrayOrScalarFunction, p: Points) -> Array: dx = p.dx dxm = p.dxm alpha = m.alpha - fx = f(p.x) # estimate second-order derivative # NOTE: this is faster than using the quadrature weights ddfx = np.empty(p.size - 1, dtype=fx.dtype) - cm = 1.0 / (dx[:-1] * dxm) - cp = 1.0 / (dx[1:] * dxm) + cm = 1.0 / dx[:-1] + cp = 1.0 / dx[1:] if m.side == Side.Left: fa = f(p.a - dx[0]) - ddfx[1:] = cm * fx[:-2] - (cm + cp) * fx[1:-1] + cp * fx[2:] + ddfx[1:] = (cm * fx[:-2] - (cm + cp) * fx[1:-1] + cp * fx[2:]) / dxm ddfx[0] = (fx[1] - 2.0 * fx[0] + fa) / dx[0] ** 2 else: fb = f(p.b + dx[-1]) - ddfx[:-1] = cm * fx[:-2] - (cm + cp) * fx[1:-1] + cp * fx[2:] + ddfx[:-1] = (cm * fx[:-2] - (cm + cp) * fx[1:-1] + cp * fx[2:]) / dxm ddfx[-1] = (fx[-2] - 2.0 * fx[-1] + fb) / dx[-1] ** 2 # FIXME: in the uniform case, we can also do an FFT, but we need different diff --git a/tests/test_diff_caputo.py b/tests/test_diff_caputo.py index 0594110..967e4f5 100644 --- a/tests/test_diff_caputo.py +++ b/tests/test_diff_caputo.py @@ -471,11 +471,14 @@ def test_caputo_vs_differint( if name in {"L2", "L2C"}: alpha += 1 + meth: caputo.CaputoMethod if name == "L1": - meth: caputo.CaputoMethod = caputo.L1(alpha=alpha) + meth = caputo.L1(alpha=alpha) differint_meth: caputo.CaputoMethod = DifferIntCaputoL1(alpha=alpha) elif name == "L2": - meth = caputo.L2(alpha=alpha) + from pycaputo.derivatives import Side + + meth = caputo.L2F(alpha=alpha, side=Side.Right) differint_meth = DifferIntCaputoL2(alpha=alpha) elif name == "L2C": meth = caputo.L2C(alpha=alpha) @@ -485,17 +488,19 @@ def test_caputo_vs_differint( from pycaputo.grid import make_points_from_name - p = make_points_from_name("uniform", 512, a=0.0, b=0.5) + # NOTE: b must be < 0.5 so that f(p.b + dx) is well defined for the test function + p = make_points_from_name("uniform", 512, a=0.0, b=0.25) df_ref = df_test(p.x, alpha=alpha) df_num = diff(meth, f_test, p) df_num_di = diff(differint_meth, f_test, p) error_vs_ref = la.norm(df_num[1:] - df_ref[1:]) / la.norm(df_ref[1:]) - error_di_vs_ref = la.norm(df_num_di[1:-1] - df_ref[1:-1]) / la.norm(df_ref[1:-1]) - error_vs_di = la.norm(df_num[4:-4] - df_num_di[4:-4]) / la.norm(df_num_di[4:-4]) + error_di_vs_ref = la.norm(df_num_di[1:] - df_ref[1:]) / la.norm(df_ref[1:]) + error_vs_di = la.norm(df_num[1:] - df_num_di[1:]) / la.norm(df_num_di[1:]) + logger.info( - "error: vs ref %.12e vs differint %.12e differint vs ref %.12e", + "error: vs ref %.12e vs differint %.12e (differint vs ref %.12e)", error_vs_ref, error_vs_di, error_di_vs_ref, @@ -523,10 +528,11 @@ def test_caputo_vs_differint( assert error_vs_ref < 1.0e-2 if name == "L1": assert error_vs_di < 1.0e-12 - else: - # NOTE: we use slightly different boundary handling for the L2 methods - # so these get larger errors compared to differint - assert error_vs_di < 1.1e-2 + elif name == "L2": + assert error_vs_di < 1.0e-11 + elif name == "L2C": + # FIXME: need to implement a L2CF similar to L2F to better compare + assert error_vs_di < 7.0e-2 # }}} From a8a973f328cc5ff20c7362f9a030553df07f4b13 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Mon, 12 Aug 2024 21:03:49 +0300 Subject: [PATCH 25/26] chore: bump dependencies --- requirements-dev.txt | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/requirements-dev.txt b/requirements-dev.txt index 17ffd0f..463a0b9 100644 --- a/requirements-dev.txt +++ b/requirements-dev.txt @@ -34,7 +34,7 @@ contourpy==1.2.1 # via matplotlib cycler==0.12.1 # via matplotlib -differint==1.0.0 +differint @ git+https://github.com/alexfikl/differint.git@ab981f6bcd9435a46d85470a1ab01cbcf536837d#egg=differint # via pycaputo (pyproject.toml) doc8==1.1.1 # via pycaputo (pyproject.toml) @@ -183,9 +183,9 @@ typos==1.23.6 # via pycaputo (pyproject.toml) urllib3==2.2.2 # via requests -uv==0.2.34 +uv==0.2.35 # via pycaputo (pyproject.toml) -zipp==3.19.2 ; python_version < '3.10' +zipp==3.20.0 ; python_version < '3.10' # via # importlib-metadata # importlib-resources From fc8d407e72ab3963af9447bd6c20edb564303be8 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Mon, 12 Aug 2024 21:05:57 +0300 Subject: [PATCH 26/26] chore: allow hatch direct references for differint --- pyproject.toml | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 900b4be..125a274 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -46,7 +46,6 @@ dependencies = [ "typing-extensions; python_version<'3.10'", ] optional-dependencies.dev = [ - "differint", "differint @ git+https://github.com/alexfikl/differint.git@fix-l2-index#egg=differint", "doc8", "matplotlib", @@ -70,6 +69,9 @@ optional-dependencies.vis = [ urls.Documentation = "https://pycaputo.readthedocs.io" urls.Repository = "https://github.com/alexfikl/pycaputo" +[tool.hatch.metadata] +allow-direct-references = true + [tool.hatch.build.targets.sdist] exclude = [ ".github",