root/applications/dual/IterativeLocal/ildp4cor2.m @ 1266

Revision 720, 14.5 kB (checked in by smidl, 16 years ago)

init

Line 
1function ildp4cor2
2    tic
3   
4    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5    %pocatecni konstanty
6    Iterace = 2; %iterace
7    K = 20; %casy
8    N = 250; %vzorky
9   
10    sigma = 1;
11    Sigmas = [[sigma^2 0 0]; [0 0 0]; [0 0 0]];
12    h = 0;
13    hdx = [0; 0; 0];
14    hdxdx = [[0 0 0]; [0 0 0]; [0 0 0]];
15    Rk = ones(1,K);
16    x0 = [0; 0; 1];
17
18    Wv = zeros(19,K);
19%     Wpi = zeros(19,K);
20    Kpi = ones(4,K);
21        Kpi(4,:) = zeros(1,K);
22%     Wpit = zeros(19,K);
23    Kpit = ones(4,K);
24%     Wpin = zeros(19,K);
25    Kpin = ones(4,K);
26
27    Epsilon = zeros(3,N);
28        Epsilon(1,:) = 1*randn(1,N);
29        Epsilon(2,:) = 1*randn(1,N);
30        Epsilon(3,:) = abs(1*randn(1,N));
31   
32    mFI = [ones(1,N);...
33           Epsilon(1,:);...
34           Epsilon(2,:);...
35           Epsilon(3,:);...
36           Epsilon(1,:).^2;...
37           Epsilon(1,:).*Epsilon(2,:);...
38           Epsilon(1,:).*Epsilon(3,:);...
39           Epsilon(2,:).^2;...
40           Epsilon(2,:).*Epsilon(3,:);...
41           Epsilon(3,:).^2;...
42
43           Epsilon(1,:).^3;...
44           Epsilon(1,:).^2.*Epsilon(2,:);...
45           Epsilon(1,:).^2.*Epsilon(3,:);...
46           Epsilon(2,:).^2.*Epsilon(1,:);...
47           Epsilon(2,:).^3;...
48           Epsilon(2,:).^2.*Epsilon(3,:);...
49           Epsilon(3,:).^2.*Epsilon(1,:);...
50           Epsilon(3,:).^2.*Epsilon(2,:);...
51           Epsilon(3,:).^3];
52    FiFiTInvFi = (mFI*mFI')\mFI;
53%     max(max(abs(FiFiTInvFi)))
54
55    Xall = zeros(3,K,N);
56
57    %globalni
58    kappa = 0;
59    nu = 0;
60   
61    for i = 1:Iterace,
62       
63        for n = 1:N,
64           Xall(:,1,n) = x0 + [sigma*randn(); 0; 0];
65           for k = 1:K-1,
66              Upi = nPi(k,Xall(:,:,n));
67              Xall(1,k+1,n) = Xall(1,k,n) + Xall(2,k,n)*Upi + sigma*randn();
68              Xall(2,k+1,n) = Xall(2,k,n) + (Xall(1,k+1,n) - Xall(1,k,n) - Xall(2,k,n)*Upi)*Upi*Xall(3,k,n)/(Upi^2*Xall(3,k,n) + sigma^2);
69              Xall(3,k+1,n) = Xall(3,k,n) - (Upi^2*Xall(3,k,n)^2)/(Upi^2*Xall(3,k,n) + sigma^2);
70           end
71        end
72       
73%         Xstripe = zeros(3,K);
74        Xstripe = mean(Xall, 3);
75
76        if(i == 1)
77           Xstripe(1,:) = zeros(1,K);
78        end       
79       
80%         subplot(2,1,1)
81%         plot(1:K,Xall(1,:,1),1:K,Xall(1,:,2),1:K,Xall(1,:,3),1:K,Xall(1,:,4),1:K,Xall(1,:,5),1:K,Xall(1,:,6),...
82%             1:K,Xall(1,:,7),1:K,Xall(1,:,8),1:K,Xall(1,:,9),1:K,Xall(1,:,10),1:K,Xall(1,:,11),1:K,Xall(1,:,12),...
83%             1:K,Xall(1,:,13),1:K,Xall(1,:,14),1:K,Xall(1,:,15),1:K,Xall(1,:,16),1:K,Xall(1,:,17),1:K,Xall(1,:,18),...
84%             1:K,Xall(1,:,19),1:K,Xall(1,:,20),1:K,Xall(1,:,21),1:K,Xall(1,:,22),1:K,Xall(1,:,23),1:K,Xall(1,:,24),...
85%             1:K,Xall(1,:,25))
86%         subplot(2,1,2)
87%         plot(1:K,Xstripe(1,:))
88
89        for n = 1:N,
90            Xall(:,K,n) = Xstripe(:,K) + Epsilon(:,n);
91        end
92       
93        for k = K-1:-1:1,
94            [i k]
95%             1]
96            for n = 1:N,
97                Xall(:,k,n) = Xstripe(:,k) + Epsilon(:,n);
98            end
99%             2]
100            for n = 1:N,
101                kappa = k;
102                nu = n;
103                [Uop(n), Hmin(n)] = fminunc(@Hamilt, nPi(k,Xall(:,:,n)), optimset('GradObj','on'));
104            end
105%             3]
106            for n = 1:N,
107                Vn(n) = Hmin(n) + Vtilde(k+1,Xall(:,k,n)-Xstripe(:,k+1));
108            end
109%             4]
110            Wv(:,k) = FiFiTInvFi*Vn';
111%             Wpit(:,k) = FiFiTInvFi*Uop';
112            yt = zeros(1,N);
113            bt = zeros(1,N);
114            pt = zeros(1,N);
115            for n = 1:N,
116%                 ytm1(n) = 0;
117%                 if(k > 1)
118%                     ytm1(n) = Xall(1,k-1,n);
119%                 end
120                yt(n) = Xall(1,k,n);
121                bt(n) = Xall(2,k,n);
122                pt(n) = Xall(3,k,n);
123            end         
124
125            mPsi = [yt',...
126                    bt'.*Uop',...
127                    pt'.*Uop',...
128                    Uop'];
129            PsiPsiTInvPsi = (mPsi'*mPsi)\mPsi';
130           
131            Kpi(:,k) = PsiPsiTInvPsi * (Rk(k)*ones(N,1));
132            Kpi(:,k)
133        end
134       
135%         linesearch
136            X = zeros(3, K);
137            Ual = zeros(1,K);
138            X(:,1) = x0 + [sigma*randn(); 0; 0]; 
139            for k = 1:K-1,                   
140                    Upi = nPiex(k, X, Kpi);
141                    Ual(k) = Upi;
142                    Ktmp = Upi*X(3,k)/(Upi^2*X(3,k) + sigma^2);
143                    X(1,k+1) = X(1,k)+X(2,k)*Upi + sigma*randn(1);
144                    X(2,k+1) = X(2,k) + Ktmp*(X(1,k+1) - X(1,k) - X(2,k)*Upi);
145                    X(3,k+1) = X(3,k)*(1-Ktmp*Upi);
146            end
147            oldloss = sum( (X(1,:) - Rk).^2 + X(3,:).*Ual.^2 );
148
149            alfa = 0;
150            loss = Inf;
151        while((loss > oldloss)||(alfa<1))
152%             Wpin = (1-alfa)*Wpit + alfa*Wpi;
153             Kpin = (1-alfa)*Kpit + alfa*Kpi;             
154                       
155            X = zeros(3, K);
156            Ual = zeros(1,K);
157            X(:,1) = x0 + [sigma*randn(); 0; 0]; 
158            for k = 1:K-1,                   
159                    Upi = nPiex(k, X, Kpin);                   
160                    Ual(k) = Upi;
161                    Ktmp = Upi*X(3,k)/(Upi^2*X(3,k) + sigma^2);
162                    X(1,k+1) = X(1,k)+X(2,k)*Upi + sigma*randn(1);
163                    X(2,k+1) = X(2,k) + Ktmp*(X(1,k+1) - X(1,k) - X(2,k)*Upi);
164                    X(3,k+1) = X(3,k)*(1-Ktmp*Upi);
165            end
166            loss = sum( (X(1,:) - Rk).^2 + X(3,:).*Ual.^2 );
167           
168            alfa = alfa + 0.1;
169            if(alfa > 1)
170               disp('zadne zlepseni')
171%                loss = -1;
172            end
173        end
174%         if(loss > -1)
175%            Wpi = Wpin;
176           Kpi = Kpin;
177%         end
178       
179    end
180   
181    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
182    Kpi
183    %graficky vystup
184            X = zeros(3, K);
185            UU = zeros(1,K);
186           
187                X(:,1) = x0 + [sigma*randn(1); 0; 0]; 
188                for k = 1:K-1,                   
189                    Upi = nPi(k, X);
190                    UU(k) = Upi;
191                    Ktmp = Upi*X(3,k)/(Upi^2*X(3,k) + sigma^2);
192                    X(1,k+1) = X(1,k)+X(2,k)*Upi + sigma*randn(1);
193                    X(2,k+1) = X(2,k) + Ktmp*(X(1,k+1) - X(1,k) - X(2,k)*Upi);
194                    X(3,k+1) = X(3,k)*(1-Ktmp*Upi);
195                end           
196           
197                subplot(5,1,1);
198                plot(1:K,X(1,:))
199                subplot(5,1,2);
200                plot(1:K,X(2,:))
201                subplot(5,1,3);
202                plot(1:K,X(3,:))
203                subplot(5,1,4);
204                plot(1:K,UU)
205     
206    %napoctena chyba
207    yy = zeros(1,K);
208    losses = zeros(1,100);
209    for n=1:100,
210        bb = randn();
211        yy(1) = x0(1) + sigma*randn();
212        for k = 1:K-1,                   
213                    Upi = nPi(k,X);               
214                    yy(k+1) = yy(k)+bb*Upi + sigma*randn(1);                   
215        end
216        losses(n) = sum( (yy - Rk).^2 );
217    end
218   
219    [min(losses) median(losses) max(losses)]
220    subplot(5,1,5);
221    hist(log(losses))
222    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
223   
224    function [valpi] = nPi(kf,Xknf)%nPi(kf,Epsil)%nPi(kf,Xknf)
225        if(kf > 1)
226            yti1 = Xknf(1,kf-1);
227        else
228            yti1 = 0;
229        end
230        valpi = (Rk(kf)-Kpi(1,kf)*Xknf(1,kf))/(Kpi(2,kf)*Xknf(2,kf) + Kpi(3,kf)*Xknf(3,kf) + Kpi(4,kf));
231%         valpi = [1;...
232%                        Epsil(1);...
233%                        Epsil(2);...
234%                        Epsil(3);...
235%                        Epsil(1)^2;...
236%                        Epsil(1)*Epsil(2);...
237%                        Epsil(1)*Epsil(3);...
238%                        Epsil(2)^2;...
239%                        Epsil(2)*Epsil(3);...
240%                        Epsil(3)^2;...
241%                        Epsil(1)^3;...
242%                        Epsil(1)^2*Epsil(2);...
243%                        Epsil(1)^2*Epsil(3);...
244%                        Epsil(2)^2*Epsil(1);...
245%                        Epsil(2)^3;...
246%                        Epsil(2)^2*Epsil(3);...
247%                        Epsil(3)^2*Epsil(1);...
248%                        Epsil(3)^2*Epsil(2);...
249%                        Epsil(3)^3]'*Wpi(:,kf);
250    end
251
252    function [valpi] = nPiex(kf,Xknf,Kpix)%nPiex(kf,Epsil,Ww)%nPi(kf,Xknf,Kpix)
253        if(kf > 1)
254            yti1 = Xknf(1,kf-1);
255        else
256            yti1 = 0;
257        end
258        valpi = (Rk(kf)-Kpi(1,kf)*Xknf(1,kf))/(Kpix(2,kf)*Xknf(2,kf) + Kpix(3,kf)*Xknf(3,kf) + Kpix(4,kf));
259%         valpi = [1;...
260%                        Epsil(1);...
261%                        Epsil(2);...
262%                        Epsil(3);...
263%                        Epsil(1)^2;...
264%                        Epsil(1)*Epsil(2);...
265%                        Epsil(1)*Epsil(3);...
266%                        Epsil(2)^2;...
267%                        Epsil(2)*Epsil(3);...
268%                        Epsil(3)^2;...
269%                        Epsil(1)^3;...
270%                        Epsil(1)^2*Epsil(2);...
271%                        Epsil(1)^2*Epsil(3);...
272%                        Epsil(2)^2*Epsil(1);...
273%                        Epsil(2)^3;...
274%                        Epsil(2)^2*Epsil(3);...
275%                        Epsil(3)^2*Epsil(1);...
276%                        Epsil(3)^2*Epsil(2);...
277%                        Epsil(3)^3]'*Ww(:,kf);
278    end
279
280    function [valHam, valGrad] = Hamilt(Uin)
281        valHam = (Xall(1,kappa+1,nu)-Rk(kappa))^2 + Xall(3,kappa,nu)*Uin^2 ...
282               + [Xall(2,kappa,nu)*Uin;...
283                  (Xall(1,kappa+1,nu) - Xall(1,kappa,nu) - Xall(2,kappa,nu)*Uin)*Uin*Xall(3,kappa,nu)/(Uin^2*Xall(3,kappa,nu) + sigma^2);...
284                  -(Uin^2*Xall(3,kappa,nu)^2)/(Uin^2*Xall(3,kappa,nu) + sigma^2)]'*Vtildedx(kappa+1,Epsilon(:,nu)+Xstripe(:,kappa)-Xstripe(:,kappa+1)) + 1/2*trace(Sigmas*Vtildedxdx(kappa+1,Epsilon(:,nu)+Xstripe(:,kappa)-Xstripe(:,kappa+1)));
285        valGrad = 2*Xall(3,kappa,nu)*Uin ...
286                + [Xall(2,kappa,nu);...
287                   (2*Uin^2*Xall(3,kappa,nu)^2*(Xall(1,kappa,nu) - Xall(1,kappa+1,nu) + Uin*Xall(2,kappa,nu)))/(sigma^2 + Uin^2*Xall(3,kappa,nu))^2 - (Uin*Xall(2,kappa,nu)*Xall(3,kappa,nu))/(sigma^2 + Uin^2*Xall(3,kappa,nu)) - (Xall(3,kappa,nu)*(Xall(1,kappa,nu) - Xall(1,kappa+1,nu) + Uin*Xall(2,kappa,nu)))/(sigma^2 + Uin^2*Xall(3,kappa,nu));...
288                   (2*Uin^3*Xall(3,kappa,nu)^2)/(sigma^2 + Uin^2*Xall(3,kappa,nu))^2 - (2*Uin*Xall(3,kappa,nu))/(sigma^2 + Uin^2*Xall(3,kappa,nu))]'*Vtildedx(kappa+1,Epsilon(:,nu)+Xstripe(:,kappa)-Xstripe(:,kappa+1));
289               
290    end
291
292    function [valVtl] = Vtilde(kin,Epsil)
293        if(kin == K)
294            valVtl = h;
295        else           
296             valVtl = [1;...
297                       Epsil(1);...
298                       Epsil(2);...
299                       Epsil(3);...
300                       Epsil(1)^2;...
301                       Epsil(1)*Epsil(2);...
302                       Epsil(1)*Epsil(3);...
303                       Epsil(2)^2;...
304                       Epsil(2)*Epsil(3);...
305                       Epsil(3)^2;...
306                       Epsil(1)^3;...
307                       Epsil(1)^2*Epsil(2);...
308                       Epsil(1)^2*Epsil(3);...
309                       Epsil(2)^2*Epsil(1);...
310                       Epsil(2)^3;...
311                       Epsil(2)^2*Epsil(3);...
312                       Epsil(3)^2*Epsil(1);...
313                       Epsil(3)^2*Epsil(2);...
314                       Epsil(3)^3]'*Wv(:,kin);
315                     
316        end
317    end
318
319    function [valVtlx] = Vtildedx(kin,Epsil)
320        if(kin == K)
321            valVtlx = hdx;
322        else
323             valVtlx = [0                         0                       0;...
324                       1                          0                       0;...
325                       0                          1                       0;...
326                       0                          0                       1;...
327                       2*Epsil(1)                 0                       0;...
328                       Epsil(2)                   Epsil(1)                0;...
329                       Epsil(3)                   0                       Epsil(1);...
330                       0                          2*Epsil(2)              0;...
331                       0                          Epsil(3)                Epsil(2);...
332                       0                          0                       2*Epsil(3);...
333                       3*Epsil(1)^2               0                       0;...
334                       2*Epsil(1)*Epsil(2)        Epsil(1)^2              0;...
335                       2*Epsil(1)*Epsil(3)        0                       Epsil(1)^2;...
336                       Epsil(2)^2                 2*Epsil(2)*Epsil(1)     0;...
337                       0                          3*Epsil(2)^2            0;...
338                       0                          2*Epsil(2)*Epsil(3)     Epsil(2)^2;...
339                       Epsil(3)^2                 0                       2*Epsil(3)*Epsil(1);...
340                       0                          Epsil(3)^2              2*Epsil(3)*Epsil(2);...
341                       0                          0                       3*Epsil(3)^2]'...
342                       *Wv(:,kin);
343                     
344        end
345    end
346
347    function [valVtlxx] = Vtildedxdx(kin,Epsil)
348        if(kin == K)
349            valVtlxx = hdxdx;
350        else
351            valVtlxx = zeros(3,3);
352            valVtlxx(1,1) = 2*Wv(5,kin) + 6*Epsil(1)*Wv(11,kin) + 2*Epsil(2)*Wv(12,kin) + 2*Epsil(3)*Wv(13,kin);
353            valVtlxx(2,2) = 2*Wv(8,kin) + 2*Epsil(1)*Wv(14,kin) + 6*Epsil(2)*Wv(15,kin) + 2*Epsil(3)*Wv(16,kin);
354            valVtlxx(3,3) = 2*Wv(10,kin) + 2*Epsil(1)*Wv(17,kin) + 2*Epsil(2)*Wv(18,kin) + 6*Epsil(3)*Wv(19,kin);
355
356            valVtlxx(1,2) = Wv(6,kin) + 2*Epsil(1)*Wv(12,kin) + 2*Epsil(2)*Wv(14,kin);
357            valVtlxx(1,3) = Wv(7,kin) + 2*Epsil(1)*Wv(13,kin) + 2*Epsil(3)*Wv(17,kin);
358            valVtlxx(2,3) = Wv(4,kin) + 2*Epsil(2)*Wv(16,kin) + 2*Epsil(3)*Wv(18,kin);
359
360            valVtlxx(2,1) = valVtlxx(1,2);
361            valVtlxx(3,1) = valVtlxx(1,3);
362            valVtlxx(3,2) = valVtlxx(2,3);
363        end
364    end
365   
366keyboard
367end
Note: See TracBrowser for help on using the browser.