1% Compute the expected error variance of the analysis.
2%
3% err = divand_error(s)
4%
5% Compute the expected error variance of the analysis based on the
6% background error covariance, observation error covariance and data
7% distribution.
8%
9% Input:
10%   s: structure created by divand_factorize
11%
12% Output:
13%   err: expected error variance of the analysis
14
15
16function err = divand_error(s)
17
18H = s.H;
19sv = s.sv;
20N = sum(s.mask(:)==1);
21n = s.n;
22
23errp = zeros(N,1);
24
25if s.primal
26  errp = diag(s.P);
27
28  %global P_CL
29  %errp = diag(P_CL);
30else
31  %dual
32  B = s.B;
33  R = s.R;
34
35  % fun(x) computes (H B H' + R)*x
36  fun = @(x) H * (B * (H'*x)) + R*x;
37  t0 = cputime;
38
39  for i=1:N
40    % show progress every 5 seconds
41    if (t0 + 5 < cputime)
42      t0 = cputime;
43      fprintf('%d/%d\n',i,N);
44    end
45
46    e = zeros(N,1);
47    e(i) = 1;
48    % analyzed covariance
49    % P = B - B * H' * inv(H * B * H' + R) * H * B
50    % analyzed variance
51    % P = e' B e - e' * B * H' * inv(H * B * H' + R) * H * B * e
52
53    Be = B * e;
54    errp(i) = e' * Be;
55    HBe = H*Be;
56    %keyboard
57
58    % tmp is inv(H * B * H' + R) * H * B * e
59    %tic
60    [tmp,flag,relres,iter] = pcg(fun, HBe,s.tol,s.maxit,s.funPC);
61    %toc
62
63    if (flag ~= 0)
64        error('divand:pcg', ['Preconditioned conjugate gradients method'...
65            ' did not converge %d %g %g'],flag,relres,iter);
66    end
67
68    errp(i) = errp(i) - HBe' * tmp;
69  end
70end
71
72err = statevector_unpack(sv,errp);
73err(~s.mask) = NaN;
74
75% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be>
76%
77% This program is free software; you can redistribute it and/or modify it under
78% the terms of the GNU General Public License as published by the Free Software
79% Foundation; either version 2 of the License, or (at your option) any later
80% version.
81%
82% This program is distributed in the hope that it will be useful, but WITHOUT
83% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
84% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
85% details.
86%
87% You should have received a copy of the GNU General Public License along with
88% this program; if not, see <http://www.gnu.org/licenses/>.
89