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

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

新規登録して質問してみよう
ただいま回答率
85.50%
Python 3.x

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

Q&A

解決済

1回答

1275閲覧

python TypeErrorが解決できない

astromelt0416

総合スコア15

Python 3.x

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

0グッド

0クリップ

投稿2018/10/19 07:01

編集2018/10/24 13:57

仕様

コードにもコメント文で書いていますが、以下に具体的に書きたいと思います。
メトロポリス法という物理学の中の統計力学におけるシミュレーションをを行う際の計算方法です。

20181024 追記メトロポリスアルゴリズムにミスがあったため、全面的に書きなおしました。

初期状態の形成

(1)Nx*Nyの行列を生成し、各成分を(i,j)であらわす。

(2)行列の各要素は0~2*np.piの乱数を取るものとする。この値をS[i,j]とする。

(3)E[i,j]を以下のように定義し、すべてのS[i,j]に対して計算を行う。
E[i,j] =cos(S(i,j)-S(i+1,j))-cos(S(i,j)-S(i-1,j))-cos(S(i,j)-S(i,j+1))-cos(S(i,j)-S(i,j-1)))
しかし、このままでは行列の端では計算ができなくなる(例えばE[1.1]を計算するとき、S[0,1]やS[1,0]が定義されていないので計算ができない)ので、以下のような周期的境界条件を課します。
~~周期的境界条件~~
例えばNx=8,Ny=8の場合、S[0,1]=S[8,1]、S[1,0]=S[1,8]というように、行列の各端で値を折り返します。

(4)行列Sのある要素S(i,j)をランダムに抽出し、その要素だけを0~2πの乱数で変動させる。

(5)変動させたS(i,j)をSal(i,j)として、以下の計算を行い、Eal(i,j)を計算する。

つまり、変動させた(i,j)以外の要素はそのままにして、手順(3)で行ったものと同じ計算をするということです。この時、行列Sの要素S[i,j]の変動に伴い、行列Eの要素E[i,j],E[i+1,j],E[i-1,j],E[i,j+1],E[1,j-1]も変動することに注意してください。

(6)変動前の行列Eおよび変動後の行列Eの行列式をそれぞれE、Ealとします。
(6-1)E > Eal の場合、要素S[i,j]→Sal[i,j]への更新に成功したとし、(7)以降の過程へ。

(6-2) E < Eal の場合、deltaE=Eal-Eを計算します。そして確率exp(-deltaE/0.001)でもともとの要素S[i,j]→Sal[i,j]への更新に成功したとします。ここで(6)~(6-2)においてS[i,j]の更新に成功した場合は行列E[i,j],E[i+1,j],E[i-1,j],E[i,j+1],E[1,j-1]の要素も当然変化します。

(7)(4)~(6-2)を1回の試行として、1000回繰り返す。

(8)(7)の後(4)~(6-2)を10000回繰り返す。しかし、この時1回ごとにS_sum = np.sum(S)、E_sum = np.sun(E)を計算し、足し合わせていく。つまりi回目の更新で得られたS_sumをS_sum(i)とすると、S_end(10000)=S_sum(1)+S_sum(2)+S_sum(3)...+S_sum(10000)を求めるということです。E_end(10000)についても同様です。

(9)(8)で求めた値を用いて、S_ave = S_end(10000)/10000,E_ave = E_end(10000)/10000 を出力できればOKです。

実際に書いてみたコード

python

