From 8567c11ca76ce8a8853bf27cb0ebc31b7691febb Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Wed, 7 Aug 2024 17:11:58 +0300 Subject: [PATCH] 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 febef25..07ed972 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 a7eb120..08a79d7 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 7d38663..15a3757 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.utils 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 a9144e7..58b71b2 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))