回答編集履歴

7

d

2018/10/24 11:42

投稿

tiitoi
tiitoi

スコア21956

test CHANGED
@@ -62,24 +62,30 @@
62
62
 
63
63
  return calc_Eal(S, S)
64
64
 
65
+
66
+
67
+ def Metropolis(Nx=8, Ny=8, iterations=100):
68
+
69
+ # (1)Nx * Ny の行列を生成し、各成分を (i,j) であらわす。
70
+
71
+ # (2)行列の各要素は 0 ~ 2 * np.pi の乱数を取るものとする。この値を S[i,j] とする。
72
+
73
+ S = np.random.uniform(0, 2 * np.pi, (Nx, Ny))
74
+
75
+ # (3) E[i,j]を以下のように定義し、すべてのS[i,j]に対して計算を行う。
76
+
77
+ E = calc_E(S)
78
+
79
+
80
+
81
+ indices = np.dstack(np.mgrid[0:Nx, 0:Ny]).reshape(-1, 2)
82
+
83
+ indices = [tuple(idx) for idx in indices]
84
+
85
+ S_sum, E_sum = 0, 0
86
+
65
87
 
66
88
 
67
- def Metropolis(Nx=8, Ny=8, iterations=100):
68
-
69
- # (1)Nx * Ny の行列を生成し、各成分を (i,j) であらわす。
70
-
71
- # (2)行列の各要素は 0 ~ 2 * np.pi の乱数を取るものとする。この値を S[i,j] とする。
72
-
73
- S = np.random.uniform(0, 2 * np.pi, (Nx, Ny))
74
-
75
- # (3) E[i,j]を以下のように定義し、すべてのS[i,j]に対して計算を行う。
76
-
77
- E = calc_E(S)
78
-
79
-
80
-
81
- indices = np.dstack(np.mgrid[0:Nx, 0:Ny]).reshape(-1, 2)
82
-
83
89
  # (7)(4)~(6)を10000回繰り返すことでおおよそすべての要素S[i,j]がスピン更新を経験するようにします。
84
90
 
85
91
  for k in range(iterations):
@@ -90,7 +96,7 @@
90
96
 
91
97
  Sal = S.copy()
92
98
 
93
- index = tuple(indices[np.random.randint(len(indices))])
99
+ index = indices[np.random.randint(len(indices))]
94
100
 
95
101
  Sal[index] = np.random.uniform(0, 2 * np.pi)
96
102
 
@@ -102,43 +108,51 @@
102
108
 
103
109
 
104
110
 
105
- # (6)ここで(3)で計算したE[i,j]と(5)で計算したEal[i,j]とを比較します。
111
+ # (6) ここで (3) で計算した E[i, j] (5) で計算した Eal[i, j] とを比較します。
112
+
106
-
113
+ delta = Eal - E
114
+
115
+ for idx in indices: # 任意の i, j について
116
+
107
- delta = Eal[index] - E[index]
117
+ delta = Eal[idx] - E[idx]
108
-
109
- # (6-1)E[i,j] > Eal[i,j] の場合、もともとの要素S[i,j]をSal[i,j]に更新し、(7)以降の過程へ。
118
+
110
-
111
- if delta < 0:
119
+ if delta < 0: # Eal[idx] < E[idx]
112
-
120
+
113
- S[index] = Sal[index] # 更新する
121
+ S[idx] = Sal[idx] # 更新する
114
-
115
- # (6-2) E[i,j] < Eal[i,j] の場合、deltaE[i,j]=Eal[i,j]-E[i,j]を計算します。
122
+
116
-
117
- # そして確率exp(-deltaE[i,j]/0.001)でもともとの要素S[i,j]をSal[i,j]に更新します。
118
-
119
- elif delta >= 0 and np.random.random() < np.exp(-delta / 0.001):
123
+ elif delta >= 0 and np.random.random() < np.exp(-delta / 0.001):
124
+
120
-
125
+ # Eal[idx] >= E[idx] かつ np.exp(-delta / 0.001) の確率
126
+
121
- S[index] = Sal[index] # 更新する
127
+ S[idx] = Sal[idx] # 更新する
128
+
129
+
130
+
122
-
131
+ # (9)
123
-
124
-
132
+
125
- # (8)(7)の操作後の行列Sを用いて行列Eを計算する。
133
+ S_sum += S.sum()
126
-
134
+
127
- E = calc_E(S)
135
+ E_sum += E.sum()
136
+
137
+
138
+
128
-
139
+ S_mean = S_sum / iterations
140
+
141
+ E_mean = E_sum / iterations
142
+
129
- return S, E
143
+ return S_mean, E_mean
130
144
 
