regression in lambdify: module=['scipy'] ignored for sqrt

[dblaschke-LANL]

sympy 1.11 has a regression where despite calling lambdify with the module=['scipy'] option, the sqrt function is nonetheless taken from numpy, thus breaking code for complex numbers (sympy <=1.10.x worked just fine in this regard).

Here is a mwe:

---

import sympy as sp
x = sp.symbols('x')
f = sp.lambdify(x,sp.sqrt(x),modules=["scipy"])
f(-2)
<lambdifygenerated-3>:2: RuntimeWarning: invalid value encountered in sqrt
return sqrt(x)

---

The problem is seen here:

---

?f
Signature: f(x)
Docstring:
Created with lambdify. Signature:

func(x)

Expression:

sqrt(x)

Source code:

def _lambdifygenerated(x):
return sqrt(x)

Imported modules:

from numpy import sqrt

---

Note how sqrt is imported from numpy under "imported modules" (this was not the case in sympy 1.10).
This should not happen with the modules=["scipy"] option.

Current workaround:

---

import scipy
f = sp.lambdify(x,sp.sqrt(x),modules=["scipy",{'sqrt':scipy.sqrt}])
f(-2)
Out[7]: 1.4142135623730951j

---

[oscarbenjamin]

this was not the case in sympy 1.10

Are you sure about this? I just tried with SymPy 1.10 and I see the same output as with SymPy 1.11 as well as current master.

In [3]: import sympy as sp
   ...: x = sp.symbols('x')
   ...: f = sp.lambdify(x,sp.sqrt(x),modules=["scipy"])
   ...: f(-2)
<lambdifygenerated-2>:2: RuntimeWarning: invalid value encountered in sqrt
  return sqrt(x)
Out[3]: nan

In [4]: sp.__version__
Out[4]: '1.10'

Also when I use the scipy.sqrt function I get a clear deprecation warning saying to use a (different) numpy function instead:

In [1]: from scipy import sqrt

In [2]: sqrt(-2)
<ipython-input-2-84fde6a6eea1>:1: DeprecationWarning: scipy.sqrt is deprecated and will be removed in SciPy 2.0.0, use numpy.lib.scimath.sqrt instead
  sqrt(-2)
Out[2]: 1.4142135623730951j

[dblaschke-LANL]

it was in fact sympy 1.10.1 from an anaconda python 3.10 environment that was working fine.

I didn't know about the deprecation of scipy.sqrt;

Perhaps it would make sense to replace

from numpy import sqrt

with

from numpy.lib.scimath import sqrt

in the case where lambidfy is called with just the modules=['scipy'] option to emulate the 'old' behavior?

Thanks for the quick response!

[oscarbenjamin]

I see the same with SymPy 1.10.1 and also many older versions of SymPy:

In [1]: import sympy as sp
   ...: x = sp.symbols('x')
   ...: f = sp.lambdify(x,sp.sqrt(x),modules=["scipy"])
   ...: f(-2)
<lambdifygenerated-1>:2: RuntimeWarning: invalid value encountered in sqrt
  return sqrt(x)
Out[1]: nan

In [2]: sp.__version__
Out[2]: '1.10.1'

I don't think that any change you are seeing is because of any change in SymPy. Have you also updated other libraries? Are you sure that it is the code shown that was working before?

[dblaschke-LANL]

yes, libraries like numpy, scipy have also been upgraded.
in older sympy versions (<=1.10.1) I see this behavior:

import sympy as sp

In [2]: sp.version
Out[2]: '1.10.1'

In [3]: import numpy as np

In [4]: np.version
Out[4]: '1.21.5'

In [5]: import scipy

In [6]: scipy.version
Out[6]: '1.7.3'

In [7]: x = sp.symbols('x')

In [8]: f = sp.lambdify(x,sp.sqrt(x),modules=["scipy"])

In [9]: f(-2)
Out[9]: 1.4142135623730951j

In [10]: import sys

In [11]: sys.version_info
Out[11]: sys.version_info(major=3, minor=10, micro=4, releaselevel='final', serial=0)

Above was an anaconda environment on macos, but I also see the same behaviour on ubuntu 22.04 (with sympy 1.9, numpy 1.21.5 and scipy 1.8)

Nonetheless, going forward I believe it would make sense for lambdify to use a sqrt implementation that supports complex numbers if the user has already requested module=['scipy'], would you agree?

[asmeurer]

Are you sure about this? I just tried with SymPy 1.10 and I see the same output as with SymPy 1.11 as well as current master.

This seems pretty clearly to be the result of #23945. scipy.sqrt apparently does automatically promote negative real arguments to complex numbers, but it's deprecated for numpy.sqrt which doesn't. We removed the usage of scipy.sqrt in favor of numpy.sqrt because of this upstream deprecation.

Generally speaking, NumPy/SciPy functions aren't going to automatically promote real dtypes (or real Python scalars) to complex scalars. If you want a complex answer, you should use -1+0j.

