teratail
回答率
85.30%
teratail header banner
teratail header banner
質問するログイン新規登録
トップPython 3.xに関する質問python TypeErrorが解決できない

編集履歴

回答編集履歴

7

d

2018/10/24 11:42

投稿

tiitoi
tiitoi

スコア21960

answer CHANGED
@@ -30,7 +30,7 @@
30
30
 
31
31
  def calc_E(S):
32
32
  return calc_Eal(S, S)
33
-
33
+
34
34
  def Metropolis(Nx=8, Ny=8, iterations=100):
35
35
  # （1）Nx * Ny の行列を生成し、各成分を (i,j) であらわす。
36
36
  # （2）行列の各要素は 0 ～ 2 * np.pi の乱数を取るものとする。この値を S[i,j] とする。
@@ -39,36 +39,43 @@
39
39
  E = calc_E(S)
40
40
 
41
41
  indices = np.dstack(np.mgrid[0:Nx, 0:Ny]).reshape(-1, 2)
42
+ indices = [tuple(idx) for idx in indices]
43
+ S_sum, E_sum = 0, 0
44
+
42
45
  # (7)(4)~(6)を10000回繰り返すことでおおよそすべての要素S[i,j]がスピン更新を経験するようにします。
43
46
  for k in range(iterations):
44
47
 
45
48
  # （4）行列Sのある要素S（i,j)をランダムに抽出し、その要素だけを0～2πの乱数で変動させる。
46
49
  Sal = S.copy()
47
- index = tuple(indices[np.random.randint(len(indices))])
50
+ index = indices[np.random.randint(len(indices))]
48
51
  Sal[index] = np.random.uniform(0, 2 * np.pi)
49
52
 
50
53
  #（5）変動させたS(i,j)をSal(i,j)として、以下の計算を行い、Eal(i,j)を計算する。
51
54
  Eal = calc_Eal(S, Sal)
52
55
 
53
- # (6)ここで(3)で計算したE[i,j]と(5)で計算したEal[i,j]とを比較します。
56
+ # (6) ここで (3) で計算した E[i, j] と (5) で計算した Eal[i, j] とを比較します。
57
+ delta = Eal - E
58
+ for idx in indices: # 任意の i, j について
54
- delta = Eal[index] - E[index]
59
+ delta = Eal[idx] - E[idx]
55
- # (6-1）E[i,j] > Eal[i,j] の場合、もともとの要素S[i,j]をSal[i,j]に更新し、(7)以降の過程へ。
56
- if delta < 0:
60
+ if delta < 0: # Eal[idx] < E[idx]
57
- S[index] = Sal[index] # 更新する
61
+ S[idx] = Sal[idx] # 更新する
58
- # (6-2) E[i,j] < Eal[i,j] の場合、deltaE[i,j]=Eal[i,j]-E[i,j]を計算します。
59
- # そして確率exp(-deltaE[i,j]/0.001)でもともとの要素S[i,j]をSal[i,j]に更新します。
60
- elif delta >= 0 and np.random.random() < np.exp(-delta / 0.001):
62
+ elif delta >= 0 and np.random.random() < np.exp(-delta / 0.001):
63
+ # Eal[idx] >= E[idx] かつ np.exp(-delta / 0.001) の確率
61
- S[index] = Sal[index] # 更新する
64
+ S[idx] = Sal[idx] # 更新する
62
-
63
- # (8)(7)の操作後の行列Sを用いて行列Eを計算する。
64
- E = calc_E(S)
65
- return S, E
66
65
 
66
+ # (9)
67
+ S_sum += S.sum()
68
+ E_sum += E.sum()
69
+
70
+ S_mean = S_sum / iterations
71
+ E_mean = E_sum / iterations
72
+ return S_mean, E_mean
73
+
67
74
  Nx, Ny = 8, 8
68
- S, E = Metropolis(Nx, Ny, iterations=10000)
75
+ S_mean, E_mean = Metropolis(Nx, Ny, iterations=10000)
69
76
  # それぞれの行列の行列式の値Ssum＝np.sum(S)、Esum=np.sum(E)をもとめることができれば計算終了です。
