# source:projects/ckbs/trunk/test/newton_step_ok.m@95

Last change on this file since 95 was 35, checked in by bradbell, 11 years ago

File size: 2.5 KB
Line
1% -------------------------------------------------------------------
2% ckbs: Constrained Kalman-Bucy Smoother Program: Copyright (C) 2006
4%          Gianluigi Pillonetto: giapi at dei dot unipd dot it
6% -------------------------------------------------------------------
7% \$begin newton_step_ok.m\$\$ \$newlinech %\$\$
8% \$spell
9%       ckbs
10%       Bdia
11%       blk
12%       ds
13%       du
14%       dy
15%       end end
16%       hdia
17%       hlow
18%       mu
19%       tmp
20% \$\$
21%
22% \$section ckbs_newton_step Example and Test\$\$
23%
24% \$index ckbs_newton_step, example and test\$\$
25% \$index newton_step, example and test\$\$
26% \$index example, newton_step\$\$
27% \$index test, newton_step\$\$
28%
30% \$newlinech \$\$ \$codep
31function [ok] = newton_step_ok()
32ok = true;
33% --------------------------------------------------------
34% You can change these parameters
35m    = 2;   % number of measurements per time point
36n    = 3;   % number of state vector components per time point
37N    = 4;   % number of time points
38% ---------------------------------------------------------
39%  Define the problem
40rand('seed', 123);
41mu    = .5;
42r     = m * N;
43p     = n * N;
44s     = rand(r, 1) + 1;
45y     = rand(p, 1);
46u     = rand(r, 1) + 1;
47b     = rand(r, 1);
48d     = rand(p, 1);
49H     = zeros(p, p);
50B     = zeros(r, r);
51Hdia  = zeros(n, n, N);
52Hlow  = zeros(n, n, N);
53Bdia  = zeros(m, n, N);
54blk_m = m * N + (1 : m);
55blk_n = n * N + (1 : n);
56for k = N : -1 : 1
57        blk_n = blk_n - n;
58        blk_m = blk_m - m;
59        %
60        Bdia(:,:, k)    = rand(m, n);
61        tmp             = rand(n, n);
62        Hlow(:,:, k)    = tmp;
63        Hdia(:,:, k)   = (tmp * tmp') + 4 * eye(n);
64        %
65        B(blk_m, blk_n) = Bdia(:,:, k);
66        H(blk_n, blk_n) = Hdia(:,:, k);
67        if k > 1
68                H(blk_n, blk_n - n) = Hlow(:,:, k);
69        end
70        if k < N
71                H(blk_n, blk_n + n) = Hlow(:,:, k + 1)';
72        end
73end
74% -------------------------------------------------------------------
75[ds, dy, du] = ckbs_newton_step(mu, s, y, u, b, d, Bdia, Hdia, Hlow);
76% -------------------------------------------------------------------
77F      = [ ...
78        s + b + B * y
79        H * y + B' * u + d
80        s .* u - mu ...
81];
82dF     = [ ...
83        eye(r)      , B           , zeros(r, r)
84        zeros(p, r) , H           , B'
85        diag(u)     , zeros(r, p) , diag(s) ...
86];
87delta        = [ ds ; dy ; du ];
88ok           = ok & ( max( abs( F + dF * delta ) ) <= 1e-10 );
89return
90end
91% \$\$ \$newlinech %\$\$
92% \$end
Note: See TracBrowser for help on using the repository browser.