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 7097 2024-12-23 16:17:37Z 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');
0123 
0124    calc_model_prior_test;
0125 
0126 
0127 
0128 function calc_model_prior_test;
0129 % Verify model prior calcs
0130 % Imported from: calc_model_prior_test.m 3127 2012-06-08 16:19:25Z bgrychtol
0131 
0132 % TODO: also test the various inverse prior calls
0133 
0134    imdl= mk_common_model('c2c2',16);
0135    try; imdl= rmfield(imdl,'RtR_prior'); end
0136    try; imdl= rmfield(imdl,'R_prior');   end
0137 
0138    any_priors= {@prior_tikhonov, ...
0139                 @prior_noser, ...
0140                 @prior_gaussian_HPF, ...
0141                 @prior_laplace};
0142 
0143   test_R=[ 0.041630544712308
0144            6.773110974081427e-04
0145            0.036206173920553
0146            2.592951965329636
0147            0.002008656225401];
0148 
0149    R_priors=   {any_priors{:}, ...
0150                 @prior_TV};
0151 
0152    % Call R_priors as R_priors
0153    eidors_cache clear
0154    for i = 1:length(R_priors); p = R_priors{i};
0155       inv_mdl= imdl;
0156       inv_mdl.R_prior= p;
0157       R= calc_R_prior(inv_mdl);
0158 %     test = std(R'*R,0,'all');
0159       test = std(reshape(R'*R,[],1),0); % compatible with R2017b
0160       unit_test_cmp(['R_prior (R''*R): ', func2str(p)], test, test_R(i), 1e-12)
0161    end
0162       
0163    % Call R_priors as RtR_priors
0164    eidors_cache clear
0165    for i = 1:length(R_priors); p = R_priors{i};
0166       inv_mdl= imdl;
0167       inv_mdl.R_prior= p;
0168       RtR= calc_RtR_prior(inv_mdl);
0169       test = std(RtR(:));
0170       unit_test_cmp(['R_prior (R''*R): ', func2str(p)], test, test_R(i), 1e-12)
0171    end
0172 
0173   test_RtR=[ 0.041630544712308
0174              0.003569905276849
0175              0.037777196349047
0176              0.282597508516543
0177              0.006965047087020];
0178       
0179    % Call RtR_priors as RtR_priors
0180    eidors_cache clear
0181    for i = 1:length(R_priors); p = R_priors{i};
0182       inv_mdl= imdl;
0183       inv_mdl.RtR_prior= p;
0184       RtR= calc_RtR_prior(inv_mdl);
0185       if diff(size(RtR))~=0  % non-square
0186          fprintf('RtR_prior: %20s  RtR_condest= NON-SQUARE\n', func2str(p) );
0187       end
0188       test = std(RtR(:));
0189       unit_test_cmp(['RtR_prior (RtR): ', func2str(p)], test, test_RtR(i), 1e-12)
0190    end
0191 
0192    test_R_= [ 0.041630544712308
0193               0.010392960833110
0194               0.040302448507613
0195               0.100634712398168];
0196    % Call RtR_priors as R_priors
0197    eidors_cache clear
0198    for i = 1:length(R_priors); p = R_priors{i};
0199       inv_mdl= imdl;
0200       inv_mdl.RtR_prior= p;
0201       if strcmp(func2str(p), 'prior_TV')
0202          continue; % not fair to ask it to ichol a non-square matrix
0203       end
0204       R= calc_R_prior(inv_mdl);
0205       test = std(R(:));
0206 %     fprintf('RtR_prior: %20s  R_condest= %5.4g\n', func2str(p), condest(R'*R));
0207       if exist('OCTAVE_VERSION') && strcmp(func2str(p),'prior_gaussian_HPF')
0208          eidors_msg('prior_gaussian_HPF test not valid in octave (need ldl)');
0209       else
0210          unit_test_cmp(['RtR_prior (R): ', func2str(p)], test, test_R_(i), 1e-5)
0211       end
0212    end