Home > atmlab > covmat > covmat3d_from_cfun.m

covmat3d_from_cfun

PURPOSE ^

COVMAT3D_FROM_CFUN Correlation function based covariance matrix for 1D-3D

SYNOPSIS ^

function S = covmat3d_from_cfun(dim,D,grid1,grid2,grid3)

DESCRIPTION ^

 COVMAT3D_FROM_CFUN   Correlation function based covariance matrix for 1D-3D

   Converts specification data to a covariance matrix. The covariance is
   described by the structure *D*. The function is more general than
   *covmat_partstat_corr* (for dimensions <= 3D), but is (much) slower,
   particularly for higher data dimensions or long data grids. 

   Correlation up to 3 dimensions are handled. The dimensions are here 
   denoted as 1, 2 and 3, in order to be general. The physical meaning
   of these dimensions differ. For e.g. gas species they correspond to
   pressure/altitude, latitude and longitude. 

   The data are assumed to be stored with dimension 1 as innermost "loop",
   dimension 2 next etc. For 3D case, the data order is as follows:
 
       [(x1,y1,z1) (x2,y1,z2) ...(xn,y1,z1) (x1,y2,z) ... (x1,y1,z2) ...]'

   The covariance is here specified in a parameterised way. Values are
   calculated as S(i,j) = si(i)*si(j)*c1(i,j)*c2(i,j)*c3(i,j),
   where si(i) is the standard deviation at point i, c1(i,j) is
   the dim1-correlation between point i and j, etc. 
   Correlation structures are seperated in each dimension (c1, c2 and c3).
   This means that only the distance in dimension 1 is considered when
   calculating c1 etc. Correlation lengths can then specified seperately
   for each dimension. The correlation for dimensions not included is set
   to 1.

   Correlation structures are specified by a function and correlation 
   lengths. The correlation length is the distance where the correlation
   has dropped to exp(-1). These lengths are averaged between the two
   involved positions. For example, the exponential correlation function
   is 
      c(i,j) = exp(-abs(x(i)-x(j))/((cl(i)+cl(j))/2))
   where x and cl is the position and correlation length, respectively.
   Remaining correlation functions are implemented in corresponding 
   manner.
   Setting SI, CL1, CL2 or CL3 to a scalar is shorthand for a constant
   standard deviation/correlation length. No grids need to be specified
   in this case. Otherwise, SI, CL1, CL2 and CL3 are interpolated
   linearly to selected retrieval positions (no extrapolation).
   No checks of size are done. Inconsistency in size between variables 
   will result in error from interpolation functions.

   Recognised fields of D are:
   ---
      SI         A priori standard deviation. Dimensionality must match *dim*.
                 Dimension order (dim1,dim2,dim3). 
      SI_GRID1   Grid of SI in dimension 1. 
      SI_GRID2   Grid of SI in dimension 1. Not needed if *dim* < 2.
      SI_GRID3   Grid of SI in dimension 1. Not needed if *dim* < 3.
      CCO        Correlation cut-off value. All correlation values, c,
                 -CCO < x < CCO are set to zero. This in order to avoid
                 non-significant correlation and make the covariance matrix
                 as sparse as possible.
      CFUN1      Correlation function in dimension 1. The options exist:
                  'drc'  Dirac. No correlation. The correlation
                         length is of no importance here.
                         Note that this applies only in the considered
                         dimension.
                  'lin'  Linearly decreasing to 0 (tenth function).
                  'exp'  Exponential.
                  'gau'  Gaussian. 
      CL1        Correlation length along dimension 1. Size as for SI.
      CL1_GRID1  As corresponding grid for SI.
      CL1_GRID2  As corresponding grid for SI.
      CL1_GRID3  As corresponding grid for SI.
      CFUN2      All fields below as for corresponding field for dim 1.
      CL2        
      CL2_GRID1  
      CL2_GRID2  
      CL2_GRID3  
      CFUN3
      CL3        
      CL3_GRID1  
      CL3_GRID2  
      CL3_GRID3  

 FORMAT   S = covmat3d_from_cfun(dim,D,grid1[,grid2,grid3])

 OUT   S       Covariance matrix.
 IN    dim     Dimensionality.
       D       Specification data.
       grid1   Retrieval grid in dimension 1.
 OPT   grid2   Retrieval grid in dimension 2. Not needed if *dim* < 2.
       grid3   Retrieval grid in dimension 3. Not needed if *dim* < 3.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

