Sympy는 dtype 객체 형식으로 숫자로 numpy 출력을 제공합니다.

Question

다음 코드를 사용하여 spherical harmonics functions Y_l^m (전체 구에 대해 4-pi- 정규화 됨) 및 해당 theta 파생물에 대한 상징적 인 Sympy 표현식을 만든 다음 세타 일부 균등 격자들을 평가하고, 피 좌표 : scipy.special.sph_harm 같은 상징적 표현없이 숫자 Y_l^m을 생성하는 기능을 제공하는 여러 다른 패키지있다

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) 

. 그러나 "정확한"파생어를 얻는 것이 중요합니다. 즉, 예를 들어 숫자 미분법을 사용하지 않아야합니다. 유한 차이 (np.gradient). 따라서 Y_l^m에 대한 기호식을 가져 와서 "가능한 한 많이"단순화하면 람다 함수가 numpy 백엔드 (벡터화 된 계산을 수행 할 수 있음)를 사용하여 만들어지고 그리드에서 평가됩니다. 마지막으로 구형 고조파의 실제 부분 만 필요합니다 (Ynm 대신 Znm을 사용하여 실제 구형 고조파를 생성 할 수도 있음을 알고 있지만 ...).

두 질문 :

1. 대부분, 수치 출력은 DTYPE 복잡하거나 np.complex128의 일반적인 2D-NumPy와 배열로 제공됩니다. 그러나 어떤 경우에는 Sympy가 dtype 객체를 갖는 배열을 생성합니다. 이것은 특히 높은 l 구면 고조파에 영향을 미칩니다. 배열 항목은 복소수가 아닌 복잡한 1-tuples로 표시됩니다. 그러나 문제는 실제 dtype이있는 배열로 브로드 캐스팅되기 때문에 해당 배열의 실수 부분을 가져 오는 것이 아무런 영향을 미치지 않으므로 오류가 발생한다는 것입니다. 특별한 이유가 있습니까? 출력이 불균일하지 않기 때문에 어떤 즉각적인 것을 보지 못합니다. np.asarray을 사용하여 dtype 컴플렉스에 추가로 캐스팅하지 않고이를 변경할 수있는 방법은 없습니까? 단지 계산 시간이 추가로 걸리므로 프로그램이 약간 복잡해 지지만 더 중요한 것은 혼란 스럽습니다.
2. 람다 함수를 생성하기 전에 이미 sympy.N을 사용하여 표현식을 평가했다는 것을 알았을 수도 있습니다. 그 이유는 구형 고조파 앞에있는 프리 팩터가 어떤 경우에는 긴 형식과 numpy인데, 이유를 아는 사람이라면 그 수의 sqrt를 계산할 수 없기 때문입니다. 이것은 일반적으로 사실이 아니며 (np.sqrt(9L) = 3.0),이 경우에는 긴 객체에 sqrt 속성이 없음이라는 오류 메시지가 표시됩니다. 나는 이것이 람다 함수 생성과 관련이 있다고 가정한다. Sympy가 이미 매번 float 형식으로 심볼릭 표현을하도록 지시하는 방법이 있습니까? 아니면 lambdify 호출을 수정하는 것이 더 나은가?

이러한 문제를 확인하려면 코드 블록을 독립 실행 형으로 테스트 할 수 있어야합니다. sympy.N 및 np.asarray 표현식을 제거하기 만하면됩니다. 첫 번째 질문은 이전에 나타난 오류와 관련이 있습니다. Y_l^m 세대에서 lmax까지 35 분이 소요됩니다. 대략 10-15 분이 걸립니다.

미리 도움 주셔서 감사합니다.

---

UPDATE

: 여기 몇 가지 최소한의 완전하고 검증 가능한 예입니다.

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 : = 31 일 객체 DTYPE 문제, 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 : 긴 SQRT 속성 문제가 모두 필요한 패키지를 가져하시기 바랍니다 n = 32에서 m = 29 :

# 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

왜 코드가 개체 배열을 생성하는지는 somethi가 아닙니다. 우리는 MCV 없이는 쉽게 대답 할 수 없습니다. 때때로 일어날 수있는 일은 아니며 재현 가능해야합니다. 배열 객체 인 경우

그러나, 그것은 쉽게 많은 계산 비용없이 모든 결과에 적용 할 수있는 copy=False 매개 변수로

arr.astype(np.complex) 

복잡한로 변환 할 수 있습니다.

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)은 원하는 것을 만들어 내지 못합니다.

---

난 당신이 사용하는 참조 화면 오프 스크롤 물건을 포함하여 코드를 더 상대 : 내 astype(complex, copy=False)와 동일합니다

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

합니다.

이미 복잡한 배열의 경우 계산 비용이 최소화됩니다. 객체 배열의 경우 새로운 배열을 만들어야하고 비용은 더 큽니다. 그러나 sympy으로 객체 배열을 만드는 것으로 알 수 없다면 비용을 지불해야합니다.