Sympy 以 dtype 对象格式提供数字 numpy 输出

Question

我正在使用以下代码为spherical harmonics functions Y_l^m（在整个球面上标准化 4-pi）及其 theta 导数创建符号 Sympy 表达式，然后想在一些均匀间隔的网格上评估它们在 theta 和 phi 坐标中：

import numpy as np
from math import pi, cos, sin
import sympy
from sympy import Ynm, simplify, diff, lambdify
from sympy.abc import n,m,theta,phi

resol = 2.5
dtheta_rad_ylm = -resol * pi/180.0
dphi_rad_ylm = resol * pi/180.0

thetaarr_rad_ylm_symm = np.arange(pi+dtheta_rad_ylm/2.0,dtheta_rad_ylm/2.0,dtheta_rad_ylm)
phiarr_rad_ylm = np.arange(0.0,2*pi,dphi_rad_ylm)
phi_grid_rad_ylm, theta_grid_rad_ylm_symm = np.meshgrid(phiarr_rad_ylm, thetaarr_rad_ylm_symm)

lmax = len(thetaarr_rad_ylm_symm)/2 - 1
nmax = (lmax+1)*(lmax+2)/2

ylms_symm_full = np.zeros((lmax+1, lmax+1, len(thetaarr_rad_ylm_symm), len(phiarr_rad_ylm)))
dylms_symm_full = np.zeros((lmax+1, lmax+1, len(thetaarr_rad_ylm_symm), len(phiarr_rad_ylm)))

for n in np.arange(0,lmax+1):
  for m in np.arange(0,n+1):
    print "generating resol %s, y_%d_%d" % (resol,n,m)

    ylm_symbolic = simplify(2 * sympy.sqrt(sympy.pi) * Ynm(n,m,theta,phi).expand(func=True))
    dylm_symbolic = simplify(diff(ylm_symbolic, theta))

    # activate and deactivate comments for second-question-related error
    # error appears later than the first-question-related error!
    ylm_lambda = lambdify((theta,phi), sympy.N(ylm_symbolic), "numpy")
    dylm_lambda = lambdify((theta,phi), sympy.N(dylm_symbolic), "numpy")
#    ylm_lambda = lambdify((theta,phi), ylm_symbolic, "numpy")
#    dylm_lambda = lambdify((theta,phi), dylm_symbolic, "numpy")

    # activate and deactivate comments for first-question-related error
    ylm_symm_full = np.asarray(ylm_lambda(theta_grid_rad_ylm_symm, phi_grid_rad_ylm), dtype=complex)
    dylm_symm_full = np.asarray(dylm_lambda(theta_grid_rad_ylm_symm, phi_grid_rad_ylm), dtype=complex)
#    ylm_symm_full = ylm_lambda(theta_grid_rad_ylm_symm, phi_grid_rad_ylm)
#    dylm_symm_full = dylm_lambda(theta_grid_rad_ylm_symm, phi_grid_rad_ylm)

    if n == 0 and m == 0:
      ylm_symm_full = np.tile(ylm_symm_full, (len(thetaarr_rad_ylm_symm), len(phiarr_rad_ylm)))
      dylm_symm_full = np.tile(dylm_symm_full, (len(thetaarr_rad_ylm_symm), len(phiarr_rad_ylm)))

    ylms_symm_full[n,m,:,:] = np.real(ylm_symm_full)
    dylms_symm_full[n,m,:,:] = np.real(dylm_symm_full)

还有其他几个包提供了在没有符号表达式的情况下生成数字 Y_l^m 的功能，例如 scipy.special.sph_harm。但是，获得“精确”导数对我来说至关重要，即不使用任何数值微分方法，例如有限差分 (np.gradient)。因此，在获得 Y_l^m 的符号公式并“尽可能”简化这些公式后，使用 numpy 后端创建 lambda 函数（以便能够进行矢量化计算），然后在网格上评估这些函数。最后我只需要球谐函数的实部（我知道我也可以用 Znm 而不是 Ynm 创建真正的球谐函数，但是......）。

两个问题：

1. 大多数情况下，数值输出随后作为 dtype complex 或 np.complex128 的通常 2d-numpy 数组给出。然而，在某些情况下，Sympy 会生成具有 dtype 对象的数组，这尤其会影响高 l 球谐函数。数组条目显示为复数 1 元组，而不仅仅是复数。然而，问题是在该数组上取实部没有效果，从而导致错误，因为它被广播到具有真实 dtype 的数组中。这有什么特别的原因吗？我没有看到任何直接的，因为输出不是不均匀的。有什么方法可以改变它，而不必使用np.asarray 将其额外转换为 dtype complex？这只需要额外的计算时间，使程序稍微复杂一些，但更重要的是令人困惑。
2. 您可能还注意到，在创建 lambda 函数之前，我已经使用 sympy.N 来计算表达式。原因是球谐函数前面的前置因子在某些情况下是长格式和 numpy 的，因为无论谁知道是什么原因，都无法计算该数字的 sqrt。请注意，这通常不是真的 (np.sqrt(9L) = 3.0)，但在这种情况下，会出现一条错误消息，指出 long 对象没有属性 sqrt。我想这也与 lambda 函数的生成有关。有什么方法可以告诉 Sympy 每次都以浮点格式给出符号表达式吗？或者，更好的是，以某种方式修改lambdify 调用？

