Sympy以dtype对象形式给出数值numpy输出

Question

我使用下面的代码为spherical harmonics functionsY_l^m（在整个球体上规范化4-pi）及其θ导数创建符号Sympy表达式，然后在θ和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μm的功能，比如scipy.special.sph_harm。然而，对于我来说，得到一个“精确”的导数是至关重要的，即不使用任何数值微分方法，例如有限差分（np.gradient）。因此，在得到Yμl^m的符号公式并“尽可能地”简化之后，使用numpy后端创建lambda函数（以便能够进行矢量化计算），然后在网格上计算这些函数。最后我只需要球谐函数的实部（我知道我也可以用Znm而不是Ynm来创建真正的球谐函数，但是…）。在

两个问题：

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

如果您想检查这些问题，代码块应该是独立的和可测试的。只要去掉sympy.N和np.A阵列表达。第一个问题与前面出现的错误有关。你的最大发电量是35，大约需要10-15分钟。在

提前感谢您的帮助！在

---

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

^{pr2}$

错误1:an=31，m=1处的对象数据类型问题：

# 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就很难回答的问题——它不能只是偶尔发生，它必须是可重复的。在

但是，如果数组是object，则很容易将其转换为complex with

^{1}$

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

^{pr2}$

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

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来确定是在创建对象数组，那么您必须承受这样的代价。在