feat(diff): adapt ModifiedL1 to new interface

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)

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.utils 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)
 
 
 # }}}