calc_TSVD_RM

 CALC_TSVD_RM: Calculated truncated Jacobian SVD reconstruction matrix
   RM= calc_TSVD_RM(mdl, hp)
 
   RM    => reconstruction matrix
   mdl   => an inverse or forward model structure 
   hp    => hyperparameter. Ratio between the last kept and the first SV
            in percent (default: .1)
 
 SEE ALSO: inv_solve_TSVD, solve_use_matrix

0001 function RM= calc_TSVD_RM(mdl,hp)
0002 % CALC_TSVD_RM: Calculated truncated Jacobian SVD reconstruction matrix
0003 %   RM= calc_TSVD_RM(mdl, hp)
0004 %
0005 %   RM    => reconstruction matrix
0006 %   mdl   => an inverse or forward model structure
0007 %   hp    => hyperparameter. Ratio between the last kept and the first SV
0008 %            in percent (default: .1)
0009 %
0010 % SEE ALSO: inv_solve_TSVD, solve_use_matrix
0011 
0012 
0013 % (C) 2011 Bartlomiej Grychtol. Licenced under GPL v2 or v3
0014 % $Id: calc_TSVD_RM.m 5431 2017-04-27 15:29:35Z aadler $
0015 
0016 if ischar(mdl) && strcmp(mdl, 'UNIT_TEST'), do_unit_test, return,end
0017 
0018 
0019 
0020 switch (mdl.type)
0021     case 'inv_model'
0022         fmdl = mdl.fwd_model;
0023         bkgnd_img = calc_jacobian_bkgnd(mdl);
0024     case 'fwd_model'
0025         fmdl = mdl;
0026         bkgnd_img = mk_image( fmdl,1) ;
0027     otherwise
0028         eidors_msg('calc_TSVD_RM: Object type not supported');
0029 end
0030 
0031 if nargin < 2
0032     hp = .1;
0033 end
0034 J= calc_jacobian(bkgnd_img);
0035 copt.fstr = 'calc_TSVD_RM';
0036 RM = eidors_cache(@calc_RM,{J, hp}, copt);
0037 
0038 function RM = calc_RM(J, hp)
0039 
0040 imdl.solve = @solve_use_matrix;
0041 copt.cache_obj = J;
0042 copt.fstr = 'svd';
0043 
0044 [U,S,V] = eidors_cache(@svd,{J, 'econ'},copt);
0045 
0046 s = diag(S);
0047 s = s./s(1);
0048 N = find(s >= hp/100,1,'last'); %disp(N);
0049 t= 1:N; % tsvd
0050 RM  = (V(:,t)/S(t,t))*U(:,t)';
0051 
0052 
0053 
0054 function do_unit_test
0055    mdl = mk_common_model('b3t2r',[16,1]);
0056    RM = calc_TSVD_RM(mdl);
0057    J  = calc_jacobian(mk_image(mdl));
0058    unit_test_cmp('RM size', size(RM), size(J'));
0059 
0060 % TODO (AA), May2015): Add better tests