Sympyはdtypeオブジェクト形式でnumpyを出力します

Question

次のコードを使用して、spherical harmonics functions Y_l^m（全球で4π正規化された）とその派生派生関数のシンボリックシンピー表現を作成し、シータにおけるいくつかの等間隔のグリッド上でそれらを評価し、phiが座標：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オブジェクトを持つ配列を生成する場合がありますが、これは特に高倍音高調波に影響します。配列の項目は、複素数の代わりに複雑な1タプルとして表示されます。しかし、問題は、実際のdtypeを持つ配列にブロードキャストされるため、その配列の実数部を取ることは効果がなく、エラーが発生するということです。これには特別な理由はありますか？出力が不均質ではないので、私はすぐには見ません。 np.asarrayを使用してdtype複合体に追加キャストすることなくこれを変更する方法はありますか？追加の計算時間がかかるだけで、プログラムは少し複雑になりますが、もっと重要なのは混乱します。
2. ラムダ関数を作成する前にsympy.Nを使って式を評価していることに気づかれたかもしれません。その理由は、球面調和関数の前のプリファクターは、その理由のどれがその数のsqrtを計算できないかを知っている人にとって、長い形式とnumpyのいくつかの場合にあります。これは一般的には真実ではありませんが（np.sqrt(9L) = 3.0）、この場合、オブジェクトにはsqrtという属性がないことを示すエラーメッセージが表示されます。これはラムダ関数の生成にも関係していると思います。 sympyに、毎回float形式の記号表現を与えるように指示する方法はありますか？あるいは、何とかlambdifyコールを変更するのがよいでしょうか？

これらの問題を確認する場合は、コードブロックをスタンドアロンでテスト可能にする必要があります。 sympy.Nとnp.asarray式を削除するだけです。最初の質問は、以前に出現したエラーに関連しています。ここで35となるlmaxまでのY_l^m世代はおよそ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

なぜオブジェクトコードがオブジェクト配列を生成するのかという疑問は、私たちはMCVeなしでは簡単には答えられません。ちょうど起こることはありません。再現可能でなければなりません。配列はオブジェクトである場合

はしかし、それは簡単にあなたは多くの計算コストなしですべての結果に適用することができ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)はあなたが望むものを生み出していないことに注意してください。

私はあなたが使用している参照画面の外にスクロールするものを含め、あなたのコードでは、より探し

---

：

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

私astype(complex, copy=False)と同じです。

既に複雑な配列の場合、計算コストは​​最小です。オブジェクト配列の場合は、新しい配列を作成しなければならず、コストもより大きくなります。しかし、sympyでオブジェクト配列を作成しているとわからない場合は、コストをかけて生きなければなりません。