docs: improve docs for new diff functionality

docs/guide_differentiation.rst

@@ -1,17 +1,17 @@
 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
 
  df = diff(f, p, alpha)
 
 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

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)

src/pycaputo/__init__.py

@@ -29,18 +29,17 @@ 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")
 
  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,18 +55,17 @@ 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")
 
  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
 

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

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))