131
145
 
132
146
 
133
147
  Nx, Ny = 8, 8
134
148
 
135
- S, E = Metropolis(Nx, Ny, iterations=10000)
149
+ S_mean, E_mean = Metropolis(Nx, Ny, iterations=10000)
136
150
 
137
151
  # それぞれの行列の行列式の値Ssum=np.sum(S)、Esum=np.sum(E)をもとめることができれば計算終了です。
138
152
 
139
- print('Size: ({}, {}), S.sum(): {}, E.sum(): {}'.format(Nx, Ny, S.sum(), E.sum()))
153
+ print('Size: ({}, {}), S_mean: {}, E_mean: {}'.format(Nx, Ny, S_mean, E_mean))
140
-
154
+
141
- # Size: (8, 8), S.sum(): 148.1828022790117, E.sum(): -108.97739996283
155
+ # Size: (8, 8), S_mean: 222.99301060088746, E_mean: 2.9429925289412937
142
156
 
143
157
  ```
144
158
 

6

s

2018/10/24 11:42

投稿

tiitoi
tiitoi

スコア21956

test CHANGED
@@ -26,17 +26,23 @@
26
26
 
27
27
 
28
28
 
29
- def update(S, Sal):
29
+ def calc_Eal(S, Sal):
30
-
30
+
31
- E = np.empty_like(S) # Sと同じ形の行列Eを作る
31
+ E = np.empty_like(S)
32
32
 
33
33
 
34
34
 
35
35
  Nx, Ny = S.shape
36
36
 
37
- S_pad = np.pad(S, (1, 1), "wrap") # 行列Sの両端で折り返す周期的境界条件
37
+ S_pad = np.pad(S, (1, 1), "wrap") # 周期的境界条件
38
-
38
+
39
- Sal_pad = np.pad(Sal, (1, 1), "wrap") # 行列Sの両端で折り返す周期的境界条件
39
+ Sal_pad = np.pad(Sal, (1, 1), "wrap") # 周期的境界条件
40
+
41
+
42
+
43
+ # 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))
44
+
45
+ # E を計算するときは、S = Sal
40
46
 
41
47
  for i, j in np.dstack(np.mgrid[1:Nx + 1, 1:Ny + 1]).reshape(-1, 2): # Eの各要素(i,j)を計算する
42
48
 
@@ -52,65 +58,87 @@
52
58
 
53
59
 
54
60
 
61
+ def calc_E(S):
62
+
63
+ return calc_Eal(S, S)
64
+
65
+
66
+
55
67
  def Metropolis(Nx=8, Ny=8, iterations=100):
56
68
 
69
+ # (1)Nx * Ny の行列を生成し、各成分を (i,j) であらわす。
70
+
71
+ # (2)行列の各要素は 0 ~ 2 * np.pi の乱数を取るものとする。この値を S[i,j] とする。
72
+
57
- S = np.random.uniform(0, 2 * np.pi, (Nx, Ny)) # Nx * Ny の行列Sを形成。各要素は0~2πの乱数
73
+ S = np.random.uniform(0, 2 * np.pi, (Nx, Ny))
74
+
58
-
75
+ # (3) E[i,j]を以下のように定義し、すべてのS[i,j]に対して計算を行う。
76
+
59
- E = update(S, S)
77
+ E = calc_E(S)
60
78
 
61
79
 
62
80
 
63
81
  indices = np.dstack(np.mgrid[0:Nx, 0:Ny]).reshape(-1, 2)
64
82
 
83
+ # (7)(4)~(6)を10000回繰り返すことでおおよそすべての要素S[i,j]がスピン更新を経験するようにします。
84
+
65
- for k in range(iterations): # モンテカルロステップ数を 10000 とする
85
+ for k in range(iterations):
66
-
67
-
68
-
86
+
87
+
88
+
69
- # Sal 作成する。
89
+ # (4)行列Sのある要素S(i,j)ランダムに抽出し、その要素だけを0~2πの乱数で変動させる。
70
90
 
71
91
  Sal = S.copy()
72
92
 
73
- # ランダムに要素を選択し、乱数に置き換える。
74
-
75
93
  index = tuple(indices[np.random.randint(len(indices))])
76
94
 
77
95
  Sal[index] = np.random.uniform(0, 2 * np.pi)
78
96
 
79
97
 
80
98
 
81
- # Eal を計算する。
99
+ #(5)変動させたS(i,j)をSal(i,j)として、以下の計算を行い、Eal(i,j)を計算する。
82
-
100
+
83
- Eal = update(S, Sal)
101
+ Eal = calc_Eal(S, Sal)
84
-
85
-
86
-
102
+
103
+
104
+
87
- # 更新 Eal[index] - E[index]
105
+ # (6)ここで(3)で計算したE[i,j]と(5)で計算したEal[i,j]とを比較します。
88
106
 
89
107
  delta = Eal[index] - E[index]
90
108
 
91
- # Eal[index] < E[index]
109
+ # (6-1)E[i,j] > Eal[i,j] の場合、もともとの要素S[i,j]をSal[i,j]に更新し、(7)以降の過程へ。
92
110
 
93
111
  if delta < 0:
94
112
 
95
113
  S[index] = Sal[index] # 更新する
96
114
 
97
- # Eal[index] >= E[index] かつ確率 np.exp(-delta / 0.001) で更新
115
+ # (6-2) E[i,j] < Eal[i,j] の場合、deltaE[i,j]=Eal[i,j]-E[i,j]を計算します。
116
+
117
+ # そして確率exp(-deltaE[i,j]/0.001)でもともとの要素S[i,j]をSal[i,j]に更新します。
98
118
 
99
119
  elif delta >= 0 and np.random.random() < np.exp(-delta / 0.001):
100
120
 
101
121
  S[index] = Sal[index] # 更新する
102
122
 
123
+
124
+
103
-
125
+ # (8)(7)の操作後の行列Sを用いて行列Eを計算する。
126
+
127
+ E = calc_E(S)
104
128
 
105
129
  return S, E
106
130
 
107
131
 
108
132
 
133
+ Nx, Ny = 8, 8
134
+
109
- S, E = Metropolis(Nx=8, Ny=8, iterations=10000)
135
+ S, E = Metropolis(Nx, Ny, iterations=10000)
136
+
110
-
137
+ # それぞれの行列の行列式の値Ssum=np.sum(S)、Esum=np.sum(E)をもとめることができれば計算終了です。
138
+
111
- print('S', S)
139
+ print('Size: ({}, {}), S.sum(): {}, E.sum(): {}'.format(Nx, Ny, S.sum(), E.sum()))
112
-
140
+
113
- print('E', E)
141
+ # Size: (8, 8), S.sum(): 148.1828022790117, E.sum(): -108.97739996283
114
142
 
115
143
  ```