70
- print('Size: ({}, {}), S.sum(): {}, E.sum(): {}'.format(Nx, Ny, S.sum(), E.sum()))
77
+ print('Size: ({}, {}), S_mean: {}, E_mean: {}'.format(Nx, Ny, S_mean, E_mean))
71
- # Size: (8, 8), S.sum(): 148.1828022790117, E.sum(): -108.97739996283
78
+ # Size: (8, 8), S_mean: 222.99301060088746, E_mean: 2.9429925289412937
72
79
  ```
73
80
 
74
81
  ## 追記

6

s

2018/10/24 11:42

投稿

tiitoi
tiitoi

スコア21960

answer CHANGED
@@ -12,12 +12,15 @@
12
12
  ```
13
13
  import numpy as np
14
14
 
15
- def update(S, Sal):
15
+ def calc_Eal(S, Sal):
16
- E = np.empty_like(S) # Sと同じ形の行列Eを作る
16
+ E = np.empty_like(S)
17
17
 
18
18
  Nx, Ny = S.shape
19
- S_pad = np.pad(S, (1, 1), "wrap") # 行列Sの両端で折り返す周期的境界条件
19
+ S_pad = np.pad(S, (1, 1), "wrap") # 周期的境界条件
20
- Sal_pad = np.pad(Sal, (1, 1), "wrap") # 行列Sの両端で折り返す周期的境界条件
20
+ Sal_pad = np.pad(Sal, (1, 1), "wrap") # 周期的境界条件
21
+
22
+ # 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))
23
+ # E を計算するときは、S = Sal
21
24
  for i, j in np.dstack(np.mgrid[1:Nx + 1, 1:Ny + 1]).reshape(-1, 2): # Eの各要素（i,j）を計算する
22
25
  E[i - 1, j - 1] = -np.cos(Sal_pad[i, j] - S_pad[i - 1, j]) \
23
26
  - np.cos(Sal_pad[i, j] - S_pad[i + 1, j]) \
@@ -25,36 +28,47 @@
25
28
  - np.cos(Sal_pad[i, j] - S_pad[i, j - 1])
26
29
  return E
27
30
 
31
+ def calc_E(S):
32
+ return calc_Eal(S, S)
33
+
28
34
  def Metropolis(Nx=8, Ny=8, iterations=100):
35
+ # （1）Nx * Ny の行列を生成し、各成分を (i,j) であらわす。
36
+ # （2）行列の各要素は 0 ～ 2 * np.pi の乱数を取るものとする。この値を S[i,j] とする。
29
- S = np.random.uniform(0, 2 * np.pi, (Nx, Ny)) # Nx * Ny の行列Sを形成。各要素は０～２πの乱数
37
+ S = np.random.uniform(0, 2 * np.pi, (Nx, Ny))
38
+ # (3) E[i,j]を以下のように定義し、すべてのS[i,j]に対して計算を行う。
30
- E = update(S, S)
39
+ E = calc_E(S)
31
40
 
32
41
  indices = np.dstack(np.mgrid[0:Nx, 0:Ny]).reshape(-1, 2)
42
+ # (7)(4)~(6)を10000回繰り返すことでおおよそすべての要素S[i,j]がスピン更新を経験するようにします。
33
- for k in range(iterations): # モンテカルロステップ数を 10000 とする
43
+ for k in range(iterations):
34
44
 
35
- # Sal を作成する。
45
+ # （4）行列Sのある要素S（i,j)をランダムに抽出し、その要素だけを0～2πの乱数で変動させる。
36
46
  Sal = S.copy()
37
- # ランダムに要素を選択し、乱数に置き換える。
38
47
  index = tuple(indices[np.random.randint(len(indices))])
39
48
  Sal[index] = np.random.uniform(0, 2 * np.pi)
40
49
 
41
- # Eal を計算する。
50
+ #（5）変動させたS(i,j)をSal(i,j)として、以下の計算を行い、Eal(i,j)を計算する。
42
- Eal = update(S, Sal)
51
+ Eal = calc_Eal(S, Sal)
43
52
 
