doublets, which are either poles with vanishingly small residues or
pole-zero pairs that are close enough together to nearly cancel, see [2]_. The
greedy nature of the AAA algorithm means Froissart doublets are rare. However, if
`rtol` is set too tight then the approximation will stagnate and many Froissart
doublets will appear. Froissart doublets can usually be removed by removing support
points and then resolving the least squares problem. The support point :math:`z_j`,
which is the closest support point to the pole :math:`a` with residue
:math:`\alpha`, is removed if the following is satisfied

.. math::

    |\alpha| / |z_j - a| < \verb|clean_up_tol| \cdot \tilde{f},

where :math:`\tilde{f}` is the geometric mean of `support_values`.


References
----------
.. [1] Y. Nakatsukasa, O. Sete, and L. N. Trefethen, "The AAA algorithm for
        rational approximation", SIAM J. Sci. Comp. 40 (2018), A1494-A1522.
        :doi:`10.1137/16M1106122`
.. [2] J. Gilewicz and M. Pindor, Pade approximants and noise: rational functions,
       J. Comp. Appl. Math. 105 (1999), pp. 285-297.
       :doi:`10.1016/S0377-0427(02)00674-X`

Examples
--------

Here we reproduce a number of the numerical examples from [1]_ as a demonstration
of the functionality offered by this method.

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from scipy.interpolate import AAA
>>> import warnings

For the first example we approximate the gamma function on ``[-3.5, 4.5]`` by
extrapolating from 100 samples in ``[-1.5, 1.5]``.

>>> from scipy.special import gamma
>>> sample_points = np.linspace(-1.5, 1.5, num=100)
>>> r = AAA(sample_points, gamma(sample_points))
>>> z = np.linspace(-3.5, 4.5, num=1000)
>>> fig, ax = plt.subplots()
>>> ax.plot(z, gamma(z), label="Gamma")
>>> ax.plot(sample_points, gamma(sample_points), label="Sample points")
>>> ax.plot(z, r(z).real, '--', label="AAA approximation")
>>> ax.set(xlabel="z", ylabel="r(z)", ylim=[-8, 8], xlim=[-3.5, 4.5])
>>> ax.legend()
>>> plt.show()

We can also view the poles of the rational approximation and their residues:

>>> order = np.argsort(r.poles())
>>> r.poles()[order]
array([-3.81591039e+00+0.j        , -3.00269049e+00+0.j        ,
       -1.99999988e+00+0.j        , -1.00000000e+00+0.j        ,
        5.85842812e-17+0.j        ,  4.77485458e+00-3.06919376j,
        4.77485458e+00+3.06919376j,  5.29095868e+00-0.97373072j,
        5.29095868e+00+0.97373072j])
>>> r.residues()[order]
array([ 0.03658074 +0.j        , -0.16915426 -0.j        ,
        0.49999915 +0.j        , -1.         +0.j        ,
        1.         +0.j        , -0.81132013 -2.30193429j,
       -0.81132013 +2.30193429j,  0.87326839+10.70148546j,
        0.87326839-10.70148546j])

For the second example, we call `AAA` with a spiral of 1000 points that wind 7.5
times around the origin in the complex plane.

>>> z = np.exp(np.linspace(-0.5, 0.5 + 15j*np.pi, 1000))
>>> r = AAA(z, np.tan(np.pi*z/2), rtol=1e-13)

We see that AAA takes 12 steps to converge with the following errors:

>>> r.errors.size
12
>>> r.errors
array([2.49261500e+01, 4.28045609e+01, 1.71346935e+01, 8.65055336e-02,
       1.27106444e-02, 9.90889874e-04, 5.86910543e-05, 1.28735561e-06,
       3.57007424e-08, 6.37007837e-10, 1.67103357e-11, 1.17112299e-13])

We can also plot the computed poles:

>>> fig, ax = plt.subplots()
>>> ax.plot(z.real, z.imag, '.', markersize=2, label="Sample points")
>>> ax.plot(r.poles().real, r.poles().imag, '.', markersize=5,
...         label="Computed poles")
>>> ax.set(xlim=[-3.5, 3.5], ylim=[-3.5, 3.5], aspect="equal")
>>> ax.legend()
>>> plt.show()

We now demonstrate the removal of Froissart doublets using the `clean_up` method
using an example from [1]_. Here we approximate the function
:math:`f(z)=\log(2 + z^4)/(1 + 16z^4)` by sampling it at 1000 roots of unity. The
algorithm is run with ``rtol=0`` and ``clean_up=False`` to deliberately cause
Froissart doublets to appear.

>>> z = np.exp(1j*2*np.pi*np.linspace(0,1, num=1000))
>>> def f(z):
...     return np.log(2 + z**4)/(1 - 16*z**4)
>>> with warnings.catch_warnings():  # filter convergence warning due to rtol=0
...     warnings.simplefilter('ignore', RuntimeWarning)
...     r = AAA(z, f(z), rtol=0, max_terms=50, clean_up=False)
>>> mask = np.abs(r.residues()) < 1e-13
>>> fig, axs = plt.subplots(ncols=2)
>>> axs[0].plot(r.poles().real[~mask], r.poles().imag[~mask], '.')
>>> axs[0].plot(r.poles().real[mask], r.poles().imag[mask], 'r.')

Now we call the `clean_up` method to remove Froissart doublets.

>>> with warnings.catch_warnings():
...     warnings.simplefilter('ignore', RuntimeWarning)
...     r.clean_up()
4
>>> mask = np.abs(r.residues()) < 1e-13
>>> axs[1].plot(r.poles().real[~mask], r.poles().imag[~mask], '.')
>>> axs[1].plot(r.poles().real[mask], r.poles().imag[mask], 'r.')
>>> plt.show()

The left image shows the poles prior of the approximation ``clean_up=False`` with
poles with residue less than ``10^-13`` in absolute value shown in red. The right
image then shows the poles after the `clean_up` method has been called.
Néd