calc_RtR_prior

 RtR = calc_RtR_prior( inv_model )
 CALC_RtR_PRIOR: calculate image regularization prior
   R'*R (which is an estimate of the inverse of the covariance)

   Typically, the image prior is matrix n_elem x n_elem of the
   normalized a priori crosscorrelation FEM element values
 
 calc_RtR_prior can be called as
    RtR_prior= calc_RtR_prior( inv_model, ... )

 and will call the function inv_model.RtR_prior
 parameters to RtR_prior should be passed in the field
 inv_model.RtR_prior_function_name.parameters

 If inv_model.RtR_prior is a matrix, calc_RtR_prior will return that matrix,
 possibly correcting for coarse2fine

 if there exists a field inv_model.rec_model, then
   the prior is calculated on the rec_model rather than
   the fwd_model. This will not be done if 
 inv_model.prior_use_fwd_not_rec= 1;

 RtR_prior    the calculated RtR regularization prior
 inv_model    is an inv_model structure

 If a function to calculate RtR_prior is not provided,
 RtR = R_prior' * R_prior;

0001 function RtR_prior = calc_RtR_prior( inv_model )
0002 % RtR = calc_RtR_prior( inv_model )
0003 % CALC_RtR_PRIOR: calculate image regularization prior
0004 %   R'*R (which is an estimate of the inverse of the covariance)
0005 %
0006 %   Typically, the image prior is matrix n_elem x n_elem of the
0007 %   normalized a priori crosscorrelation FEM element values
0008 %
0009 % calc_RtR_prior can be called as
0010 %    RtR_prior= calc_RtR_prior( inv_model, ... )
0011 %
0012 % and will call the function inv_model.RtR_prior
0013 % parameters to RtR_prior should be passed in the field
0014 % inv_model.RtR_prior_function_name.parameters
0015 %
0016 % If inv_model.RtR_prior is a matrix, calc_RtR_prior will return that matrix,
0017 % possibly correcting for coarse2fine
0018 %
0019 % if there exists a field inv_model.rec_model, then
0020 %   the prior is calculated on the rec_model rather than
0021 %   the fwd_model. This will not be done if
0022 % inv_model.prior_use_fwd_not_rec= 1;
0023 %
0024 % RtR_prior    the calculated RtR regularization prior
0025 % inv_model    is an inv_model structure
0026 %
0027 % If a function to calculate RtR_prior is not provided,
0028 % RtR = R_prior' * R_prior;
0029 
0030 % (C) 2005-2018 Andy Adler. License: GPL version 2 or version 3
0031 % $Id: calc_RtR_prior.m 5683 2018-02-24 21:53:16Z aadler $
0032 
0033 if ischar(inv_model) && strcmp(inv_model,'UNIT_TEST'); do_unit_test; return; end
0034 
0035 inv_model = rec_or_fwd_model( inv_model);
0036 
0037 if isfield(inv_model,'RtR_prior')
0038    if isnumeric(inv_model.RtR_prior)
0039       RtR_prior = inv_model.RtR_prior;
0040    else
0041       try inv_model.RtR_prior = str2func(inv_model.RtR_prior); end
0042       RtR_prior= feval( inv_model.RtR_prior, inv_model );
0043    end
0044    % cache results of c2f call
0045    RtR_prior= eidors_cache(@c2f_RtR_prior, ...
0046          {inv_model,RtR_prior}, 'calc_RtR_prior' );
0047 elseif isfield(inv_model,'R_prior')
0048    RtR_prior = eidors_cache(@calc_from_R_prior, inv_model, 'calc_RtR_prior');
0049 else
0050    error('calc_RtR_prior: neither R_prior nor RtR_prior provided');
0051 end
0052 
0053 function RtR_prior = c2f_RtR_prior( inv_model, RtR_prior);
0054    if isfield(inv_model.fwd_model,'coarse2fine')
0055       c2f= inv_model.fwd_model.coarse2fine;
0056       if size(RtR_prior,1)==size(c2f,1)
0057    %     we need to take into account coarse2fine - using a reasonable tol
0058          eidors_msg('calc_RtR_prior: using coarse2fine to model RtR');
0059          f2c= c2f'; %pinv(c2f,1e-6);
0060          RtR_prior=f2c*RtR_prior*c2f;
0061 
0062          % now update the object in the cache
0063       end
0064    end
0065 
0066 function RtR_prior = calc_from_R_prior(inv_model)
0067 
0068    % The user has provided an R prior. We can use this to
0069    % calculate RtR= R'*R;
0070    if isnumeric(inv_model.R_prior)
0071       R = inv_model.R_prior;
0072    else
0073       try inv_model.R_prior = str2func(inv_model.R_prior); end
0074       R= eidors_cache( inv_model.R_prior, inv_model );
0075    end
0076 
0077    RtR_prior = c2f_RtR_prior(inv_model, R'*R );
0078 
0079    
0080    
0081 
0082 function inv_model = rec_or_fwd_model( inv_model);
0083 
0084    if isfield(inv_model,'rec_model');
0085       use_rec_model = 1;
0086       try if inv_model.prior_use_fwd_not_rec== 1;
0087          use_rec_model = 0;
0088       end; end
0089 
0090       if use_rec_model
0091           % copy the normalize flag from the fwd_model to prevent warnings
0092          inv_model.rec_model = mdl_normalize(inv_model.rec_model, ...
0093                                        mdl_normalize(inv_model.fwd_model));
0094          inv_model.fwd_model= inv_model.rec_model;
0095          inv_model= rmfield(inv_model,'rec_model');
0096       end
0097    end
0098 
0099 
0100 function do_unit_test
0101    imdl = mk_common_model('a2c2',16);;
0102    imdl.RtR_prior = @prior_tikhonov;
0103    RtR = calc_RtR_prior(imdl);
0104    unit_test_cmp('RtR=prior_tik',RtR,speye(size(RtR,1)),1e-14);
0105    imdl.RtR_prior = @prior_laplace;
0106    RtR = calc_RtR_prior(imdl);
0107    unit_test_cmp('RtR=prior_lapl',RtR(1:3,1:3),[6,-2,0;-2,6,-2;0,-2,6],1e-14);
0108 
0109    imdl = rmfield(imdl,'RtR_prior');
0110    imdl.R_prior = @prior_tikhonov;
0111    RtR = calc_RtR_prior(imdl);
0112    unit_test_cmp('R=prior_tik',RtR,speye(size(RtR,1)),1e-14);
0113    imdl.R_prior = @prior_TV;
0114    RtR = calc_RtR_prior(imdl);
0115    unit_test_cmp('R=prior_TV',RtR(1:3,1:3)*16,[4,-1,0;-1,4,-1;0,-1,4],1e-14);
0116 
0117    imdl = mk_common_model('b2c2',16);
0118    grid = linspace(-1,1,16);
0119    [imdl.rec_model, imdl.fwd_model.coarse2fine] = ...
0120         mk_grid_model( imdl.fwd_model, grid, grid);
0121    RtR = calc_RtR_prior(imdl);
0122    disp('cache should only be called once');