Derivative with numpy 1.25.2 errors with complex valued function

[rparini]

This happens on the latest version of numpy (1.25.2), I guess something changed about the percentile function that numdifftools is using. It works fine on numpy 1.24.4.

>>> import numpy as np
>>> np.__version__
'1.25.2'
>>> import numdifftools
>>> numdifftools.__version__
'0.9.41'
>>> f = lambda z: 1j*z**2
>>> df = numdifftools.Derivative(f)
>>> df(2)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/rparini/.pyenv/versions/3.9.13/lib/python3.9/site-packages/numdifftools/core.py", line 289, in __call__
    derivative, info = self._extrapolate(*results)
  File "/Users/rparini/.pyenv/versions/3.9.13/lib/python3.9/site-packages/numdifftools/limits.py", line 206, in _extrapolate
    der, info = self._get_best_estimate(der1, errors1, steps, shape)
  File "/Users/rparini/.pyenv/versions/3.9.13/lib/python3.9/site-packages/numdifftools/limits.py", line 186, in _get_best_estimate
    errors += _Limit._add_error_to_outliers(der)
  File "/Users/rparini/.pyenv/versions/3.9.13/lib/python3.9/site-packages/numdifftools/limits.py", line 170, in _add_error_to_outliers
    p25, median, p75 = np.percentile(der, [25,50, 75], axis=0)
  File "/Users/rparini/.pyenv/versions/3.9.13/lib/python3.9/site-packages/numpy/lib/function_base.py", line 4277, in percentile
    raise TypeError("a must be an array of real numbers")
TypeError: a must be an array of real numbers

[PaulKGrimes]

This is still an issue with numpy 1.26.1

[MaxMelching]

In case people still come back to this issue, here is my quick fix for this. Though I have no idea how rigorous this is from an error propagation point of view, this makes the code run again for complex input:

def _add_error_to_outliers_fixed(der, trim_fact=10):
    """
    discard any estimate that differs wildly from the
    median of all estimates. A factor of 10 to 1 in either
    direction is probably wild enough here. The actual
    trimming factor is defined as a parameter.
    """
    if np.iscomplexobj(der):    
        return np.sqrt(
            _add_error_to_outliers_fixed(np.real(der), trim_fact)**2
            + _add_error_to_outliers_fixed(np.imag(der), trim_fact)**2
        )
    
    try:
        if np.any(np.isnan(der)):
            p25, median, p75 = np.nanpercentile(der, [25,50, 75], axis=0) 
        else:
            p25, median, p75 = np.percentile(der, [25,50, 75], axis=0)

        iqr = np.abs(p75 - p25)
    except ValueError as msg:
        warnings.warn(str(msg))
        return 0 * der

    a_median = np.abs(median)
    outliers = (((abs(der) < (a_median / trim_fact)) +
                (abs(der) > (a_median * trim_fact))) * (a_median > 1e-8) +
                ((der < p25 - 1.5 * iqr) + (p75 + 1.5 * iqr < der)))
    errors = outliers * np.abs(der - median)
    return errors

_Limit._add_error_to_outliers = staticmethod(_add_error_to_outliers_fixed)

The results seem to be sensible.

Also, something like this is required because looking through the corresponding numpy changes, it is not a bug that np.percentile does not accept complex numbers anymore. They argue that the percentile of a complex number is not even something that can be well-defined, so I do not think support for this will be brought back in the future.