forked from BarakeelFanseu/2D-GEM
-
Notifications
You must be signed in to change notification settings - Fork 0
/
zoomOut_refine.m
49 lines (42 loc) · 1.34 KB
/
zoomOut_refine.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
function [T12, C21, all_T12, all_C21] = zoomOut_refine(M, N, options)
try
[evecs, evals] = eigs(M.Phi.W, M.Phi.A, options.spec_dim, 1e-6);
catch
% In case of trouble make the laplacian definite
[evecs, evals] = eigs(M.Phi.W - 1e-8*speye(M.n), M.Phi.A, options.spec_dim, 'sm');
end
evals = diag(evals);
if ~isreal(evecs)
evecs(1:2:end) = real(evecs(1:2:end));
evecs(2:2:end) = imag(evecs(1:2:end));
end
[evals, order] = sort(abs(evals),'ascend');
evecs = evecs(:,order);
B1_all = evecs;
try
[evecs, evals] = eigs(N.Phi.W, N.Phi.A, options.spec_dim, 1e-6);
catch
% In case of trouble make the laplacian definite
[evecs, evals] = eigs(N.Phi.W - 1e-8*speye(N.n), N.Phi.A, options.spec_dim, 'sm');
end
evals = diag(evals);
if ~isreal(evecs)
evecs(1:2:end) = real(evecs(1:2:end));
evecs(2:2:end) = imag(evecs(1:2:end));
end
[evals, order] = sort(abs(evals),'ascend');
evecs = evecs(:,order);
B2_all = evecs;
% T12 = N.shots*M.shots';
% T12 = greedy_match(T12);
% [T12, ~] = find(T12);
T12 = knnsearch(N.shots, M.shots,'NSMethod','kdtree');
if nargout > 2, all_T12 = {}; all_C21 = {}; end
for k = options.k_init : options.k_step : options.k_final
B1 = B1_all(:, 1:k);
B2 = B2_all(:, 1:k);
C21 = B1\B2(T12,:);
T12 = knnsearch(B2*C21', B1);
if nargout > 2, all_T12{end+1} = T12; all_C21{end+1} = C21;
end
end