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
 '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 electrode movement
0019 % 'Choose NF=1.0';  Choose hyperparameter value appropriate for specified noise figure (NF)
0020 
0021 % (C) 2010 Andy Adler. License: GPL version 2 or version 3
0022 % $Id: select_imdl.html 2819 2011-09-07 16:43:11Z aadler $
0023 
0024 if isstr(mdl) && strcmp(mdl,'UNIT_TEST'); do_unit_test; return; end
0025 
0026 if nargin == 1; options = {}; end
0027 
0028 switch mdl.type
0029   case 'inv_model'; inv_mdl = mdl;
0030   case 'fwd_model'; inv_mdl = basic_imdl( mdl );
0031   otherwise;        error('select_imdl: expects inv_model or fwd_model input');
0032 end
0033 
0034 for i=1:length(options);
0035   % split the option on an equals sign
0036   [s,f,tok]= regexp(options{i},'(.[^=]*)=?(.*)');
0037   tok = tok{1};
0038   for t=1:size(tok,1);
0039      opt{t} = options{i}(tok(t,1):tok(t,2));
0040   end
0041 % This is the Mat7 only code:
0042 %   opt= regexp(options{i},'(.[^=]*)=?(.*)','tokens');
0043 %   opt= opt{1};
0044   switch opt{1}
0045     case 'NOSER dif';       inv_mdl = NOSER_dif( inv_mdl );
0046     case 'Basic GN dif';    inv_mdl = Basic_GN_Dif( inv_mdl );
0047     case 'Basic GN abs';    inv_mdl = Basic_GN_Abs( inv_mdl );
0048     case 'Nodal GN dif';    inv_mdl = Nodal_GN_Dif( inv_mdl );
0049     case 'TV solve dif';    inv_mdl = TV_solve_Dif( inv_mdl );
0050     case 'Elec Move GN';    inv_mdl = Elec_Move_GN( inv_mdl );
0051     case 'Choose NF';       inv_mdl = Choose_NF( inv_mdl, str2num(opt{2}) );
0052     
0053     otherwise; error('option {%s} not understood', options{i});
0054   end
0055 end
0056 
0057 function imdl = basic_imdl( fmdl );
0058    imdl.name= 'Basic imdl from select_imdl';
0059    imdl.type= 'inv_model';
0060 
0061    imdl.solve= @aa_inv_solve;
0062    imdl.hyperparameter.value = .01;
0063    imdl.RtR_prior = @laplace_image_prior;
0064    imdl.jacobian_bkgnd.value = 1;
0065    imdl.reconst_type= 'difference';
0066    imdl.fwd_model = fmdl;
0067 
0068 function imdl = NOSER_dif( imdl );
0069    imdl.RtR_prior = @noser_image_prior;
0070    try; imdl = rmfield(imdl,'R_prior'); end
0071    imdl.hyperparameter.value = .03;
0072    imdl.solve= @aa_inv_solve;
0073    imdl.reconst_type= 'difference';
0074 
0075 function imdl = Basic_GN_Dif( imdl );
0076    imdl.RtR_prior = @laplace_image_prior;
0077    try; imdl = rmfield(imdl,'R_prior'); end
0078    imdl.solve= @aa_inv_solve;
0079    imdl.reconst_type= 'difference';
0080 
0081 function imdl = Basic_GN_Abs( imdl );
0082    imdl.RtR_prior = @laplace_image_prior;
0083    try; imdl = rmfield(imdl,'R_prior'); end
0084    imdl.solve= @GN_abs_solve;
0085    imdl.parameters.max_iterations= 10;
0086    imdl.reconst_type= 'absolute';
0087 
0088 function imdl = TV_solve_Dif( imdl );
0089    imdl.R_prior = @ab_calc_tv_prior;
0090    try; imdl = rmfield(imdl,'RtR_prior'); end
0091    imdl.solve= @ab_tv_diff_solve;
0092    imdl.reconst_type= 'difference';
0093    imdl.parameters.max_iterations= 15;
0094    imdl.hyperparameter.value = 1e-1;
0095    imdl.parameters.term_tolerance = 1e-3;
0096 
0097 function imdl = Elec_Move_GN( imdl );
0098    % keep previous model as conductivity jacobian, so it should be ok
0099    imdl.fwd_model.conductivity_jacobian = imdl.fwd_model.jacobian; 
0100    imdl.fwd_model.jacobian = @aa_e_move_jacobian;
0101    imdl.RtR_prior =          @aa_e_move_image_prior;
0102    imdl.solve= @aa_inv_solve;
0103 
0104    MV_prior = 1./mean( std( imdl.fwd_model.nodes ));
0105    imdl.aa_e_move_image_prior.parameters = MV_prior;
0106    imdl.aa_e_move_image_prior.RegC.func = @laplace_image_prior;
0107 %  imdl.aa_e_move_image_prior.RegC.func = @gaussian_HPF_prior; % not OK for 3D
0108    imdl.hyperparameter.value = .03;
0109 
0110    n_elems = size(imdl.fwd_model.elems,1);
0111    imdl.inv_solve.select_parameters = 1:n_elems;
0112 
0113    imdl.prior_use_fwd_not_rec = 1; % for c2f mapping
0114 
0115 function imdl = Nodal_GN_Dif( imdl );
0116    imdl.solve = @nodal_solve;
0117 
0118 
0119 function imdl = Choose_NF( imdl, NF_req );
0120    if ~strcmp(imdl.reconst_type, 'difference');
0121       error('Choose NF only works for difference solvers right now');
0122    end
0123 
0124 % Find 4 elems in mesh ctr to be the NF target elems
0125    xyz_elems = interp_mesh( imdl.fwd_model );
0126    ctr_elems = mean(xyz_elems, 1);
0127    xyz_elems = xyz_elems - ones(size(xyz_elems,1),1)*ctr_elems;
0128    d_elems   = sqrt( sum( xyz_elems.^2, 2 ));
0129    [jnk, e_idx] = sort(d_elems);
0130 
0131    imdl.hyperparameter.tgt_elems = e_idx(1:4);
0132    imdl.hyperparameter.noise_figure = NF_req;
0133 
0134    sv_log = eidors_msg('log_level'); eidors_msg('log_level', 1);
0135    HP = choose_noise_figure( imdl );
0136    eidors_msg('log_level', sv_log);
0137 
0138    imdl.hyperparameter.value = HP;
0139 
0140 
0141 function do_unit_test
0142 % Test difference solvers on lung images
0143    load montreal_data_1995;
0144    imdl = mk_common_model('b2t3',16); 
0145 
0146    for i=1:100; % Break when finished
0147       vh = zc_resp(:,1); vi= zc_resp(:,23);
0148       switch i
0149          case 01;
0150             imdl0 = select_imdl( imdl );
0151          case 02;
0152             imdl0 = select_imdl( imdl.fwd_model );
0153          case 03;
0154             imdl0 = select_imdl( imdl, {'NOSER dif'} );
0155          case 04;
0156             imdl0 = select_imdl( imdl, {'NOSER dif','Choose NF=1.1'});
0157          case 05;
0158             imdl0 = select_imdl( imdl, {'Basic GN dif'} );
0159          case 06;
0160             imdl0 = select_imdl( imdl, {'TV solve dif'} );
0161          case 07;
0162             imdl0 = select_imdl( imdl.fwd_model, {'Basic GN dif','Elec Move GN','Choose NF=0.5'} );
0163          case 08;
0164             imdl0 = select_imdl( imdl, {'Elec Move GN'} );
0165          case 09;
0166             imdl0 = mk_common_model('b2C2',16); 
0167             imdl0 = select_imdl( imdl0, {'Elec Move GN'} );
0168          case 10;
0169             imdl0 = select_imdl( imdl, {'Nodal GN dif'} );
0170          case 11;
0171             imdl0 = select_imdl( imdl, {'Nodal GN dif', 'Choose NF=0.50'} );
0172          case 12;
0173             imdl0 = mk_common_model('b2C2',16); 
0174             imdl0 = select_imdl( imdl0, {'Basic GN dif', 'TV solve dif'} );
0175             imdl0.parameters.max_iterations= 2;
0176             imdl0 = select_imdl( imdl0, {'Choose NF=0.8'} );
0177             [vh,vi] = simulate_movement(mk_image(imdl0), [0;0.5;0.35]);
0178          case 13;
0179             imdl0 = mk_common_model('b2C2',16); 
0180             imdl0 = select_imdl( imdl0, {'Nodal GN dif', 'Choose NF=0.50'} );
0181             [vh,vi] = simulate_movement(mk_image(imdl0), [0;0.5;0.05]);
0182          case 14;
0183             imdl0 = mk_common_model('b2C2',16); 
0184             imdl0 = select_imdl( imdl0, {'Basic GN abs'} );
0185             [vh,vi] = simulate_movement(mk_image(imdl0), [0;0.5;0.05]);
0186          case 15; break
0187       end;
0188 
0189 %     disp([i,imdl0.hyperparameter.value]);
0190       if strcmp( imdl0.reconst_type, 'absolute')
0191          imgr = inv_solve( imdl0, vi);
0192       else
0193          imgr = inv_solve( imdl0, vh, vi);
0194       end
0195       subplot(4,4,i); show_slices( imgr );
0196    end