44
- # 更新 Eal[index] - E[index]
53
+ # (6)ここで(3)で計算したE[i,j]と(5)で計算したEal[i,j]とを比較します。
45
54
  delta = Eal[index] - E[index]
46
- # Eal[index] < E[index]
55
+ # (6-1）E[i,j] > Eal[i,j] の場合、もともとの要素S[i,j]をSal[i,j]に更新し、(7)以降の過程へ。
47
56
  if delta < 0:
48
57
  S[index] = Sal[index] # 更新する
49
- # Eal[index] >= E[index] かつ確率 np.exp(-delta / 0.001) で更新
58
+ # (6-2) E[i,j] < Eal[i,j] の場合、deltaE[i,j]=Eal[i,j]-E[i,j]を計算します。
59
+ # そして確率exp(-deltaE[i,j]/0.001)でもともとの要素S[i,j]をSal[i,j]に更新します。
50
60
  elif delta >= 0 and np.random.random() < np.exp(-delta / 0.001):
51
61
  S[index] = Sal[index] # 更新する
62
+
52
-
63
+ # (8)(7)の操作後の行列Sを用いて行列Eを計算する。
64
+ E = calc_E(S)
53
65
  return S, E
54
66
 
67
+ Nx, Ny = 8, 8
55
- S, E = Metropolis(Nx=8, Ny=8, iterations=10000)
68
+ S, E = Metropolis(Nx, Ny, iterations=10000)
69
+ # それぞれの行列の行列式の値Ssum＝np.sum(S)、Esum=np.sum(E)をもとめることができれば計算終了です。
56
- print('S', S)
70
+ print('Size: ({}, {}), S.sum(): {}, E.sum(): {}'.format(Nx, Ny, S.sum(), E.sum()))
57
- print('E', E)
71
+ # Size: (8, 8), S.sum(): 148.1828022790117, E.sum(): -108.97739996283
58
72
  ```
59
73
 
60
74
  ## 追記

5

a

2018/10/20 04:14

投稿

tiitoi
tiitoi

スコア21960

answer CHANGED
@@ -12,41 +12,43 @@
12
12
  ```
13
13
  import numpy as np
14
14
 
15
- def update(S):
15
+ def update(S, Sal):
16
+ E = np.empty_like(S) # Sと同じ形の行列Eを作る
17
+
16
18
  Nx, Ny = S.shape
17
- E = np.empty_like(S) # Sと同じ形の行列Eを作る
18
19
  S_pad = np.pad(S, (1, 1), "wrap") # 行列Sの両端で折り返す周期的境界条件
20
+ Sal_pad = np.pad(Sal, (1, 1), "wrap") # 行列Sの両端で折り返す周期的境界条件
19
21
  for i, j in np.dstack(np.mgrid[1:Nx + 1, 1:Ny + 1]).reshape(-1, 2): # Eの各要素（i,j）を計算する
20
- E[i - 1, j - 1] = -np.cos(S_pad[i, j] - S_pad[i - 1, j]) \
22
+ E[i - 1, j - 1] = -np.cos(Sal_pad[i, j] - S_pad[i - 1, j]) \
21
- - np.cos(S_pad[i, j] - S_pad[i + 1, j]) \
23
+ - np.cos(Sal_pad[i, j] - S_pad[i + 1, j]) \
22
- - np.cos(S_pad[i, j] - S_pad[i, j + 1]) \
24
+ - np.cos(Sal_pad[i, j] - S_pad[i, j + 1]) \
23
- - np.cos(S_pad[i, j] - S_pad[i, j - 1])
25
+ - np.cos(Sal_pad[i, j] - S_pad[i, j - 1])
24
26
  return E
25
27
 
26
28
  def Metropolis(Nx=8, Ny=8, iterations=100):
27
29
  S = np.random.uniform(0, 2 * np.pi, (Nx, Ny)) # Nx * Ny の行列Sを形成。各要素は０～２πの乱数
28
- E = update(S)
30
+ E = update(S, S)
29
31
 
30
32
  indices = np.dstack(np.mgrid[0:Nx, 0:Ny]).reshape(-1, 2)
