calc_R_prior

 R = calc_R_prior( inv_model, varargin )
 CALC_R_PRIOR: calculate regularization matrix R
   The image prior is matrix n_elem x ??? 
 
 NOTE: This function is mostly just a hack. It's not obvious
    how to to this correctly.

 calc_R_prior can be called as
    R_prior= calc_R_prior( inv_model, ... )

 and will call the function inv_model.R_prior
 parameters to R_prior should be passed in the field
 inv_model.R_prior_function_name.parameters
 
 If inv_model.R_prior is a matrix, calc_R_prior will return that matrix,
 possibly correcting for coarse2fine

 R_prior      calculated regularization prior R
 inv_model    is an inv_model structure

0001 function R_prior = calc_R_prior( inv_model, varargin )
0002 % R = calc_R_prior( inv_model, varargin )
0003 % CALC_R_PRIOR: calculate regularization matrix R
0004 %   The image prior is matrix n_elem x ???
0005 %
0006 % NOTE: This function is mostly just a hack. It's not obvious
0007 %    how to to this correctly.
0008 %
0009 % calc_R_prior can be called as
0010 %    R_prior= calc_R_prior( inv_model, ... )
0011 %
0012 % and will call the function inv_model.R_prior
0013 % parameters to R_prior should be passed in the field
0014 % inv_model.R_prior_function_name.parameters
0015 %
0016 % If inv_model.R_prior is a matrix, calc_R_prior will return that matrix,
0017 % possibly correcting for coarse2fine
0018 %
0019 % R_prior      calculated regularization prior R
0020 % inv_model    is an inv_model structure
0021 
0022 % (C) 2005-2008 Andy Adler. License: GPL version 2 or version 3
0023 % $Id: calc_R_prior.m 6035 2019-12-30 19:17:12Z aadler $
0024 
0025 if ischar(inv_model) && strcmp(inv_model,'UNIT_TEST'); do_unit_test; return; end
0026 
0027 inv_model = rec_or_fwd_model( inv_model);
0028 
0029 
0030 if isfield(inv_model,'R_prior')
0031    if isnumeric(inv_model.R_prior)
0032       R_prior = inv_model.R_prior;
0033    else
0034       R_prior= eidors_cache( inv_model.R_prior, inv_model );
0035    end
0036 elseif isfield(inv_model,'RtR_prior')
0037    R_prior = eidors_cache(@calc_from_RtR_prior, inv_model,'calc_R_prior');
0038 else
0039    error('calc_R_prior: neither R_prior or RtR_prior func provided');
0040 end
0041 
0042 if isfield(inv_model.fwd_model,'coarse2fine')
0043    c2f= inv_model.fwd_model.coarse2fine;
0044    if size(R_prior,1)==size(c2f,1)
0045 %     we need to take into account coarse2fine - using a reasonable tol
0046       f2c = c2f';         
0047       R_prior = f2c*R_prior*c2f;
0048    end
0049 end
0050 
0051 function R_prior = calc_from_RtR_prior(inv_model)
0052    % The user has provided an RtR prior. We can use this to
0053    % get R =RtR^(1/2). Not that this is non unique
0054    if isnumeric(inv_model.RtR_prior)
0055       RtR_prior = inv_model.RtR_prior;
0056    else
0057       RtR_prior= eidors_cache( inv_model.RtR_prior, inv_model );
0058    end
0059    
0060    if exist('ldl') % not available for octave
0061       [L,D,P] = ldl(RtR_prior);
0062       R_prior = sqrt(D)*L'*P';
0063    else
0064       R_prior = sparse(size(RtR_prior,1),size(RtR_prior,2));
0065       [R,f,P] = chol(RtR_prior); % will often be sigular
0066       if f==0; f=size(RtR_prior,1); end % successful
0067       if f==1; f=size(R,1); end
0068       R_prior(1:f,:) = R*P';
0069    end
0070 
0071 function inv_model = rec_or_fwd_model( inv_model);
0072 
0073    if isfield(inv_model,'rec_model');
0074       use_rec_model = 1;
0075       try if inv_model.prior_use_fwd_not_rec== 1;
0076          use_rec_model = 0;
0077       end; end
0078 
0079       if use_rec_model
0080           % copy the normalize flag from the fwd_model to prevent warnings
0081          inv_model.rec_model = mdl_normalize(inv_model.rec_model, ...
0082                                        mdl_normalize(inv_model.fwd_model));
0083          inv_model.fwd_model= inv_model.rec_model;
0084          inv_model= rmfield(inv_model,'rec_model');
0085       end
0086    end
0087 
0088 function do_unit_test
0089 
0090    priors={'prior_tikhonov',
0091            'prior_laplace',
0092            'prior_gaussian_HPF',
0093            'prior_movement'};
0094    models= {'a2c0','f2c2','n3r2'};
0095 
0096    for j= 1:length(models);
0097       imdl = mk_common_model(models{j},16);
0098       for i = 1:length(priors);
0099          imdl.RtR_prior = feval(priors{i},imdl);
0100 
0101          R   = calc_R_prior(imdl);
0102          RtR = calc_RtR_prior(imdl);  
0103          unit_test_cmp(['P:',priors{i}], R'*R, RtR, 1e-10);
0104       end
0105    end
0106 
0107