Q&A
質問の意図がよくわかりませんが、 prob関数の中身がごちゃってるので似た処理をまとめる関数を作って整理したいということですか？
以下の時間発展の計算をひとまとめに関数として定義して、欲しい値だけを抽出するには、どのようにすれば良いのでしょうか？
（少しdef関数に苦手意識あり）
一部のコードを掲載しますが、要望があれば全部掲載します。
python
1import numpy as np 2import matplotlib.pyplot as plt 3import math 4from mpl_toolkits.mplot3d import Axes3D 5from matplotlib import cm 6import matplotlib.colors as colors 7 8n=10 #tの範囲 9m=15 #偶数 10 11p=1/2 12q=1-p 13 14#quantum coin 15P = [[-p, q, math.sqrt(p*q), math.sqrt(p*q)],[0,0,0,0],[0,0,0,0],[0,0,0,0]] 16Q = [[0,0,0,0],[q, -p, math.sqrt(p*q), math.sqrt(p*q)],[0,0,0,0],[0,0,0,0]] 17R = [[0,0,0,0],[0,0,0,0],[math.sqrt(p*q), math.sqrt(p*q), -q, p],[0,0,0,0]] 18S = [[0,0,0,0],[0,0,0,0],[0,0,0,0],[math.sqrt(p*q), math.sqrt(p*q), p, -q]] 19 20t_list = [] 21x_list = [] 22y_list = [] 23 24phi_map = np.zeros((2*m+1, 2*m+1,4)) #縦横2m+1の座標平面上にいる波動状態をセット np.zeros((行,列,[]の中身の数)) 25phi_map[0,0]= np.array([1/2,1/2,-1/2,-1/2]) 26p_map=np.zeros([2*m+1,2*m+1]) 27 28for i in range(0,2*m+1): 29 p = np.dot(phi_map[i,i], np.conj(phi_map[i,i])) 30 p_map[i,i]=p 31 x_list.append(i) 32 y_list.append(i) 33 34def prob(t,x,y): 35 for t in range(0,n+1): 36 t_list.append(t) 37 if t == 0: 38 phi_map 39 p_map 40 else: 41 next_phi_map = np.zeros((2*m+1,2*m+1, 4)) 42 for x in range(0,2*m+1): 43 if x == 0: 44 for y in range(0,2*m+1): 45 if y == 0: 46 next_phi_map[x,y] = np.array([np.inner(P, phi_map[x+1,y]) + np.inner(R, phi_map[x,y+1])]) 47 elif y == 2*m: 48 next_phi_map[x,y] = np.array([np.inner(P, phi_map[x+1,y]) + np.inner(S, phi_map[x,y-1])]) 49 else: 50 next_phi_map[x,y] = np.array([np.inner(P, phi_map[x+1,y]) + np.inner(S, phi_map[x,y-1]) + np.inner(R, phi_map[x,y+1])]) 51 elif x == 2*m: 52 for y in range(0,2*m+1): 53 if y == 0: 54 next_phi_map[x,y] = np.array([np.inner(Q, phi_map[x-1,y]) + np.inner(R, phi_map[x,y+1])]) 55 elif y == 2*m: 56 next_phi_map[x,y] = np.array([np.inner(Q, phi_map[x-1,y]) + np.inner(S, phi_map[x,y-1])]) 57 else: 58 next_phi_map[x,y] = np.array([np.inner(Q, phi_map[x-1,y]) + np.inner(S, phi_map[x,y-1]) + np.inner(R, phi_map[x,y+1])]) 59 else: 60 for y in range(0,2*m+1): 61 if y == 0: 62 next_phi_map[x,y] = np.array([np.inner(P, phi_map[x+1,y]) + np.inner(Q, phi_map[x-1,y]) + np.inner(R, phi_map[x,y+1])]) 63 elif y == 2*m: 64 next_phi_map[x,y] = np.array([np.inner(P, phi_map[x+1,y]) + np.inner(Q, phi_map[x-1,y]) + np.inner(S, phi_map[x,y-1])]) 65 else: 66 next_phi_map[x,y] = np.array([np.inner(P, phi_map[x+1,y]) + np.inner(Q, phi_map[x-1,y]) + np.inner(R, phi_map[x,y+1]) + np.inner(S, phi_map[x,y-1])]) 67 p_map[x,y] = np.dot(next_phi_map[x,y], np.conj(next_phi_map[x,y])) 68 phi_map = next_phi_map 69 70 return (p_map[x,y])#←ここは違うので無視してください。 71 72 73fig = plt.figure() 74ax = Axes3D(fig) 75X,Y = np.meshgrid(x_list, y_list) 76 77ax.set_xlabel("x") 78ax.set_ylabel("y") 79ax.set_zlabel("probability") 80 81 82ax.set_xlim(2*m,0) 83ax.set_ylim(0,2*m) 84ax.set_zlim(0,0.015) 85 86#Z軸の色を設定 87offset = p_map.ravel() + np.abs(p_map.min()) 88fracs = offset.astype(float)/offset.max() 89norm = colors.Normalize(fracs.min(), fracs.max()) 90clrs = cm.cool(norm(fracs)) 91 92ax.bar3d(X.ravel(), Y.ravel(), p_map.ravel() ,0.1, 0.1, -p_map.ravel(),color =clrs)#,cmap=cm.hot) 93ax.w_xaxis.set_pane_color((0, 0, 0, 0)) 94ax.w_yaxis.set_pane_color((0, 0, 0, 0)) 95ax.w_zaxis.set_pane_color((0, 0, 0, 0)) 96#plt.show() 97 98
言葉足らずで申し訳ないです。コードをすべて掲載しました。
以下のような、座標にそれぞれの確率の値がでます。この時間、座標に依存した、特定の箇所の確率だけを取り出したい場合、def関数、prob(t,x,y)として行うと良いのでしょうか？
（すいません、これで良いでしょうか？）
上図は一例です（コードの結果とは少し違うかもしれません）。
「欲しい値」がいったい何なのかを明確に述べてないところが伝わらない原因ではないでしょうか？t=0からt=Nまでの全ての確率がほしいのですか？t=Nの確率がほしい？それともt=0,1のそれぞれを順番にほしいのですか？明確に述べてください。
https://teratail.com/questions/126024 の続きですね！コード全部載せたほうが回答を得やすいと思います（動かせるから）。
すいません、説明足らずで、コードと一例の結果を掲載しました。
上のコードではtが0~10秒間に正方格子上を物体が動いた時の、各位置にいる確率分布を計算しています。この0～10秒間計算した状態で、例えば、t=5秒の時の(1,5)の座標の確率が欲しい！みたいな自分が欲しいと思った特定の座標の確率を抽出したいということです。
回答2件
0
ベストアンサー
Pythonの関数として「一度tに関する全てのx, yの確率を計算したらprob(t, x, y)で個々の確率を一つ返せるようにする」ということなら
Python
1g_t = None 2g_p = None 3 4def prob(t, x, y): 5 if global_t != t: 6 global g_t, g_p 7 g_t = t 8 g_p = ... # t時点での全てのx, yの確率を2次元配列で求める計算 9 return g_p[x, y]
とすれば何度も同じ計算をせずに済みますね。しかし・・・こんな関数必要なんでしょうか？
prob(t)で時刻tにおける確率の2次元配列を返すようにすればプロットできますので必要ないと思います。
Python
1n, w = 20, 30 2xs = np.array(range(w)) 3ys = np.array(range(w)) 4p0 = np.zeros((w, w)) 5p0[w//2, w//2] = 1 6... 7 8def prob(n): 9 p_map = p0 10 for _ in range(n): 11 p_map = ... # 次の時刻の確率を直前の確率から計算する 12 return p_map # 時刻tにおけるw*wの確率の2次元配列を返す 13 14Z = prob(T) 15X, Y = np.meshgrid(xs, ys) 16... 17ax.bar3d(X.ravel(), Y.ravel(), Z.ravel(), ...)
でよいのではないでしょうか？またさらにt=1, 2, ..., nまで連続的にプロットするなら
Python
1... 2 3def prob_next(p_map): 4 p_map = ... # 次の時刻の確率の2次元配列を直前の確率の2次元配列から計算する 5 return p_map 6 7Z = p0 # 失礼しました。最初のコードではfor文の中に入れてました orz 8for _ in range(n): 9 Z = prob_next(Z) 10 X, Y = np.meshgrid(xs, ys) 11 ... 12 ax.bar3d(X.ravel(), Y.ravel(), Z.ravel(), ...) 13 plt.show() # あるいはアニメーションにするとか
にしておくと無駄な計算をせずに済むと思います。
投稿2018/05/22 04:05
編集2018/05/22 06:13総合スコア18394
ありがとうございます！
2018/05/22 23:30の質問コメントですが、つまりは一旦t=10を計算済みでもいつでも任意のtの確率が求めたいということですよね？そういうことならprobの定義で参照するglobal変数を一切なくすか、参照しているglobalを変更しないようなプログラムにすればよいという話だと思います。
はい、そういうことです。
アドバイスありがとうございます、ちょっとやってみます！
0
python
1def prob(t,x,y): 2 for t in range(0,n+1):
tを引き数にしたはいいのですが、いきなりfor文の変数として書き換えられています。
x,yも下で書き換わっています。
もっと関数を分解したいということですか？
python
1def evolve(): 2 pass 3 #get t+1 with t 4 5def get_evolved_history(): 6 pass 7 #for t in T 8 # ans.append(evolve()) 9 10def prob(): 11 pass 12 #hist = get_evolved_history() 13 #return hist(t,x,y)
投稿2018/05/21 21:37
総合スコア8560
for文で行う、確率計算をすべてdef関数でまとめ上げてしまって、defの引数を指定すれば、特定の欲しい確率だけを取り出せないかな～と思いました。
history[t,x,y,:]という4つ足のテンソルにデータを収めれば関数をわざわざ定義する必要がありません。
https://teratail.com/questions/125128
の回答に1次元の例があるので、それを2次元に拡張するのはいかがですか？
ありがとうございます。やってみます
あなたの回答
tips
プレビュー