31
33
  for k in range(iterations): # モンテカルロステップ数を 10000 とする
34
+
35
+ # Sal を作成する。
36
+ Sal = S.copy()
32
- # ランダムに要素を選択する。
37
+ # ランダムに要素を選択し、乱数に置き換える。
33
38
  index = tuple(indices[np.random.randint(len(indices))])
34
-
35
- # 乱数を生成する。
36
- E_tmp = E.copy()
37
- E_tmp[index] = np.random.uniform(0, 2 * np.pi)
39
+ Sal[index] = np.random.uniform(0, 2 * np.pi)
38
-
40
+
39
- # 計算
41
+ # Eal を計算する。
40
- E_tmp = update(E_tmp)
42
+ Eal = update(S, Sal)
41
-
43
+
42
- # 更新 E_tmp[index] - E[index]
44
+ # 更新 Eal[index] - E[index]
43
- delta = E_tmp[index] - E[index]
45
+ delta = Eal[index] - E[index]
44
- # E_tmp[index] < E[index]
46
+ # Eal[index] < E[index]
45
47
  if delta < 0:
46
- E[index] = E_tmp[index] # 更新する
48
+ S[index] = Sal[index] # 更新する
47
- # E_tmp[index] <= E[index] かつ確率 np.exp(-delta / 0.001) で更新
49
+ # Eal[index] >= E[index] かつ確率 np.exp(-delta / 0.001) で更新
48
50
  elif delta >= 0 and np.random.random() < np.exp(-delta / 0.001):
49
- E[index] = E_tmp[index] # 更新する
51
+ S[index] = Sal[index] # 更新する
50
52
 
51
53
  return S, E
52
54
 

4

c

2018/10/19 16:40

投稿

tiitoi
tiitoi

スコア21960

answer CHANGED
@@ -11,8 +11,6 @@
11
11
 
12
12
  ```
13
13
  import numpy as np
14
- import random
15
- np.set_printoptions(linewidth=200)
16
14
 
17
15
  def update(S):
18
16
  Nx, Ny = S.shape
@@ -41,13 +39,15 @@
41
39
  # 計算
42
40
  E_tmp = update(E_tmp)
43
41
 
44
- # 更新
42
+ # 更新 E_tmp[index] - E[index]
45
43
  delta = E_tmp[index] - E[index]
44
+ # E_tmp[index] < E[index]
45
+ if delta < 0:
46
+ E[index] = E_tmp[index] # 更新する
47
+ # E_tmp[index] <= E[index] かつ確率 np.exp(-delta / 0.001) で更新
46
- if delta < 0 or np.random.random() > np.exp(-delta / 0.001):
48
+ elif delta >= 0 and np.random.random() < np.exp(-delta / 0.001):
47
- continue # 更新しない
49
+ E[index] = E_tmp[index] # 更新する
48
50
 
49
- E[index] = E_tmp[index] # 更新する
50
-
51
51
  return S, E
52
52
 
53
53
  S, E = Metropolis(Nx=8, Ny=8, iterations=10000)

3

s

2018/10/19 11:37

投稿

tiitoi
tiitoi

スコア21960

answer CHANGED
@@ -42,11 +42,11 @@
42
42
  E_tmp = update(E_tmp)
43
43
 
44
44
  # 更新
45
- if E[index] > E_tmp[index]:
46
- continue
47
- else:
48
- delta = E_tmp[index] - E[index]
45
+ delta = E_tmp[index] - E[index]
46
+ if delta < 0 or np.random.random() > np.exp(-delta / 0.001):
47
+ continue # 更新しない
48
+
49
- E[index] = np.exp(-delta / 0.001)
49
+ E[index] = E_tmp[index] # 更新する
50
50
 
51
51
  return S, E
52
52
 

2

d

2018/10/19 11:06

投稿

tiitoi
tiitoi

スコア21960

answer CHANGED
@@ -53,4 +53,39 @@
53
53
  S, E = Metropolis(Nx=8, Ny=8, iterations=10000)
54
54
  print('S', S)
55
55
  print('E', E)
56
+ ```
57
+
58
+ ## 追記
59
+
60
+ indices には配列 E の各要素のインデックスが入っています。
61
+ ```
62
+ indices = np.array([[0, 0],
63
+ [0, 1],
64
+ [0, 2],
65
+ [0, 3],
66
+ [0, 4],
67
+ [0, 5],
68
+ [0, 6],
69
+ [0, 7]])
70
+ # 実際はもっとあるが、一部だけ
71
+ ```
72
+
73
+ 配列の長さを len() で取得すると、8 となります。
74
+ ```
75
+ print(len(indices)) # 8
76
+ ```
77
+
78
+ また、indices の4つ目の値を取得するには、以下のようになるかと思います。
79
+
80
+ ```
81
+ print(indices[4]) # [0 4]
82
+ ```
83
+
84
+ なので、[0, 8) の範囲の乱数を生成して、そのインデックスの indices の要素を得ることで、
85
+ E の要素をランダムに選ぶことになります。
86
+
87
+ ```
88
+ index_to_pick = np.random.randint(len(indices)) # [0, 8) の範囲の乱数を生成
89
+ print(index_to_pick) # 2
90
+ print(indices[index_to_pick]) # [0 2]
56
91
  ```

