Rework derivative discretizations

docs/changelog.rst

@@ -1,6 +1,35 @@
 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.
+
+Changes
+^^^^^^^
+
+* Renamed `pycaputo.differentiation.caputo.SpectralJacobi` to
+  :class:`~pycaputo.differentiation.caputo.Jacobi`.
+
 pycaputo 0.7.0 (July 13, 2024)
 ------------------------------
 

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,37 @@ 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.differentiation_matrix`: a function that
+  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.
+These methods are all registered with the :func:`~functools.singledispatch` mechanism.
+An example setup is provided below.
 
 .. literalinclude:: ../examples/guide-differentiation.py
     :language: python

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

examples/caputo-derivative-l1-spectrum.py

@@ -0,0 +1,85 @@
+# SPDX-FileCopyrightText: 2023 Alexandru Fikl <alexfikl@gmail.com>
+# 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.typing 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$")
+
+# }}}

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)

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)

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,42 @@ def d(self) -> RiemannLiouvilleDerivative:
 # {{{
 
 # [register-start]
-from pycaputo.differentiation import diff
+from pycaputo.differentiation import (
+    diff,
+    differentiation_matrix,
+    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)
+
+
+@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,
+    n: int,
+) -> Scalar:
+    # ... add an actual implementation here ...
+    return np.array(0.0)
 
 
 @diff.register(RiemannLiouvilleDerivativeMethod)

pyproject.toml

@@ -46,7 +46,7 @@ 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",
     "mypy",
@@ -69,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",

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

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