Home > atmlab > retrieval > oem.m

oem

PURPOSE ^

OEM Data inversion by "optimal estimation method"

SYNOPSIS ^

function [X,R] = oem(O,Q,R,comfun,Sx,Se,Sxinv,Seinv,xa,y)

DESCRIPTION ^

 OEM   Data inversion by "optimal estimation method"

   The name OEM is here kept as it is so well established in the atmospheric
   community. A better name would be "the maximum a posteriori solution". See
   "Inverse methods for atmospheric sounding: Theory and practice" by C.D.
   Rodgers for background and theory.

   The OEM calculations are controlled by the structure O. The fields of O,
   and associated descriptions, are listed by
      qinfo(@oem)

   The Marquardt-Levenberg method is implemented following Eq. 5.36 of
   Rodgers' book. The gamma factor, here denoted as ga, is not allowed to
   exceed O.ga_max. When a succesful iteration has been performed, ga is
   decreased following O.ga_factor_ok. For cases when a lower cost is not
   obtained, ga is increased by O.ga_factor_not_ok. Values of ga below 1 will
   be treated as zero. If not a lower cost is obtained for ga=0, then ga=1 is
   tested. O.ga_start is used as start value for ga.

   The results are gathered in a structure X, with fields:
      x         : Retrieved state.
      converged : Convergence flag. 1 if convergency criterion has been 
                  reached. Otherwise 0 for non-linear inversions when the
                  maximum number of iterations is reached, or -1 for 
                  Marquardt-Levenberg when the maximum value for gamma is
                  passed. Set to NaN for linear retrievals.
      cost      : Total cost, normalised by length of y. A succesful retrieval
                  shall then have a cost in the order of 1.
      cost_y    : Cost corresponding to the measurement vector.
      cost_x    : Cost corresponding to the state vector.
      dx        : Change in the state vector x between each iteration.
                  A normalisation to the length of x is applied.
      ga        : The Marquardt-Levenberg parameter for each iteration.
      e         : Total retrieval error (square root of diagonal elements 
                  of S).
      eo        : Observation error (square root of diagonal elements of So).
      es        : Smoothing error (square root of diagonal elements of Ss).
      ex        : A priori uncertainty (square root of diag elements of Sx).
      yf        : Fitted spectrum.
      A         : Averaging kernel matrix.
      G         : Gain matrix (also denoted as Dy).
      J         : Jacobian (weighting function matrix).
      S         : Covariance matrix for total error. That is, the sum
                  of So and Ss.
      So        : Covariance matrix for observation error.
      Ss        : Covariance matrix for smoothing error.
      Xiter     : The x-vector after each iteration.
   The fields *x* and *converged* are always included (if an actual inversion 
   is performed). Other fields are included if corresponding O field is set.
   Note that the case x=x0, where x0 is the start value of x, is treated as
   iteration 1, and the first value in variables covering all iterations 
   corresponds to the sitatuation when no inversion/iteration has been
   performed. 

   The communication with the forward model is handled by the function
   handle comfun. If comfun is set to @arts_oem, the following calls will 
   be made:
      [R,yf,J]  = arts_oem( Q, R, x, iter );
      [R,yf]    = arts_oem( Q, R, x, iter );
   That is, this function expects the comfun function to provide spectrum 
   and jacobian (only spectrum in the second case) for present state (the x 
   vector). The oem function does not use any information in Q and R, and
   the comfun function is free to use these variables freely for bookkeeping
   etc. This function handles y, x, Sx and Se as totally generic variables
   and unit conversion and similar tasks are left to the comfun and
   post-processing. 

   Use the first format to check default values for O. The second format
   is intended for extracting retrieval characteristics (corresponding to
   a linear inversion). The third format performs a complete inversion.

 FORMAT   O = oem
            or
          [X,R] = oem(O,Q,R,comfun,Sx,Se,Sxinv,Seinv,xa)
            or
          [X,R] = oem(O,Q,R,comfun,Sx,Se,Sxinv,Seinv,xa,y)
        
 OUT   O        O structure with default value for all recognised fields.
       X        Retrieved state and additional information.
       R        Retreival variable structure from last calculation.
 IN    O        OEM setting structure.
       Q        Forward model setting structure.
       R        Retreival variable structure.
       comfun   Handle to function for communication with forward model.
       Sx       Covariance matrix for a priori uncertainty.
       Se       Covariance matrix for observation uncertainties.
       Sxinv    The inverse of Sx. If set to [], the inverse is calculated 
                internally.
       Seinv    The inverse of Se. If set to [], the inverse is calculated 
                internally.
       xa       A priori profile.
       y        Measurement vector.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

DOWNLOAD ^

oem.m

SOURCE CODE ^

