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

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

ただいまの
回答率

87.95%

python TypeErrorが解決できない

解決済

回答 1

投稿 編集

  • 評価
  • クリップ 0
  • VIEW 669

score 15

 仕様

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

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です。 

 実際に書いてみたコード

#!/usr/bin/env python3
# -*- coding: utf-8 -*-


import numpy as np
from random import random,randrange

Nx = int(input("Nx="))
Ny = int(input("Ny="))

np.set_printoptions(linewidth=200)

def Initial(): #ランダムな初期状態の形成
    S = np.random.uniform(0,2*np.pi,(Nx,Ny))#Nx*Nyの行列Sを形成。各要素は0~2πの乱数
    print(S)

    Einitial = np.empty_like(S)#Sと同じ形の行列Eを作る
    S_pad = np.pad(S,(1,1),"wrap")#行列Sの両端で折り返す周期的境界条件

    for i, j in np.dstack(np.mgrid[1:Nx+1, 1:Ny+1]).reshape(-1, 2):#Eの各要素(i,j)を計算する
        Einitial[i - 1, j - 1] = -np.cos(S_pad[i, j] - S_pad[i - 1, j]) \
                          - np.cos(S_pad[i, j] - S_pad[i + 1, j]) \
                          - np.cos(S_pad[i, j] - S_pad[i, j + 1]) \
                          - np.cos(S_pad[i, j] - S_pad[i, j - 1])

    Einitialsum = np.sum(Einitial)
    print("Einitial=", Einitial)
    print("Einitialsum=", Einitialsum)


def Metropolis():#メトロポリスアルゴリズム
    for k in range(10000):#モンテカルロステップ数を10000とする
        i = randrange(1,Nx-1)
        j = randrange(1,Ny-1)
        #ランダムに抽出した要素S(i,j)の値をランダムで決定する
        Sal = np.random.uniform(0,2*np.pi)

        #変動させたS(i,j)を用いて計算。S(i,j)以外はもともとの値
        Eal[i - 1, j - 1] = -np.cos(Sal[i, j] - S_pad[i - 1, j]) \
                            - np.cos(Sal[i, j] - S_pad[i + 1, j]) \
                            - np.cos(Sal[i, j] - S_pad[i, j + 1]) \
                            - np.cos(Sal[i, j] - S_pad[i, j - 1])

        if E[i,j] > Eal[i,j]:
            S[i,j] = Sal[i,j]
            E[i,j] = Eal[i,j]
        else:
            deltaE[i,j] = Eal[i,j] - E[i,j]
            if random() < np.exp(-deltaE[i,j]/0.001):
                S[i,j] = Sal[i,j]
                E[i,j] = Eal[i,j]

    Send = np.sum(S)
    Eend = np.sum(E)

    print(Send)
    print(Eend)



if __name__ == '__main__':
    Initial()
    Metropolis()


Traceback (most recent call last):
  File "C:/Users/user2/Desktop/2DXYTEST/testtesttesttest9.py", line 63, in <module>
    Metropolis()
  File "C:/Users/user2/Desktop/2DXYTEST/testtesttesttest9.py", line 42, in Metropolis
    - np.cos(Sal[i, j] - S_pad[i, j - 1])
TypeError: '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の分岐以降でエラーが出ています)

def Metropolis(Nx=8, Ny=8, iterations=10000):
    S = np.random.uniform(0, 2 * np.pi, (Nx, Ny)) # Nx * Ny の行列Sを形成。各要素は0~2πの乱数
    E = update(S)

    indices = np.dstack(np.mgrid[0:Nx, 0:Ny]).reshape(-1, 2)
    for k in range(iterations): # モンテカルロステップ数を 10000 とする
        # ランダムに要素を選択する。
        index = tuple(indices[np.random.randint(len(indices))])

        # 乱数を生成する。
        E_tmp = E.copy()
        E_tmp[index] = np.random.uniform(0, 2 * np.pi)

        # 計算
        E_tmp = update(E_tmp)

        # 更新
        if E[index] > E_tmp[index]:
            continue
        else:
            delta = E_tmp[index] - E[index]
            if random() < np.exp(-delta / 0.001):#変更部分です
                continue

    return S, E


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

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

def Metropolis(Nx=8, Ny=8, iterations=10000):
    S = np.random.uniform(0, 2 * np.pi, (Nx, Ny)) # Nx * Ny の行列Sを形成。各要素は0~2πの乱数
    E = update(S)

    indices = np.dstack(np.mgrid[0:Nx, 0:Ny]).reshape(-1, 2)
    for k in range(iterations): # モンテカルロステップ数を 10000 とする
        # ランダムに要素を選択する。
        index = tuple(indices[np.random.randint(len(indices))])

        # 乱数を生成する。
        E_tmp = E.copy()
        E_tmp[index] = np.random.uniform(0, 2 * np.pi)

        # 計算
        E_tmp = update(E_tmp)

        # 更新
        if E[index] > E_tmp[index]:
            E[index] = E_tmp[index]

        else:
            delta = E_tmp[index] - E[index]
            if np.random.random() < np.exp(-delta / 0.001):
                E[index] = E_tmp[index]


    return S, E
  • 気になる質問をクリップする

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

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

    クリップを取り消します

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

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

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

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

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

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

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

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

    質問の評価を下げる

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

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

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

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

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

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

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

    詳細な説明はこちら

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

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

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

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

  • hayataka2049

    2018/10/19 16:03

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

    キャンセル

回答 1

checkベストアンサー

+1

#ランダムに抽出した要素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)行列の各要素は 02 * 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)をランダムに抽出し、その要素だけを02πの乱数で変動させる。
        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/24 20:45

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

    キャンセル

  • 2018/10/24 23:01

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

    キャンセル

  • 2018/10/24 23:54 編集

    > 行列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自体は変化しないという話だったと思うのですが、そうではなかったということでしょうか?

    -------------------------------------

    申し訳ありませんが、メトロポリス法自体を知らないため、質問欄に記載のアルゴリズムをこちらで解釈して記述したコードでなにか計算値が算出されたとして、それが合っているのかどうかこちらではわからないので、恐縮ですが、正しいアルゴリズムにコードをこちらで修正するのは限界があります。
    アルゴリズムの解説等が載っているサイトはありますか?
    質問者さんの方で正しくコードを修正するのは難しいのでしょうか?

    キャンセル

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

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

関連した質問

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