select_imdl

 SELECT_IMDL: select pre-packaged inverse model features
 inv_mdl= select_imdl( mdl, options )

   mdl - inv_model structure - parameters are replaced as specified
    OR
   mdl - fwd_model - a basic GN difference model is created, and parameters replaced

 OPTIONS => {'opt1','opt2'} options are processed in the order specified

 Available options are:

 'Basic GN dif';   Basic GN one step difference solver with Laplace prior
 'Basic GN abs';   Basic Gauss-Newton absolute solver with Laplace prior
 'NOSER dif';      Basic GN one step difference solver with NOSER prior 
 'Nodal GN dif';   Basic GN solver, solves onto nodes
 'TV solve dif';   Total Variation PDIPM difference solver 
 'Elec Move GN';   One step GN difference solver with compensation for
                   electrode movement. Conductivity prior, hyperparameter
                   and Jacobian will be preserved if specified in the
                   model
 'Choose NF=1.0';  Choose hyperparameter value appropriate for specified noise figure (NF)

0001 function inv_mdl= select_imdl( mdl, options )
0002 % SELECT_IMDL: select pre-packaged inverse model features
0003 % inv_mdl= select_imdl( mdl, options )
0004 %
0005 %   mdl - inv_model structure - parameters are replaced as specified
0006 %    OR
0007 %   mdl - fwd_model - a basic GN difference model is created, and parameters replaced
0008 %
0009 % OPTIONS => {'opt1','opt2'} options are processed in the order specified
0010 %
0011 % Available options are:
0012 %
0013 % 'Basic GN dif';   Basic GN one step difference solver with Laplace prior
0014 % 'Basic GN abs';   Basic Gauss-Newton absolute solver with Laplace prior
0015 % 'NOSER dif';      Basic GN one step difference solver with NOSER prior
0016 % 'Nodal GN dif';   Basic GN solver, solves onto nodes
0017 % 'TV solve dif';   Total Variation PDIPM difference solver
0018 % 'Elec Move GN';   One step GN difference solver with compensation for
0019 %                   electrode movement. Conductivity prior, hyperparameter
0020 %                   and Jacobian will be preserved if specified in the
0021 %                   model
0022 % 'Choose NF=1.0';  Choose hyperparameter value appropriate for specified noise figure (NF)
0023 
0024 % (C) 2010 Andy Adler. License: GPL version 2 or version 3
0025 % $Id: select_imdl.m 5554 2017-06-18 16:47:36Z aadler $
0026 
0027 if ischar(mdl) && strcmp(mdl,'UNIT_TEST'); do_unit_test; return; end
0028 
0029 if nargin == 1; options = {}; end
0030 
0031 switch mdl.type
0032   case 'inv_model'; inv_mdl = mdl;
0033   case 'fwd_model'; inv_mdl = basic_imdl( mdl );
0034   otherwise;        error('select_imdl: expects inv_model or fwd_model input');
0035 end
0036 
0037 for i=1:length(options);
0038   % split the option on an equals sign
0039   [s,f,tok]= regexp(options{i},'(.[^=]*)=?(.*)');
0040   tok = tok{1};
0041   for t=1:size(tok,1);
0042      opt{t} = options{i}(tok(t,1):tok(t,2));
0043   end
0044 % This is the Mat7 only code:
0045 %   opt= regexp(options{i},'(.[^=]*)=?(.*)','tokens');
0046 %   opt= opt{1};
0047   switch opt{1}
0048     case 'NOSER dif';       inv_mdl = NOSER_dif( inv_mdl );
0049     case 'Basic GN dif';    inv_mdl = Basic_GN_Dif( inv_mdl );
0050     case 'Basic GN abs';    inv_mdl = Basic_GN_Abs( inv_mdl );
0051     case 'Nodal GN dif';    inv_mdl = Nodal_GN_Dif( inv_mdl );
0052     case 'TV solve dif';    inv_mdl = TV_solve_Dif( inv_mdl );
0053     case 'Elec Move GN';    inv_mdl = Elec_Move_GN( inv_mdl );
0054     case 'Choose NF';       inv_mdl = Choose_NF( inv_mdl, str2num(opt{2}) );
0055     
0056     otherwise; error('option {%s} not understood', options{i});
0057   end
0058 end
0059 
0060 function imdl = basic_imdl( fmdl );
0061    imdl.name= 'Basic imdl from select_imdl';
0062    imdl.type= 'inv_model';
0063 
0064    imdl.solve= 'eidors_default';
0065    imdl.hyperparameter.value = .01;
0066    imdl.RtR_prior = 'eidors_default';
0067    imdl.jacobian_bkgnd.value = 1;
0068    imdl.reconst_type= 'difference';
0069    imdl.fwd_model = fmdl;
0070 
0071 function imdl = NOSER_dif( imdl );
0072    imdl.RtR_prior = @prior_noser;
0073    try; imdl = rmfield(imdl,'R_prior'); end
0074    imdl.hyperparameter.value = .03;
0075    imdl.solve= @inv_solve_diff_GN_one_step;
0076    imdl.reconst_type= 'difference';
0077 
0078 function imdl = Basic_GN_Dif( imdl );
0079    imdl.RtR_prior = @prior_laplace;
0080    try; imdl = rmfield(imdl,'R_prior'); end
0081    imdl.solve= @inv_solve_diff_GN_one_step;
0082    imdl.reconst_type= 'difference';
0083 
0084 function imdl = Basic_GN_Abs( imdl );
0085    imdl.RtR_prior = @prior_laplace;
0086    try; imdl = rmfield(imdl,'R_prior'); end
0087    imdl.solve= @inv_solve_gn;
0088    imdl.inv_solve_gn.max_iterations= 10;
0089    imdl.reconst_type= 'absolute';
0090 
0091 function imdl = TV_solve_Dif( imdl );
0092    imdl.R_prior = @prior_TV;
0093    try; imdl = rmfield(imdl,'RtR_prior'); end
0094    imdl.solve= @inv_solve_TV_pdipm;
0095    imdl.reconst_type= 'difference';
0096    imdl.inv_solve.max_iterations= 15;
0097    imdl.hyperparameter.value = 1e-1;
0098    imdl.inv_solve.term_tolerance = 1e-3;
0099 
0100 function imdl = Elec_Move_GN( imdl );
0101    % keep previous model as conductivity jacobian, so it should be ok
0102    imdl.fwd_model.conductivity_jacobian = imdl.fwd_model.jacobian; 
0103    imdl.fwd_model.jacobian = @jacobian_movement;
0104    
0105    % keep previous prior
0106    try
0107        imdl.prior_movement.RegC.func = imdl.RtR_prior;
0108    catch
0109        imdl.prior_movement.RegC.func = @prior_laplace;
0110        %  imdl.prior_movement.RegC.func = @prior_gaussian_HPF; % not OK for 3D
0111    end
0112    
0113    imdl.RtR_prior =          @prior_movement;
0114    MV_prior = 1./mean( std( imdl.fwd_model.nodes ));
0115    imdl.prior_movement.parameters = MV_prior;
0116    
0117    imdl.solve= @inv_solve_diff_GN_one_step;
0118    
0119    % keep hyperparameter
0120    try 
0121        hp = imdl.hyperparameter.value;
0122    catch
0123        hp = 0.03;
0124    end
0125    imdl.hyperparameter.value = hp;
0126    
0127    try
0128       n_elems = size(imdl.rec_model.coarse2fine,2);
0129    catch
0130       n_elems = size(imdl.fwd_model.elems,1);
0131    end
0132    imdl.inv_solve.select_parameters = 1:n_elems;
0133 
0134    imdl.prior_use_fwd_not_rec = 1; % for c2f mapping
0135 
0136 function imdl = Nodal_GN_Dif( imdl );
0137    imdl.solve = @nodal_solve;
0138 
0139 
0140 function imdl = Choose_NF( imdl, NF_req );
0141    if ~strcmp(imdl.reconst_type, 'difference');
0142       error('Choose NF only works for difference solvers right now');
0143    end
0144 
0145 % Find 4 elems in mesh ctr to be the NF target elems
0146    xyz_elems = interp_mesh( imdl.fwd_model );
0147    ctr_elems = mean(xyz_elems, 1);
0148    xyz_elems = xyz_elems - ones(size(xyz_elems,1),1)*ctr_elems;
0149    d_elems   = sqrt( sum( xyz_elems.^2, 2 ));
0150    [jnk, e_idx] = sort(d_elems);
0151 
0152    imdl.hyperparameter.tgt_elems = e_idx(1:4);
0153    imdl.hyperparameter.noise_figure = NF_req;
0154 
0155    sv_log = eidors_msg('log_level'); eidors_msg('log_level', 1);
0156    HP = choose_noise_figure( imdl );
0157    eidors_msg('log_level', sv_log);
0158 
0159    imdl.hyperparameter.value = HP;
0160 
0161 
0162 function do_unit_test
0163 % Test difference solvers on lung images
0164    load montreal_data_1995;
0165    imdl = mk_common_model('b2t3',16); 
0166 
0167    for i=1:100; % Break when finished
0168       vh = zc_resp(:,1); vi= zc_resp(:,23);
0169       fprintf('solver #%d\n',i);
0170       switch i
0171          case 01;
0172             imdl0 = select_imdl( imdl );
0173          case 02;
0174             imdl0 = select_imdl( imdl.fwd_model );
0175          case 03;
0176             imdl0 = select_imdl( imdl, {'NOSER dif'} );
0177          case 04;
0178             imdl0 = select_imdl( imdl, {'NOSER dif','Choose NF=1.1'});
0179          case 05;
0180             imdl0 = select_imdl( imdl, {'Basic GN dif'} );
0181          case 06;
0182             imdl0 = select_imdl( imdl, {'TV solve dif'} );
0183          case 07;
0184             imdl0 = select_imdl( imdl.fwd_model, {'Basic GN dif','Elec Move GN'} );
0185          case 08;
0186             imdl0 = select_imdl( imdl.fwd_model, {'Basic GN dif','Elec Move GN','Choose NF=0.5'} );
0187          case 09;
0188             imdl0 = select_imdl( imdl, {'Elec Move GN'} );
0189          case 10;
0190             imdl0 = mk_common_model('b2C2',16); 
0191             imdl0 = select_imdl( imdl0, {'Elec Move GN'} );
0192          case 11;
0193             imdl0 = select_imdl( imdl, {'Nodal GN dif'} );
0194          case 12;
0195             imdl0 = select_imdl( imdl, {'Nodal GN dif', 'Choose NF=0.50'} );
0196          case 13;
0197             imdl0 = mk_common_model('b2C2',16); 
0198             imdl0 = select_imdl( imdl0, {'Basic GN dif', 'TV solve dif'} );
0199             imdl0.inv_solve.max_iterations= 2;
0200             imdl0 = select_imdl( imdl0, {'Choose NF=0.8'} );
0201             [vh,vi] = simulate_movement(mk_image(imdl0), [0;0.5;0.35]);
0202          case 14;
0203             imdl0 = mk_common_model('b2C2',16); 
0204             imdl0 = select_imdl( imdl0, {'Nodal GN dif', 'Choose NF=0.50'} );
0205             [vh,vi] = simulate_movement(mk_image(imdl0), [0;0.5;0.05]);
0206          case 15;
0207             imdl0 = mk_common_model('b2C2',16); 
0208             imdl0 = select_imdl( imdl0, {'Basic GN abs'} );
0209             [vh,vi] = simulate_movement(mk_image(imdl0), [0;0.5;0.05]);
0210          case 16; break
0211       end;
0212 
0213 %     disp([i,imdl0.hyperparameter.value]);
0214       if strcmp( imdl0.reconst_type, 'absolute')
0215          imgr = inv_solve( imdl0, vi);
0216       else
0217          imgr = inv_solve( imdl0, vh, vi);
0218       end
0219       subplot(4,4,i); show_slices( imgr );
0220    end