1#!/usr/bin/env python3 2# -*- coding: utf-8 -*- 3 4 5import numpy as np 6from random import random,randrange 7 8Nx = int(input("Nx=")) 9Ny = int(input("Ny=")) 10 11np.set_printoptions(linewidth=200) 12 13def Initial(): #ランダムな初期状態の形成 14 S = np.random.uniform(0,2*np.pi,(Nx,Ny))#Nx*Nyの行列Sを形成。各要素は0~2πの乱数 15 print(S) 16 17 Einitial = np.empty_like(S)#Sと同じ形の行列Eを作る 18 S_pad = np.pad(S,(1,1),"wrap")#行列Sの両端で折り返す周期的境界条件 19 20 for i, j in np.dstack(np.mgrid[1:Nx+1, 1:Ny+1]).reshape(-1, 2):#Eの各要素(i,j)を計算する 21 Einitial[i - 1, j - 1] = -np.cos(S_pad[i, j] - S_pad[i - 1, j]) \ 22 - np.cos(S_pad[i, j] - S_pad[i + 1, j]) \ 23 - np.cos(S_pad[i, j] - S_pad[i, j + 1]) \ 24 - np.cos(S_pad[i, j] - S_pad[i, j - 1]) 25 26 Einitialsum = np.sum(Einitial) 27 print("Einitial=", Einitial) 28 print("Einitialsum=", Einitialsum) 29 30 31def Metropolis():#メトロポリスアルゴリズム 32 for k in range(10000):#モンテカルロステップ数を10000とする 33 i = randrange(1,Nx-1) 34 j = randrange(1,Ny-1) 35 #ランダムに抽出した要素S(i,j)の値をランダムで決定する 36 Sal = np.random.uniform(0,2*np.pi) 37 38 #変動させたS(i,j)を用いて計算。S(i,j)以外はもともとの値 39 Eal[i - 1, j - 1] = -np.cos(Sal[i, j] - S_pad[i - 1, j]) \ 40 - np.cos(Sal[i, j] - S_pad[i + 1, j]) \ 41 - np.cos(Sal[i, j] - S_pad[i, j + 1]) \ 42 - np.cos(Sal[i, j] - S_pad[i, j - 1]) 43 44 if E[i,j] > Eal[i,j]: 45 S[i,j] = Sal[i,j] 46 E[i,j] = Eal[i,j] 47 else: 48 deltaE[i,j] = Eal[i,j] - E[i,j] 49 if random() < np.exp(-deltaE[i,j]/0.001): 50 S[i,j] = Sal[i,j] 51 E[i,j] = Eal[i,j] 52 53 Send = np.sum(S) 54 Eend = np.sum(E) 55 56 print(Send) 57 print(Eend) 58 59 60 61if __name__ == '__main__': 62 Initial() 63 Metropolis()

output

1Traceback (most recent call last): 2 File "C:/Users/user2/Desktop/2DXYTEST/testtesttesttest9.py", line 63, in <module> 3 Metropolis() 4 File "C:/Users/user2/Desktop/2DXYTEST/testtesttesttest9.py", line 42, in Metropolis 5 - np.cos(Sal[i, j] - S_pad[i, j - 1]) 6TypeError: 'float' object is not subscriptable

type Errorが出るのですが、現状解決方法がわかりません。また、メトロポリスアルゴリズムにおいても、Sal[i,j]やE[i,j]のように書いていますがこれで最後まで実行できるかわかりません。pythonに詳しい方、よろしくお願いします。

参考

以前にも以下のような質問をし、本プログラム作成の際にお世話になりました。
https://teratail.com/questions/152409
https://teratail.com/questions/152518

###追記(elseの分岐以降でエラーが出ています)

python

1def Metropolis(Nx=8, Ny=8, iterations=10000): 2 S = np.random.uniform(0, 2 * np.pi, (Nx, Ny)) # Nx * Ny の行列Sを形成。各要素は0~2πの乱数 3 E = update(S) 4 5 indices = np.dstack(np.mgrid[0:Nx, 0:Ny]).reshape(-1, 2) 6 for k in range(iterations): # モンテカルロステップ数を 10000 とする 7 # ランダムに要素を選択する。 8 index = tuple(indices[np.random.randint(len(indices))]) 9 10 # 乱数を生成する。 11 E_tmp = E.copy() 12 E_tmp[index] = np.random.uniform(0, 2 * np.pi) 13 14 # 計算 15 E_tmp = update(E_tmp) 16 17 # 更新 18 if E[index] > E_tmp[index]: 19 continue 20 else: 21 delta = E_tmp[index] - E[index] 22 if random() < np.exp(-delta / 0.001):#変更部分です 23 continue 24 25 return S, E 26

output

1Traceback (most recent call last): 2 File "C:/Users/user2/Desktop/2DXYTEST/testtesttesttest13.py", line 44, in <module> 3 S, E = Metropolis(Nx=8, Ny=8, iterations=10000) 4 File "C:/Users/user2/Desktop/2DXYTEST/testtesttesttest13.py", line 37, in Metropolis 5 if random() < np.exp(-delta / 0.001): 6TypeError: 'module' object is not callable 7

###追記2(#更新部分について、教えていただいた内容をもとに条件分岐を訂正しました)

python