1

d

2018/10/19 09:51

投稿

tiitoi
tiitoi

スコア21960

answer CHANGED
@@ -6,8 +6,51 @@
6
6
  これは [0, 2pi) に従う乱数を1つ返す関数です。
7
7
  float の変数に対して、Sal[i, j] のようにしているのでエラーになっています。
8
8
 
9
+
9
- S と同じ形状の乱数を作成する場合は第3引数で size を指定してください。
10
+ 上記のアルゴリズムに従うなら、以下のような感じですかね。
11
+
10
12
  ```
13
+ import numpy as np
14
+ import random
15
+ np.set_printoptions(linewidth=200)
16
+
17
+ def update(S):
18
+ Nx, Ny = S.shape
19
+ E = np.empty_like(S) # Sと同じ形の行列Eを作る
20
+ S_pad = np.pad(S, (1, 1), "wrap") # 行列Sの両端で折り返す周期的境界条件
21
+ for i, j in np.dstack(np.mgrid[1:Nx + 1, 1:Ny + 1]).reshape(-1, 2): # Eの各要素（i,j）を計算する
22
+ E[i - 1, j - 1] = -np.cos(S_pad[i, j] - S_pad[i - 1, j]) \
23
+ - np.cos(S_pad[i, j] - S_pad[i + 1, j]) \
24
+ - np.cos(S_pad[i, j] - S_pad[i, j + 1]) \
25
+ - np.cos(S_pad[i, j] - S_pad[i, j - 1])
26
+ return E
27
+
28
+ def Metropolis(Nx=8, Ny=8, iterations=100):
29
+ S = np.random.uniform(0, 2 * np.pi, (Nx, Ny)) # Nx * Ny の行列Sを形成。各要素は０～２πの乱数
30
+ E = update(S)
31
+
32
+ indices = np.dstack(np.mgrid[0:Nx, 0:Ny]).reshape(-1, 2)
33
+ for k in range(iterations): # モンテカルロステップ数を 10000 とする
11
- #ランダムに抽出した要素S（i,j）の値をランダムで決定する
34
+ # ランダムに要素を選択する。
35
+ index = tuple(indices[np.random.randint(len(indices))])
36
+
37
+ # 乱数を生成する。
38
+ E_tmp = E.copy()
12
- Sal = np.random.uniform(0,2*np.pi, S.shape)
39
+ E_tmp[index] = np.random.uniform(0, 2 * np.pi)
40
+
41
+ # 計算
42
+ E_tmp = update(E_tmp)
43
+
44
+ # 更新
45
+ if E[index] > E_tmp[index]:
46
+ continue
47
+ else:
48
+ delta = E_tmp[index] - E[index]
49
+ E[index] = np.exp(-delta / 0.001)
50
+
51
+ return S, E
52
+
53
+ S, E = Metropolis(Nx=8, Ny=8, iterations=10000)
54
+ print('S', S)
55
+ print('E', E)
13
56
  ```