This sort of thing is going to be even more streamlined than it currently is in future versions of NumPy (see https://numpy.org/neps/nep-0050-scalar-promotion.html). The main motivation is that the resulting dtype of an operation should depend only on the type of the input, not its value. For sqrt, the value is complex depending on the value, so the only way it could give complex outputs for negative real (float or int) inputs would be to always return a complex dtype result even for positive inputs, which would not be desirable.

Nonetheless, going forward I believe it would make sense for lambdify to use a sqrt implementation that supports complex numbers if the user has already requested module=['scipy'], would you agree?

Just to be clear, np.sqrt does support complex numbers:

>>> np.sqrt(-1+0j)
1j

According to the deprecation warning:

<stdin>:1: DeprecationWarning: scipy.sqrt is deprecated and will be removed in SciPy 2.0.0, use numpy.lib.scimath.sqrt instead

You can use numpy.lib.scimath.sqrt. So if you really want a function that promote to complex based on values you can use this:

>>> lambdify(x, sqrt(x), [numpy.lib.scimath])(-1)
1j

I'm rather inclined for the lambdify default to remain inline with the NumPy defaults, though.

[oscarbenjamin]

Are you sure about this? I just tried with SymPy 1.10 and I see the same output as with SymPy 1.11 as well as current master.

This seems pretty clearly to be the result of #23945.

Have you actually tested different versions of sympy? I can't reproduce the situation where the behaviour was previously different. Here's a fresh venv with sympy 1.10.1 (and numpy and scipy) all installed with pip from pypi:

>>> import sympy as sp
>>> x = sp.symbols('x')
>>> f = sp.lambdify(x,sp.sqrt(x),modules=["scipy"])
>>> f(-2)
<lambdifygenerated-1>:2: RuntimeWarning: invalid value encountered in sqrt
  return sqrt(x)
nan
>>> sp.__version__
'1.10.1'

Here's what's installed:

$ pip list
Package    Version
---------- -------
mpmath     1.2.1
numpy      1.23.3
pip        20.2.3
scipy      1.9.1
setuptools 49.2.1
sympy      1.10.1

[dblaschke-LANL]

OK, did some more testing: the code works with scipy<=1.8 and sympy <=1.10.1, but returns nan with scipy 1.9 or sympy 1.11 (or both)

@ Aaron: I'm not suggesting to change any defaults, just suggesting to fall back to numpy.lib.scimath for those functions that have been deprecated in scipy, and only when the user explicitly requests scipy functions instead of numpy funcions via
sp.lambdify(x,f(x),modules=["scipy"])

(and I agree that passing modules=[numpy.lib.scimath] is another option, but the above suggestion avoids breaking older user-code)

[asmeurer]

Oh I apparently had scipy 1.5 installed on my machine. So this behavior is related to older scipy versions.

>>> lambdify(x, sqrt(x), ['scipy'])(-1)
1j
>>> import scipy
>>> scipy.__version__
'1.5.3'
>>> sympy.__version__
'1.10.1'

I don't quite get exactly what changed in the newer scipy versions. I guess it doesn't really matter. The point is that upgrading either SymPy or SciPy from these versions changes this behavior.

@ Aaron: I'm not suggesting to change any defaults, just suggesting to fall back to numpy.lib.scimath for those functions that have been deprecated in scipy, and only when the user explicitly requests scipy functions instead of numpy functions via
sp.lambdify(x,f(x),modules=["scipy"])

So the feature request here is to automatically use numpy.lib.scimath (also called numpy.emath) in lambdify when scipy is in the modules.

I think this is a bad idea. We should definitely support explicitly adding that module, but I don't like automatically injecting different modules when certain modules are present. We do inject 'numpy' with 'scipy', which is a nice convenience (see also #20294), but in general I think it's better to be explicit about this. It would be confusing for users if adding scipy also changed which NumPy functions are used.

Also I do think that the normal NumPy namespace behavior of functions that don't have value-dependent result data types is better for most applications. It's much more consistent and easier to reason about. And it's also faster, because a value-based sqrt has to read every element of the input array to check if any elements are negative before it can allocate the result array:

In [1]: import numpy as np

In [2]: a = np.arange(10000, dtype=np.float64)

In [3]: %timeit b = np.sqrt(a)
8.38 µs ± 77.1 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

In [4]: %timeit b = np.emath.sqrt(a)
30.6 µs ± 1.32 µs per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

I get that there's a backwards compatibility break here, but to me, that break is due to changes in SciPy, not in SymPy. Namely, that functions like sqrt are no longer included with SciPy. SymPy only kept up with this change ahead of the deprecation.

[dblaschke-LANL]

OK, I see your point

[alexbatgithub]

For those who are interested in an interim workaround:

f = sp.lambdify(x,sp.sqrt(x),modules=[{'sqrt':np.emath.sqrt},"scipy"])

[oscarbenjamin]

The real problem here is that lambdify's modules argument is a bad idea.

It would be better if when lambdfying you would specify what sort of type/domain you want to operate with like lambdify(..., type='complex128') or something like that. Just pulling function names off of a module object is not a good way of controlling the desired behaviour.