1def Metropolis(Nx=8, Ny=8, iterations=10000): 2 S = np.random.uniform(0, 2 * np.pi, (Nx, Ny)) # Nx * Ny の行列Sを形成。各要素は0~2πの乱数 3 E = update(S) 4 5 indices = np.dstack(np.mgrid[0:Nx, 0:Ny]).reshape(-1, 2) 6 for k in range(iterations): # モンテカルロステップ数を 10000 とする 7 # ランダムに要素を選択する。 8 index = tuple(indices[np.random.randint(len(indices))]) 9 10 # 乱数を生成する。 11 E_tmp = E.copy() 12 E_tmp[index] = np.random.uniform(0, 2 * np.pi) 13 14 # 計算 15 E_tmp = update(E_tmp) 16 17 # 更新 18 if E[index] > E_tmp[index]: 19 E[index] = E_tmp[index] 20 21 else: 22 delta = E_tmp[index] - E[index] 23 if np.random.random() < np.exp(-delta / 0.001): 24 E[index] = E_tmp[index] 25 26 27 return S, E

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

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

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

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

hayataka2049

2018/10/19 07:03

質問のタイトルが不適当なのでTypeErrorに変えておいてください
guest

回答1

0

ベストアンサー

#ランダムに抽出した要素S(i,j)の値をランダムで決定する Sal = np.random.uniform(0,2*np.pi)

これは [0, 2pi) に従う乱数を1つ返す関数です。
float の変数に対して、Sal[i, j] のようにしているのでエラーになっています。

上記のアルゴリズムに従うなら、以下のような感じですかね。

import numpy as np def calc_Eal(S, Sal): E = np.empty_like(S) Nx, Ny = S.shape S_pad = np.pad(S, (1, 1), "wrap") # 周期的境界条件 Sal_pad = np.pad(Sal, (1, 1), "wrap") # 周期的境界条件 # Eal(i,j)=cos(Sal(i,j)-S(i+1,j))-cos(Sal(i,j)-S(i-1,j))-cos(Sal(i,j)-S(i,j+1))-cos(Sal(i,j)-S(i,j-1)) # E を計算するときは、S = Sal for i, j in np.dstack(np.mgrid[1:Nx + 1, 1:Ny + 1]).reshape(-1, 2): # Eの各要素(i,j)を計算する E[i - 1, j - 1] = -np.cos(Sal_pad[i, j] - S_pad[i - 1, j]) \ - np.cos(Sal_pad[i, j] - S_pad[i + 1, j]) \ - np.cos(Sal_pad[i, j] - S_pad[i, j + 1]) \ - np.cos(Sal_pad[i, j] - S_pad[i, j - 1]) return E def calc_E(S): return calc_Eal(S, S) def Metropolis(Nx=8, Ny=8, iterations=100): # (1)Nx * Ny の行列を生成し、各成分を (i,j) であらわす。 # (2)行列の各要素は 0 ~ 2 * np.pi の乱数を取るものとする。この値を S[i,j] とする。 S = np.random.uniform(0, 2 * np.pi, (Nx, Ny)) # (3) E[i,j]を以下のように定義し、すべてのS[i,j]に対して計算を行う。 E = calc_E(S) indices = np.dstack(np.mgrid[0:Nx, 0:Ny]).reshape(-1, 2) indices = [tuple(idx) for idx in indices] S_sum, E_sum = 0, 0 # (7)(4)~(6)を10000回繰り返すことでおおよそすべての要素S[i,j]がスピン更新を経験するようにします。 for k in range(iterations): # (4)行列Sのある要素S(i,j)をランダムに抽出し、その要素だけを0~2πの乱数で変動させる。 Sal = S.copy() index = indices[np.random.randint(len(indices))] Sal[index] = np.random.uniform(0, 2 * np.pi) #(5)変動させたS(i,j)をSal(i,j)として、以下の計算を行い、Eal(i,j)を計算する。 Eal = calc_Eal(S, Sal) # (6) ここで (3) で計算した E[i, j] と (5) で計算した Eal[i, j] とを比較します。 delta = Eal - E for idx in indices: # 任意の i, j について delta = Eal[idx] - E[idx] if delta < 0: # Eal[idx] < E[idx] S[idx] = Sal[idx] # 更新する elif delta >= 0 and np.random.random() < np.exp(-delta / 0.001): # Eal[idx] >= E[idx] かつ np.exp(-delta / 0.001) の確率 S[idx] = Sal[idx] # 更新する # (9) S_sum += S.sum() E_sum += E.sum() S_mean = S_sum / iterations E_mean = E_sum / iterations return S_mean, E_mean Nx, Ny = 8, 8 S_mean, E_mean = Metropolis(Nx, Ny, iterations=10000) # それぞれの行列の行列式の値Ssum=np.sum(S)、Esum=np.sum(E)をもとめることができれば計算終了です。 print('Size: ({}, {}), S_mean: {}, E_mean: {}'.format(Nx, Ny, S_mean, E_mean)) # Size: (8, 8), S_mean: 222.99301060088746, E_mean: 2.9429925289412937

