Home > atmlab > arts > arts_oem_init.m

arts_oem_init

PURPOSE ^

ARTS_OEM_INIT Initialization of OEM inversions by ARTS

SYNOPSIS ^

function [Q,O,R,xa] = arts_oem_init(Q,O,workfolder)

DESCRIPTION ^

 ARTS_OEM_INIT   Initialization of OEM inversions by ARTS

    This function together with *arts_oem* handle bookkeeping and
    communication to perform OEM inversion with ARTS. Not all possible
    retrieval set-ups are handled. However, it should be possible to apply
    this function for most cases, while *arts_oem* is less general (where
    handling of "non-standard" sensor retrievals are most problematic).
    The function ignores any ATMDATA, and atmospheric input shall be
    provided by RAW_ATMOSPHERE and/or ready fields (such as T_FIELD and
    WIND_V_FIELD).

    The standard sequence of function calls is:

       [Qoem,O,R,xa] = arts_oem_init( Q, O, workfolder );
       Sx            = arts_sx( Qoem, R );
       Se            = covmat ...
       X             = oem( O, Qoem, R, @arts_oem, Sx, Se, [], [], xa, y );
       %
       clear Qoem

    For "non-standard" cases, you have to make a function that can replace
    *arts_oem* and give its function handle as input to *oem*.

    Note that *arts_oem_init* can modify Q in such way that it can not be
    used for further calculations. 

    The implementation strategy is to avoid unnecessery reading and
    writing of files, in order to obtain high calculation speed. The drawback
    is that more variables must be kept in memory.

    All pure forward model variables are stored in Q, while all extra data 
    needed to handle the inversions are put into R. This data can be used
    for further processing or "re-packing" of inversions results.
    These fields always exist:
       workfolder: Path to the temporary working folder.
       jq        : ARTS description of jacobian quantities.
       ji        : Start and stop index in J for each jacobian quantities.
       cfile_y   : Name of control file to generate only y
       cfile_yj  : Name of control to generate y and J
       i_rel     : Index of species retrieved as rel (not including logrel)
       p_grid    : As the ARTS WSM with same name.
       lat_grid  : As the ARTS WSM with same name.
       lon_grid  : As the ARTS WSM with same name.
       minmax    : A matrix with min and mix limits. 

    R can further contain these fields, depending on retrieval settings:
       vmr_field0      : A priori vmr field
       sensor_response : A copy of Q.SENSOR_RESPONSE 

 FORMAT   [Q,O,R,xa] = arts_oem_init(Q,O,workfolder)
        
 Out   Q            Modified Q structure
       O            Modified O structure
       R            Structure with retrieval variables. See *arts_oem*.
       xa           A priori state vector
 In    O            OEM settings.
       Q            Qarts settings.
       workfolder   Use this as "workfolder".

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

DOWNLOAD ^

arts_oem_init.m

SOURCE CODE ^

