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 OR image

 OPTIONS => {'opt1','opt2'} options are processed in the order specified
 if OPTIONS is char, treated as first option

 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 
 'GREIT':          Algorithm of mk_GREIT_model
   parameters provided as select_imdl(fmdl,{'GREIT:NF=0.3 64x64 rad=0.25'});
   defaults NF=0.5 32x32 rad=0.2
 'GREIT':          3D GREIT if parameters:
     {'GREIT:NF=0.3 32x32x6'});   % default vrange=.2
     {'GREIT:NF=0.3 vrange=0.3 32x32x6'});
   Here 6 slices cover +/-vrange*max(diam) from verticle centre

 '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,opt_out]= 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 OR image
0008 %
0009 % OPTIONS => {'opt1','opt2'} options are processed in the order specified
0010 % if OPTIONS is char, treated as first option
0011 %
0012 % Available options are:
0013 %
0014 % 'Basic GN dif';   Basic GN one step difference solver with Laplace prior
0015 % 'Basic GN abs';   Basic Gauss-Newton absolute solver with Laplace prior
0016 % 'NOSER dif';      Basic GN one step difference solver with NOSER prior
0017 % 'Nodal GN dif';   Basic GN solver, solves onto nodes
0018 % 'TV solve dif':   Total Variation PDIPM difference solver
0019 % 'GREIT':          Algorithm of mk_GREIT_model
0020 %   parameters provided as select_imdl(fmdl,{'GREIT:NF=0.3 64x64 rad=0.25'});
0021 %   defaults NF=0.5 32x32 rad=0.2
0022 % 'GREIT':          3D GREIT if parameters:
0023 %     {'GREIT:NF=0.3 32x32x6'});   % default vrange=.2
0024 %     {'GREIT:NF=0.3 vrange=0.3 32x32x6'});
0025 %   Here 6 slices cover +/-vrange*max(diam) from verticle centre
0026 %
0027 % 'Elec Move GN':   One step GN difference solver with compensation for
0028 %                   electrode movement. Conductivity prior, hyperparameter
0029 %                   and Jacobian will be preserved if specified in the
0030 %                   model
0031 % 'Choose NF=1.0';  Choose hyperparameter value appropriate for specified noise figure (NF)
0032 
0033 % (C) 2010-2024 Andy Adler. License: GPL version 2 or version 3
0034 % $Id: select_imdl.m 6972 2024-10-03 11:40:00Z aadler $
0035 
0036 if ischar(mdl) && strcmp(mdl,'UNIT_TEST'); do_unit_test; return; end
0037 
0038 if nargin == 1; options = {}; end
0039 
0040 switch mdl.type
0041   case 'inv_model'; inv_mdl = mdl;
0042   case 'fwd_model'; inv_mdl = basic_imdl( mdl );
0043   case 'image';     inv_mdl = basic_imdl( mdl.fwd_model, mdl.elem_data );
0044   otherwise;        error('select_imdl: expects inv_model or fwd_model input');
0045 end
0046 
0047 opt_out = struct(); % default setting
0048 if ischar(options); options = {options}; end
0049 for i=1:length(options);
0050 % Pre matlab v7 code
0051 % [s,f,tok]= regexp(options{i},'(.[^=:]*)[=:]?(.*)');
0052 % tok = tok{1};
0053 % for t=1:size(tok,1);
0054 %    opt{t} = options{i}(tok(t,1):tok(t,2));
0055 % end
0056   % split the option on an equals sign
0057   opt= regexp(options{i},'(.[^=:]*)[=:]?(.*)','tokens'); opt= opt{1};
0058   switch opt{1}
0059     case 'NOSER dif';       inv_mdl = NOSER_dif( inv_mdl );
0060     case 'Basic GN dif';    inv_mdl = Basic_GN_Dif( inv_mdl );
0061     case 'Basic GN abs';    inv_mdl = Basic_GN_Abs( inv_mdl );
0062     case 'Nodal GN dif';    inv_mdl = Nodal_GN_Dif( inv_mdl );
0063     case 'TV solve dif';    inv_mdl = TV_solve_Dif( inv_mdl );
0064     case 'Elec Move GN';    inv_mdl = Elec_Move_GN( inv_mdl );
0065     case 'Choose NF';       inv_mdl = Choose_NF( inv_mdl, str2num(opt{2}) );
0066     case 'GREIT';           [inv_mdl,opt_out] = GREIT(inv_mdl, opt{2});
0067     
0068     otherwise; error('option {%s} not understood', options{i});
0069   end
0070 end
0071 
0072 % standard field order
0073 inv_mdl = eidors_obj('inv_model',inv_mdl);
0074 
0075 % check we are giving back a good specimen
0076 valid_inv_model(inv_mdl);
0077 
0078 
0079 function imdl = basic_imdl( fmdl, elem_data );
0080    if nargin<=1; elem_data = 1; end
0081    imdl.name= 'Basic imdl from select_imdl';
0082    imdl.type= 'inv_model';
0083 
0084    imdl.solve= 'eidors_default';
0085    imdl.hyperparameter.value = .01;
0086    imdl.RtR_prior = 'eidors_default';
0087    imdl.jacobian_bkgnd.value = elem_data;
0088    imdl.reconst_type= 'difference';
0089    imdl.fwd_model = fmdl;
0090 
0091 function imdl = NOSER_dif( imdl );
0092    imdl.RtR_prior = @prior_noser;
0093    try; imdl = rmfield(imdl,'R_prior'); end
0094    imdl.hyperparameter.value = .03;
0095    imdl.solve= @inv_solve_diff_GN_one_step;
0096    imdl.reconst_type= 'difference';
0097 
0098 function imdl = Basic_GN_Dif( imdl );
0099    imdl.RtR_prior = @prior_laplace;
0100    try; imdl = rmfield(imdl,'R_prior'); end
0101    imdl.solve= @inv_solve_diff_GN_one_step;
0102    imdl.reconst_type= 'difference';
0103 
0104 function imdl = Basic_GN_Abs( imdl );
0105    imdl.RtR_prior = @prior_laplace;
0106    try; imdl = rmfield(imdl,'R_prior'); end
0107    imdl.solve= @inv_solve_gn;
0108    imdl.inv_solve_gn.max_iterations= 10;
0109    imdl.reconst_type= 'absolute';
0110 
0111 function imdl = TV_solve_Dif( imdl );
0112    imdl.R_prior = @prior_TV;
0113    try; imdl = rmfield(imdl,'RtR_prior'); end
0114    imdl.solve= @inv_solve_TV_pdipm;
0115    imdl.reconst_type= 'difference';
0116    imdl.inv_solve.max_iterations= 15;
0117    imdl.hyperparameter.value = 1e-1;
0118    imdl.inv_solve.term_tolerance = 1e-3;
0119 
0120 function imdl = Elec_Move_GN( imdl );
0121    % keep previous model as conductivity jacobian, so it should be ok
0122    imdl.fwd_model.conductivity_jacobian = imdl.fwd_model.jacobian; 
0123    imdl.fwd_model.jacobian = @jacobian_movement;
0124    
0125    % keep previous prior
0126    try
0127        imdl.prior_movement.RegC.func = imdl.RtR_prior;
0128    catch
0129        imdl.prior_movement.RegC.func = @prior_laplace;
0130        %  imdl.prior_movement.RegC.func = @prior_gaussian_HPF; % not OK for 3D
0131    end
0132    
0133    imdl.RtR_prior =          @prior_movement;
0134    MV_prior = 1./mean( std( imdl.fwd_model.nodes ));
0135    imdl.prior_movement.parameters = MV_prior;
0136    
0137    imdl.solve= @inv_solve_diff_GN_one_step;
0138    
0139    % keep hyperparameter
0140    try 
0141        hp = imdl.hyperparameter.value;
0142    catch
0143        hp = 0.03;
0144    end
0145    imdl.hyperparameter.value = hp;
0146    
0147    try
0148       n_elems = size(imdl.rec_model.coarse2fine,2);
0149    catch
0150       n_elems = size(imdl.fwd_model.elems,1);
0151    end
0152    imdl.inv_solve.select_parameters = 1:n_elems;
0153 
0154    imdl.prior_use_fwd_not_rec = 1; % for c2f mapping
0155 
0156 function imdl = Nodal_GN_Dif( imdl );
0157    imdl.solve = @nodal_solve;
0158 
0159 function [imdl,opt] = GREIT(inv_mdl, opts)
0160    opt.square_pixels = true;
0161    opt.noise_figure = 0.5;
0162    opt.imgsz = [32 32];
0163    % Allow error compensation afterward
0164    opt.keep_model_components = true;
0165    radius = 0.2;
0166 
0167    vrange = 0.2; % TODO: add adjust
0168    matches = regexp(opts,'\S+','match');
0169    for i=1:length(matches); mi = matches{i};
0170       tok= regexp(mi,'NF=([\d\.]+)','tokens');
0171       if ~isempty(tok);
0172          opt.noise_figure = str2num(tok{1}{1});
0173       end
0174 
0175       tok= regexp(mi,'rad=([\d\.]+)','tokens');
0176       if ~isempty(tok);
0177          radius = str2num(tok{1}{1});
0178       end
0179 
0180       tok= regexp(mi,'vrange=([\d\.]+)','tokens');
0181       if ~isempty(tok);
0182          vrange = str2num(tok{1}{1});
0183       end
0184 
0185       tok= regexp(mi,'(\d+)x(\d+)','tokens');
0186       if ~isempty(tok);
0187          opt.imgsz = [str2num(tok{1}{1}), str2num(tok{1}{2})];
0188       end
0189 
0190       % 3D GREIT
0191       tok= regexp(mi,'(\d+)x(\d+)x(\d+)','tokens');
0192       if ~isempty(tok);
0193          tokA=str2num(tok{1}{1});
0194          tokB=str2num(tok{1}{2});
0195          tokC=str2num(tok{1}{3});
0196          [minf,maxf] = bounds(inv_mdl.fwd_model.nodes);
0197          vctr = mean([minf(3),maxf(3)]);
0198          mrad = max( maxf(1:2) - minf(1:2) );
0199 
0200          vopt.imgsz = [tokA, tokB];
0201          vopt.square_pixels = true;
0202          vopt.zvec = linspace(-1,1,tokC+1)*mrad*vrange + vctr;
0203          vopt.save_memory = 2;
0204 
0205 % GREIT 3D with a 1x32 electrode layout
0206          [inv_mdl,opt.distr] = GREIT3D_distribution(inv_mdl, vopt);
0207       end
0208    end
0209 
0210    imdl= mk_GREIT_model(inv_mdl, radius, [], opt);
0211    
0212 
0213 function imdl = Choose_NF( imdl, NF_req );
0214    if ~strcmp(imdl.reconst_type, 'difference');
0215       error('Choose NF only works for difference solvers right now');
0216    end
0217 
0218 % Find 4 elems in mesh ctr to be the NF target elems
0219    xyz_elems = interp_mesh( imdl.fwd_model );
0220    ctr_elems = mean(xyz_elems, 1);
0221    xyz_elems = xyz_elems - ones(size(xyz_elems,1),1)*ctr_elems;
0222    d_elems   = sqrt( sum( xyz_elems.^2, 2 ));
0223    [jnk, e_idx] = sort(d_elems);
0224 
0225    imdl.hyperparameter.tgt_elems = e_idx(1:4);
0226    imdl.hyperparameter.noise_figure = NF_req;
0227 
0228    sv_log = eidors_msg('log_level'); eidors_msg('log_level', 1);
0229    HP = choose_noise_figure( imdl );
0230    eidors_msg('log_level', sv_log);
0231 
0232    imdl.hyperparameter.value = HP;
0233 
0234 
0235 function do_unit_test
0236 % Test difference solvers on lung images
0237    load montreal_data_1995;
0238    imdl = mk_common_model('b2t3',16); 
0239 
0240    i=0;while true; i=i+1; % Break when finished
0241       vh = zc_resp(:,1); vi= zc_resp(:,23);
0242       fprintf('solver #%d\n',i);
0243       switch i
0244          case 01;
0245             imdl0 = select_imdl( imdl );
0246          case 02;
0247             imdl0 = select_imdl( imdl.fwd_model );
0248          case 03;
0249             imdl0 = select_imdl( imdl, {'NOSER dif'} );
0250          case 04;
0251             imdl0 = select_imdl( imdl, {'NOSER dif','Choose NF=1.1'});
0252          case 05;
0253             imdl0 = select_imdl( imdl, {'Basic GN dif'} );
0254          case 06;
0255             imdl0 = select_imdl( imdl, {'TV solve dif'} );
0256          case 07;
0257             imdl0 = select_imdl( imdl.fwd_model, {'Basic GN dif','Elec Move GN'} );
0258          case 08;
0259             imdl0 = select_imdl( imdl.fwd_model, {'Basic GN dif','Elec Move GN','Choose NF=0.5'} );
0260          case 09;
0261             imdl0 = select_imdl( imdl, {'Elec Move GN'} );
0262          case 10;
0263             imdl0 = mk_common_model('b2C2',16); 
0264             imdl0 = select_imdl( imdl0, {'Elec Move GN'} );
0265          case 11;
0266             imdl0 = select_imdl( imdl, {'Nodal GN dif'} );
0267          case 12;
0268             imdl0 = select_imdl( imdl, {'Nodal GN dif', 'Choose NF=0.50'} );
0269          case 13;
0270             fmdl = mk_library_model('adult_male_16el');
0271             [stim,msel] = mk_stim_patterns(16,1,[0,1],[0,1],{},1);
0272             fmdl.stimulation = stim; 
0273             fmdl.meas_select = msel; 
0274             imdl0 = select_imdl( fmdl, 'GREIT');
0275 
0276          case 14;
0277             imdl0 = mk_common_model('b2C2',16); 
0278             imdl0 = select_imdl( imdl0, {'Basic GN dif', 'TV solve dif'} );
0279             imdl0.inv_solve.max_iterations= 2;
0280             imdl0 = select_imdl( imdl0, {'Choose NF=0.8'} );
0281             [vh,vi] = simulate_movement(mk_image(imdl0), [0;0.5;0.35]);
0282          case 15;
0283             imdl0 = mk_common_model('b2C2',16); 
0284             imdl0 = select_imdl( imdl0, {'Nodal GN dif', 'Choose NF=0.50'} );
0285             [vh,vi] = simulate_movement(mk_image(imdl0), [0;0.5;0.05]);
0286          case 16;
0287             imdl0 = mk_common_model('b2C2',16); 
0288             imdl0 = select_imdl( imdl0, {'Basic GN abs'} );
0289             [vh,vi] = simulate_movement(mk_image(imdl0), [0;0.5;0.05]);
0290          case 17;
0291             imdl = mk_common_model('n3r2',[16,2]); fmdl = imdl.fwd_model;
0292             fmdl.nodes(:,3) = fmdl.nodes(:,3) - 1.5;
0293             imdl0= select_imdl(fmdl,{'GREIT:NF=0.3 64x64 rad=0.25'});
0294             img = mk_image(fmdl,1); vh = fwd_solve(img);
0295             img.elem_data(fmdl.demo_targ_elems.A) = 1.2;
0296             img.elem_data(fmdl.demo_targ_elems.B) = 0.8;
0297                                     vi = fwd_solve(img);
0298            
0299          case 18;
0300             imdl = mk_common_model('n3r2',[16,2]); fmdl = imdl.fwd_model;
0301             fmdl.nodes(:,3) = fmdl.nodes(:,3) - 1.5;
0302             imdl0= select_imdl(fmdl,{'GREIT:NF=0.3 32x32 rad=0.25'});
0303             img = mk_image(fmdl,1); vh = fwd_solve(img);
0304             img.elem_data(fmdl.demo_targ_elems.A) = 1.2;
0305             img.elem_data(fmdl.demo_targ_elems.B) = 0.8;
0306                                     vi = fwd_solve(img);
0307             
0308          otherwise; break
0309       end;
0310 
0311 %     disp([i,imdl0.hyperparameter.value]);
0312       if strcmp( imdl0.reconst_type, 'absolute')
0313          imgr = inv_solve( imdl0, vi);
0314       else
0315          imgr = inv_solve( imdl0, vh, vi);
0316       end
0317       subplot(4,5,i); show_slices( imgr, [inf,inf,0] );
0318    end