質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

ただいまの
回答率

90.35%

  • Python 3.x

    7319questions

    Python 3はPythonプログラミング言語の最新バージョンであり、2008年12月3日にリリースされました。

forで作り上げる計算を関数として定義し、ほしい値だけを取り出すには?

解決済

回答 2

投稿 編集

  • 評価
  • クリップ 1
  • VIEW 215

Fallout_18

score 67

以下の時間発展の計算をひとまとめに関数として定義して、欲しい値だけを抽出するには、どのようにすれば良いのでしょうか?
(少しdef関数に苦手意識あり)
一部のコードを掲載しますが、要望があれば全部掲載します。

import numpy as np
import matplotlib.pyplot as plt
import math
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import matplotlib.colors as colors

n=10  #tの範囲
m=15  #偶数

p=1/2
q=1-p

#quantum coin
P = [[-p, q, math.sqrt(p*q), math.sqrt(p*q)],[0,0,0,0],[0,0,0,0],[0,0,0,0]]
Q = [[0,0,0,0],[q, -p, math.sqrt(p*q), math.sqrt(p*q)],[0,0,0,0],[0,0,0,0]]
R = [[0,0,0,0],[0,0,0,0],[math.sqrt(p*q), math.sqrt(p*q), -q, p],[0,0,0,0]]
S = [[0,0,0,0],[0,0,0,0],[0,0,0,0],[math.sqrt(p*q), math.sqrt(p*q), p, -q]]

t_list = []
x_list = []
y_list = []

phi_map = np.zeros((2*m+1, 2*m+1,4)) #縦横2m+1の座標平面上にいる波動状態をセット np.zeros((行,列,[]の中身の数))
phi_map[0,0]= np.array([1/2,1/2,-1/2,-1/2])
p_map=np.zeros([2*m+1,2*m+1])

for i in range(0,2*m+1):
    p = np.dot(phi_map[i,i], np.conj(phi_map[i,i]))
    p_map[i,i]=p
    x_list.append(i)
    y_list.append(i)

def prob(t,x,y):
    for t in range(0,n+1):
        t_list.append(t)
        if t == 0:
            phi_map
            p_map
        else:
            next_phi_map = np.zeros((2*m+1,2*m+1, 4))
            for x in range(0,2*m+1):
                if x == 0:
                    for y in range(0,2*m+1):
                        if y == 0:
                            next_phi_map[x,y] = np.array([np.inner(P, phi_map[x+1,y]) + np.inner(R, phi_map[x,y+1])])
                        elif y == 2*m:
                            next_phi_map[x,y] = np.array([np.inner(P, phi_map[x+1,y]) + np.inner(S, phi_map[x,y-1])])
                        else:
                            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])])
                elif x == 2*m:
                    for y in range(0,2*m+1):
                        if y == 0:
                            next_phi_map[x,y] = np.array([np.inner(Q, phi_map[x-1,y]) + np.inner(R, phi_map[x,y+1])])
                        elif y == 2*m:
                            next_phi_map[x,y] = np.array([np.inner(Q, phi_map[x-1,y]) + np.inner(S, phi_map[x,y-1])])
                        else:
                            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])])
                else:
                    for y in range(0,2*m+1):
                        if y == 0:
                            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])])
                        elif y == 2*m:
                            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])])
                        else:
                            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])])
                        p_map[x,y] = np.dot(next_phi_map[x,y], np.conj(next_phi_map[x,y]))
            phi_map = next_phi_map

        return (p_map[x,y])#←ここは違うので無視してください。


fig = plt.figure()
ax = Axes3D(fig)
X,Y = np.meshgrid(x_list, y_list)

ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("probability")


ax.set_xlim(2*m,0)
ax.set_ylim(0,2*m)
ax.set_zlim(0,0.015)

#Z軸の色を設定
offset = p_map.ravel() + np.abs(p_map.min())
fracs = offset.astype(float)/offset.max()
norm = colors.Normalize(fracs.min(), fracs.max())
clrs = cm.cool(norm(fracs))

ax.bar3d(X.ravel(), Y.ravel(), p_map.ravel() ,0.1, 0.1, -p_map.ravel(),color =clrs)#,cmap=cm.hot)
ax.w_xaxis.set_pane_color((0, 0, 0, 0))
ax.w_yaxis.set_pane_color((0, 0, 0, 0))
ax.w_zaxis.set_pane_color((0, 0, 0, 0))
#plt.show()


言葉足らずで申し訳ないです。コードをすべて掲載しました。
以下のような、座標にそれぞれの確率の値がでます。この時間、座標に依存した、特定の箇所の確率だけを取り出したい場合、def関数、prob(t,x,y)として行うと良いのでしょうか?
(すいません、これで良いでしょうか?)
イメージ説明
上図は一例です(コードの結果とは少し違うかもしれません)。

  • 気になる質問をクリップする

    クリップした質問は、後からいつでもマイページで確認できます。

    またクリップした質問に回答があった際、通知やメールを受け取ることができます。

    クリップを取り消します

  • 良い質問の評価を上げる

    以下のような質問は評価を上げましょう

    • 質問内容が明確
    • 自分も答えを知りたい
    • 質問者以外のユーザにも役立つ

    評価が高い質問は、TOPページの「注目」タブのフィードに表示されやすくなります。

    質問の評価を上げたことを取り消します

  • 評価を下げられる数の上限に達しました

    評価を下げることができません

    • 1日5回まで評価を下げられます
    • 1日に1ユーザに対して2回まで評価を下げられます

    質問の評価を下げる

    teratailでは下記のような質問を「具体的に困っていることがない質問」、「サイトポリシーに違反する質問」と定義し、推奨していません。

    • プログラミングに関係のない質問
    • やってほしいことだけを記載した丸投げの質問
    • 問題・課題が含まれていない質問
    • 意図的に内容が抹消された質問
    • 広告と受け取られるような投稿

    評価が下がると、TOPページの「アクティブ」「注目」タブのフィードに表示されにくくなります。

    質問の評価を下げたことを取り消します

    この機能は開放されていません

    評価を下げる条件を満たしてません

    評価を下げる理由を選択してください

    詳細な説明はこちら

    上記に当てはまらず、質問内容が明確になっていない質問には「情報の追加・修正依頼」機能からコメントをしてください。

    質問の評価を下げる機能の利用条件

    この機能を利用するためには、以下の事項を行う必要があります。