追記

indices には配列 E の各要素のインデックスが入っています。

indices = np.array([[0, 0], [0, 1], [0, 2], [0, 3], [0, 4], [0, 5], [0, 6], [0, 7]]) # 実際はもっとあるが、一部だけ

配列の長さを len() で取得すると、8 となります。

print(len(indices)) # 8

また、indices の4つ目の値を取得するには、以下のようになるかと思います。

print(indices[4]) # [0 4]

なので、[0, 8) の範囲の乱数を生成して、そのインデックスの indices の要素を得ることで、
E の要素をランダムに選ぶことになります。

index_to_pick = np.random.randint(len(indices)) # [0, 8) の範囲の乱数を生成 print(index_to_pick) # 2 print(indices[index_to_pick]) # [0 2]

投稿2018/10/19 07:11

編集2018/10/24 11:42
tiitoi

総合スコア21954

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

astromelt0416

2018/10/19 09:42

tiitoi様 授業の都合で返信が遅くなってしまいました。申し訳ございません。 返信に載せていただいたコードで予期していた結果がしっかり出ましたので、大変助かりました。一つ質問させていただきたいです。 index = tuple(indices[np.random.randint(len(indices))]) の箇所でランダムに要素を選択しているとのことですが、len()を使ってどうしてこのような結果になるのでしょうか。len()はリストのサイズ取得ぐらいしか役割を知らなかったので、教えていただけると嬉しいです。
tiitoi

2018/10/19 09:52

追記しました。len() は配列の長さを得ているだけです。配列からランダムに要素を取得するために [0, len(array)) の範囲で乱数を作り、その乱数にあたる配列の要素を取得しています。
astromelt0416

2018/10/19 10:12

上記の内容に関しては理解できました。ありがとうございます。そして、申し訳ないのですが今再度tiitoi様のコードを見返していたら1か所気になるところがありました。♯更新のところでelse以下なのですが、私の書いた条件式と少し異なっているようです。この部分を個人的に修正して試してみたのですがエラーが出てしまいました。もしよろしければ見ていただけますと幸いです。(コードは本文に追記します。)
tiitoi

2018/10/19 10:42

random() は random モジュールの関数なので、 random.random() または numpy の np.random.random() を使うと [0, 1) の乱数が生成されます。 これでどうでしょうか?
tiitoi

2018/10/19 10:45

ちなみに random.random() > 0.6 の場合更新すると今なっていますが、確率 0.6 で更新するといった場合、 if random.random() < 0.6 更新 else continue のような気がするのですが、どうでしょうか?
astromelt0416

2018/10/19 11:02

よくよく考えてみるとtiitoi様の言っていることが正しい気がします。確かにelse以下がそうなる気がします。私が条件分岐でのcontinueを理解できていませんでした。しかし、そう考えると if E[index] > E_tmp[index]: → continue ではなく、ここでも if E[index] > E_tmp[index]: → 更新 となる気がします。
tiitoi

2018/10/19 11:12

回答を修正しました。 E[i,j] > Eal[i,j]、つまり delta = E_tmp[index] - E[index] < 0 の場合も更新しない認識だったのですが、そうではないのでしょうか?
astromelt0416

2018/10/19 11:21

私の書き方が悪かったところがあるかもしれないので再度スピン更新の箇所だけ訂正します。 まずE[i,j] < Eal[i,j]ならば更新します。E[i,j] > Eal[i,j]の場合のみdelta = E_tmp[index] - E[index](>0)を計算し確率exp(-delta/0.001))で更新します。 つまりE[i,j] > Eal[i,j]でもある確率で更新はします。戸惑わさせてしまい申し訳ございません。
tiitoi

2018/10/19 11:40

* E[i,j] > Eal[i,j] の場合、更新する。 * E[i,j] <= Eal[i,j] の場合、確率 exp(-(E_tmp[index] - E[index])/0.001)) で更新する。 つまり、dleta = E_tmp[index] - E[index] とおくと、 * dleta < 0 の場合、更新する。 * delta >= 0 の場合、確率 exp(-delta/0.001)) で更新する。 ということですね。回答を修正しました。
astromelt0416

2018/10/19 12:06

