-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathpcs_evolution.m
108 lines (88 loc) · 3.13 KB
/
pcs_evolution.m
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
L = 100;
D_object = 30;
min_maxeig = -0.5;
maxeig = -0.01;
nrepeat = 10;
r = [0:L-1];
dt = 0.05;
Lt = floor(L/dt*0.8);
subdivisions = 10;
H_reservoir = - sparse(mod(r+1,L)+1,r+1,0.5) - sparse(mod(r-1,L)+1,r+1,0.5);
H_object = rand(D_object,D_object) + i*rand(D_object,D_object);
H_object = 0.5*(H_object + H_object');
%H_object = [1 0; 0 -1];
list = sort(eig(H_object));
H_object_new = H_object-(list(1)+list(2))/2*eye(D_object);
scale = max(eig(H_object_new));
H = sparse(D_object*L,D_object*L);
for ii = 0:D_object-1
H(ii*L+[0:L-1]+1,ii*L+[0:L-1]+1) = H_reservoir;
end
vobject_init = zeros(D_object,1);
v_init = zeros(D_object*L,1);
%Algorithm
for ii = 0:D_object-1
vobject_init(ii+1) = (rand-0.5) + i*(rand-0.5);
end
H_object_new = H_object_new/scale*(min_maxeig);
for jj = 0:D_object-1
for ii = 0:D_object-1
H(ii*L+1,jj*L+1) = H(ii*L+1,jj*L+1) + H_object_new(ii+1,jj+1)/2.0;
H(ii*L+2,jj*L+1) = H(ii*L+2,jj*L+1) + H_object_new(ii+1,jj+1)/2.0;
H(ii*L+1,jj*L+2) = H(ii*L+1,jj*L+2) + H_object_new(ii+1,jj+1)/2.0;
H(ii*L+2,jj*L+2) = H(ii*L+2,jj*L+2) + H_object_new(ii+1,jj+1)/2.0;
end
end
[vv,dd] = eig(H_object);
[ee,ord] = sort(diag(dd));
index = find(abs(vv(:,ord(1))) == max(abs(vv(:,ord(1)))));
vv_exact = vv(:,ord(1))/vv(index,ord(1));
for ntrial = 1:nrepeat
for ii = 0:D_object-1
v_init(ii*L+1) = vobject_init(ii+1,1);
v_init(ii*L+2) = vobject_init(ii+1,1);
end
psi1(:,1) = exponentiate_pcs(v_init,H,-i*dt,subdivisions);
for nt = 1:Lt
maxeig_run = -2*(maxeig+min_maxeig)/Lt*abs(nt-Lt/2) + maxeig;
H_object_new = H_object_new/scale*maxeig_run;
for jj = 0:D_object-1
for ii = 0:D_object-1
H(ii*L+1,jj*L+1) = H(ii*L+1,jj*L+1) + H_object_new(ii+1,jj+1)/2.0;
H(ii*L+2,jj*L+1) = H(ii*L+2,jj*L+1) + H_object_new(ii+1,jj+1)/2.0;
H(ii*L+1,jj*L+2) = H(ii*L+1,jj*L+2) + H_object_new(ii+1,jj+1)/2.0;
H(ii*L+2,jj*L+2) = H(ii*L+2,jj*L+2) + H_object_new(ii+1,jj+1)/2.0;
end
end
[vv,dd] = eig(H_object);
[ee,ord] = sort(diag(dd));
index = find(abs(vv(:,ord(1))) == max(abs(vv(:,ord(1)))));
vv_exact = vv(:,ord(1))/vv(index,ord(1));
psi1(:,nt+1) = exponentiate_pcs(psi1(:,nt),H,-i*dt,subdivisions);
end
figure(ntrial)
hold on
for ii = 0:D_object-1
plot(dt*[1:Lt],log(abs(psi1(ii*L+1,1:Lt)+psi1(ii*L+2,1:Lt)).^2))
end
hold off
for ii = 0:D_object-1
v_PC(ii+1,1) = psi1(ii*L+1,Lt) + psi1(ii*L+2,Lt);
end
index = find(abs(v_PC) == max(abs(v_PC)));
ntrial
v_PC = v_PC/v_PC(index);
[v_PC vv_exact]
initial_overlap(ntrial,1) = (abs(vv_exact'*vobject_init))^2/((vv_exact'*vv_exact)*(vobject_init'*vobject_init));
final_overlap(ntrial,1) = (abs(v_PC'*vv_exact))^2/((vv_exact'*vv_exact)*(v_PC'*v_PC));
initial_overlap(ntrial,1)
final_overlap(ntrial,1)
vobject_init = v_PC;
end
"initial"
initial_overlap(1)
"overlaps"
final_overlap
figure(nrepeat+1)
plot(log(1-final_overlap))
eig(H_object_new)