-
Notifications
You must be signed in to change notification settings - Fork 0
/
bzd
151 lines (146 loc) · 6.46 KB
/
bzd
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
ClearAll["Global`*"]
k[l_, t_] :=
SparseArray[# -> 1 & /@ {{1, 2}, {1, l}, {l, 1}, {l, l - 1},
Sequence @@ (Table[{{i, i - 1}, {i, i + 1}}, {i, 2, l - 1}] //
Flatten[#, 1] &)}];
R[xs_, i_] := {1 + (1 -
G[[1]][[i]][[i]]) (Exp[2 *lmd* S[[xs]][[i]]] - 1),
1 + (1 - G[[2]][[i]][[i]]) (Exp[-2 *lmd * S[[xs]][[i]]] - 1)};
Upq[i_, xs_, csx_, dd_, rx_] :=
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[rx == 1,
If[Mod[csx, 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]].Inverse[Bs[[xs]]],
Bd[[xs]].G[[2]].Inverse[Bd[[xs]]]};]]]];
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];
qrd[x_] := Block[{q, d, r}, {q, r} = QRDecomposition[x];
d = DiagonalMatrix[Table[r[[i]][[i]], {i, Length@r}]];
r = Table[
r[[i]][[j]]/d[[i]][[i]], {i, Length@r}, {j,
Length@r}]; {ConjugateTranspose[q], d, r}];
ld[x_] := Block[{jg = dw}, Table[jg = jg.i, {i, x}]; jg];
fxk[x_] := Table[ld[x[[i ;; i + 3]]], {i, m/4}];
(*Svd[x_]:=Block[{u,w,v,vp},Nest[If[#\[Equal]1,{u,w,v}=\
SingularValueDecomposition[x[[-#]]];
v=Transpose@v;#+1,{u,w,vp}=SingularValueDecomposition[x[[-#]].u.w];vp=\
Transpose@vp;
v=vp.v;#+1]&,1,Length@x];{u,w,v}];*)
Svd[x_] :=
Block[{u, w, v, vp},
Nest[If[# == 1, {u, w, v} = qrd[x[[#]]]; # + 1, {u, w, vp} =
qrd[x[[#]].u.w]; v = vp.v; # + 1] &, 1, Length@x]; {u, w, v}];
gd[x_] :=
Block[{u, w, v, zj}, {u, w, v} =
Svd[fxk[{Sequence @@ sdd[x][[2 ;; -1]],
Sequence @@ dds[x][[2 ;; -1]]}]];
zj = qrd[Inverse[u].Inverse[v] + w];
Inverse[zj[[1]].v].Inverse[zj[[2]]].Inverse[u.zj[[3]]]];
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, Bd[[m - #]]}; # + 1] &, 0, m - x + 1]; jg];
gs[x_] :=
Block[{u, w, v, zj}, {u, w, v} =
Svd[fxk[{Sequence @@ sdds[x][[2 ;; -1]],
Sequence @@ ddss[x][[2 ;; -1]]}]];
zj = qrd[Inverse[u].Inverse[v] + w];
Inverse[zj[[1]].v].Inverse[zj[[2]]].Inverse[u.zj[[3]]]];
(*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_] :=
If[Times @@ R[xs, i] > RandomReal[], S[[xs]][[i]] = -S[[xs]][[i]];
rjs = Times @@ R[xs, i]; 1, 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"]; Print[G],
Gs = {gs[yxs], gd[yxs]}; error = Mean@Flatten[Abs[Gs - G]];
G = Gs]];
yr[warm_] :=
Block[{cs = 0, r},
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]];
Gcc[r_] =
If[r == 1, qscl[U, bt, m, t, l, mu];
Gcc = 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}];];
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];
cs = cs + 1;
jg = {Sequence @@ jg, {Nz[1] + Nz[-1],
Nz[1] - Nz[-1]}}; {#[[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]}}; {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]}}; {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, rjs}, qs[1, U, bt, m, t, l, mu, 1]; yr[warm]; cl[o]];