はい、更新に関しましてはそういうことになります。修正までしていただきありがとうございました。tiitoi様のコードでEの正しい行列が得られることを確認しました。本当にありがとうございます。 以下は急ぎではないので、お時間があるときに少し見ていただきたいです。 本プログラムにおいて目標はスピンを更新後のE及びSの行列を得ることです。本来、行列S[i,j]→Sal[i,j]を仮定し、そこから導かれるEal[i,j]とE[i,j]の比較により更新を行うというものでした。したがって最後の段階でreturnされるべきなのは更新後のEとSなのです。しかし、結果を見るあたりreturnされているSが更新前のSであるような気がします。このあたりをお時間があるときに考えていただければ幸いです。よろしくお願い致します。
tiitoi

2018/10/19 16:44

よく読んだら記載のアルゴリズム通りになっていないところがあったので、直しました。 ただアルゴリズム内でEの更新がかかれてなかったのですが、更新する必要はないのでしょうか?
astromelt0416

2018/10/19 22:24

回答ありがとうございます。 10000回の試行の末にまんべんなく更新された行列Sを用いて最後にEを計算しreturnするという考えでした。行列Sの要素S[i,j]の更新のために適宜E「i,j] の値は計算する必要はありますが、その都度Eの更新はする必要はないと思います。
astromelt0416

2018/10/19 22:34

質問文におきまして、”仕様”の部分に追記しました。
tiitoi

2018/10/20 04:21 編集

質問欄にアルゴリズムの手順の該当する箇所にコメントを入れました。 これまで理解している通りに計算した結果以下のようになりましたが、質問欄の値にはなりませんでした。 # Size: (8, 8), S.sum(): 148.1828022790117, E.sum(): -108.97739996283 気になる点としては、 ・「ランダムに抽出し、その要素だけを乱数で変動させる」というのはその要素を乱数の値で置き換えるという解釈でよいのか。 ・E[i,j]とEal[i,j]を比較して、E[i,j]を更新する部分はランダムに置き換えたインデックスのみ対象に行われるという解釈でよいのか。 ということですかね。 おそらく自分の質問欄の解釈が間違っているため、正しい値が出ないのだと思いますが、メトロポリス法自体を理解していないため、こちらでどこが問題なのか見つけるのは難しいです。すみません。 回答のコードを参考に解釈が間違っている箇所は修正してみてください。
astromelt0416

2018/10/24 08:58

tiitoi様 あれからコードの検証を行いあらゆることを試したのですがやはりうまくいきませんでした。tiitoi様の解釈に誤りが見られなかったのでもしやと思い、英語の文献を再度漁ったところアルゴリズムにいくつかミスがあることに気付きました。今更ですが質問文を改訂しておいたのでもしお時間あれば見ていただけると助かります。
tiitoi

2018/10/24 11:45

解答欄に記載の通りにコードを修正して試してみましたが、質問欄に記載していただいたような値にはならないですね。 E[i,j]とEal[i,j]を比較して、E[i,j]を更新する部分は、各成分ごとに比較して、更新という認識でいいのですかね?
astromelt0416

2018/10/24 14:01

tiitoi様 ずっと行列Sの成分S[i,j]を更新した際に変化するのはE[i,j]だけだと勘違いしておりました。正しく質問投稿文の方も書き直しました。比較するのは行列Eの成分ではなく、行列Eの行列式の値E=np.sum(E)同士になるとのことです。何度も申し訳ございません。
tiitoi

2018/10/24 14:58 編集

> 行列Eの行列式の値E=np.sum(E)同士 すいません。ちょっと意味がわからないです。 np.sum(E) はEのすべての要素の和が計算されますが、それを比較するとはどういうことでしょうか。 > 行列E[i,j],E[i+1,j],E[i-1,j],E[i,j+1],E[1,j-1]の要素も当然変化します。 ここもどういうことでしょうか。以前の話だとループ中はE自体は変化しないという話だったと思うのですが、そうではなかったということでしょうか? ------------------------------------- 申し訳ありませんが、メトロポリス法自体を知らないため、質問欄に記載のアルゴリズムをこちらで解釈して記述したコードでなにか計算値が算出されたとして、それが合っているのかどうかこちらではわからないので、恐縮ですが、正しいアルゴリズムにコードをこちらで修正するのは限界があります。 アルゴリズムの解説等が載っているサイトはありますか? 質問者さんの方で正しくコードを修正するのは難しいのでしょうか?
guest

あなたの回答

tips

太字

斜体

打ち消し線

見出し

引用テキストの挿入

コードの挿入

リンクの挿入

リストの挿入

番号リストの挿入

表の挿入

水平線の挿入

プレビュー

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

ただいまの回答率
85.50%

質問をまとめることで
思考を整理して素早く解決

テンプレート機能で
簡単に質問をまとめる

質問する

関連した質問