0001 % OEM   Data inversion by "optimal estimation method"
0002 %
0003 %   The name OEM is here kept as it is so well established in the atmospheric
0004 %   community. A better name would be "the maximum a posteriori solution". See
0005 %   "Inverse methods for atmospheric sounding: Theory and practice" by C.D.
0006 %   Rodgers for background and theory.
0007 %
0008 %   The OEM calculations are controlled by the structure O. The fields of O,
0009 %   and associated descriptions, are listed by
0010 %      qinfo(@oem)
0011 %
0012 %   The Marquardt-Levenberg method is implemented following Eq. 5.36 of
0013 %   Rodgers' book. The gamma factor, here denoted as ga, is not allowed to
0014 %   exceed O.ga_max. When a succesful iteration has been performed, ga is
0015 %   decreased following O.ga_factor_ok. For cases when a lower cost is not
0016 %   obtained, ga is increased by O.ga_factor_not_ok. Values of ga below 1 will
0017 %   be treated as zero. If not a lower cost is obtained for ga=0, then ga=1 is
0018 %   tested. O.ga_start is used as start value for ga.
0019 %
0020 %   The results are gathered in a structure X, with fields:
0021 %      x         : Retrieved state.
0022 %      converged : Convergence flag. 1 if convergency criterion has been
0023 %                  reached. Otherwise 0 for non-linear inversions when the
0024 %                  maximum number of iterations is reached, or -1 for
0025 %                  Marquardt-Levenberg when the maximum value for gamma is
0026 %                  passed. Set to NaN for linear retrievals.
0027 %      cost      : Total cost, normalised by length of y. A succesful retrieval
0028 %                  shall then have a cost in the order of 1.
0029 %      cost_y    : Cost corresponding to the measurement vector.
0030 %      cost_x    : Cost corresponding to the state vector.
0031 %      dx        : Change in the state vector x between each iteration.
0032 %                  A normalisation to the length of x is applied.
0033 %      ga        : The Marquardt-Levenberg parameter for each iteration.
0034 %      e         : Total retrieval error (square root of diagonal elements
0035 %                  of S).
0036 %      eo        : Observation error (square root of diagonal elements of So).
0037 %      es        : Smoothing error (square root of diagonal elements of Ss).
0038 %      ex        : A priori uncertainty (square root of diag elements of Sx).
0039 %      yf        : Fitted spectrum.
0040 %      A         : Averaging kernel matrix.
0041 %      G         : Gain matrix (also denoted as Dy).
0042 %      J         : Jacobian (weighting function matrix).
0043 %      S         : Covariance matrix for total error. That is, the sum
0044 %                  of So and Ss.
0045 %      So        : Covariance matrix for observation error.
0046 %      Ss        : Covariance matrix for smoothing error.
0047 %      Xiter     : The x-vector after each iteration.
0048 %   The fields *x* and *converged* are always included (if an actual inversion
0049 %   is performed). Other fields are included if corresponding O field is set.
0050 %   Note that the case x=x0, where x0 is the start value of x, is treated as
0051 %   iteration 1, and the first value in variables covering all iterations
0052 %   corresponds to the sitatuation when no inversion/iteration has been
0053 %   performed.
0054 %
0055 %   The communication with the forward model is handled by the function
0056 %   handle comfun. If comfun is set to @arts_oem, the following calls will
0057 %   be made:
0058 %      [R,yf,J]  = arts_oem( Q, R, x, iter );
0059 %      [R,yf]    = arts_oem( Q, R, x, iter );
0060 %   That is, this function expects the comfun function to provide spectrum
0061 %   and jacobian (only spectrum in the second case) for present state (the x
0062 %   vector). The oem function does not use any information in Q and R, and
0063 %   the comfun function is free to use these variables freely for bookkeeping
0064 %   etc. This function handles y, x, Sx and Se as totally generic variables
0065 %   and unit conversion and similar tasks are left to the comfun and
0066 %   post-processing.
0067 %
0068 %   Use the first format to check default values for O. The second format
0069 %   is intended for extracting retrieval characteristics (corresponding to
0070 %   a linear inversion). The third format performs a complete inversion.
0071 %
0072 % FORMAT   O = oem
0073 %            or
0074 %          [X,R] = oem(O,Q,R,comfun,Sx,Se,Sxinv,Seinv,xa)
0075 %            or
0076 %          [X,R] = oem(O,Q,R,comfun,Sx,Se,Sxinv,Seinv,xa,y)
0077 %
0078 % OUT   O        O structure with default value for all recognised fields.
0079 %       X        Retrieved state and additional information.
0080 %       R        Retreival variable structure from last calculation.
0081 % IN    O        OEM setting structure.
0082 %       Q        Forward model setting structure.
0083 %       R        Retreival variable structure.
0084 %       comfun   Handle to function for communication with forward model.
0085 %       Sx       Covariance matrix for a priori uncertainty.
0086 %       Se       Covariance matrix for observation uncertainties.
0087 %       Sxinv    The inverse of Sx. If set to [], the inverse is calculated
0088 %                internally.
0089 %       Seinv    The inverse of Se. If set to [], the inverse is calculated
0090 %                internally.
0091 %       xa       A priori profile.
0092 %       y        Measurement vector.
0093 
0094 function [X,R] = oem(O,Q,R,comfun,Sx,Se,Sxinv,Seinv,xa,y)
0095 
0096   
0097 %=============================================================================
0098 %=== Definition of O
0099 %=============================================================================
0100   
0101 if ~nargin
0102 
0103 %-----------------------------------------------------------------------------
0104 O.A = false;
0105 I.A = [ ...
0106 'Flag to include averaging kernel matrix in X. Default is 0.'
0107 ];
0108 %-----------------------------------------------------------------------------
0109 O.cost = false;
0110 I.cost = [ ...
0111 'Flag to include cost values in X and to calculate cost even if solution ',...
0112 'method does not require cost values. This affects also the output to ',...
0113 '*outfids*.'
0114 ];
0115 %-----------------------------------------------------------------------------
0116 O.dx = false;
0117 I.dx = [ ...
0118 'Flag to include the sequence of convergence values (defined as for *stop_dx*).'
0119 ];
0120 %-----------------------------------------------------------------------------
0121 O.e = false;
0122 I.e = [ ...
0123 'Flag to include in X the estimate of total retrieval error  (square root ',...
0124 'of diagonal elements of S). That is, the standard deviation for the sum ',...
0125 'of observation and smoothing errors.'
0126 ];
0127 %-----------------------------------------------------------------------------
0128 O.eo = false;
0129 I.eo = [ ...
0130 'Flag to include in X the estimate of observation error (square root of ',...
0131 'diagonal elements of So).'
0132 ];
0133 %-----------------------------------------------------------------------------
0134 O.es = false;
0135 I.es = [ ...
0136 'Flag to include in X the estimate of smootning error (square root of ',...
0137 'diagonal elements of Ss).'
0138 ];
0139 %-----------------------------------------------------------------------------
0140 O.ex = false;
0141 I.ex = [ ...
0142 'Flag to include in X the a priori uncertainty (square root of ',...
0143 'diagonal elements of Sx).'
0144 ];
0145 %-----------------------------------------------------------------------------
0146 O.G = false;
0147 I.G = [ ...
0148 'Flag to include gain matrix in X (alse denoted as Dy).'
0149 ];
0150 %-----------------------------------------------------------------------------
0151 O.ga = false;
0152 I.ga = [ ...
0153 'Flag to include Marquardt-Levenberg parameter in X. Default is 0.'    
0154 ];
0155 %-----------------------------------------------------------------------------
0156 O.ga_factor_not_ok = 3;
0157 I.ga_factor_not_ok = [ ...
0158 'The factor with which the Marquardt-Levenberg factor is increased when not ',...
0159 'a lower cost value is obtained. This starts a new sub-teration. This value ',...
0160 'must be > 1.'    
0161 ];
0162 %-----------------------------------------------------------------------------
0163 O.ga_factor_ok = 2;
0164 I.ga_factor_ok = [ ...
0165 'The factor with which the Marquardt-Levenberg factor is decreased after a ',...
0166 'lower cost values has been reached. This value must be > 1.'    
0167 ];
0168 %-----------------------------------------------------------------------------
0169 O.ga_max = 100;
0170 I.ga_max = [ ...
0171 'Maximum value for gamma factor for the Marquardt-Levenberg method. The ',...
0172 'stops if this value is reached and cost value is still not decreased. ',...'
0173 'This value must be > 0.'
0174 ];
0175 %-----------------------------------------------------------------------------
0176 O.ga_start = 4;
0177 I.ga_start = [ ...
0178 'Start value for gamma factor for the Marquardt-Levenberg method. Type:',...
0179 '#   help oem',...
0180 '#for a definition of the gamma factor. This value must be >= 0.'
0181 ];
0182 %-----------------------------------------------------------------------------
0183 O.itermethod = 'ML';
0184 I.itermethod = [ ...
0185 'Iteration method. Choices are ''GN'' for Gauss-Newton and ''ML'' or',...
0186 '''LM'' for Marquardt-Levenberg.'
0187 ];
0188 %-----------------------------------------------------------------------------
0189 O.J = false;
0190 I.J = [ ...
0191 'Flag to include weighting function matrix in X.'
0192 ];
0193 %-----------------------------------------------------------------------------
0194 O.jexact = false;
0195 I.jexact = [ ...
0196 'Flag to select recalculation of J after last iteration. If not set to 1, ',...
0197 'J will correspond to x before last iteration. Also used for the linear case.'
0198 ];
0199 %-----------------------------------------------------------------------------
0200 O.jfast = false;
0201 I.jfast = [ ...
0202 'Flag to always calculate the Jacobian in parallel to the spectrum. This ',...
0203 'field is only used for the Marquardt-Levenberg case. This option can save ',...
0204 'time if the calculation of the Jacobian is very fast and the convergence ',...
0205 'is smooth (few cases where ga has to be increased). The advanatge of this ',...
0206 'option is that the next iteration can be started without a call of the ',...
0207 'forward model.'
0208 ];
0209 %-----------------------------------------------------------------------------
0210 O.linear = false;
0211 I.linear = [ ...
0212 'Flag to trigger a linear inversion. Fields like itermethod are ignored if ',...
0213 'this option is selected. Default is non-linear (0).'
0214 ];
0215 %-----------------------------------------------------------------------------
0216 O.maxiter = 99;
0217 I.maxiter = [ ...
0218 'Maximum number of iterations.'
0219 ];
0220 %-----------------------------------------------------------------------------
0221 O.msg = [];
0222 I.msg = [ ...
0223 'Message to put at the start of output messages. Can include e.g. number ',...
0224 'of retrieval case.'
0225 ];
0226 %-----------------------------------------------------------------------------
0227 O.outfids = 1;
0228 I.outfids = [ ...
0229 'File identifiers for output messages. Inlcude 1 for the screen. Set to [] ',...
0230 'for no output att all.'
0231 ];
0232 %-----------------------------------------------------------------------------
0233 O.S = false;
0234 I.S = [ ...
0235 'Flag to include covariance matrix for total error in X. That is, the sum ',...
0236 'of So and Ss.'
0237 ];
0238 %-----------------------------------------------------------------------------
0239 O.So = false;
0240 I.So = [ ...
0241 'Flag to include covariance matrix for observation error in X.'
0242 ];
0243 %-----------------------------------------------------------------------------
0244 O.Ss = false;
0245 I.Ss = [ ...
0246 'Flag to include covariance matrix for smoothing error in X.'
0247 ];
0248 %-----------------------------------------------------------------------------
0249 O.stop_dx = 1e-3;
0250 I.stop_dx = [ ...
0251 'Stop criterion. The iteration is halted when the change in x ',...
0252 'is < stop_dx (see Eq. 5.29 in Rodgers'' book). A normalisation to the ',...
0253 'length of x is applied.'
0254 ];
0255 %-----------------------------------------------------------------------------
0256 O.sxnorm = false;
0257 I.sxnorm = [ ...
0258 'Flag to internally perform a normalisation of x, based on the diagonal ',...
0259 'elements of Sx. Numerical problems can occur when the retrieved values ',...
0260 'differ strongly in magnitude (due to poor condition number for matrix ',...
0261 'inversions). This flag can be used to overcome this problem.',...
0262 '#The inverse of Sx must be calculated internally if this option is used ',...
0263 'and there is no use in pre-calculating *Sxinv*.'
0264 ];
0265 %-----------------------------------------------------------------------------
0266 O.yf = false;
0267 I.yf = [ ...
0268 'Flag to include "fitted spectrum" in X. That is, the simulated ',...
0269 'measurement matching retrieved state.'
0270 ];
0271 %-----------------------------------------------------------------------------
0272 O.Xiter = false;
0273 I.Xiter = [ ...
0274 'Flag to include all iteration states in X.'
0275 ];
0276 %-----------------------------------------------------------------------------
0277  
0278 X = O;
0279 R = I; 
0280 return  
0281 end
0282   
0283 
0284 
0285 
0286 
0287 %==============================================================================
0288 %=== Start of calculation part
0289 %==============================================================================
0290 
0291 
0292 %=== Check of input ===========================================================
0293                                                                             %&%
0294   rqre_nargin( 9, nargin );                                                 %&%
0295                                                                             %&%
0296   %- O                                                                      %&%
0297   %                                                                         %&%
0298   rqre_datatype( O, @isstruct );                                            %&%
0299   qcheck( @oem, O );                                                        %&%
0300   %                                                                         %&%
0301   rqre_datatype( O.linear, @isboolean, 'O.linear' );                        %&%
0302   rqre_datatype( O.outfids, {@isempty,@isvector}, 'O.outfids' );            %&%
0303   %                                                                         %&%
0304   if ~O.linear                                                              %&%
0305     rqre_datatype( O.itermethod, @ischar, 'O.itermethod' );                 %&%
0306     rqre_datatype( O.maxiter, @istensor0, 'O.maxiter' );                    %&%
0307     rqre_datatype( O.stop_dx, @istensor0, 'O.stop_dx' );                    %&%
0308     rqre_in_range( O.maxiter, 1, [], 'O.maxiter' );                         %&%
0309     rqre_in_range( O.stop_dx, 0, [], 'O.stop_dx' );                         %&%
0310     %                                                                       %&%
0311     if strcmp( upper(O.itermethod), 'GN' )                                  %&%
0312     elseif strcmp( upper(O.itermethod), 'ML' )  |  ...                      %&%
0313            strcmp( upper(O.itermethod), 'LM' )                              %&%
0314       %                                                                     %&%
0315       rqre_datatype( O.ga_start, @istensor0, 'O.ga_start' );                %&%
0316       rqre_datatype( O.ga_max, @istensor0, 'O.ga_max' );                    %&%
0317       rqre_datatype( O.ga_factor_ok, @istensor0, 'O.ga_factor_ok' );        %&%
0318       rqre_datatype( O.ga_factor_not_ok, @istensor0, 'O.ga_factor_not_ok' );%&%
0319       %                                                                     %&%
0320       rqre_in_range( O.ga_start, 0, [], 'O.ga_start' );                     %&%
0321       if O.ga_max <= 1                                                      %&%
0322         error( 'O.ga_max must be > 1.' );                                   %&%
0323       end                                                                   %&%
0324       if O.ga_factor_ok <= 1                                                %&%
0325         error( 'O.ga_factor_ok must be > 1.' );                             %&%
0326       end                                                                   %&%
0327       if O.ga_factor_not_ok <= 1                                            %&%
0328         error( 'O.ga_factor_not_ok must be > 1.' );                         %&%
0329       end                                                                   %&%
0330     else                                                                    %&%
0331       error( ['Unknown choice for iteration method. ', ...                  %&%
0332                              'Possible choices are ''GN'' and ''ML''.' ] ); %&%
0333     end                                                                     %&%
0334   end                                                                       %&%
0335   %                                                                         %&%
0336                                                                             %&%
0337   rqre_datatype( Q, @isstruct );                                            %&%
0338   rqre_datatype( R, {@isempty,@isstruct} );                                 %&%
0339   rqre_datatype( comfun, @isfunction_handle );                              %&%
0340                                                                             %&%
0341   rqre_datatype( Se, @istensor2 );                                          %&%
0342   rqre_datatype( Sx, @istensor2 );                                          %&%
0343   if size(Se,1) ~= size(Se,2)                                               %&%
0344     error('Input argument *Se* must be a square matrix.');                  %&%
0345   end                                                                       %&%
0346   if size(Sx,1) ~= size(Sx,2)                                               %&%
0347     error('Input argument *Sx* must be a square matrix.');                  %&%
0348   end                                                                       %&%
0349                                                                             %&%
0350   if nargin >= 7  &&  ~isempty(Sxinv)                                       %&%
0351     if size(Sxinv,1) ~= size(Sxinv,2)                                       %&%
0352       error('Input argument *Sxinv* must be a square matrix.');             %&%
0353     end                                                                     %&%
0354     if size(Sx,1) ~= size(Sxinv,1)                                          %&%
0355       error('Mismatch in size between *Sxinv* and *Sx*.');                  %&%
0356     end                                                                     %&%
0357   end                                                                       %&%
0358   %                                                                         %&%
0359   if nargin >= 8  &&  ~isempty(Seinv)                                       %&%
0360     if size(Seinv,1) ~= size(Seinv,2)                                       %&%
0361       error('Input argument *Seinv* must be a square matrix.');             %&%
0362     end                                                                     %&%
0363     if size(Se,1) ~= size(Seinv,1)                                          %&%
0364       error('Mismatch in size between *Seinv* and *Se*.');                  %&%
0365     end                                                                     %&%
0366   end                                                                       %&%
0367                                                                             %&%
0368   if nargin >= 9                                                            %&%
0369     rqre_datatype( xa, @isvector );                                         %&%
0370     if size(xa,1) ~= size(Sx,1)                                             %&%
0371       error('Mismatch in size between *Sx* and *xa*.');                     %&%
0372     end                                                                     %&%
0373   end                                                                       %&%
0374                                                                             %&%
0375   if nargin < 10                                                            %&%
0376     if O.cost                                                               %&%
0377       error('O.cost is not handled for characterisation.');                 %&%
0378     end                                                                     %&%
0379     if O.yf                                                                 %&%
0380       error('O.yf is not handled for characterisation.');                   %&%
0381     end                                                                     %&%
0382   else                                                                      %&%
0383     rqre_datatype( y, @isvector );                                          %&%
0384     if size(y,1) ~= size(Se,1)                                              %&%
0385       error('Mismatch in size between *Se* and *y*.');                      %&%
0386     end                                                                     %&%
0387   end                                                                       %&%
0388 
0389   %- Some post-processing of input
0390   %
0391   if ~O.linear
0392     if strcmp( upper(O.itermethod), 'GN' )
0393       O.itermethod = 'GN';  
0394     elseif strcmp( upper(O.itermethod), 'ML' )  |  ...
0395            strcmp( upper(O.itermethod), 'LM' )
0396         O.itermethod = 'ML';
0397         O.cost       = true;  % Note default value here, cost must
0398         %                       anyhow be calculated
0399     end
0400   end
0401   %
0402   if O.sxnorm
0403     xnorm = full( sqrt( diag( Sx ) ) );  % Avoid that xnorm inherits sparse
0404     Sx    = Sx ./ (xnorm*xnorm');        % from Sx!
0405     Sxinv = [];
0406   else
0407     xnorm = NaN;
0408   end
0409   %
0410   if nargin < 7  |  isempty(Sxinv)
0411     Sxinv = full_or_sparse( Sx \ speye(size(Sx)) );
0412   end
0413   %
0414   if nargin < 8  |  isempty(Seinv)
0415     Seinv = full_or_sparse( Se \ speye(size(Se)) );
0416   end
0417 
0418     
0419 %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
0420 
0421 
0422 
0423   
0424 %=== Stuff common to all calculation options ==================================
0425 
0426   %- Print initial screen message
0427   %
0428   if ~isempty( O.msg )  |  out(2)
0429     out( 1, 1, O.outfids );
0430   end
0431   if ~isempty( O.msg )
0432     out( 1, O.msg, O.outfids );
0433     out( 3, 0, O.outfids );
0434   end
0435 
0436   %- Make sure that X is initialised as a structure
0437   %
0438   X.converged = NaN;   % Expected output for the linear case
0439 
0440   %- xa is used as initial linearisation point
0441   %
0442   x    = xa;
0443   iter = 1;
0444   %
0445   if O.Xiter
0446     X.Xiter(:,iter) = x;
0447   end
0448   %
0449   if O.dx
0450     X.dx(iter) = NaN;
0451   end
0452 
0453   %- Is J valid for last x?
0454   %
0455   j_updated = 0;
0456 
0457   %- dx and ga not always used, but always printed:
0458   %
0459   dx = NaN;
0460   ga = 0;
0461   
0462 
0463   
0464   
0465   
0466 %=== Only characterisation ====================================================
0467 
0468   if nargin < 10
0469     %
0470     out( 1, 'Performing linear characterisation ...' );
0471     J = NaN;
0472 
0473 %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
0474 
0475   
0476   
0477   
0478   
0479 %=== Make inversion ===========================================================
0480 %
0481   else
0482     
0483     % The different cases are handled separately to avoid a lot of flags, if
0484     % statements and control variables. Helps also to save memory and to make
0485     % the calculations in most efficient way. The drawback is that more or less
0486     % the same code is found at several places.
0487     
0488     
0489     %=== Linear inversion =====================================================
0490     %
0491     if O.linear
0492 
0493       %- Messages
0494       out( 3, 'Linear inversion', O.outfids );
0495       if O.cost
0496         out( 3, 0, O.outfids );
0497         nlinprint( O.outfids );
0498       end
0499 
0500       %- Calculate Jacobian
0501       [R,yf,J] = calc_y_and_j( comfun, Q, R, O, xnorm, x, iter );
0502 
0503       %- Calculate cost
0504       if O.cost
0505         X = calc_cost_for_X( y, yf, Seinv, xa, x, Sxinv, O, xnorm, X, iter );
0506         nlinprint( O.outfids, iter, ga, X.cost(iter), X.cost_x(iter), ...
0507                                                         X.cost_y(iter), dx );
0508       end
0509       
0510       %- Some matrix products (also used in common part below)
0511       JtSeinv = J' * Seinv;
0512       SJSJ    = Sxinv + JtSeinv * J;
0513 
0514       %- Calculate x
0515       xstep   = SJSJ \ ( JtSeinv * ( y - yf ) ); 
0516       %
0517       if O.sxnorm 
0518         x = x + xnorm .* xstep;
0519       else
0520         x = x + xstep;
0521       end
0522       %
0523       iter = iter + 1;      
0524       if O.Xiter
0525         X.Xiter(:,iter) = x;
0526       end
0527       
0528       %- Flag that yf and J are not valid for x
0529       yf        = NaN;
0530       j_updated = 0;
0531 
0532       % cost for last state is handled below in a common way
0533         
0534         
0535         
0536 
0537         
0538     %=== GN iteration ========================================================
0539     %
0540     elseif strcmp( O.itermethod, 'GN' )
0541       
0542       %- Messages
0543       out( 3, 'Non-linear inversion: Gauss-Newton', O.outfids );
0544       if O.cost
0545         out( 3, 0, O.outfids );
0546         nlinprint( O.outfids );
0547       end
0548 
0549       converged = 0;
0550       %
0551       while ~converged  &  iter <= O.maxiter
0552 
0553         %- Calculate Jacobian
0554         [R,yf,J] = calc_y_and_j( comfun, Q, R, O, xnorm, x, iter );
0555 
0556         %- Calculate cost
0557         if O.cost
0558           X = calc_cost_for_X( y, yf, Seinv, xa, x, Sxinv, O, xnorm, X, iter );
0559           nlinprint( O.outfids, iter, ga, X.cost(iter), X.cost_x(iter), ...
0560                                                           X.cost_y(iter), dx );
0561         end
0562         
0563         %- Some matrix products (also used in common part below)
0564         JtSeinv = J' * Seinv;
0565         SJSJ    = Sxinv + JtSeinv * J;
0566 
0567         %- Weighted difference between y and yf
0568         if iter == 1
0569           wdy = JtSeinv * ( y - yf );
0570         else
0571           if O.sxnorm 
0572             wdy = JtSeinv * ( y - yf ) - Sxinv * ( ( x - xa ) ./ xnorm );
0573           else
0574             wdy = JtSeinv * ( y - yf ) - Sxinv * ( x - xa );
0575           end
0576         end
0577       
0578         %- Calculate new x
0579         %
0580         xold  = x;  
0581         xstep = SJSJ \ wdy;
0582         %
0583         if O.sxnorm 
0584           x = x + xnorm .* xstep;
0585         else
0586           x = x + xstep;
0587         end
0588         %
0589         iter = iter + 1;   
0590         %
0591         if O.Xiter
0592           X.Xiter(:,iter) = x;
0593         end
0594 
0595         %- Test for convergence
0596         [converged,dx,X] = test_convergence( xold, x, SJSJ, O, xnorm, X, iter );
0597         
0598       end % while ~converged
0599       
0600       %- Flag that yf and J are not valid for x
0601       yf        = NaN;
0602       j_updated = 0;
0603 
0604       % cost for last state is handled below in a common way
0605       
0606       
0607       
0608 
0609         
0610     %=== ML iteration ========================================================
0611     %
0612     elseif strcmp( O.itermethod, 'ML' )
0613       
0614       %- Messages
0615       out( 3, 'Non-linear inversion: Marquardt-Levenberg', O.outfids );
0616       if O.cost
0617         out( 3, 0, O.outfids );
0618         nlinprint( O.outfids );
0619       end
0620       
0621       %- Iteration
0622       ga_treshold = 1;
0623       ga          = O.ga_start;
0624       converged   = 0;
0625       first_iter  = 1;
0626       %
0627       while ~converged  &  iter <= O.maxiter
0628 
0629         %- Calculate Jacobian
0630         if ~j_updated
0631           [R,yf,J] = calc_y_and_j( comfun, Q, R, O, xnorm, x, iter );
0632         end
0633 
0634         %- If first iteration, cost is unknown.
0635         %- For later iterations, the gamma factor shall be decreased
0636         if first_iter
0637           [cost,cost_y,cost_x] = ...
0638                               calc_cost( y, yf, Seinv, xa, x, Sxinv, O, xnorm );
0639           %
0640           nlinprint( O.outfids, iter, NaN, cost, cost_x, cost_y, dx );
0641           if O.cost
0642             X.cost(iter)   = cost;
0643             X.cost_y(iter) = cost_y;
0644             X.cost_x(iter) = cost_x;
0645             X.ga(iter) = NaN;
0646           end
0647           first_iter = 0;
0648         else
0649           if ga >= O.ga_factor_ok * ga_treshold;
0650             ga = ga / O.ga_factor_ok;
0651           else
0652             ga = 0;
0653           end   
0654         end 
0655         
0656         %- Weighted difference between y and yf (not dependent on ga)
0657         %
0658         JtSeinv  = J' * Seinv;   % Neither dependent on ga
0659         JtSeinvJ = JtSeinv * J;   
0660         %
0661         if iter == 1
0662           wdy = JtSeinv * ( y - yf );
0663         else
0664           if O.sxnorm 
0665             wdy = JtSeinv * ( y - yf ) - Sxinv * ( ( x - xa ) ./ xnorm );
0666           else
0667             wdy = JtSeinv * ( y - yf ) - Sxinv * ( x - xa );
0668           end   
0669         end
0670         
0671         %- Find new x
0672         %
0673         xfound   = 0;
0674         cost_old = cost;
0675         %
0676         while ~xfound 
0677 
0678           %- New x state
0679           if ga > 0
0680             SJSJ = (1+ga)*Sxinv + JtSeinvJ;      
0681           else  
0682             SJSJ = Sxinv + JtSeinvJ;
0683           end   
0684           %
0685           xstep = SJSJ \ wdy; 
0686           %
0687           if O.sxnorm 
0688             xnew = x + xnorm .* xstep;
0689           else
0690             xnew = x + xstep;
0691           end
0692           
0693           %- New fitted spectrum
0694           if O.jfast
0695             [R,yf,J]  = calc_y_and_j( comfun, Q, R, O, xnorm, xnew, iter+1 );
0696             j_updated = 1; 
0697           else
0698             [R,yf] = feval( comfun, Q, R, xnew, iter+1 );
0699           end
0700 
0701           %- Check if lower cost has been reached
0702           [cost,cost_y,cost_x] = ...
0703                           calc_cost( y, yf, Seinv, xa, xnew, Sxinv, O, xnorm );
0704           %
0705           if cost < cost_old
0706             xfound = 1;
0707             iter   = iter + 1;
0708             if O.cost
0709               X.cost(iter)   = cost;
0710               X.cost_y(iter) = cost_y;
0711               X.cost_x(iter) = cost_x;
0712               X.ga(iter) = ga;       
0713             end
0714             if O.Xiter
0715               X.Xiter(:,iter) = xnew;
0716             end
0717           else
0718             nlinprint( O.outfids, -99, ga, cost, cost_x, cost_y, NaN );
0719             if ga < ga_treshold
0720               ga = ga_treshold;
0721             else
0722               if ga < O.ga_max
0723                 ga = ga * O.ga_factor_not_ok;
0724                 if ga > O.ga_max
0725                   ga = O.ga_max;
0726                 end
0727               else
0728                 X.converged = -1;
0729                 yf          = NaN;
0730                 break;    % New x could not be found !!!
0731               end
0732             end   
0733           end
0734           
0735         end %while ~xfound
0736 
0737         %- Test for convergence
0738         if xfound
0739           [converged,dx,X] = ...
0740                           test_convergence( x, xnew, SJSJ, O, xnorm, X, iter );
0741           nlinprint( O.outfids, iter, ga, cost, cost_x, cost_y, dx );
0742           x = xnew;
0743         else
0744           break;
0745         end
0746       end % while ~converged
0747       
0748       %- J can not be used for error characterisation if ga > 0
0749       %
0750       if ga > 0
0751         j_updated = 0;
0752         J         = NaN;
0753       end
0754       
0755     else
0756       error( 'Unknown retrieval option.' );    
0757     end
0758   end % of retrieval option part
0759 
0760 
0761   
0762   
0763 
0764 
0765 %=== Update/calculate basic inversion variables ===============================
0766 
0767   %- Update J?
0768   %
0769   if O.J | O.G | O.A | O.S | O.So | O.Ss | O.e | O.eo | O.es    
0770     %
0771     if isnan(J)  |  ( O.jexact  &  ~j_updated ) 
0772       [R,yf,J] = calc_y_and_j( comfun, Q, R, O, xnorm, x, iter );
0773       %
0774       JtSeinv  = J' * Seinv;
0775       SJSJ     = Sxinv + JtSeinv * J;
0776     end
0777   end
0778   
0779   %- Calculate G ?
0780   %
0781   if O.G | O.A | O.S | O.So | O.Ss | O.e | O.eo | O.es
0782     G = ( SJSJ \ speye(size(Sxinv)) ) * JtSeinv;
0783     % back-normalisation with xnorm must be done after calculation of A
0784   end
0785 
0786   %- Calculate A ?
0787   %
0788   if O.A | O.S | O.Ss | O.e | O.es
0789     A = G * J;
0790   end
0791 
0792   
0793   
0794   
0795   
0796 %=== Fill X ===================================================================
0797 
0798   X.x = x;
0799     
0800   if O.cost  &  length(X.cost) < iter
0801     if isnan(yf)
0802       [R,yf] = feval( comfun, Q, R, x, iter );
0803     end
0804     X = calc_cost_for_X( y, yf, Seinv, xa, x, Sxinv, O, xnorm, X, iter );
0805     nlinprint( O.outfids, iter, ga, X.cost(iter), X.cost_x(iter), ...
0806                                                   X.cost_y(iter), dx );
0807                                               
0808     X.ga(iter) = ga;
0809   end
0810 
0811   if O.yf
0812     if isnan(yf)
0813       [R,yf] = feval( comfun, Q, R, x, iter );
0814     end
0815     %
0816     X.yf = yf;
0817   end
0818 
0819   if O.J
0820     if O.sxnorm
0821       X.J = J ./ repmat( xnorm', size(J,1) , 1 );
0822     else
0823       X.J = J;
0824     end
0825   end
0826 
0827   if O.G
0828     if O.sxnorm
0829       X.G = G .* repmat( xnorm, 1, size(G,2) );
0830     else
0831       X.G = G;
0832     end
0833   end
0834 
0835   if O.A 
0836     if O.sxnorm
0837       X.A = A .* (xnorm*(1./xnorm'));  %scaling is ni/nj, where ni is xnorm
0838     else                               %for retrieved value, and nj for x-state
0839       X.A = A; 
0840     end
0841   end  
0842 
0843   if O.S  |  O.So  |  O.e  |  O.eo
0844     So = G * Se * G';
0845     %
0846     if O.sxnorm
0847       So = So .* (xnorm*xnorm');
0848     end
0849     %
0850     if O.So
0851       X.So = So;
0852     end
0853     %
0854     if O.eo
0855       X.eo = full( sqrt( diag( So ) ) );
0856     end
0857   end  
0858 
0859   if O.S  |  O.Ss  |  O.e  |  O.es
0860     AI = A - eye(size(A,1));
0861     Ss = AI * Sx * AI';
0862     clear AI;
0863     %
0864     if O.sxnorm
0865       Ss = Ss .* (xnorm*xnorm');
0866     end
0867     %
0868     if O.S
0869       X.S = Ss + So;
0870     end
0871     %
0872     if O.Ss
0873       X.Ss = Ss;
0874     end
0875     %
0876     if O.e
0877       X.e = full( sqrt( diag( Ss ) + diag( So ) ) );
0878     end
0879     %
0880     if O.es
0881       X.es = full( sqrt( diag( Ss ) ) );
0882     end
0883   end  
0884 
0885   if O.ex
0886     if O.sxnorm
0887       X.ex = full( sqrt( diag( Sx ) ) .* xnorm );
0888     else
0889       X.ex = full( sqrt( diag( Sx ) ) );
0890     end
0891   end
0892 %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
0893 
0894 
0895 
0896 %=== End of main part =========================================================
0897   
0898   %= End screen meassages
0899   %
0900   if ~isempty( O.msg )  |  out(2)
0901     out( 1, -1, O.outfids );
0902   end
0903   
0904   return
0905 
0906 %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
0907 %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
0908 %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
0909 %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
0910 %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
0911 
0912 
0913 
0914 
0915 
0916 %= Internal sub-functions =====================================================
0917 
0918 %=== nlinprint ================================================================
0919   %
0920   function nlinprint(outfids,iter,ga,cost,xcost,ycost,dx)
0921     %
0922     if nargin == 1
0923      out( 2,'             Gamma       Total   Profile   Spectrum    Converg.',...
0924                                                                      outfids );
0925      out( 2,'Iteration   factor        cost      cost       cost     measure',...
0926                                                                      outfids );
0927     else
0928      if iter == 1, ga = NaN; end
0929      out( 2, sprintf('%9.0f%9.1f%12.3f%10.2f%11.2f%12.2f', iter, ga, ...
0930                                            cost, xcost, ycost, dx ), outfids );
0931     end
0932   return
0933 %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
0934 
0935   
0936 %=== calc_cost ================================================================
0937   %
0938   function [cost,cost_y,cost_x] = calc_cost(y,yf,Seinv,xa,x,Sxinv,O,xnorm)
0939     dd       = x - xa;
0940     if O.sxnorm
0941       dd = dd ./ xnorm;
0942     end
0943     ny       = length( y );
0944     cost_x   = (dd' * Sxinv * dd) / ny;  
0945     dd       = y - yf;
0946     cost_y   = (dd' * Seinv * dd) / ny; 
0947     cost     = cost_x + cost_y;
0948   return
0949 %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
0950 
0951     
0952 %=== calc_cost_for_X ==========================================================
0953   %
0954   function X = calc_cost_for_X(y,yf,Seinv,xa,x,Sxinv,O,xnorm,X,iter)
0955     %
0956     [cost,cost_y,cost_x] = calc_cost( y, yf, Seinv, xa, x, Sxinv, O, xnorm );
0957     %
0958     X.cost(iter)   = cost;
0959     X.cost_y(iter) = cost_y;
0960     X.cost_x(iter) = cost_x;
0961   return
0962 %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
0963     
0964 
0965 %=== calc_y_and_j =============================================================
0966   %
0967   function [R,yf,J] = calc_y_and_j(comfun,Q,R,O,xnorm,x,iter)
0968     %
0969     [R,yf,J]  = feval( comfun, Q, R, x, iter );
0970     %
0971     if O.sxnorm
0972       J = J .* repmat( xnorm', size(J,1), 1 );
0973     end
0974   return
0975 %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
0976 
0977 
0978 %=== test_convergence =========================================================
0979   %
0980   function [converged,dx,X] = test_convergence(xold,xnew,SJSJ,O,xnorm,X,iter)
0981     %
0982     dd = xnew - xold;
0983     %
0984     if O.sxnorm 
0985       dd = dd ./ xnorm;
0986     end
0987     %
0988     dx = ( dd' * SJSJ * dd ) / length(xnew);
0989     %
0990     if dx <= O.stop_dx
0991       converged   = 1;
0992     else
0993       converged   = 0;
0994     end
0995     X.converged   = converged;
0996     if O.dx
0997       X.dx(iter)  = dx;
0998     end
0999     %
1000   return  
1001 %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Generated on Mon 15-Sep-2014 13:31:28 by m2html © 2005