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