0001 % ARTS_OEM_INIT   Initialization of OEM inversions by ARTS
0002 %
0003 %    This function together with *arts_oem* handle bookkeeping and
0004 %    communication to perform OEM inversion with ARTS. Not all possible
0005 %    retrieval set-ups are handled. However, it should be possible to apply
0006 %    this function for most cases, while *arts_oem* is less general (where
0007 %    handling of "non-standard" sensor retrievals are most problematic).
0008 %    The function ignores any ATMDATA, and atmospheric input shall be
0009 %    provided by RAW_ATMOSPHERE and/or ready fields (such as T_FIELD and
0010 %    WIND_V_FIELD).
0011 %
0012 %    The standard sequence of function calls is:
0013 %
0014 %       [Qoem,O,R,xa] = arts_oem_init( Q, O, workfolder );
0015 %       Sx            = arts_sx( Qoem, R );
0016 %       Se            = covmat ...
0017 %       X             = oem( O, Qoem, R, @arts_oem, Sx, Se, [], [], xa, y );
0018 %       %
0019 %       clear Qoem
0020 %
0021 %    For "non-standard" cases, you have to make a function that can replace
0022 %    *arts_oem* and give its function handle as input to *oem*.
0023 %
0024 %    Note that *arts_oem_init* can modify Q in such way that it can not be
0025 %    used for further calculations.
0026 %
0027 %    The implementation strategy is to avoid unnecessery reading and
0028 %    writing of files, in order to obtain high calculation speed. The drawback
0029 %    is that more variables must be kept in memory.
0030 %
0031 %    All pure forward model variables are stored in Q, while all extra data
0032 %    needed to handle the inversions are put into R. This data can be used
0033 %    for further processing or "re-packing" of inversions results.
0034 %    These fields always exist:
0035 %       workfolder: Path to the temporary working folder.
0036 %       jq        : ARTS description of jacobian quantities.
0037 %       ji        : Start and stop index in J for each jacobian quantities.
0038 %       cfile_y   : Name of control file to generate only y
0039 %       cfile_yj  : Name of control to generate y and J
0040 %       i_rel     : Index of species retrieved as rel (not including logrel)
0041 %       p_grid    : As the ARTS WSM with same name.
0042 %       lat_grid  : As the ARTS WSM with same name.
0043 %       lon_grid  : As the ARTS WSM with same name.
0044 %       minmax    : A matrix with min and mix limits.
0045 %
0046 %    R can further contain these fields, depending on retrieval settings:
0047 %       vmr_field0      : A priori vmr field
0048 %       sensor_response : A copy of Q.SENSOR_RESPONSE
0049 %
0050 % FORMAT   [Q,O,R,xa] = arts_oem_init(Q,O,workfolder)
0051 %
0052 % Out   Q            Modified Q structure
0053 %       O            Modified O structure
0054 %       R            Structure with retrieval variables. See *arts_oem*.
0055 %       xa           A priori state vector
0056 % In    O            OEM settings.
0057 %       Q            Qarts settings.
0058 %       workfolder   Use this as "workfolder".
0059 
0060 % 2010-01-11   Created by Patrick Eriksson.
0061 
0062 function [Q,O,R,xa] = arts_oem_init(Q,O,workfolder)
0063 
0064 
0065 %-----------------------------------------------------------------------------
0066 %--- Checks of input and settings --------------------------------------------
0067 %-----------------------------------------------------------------------------
0068 if atmlab( 'STRICT_ASSERT' )
0069 
0070   % Input
0071   rqre_nargin( 3, nargin );
0072   rqre_datatype( Q, @isstruct ); 
0073   rqre_datatype( O, @isstruct ); 
0074   rqre_datatype( workfolder, @ischar ); 
0075 
0076   % These checks are of course not complete. They focus on demands of
0077   % arts_oem, mainly to make sure that iteration results can be mapped back
0078   % to arts correctly. Other stuff expected to be captured by qarts2cfile
0079   % and arts itself.
0080   
0081   % General settings
0082   if ~qarts_isset( Q.J_DO) | ~Q.J_DO
0083     error( 'Q.J_DO must be set to true.' );
0084   end
0085   
0086   % t_field and vmr_field
0087   if ~qarts_isset( Q.RAW_ATMOSPHERE ) | ~Q.RAW_ATMOSPHERE
0088     if isfield( Q, 'T.RETRIEVE' ) & Q.T.RETRIEVE & ~qarts_isset( Q.T_FIELD )
0089       error( ['If temperature is retrieved, T_FIELD or RAW_ATMOSPHERE ',...
0090                                                            'must be set.'] );
0091     end
0092     if isfield( Q.T, 'MINMAX' )
0093       rqre_datatype( Q.T.MINMAX, @istensor1, 'Q.T.MINMAX' );
0094       if length( Q.T.MINMAX ) > 2
0095         error( 'Max allowed length of Q.T.MINMAX is 2.' );
0096       end
0097     end
0098     %
0099     if isfield( Q, 'ABS_SPECIES.RETRIEVE' ) & ...
0100                  any([Q.ABS_SPECIES.RETRIEVE]) & ~qarts_isset( Q.VMR_FIELD )
0101       error( ['If any absorption species is retrieved, VMR_FIELD or ',...
0102                                           'RAW_ATMOSPHERE must be set.'] );
0103     end
0104     for i = 1 : length(Q.ABS_SPECIES)
0105       if isfield( Q.ABS_SPECIES(i), 'MINMAX' )
0106         rqre_datatype( Q.ABS_SPECIES(i).MINMAX, @istensor1, ...
0107                                      sprintf('Q.ABS_SPECIES(%d).MINMAX',i) );
0108         if length( Q.ABS_SPECIES(i).MINMAX ) > 2
0109           error( 'Max allowed length of Q.ABS_SPECIES(%d).MINMAX is 2.', i );
0110         end
0111       end 
0112     end 
0113   end  
0114 
0115   % vind fields
0116   if isfield( Q, 'WIND_U.RETRIEVE' ) & Q.WIND_U.RETRIEVE & ...
0117                                               ~qarts_isset( Q.WIND_U_FIELD )
0118     error( 'If wind u-component is retrieved, WIND_U_FIELD must be set.' );
0119   end  
0120   if isfield( Q, 'WIND_V.RETRIEVE' ) & Q.WIND_V.RETRIEVE & ...
0121                                               ~qarts_isset( Q.WIND_V_FIELD )
0122     error( 'If wind v-component is retrieved, WIND_V_FIELD must be set.' );
0123   end  
0124   if isfield( Q, 'WIND_W.RETRIEVE' ) & Q.WIND_W.RETRIEVE & ...
0125                                               ~qarts_isset( Q.WIND_W_FIELD )
0126     error( 'If wind w-component is retrieved, WIND_W_FIELD must be set.' );
0127   end  
0128   
0129   % Sensor stuff
0130   if qarts_isset(Q.POINTING) & Q.POINTING.RETRIEVE
0131     if ~qarts_isset( Q.SENSOR_LOS )
0132       error( ['When any pointing retrieval is performed, Q.SENSOR_LOS ',...
0133               'must be set.'] );
0134     end
0135     if ~qarts_isset( Q.SENSOR_TIME )
0136       error( ['When any pointing retrieval is performed, Q.SENSOR_TIME ',...
0137               'must be set.'] );
0138     end
0139   end
0140   %
0141   if ( qarts_isset(Q.FSHIFTFIT)   & Q.FSHIFTFIT.RETRIEVE )  | ...
0142      ( qarts_isset(Q.FSTRETCHFIT) & Q.FSTRETCHFIT.RETRIEVE )
0143     if ~isstruct( Q.SENSOR_RESPONSE )
0144       error( ['When any frequency fit is performed Q.SENSOR_RESPONSE ',...
0145               'must follow qartsSensor (i.e. be a structure).'] );
0146     end
0147   end
0148 end
0149 %-----------------------------------------------------------------------------
0150 %--- End checks --------------------------------------------------------------
0151 %-----------------------------------------------------------------------------
0152 
0153 
0154 
0155 %- Initialization of variables
0156 %
0157 R            = [];
0158 R.i_rel      = [];
0159 R.workfolder = workfolder;
0160 %
0161 [Q,R,t_field,z_field,vmr_field,wind_u,wind_v,wind_w] = init_local( Q, R );
0162 
0163 
0164 
0165 %- Some internal variables
0166 %
0167 nq          = length( R.jq );
0168 nx          = R.ji{nq}{2};
0169 i_asj       = find( [ Q.ABS_SPECIES.RETRIEVE ] );
0170 store_vmr0  = false;
0171 
0172 
0173 %- xa
0174 %
0175 xa       = zeros( nx, 1 );
0176 
0177 
0178 %- R.minmax
0179 %
0180 R.minmax = repmat( [-Inf Inf], nx, 1 );
0181 
0182 
0183 
0184 %--- Loop retrieval quantities and fill xa and R fields ----------------------
0185 %-----------------------------------------------------------------------------
0186 
0187 for i = 1 : nq
0188 
0189   ind = R.ji{i}{1} : R.ji{i}{2};
0190 
0191   switch R.jq{i}.maintag
0192 
0193 
0194   case 'Absorption species'   %-----------------------------------------------
0195     %
0196     store_vmr0 = true;
0197     ig         = i_asj(i);    % Gas species index
0198     %
0199     if strcmp( R.jq{i}.mode, 'rel' )
0200       xa(ind)    = 1;
0201       R.i_rel    = [ R.i_rel ind ];
0202     elseif strcmp( R.jq{i}.mode, 'vmr' )
0203       xa(ind)    = mat2col( arts_regrid( Q.ATMOSPHERE_DIM, Q, ...
0204                           get_vmrfield_local(vmr_field,ig), R.jq{i}.grids ) );
0205     elseif strcmp( R.jq{i}.mode, 'nd' )
0206       xa(ind)    = vmr2nd( ...
0207                      mat2col( arts_regrid( Q.ATMOSPHERE_DIM, Q, ...   % VMR
0208                        get_vmrfield_local(vmr_field,ig), R.jq{i}.grids ) ), ...
0209                      repmat( R.jq{i}.grids{1}, ...                    % P
0210                        length(ind)/length(R.jq{i}.grids{1}), 1 ), ...
0211                      mat2col( arts_regrid( Q.ATMOSPHERE_DIM, Q, ...   % T
0212                        t_field, R.jq{i}.grids ) ) ...
0213                    );
0214     elseif strcmp( R.jq{i}.mode, 'logrel' )
0215       xa(ind)    = 0;
0216     else
0217       error( sprintf('Unknown gas species retrieval unit (%s).', ...
0218                                                               R.jq{i}.mode ) ); 
0219     end
0220     if isfield( Q.ABS_SPECIES(ig), 'MINMAX' )  &  ...
0221                                              ~isempty(Q.ABS_SPECIES(ig).MINMAX)
0222       
0223       R.minmax(ind,1:length(Q.ABS_SPECIES(ig).MINMAX)) = ...
0224                            repmat( Q.ABS_SPECIES(ig).MINMAX', length(ind), 1 );
0225     end
0226     
0227   case 'Atmospheric temperatures'   %-----------------------------------------
0228     %
0229     xa(ind) = mat2col( arts_regrid( Q.ATMOSPHERE_DIM, Q, t_field, ...
0230                                                              R.jq{i}.grids ) );
0231     if isfield( Q.T, 'MINMAX' )  & ~isempty(Q.T.MINMAX)
0232       R.minmax(ind,1:length(Q.T.MINMAX)) = ...
0233                                          repmat( Q.T.MINMAX', length(ind), 1 );
0234     end
0235 
0236   case 'Wind'   %-------------------------------------------------------------
0237     %
0238     if strcmp( R.jq{i}.subtag, 'u' )
0239       Q = force_file_local( Q, 'WIND_U_FIELD', 'Tensor3', R );
0240       if isempty( wind_u )
0241         wind_u = zeros( size( t_field ) );
0242       end
0243       xa(ind) = mat2col( arts_regrid( Q.ATMOSPHERE_DIM, Q, wind_u, ...
0244                                                              R.jq{i}.grids ) );
0245     elseif strcmp( R.jq{i}.subtag, 'v' )
0246       Q = force_file_local( Q, 'WIND_V_FIELD', 'Tensor3', R );
0247       if isempty( wind_v )
0248         wind_v = zeros( size( t_field ) );
0249       end
0250       xa(ind) = mat2col( arts_regrid( Q.ATMOSPHERE_DIM, Q, wind_v, ...
0251                                                              R.jq{i}.grids ) );
0252     elseif strcmp( R.jq{i}.subtag, 'w' )
0253       Q = force_file_local( Q, 'WIND_W_FIELD', 'Tensor3', R );
0254       if isempty( wind_w )
0255         wind_w = zeros( size( t_field ) );
0256       end
0257       xa(ind) = mat2col( arts_regrid( Q.ATMOSPHERE_DIM, Q, wind_w, ...
0258                                                              R.jq{i}.grids ) );
0259     else
0260       error( 'Unknown wind subtag.' );
0261     end
0262     
0263    %---------------------------------------------------------------------------
0264    case { 'Frequency', 'Sensor pointing', 'Polynomial baseline fit', ...
0265                                           'Sinusoidal baseline fit' }
0266     %
0267     xa(ind) = 0;    
0268 
0269   otherwise   %---------------------------------------------------------------
0270     error('Unknown retrieval quantitity.'); 
0271   end
0272 end 
0273 
0274 
0275 
0276 %--- More to put into R ? -----------------------------------------------------
0277 
0278 if store_vmr0
0279   R.vmr_field0 = vmr_field;
0280 end
0281 
0282 
0283 
0284 %--- Create control files -----------------------------------------------------
0285 %
0286 %- Only y
0287 Q.J_DO    = false;
0288 parts     = qarts2cfile( 'y_after_init' );
0289 S         = qarts2cfile( Q, parts, R.workfolder );
0290 cfile     = fullfile( R.workfolder, 'cfile_y.arts' );
0291 strs2file( cfile, S );
0292 R.cfile_y = cfile;
0293 %
0294 %- y+j
0295 Q.J_DO     = true;
0296 S          = qarts2cfile( Q, parts, R.workfolder );
0297 cfile      = fullfile( R.workfolder, 'cfile_yj.arts' );
0298 strs2file( cfile, S );
0299 R.cfile_yj = cfile;
0300 
0301 return
0302 
0303 
0304 
0305 
0306 
0307 %-----------------------------------------------------------------------------
0308 %--- Internal sub-functions
0309 %-----------------------------------------------------------------------------
0310 
0311 function Q = force_file_local( Q, fieldname, datatype, R )
0312   if ~ischar( getfield( Q, fieldname ) )
0313     filename = fullfile( R.workfolder, [lower(fieldname),'.xml'] );
0314     xmlStore( filename, Q.(fieldname), datatype );
0315     Q.(fieldname) = filename;
0316   end
0317 return
0318 
0319 
0320 
0321 function V = get_vmrfield_local( vmr_field, i )
0322   s               = ones( 1, 4 );
0323   s2              = size( vmr_field );
0324   s(1:length(s2)) = s2;
0325   V = zeros( s(2:end) );
0326   V(:,:,:) = vmr_field(i,:,:,:);
0327 return
0328 
0329 
0330 
0331 function [Q,R,t_field,z_field,vmr_field,wind_u,wind_v, wind_w] = ...
0332                                                              init_local( Q, R )
0333 
0334   %- Copy some data to R
0335   %
0336   % This is done here to possible save time; the data are forced to files below
0337   %
0338   R.p_grid   = qarts_get( Q.P_GRID );
0339   R.lat_grid = qarts_get( Q.LAT_GRID );
0340   R.lon_grid = qarts_get( Q.LON_GRID );
0341 
0342   %- Force some ARTS variables to be files
0343   %
0344   % This is done to avoid creating the same file repeatedly in qarts2cfile.
0345   % This can be done for variables that never are part of the retrieval
0346   % quantities. Some variables, beside the ones included here, are handled
0347   % below.
0348   Q = force_file_local( Q, 'F_GRID', 'Vector', R );
0349   %
0350   Q = force_file_local( Q, 'P_GRID', 'Vector', R );
0351   if Q.ATMOSPHERE_DIM > 1
0352     Q = force_file_local( Q, 'LAT_GRID', 'Vector', R );
0353     if Q.ATMOSPHERE_DIM > 2
0354       Q = force_file_local( Q, 'LON_GRID', 'Vector', R );
0355     end
0356   end
0357   %
0358   Q = force_file_local( Q, 'SENSOR_POS', 'Matrix', R );
0359   Q = force_file_local( Q, 'SENSOR_LOS', 'Matrix', R );
0360 
0361   %- Atmospheric fields
0362   %
0363   if qarts_isset( Q.RAW_ATMOSPHERE )  &  Q.RAW_ATMOSPHERE
0364     [t_field,z_field,vmr_field]       = arts_atmfields( Q );
0365     [Q.T_FIELD,Q.Z_FIELD,Q.VMR_FIELD] = deal( t_field, z_field, vmr_field );   
0366     Q.RAW_ATMOSPHERE                  = {};
0367   else
0368     t_field   = qarts_get( Q.T_FIELD );
0369     z_field   = qarts_get( Q.Z_FIELD );
0370     vmr_field = qarts_get( Q.VMR_FIELD );
0371   end
0372   %
0373   Q = force_file_local( Q, 'T_FIELD', 'Tensor3', R );
0374   Q = force_file_local( Q, 'Z_FIELD', 'Tensor3', R );
0375   Q = force_file_local( Q, 'VMR_FIELD', 'Tensor4', R );
0376   %
0377   % Winds (some stuff done, further above if wind is retrieved):
0378   if qarts_isset( Q.WIND_U_FIELD )
0379     wind_u = qarts_get( Q.WIND_U_FIELD );
0380   else
0381     wind_u = NaN;
0382   end
0383   if qarts_isset( Q.WIND_V_FIELD )
0384     wind_v = qarts_get( Q.WIND_V_FIELD );
0385   else
0386     wind_v = NaN;
0387   end
0388   if qarts_isset( Q.WIND_W_FIELD )
0389     wind_w = qarts_get( Q.WIND_W_FIELD );
0390   else
0391     wind_w = NaN;
0392   end
0393   
0394   %- Absorption
0395   %
0396   if strcmp( Q.ABSORPTION, 'LoadTable' )
0397     Q = force_file_local( Q, 'ABS_LOOKUP', 'GasAbsLookup', R );
0398   elseif strcmp( Q.ABSORPTION, 'CalcTable' )
0399     Q            = qarts_abstable( Q );
0400     Q.ABS_LOOKUP = arts_abstable( Q, R.workfolder );
0401     Q.ABSORPTION = 'LoadTable';
0402   end
0403   
0404   %- Run arts to create jacobian and sensor variables
0405   %
0406   parts = qarts2cfile( 'y_init' );
0407   %
0408   S     = qarts2cfile( Q, parts, R.workfolder );
0409   cfile = fullfile( R.workfolder, 'cfile.arts' );
0410   strs2file( cfile, S );
0411   arts( cfile );
0412 
0413   %- Jacobian
0414   %
0415   R.jq = xmlLoad( fullfile( R.workfolder, 'jacobian_quantities.xml' ) );
0416   R.ji = xmlLoad( fullfile( R.workfolder, 'jacobian_indices.xml' ) );
0417   %
0418   if length( R.jq ) == 0
0419     error( 'No retrieval quantities defined!?' );
0420   end
0421   %
0422   for i = 1 : length(R.ji)           % Fix different indexing
0423     for j = 1 : length(R.ji{i})
0424       R.ji{i}{j} = R.ji{i}{j} + 1;
0425     end
0426   end
0427 
0428   %- Sensor
0429   %
0430   if qarts_isset(Q.POINTING) & Q.POINTING.RETRIEVE
0431     R.los = qarts_get( Q.SENSOR_LOS );
0432   end
0433   %
0434   if ( qarts_isset(Q.FSHIFTFIT)   & Q.FSHIFTFIT.RETRIEVE )  | ...
0435      ( qarts_isset(Q.FSTRETCHFIT) & Q.FSTRETCHFIT.RETRIEVE )
0436     %- A copy of Q.SENSOR_RESPONSE is needed here
0437     R.sensor_response = Q.SENSOR_RESPONSE;
0438   end
0439   %
0440   Q.SENSOR_RESPONSE     = fullfile( R.workfolder, 'sensor_response.xml' );
0441   %
0442   Q.SENSOR_RESPONSE     = fullfile( R.workfolder, 'sensor_response.xml' );
0443   Q.SENSOR_RESPONSE_F   = fullfile( R.workfolder, 'sensor_response_f.xml' );
0444   Q.SENSOR_RESPONSE_ZA  = fullfile( R.workfolder, 'sensor_response_za.xml' );
0445   Q.SENSOR_RESPONSE_AA  = fullfile( R.workfolder, 'sensor_response_aa.xml' );
0446   Q.SENSOR_RESPONSE_POL = fullfile( R.workfolder, 'sensor_response_pol.xml' );
0447   Q.SENSOR_RESPONSE_F_GRID = ...
0448                         fullfile( R.workfolder, 'sensor_response_f_grid.xml' );
0449   Q.SENSOR_RESPONSE_ZA_GRID = ...
0450                        fullfile( R.workfolder, 'sensor_response_za_grid.xml' );
0451   Q.SENSOR_RESPONSE_AA_GRID  = ...
0452                        fullfile( R.workfolder, 'sensor_response_aa_grid.xml' );
0453   Q.SENSOR_RESPONSE_POL_GRID = ...
0454                       fullfile( R.workfolder, 'sensor_response_pol_grid.xml' );
0455   Q.MBLOCK_ZA_GRID      = fullfile( R.workfolder, 'mblock_za_grid.xml' );
0456   Q.MBLOCK_AA_GRID      = fullfile( R.workfolder, 'mblock_aa_grid.xml' );
0457 return

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