-
Notifications
You must be signed in to change notification settings - Fork 0
/
xxxxxx
206 lines (195 loc) · 8.09 KB
/
xxxxxx
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
ClearAll["Global`*"]
k[l_, t_] :=
SparseArray[# -> -t & /@ {{1, 2}, {1, l}, {l, 1}, {l, l - 1},
Sequence @@ (Table[{{i, i - 1}, {i, i + 1}}, {i, 2, l - 1}] //
Flatten[#, 1] &)}];
ld[x_] := Block[{jg = dw}, Table[jg = jg.i, {i, x}]; jg];
fxk[x_] := Table[ld[x[[i ;; i + 3]]], {i, m/2}];
Svd[x_] :=
Block[{u, w, v, vp},
Nest[If[# == 1, {u, w, v} = SingularValueDecomposition[x[[-#]]];
v = Conjugate[Transpose@v]; # + 1, {u, w, vp} =
SingularValueDecomposition[x[[-#]].u.w];
vp = Conjugate[Transpose@vp]; v = vp.v; # + 1] &, 1,
Length@x]; {u, w, v}];
sdd[x_] :=
Block[{jg = {dw}},
Nest[If[# == 0, # + 1, jg = {Sequence @@ jg, Bd[[#]]}; # + 1] &, 0,
x]; jg];
dds[x_] :=
Block[{jg = {dw}},
Nest[If[# == 0, jg = {Sequence @@ jg, Bd[[m - #]]}; # + 1,
jg = {Sequence @@ jg, Bd[[m - #]]}; # + 1] &, 0, m - x + 1]; jg];
gd[x_] :=
Block[{u, w, v, zj}, {u, w, v} =
Svd[{Sequence @@ sdd[x][[2 ;; -1]],
Sequence @@ dds[x][[2 ;; -1]]}];
zj = SingularValueDecomposition[
Transpose[u].Transpose[Conjugate[v]] +
w]; (Transpose[Conjugate[v]].zj[[3]]).Table[
If[zj[[2]][[i]][[j]] != 0 && i == j,
zj[[2]][[i]][[j]] = 1/zj[[2]][[i]][[j]], zj[[2]][[i]][[j]]], {i,
Length@zj[[2]]}, {j,
Length@zj[[2]]}].(Transpose[zj[[1]]].Transpose[u])];
gs[x_] :=
Block[{u, w, v, zj}, {u, w, v} =
Svd[{Sequence @@ sdds[x][[2 ;; -1]],
Sequence @@ ddss[x][[2 ;; -1]]}];
zj = SingularValueDecomposition[
Transpose[u].Transpose[Conjugate[v]] +
w]; (Transpose[Conjugate[v]].zj[[3]]).Table[
If[zj[[2]][[i]][[j]] != 0 && i == j,
zj[[2]][[i]][[j]] = 1/zj[[2]][[i]][[j]], zj[[2]][[i]][[j]]], {i,
Length@zj[[2]]}, {j,
Length@zj[[2]]}].(Transpose[zj[[1]]].Transpose[u])];
sdds[x_] :=
Block[{jg = {dw}},
Nest[If[# == 0, # + 1, jg = {Sequence @@ jg, Bs[[#]]}; # + 1] &, 0,
x]; jg];
ddss[x_] :=
Block[{jg = {dw}},
Nest[If[# == 0, jg = {Sequence @@ jg, Bs[[m - #]]}; # + 1,
jg = {Sequence @@ jg, Bs[[m - #]]}; # + 1] &, 0, m - x + 1];
jg];
R[xs_, i_] := {1 + (1 -
G[[1]][[i]][[i]]) (Exp[-2*lmd*-1*S[[xs]][[i]]] - 1),
1 + (1 - G[[2]][[i]][[i]]) (Exp[-2 *lmd*S[[xs]][[i]]] - 1)};
Upq[i_, xs_, cs_, dd_, r_] :=
Block[{gmd = Exp[-2 *lmd*1*S[[xs]][[i]]] - 1,
gms = Exp[2 *lmd* 1* S[[xs]][[i]]] - 1,
dww = Table[If[j == h == i, 1, 0], {h, l}, {j, l}]},
If[r == 1,
If[Mod[cs, gxjg] == 0, qs[xs, U, bt, m, t, l, mu, 0],
If[dd == 1,
G[[1]] =
Table[G[[1]][[h]][[
n]] - ((dww[[h]][[n]] - G[[1]][[h]][[i]]) gms G[[1]][[i]][[
n]])/(1 + (1 + G[[1]][[i]][[i]])*gms), {h, l}, {n, l}];
G[[2]] =
Table[G[[2]][[h]][[
n]] - ((dww[[h]][[n]] - G[[2]][[h]][[i]]) gms G[[2]][[i]][[
n]])/(1 + (1 + G[[1]][[i]][[i]])*gms), {h, l}, {n, l}],
G = {Bs[[xs]].G[[1]].PseudoInverse[Bs[[xs]]],
Bd[[xs]].G[[2]].PseudoInverse[Bd[[xs]]]}]]]];
(*Nz[sx_]:=If[sx\[Equal]1,Mean@Flatten[Table[1-gs[i][[j]][[j]],{i,m},{\
j,l}],1],Mean@Flatten[Table[1-gd[i][[j]][[j]],{i,m},{j,l}],1]];*)
accpet[i_, xs_] :=
Block[{rrr = Times @@ R[xs, i]},
If[rrr > RandomReal[],
If[rrr > 1 || rrr < 0, rjs = {Sequence @@ rjs, rrr}];
S[[xs]][[i]] = -S[[xs]][[i]]; 1, 0]];
(*qs[xs,U,bt,m,t,l,mu,0]*)
qs[yxs_, U0_, bt0_, m0_, t0_, l0_, mu0_, DT0_] :=
Block[{U = U0, bt = bt0, l = l0, m = m0, mu = mu0, t = t0, Gs},
lmd = Log[ (Sqrt[Exp[Abs[U] dt] - 1] + Exp[Abs[U] dt .5])];
vs = Table[lmd/dt S[[xs]][[i]] - mu + U 0.5, {xs, m}, {i, l}];
vd = Table[(-1 lmd)/dt S[[xs]][[i]] - mu + U 0.5, {xs, m}, {i,
l}];
Kma = Table[If[EvenQ@i, k[l, 1][[i]], Table[0, l]], {i, l}] //
Normal;
Kmb = Table[If[OddQ@i, k[l, 1][[i]], Table[0, l]], {i, l}] //
Normal; Bd =
Table[(Kma*Sinh[dt t] + dw Cosh[dt t]).(Kmb*Sinh[dt t] +
dw Cosh[dt t]).(DiagonalMatrix@
Table[Exp[-dt vd[[xs]][[i]]], {i, l}]), {xs, m}];
Bs = Table[(Kma*Sinh[dt t] + dw Cosh[dt t]).(Kmb*Sinh[dt t] +
dw Cosh[dt t]).(DiagonalMatrix@
Table[Exp[-dt vs[[xs]][[i]]], {i, l}]), {xs, m}];
If[DT0 == 1, G = {gs[yxs], gd[yxs]}; Print["start"];
error = {error, "start"}, Gs = {gs[yxs], gd[yxs]};
error = {Sequence @@ error, Mean@Flatten[Abs[Gs - G], 1]};
G = Gs;]];
yr[warm_] :=
Block[{cs = 0, r}, Print["warm"];
Nest[If[#[[1]] != l, r = accpet[Sequence @@ #];
Upq[#[[1]], #[[2]], cs, 1, r];
cs = cs + 1; {#[[1]] + 1, #[[2]]},
If[#[[2]] != m, r = accpet[Sequence @@ #];
Upq[#[[1]], #[[2]], cs, 2, r]; cs = cs + 1;
{1, #[[2]] + 1}, r = accpet[Sequence @@ #];
Upq[#[[1]], #[[2]], cs, 1, r]; cs = cs + 1;
G = {gs[1], gd[1]}; {1, 1}]] &, {1, 1}, warm];
Print["warm over"]];
Gcc[r_] = If[r == 1, qscl[U, bt, m, t, l, mu];
Gc = Table[{gs[i], gd[i]}, {i, m}]];
qscl[U0_, bt0_, m0_, t0_, l0_, mu0_] :=
Block[{U = U0, bt = bt0, l = l0, m = m0, mu = mu0, t = t0, Gs},
lmd = Log[(Sqrt[Exp[Abs[U] dt] - 1] + Exp[Abs[U] dt .5])];
vs = Table[lmd/dt*S[[xs]][[i]] - mu + U 0.5, {xs, m}, {i, l}];
vd = Table[(-1 lmd)/dt*S[[xs]][[i]] - mu + U 0.5, {xs, m}, {i,
l}];
Kma = Table[If[EvenQ@i, k[l, 1][[i]], Table[0, l]], {i, l}] //
Normal;
Kmb = Table[If[OddQ@i, k[l, 1][[i]], Table[0, l]], {i, l}] //
Normal;
Bd = Table[(Kma*Sinh[dt t] + dw Cosh[dt t]).(Kmb*Sinh[dt t] +
dw Cosh[dt t]).(DiagonalMatrix@
Table[Exp[-dt vd[[xs]][[i]]], {i, l}]), {xs, m}];
Bs = Table[(Kma*Sinh[dt t] + dw Cosh[dt t]).(Kmb*Sinh[dt t] +
dw Cosh[dt t]).(DiagonalMatrix@
Table[Exp[-dt vs[[xs]][[i]]], {i, l}]), {xs, m}];];
mz[x_] :=
Which[x == 1,
Mean@Flatten[
Table[Abs[(1 - Gc[[i]][[1]][[j]][[j]])^2 + (1 -
Gc[[i]][[2]][[j]][[j]])^2 -
2*(1 - Gc[[i]][[2]][[j]][[j]])*(1 -
Gc[[i]][[1]][[j]][[j]]) + (1 - Gc[[i]][[2]][[j]][[j]])*
Gc[[i]][[1]][[j]][[j]] + (1 - Gc[[i]][[1]][[j]][[j]])*
Gc[[i]][[2]][[j]][[j]]], {i, m}, {j, l}], 1], x == 2
, Mean@
Flatten[Table[
Abs[(1 - Gc[[i]][[2]][[j]][[j]])*
Gc[[i]][[1]][[j]][[j]] + (1 - Gc[[i]][[2]][[j]][[j]])*
Gc[[i]][[1]][[j]][[j]]], {i, m}, {j, l}], 1]];
Nz[sx_] :=
If[sx == 1,
Mean@Flatten[Table[1 - Gc[[i]][[1]][[j]][[j]], {i, m}, {j, l}], 1],
Mean@Flatten[Table[1 - Gc[[i]][[2]][[j]][[j]], {i, m}, {j, l}],
1]];
cl[o_] :=
Block[{jg = {{ks, ks}}, cs = 0, r, Gc}, Gcc[1];
Gcs = {Sequence @@ Gcs, Gc}; Print["cl"];
G = {Gc[[1]][[1]], Gc[[1]][[2]]};
Nest[If[#[[1]] != l, r = accpet[Sequence @@ #]; Gcc[r];
G = {Gc[[#[[2]]]][[1]], Gc[[#[[2]]]][[2]]};
cs = cs + 1;
jg = {Sequence @@ jg, {Nz[1] + Nz[-1], Nz[1] - Nz[-1], mz[1],
mz[2]}}; {#[[1]] + 1, #[[2]]},
If[#[[2]] != m, r = accpet[Sequence @@ #];
Gcc[r]; G = {Gc[[#[[2]]]][[1]], Gc[[#[[2]]]][[2]]}; cs = cs + 1;
jg = {Sequence @@ jg, {Nz[1] + Nz[-1], Nz[1] - Nz[-1], mz[1],
mz[2]}}; {1, #[[2]] + 1}, r = accpet[Sequence @@ #]; Gcc[r];
cs = cs + 1;
G = {Gc[[1]][[1]], Gc[[1]][[2]]};
jg = {Sequence @@ jg, {Nz[1] + Nz[-1], Nz[1] - Nz[-1], mz[1],
mz[2]}}; {1, 1}]] &, {1, 1}, o]; jg];
Qmc[U0_, bt0_, m0_, t0_, l0_, mu0_, warm_, o_] :=
Block[{U = U0, bt = bt0, l = l0, m = m0, mu = mu0, t = t0,
dt = bt0/m0, dw = IdentityMatrix[l0],
S = Table[RandomChoice[{1, -1}], m0, l0], Bd, Bs, Kmb, Kma, vs,
vd, G}, qs[1, U, bt, m, t, l, mu, 1]; yr[warm]; cl[o]];
gxjg = 3;
U = 4;
Gcs = {{0, 0}};
rjs = {"ks"};
error = {0};
bt = {10, 8, 3, 1, 1/4, 4, 3, 2, 1/2};
l = 6;
m = 8*bt;
mu = U/2;
t = 1;
warm = 5000;
o = 5000;
obs = Table[
Qmc[U, bt[[p]], m[[p]], t, l, mu, warm, o], {p, Length@bt}] //
AbsoluteTiming
jgs = Table[obs[[2]][[i]][[2 ;; -1]], {i, Length@bt}]
nn = Table[{1/bt[[i]], jgs[[i]][[All, 1]] // Mean}, {i, Length@bt}];
mmm = Table[{1/bt[[i]], jgs[[i]][[All, 2]] // Mean}, {i, Length@bt}];
mmmm = Table[{1/bt[[i]], Abs[jgs[[i]][[All, 2]]] // Mean}, {i,
Length@bt}];
mzz = Table[{1/bt[[i]], Abs[jgs[[i]][[All, 3]]] // Mean}, {i,
Length@bt}];
mxx = Table[{1/bt[[i]], Abs[jgs[[i]][[All, 4]]] // Mean}, {i,
Length@bt}];