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.html 2819 2011-09-07 16:43:11Z aadler $
0015 
0016 if isstr(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 
0035 J= calc_jacobian(fmdl,bkgnd_img);
0036 imdl.solve = @solve_use_matrix;
0037 S = eidors_obj('get-cache',mdl,'TSVD_S');
0038 if isempty(S)
0039     [U,S,V] = svd(J,'econ');
0040     eidors_obj('set-cache',mdl,'TSVD_S',S);
0041     eidors_obj('set-cache',mdl,'TSVD_V',V);
0042     eidors_obj('set-cache',mdl,'TSVD_U',U);
0043 else
0044     U = eidors_obj('get-cache',mdl,'TSVD_U');
0045     V = eidors_obj('get-cache',mdl,'TSVD_V');
0046 end
0047 s = diag(S);
0048 s = s./s(1);
0049 N = find(s >= hp/100,1,'last'); %disp(N);
0050 t= 1:N; % tsvd
0051 RM = eidors_obj('get-cache',mdl,'TSVD_RM',N);
0052 if isempty(RM)
0053     RM  = (V(:,t)/S(t,t))*U(:,t)';
0054     eidors_obj('set-cache',mdl,'TSVD_RM',RM,N);
0055 end
0056 
0057 
0058 
0059 function do_unit_test
0060 mdl = mk_common_model('b3t2r',[16,1]);
0061 RM = calc_TSVD_RM(mdl);