fode: add estimate for the Lipschitz constant

docs/fode.rst

@@ -45,3 +45,7 @@ Caputo Derivative FODEs
 .. autoclass:: CaputoPECEMethod
 .. autoclass:: CaputoPECMethod
 .. autoclass:: CaputoModifiedPECEMethod
+
+.. class:: Array
+
+   See :class:`pycaputo.utils.Array`.

docs/misc.rst

@@ -21,6 +21,11 @@ Generating Functions
 
 .. automodule:: pycaputo.generating_functions
 
+Lipschitz Constant Estimation
+-----------------------------
+
+.. automodule:: pycaputo.lipschitz
+
 Logging
 -------
 

docs/references.rst

@@ -6,6 +6,11 @@ References
     SIAM Journal on Mathematical Analysis, Vol. 17, pp. 704--719, 1986,
     `DOI <https://doi.org/10.1137/0517050>`__.
 
+.. [Wood1996] G. R. Wood, B. P. Zhang,
+    *Estimation of the Lipschitz Constant of a Function*,
+    Journal of Global Optimization, Vol. 8, 1996,
+    `DOI <https://doi.org/10.1007/bf00229304>`__.
+
 .. [Diethelm2002] K. Diethelm, N. J. Ford, A. D. Freed,
     *A Predictor-Corrector Approach for the Numerical Solution of
     Fractional Differential Equations*,

examples/example-lipschitz-weibull.py

@@ -0,0 +1,60 @@
+# SPDX-FileCopyrightText: 2023 Alexandru Fikl <alexfikl@gmail.com>
+# SPDX-License-Identifier: MIT
+
+import numpy as np
+import scipy.stats as ss
+
+from pycaputo.lipschitz import fit_reverse_weibull, uniform_sample_maximum_slopes
+from pycaputo.utils import Array, figure, set_recommended_matplotlib
+
+
+def f(x: Array) -> Array:
+    return np.array(x - x**3 / 3)
+
+
+# parameters
+
+a, b = -1.0, 1.0
+c, loc, scale = 1.0, 1.0, 0.002
+
+nslopes = 11
+nbatches = 100
+delta = 0.05
+
+# fit distribution
+
+rng = np.random.default_rng(seed=42)
+smax = uniform_sample_maximum_slopes(
+    f,
+    a,
+    b,
+    delta=delta,
+    nslopes=nslopes,
+    nbatches=nbatches,
+    rng=rng,
+)
+smax = np.sort(smax)
+
+cdf_opt = fit_reverse_weibull(smax)
+cdf_ref = ss.weibull_max(c=c, loc=loc, scale=scale)
+cdf_empirical = ss.ecdf(smax).cdf
+
+print(
+    "Optimal paramters: c {c:.12e} loc {loc:.12e} scale {scale:.12e}".format(
+        **cdf_opt.kwds
+    )
+)
+
+# plot
+
+set_recommended_matplotlib()
+
+with figure("example-lipschitz-weibull") as fig:
+    ax = fig.gca()
+
+    ax.step(smax, cdf_opt.cdf(smax), "k--", label="$Optimal$")
+    ax.step(smax, cdf_ref.cdf(smax), "k:", label="$Example$")
+    ax.step(smax, cdf_empirical.probabilities, label="$Empirical$")
+    ax.set_xlabel("$s_{max}$")
+
+    ax.legend()

pycaputo/fode/base.py

@@ -16,6 +16,9 @@
 logger = get_logger(__name__)
 
 
