0001 function R_prior = calc_R_prior( inv_model, varargin )
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022 if ischar(inv_model) && strcmp(inv_model,'UNIT_TEST'); do_unit_test; return; end
0023
0024 inv_model = rec_or_fwd_model( inv_model);
0025
0026
0027 if isfield(inv_model,'R_prior')
0028 if isnumeric(inv_model.R_prior)
0029 R_prior = inv_model.R_prior;
0030 else
0031 R_prior= eidors_cache( inv_model.R_prior, inv_model );
0032 end
0033 elseif isfield(inv_model,'RtR_prior')
0034 R_prior = eidors_cache(@calc_from_RtR_prior, inv_model,'calc_R_prior');
0035 else
0036 error('calc_R_prior: neither R_prior or RtR_prior func provided');
0037 end
0038
0039 if isfield(inv_model.fwd_model,'coarse2fine')
0040 c2f= inv_model.fwd_model.coarse2fine;
0041 if size(R_prior,1)==size(c2f,1)
0042
0043 f2c = c2f';
0044 R_prior = f2c*R_prior*c2f;
0045 end
0046 end
0047
0048 function R_prior = calc_from_RtR_prior(inv_model)
0049
0050
0051 if isnumeric(inv_model.RtR_prior)
0052 RtR_prior = inv_model.RtR_prior;
0053 else
0054 RtR_prior= eidors_cache( inv_model.RtR_prior, inv_model );
0055 end
0056
0057 [L,D,P] = ldl(RtR_prior);
0058 R_prior = sqrt(D)*L'*P';
0059 [L,D,P] = ldl(RtR_prior);
0060
0061 function inv_model = rec_or_fwd_model( inv_model);
0062
0063 if isfield(inv_model,'rec_model');
0064 use_rec_model = 1;
0065 try if inv_model.prior_use_fwd_not_rec== 1;
0066 use_rec_model = 0;
0067 end; end
0068
0069 if use_rec_model
0070
0071 inv_model.rec_model = mdl_normalize(inv_model.rec_model, ...
0072 mdl_normalize(inv_model.fwd_model));
0073 inv_model.fwd_model= inv_model.rec_model;
0074 inv_model= rmfield(inv_model,'rec_model');
0075 end
0076 end
0077
0078 function do_unit_test
0079
0080 priors={'prior_tikhonov',
0081 'prior_laplace',
0082 'prior_gaussian_HPF',
0083 'prior_movement'};
0084 models= {'a2c0','f2c2','n3r2'};
0085
0086 for j= 1:length(models);
0087 imdl = mk_common_model(models{j},16);
0088 for i = 1:length(priors);
0089 imdl.RtR_prior = feval(priors{i},imdl);
0090
0091 R = calc_R_prior(imdl);
0092 RtR = calc_RtR_prior(imdl);
0093 unit_test_cmp(['P:',priors{i}], R'*R, RtR, 1e-10);
0094 end
0095 end
0096
0097