116
144
 

5

a

2018/10/20 04:14

投稿

tiitoi
tiitoi

スコア21956

test CHANGED
@@ -26,23 +26,27 @@
26
26
 
27
27
 
28
28
 
29
- def update(S):
29
+ def update(S, Sal):
30
+
31
+ E = np.empty_like(S) # Sと同じ形の行列Eを作る
32
+
33
+
30
34
 
31
35
  Nx, Ny = S.shape
32
36
 
33
- E = np.empty_like(S) # Sと同じ形行列Eを作る
37
+ S_pad = np.pad(S, (1, 1), "wrap") # 行列Sの両端で折り返す周期的境界条件
34
38
 
35
- S_pad = np.pad(S, (1, 1), "wrap") # 行列Sの両端で折り返す周期的境界条件
39
+ Sal_pad = np.pad(Sal, (1, 1), "wrap") # 行列Sの両端で折り返す周期的境界条件
36
40
 
37
41
  for i, j in np.dstack(np.mgrid[1:Nx + 1, 1:Ny + 1]).reshape(-1, 2): # Eの各要素(i,j)を計算する
38
42
 
39
- E[i - 1, j - 1] = -np.cos(S_pad[i, j] - S_pad[i - 1, j]) \
43
+ E[i - 1, j - 1] = -np.cos(Sal_pad[i, j] - S_pad[i - 1, j]) \
40
44
 
41
- - np.cos(S_pad[i, j] - S_pad[i + 1, j]) \
45
+ - np.cos(Sal_pad[i, j] - S_pad[i + 1, j]) \
42
46
 
43
- - np.cos(S_pad[i, j] - S_pad[i, j + 1]) \
47
+ - np.cos(Sal_pad[i, j] - S_pad[i, j + 1]) \
44
48
 
45
- - np.cos(S_pad[i, j] - S_pad[i, j - 1])
49
+ - np.cos(Sal_pad[i, j] - S_pad[i, j - 1])
46
50
 
47
51
  return E
48
52
 
@@ -52,7 +56,7 @@
52
56
 
53
57
  S = np.random.uniform(0, 2 * np.pi, (Nx, Ny)) # Nx * Ny の行列Sを形成。各要素は0~2πの乱数
54
58
 
55
- E = update(S)
59
+ E = update(S, S)
56
60
 
