GREIT_errors

 GREIT_errors: Add Error Compensation to GREIT-type model
 imdl= imdl_errors(imdli, options )
  imdl - inv_model structure

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

 Available options are:

 'Channel:[list]' or 'Channel:[list] SNR=100'  # remove individual channels
 'Electrode:[list]' or 'Electrode:[list] SNR=100' # remove electrodes
 'RecipError:' or 'RecipError: Tau=1e-3'       # use calc_reciproc_error

 Data is a sequence of measured data to model electrode errors with

 See also: select_imdl

 Model must have a imdl.solve_use_matrix field with RM, M, PJt fields
  these are created with keep_model_components option to mk_GREIT_model

0001 function imdl= GREIT_errors(imdli, opts, data )
0002 % GREIT_errors: Add Error Compensation to GREIT-type model
0003 % imdl= imdl_errors(imdli, options )
0004 %  imdl - inv_model structure
0005 %
0006 % OPTIONS => {'opt1','opt2'} options are processed in the order specified
0007 % if OPTIONS is char, treated as first option
0008 %
0009 % Available options are:
0010 %
0011 % 'Channel:[list]' or 'Channel:[list] SNR=100'  # remove individual channels
0012 % 'Electrode:[list]' or 'Electrode:[list] SNR=100' # remove electrodes
0013 % 'RecipError:' or 'RecipError: Tau=1e-3'       # use calc_reciproc_error
0014 %
0015 % Data is a sequence of measured data to model electrode errors with
0016 %
0017 % See also: select_imdl
0018 %
0019 % Model must have a imdl.solve_use_matrix field with RM, M, PJt fields
0020 %  these are created with keep_model_components option to mk_GREIT_model
0021 
0022 % (C) 2024 Andy Adler. License: GPL version 2 or version 3
0023 % $Id: GREIT_errors.m 6682 2024-03-14 19:46:10Z aadler $
0024 
0025 if ischar(imdli) && strcmp(imdli,'UNIT_TEST'); do_unit_test; return; end
0026 if nargin<3; data=[]; end
0027 
0028 imdl = imdli;
0029 if ischar(opts); opts = {opts}; end
0030 for i=1:length(opts);
0031    if iscell(opts); opti = opts{i};
0032    else           ; opti = opts; end
0033    imdl = process(imdl, opti, data);
0034 end
0035 
0036 function imdl = process(imdli, opti, data);
0037    if ~ischar(opti); error('Expecting string options'); end
0038    SNR = 1000; % Default values
0039    Tau = -3e-4;
0040    imdl = imdli;
0041    PJt        = imdli.solve_use_matrix.PJt;
0042    M          = imdli.solve_use_matrix.M;
0043    noiselev   = imdli.solve_use_matrix.noiselev;
0044 
0045    % split the option on an equals sign
0046    opt= regexp(opti,'(.[^:]*):?(.\S*)\s*(.*)','tokens'); opt= opt{1};
0047    matches = regexp(opt{3},'(\S+)=(\S+)','tokens');
0048    for i=1:length(matches); switch matches{i}{1}
0049       case 'SNR'; SNR = str2num(matches{i}{2});
0050       case 'Tau'; Tau =-abs(str2num(matches{i}{2})) 
0051       otherwise;  error('unknown parameter');
0052    end; end
0053 
0054    switch lower(opt{1});
0055       case 'channel';
0056          noise_covar = sparse(size(M,1),size(M,2));
0057          for ch= str2num(opt{2});
0058             noise_covar(ch,ch) = SNR;
0059          end
0060       case 'electrode';
0061          imdli.meas_icov_rm_elecs.elec_list = str2num(opt{2});
0062          imdli.meas_icov_rm_elecs.exponent  = -1;
0063          imdli.meas_icov_rm_elecs.SNR       = SNR;
0064          noise_covar = meas_icov_rm_elecs(imdli);
0065       case 'reciperror';
0066          imdl.calc_reciproc_error.tau = Tau;
0067          noise_covar= calc_reciproc_error( imdl, data );
0068       otherwise error('option not understood');
0069    end
0070 
0071    M   = M + noiselev^2 * noise_covar;
0072    imdl.solve_use_matrix.RM = left_divide(M.',PJt.').';    %PJt/M;
0073    
0074 
0075 function do_unit_test
0076    for sp=5:9; switch sp
0077       case 1; [vv,imdl] = test_data(1);
0078       case 2; [vv,imdl] = test_data(2);
0079       case 3; [vv,imdl] = test_data(2);
0080          imdl = GREIT_errors(imdl,'Electrode:[2]');
0081       case 4; [vv,imdl] = test_data(2);
0082          imdl = GREIT_errors(imdl,'Channel:[3]');
0083       case 5; [vv,imdl] = test_data(2);
0084          imdl = GREIT_errors(imdl,'Channel:[3] SNR=10000');
0085       case 6; [vv,imdl] = test_data(2);
0086          imdl = GREIT_errors(imdl,'RecipError: Tau=.0001',vv);
0087       case 7; [vv,imdl] = test_data(2);
0088          imdl = GREIT_errors(imdl,'RecipError:',vv);
0089       otherwise; break;
0090   
0091       end
0092       imgr = inv_solve(imdl,vv(:,1),vv(:,end));
0093       imgr.calc_colours.ref_level = 0;
0094       subplot(3,3,sp); hh=show_fem(imgr);
0095       set(hh,'EdgeColor','none');
0096    end
0097 
0098 function [vv,imdl] = test_data(mode);
0099    fmdl = mk_library_model('adult_male_16el_lungs');
0100    fmdl.stimulation = mk_stim_patterns(16,1,[0,1],[0,1],{});
0101    img = mk_image(fmdl,1);
0102    k=0;for sig = linspace(0.2,0.1,50); k=k+1;
0103       img.elem_data(vertcat(fmdl.mat_idx{2:3}))=sig;
0104       vv(:,k) = getfield(fwd_solve(img),'meas');
0105    end
0106    imdl = select_imdl(fmdl,{'GREIT:NF=0.3 64x64 rad=0.25'});
0107 
0108    rng('default');
0109    switch mode
0110       case 1;
0111       case 2; vv(3,:) = 0.1*randn(size( vv(3,:) ));
0112       otherwise; error 'huh?'
0113    end