如果您想检查这些问题，代码块应该是独立且可测试的。只需删除 sympy.N 和 np.asarray 表达式。第一个问题与之前出现的错误有关。 Y_l^m 生成到 lmax 这里是 35 大约需要 10-15 分钟。

提前感谢您的帮助！

---

更新：以下是一些最小、完整且可验证的示例。对于两者，请导入所需的包：

import numpy as np
from math import pi, cos, sin
import sympy
from sympy import Ynm, simplify, diff, lambdify
from sympy.abc import n,m,theta,phi

错误 #1： an = 31, m = 1 处的对象 dtype 问题：

# minimal, complete and verifiable example (MCVe) #1
# error message:

#---> 43     dylms_symm_full[n,m,:,:] = np.real(dylm_symm_full)
#TypeError: can't convert complex to float

ylm_symbolic = simplify(2 * sympy.sqrt(sympy.pi) * Ynm(31,1,theta,phi).expand(func=True))
dylm_symbolic = simplify(diff(ylm_symbolic, theta))

ylm_lambda = lambdify((theta,phi), ylm_symbolic, "numpy")
dylm_lambda = lambdify((theta,phi), dylm_symbolic, "numpy")

ylm_symm_full = ylm_lambda(theta_grid_rad_ylm_symm, phi_grid_rad_ylm)
dylm_symm_full = dylm_lambda(theta_grid_rad_ylm_symm, phi_grid_rad_ylm)

ylms_symm_full = np.zeros((len(thetaarr_rad_ylm_symm), len(phiarr_rad_ylm)))
dylms_symm_full = np.zeros((len(thetaarr_rad_ylm_symm), len(phiarr_rad_ylm)))

ylms_symm_full[:,:] = np.real(ylm_symm_full)
dylms_symm_full[:,:] = np.real(dylm_symm_full)

print ylm_symm_full
print dylm_symm_full

错误 #2： 在 n = 32，m = 29 时出现长 sqrt 属性问题：

# minimal, complete and verifiable example (MCVe) #2
# error message:

#---> 33     ylm_symm_full = np.asarray(ylm_lambda(theta_grid_rad_ylm_symm, phi_grid_rad_ylm), dtype=complex)
#/opt/local/anaconda/anaconda-2.2.0/lib/python2.7/site-packages/numpy/__init__.pyc in <lambda>(_Dummy_4374, _Dummy_4375)
#AttributeError: 'long' object has no attribute 'sqrt'

ylm_symbolic = simplify(2 * sympy.sqrt(sympy.pi) * Ynm(32,29,theta,phi).expand(func=True))
dylm_symbolic = simplify(diff(ylm_symbolic, theta))

ylm_lambda = lambdify((theta,phi), ylm_symbolic, "numpy")
dylm_lambda = lambdify((theta,phi), dylm_symbolic, "numpy")

ylm_symm_full = np.asarray(ylm_lambda(theta_grid_rad_ylm_symm, phi_grid_rad_ylm), dtype=complex)
dylm_symm_full = np.asarray(dylm_lambda(theta_grid_rad_ylm_symm, phi_grid_rad_ylm), dtype=complex)

ylms_symm_full = np.zeros((len(thetaarr_rad_ylm_symm), len(phiarr_rad_ylm)))
dylms_symm_full = np.zeros((len(thetaarr_rad_ylm_symm), len(phiarr_rad_ylm)))

ylms_symm_full[:,:] = np.real(ylm_symm_full)
dylms_symm_full[:,:] = np.real(dylm_symm_full)

print ylm_symbolic                # the symbolic Y_32^29 expression
print type(175844649714253329810) # the number that causes the problem

Answer 1

为什么您的代码有时会生成一个对象数组，这不是我们在没有 MCVe 的情况下无法轻易回答的问题 - 它不能只是偶尔发生，它必须是可重现的。

但是如果数组是对象，它可能很容易转换为复杂的

arr.astype(np.complex)

使用copy=False 参数，您可以将其应用于所有结果，而无需太多计算成本。

arr.astype(np.complex, copy=False).real

对象版本的元素不是元组；它们是标量复数值，就这样打印。

In [187]: arr=np.random.rand(3)+np.random.rand(3)*1j
In [188]: arrO=arr.astype(object)
In [189]: arrO
Out[189]: 
array([(0.6129476673822528+0.09323924558124808j),
       (0.9540542895309456+0.81929476753951j),
       (0.8068967867200485+0.9494305517611881j)], dtype=object)
In [190]: type(arrO[0])
Out[190]: complex
In [191]: arr.real
Out[191]: array([ 0.61294767,  0.95405429,  0.80689679])
In [193]: arrO[0]
Out[193]: (0.6129476673822528+0.09323924558124808j)
In [194]: arrO.astype(np.complex).real
Out[194]: array([ 0.61294767,  0.95405429,  0.80689679])

一些数学运算会对对象数组的元素进行“渗透”，但real 不是其中之一。所以你注意到np.real(arrO) 不会产生你想要的东西。

---

更多地查看你的代码，包括我看到你正在使用的从屏幕上滚动出来的东西：

np.asarray(dylm_lambda(...), dtype=complex)

这和我的astype(complex, copy=False)一样。

对于已经很复杂的数组，计算成本是最小的。对于对象数组，它必须创建一个新数组，并且成本更高。但是，如果您无法通过sympy 确定对象数组是在创建对象数组，那么您就必须忍受成本。