57
61
 
58
62
 
@@ -60,41 +64,41 @@
60
64
 
61
65
  for k in range(iterations): # モンテカルロステップ数を 10000 とする
62
66
 
67
+
68
+
69
+ # Sal を作成する。
70
+
71
+ Sal = S.copy()
72
+
63
- # ランダムに要素を選択る。
73
+ # ランダムに要素を選択し、乱数に置き換える。
64
74
 
65
75
  index = tuple(indices[np.random.randint(len(indices))])
66
76
 
67
-
77
+ Sal[index] = np.random.uniform(0, 2 * np.pi)
68
78
 
69
- # 乱数を生成する。
70
79
 
71
- E_tmp = E.copy()
72
80
 
73
- E_tmp[index] = np.random.uniform(0, 2 * np.pi)
81
+ # Eal を計算する。
74
82
 
75
-
83
+ Eal = update(S, Sal)
76
84
 
77
- # 計算
78
85
 
79
- E_tmp = update(E_tmp)
80
86
 
81
-
87
+ # 更新 Eal[index] - E[index]
82
88
 
83
- # 更新 E_tmp[index] - E[index]
89
+ delta = Eal[index] - E[index]
84
90
 
85
- delta = E_tmp[index] - E[index]
86
-
87
- # E_tmp[index] < E[index]
91
+ # Eal[index] < E[index]
88
92
 
89
93
  if delta < 0:
90
94
 
91
- E[index] = E_tmp[index] # 更新する
95
+ S[index] = Sal[index] # 更新する
92
96
 
93
- # E_tmp[index] <= E[index] かつ確率 np.exp(-delta / 0.001) で更新
97
+ # Eal[index] >= E[index] かつ確率 np.exp(-delta / 0.001) で更新
94
98
 
95
99
  elif delta >= 0 and np.random.random() < np.exp(-delta / 0.001):
96
100
 
97
- E[index] = E_tmp[index] # 更新する
101
+ S[index] = Sal[index] # 更新する
98
102
 
99
103
 
100
104
 

4

c

2018/10/19 16:40

投稿

tiitoi
tiitoi

スコア21956

test CHANGED
@@ -23,10 +23,6 @@
23
23
  ```
24
24
 
25
25
  import numpy as np
26
-
27
- import random
28
-
29
- np.set_printoptions(linewidth=200)
30
26
 
31
27
 
32
28
 
@@ -84,19 +80,23 @@
84
80
 
85
81
 
86
82
 
87
- # 更新
83
+ # 更新 E_tmp[index] - E[index]
88
84
 
89
85
  delta = E_tmp[index] - E[index]
90
86
 
91
- if delta < 0 or np.random.random() > np.exp(-delta / 0.001):
87
+ # E_tmp[index] < E[index]
92
88
 
89
+ if delta < 0:
90
+
93
- continue # 更新しない
91
+ E[index] = E_tmp[index] # 更新する
92
+
93
+ # E_tmp[index] <= E[index] かつ確率 np.exp(-delta / 0.001) で更新
94
+
95
+ elif delta >= 0 and np.random.random() < np.exp(-delta / 0.001):
96
+
97
+ E[index] = E_tmp[index] # 更新する
94
98
 
95
99
 
96
-
97
- E[index] = E_tmp[index] # 更新する
98
-
99
-
100
100
 
101
101
  return S, E
102
102
 

3

s

2018/10/19 11:37

投稿

tiitoi
tiitoi

スコア21956

test CHANGED
@@ -86,15 +86,15 @@
86
86
 
87
87
  # 更新
88
88
 
89
- if E[index] > E_tmp[index]:
89
+ delta = E_tmp[index] - E[index]
90
90
 
91
- continue
91
+ if delta < 0 or np.random.random() > np.exp(-delta / 0.001):
92
92
 
93
- else:
93
+ continue # 更新しない
94
94
 
95
- delta = E_tmp[index] - E[index]
96
95
 
96
+
97
- E[index] = np.exp(-delta / 0.001)
97
+ E[index] = E_tmp[index] # 更新する
98
98
 
99
99
 
100
100
 

2

d

2018/10/19 11:06

投稿

tiitoi
tiitoi

スコア21956

test CHANGED
@@ -109,3 +109,73 @@
109
109
  print('E', E)
110
110
 
111
111
  ```