質問への追記・修正、ベストアンサー選択の依頼

  • set0gut1

    2018/05/22 00:55

    https://teratail.com/questions/126024 の続きですね!コード全部載せたほうが回答を得やすいと思います(動かせるから)。

    キャンセル

  • Fallout_18

    2018/05/22 12:15

    すいません、説明足らずで、コードと一例の結果を掲載しました。

    キャンセル

  • Fallout_18

    2018/05/22 23:30

    上のコードではtが0~10秒間に正方格子上を物体が動いた時の、各位置にいる確率分布を計算しています。この0~10秒間計算した状態で、例えば、t=5秒の時の(1,5)の座標の確率が欲しい!みたいな自分が欲しいと思った特定の座標の確率を抽出したいということです。

    キャンセル

回答 2

checkベストアンサー

+2

Pythonの関数として「一度tに関する全てのx, yの確率を計算したらprob(t, x, y)で個々の確率を一つ返せるようにする」ということなら

g_t = None
g_p = None

def prob(t, x, y):
    if global_t != t:
        global g_t, g_p
        g_t = t
        g_p = ... # t時点での全てのx, yの確率を2次元配列で求める計算
    return g_p[x, y]


とすれば何度も同じ計算をせずに済みますね。しかし・・・こんな関数必要なんでしょうか?
prob(t)で時刻tにおける確率の2次元配列を返すようにすればプロットできますので必要ないと思います。

n, w = 20, 30
xs = np.array(range(w))
ys = np.array(range(w))
p0 = np.zeros((w, w))
p0[w//2, w//2] = 1
...

def prob(n):
     p_map = p0
     for _ in range(n):
         p_map = ...  # 次の時刻の確率を直前の確率から計算する
     return p_map  # 時刻tにおけるw*wの確率の2次元配列を返す

Z = prob(T)
X, Y = np.meshgrid(xs, ys)
...
ax.bar3d(X.ravel(), Y.ravel(), Z.ravel(), ...)


でよいのではないでしょうか?またさらにt=1, 2, ..., nまで連続的にプロットするなら

...

def prob_next(p_map):
    p_map = ...  # 次の時刻の確率の2次元配列を直前の確率の2次元配列から計算する
    return p_map

Z = p0  # 失礼しました。最初のコードではfor文の中に入れてました orz
for _ in range(n):
    Z = prob_next(Z)
    X, Y = np.meshgrid(xs, ys)
    ...
    ax.bar3d(X.ravel(), Y.ravel(), Z.ravel(), ...)
    plt.show()   # あるいはアニメーションにするとか


にしておくと無駄な計算をせずに済むと思います。

投稿

編集

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

  • 2018/05/22 23:33

    ありがとうございます!

    キャンセル

  • 2018/05/22 23:38

    2018/05/22 23:30の質問コメントですが、つまりは一旦t=10を計算済みでもいつでも任意のtの確率が求めたいということですよね?そういうことならprobの定義で参照するglobal変数を一切なくすか、参照しているglobalを変更しないようなプログラムにすればよいという話だと思います。

    キャンセル

  • 2018/05/22 23:41

    はい、そういうことです。
    アドバイスありがとうございます、ちょっとやってみます!

    キャンセル

+2

def prob(t,x,y):
    for t in range(0,n+1):

tを引き数にしたはいいのですが、いきなりfor文の変数として書き換えられています。
x,yも下で書き換わっています。


もっと関数を分解したいということですか?

def evolve():
    pass
    #get t+1 with t

def get_evolved_history():
    pass
    #for t in T
    #    ans.append(evolve())

def prob():
    pass
    #hist = get_evolved_history()
    #return hist(t,x,y)

投稿

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

  • 2018/05/22 12:18

    for文で行う、確率計算をすべてdef関数でまとめ上げてしまって、defの引数を指定すれば、特定の欲しい確率だけを取り出せないかな~と思いました。

    キャンセル

  • 2018/05/22 18:13

    history[t,x,y,:]という4つ足のテンソルにデータを収めれば関数をわざわざ定義する必要がありません。

    https://teratail.com/questions/125128
    の回答に1次元の例があるので、それを2次元に拡張するのはいかがですか?

    キャンセル

  • 2018/05/22 23:33

    ありがとうございます。やってみます

    キャンセル

15分調べてもわからないことは、teratailで質問しよう!

  • ただいまの回答率 90.35%
  • 質問をまとめることで、思考を整理して素早く解決
  • テンプレート機能で、簡単に質問をまとめられる

同じタグがついた質問を見る

  • Python 3.x

    7319questions

    Python 3はPythonプログラミング言語の最新バージョンであり、2008年12月3日にリリースされました。