root/applications/dual/IterativeLocal/ildp.m @ 1179

Revision 755, 10.3 kB (checked in by vahalam, 15 years ago)
Line 
1function ildp
2
3tic
4    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5    %pocatecni konstanty
6   
7%     [b=1, yr=1] - to jsme zkouseli spolu  sigma = 0.1; rho = 0.5;
8%     [b=-1, yr=1] sigma = 0.1; rho = 0.5; !aprior
9%     [b=10, yr=5] sigma = 0.5; rho = 2.0;
10%     [b=0.6, yr=10] sigma = 0.1; rho = 1.0;
11   
12    Iterace = 10; %iterace
13    K = 5; %casy
14    N = 100; %vzorky
15   
16    sigma = 0.1;
17    Sigmas = [[sigma^2 0 0]; [0 0 0]; [0 0 0]];
18   
19    h = 0;
20    hdx = [0; 0; 0];
21    hdxdx = [[0 0 0]; [0 0 0]; [0 0 0]];
22   
23    Rk = 1*ones(1, K);
24    x0 = [0; 1; 1];
25   
26    %velikost okoli pro lokalni metodu
27    rho = 0.5;
28    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
29    %globalni promenne
30
31    Kpi = ones(4, K);
32        Kpi(4, :) = zeros(1, K);
33%         Kpi(3, :) = zeros(1, K);
34       
35    Wv = zeros(9, K);
36   
37    Xkn = zeros(3, K, N);
38    Xstripe = zeros(3, K);   
39 
40    gka = 0;
41    gnu = 0;
42   
43    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
44    %iteracni smycka
45
46    for i = 1:Iterace,
47       %generovani stavu
48       for n = 1:N,
49            Xkn(:, 1, n) = x0 + [sigma*randn(); 0; 0];
50            for k = 1:K-1,
51               Uk = uPi(k, Xkn(:, k, n));
52               Kk = Uk*Xkn(3, k, n)/(Uk^2*Xkn(3, k, n) + sigma^2);
53               Xkn(1, k+1, n) = Xkn(1, k, n) + Xkn(2, k, n)*Uk + sigma*randn();
54               Xkn(2, k+1, n) = Xkn(2, k, n) + Kk*(Xkn(1, k+1, n) - Xkn(1, k, n) - Xkn(2, k, n)*Uk);
55               Xkn(3, k+1, n) = (1 - Kk*Uk)*Xkn(3, k, n);
56            end           
57       end
58       Xstripe = mean(Xkn, 3);
59                 
60       for k = K-1:-1:1,
61            gka = k;
62%             1]
63            for n = 1:N,
64                %krive okoli
65               Xkn(1, k, n) = Xstripe(1, k) + rho*randn();
66               Xkn(2, k, n) = Xstripe(2, k) + rho*randn(); 
67               Xkn(3, k, n) = Xstripe(3, k)*exp(rho*randn());
68            end           
69%             2]
70            for n = 1:N,
71               gnu = n;
72                [Uopt(n), Hmin(n)] = fminunc(@Hamilt, uPi(k, Xkn(:, k, n)), optimset('GradObj','on','Display','notify'));
73            end
74%             3]
75            for n = 1:N,               
76               Vn(n) = Hmin(n) + Vtilde(k+1, Xkn(:, k, n));
77            end           
78%             4]
79            Epsilon = zeros(3, N);
80            for n = 1:N,
81                Epsilon(1, n) = Xkn(1, k, n) - Xstripe(1, k);
82                Epsilon(2, n) = Xkn(2, k, n) - Xstripe(2, k);
83                Epsilon(3, n) = Xkn(3, k, n)/Xstripe(3, k);
84            end
85            mFi = matrixFi(Epsilon);
86            FiFiTInvFi = (mFi*mFi')\mFi;
87            Wv(:,k) = FiFiTInvFi * Vn';
88            for n = 1:N,
89               yt(n) = Xkn(1, k, n);
90               bt(n) = Xkn(2, k, n);
91               pt(n) = Xkn(3, k, n);
92            end
93            mPsi = [yt',...
94                    bt'.*Uopt',...
95                    pt'.*Uopt',...
96                    Uopt'];
97            PsiPsiTInvPsi = (mPsi'*mPsi)\mPsi';             
98            Kpi(:,k) = PsiPsiTInvPsi * (Rk(k)*ones(N,1));           
99       end
100    end
101    %%%%%%%%%%%
102    toc
103   
104    Kpi
105   
106                for n = 1:N,
107                    Xkn(:, 1, n) = x0 + [sigma*randn(); 0; 0];
108                    for k = 1:K-1,
109                       Uk = uPi(k, Xkn(:, k, n));
110                            UU(k,n) = Uk;
111                       Kk = Uk*Xkn(3, k, n)/(Uk^2*Xkn(3, k, n) + sigma^2);
112                       Xkn(1, k+1, n) = Xkn(1, k, n) + Xkn(2, k, n)*Uk + sigma*randn();
113                       Xkn(2, k+1, n) = Xkn(2, k, n) + Kk*(Xkn(1, k+1, n) - Xkn(1, k, n) - Xkn(2, k, n)*Uk);
114                       Xkn(3, k+1, n) = (1 - Kk*Uk)*Xkn(3, k, n);       
115                    end
116                        hold all
117                        subplot(4,1,1);
118                        plot(1:K,Xkn(1,:,n))
119                        hold all
120                        subplot(4,1,2);
121                        plot(1:K,Xkn(2,:,n))
122                        hold all
123                        subplot(4,1,3);
124                        plot(1:K,Xkn(3,:,n))
125                        hold all
126                        subplot(4,1,4);
127                        plot(1:K-1,UU(:,n))
128
129                end
130                title('iLDP')
131                figure
132               for n = 1:N,
133                    Xkn(:, 1, n) = x0 + [sigma*randn(); 0; 0];
134                    for k = 1:K-1,
135                       Uk = (Rk(k) - Xkn(1, k, n))/(Xkn(2, k, n) + Xkn(3, k, n));
136                            UU(k,n) = Uk;
137                       Kk = Uk*Xkn(3, k, n)/(Uk^2*Xkn(3, k, n) + sigma^2);
138                       Xkn(1, k+1, n) = Xkn(1, k, n) + Xkn(2, k, n)*Uk + sigma*randn();
139                       Xkn(2, k+1, n) = Xkn(2, k, n) + Kk*(Xkn(1, k+1, n) - Xkn(1, k, n) - Xkn(2, k, n)*Uk);
140                       Xkn(3, k+1, n) = (1 - Kk*Uk)*Xkn(3, k, n);     
141                    end
142                        hold all
143                        subplot(4,1,1);
144                        plot(1:K,Xkn(1,:,n))
145                        hold all
146                        subplot(4,1,2);
147                        plot(1:K,Xkn(2,:,n))
148                        hold all
149                        subplot(4,1,3);
150                        plot(1:K,Xkn(3,:,n))
151                        hold all
152                        subplot(4,1,4);
153                        plot(1:K-1,UU(:,n))
154
155               end
156               title('CE')
157               
158               for vzorek = 1:100,
159                   loss(vzorek) = 0;                   
160                   bb = randn() + x0(2);
161                   yy(1) = x0(1);
162                   for k=1:K-1,
163                       yy(k+1)=yy(k)+bb*uPi(k,[yy(k); bb; 0]);
164                       loss(vzorek) = loss(vzorek) + (yy(k+1) - Rk(k+1))^2;
165                   end                 
166               end
167               figure
168               hist(log(loss))
169               
170    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
171    %pomocne funkce
172   
173    function [val_uPi] = uPi(k_uPi, x_uPi)
174        val_uPi = (Rk(k_uPi) - Kpi(1, k_uPi)*x_uPi(1))/(Kpi(2, k_uPi)*x_uPi(2) + Kpi(3, k_uPi)*x_uPi(3) + Kpi(4, k_uPi));
175
176    end
177
178    function [val_ham, val_grad] = Hamilt(u_ham)
179%         Vtddx = Vtilde_dx(gka+1, Xkn(:, gka, gnu));
180       val_ham = (Xkn(1, gka, gnu) + Xkn(2, gka, gnu)*u_ham - Rk(gka+1))^2 + Xkn(3, gka, gnu)*u_ham^2 ... + sigma^2 ...  %ztrata l
181               + [Xkn(2, gka, gnu)*u_ham; ...                                   
182                  Xkn(3, gka, gnu)*u_ham*(Xkn(1, gka+1, gnu) - Xkn(1, gka, gnu) - Xkn(2, gka, gnu)*u_ham)/(Xkn(3, gka, gnu)*u_ham^2 + sigma^2); ...
183                  -Xkn(3, gka, gnu)^2*u_ham^2/(Xkn(3, gka, gnu)*u_ham^2 + sigma^2)]' ... %fce f               
184                *Vtilde_dx(gka+1, Xkn(:, gka, gnu)) ...
185                + Wv(5, gka+1)*sigma;%+ Wv(4, gka+1)*sigma;
186           
187       val_grad = 2*(Xkn(1, gka, gnu) + Xkn(2, gka, gnu)*u_ham - Rk(gka+1))*Xkn(2, gka, gnu) + 2*Xkn(3, gka, gnu)*u_ham ...  %ztrata du
188                + [Xkn(2, gka, gnu);...                                                                                                                                                                                                                                                     
189                   (2*u_ham^2*Xkn(3, gka, gnu)^2*(Xkn(1, gka, gnu) - Xkn(1, gka+1, gnu) + u_ham*Xkn(2, gka, gnu)))/(sigma^2 + u_ham^2*Xkn(3, gka, gnu))^2 - (u_ham*Xkn(2, gka, gnu)*Xkn(3, gka, gnu))/(sigma^2 + u_ham^2*Xkn(3, gka, gnu)) - (Xkn(3, gka, gnu)*(Xkn(1, gka, gnu) - Xkn(1, gka+1, gnu) + u_ham*Xkn(2, gka, gnu)))/(sigma^2 + u_ham^2*Xkn(3, gka, gnu));...
190                   (2*u_ham^3*Xkn(3, gka, gnu)^3)/(sigma^2 + u_ham^2*Xkn(3, gka, gnu))^2 - (2*u_ham*Xkn(3, gka, gnu)^2)/(sigma^2 + u_ham^2*Xkn(3, gka, gnu))]' ... %fce f du
191                * Vtilde_dx(gka+1, Xkn(:, gka, gnu)); %derivace Bellman fce               
192    end
193
194    function [val_Vt] = Vtilde(k_Vt, x_Vt)
195        if(k_Vt == K)
196            val_Vt = h;
197        else
198            Epsl = zeros(3, 1);           
199                Epsl(1) = x_Vt(1) - Xstripe(1, k_Vt);
200                Epsl(2) = x_Vt(2) - Xstripe(2, k_Vt);
201                Epsl(3) = x_Vt(3)/Xstripe(3, k_Vt);
202           
203            val_Vt = vectFi(Epsl)' * Wv(:,k_Vt);
204        end
205    end
206
207    function [val_Vt] = Vtilde_dx(k_Vt, x_Vt)
208        if(k_Vt == K)
209            val_Vt = hdx;
210        else
211             Epsl = zeros(3, 1);           
212                Epsl(1) = x_Vt(1) - Xstripe(1, k_Vt);
213                Epsl(2) = x_Vt(2) - Xstripe(2, k_Vt);
214                Epsl(3) = x_Vt(3)/Xstripe(3, k_Vt);
215               
216            val_Vt = difFi(Epsl)' * Wv(:,k_Vt);
217        end
218    end
219
220    function [val_Fi] = vectFi(x_Fi)
221        val_Fi = [ ...
222                   1; ...
223                   x_Fi(1); ...
224                   x_Fi(2); ...
225                   log(x_Fi(3)); ...
226                   x_Fi(1)^2; ...
227                   x_Fi(1)*x_Fi(2); ...
228                   x_Fi(1)*log(x_Fi(3)); ...
229                   x_Fi(2)^2; ...
230                   x_Fi(2)*log(x_Fi(3)); ...
231%                    2*ln(x_Fi(3)); ...
232                 ];
233    end
234
235    function [val_Fi] = matrixFi(x_Fi)
236        val_Fi = [ ...
237                   ones(1, N); ...
238                   x_Fi(1, :); ...
239                   x_Fi(2, :); ...
240                   log(x_Fi(3, :)); ...
241                   x_Fi(1, :).^2; ...
242                   x_Fi(1, :).*x_Fi(2, :); ...
243                   x_Fi(1, :).*log(x_Fi(3, :)); ...
244                   x_Fi(2, :).^2; ...
245                   x_Fi(2, :).*log(x_Fi(3, :)); ...
246                 ];
247
248    end
249
250    function [val_Fi] = difFi(x_Fi)
251        val_Fi = [ ...
252                   0            0           0; ...
253                   1            0           0; ...
254                   0            1           0; ...
255                   0            0           1/(x_Fi(3)); ...
256                   2*x_Fi(1)    0           0; ...
257                   x_Fi(2)      x_Fi(1)     0; ...
258                   log(x_Fi(3))    0           x_Fi(1)/(x_Fi(3)); ...
259                   0            2*x_Fi(2)   0; ...
260                   0            log(x_Fi(3))   x_Fi(2)/(x_Fi(3)); ...
261                 ];
262    end
263end
Note: See TracBrowser for help on using the browser.