112
+
113
+
114
+
115
+ ## 追記
116
+
117
+
118
+
119
+ indices には配列 E の各要素のインデックスが入っています。
120
+
121
+ ```
122
+
123
+ indices = np.array([[0, 0],
124
+
125
+ [0, 1],
126
+
127
+ [0, 2],
128
+
129
+ [0, 3],
130
+
131
+ [0, 4],
132
+
133
+ [0, 5],
134
+
135
+ [0, 6],
136
+
137
+ [0, 7]])
138
+
139
+ # 実際はもっとあるが、一部だけ
140
+
141
+ ```
142
+
143
+
144
+
145
+ 配列の長さを len() で取得すると、8 となります。
146
+
147
+ ```
148
+
149
+ print(len(indices)) # 8
150
+
151
+ ```
152
+
153
+
154
+
155
+ また、indices の4つ目の値を取得するには、以下のようになるかと思います。
156
+
157
+
158
+
159
+ ```
160
+
161
+ print(indices[4]) # [0 4]
162
+
163
+ ```
164
+
165
+
166
+
167
+ なので、[0, 8) の範囲の乱数を生成して、そのインデックスの indices の要素を得ることで、
168
+
169
+ E の要素をランダムに選ぶことになります。
170
+
171
+
172
+
173
+ ```
174
+
175
+ index_to_pick = np.random.randint(len(indices)) # [0, 8) の範囲の乱数を生成
176
+
177
+ print(index_to_pick) # 2
178
+
179
+ print(indices[index_to_pick]) # [0 2]
180
+
181
+ ```

1

d

2018/10/19 09:51

投稿

tiitoi
tiitoi

スコア21956

test CHANGED
@@ -14,12 +14,98 @@
14
14
 
15
15
 
16
16
 
17
+
18
+
17
- S と同形状の乱数を作成る場合は第3引数で size を指定してください
19
+ 上記のアルゴリズムに従うなら、以下のような感かね
20
+
21
+
18
22
 
19
23
  ```
20
24
 
21
- #ランダムに抽出した要素S(i,j)の値をランダムで決定する
25
+ import numpy as np
22
26
 
27
+ import random
28
+
29
+ np.set_printoptions(linewidth=200)
30
+
31
+
32
+
33
+ def update(S):
34
+
35
+ Nx, Ny = S.shape
36
+
37
+ E = np.empty_like(S) # Sと同じ形の行列Eを作る
38
+
39
+ S_pad = np.pad(S, (1, 1), "wrap") # 行列Sの両端で折り返す周期的境界条件
40
+
41
+ for i, j in np.dstack(np.mgrid[1:Nx + 1, 1:Ny + 1]).reshape(-1, 2): # Eの各要素(i,j)を計算する
42
+
43
+ E[i - 1, j - 1] = -np.cos(S_pad[i, j] - S_pad[i - 1, j]) \
44
+
45
+ - np.cos(S_pad[i, j] - S_pad[i + 1, j]) \
46
+
47
+ - np.cos(S_pad[i, j] - S_pad[i, j + 1]) \
48
+
49
+ - np.cos(S_pad[i, j] - S_pad[i, j - 1])
50
+
51
+ return E
52
+
53
+
54
+
55
+ def Metropolis(Nx=8, Ny=8, iterations=100):
56
+
57
+ S = np.random.uniform(0, 2 * np.pi, (Nx, Ny)) # Nx * Ny の行列Sを形成。各要素は0~2πの乱数
58
+
59
+ E = update(S)
60
+
61
+
62
+
63
+ indices = np.dstack(np.mgrid[0:Nx, 0:Ny]).reshape(-1, 2)
64
+
65
+ for k in range(iterations): # モンテカルロステップ数を 10000 とする
66
+
67
+ # ランダムに要素を選択する。
68
+
69
+ index = tuple(indices[np.random.randint(len(indices))])
70
+
71
+
72
+
73
+ # 乱数を生成する。
74
+
75
+ E_tmp = E.copy()
76
+
23
- Sal = np.random.uniform(0,2*np.pi, S.shape)
77
+ E_tmp[index] = np.random.uniform(0, 2 * np.pi)
78
+
79
+
80
+
81
+ # 計算
82
+
83
+ E_tmp = update(E_tmp)
84
+
85
+
86
+
87
+ # 更新
88
+
89
+ if E[index] > E_tmp[index]:
90
+
91
+ continue
92
+
93
+ else:
94
+
95
+ delta = E_tmp[index] - E[index]
96
+
97
+ E[index] = np.exp(-delta / 0.001)
98
+
99
+
100
+
101
+ return S, E
102
+
103
+
104
+
105
+ S, E = Metropolis(Nx=8, Ny=8, iterations=10000)
106
+
107
+ print('S', S)
108
+
109
+ print('E', E)
24
110
 
25
111
  ```