1% Compute a variational analysis of arbitrarily located observations. 2% 3% [fi,err,s] = divand(mask,pmn,xi,x,f,len,lambda,...); 4% [fi,err] = divand(mask,pmn,xi,x,f,len,lambda,...); 5% [fi,s] = divand(mask,pmn,xi,x,f,len,lambda,...); 6% 7% Perform an n-dimensional variational analysis of the observations f located at 8% the coordinates x. The array fi represent the interpolated field at the grid 9% defined by the coordinates xi and the scales factors pmn. 10% 11% Input: 12% mask: binary mask delimiting the domain. 1 is inside and 0 outside. 13% For oceanographic application, this is the land-sea mask. 14% 15% pmn: scale factor of the grid. pmn is a cell array with n elements. Every 16% element represents the scale factor of the corresponding dimension. Its 17% inverse is the local resolution of the grid in a particular dimension. 18% 19% xi: cell array with n elements. Every element represents a coordinate 20% of the final grid on which the observations are interpolated 21% x: cell array with n elements. Every element represents a coordinate of 22% the observations 23% f: value of the observations *minus* the background estimate (m-by-1 array). 24% (see note) 25% len: correlation length 26% lambda: signal-to-noise ratio of observations (if lambda is a scalar). 27% The larger this value is, the closer is the field fi to the 28% observation. 29% if lambda is a scalar: 30% R = 1/lambda I, where R is the observation error covariance 31% matrix), 32% if lambda is a vector 33% a vector (R = diag(lambda)) or a matrix or a CovarParam object (R = 34% lambda). 35% 36% 37% Optional input arguments specified as pairs of keyword and values: 38% 'velocity', vel: velocity of advection constraint. The default is 39% no-advection constraint 40% 41% 'alpha': alpha is vector of coefficients multiplying various terms in the 42% cost function. The first element multiplies the norm. 43% The other i-th element of alpha multiplies the (i+1)-th derivative. 44% Per default, the highest derivative is m = ceil(1+n/2) where n is the 45% dimension of the problem. 46% 47% The values of alpha is the (m+1)th row of the Pascal triangle: 48% m=0 1 49% m=1 1 1 50% m=1 1 2 1 (n=1,2) 51% m=2 1 3 3 1 (n=3,4) 52% ... 53% 54% 'diagnostics': 0 or 1 turns diagnostic and debugging information on (1) or 55% off (0, default). If on, they will be returned as the last output 56% argument 57% 58% 'EOF', EOF: sub-space constraint. Orthogonal (EOF' WE^2 EOF = I) (units of 59% EOF: m^(-n/2)) 60% 61% 'EOF_scaling', EOF_scaling: (dimensional) 62% 63% 'constraint': a structure with user specified constrain 64% 65% 'moddim': modulo for cyclic dimension (vector with n elements). 66% Zero is used for non-cyclic dimensions. Halo points should 67% not be included for cyclic dimensions. For example if the first dimension 68% is cyclic, then the grid point corresponding to mask(1,j) should be 69% between mask(end,1) (left neighbor) and mask(2,j) (right neighbor) 70% 71% 'fracdim': fractional indices (n-by-m array). If this array is specified, 72% then x and xi are not used. 73% 74% 'inversion': direct solver ('chol' for Cholesky factorization) or a 75% interative solver ('pcg' for preconditioned conjugate gradient) can be 76% used. 77% 78% 'compPC': function that returns a preconditioner for the primal formulation 79% if inversion is set to 'pcg'. The function has the following arguments: 80% 81% [M1,M2] = compPC(iB,H,R) 82% 83% where iB is the inverse background error covariance, H the observation 84% operator and R the error covariance of the observation. The used 85% preconditioner M will be M = M1 * M2 (or just M = M1 if M2 is empty). 86% Per default a modified incomplete Cholesky factorization will be used a 87% preconditioner. 88% 89% Note: 'velocity' and 'constraint' may appear multiple times 90% 91% Output: 92% fi: the analysed field 93% err: error variance of the analysis field relative to the error variance of 94% the background 95% s: structure 96% s.iB: adimensional 97% s.E: scaled EOF (dimensional) 98% 99% Note: 100% If zero is not a valid first guess for your variable (as it is the case for 101% e.g. ocean temperature), you have to subtract the first guess from the 102% observations before calling divand and then add the first guess back in. 103% 104% Example: 105% see divand_simple_example.m 106% 107 108 109function varargout = divand(mask,pmn,xi,x,f,len,lambda,varargin) 110 111% default values 112 113velocity = {}; 114EOF = []; 115diagnostics = 0; 116EOF_lambda = 0; 117primal = 1; 118factorize = 1; 119tol = 1e-6; 120maxit = 100; 121minit = 10; 122constraints = {}; 123inversion = 'chol'; 124moddim = []; 125fracindex = []; 126alpha = []; 127keepLanczosVectors = 0; 128compPC = @divand_pc_sqrtiB; 129 130prop = varargin; 131for i=1:length(prop)/2 132 if strcmp(prop{2*i-1},'velocity') 133 velocity{end+1} = prop{2*i}; 134 elseif strcmpi(prop{2*i-1},'EOF') 135 EOF = prop{2*i}; 136 elseif strcmpi(prop{2*i-1},'EOF_scaling') 137 EOF_lambda = prop{2*i}; 138 elseif strcmpi(prop{2*i-1},'primal') 139 primal = prop{2*i}; 140 elseif strcmpi(prop{2*i-1},'factorize') 141 factorize = prop{2*i}; 142 elseif strcmpi(prop{2*i-1},'inversion') 143 inversion = prop{2*i}; 144 elseif strcmpi(prop{2*i-1},'tol') 145 tol = prop{2*i}; 146 elseif strcmpi(prop{2*i-1},'maxit') 147 maxit = prop{2*i}; 148 elseif strcmpi(prop{2*i-1},'minit') 149 minit = prop{2*i}; 150 elseif strcmpi(prop{2*i-1},'diagnostics') 151 diagnostics = prop{2*i}; 152 elseif strcmpi(prop{2*i-1},'constraint') 153 constraints{end+1} = prop{2*i}; 154 elseif strcmpi(prop{2*i-1},'fracindex') 155 fracindex = prop{2*i}; 156 elseif strcmpi(prop{2*i-1},'moddim') 157 moddim = prop{2*i}; 158 elseif strcmpi(prop{2*i-1},'alpha') 159 alpha = prop{2*i}; 160 elseif strcmpi(prop{2*i-1},'keepLanczosVectors') 161 keepLanczosVectors = prop{2*i}; 162 elseif strcmpi(prop{2*i-1},'compPC') 163 compPC = prop{2*i}; 164 else 165 warning('divand:unknownProperty','unknown property "%s" is irgnored',prop{2*i-1}); 166 end 167end 168 169% check inputs 170 171if ~any(mask(:)) 172 error('no sea points in mask'); 173end 174 175s = divand_background(mask,pmn,len,alpha,moddim); 176s.betap = 0; 177s.EOF_lambda = EOF_lambda; 178s.primal = primal; 179s.factorize = factorize; 180s.tol = tol; 181s.maxit = maxit; 182s.minit = minit; 183s.inversion = inversion; 184s.keepLanczosVectors = keepLanczosVectors; 185s.compPC = compPC; 186 187% remove non-finite elements from observations 188f = f(:); 189valid = isfinite(f); 190x = cat_cell_array(x); 191 192if ~all(valid) 193 fprintf(1,'remove %d (out of %d) non-finite elements from observation vector\n',sum(~valid),numel(f)); 194 x = reshape(x,[length(f) s.n]); 195 f = f(valid); 196 x = reshape(x(repmat(valid,[1 s.n])),[length(f) s.n]); 197 198 if ~isempty(fracindex) 199 fracindex = fracindex(:,valid); 200 end 201 202 if isscalar(lambda) 203 % do nothing 204 elseif isvector(lambda) 205 lambda = lambda(valid); 206 elseif ismatrix(lambda) 207 lambda = lambda(valid,valid); 208 end 209end 210 211apply_EOF_contraint = ~(isempty(EOF) | all(EOF_lambda == 0)); 212 213s.mode = 1; 214 215if ~apply_EOF_contraint 216 s.betap = 0; 217else 218 if s.mode==0 219 s.betap = max(EOF_lambda)/s.coeff; % units m^(-n) 220 elseif s.mode==1 221 s.betap = max(max(EOF_lambda)-1,0)/s.coeff; 222 end 223end 224 225%assert(s.betap,0,1e-8) 226% increase contraint on total enegery to ensure system is still positive defined 227%s.betap 228s.iB = s.iB + s.betap * s.WE'*s.WE; 229 230% add observation constrain to cost function 231s = divand_addc(s,divand_obs(s,xi,x,f,lambda,fracindex)); 232 233% add advection constraint to cost function 234for i=1:length(velocity) 235 %s = divand_advection(s,velocity); 236 s = divand_addc(s,divand_constr_advec(s,velocity{i})); 237end 238 239 240% add all additional constrains 241for i=1:length(constraints) 242 s = divand_addc(s,constraints{i}); 243end 244 245 246if apply_EOF_contraint 247 s = divand_eof_contraint(s,EOF_lambda,EOF); 248end 249 250% factorize a posteori error covariance matrix 251% or compute preconditioner 252s = divand_factorize(s); 253 254if ~apply_EOF_contraint 255 [fi,s] = divand_solve(s,f); 256else 257 [fi,s] = divand_solve_eof(s,f); 258end 259 260varargout{1} = fi; 261 262if nargout-diagnostics >= 2 263 err = divand_error(s); 264 varargout{2} = err; 265end 266 267if diagnostics 268 s.B = CovarIS(s.iB); 269 [s.Jb,s.Jo,s.Jeof,s.J] = divand_diagnose(s,fi,f); 270 s.valid = valid; 271 varargout{nargout} = s; 272end 273 274 275% Copyright (C) 2008-2014 Alexander Barth <barth.alexander@gmail.com> 276% 277% This program is free software; you can redistribute it and/or modify it under 278% the terms of the GNU General Public License as published by the Free Software 279% Foundation; either version 2 of the License, or (at your option) any later 280% version. 281% 282% This program is distributed in the hope that it will be useful, but WITHOUT 283% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 284% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 285% details. 286% 287% You should have received a copy of the GNU General Public License along with 288% this program; if not, see <http://www.gnu.org/licenses/>. 289 290% LocalWords: fi divand pmn len diag CovarParam vel ceil moddim fracdim