+# {{{ time step
+
+
 def make_predict_time_step_fixed(dt: float) -> ScalarStateFunction:
     """
     :returns: a callable returning a fixed time step *dt*.
@@ -73,6 +76,67 @@ def predict_time_step(t: float, y: Array) -> float:
     return predict_time_step
 
 
+def make_predict_time_step_lipschitz(
+    f: StateFunction,
+    alpha: float,
+    yspan: tuple[float, float],
+    *,
+    theta: float = 0.9,
+) -> ScalarStateFunction:
+    r"""Estimate a time step based on the Lipschitz constant.
+
+    This method estimates the time step using
+
+    .. math::
+
+        \Delta t = \theta \frac{1}{(\Gamma(2 - \alpha) L)^{\frac{1}{\alpha}}},
+
+    where :math:`L` is the Lipschitz constant estimated by
+    :func:`estimate_lipschitz_constant` and :math:`\theta` is a coefficient that
+    should be smaller than 1.
+
+    .. warning::
+
+        This method requires many evaluations of the right-hand side function
+        *f* and will not be efficient in practice. Ideally, the Lipschitz
+        constant is approximated by some theoretical result.
+
+    :arg f: the right-hand side function used in the differential equation.
+    :arg alpha: fractional order of the derivative.
+    :arg yspan: an expected domain of the state variable. The Lipschitz constant
+        is estimated by sampling in this domain.
+    :arg theta: a coefficient used to control the time step.
+
+    :returns: a callable returning an estimate of the time step based on the
+        Lipschitz constant.
+    """
+    from functools import partial
+    from math import gamma
+
+    if not 0.0 < alpha < 1.0:
+        raise ValueError(f"Derivative order must be in (0, 1): {alpha}")
+
+    if theta <= 0.0:
+        raise ValueError(f"Theta constant cannot be negative: {theta}")
+
+    if yspan[0] > yspan[1]:
+        yspan = (yspan[1], yspan[0])
+
+    from pycaputo.lipschitz import estimate_lipschitz_constant
+
+    def predict_time_step(t: float, y: Array) -> float:
+        if y.shape != (1,):
+            raise ValueError(f"Only scalar functions are supported: {y.shape}")
+
+        L = estimate_lipschitz_constant(partial(f, t), yspan[0], yspan[1])
+        return float((gamma(2 - alpha) * L) ** (-1.0 / alpha))
+
+    return predict_time_step
+
+
+# }}}
+
+
 # {{{ events
 
 

pycaputo/lipschitz.py

@@ -0,0 +1,160 @@
+# SPDX-FileCopyrightText: 2023 Alexandru Fikl <alexfikl@gmail.com>
+# SPDX-License-Identifier: MIT
+
+from __future__ import annotations
+
+import numpy as np
+import scipy.stats
+
+from pycaputo.logging import get_logger
+from pycaputo.utils import Array, ScalarFunction
+
+logger = get_logger(__name__)
+
+
+# {{{ estimate_lipschitz_constant
+
+
+def uniform_diagonal_sample(
+    a: float,
+    b: float,
+    n: int,
+    *,
+    delta: float = 0.05,
+    rng: np.random.Generator | None = None,
+) -> tuple[Array, Array]:
+    r"""Sample points uniformly around the diagonal of :math:`[a, b] \times [a, b]`.
+
+    This function samples point in
+
+    .. math::
+
+        \{(x, y) \in [a, b] \times [a, b] \mid |x - y| < \delta\}.
+
+    :arg n: number of points to sample.
+    :arg rng: random number generator to use, which defaults to
+        :func:`numpy.random.default_rng`.
+    :returns: an array of shape ``(2, n)``.
+    """
+    if rng is None:
+        rng = np.random.default_rng()
+
+    # TODO: is this actually uniformly distributed on the domain?
+    # FIXME: any nice vectorized way to do this?
+    x, y = np.empty((2, n))
+    for i in range(n):
+        x[i] = rng.uniform(a, b)
+        y[i] = rng.uniform(max(a, x[i] - delta), min(x[i] + delta, b))
+
+    return x, y
+
+
+def uniform_sample_maximum_slopes(
+    f: ScalarFunction,
+    a: float,
+    b: float,
+    *,
+    delta: float | None = None,
+    nslopes: int = 8,
+    nbatches: int = 64,
+    rng: np.random.Generator | None = None,
+) -> Array:
+    r"""Uniformly sample slopes of *f* in a neighborhood the diagonal of
+    :math:`[a, b] \times [a, b]`.
+
+    This function uses :func:`uniform_diagonal_sample` and takes the same
+    arguments as :func:`estimate_lipschitz_constant`.
+
+    :returns: an array of shape ``(nbatches,)`` of maximum slopes for the function
+        *f* on the interval :math:`[a, b]`.
+    """
+    if a > b:
+        raise ValueError(f"Invalid interval: [{a}, {b}]")
+
+    d = np.sqrt(2) * (b - a)
+    if delta is None:
+        # NOTE: in the numerical tests from [Wood1996] they take `delta = 0.05`
+        # regardless of the interval, but here we take a scaled version by default
+        # that's about `0.044` on `[-1, 1]` and grows on bigger intervals
+        delta = 1 / 64 * d
+
+    if delta <= 0:
+        raise ValueError(f"delta should be positive: {delta}")
+
+    if rng is None:
+        rng = np.random.default_rng()
+
+    smax = np.empty(nbatches)
+    for m in range(nbatches):
+        x, y = uniform_diagonal_sample(a, b, nslopes, delta=delta, rng=rng)
+        s = np.abs(f(x) - f(y)) / (np.abs(x - y) + 1.0e-15)
+
+        smax[m] = np.max(s)
+
+    return smax
+
+
+def fit_reverse_weibull(x: Array) -> scipy.stats.rv_continuous:
+    """Fits a reverse Weibull distribution to the CDF of *x*.
+
+    See :data:`scipy.stats.weibull_max` for details on the distribution.
+
+    :returns: the distribution with the optimal parameters.
+    """
+    # NOTE: MLE seems to work better than MM
+    c, loc, scale = scipy.stats.weibull_max.fit(x, method="MLE")
+    return scipy.stats.weibull_max(c=c, loc=loc, scale=scale)
+
+
+def estimate_lipschitz_constant(
+    f: ScalarFunction,
+    a: float,
+    b: float,
+    *,
+    delta: float | None = None,
+    nslopes: int = 8,
+    nbatches: int = 64,
+    rng: np.random.Generator | None = None,
+) -> float:
+    r"""Estimate the Lipschitz constant of *f* based on [Wood1996]_.
+
+    The Lipschitz constant is defined, for a function
+    :math:`f: [a, b] \to \mathbb{R}`, by
+
+    .. math::
+
+        |f(x) - f(y)| \le L |x - y|, \qquad \forall x, y \in [a, b]
+
+    .. warning::
+
+        The current implementation of this function seems to underpredict the
+        Lipschitz constant, unlike the results from [Wood1996]_. This is likely
+        a consequence of the method used to fit the data to the Weibull
+        distribution.
+
+    :arg a: left-hand side of the domain.
+    :arg b: right-hand side of the domain.
+    :arg delta: a width of a strip around the diagonal of the
+        :math:`[a, b] \times [a, b]` domain for :math:`(x, y)` sampling.
+    :arg nslopes: number of slopes to sample.
+    :arg nbatches: number of batches of slopes to sample.
+
+    :returns: an estimate of the Lipschitz constant.
+    """
+
+    smax = uniform_sample_maximum_slopes(
+        f,
+        a,
+        b,
+        delta=delta,
+        nslopes=nslopes,
+        nbatches=nbatches,
+        rng=rng,
+    )
+    cdf = fit_reverse_weibull(smax)
+
+    # NOTE: the estimate for the Lipschitz constant is the location parameter
+    return float(cdf.kwds["loc"])
+
+
+# }}}

pycaputo/logging.py

@@ -6,6 +6,21 @@
 import logging
 import os
 
+import rich
+
+
+def stringify_table(table: rich.table.Table) -> str:
+    """Stringify a rich table."""
+    import io
+
+    from rich.console import Console
+
+    file = io.StringIO()
+    console = Console(file=file)
+    console.print(table)
+
+    return str(file.getvalue())
+
 
 def get_logger(
     module: str,

pyproject.toml

@@ -98,6 +98,7 @@ ignore = [
     "E402",     # module-import-not-at-top-of-file
     "I001",     # unsorted-imports
     "N806",     # non-lowercase-variable-in-function
+    "N803",     # non-lowercase-argument
     "PLR0911",  # too-many-return-statements
     "PLR0912",  # too-many-branches
     "PLR0913",  # too-many-arguments

tests/test_fode.py

@@ -15,7 +15,7 @@
 from pycaputo.logging import get_logger
 from pycaputo.utils import Array, set_recommended_matplotlib
 
-logger = get_logger("pycaputo.test_caputo")
+logger = get_logger("pycaputo.test_fode")
 set_recommended_matplotlib()