source: projects/ckbs/trunk/test/sumsq_grad_ok.m @ 95

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

Revmove reference to gpl2.txt and further shorten copyright and license statement

File size: 2.2 KB
Line 
1% -------------------------------------------------------------------
2% ckbs: Constrained Kalman-Bucy Smoother Program: Copyright (C) 2006
3% Authors: Bradlely Bell:        bradbell at washington dot edu
4%          Gianluigi Pillonetto: giapi at dei dot unipd dot it
5% License: GNU General Public License Version 2
6% -------------------------------------------------------------------
7% $begin sumsq_grad_ok.m$$ $newlinech %$$
8% $spell
9%       ckbs
10%       dg
11%       dh
12%       diff
13%       end end
14%       grad grad
15%       obj
16%       qinv
17%       rinv
18%       sm
19%       sumsq
20%       tmp
21%       xm
22%       xp
23% $$
24%
25% $section ckbs_sumsq_grad Example and Test$$
26%
27% $index ckbs_sumsq_grad, example and test$$
28% $index sumsq_grad, example and test$$
29% $index example, sumsq_grad$$
30% $index test, sumsq_grad$$
31%
32% $head Source Code$$
33% $newlinech $$ $codep
34function [ok] = sumsq_grad_ok()
35ok = true;
36% --------------------------------------------------------
37% You can change these parameters
38m    = 1;   % number of measurements per time point
39n    = 2;   % number of state vector components per time point
40N    = 3;   % number of time points
41% ---------------------------------------------------------
42%  Define the problem
43rand('seed', 123);
44x    = rand(n, N);
45z    = rand(m, N);
46h    = rand(m, N);
47g    = rand(n, N);
48dg   = zeros(n, n, N);
49dh   = zeros(m, n, N);
50qinv = zeros(n, n, N);
51rinv = zeros(m, m, N);
52for k = 1 : N
53        dh(:, :, k)   = rand(m, n);
54        dg(:, :, k)   = rand(n, n);
55        tmp           = rand(m, m);
56        rinv(:, :, k) = (tmp + tmp') / 2 + 2 * eye(m);
57        tmp           = rand(n, n);
58        qinv(:, :, k) = (tmp + tmp') / 2 + 2 * eye(n);
59end
60% ---------------------------------------------------------
61% Compute the gradient using ckbs_sumsq_grad
62grad = ckbs_sumsq_grad(x, z, g, h, dg, dh, qinv, rinv);
63% ---------------------------------------------------------
64% Use finite differences to check gradient
65step   = 1;
66for k  = 1 : N
67for i  = 1 : n
68        % Check second partial w.r.t x(i,k)
69        xm        = x;
70        xm(i, k)  = xm(i, k) - step;
71        Sm        = ckbs_sumsq_obj(xm, z, g, h, dg, dh, qinv, rinv);
72        %
73        xp        = x;
74        xp(i, k)  = xp(i, k) + step;
75        Sp        = ckbs_sumsq_obj(xp, z, g, h, dg, dh, qinv, rinv);
76        %
77        check  = (Sp - Sm) / ( 2 * step);
78        diff   = grad(i, k) - check;
79        ok     = ok & ( abs(diff) < 1e-10 );
80end
81end
82return
83end
84% $$ $newlinech %$$
85% $end
Note: See TracBrowser for help on using the repository browser.