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