-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathModInv.m
129 lines (108 loc) · 3.1 KB
/
ModInv.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
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
function [ Xmc, F7, F5, F6 ] = ModInv( finv, r, Vp, OBS, PAR, per, M2, TPSDR, Opcion )
switch Opcion
case 1
F5 = Fig5( finv, M2, TPSDR, r);
F6 = Fig6;
Xmc = Vp';
i = 0;
RMS = 1;
Z = Jacobiano( finv, r, Vp, OBS, PAR, per, TPSDR );
disp(' ')
pause(1)
while(RMS>0.05)
i = i+1;
RMS = sqrt(sum((M2 - TPSDR).^2)/length(M2));
figure(3)
hold on
bar(i,RMS)
set(gca, 'XLim', [0.5, i+.5], 'XTick', 1:1:i)
INV_ZtZ = inv(Z'*Z);
TPSDmc = DirectoCCA(finv,r,Xmc)'; %transpuesto solo para visualizacion
Xmc = Xmc + INV_ZtZ * Z' * ( M2 - TPSDmc );
TPSDcal = DirectoCCA(finv,r,Xmc)';
Z = Jacobiano( finv, r, Xmc, OBS, PAR, per, TPSDcal );
figure(2);
hold on
FF2 = loglog(finv,TPSDcal,'--r','LineWidth',1);
legend('M_{Obs}','PSD_{0}',strcat('PSD_{iter:', num2str(i),'}'))
pause(0.5)
% Determinante=det(Z'*Z);
% DVS_ZtZ=svd(Z'*Z);
% Ind=DVS_ZtZ(length(Z))/DVS_ZtZ(1);
% % Ind=DVS_ZtZ(1)/DVS_ZtZ(length(Z));
%
%
% if (Ind < 10^3)
% disp(['Inestabilidad en la inversa de la matriz ZtZ segun el indice DVS=', ' ', num2str(Ind)])
% break
% elseif(Determinante < 0)
% disp(['Inestabilidad en la inversa de la matriz ZtZ. Determinante=', ' ' ,num2str(Determinante)])
% break
% else
%
% end
if RMS <= 0.05
break
else
delete(FF2)
end
legend('M_{Obs}','PSD_{0}')
TPSDR=TPSDcal;
% disp(['Iteracion: ',num2str(i)])
end
disp(['Iteracion: ',num2str(i)])
F7 = Fig7( finv, Xmc);
RMS=RMS*100;
disp(['RMS = ',num2str(RMS),'%'])
VpCal = Xmc
case 2
F5 = Fig5( finv, M2, TPSDR, r );
F6 = Fig6;
Xmc = Vp';
i = 0;
RMS = 1;
I = diag(ones(1,PAR));
% alfa = 0.000001;
disp(' ')
disp('Define alfa para la regularizacion')
disp('alfa muy pequenio --> Minimo cuadrado')
disp('alfa grande --------> Sobre amortiguamiento')
disp('alfa medio ---------> Regularizado')
alfa = input('alfa = ');
disp(' ')
Z = Jacobiano( finv, r, Vp, OBS, PAR, per, TPSDR );
% pause(2)
while(RMS > 0.05)
i = i+1;
RMS = sqrt(sum((M2 - TPSDR).^2)/length(M2));
figure(3)
hold on
bar(i,RMS)
set(gca, 'XLim', [0.5, i+.5], 'XTick', 1:1:i)
TPSDmc = DirectoCCA(finv,r,Xmc)'; %transpuesto solo para visualizacion
Vp_reg = Xmc + inv(Z'*Z + (alfa*I)) * Z' * (M2 - TPSDmc);
TPSDcal = DirectoCCA(finv,r,Vp_reg)'; %transpuesto solo para visualizacion
Z = Jacobiano( finv, r, Vp_reg, OBS, PAR, per, TPSDcal );
figure(2);
hold on
FF2 = loglog(finv,TPSDcal,'--k','LineWidth',1);
legend('M_{Obs}','PSD_{0}',strcat('PSD_{iter:', num2str(i),'}'))
pause(0.5)
% disp(['Iteracion: ',num2str(i)])
% Residual = abs(sum(M2 - TPSDcal));
if RMS<=0.05
break
else
delete(FF2)
end
legend('M_{Obs}','PSD_{0}')
TPSDR = TPSDcal;
Xmc = Vp_reg;
end
disp(['Iteracion: ',num2str(i)])
F7 = Fig7( finv, Vp_reg);
RMS=RMS*100;
disp(['RMS = ',num2str(RMS),'%'])
VpCal = Xmc
end
end