Wynn epsilon table algorithm

epsilon.py

# -*- coding:utf-8 -*-
from __future__ import division, absolute_import, print_function

import sys
import math
import collections

try:
    from math import inf, nan
except ImportError:
    inf = float('inf')
    nan = float('nan')

eps = sys.float_info.epsilon


class WynnEpsilon(object):
    r"""Wynn's epsilon algorithm.

    Extrapolates the limit of the sequence::

       {s0, s1, s2, ...}

    Parameters
    ----------
    sequence : iterable of array_like, optional
        The sequence or its initial part.
    max_rows : int, optional
        Maximum number of table rows to retain.
        Must be >= 2. Default: all
    max_cols : int, optional
        Maximum number of table columns to compute.
        Default: all
    irregular_tol : float, optional
        Tolerance threshold for detecting irregular behavior
        in the table.

    Attributes
    ----------
    max_rows
    max_cols
    irregular_tol
    limit_estimate
    error_estimate
    epsilon
    
    Methods
    -------
    append
    extend

    Notes
    -----
    The algorithm computes the even-epsilon table

    .. math::

       \begin{array}{c|ccccc}
            p   &  k=-1  &  k=0  &       k=1        &       k=2        & \cdots \\
         \hline%& -------&-------&------------------&------------------&----------
            0   & \infty &  s_0  &                  &                  &        \\
            1   & \infty &  s_1  &                  &                  &        \\
            2   & \infty &  s_2  & \epsilon_2^{(0)} &                  &        \\
            3   & \infty &  s_3  & \epsilon_2^{(1)} & \epsilon_4^{(0)} &        \\
         \vdots
       \end{array}

    Note that usually in literature this is expressed in a shifted
    form, where the rows above are diagonals.

    The values :math:`\epsilon` are associated with Padé approximants
    related to the sequence. [1]_

    Only last two rows are required for the recurrence relation::

         e3   e0   ##
         ##   e1   ##
         ##   e2   new

         d1 = e1 - e3
         d2 = e2 - e1
         d3 = e1 - e0

         new = e1 + 1/(1/d1 + 1/d2 - 1/d3)

    When divisors become small, this can signal convergence or
    irregular behavior of the extrapolant. [2]_ Handling of these
    situations in this implementation follows closely that in
    QUADPACK/QELG. [3]_

    References
    ----------
    [1] P. Wynn, "Upon Systems of Recursions which Obtain among
        the Quotients of the Pade Table,"
        Numerische Mathematik 8 (3), 264 (1966).
    [2] R. P. Eddy, "The Even-Rho and Even-Epsilon Algorithms
        for Accelerating Convergence of a Numerical Sequence,"
        DTNSRDC-81/083, David W Taylor Naval Ship Research and
        development center (Bethesda MD, 1981).
    [3] R. Piessens, E. de Doncker, QUADPACK (1983).

    """
    def __init__(self, sequence=None, max_rows=None, max_cols=None, irregular_tol=1e-4):
        self.limit_estimate = nan
        self.error_estimate = nan

        self._prev_results = collections.deque([], 3)

        self.irregular_tol = irregular_tol

        if max_rows is not None and max_rows != inf:
            if not (max_rows >= 2):
                raise ValueError("max_rows must be >= 2")
            self.epsilon = collections.deque([], max_rows)
        else:
            self.epsilon = collections.deque([])

        self.max_cols = max_cols

        if sequence is not None:
            self.extend(sequence)

    def append(self, s):
        """
        Add the next known element in the sequence.
        """
        # Compute new elements
        if len(self.epsilon) < 2:
            n = 0
        else:
            n = min(len(self.epsilon[-1]), len(self.epsilon[-2]))

        if self.max_cols is not None:
            n = min(n, self.max_cols)

        best_error = 0*s + inf
        best_result = s

        new_row = [s]

        for k in range(n):
            if k == 0:
                e3 = inf
            else:
                e3 = self.epsilon[-2][k-1]
            e0 = self.epsilon[-2][k]
            e1 = self.epsilon[-1][k]
            e2 = new_row[k]

            d1 = e1 - e3   # = 1 / H_k^{p-1}
            d2 = e2 - e1   # = 1 / V_k^p
            d3 = e1 - e0   # = 1 / V_k^{p-1}

            err1 = abs(d1)
            tol1 = max(abs(e1), abs(e3)) * eps

            err2 = abs(d2)
            tol2 = max(abs(e2), abs(e1)) * eps

            err3 = abs(d3)
            tol3 = max(abs(e1), abs(e0)) * eps

            if err2 <= tol2 and err3 <= tol3:
                # Close to convergence
                best_result = e2
                best_error = err2 + err3
                break

            if (err1 <= tol1 and k > 0) or err2 <= tol2 or err3 <= tol3:
                # Two elements close to each other: drop columns
                break

            ss = 1/d1 + 1/d2 - 1/d3
            epsinf = abs(ss * e1)

            if epsinf < self.irregular_tol:
                # Irregular behavior in the table: drop columns
                break

            new = e1 + 1/ss
            error = err2 + abs(new - e2) + abs(err3)

            if error < best_error:
                best_result = new
                best_error = error

            new_row.append(new)

        self.epsilon.append(new_row)

        self.limit_estimate = best_result
        self._prev_results.append(best_result)

        if len(self._prev_results) < 3:
            self.error_estimate = inf
        else:
            self.error_estimate = max(abs(best_result - x)
                                      for x in self._prev_results)
            self.error_estimate = max(self.error_estimate,
                                      eps*abs(self.limit_estimate))

    def extend(self, terms):
        for term in terms:
            self.append(term)

    def __str__(self):
        ncols = max(len(row) for row in self.epsilon)
        header = "  ".join(["="*15]*ncols)
        text = header + "\n"
        text += "\n".join("  ".join("{:^15g}".format(x) for x in row)
                          for row in self.epsilon)
        text += "\n" + header
        text += "\n-> {} +/- {}".format(self.limit_estimate, self.error_estimate)
        return text


def test_pi_limit():
    import pytest

    def s(n):
        # s(n) -> 4 atan(1)
        return sum((-1)**j * 4 / (2 * j + 1) for j in range(n + 1))

    etab = WynnEpsilon()
    etab.extend([s(j) for j in range(5)])

    # Compare to [Graves-Morris et al., J. Comput. Appl. Math. 122, 51 (2000)]
    # for diagonals.
    expected = [[4.0],
                [2.667],
                [3.467, 3.167],
                [2.895, 3.133],
                [3.340, 3.145, 3.142]]
    for j, row in enumerate(expected):
        assert etab.epsilon[j] == pytest.approx(row, abs=1e-3)

    assert etab.limit_estimate == pytest.approx(math.pi, abs=etab.error_estimate)

    print(etab)