DOWNLOAD ^

covmat3d_from_cfun.m

SOURCE CODE ^

0001 % COVMAT3D_FROM_CFUN   Correlation function based covariance matrix for 1D-3D
0002 %
0003 %   Converts specification data to a covariance matrix. The covariance is
0004 %   described by the structure *D*. The function is more general than
0005 %   *covmat_partstat_corr* (for dimensions <= 3D), but is (much) slower,
0006 %   particularly for higher data dimensions or long data grids.
0007 %
0008 %   Correlation up to 3 dimensions are handled. The dimensions are here
0009 %   denoted as 1, 2 and 3, in order to be general. The physical meaning
0010 %   of these dimensions differ. For e.g. gas species they correspond to
0011 %   pressure/altitude, latitude and longitude.
0012 %
0013 %   The data are assumed to be stored with dimension 1 as innermost "loop",
0014 %   dimension 2 next etc. For 3D case, the data order is as follows:
0015 %
0016 %       [(x1,y1,z1) (x2,y1,z2) ...(xn,y1,z1) (x1,y2,z) ... (x1,y1,z2) ...]'
0017 %
0018 %   The covariance is here specified in a parameterised way. Values are
0019 %   calculated as S(i,j) = si(i)*si(j)*c1(i,j)*c2(i,j)*c3(i,j),
0020 %   where si(i) is the standard deviation at point i, c1(i,j) is
0021 %   the dim1-correlation between point i and j, etc.
0022 %   Correlation structures are seperated in each dimension (c1, c2 and c3).
0023 %   This means that only the distance in dimension 1 is considered when
0024 %   calculating c1 etc. Correlation lengths can then specified seperately
0025 %   for each dimension. The correlation for dimensions not included is set
0026 %   to 1.
0027 %
0028 %   Correlation structures are specified by a function and correlation
0029 %   lengths. The correlation length is the distance where the correlation
0030 %   has dropped to exp(-1). These lengths are averaged between the two
0031 %   involved positions. For example, the exponential correlation function
0032 %   is
0033 %      c(i,j) = exp(-abs(x(i)-x(j))/((cl(i)+cl(j))/2))
0034 %   where x and cl is the position and correlation length, respectively.
0035 %   Remaining correlation functions are implemented in corresponding
0036 %   manner.
0037 %   Setting SI, CL1, CL2 or CL3 to a scalar is shorthand for a constant
0038 %   standard deviation/correlation length. No grids need to be specified
0039 %   in this case. Otherwise, SI, CL1, CL2 and CL3 are interpolated
0040 %   linearly to selected retrieval positions (no extrapolation).
0041 %   No checks of size are done. Inconsistency in size between variables
0042 %   will result in error from interpolation functions.
0043 %
0044 %   Recognised fields of D are:
0045 %   ---
0046 %      SI         A priori standard deviation. Dimensionality must match *dim*.
0047 %                 Dimension order (dim1,dim2,dim3).
0048 %      SI_GRID1   Grid of SI in dimension 1.
0049 %      SI_GRID2   Grid of SI in dimension 1. Not needed if *dim* < 2.
0050 %      SI_GRID3   Grid of SI in dimension 1. Not needed if *dim* < 3.
0051 %      CCO        Correlation cut-off value. All correlation values, c,
0052 %                 -CCO < x < CCO are set to zero. This in order to avoid
0053 %                 non-significant correlation and make the covariance matrix
0054 %                 as sparse as possible.
0055 %      CFUN1      Correlation function in dimension 1. The options exist:
0056 %                  'drc'  Dirac. No correlation. The correlation
0057 %                         length is of no importance here.
0058 %                         Note that this applies only in the considered
0059 %                         dimension.
0060 %                  'lin'  Linearly decreasing to 0 (tenth function).
0061 %                  'exp'  Exponential.
0062 %                  'gau'  Gaussian.
0063 %      CL1        Correlation length along dimension 1. Size as for SI.
0064 %      CL1_GRID1  As corresponding grid for SI.
0065 %      CL1_GRID2  As corresponding grid for SI.
0066 %      CL1_GRID3  As corresponding grid for SI.
0067 %      CFUN2      All fields below as for corresponding field for dim 1.
0068 %      CL2
0069 %      CL2_GRID1
0070 %      CL2_GRID2
0071 %      CL2_GRID3
0072 %      CFUN3
0073 %      CL3
0074 %      CL3_GRID1
0075 %      CL3_GRID2
0076 %      CL3_GRID3
0077 %
0078 % FORMAT   S = covmat3d_from_cfun(dim,D,grid1[,grid2,grid3])
0079 %
0080 % OUT   S       Covariance matrix.
0081 % IN    dim     Dimensionality.
0082 %       D       Specification data.
0083 %       grid1   Retrieval grid in dimension 1.
0084 % OPT   grid2   Retrieval grid in dimension 2. Not needed if *dim* < 2.
0085 %       grid3   Retrieval grid in dimension 3. Not needed if *dim* < 3.
0086 
0087 % 2006-08-21   Created by Patrick Eriksson.
0088 
0089 
0090 function S = covmat3d_from_cfun(dim,D,grid1,grid2,grid3)
0091                                                                           %&%
0092                                                                           %&%
0093 %--- Simple checks of input                                               %&%
0094 %                                                                         %&%
0095 rqre_nargin( 3, nargin );                                                 %&%
0096 rqre_alltypes( dim, {@istensor0,@iswhole} );                              %&%
0097 rqre_in_range( dim, 1, 3 );                                               %&%
0098 rqre_datatype( D, @isstruct );                                            %&%
0099 %                                                                         %&%
0100 if nargin ~= 2+dim                                                        %&%
0101   error( ...                                                              %&%
0102   'Dimensionality (*dim*) and number of input arguments do not match.');  %&%
0103 end                                                                       %&%
0104 
0105 
0106 %--- Determine "x, y and z" position for each point
0107 %
0108 x  = vec2col(grid1);
0109 %
0110 if dim >= 2
0111   x = repmat( x, length(grid2), 1 );
0112   y = repmat( vec2row(grid2), length(grid1), 1 );
0113   y = y(:);
0114   if dim == 3
0115     x = repmat( x, length(grid3) );
0116     y = repmat( y, length(grid3) );
0117     z = repmat( vec2row(grid3), length(grid1)*length(grid2), 1 );
0118     z = z(:);
0119   end
0120 end
0121 %
0122 n = length(x);
0123 
0124 
0125 
0126 %--- Correlation part
0127 %
0128 C = 1;
0129 %
0130 %- x-dimension
0131 if strcmp( D.CFUN1, 'drc' )
0132   cl = NaN;
0133 else 
0134   if isscalar( D.CL1 )
0135     cl = repmat( D.CL1, n, 1 );
0136   else
0137     if dim == 1
0138       cl = interpd( x, D.CL1, D.CL1_GRID1, 'linear' );
0139     elseif dim == 2
0140       cl = interpd( [x,y], D.CL1, D.CL1_GRID1, D.CL1_GRID2, 'linear' );
0141     else
0142       cl = interpd( [x,y,z], D.CL1, D.CL1_GRID1, D.CL1_GRID2, D.CL1_GRID3,...
0143                                                                   'linear' );
0144     end
0145     if any( isnan(cl) )
0146       error( ['NaN obtained for dim1 correlation length. ',...
0147               'Probably caused by extrapolation at re-gridding.'] )
0148     end
0149   end
0150 end
0151 C = C .* make_cmatrix( x, cl, D.CFUN1 );
0152 %
0153 %- y-dimension
0154 if dim >= 2
0155   if strcmp( D.CFUN2, 'drc' )
0156     cl = NaN;
0157   else 
0158     if isscalar( D.CL2 )
0159       cl = repmat( D.CL2, n, 1 );
0160     else
0161       if dim == 2
0162         cl = interpd( [x,y], D.CL2, D.CL2_GRID1, D.CL2_GRID2, 'linear' );
0163       else
0164         cl = interpd( [x,y,z], D.CL2, D.CL2_GRID1, D.CL2_GRID2, ...
0165                                                      D.CL2_GRID3, 'linear' );
0166       end
0167       if any( isnan(cl) )
0168         error( ['NaN obtained for dim2 correlation length. ',...
0169                 'Probably caused by extrapolation at re-gridding.'] )
0170       end
0171     end
0172   end
0173   C = C .* make_cmatrix( y, cl, D.CFUN2 );
0174 end
0175 %
0176 %- z-dimension
0177 if dim == 3
0178   if strcmp( D.CFUN3, 'drc' )
0179     cl = NaN;
0180   else 
0181     if isscalar( D.CL3 )
0182       cl = repmat( D.CL3, n, 1 );
0183     else
0184       cl = interpd( [x,y,z], D.CL3, D.CL3_GRID1, D.CL3_GRID2, ...
0185                                                      D.CL3_GRID3, 'linear' );
0186     end
0187     if any( isnan(cl) )
0188       error( ['NaN obtained for dim3 correlation length. ',...
0189               'Probably caused by extrapolation at re-gridding.'] )
0190     end
0191   end
0192   C = C .* make_cmatrix( z, cl, D.CFUN3 );
0193 end
0194 
0195 
0196 %--- Apply correlation cut-off
0197 %
0198 if D.CCO > 0
0199   ind = find( C  &  abs(C) < D.CCO );
0200   %
0201   if isempty( ind )
0202     C(ind) = 0;
0203   end
0204 end
0205 
0206 
0207 %--- Include standard deviation
0208 %
0209 if isscalar( D.SI )
0210   si = repmat( D.SI, n, 1 );
0211 else
0212   if dim == 1
0213     si = interpd( x, D.SI, D.SI_GRID1, 'linear' );
0214   elseif dim == 2
0215     si = interpd( [x,y], D.SI, D.SI_GRID1, D.SI_GRID2, 'linear' );
0216   else
0217     si = interpd( [x,y,z], D.SI, D.SI_GRID1, D.SI_GRID2, D.SI_GRID3, ...
0218                                                                   'linear' );
0219   end
0220 end
0221 %
0222 if any( isnan(si) )
0223   error( ['NaN obtained for standard deviation. ',...
0224           'Probably caused by extrapolation at re-gridding.'] )
0225 end
0226 %
0227 S = (si * si') .* C;
0228 
0229 
0230 return
0231 
0232 %------------------------------------------------------------------------------
0233 
0234 
0235 function C = make_cmatrix( x, cl, cfun )
0236   %
0237   n = length( x );
0238   C = speye( n, n );
0239   %
0240   if strcmp( cfun, 'drc' )
0241     for row=1:(n-1)
0242       for col=(row+1):n
0243         if x(row) == x(col)
0244           C(row,col) = 1; 
0245           C(col,row) = 1;
0246         end
0247       end
0248     end
0249 
0250   elseif strcmp( cfun, 'lin' )
0251     for row=1:(n-1)
0252       for col=(row+1):n
0253         c = 1 - 2 * (1-exp(-1)) * abs( (x(row)-x(col))/(cl(row)+cl(col)) );
0254         if c>0 
0255           C(row,col) = c; 
0256           C(col,row) = c;
0257         end
0258       end
0259     end
0260 
0261   elseif strcmp( cfun, 'exp' )
0262     for row=1:(n-1)
0263       for col=(row+1):n
0264         c = exp(-abs(x(row)-x(col))/((cl(row)+cl(col))/2));
0265         C(row,col) = c; 
0266         C(col,row) = c;
0267       end
0268     end
0269 
0270   elseif strcmp( cfun, 'gau' )
0271     for row=1:(n-1)
0272       for col=(row+1):n
0273         c = exp( -4 * ( (x(row)-x(col))/(cl(row)+cl(col)) )^2  );
0274         C(row,col) = c; 
0275         C(col,row) = c;
0276       end
0277     end
0278 
0279   else
0280      error(sprintf('Unknown selection for correlation function (%s).',cfun));
0281   end
0282   %
0283 return

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