Description of inv

0001 function img = inv_solve( inv_model, data1, data2)
0002 % INV_SOLVE: calculate imag from an inv_model and data
0003 %
0004 % inv_solve can be called as
0005 %     img= inv_solve( inv_model, data1, data2)
0006 %   if inv_model.reconst_type = 'difference'
0007 % or
0008 %     img= inv_solve( inv_model, data )
0009 %   if inv_model.reconst_type = 'static'
0010 %
0011 %   if inv_model.reconst_to = 'nodes' then output
0012 %      img.node_data has output data
0013 %   else   reconst_to = 'elems' (DEFAULT) then output to
0014 %      img.elem_data has output data
0015 %
0016 % in each case it will call the inv_model.solve
0017 %
0018 % data      is a measurement data structure
0019 % inv_model is a inv_model structure
0020 % img       is an image structure
0021 %           or a vector of images is data1 or data2 are vectors
0022 %
0023 % For difference EIT:
0024 % data1      => difference data at earlier time (ie homogeneous)
0025 % data2      => difference data at later time   (ie inhomogeneous)
0026 %
0027 % data can be:
0028 %   - an EIDORS data object
0029 %
0030 %   - an M x S matrix, where M is the total number
0031 %         of measurements expected by inv_model
0032 %
0033 %   - an M x S matrix, where M is n_elec^2
0034 %        when not using data from current injection
0035 %        electrodes, it is common to be given a full
0036 %        measurement set.  For example, 16 electrodes give
0037 %        208 measures, but 256 measure sets are common.
0038 %        Data will be selected based on fwd_model.meas_select.
0039 %
0040 % If S > 1 for both data1 and data2 then the matrix sizes must be equal
0041 %
0042 % Parameters:
0043 %   inv_model.inv_solve.select_parameters: indices of parameters to return
0044 %                         DEFAULT: return all paramteres
0045 %  Scale solution (to correct for amplitude or other defects)
0046 %   inv_model.inv_solve.scale_solution.offset
0047 %   inv_model.inv_solve.scale_solution.scale
0048 %  Disable solution error calculations
0049 %   inv_model.inv_solve.calc_solution_error = 0
0050 %
0051 % img.time is frame sequence number. It is in seconds if
0052 %   imdl.fwd_model.frame_rate is specified
0053 
0054 
0055 % (C) 2005-2010 Andy Adler. License: GPL version 2 or version 3
0056 % $Id: inv_solve.m 6665 2024-03-08 18:41:02Z aadler $
0057 
0058 if ischar(inv_model) && strcmp(inv_model,'UNIT_TEST'); do_unit_test; return; end
0059 
0060 inv_model = prepare_model( inv_model );
0061 opts = parse_parameters( inv_model );
0062 
0063 print_hello(inv_model.solve);
0064 try inv_model.solve = str2func(inv_model.solve); end
0065 
0066 if opts.abs_solve
0067    if nargin~=2;
0068       error('only one data set is allowed for a static reconstruction');
0069    end
0070    
0071    fdata = filt_data(inv_model,data1);
0072 
0073    imgc= eidors_cache( inv_model.solve, {inv_model, fdata}, 'inv_solve');
0074 
0075 else
0076    if nargin~=3;
0077       error('two data sets are required for a difference reconstruction');
0078    end
0079 
0080    % expand data sets if one is provided that is longer
0081    data_width= max(num_frames(data1), num_frames(data2));
0082 
0083    fdata1 = filt_data( inv_model, data1, data_width );
0084    fdata2 = filt_data( inv_model, data2, data_width );
0085    % TODO: Check if solver can handle being called with multiple data
0086    imgc= eidors_cache( inv_model.solve, {inv_model, fdata1, fdata2},'inv_solve');
0087    
0088 
0089 end
0090 
0091 check_parametrization_handling(inv_model,imgc);
0092      
0093 img = eidors_obj('image', imgc );
0094 % img = data_mapper(img,1);
0095 if ~isfield(img,'current_params') % current_params is Bartek's units conversion control for img.elem_data.
0096 
0097   img.current_params = [];
0098 end
0099 
0100 % If we reconstruct with a different 'rec_model' then
0101 %  put this into the img
0102 if isfield(inv_model,'rec_model')
0103    img.fwd_model= inv_model.rec_model;
0104 end
0105 
0106 
0107 if opts.select_parameters;
0108    img.elem_data = img.elem_data( opts.select_parameters, :);
0109 end
0110 
0111 if isfield(img,'elem_data')
0112    img.time = 1:size(img.elem_data,2);
0113 else
0114    img.time = 1:size(img.node_data,2);
0115 end
0116 
0117 if isfield(inv_model.fwd_model,'frame_rate');
0118    img.time = img.time / inv_model.fwd_model.frame_rate;
0119 end
0120 
0121 if ~opts.reconst_to_elems
0122   img.node_data= img.elem_data;
0123   img = rmfield(img,'elem_data');
0124 end
0125 
0126 % Scale if required
0127 try; img.elem_data = opts.offset + opts.scale * img.elem_data; end
0128 try; img.node_data = opts.offset + opts.scale * img.node_data; end
0129 
0130 % MATLAB IS SUPER SUPER STUPID HERE. YOU CAN'T ASSIGN IF FIELDS ARE IN
0131 %  A DIFFERENT ORDER. Example
0132 %>> A.a= 1; A.b= 2; B.b= 3; B.a = 3; A(2) = B
0133 %??? Subscripted assignment between dissimilar structures.
0134 img.info.error = NaN;
0135 img = orderfields(img);
0136 
0137 % calculate residuals
0138 try 
0139    do_calc = inv_model.inv_solve.calc_solution_error;
0140 catch
0141    eidors_msg('inv_solve: not calculting solution residual',3);
0142    do_calc = false;
0143 end
0144 if ~do_calc;
0145    return;
0146 end
0147 eidors_msg('inv_solve: Calculating solution residual (inv_model.inv_solve.calc_solution_error = 0 to disable)',2);
0148 try
0149    if opts.abs_solve
0150       img.info.error = calc_solution_error( imgc, inv_model, fdata);
0151    else
0152       img.info.error = calc_solution_error( imgc, inv_model, fdata1, fdata2 );
0153    end
0154    eidors_msg('inv_solve: Solution Error: %f', img.info.error,  2);
0155 catch
0156    eidors_msg('inv_solve: Solution Error calculation failed.',2);
0157 end
0158 
0159 function print_hello(solver)
0160     if isa(solver,'function handle')
0161         solver = func2str(solver);
0162     end
0163     if strcmp(solver, 'eidors_default');
0164         solver = eidors_default('get','inv_solve');
0165     end
0166     eidors_msg('inv_solve: %s', solver,3);
0167 
0168 
0169 function opts = parse_parameters( imdl )
0170    if  strcmp(imdl.reconst_type,'static') || ...
0171        strcmp(imdl.reconst_type,'absolute')
0172       opts.abs_solve = 1;
0173    elseif strcmp(imdl.reconst_type,'difference')
0174       opts.abs_solve = 0;
0175    else
0176       error('inv_model.reconst_type (%s) not understood', imdl.reconst_type); 
0177    end
0178 
0179    opts.select_parameters = [];
0180    try
0181       opts.select_parameters = imdl.inv_solve.select_parameters;
0182    end
0183 
0184    opts.reconst_to_elems = 1;
0185    try if strcmp( imdl.reconst_to, 'nodes' )
0186       opts.reconst_to_elems = 0;
0187    end; end
0188    
0189    opts.scale  = 1;
0190    try opts.scale = imdl.inv_solve.scale_solution.scale; end
0191 
0192    opts.offset = 0;
0193    try opts.offset = imdl.inv_solve.scale_solution.offset; end
0194 
0195 function imdl = deprecate_imdl_parameters(imdl)
0196    if isfield(imdl, 'parameters')
0197       if isstruct(imdl.parameters)
0198 %          disp(imdl)
0199 %          disp(imdl.parameters)
0200 %          warning('EIDORS:deprecatedParameters','INV_SOLVE inv_model.parameters.* are deprecated in favor of inv_model.inv_solve.* as of 30-Apr-2014.');
0201 
0202          if ~isfield(imdl, 'inv_solve')
0203             imdl.inv_solve = imdl.parameters;
0204          else % we merge
0205             % merge struct trick from:
0206             %  http://stackoverflow.com/questions/38645
0207             A = imdl.parameters;
0208             B = imdl.inv_solve;
0209             M = [fieldnames(A)' fieldnames(B)'; struct2cell(A)' struct2cell(B)'];
0210             try % assumes no collisions
0211                imdl.inv_solve=struct(M{:});
0212             catch % okay, collisions - do unique to resolve them
0213                [tmp, rows] = unique(M(1,:), 'last');
0214                M=M(:,rows);
0215                imdl.inv_solve=struct(M{:});
0216             end
0217          end
0218          imdl = rmfield(imdl, 'parameters');
0219          imdl.parameters = imdl.inv_solve; % backwards compatible!
0220       else
0221          error('unexpected inv_model.parameters where parameters is not a struct... i do not know what to do');
0222       end
0223    end
0224 
0225 function mdl = prepare_model( mdl )
0226     fmdl = mdl.fwd_model;
0227     fmdl = mdl_normalize(fmdl,mdl_normalize(fmdl));
0228     if ~isfield(fmdl,'elems');
0229         return;
0230     end
0231 
0232     fmdl.elems  = double(fmdl.elems);
0233     fmdl.n_elem = size(fmdl.elems,1);
0234     fmdl.n_node = size(fmdl.nodes,1);
0235     if isfield(fmdl,'electrode');
0236         fmdl.n_elec = length(fmdl.electrode);
0237     else
0238         fmdl.n_elec = 0;
0239     end
0240 
0241     mdl.fwd_model= fmdl;
0242     if ~isfield(mdl,'reconst_type');
0243         mdl.reconst_type= 'difference';
0244     end
0245 
0246     % warn if we have deprecated inv_model.parameters laying about
0247     mdl = deprecate_imdl_parameters(mdl);
0248 
0249 function check_parametrization_handling(inv_model,imgc)
0250 if isfield(inv_model, 'jacobian_bkgnd') && ... 
0251     has_params(inv_model.jacobian_bkgnd) && ~has_params(imgc)
0252    if isa(inv_model.solve,'function_handle')
0253       solver = func2str(inv_model.solve);
0254    else
0255       solver = inv_model.solve;
0256    end
0257    if strcmp(solver,'eidors_default');
0258       solver = eidors_default('get','inv_solve');
0259    end
0260    warning('EIDORS:ParamsObliviousSolver',...
0261       ['The solver %s did not handle the parametrization data properly.\n'...
0262        'The results may be incorrect. Please check the code to verify.'], ...
0263        solver);
0264 end
0265          
0266     
0267 function b = has_params(s)
0268 b = false;
0269 if isstruct(s)
0270    b = any(ismember(fieldnames(s),supported_params));
0271 end
0272 
0273 
0274 % TODO: this code really needs to be cleaned, but not before eidors 3.4
0275 function nf= num_frames(d0)
0276    if isnumeric( d0 )
0277       nf= size(d0,2);
0278    elseif d0(1).type == 'data';
0279       nf= size( horzcat( d0(:).meas ), 2);
0280    else
0281       error('Problem calculating number of frames. Expecting numeric or data object');
0282    end
0283    
0284 % test for existance of meas_select and filter data
0285 function d2= filt_data(inv_model, d0, data_width )
0286    if ~isnumeric( d0 )
0287        % we probably have a 'data' object
0288 
0289        d1 = [];
0290        for i=1:length(d0)
0291           if strcmp( d0(i).type, 'data' )
0292               d1 = [d1, d0(i).meas];
0293           else
0294               error('expecting an object of type data');
0295           end
0296        end
0297 
0298    else
0299       % we have a matrix of data. Hope for the best
0300       d1 = d0;
0301    end
0302 
0303    d1= double(d1); % ensure we can do math on our object
0304 
0305    if isfield(inv_model.fwd_model,'meas_select') && ...
0306      ~isempty(inv_model.fwd_model.meas_select)
0307       % we have a meas_select parameter that isn []
0308 
0309       meas_select= inv_model.fwd_model.meas_select;
0310       if     size(d1,1) == length(meas_select)
0311          d2= d1(meas_select,:);
0312       elseif size(d1,1) == sum(meas_select==1)
0313          d2= d1;
0314       else
0315          error('inconsistent difference data: (%d ~= %d). Maybe check fwd_model.meas_select',  ...
0316                size(d1,1), length(meas_select));
0317       end
0318    else
0319       d2= d1;
0320    end
0321 
0322    if nargin==3 % expand to data width
0323       d2_width= size(d2,2);
0324       if d2_width == data_width
0325          % ok
0326       elseif d2_width == 1
0327          d2= d2(:,ones(1,data_width));
0328       else
0329          error('inconsistent difference data: (%d ~= %d)',  ...
0330                d2_width, data_width);
0331       end
0332    end
0333 
0334 % Test code
0335 function do_unit_test
0336    k=0; N=5; nd = 5;
0337 
0338    imdl = mk_common_model('d2c2',16);
0339    imdl = select_imdl( imdl, {'Choose NF=2.5'});
0340    mvx = linspace(-0.8,0.8,nd);
0341    [vh,vi] = simulate_movement(mk_image(imdl), [mvx;0*mvx;0.05+0*mvx]);
0342 
0343    img= inv_solve(imdl,vh,vi); img.show_slices.img_cols = 5;
0344    k=k+1; subplot(N,1,k); show_slices(img);
0345    unit_test_cmp('inv_solve: 1a', mean(img.elem_data), 3e-3, 1e-3);
0346    unit_test_cmp('inv_solve: 1b',  std(img.elem_data), 1.5e-2, 1e-2);
0347 
0348    vhm= eidors_obj('data','nn','meas',vh);
0349    img= inv_solve(imdl,vhm,vi);
0350    unit_test_cmp('inv_solve: 2a', mean(img.elem_data), 3e-3, 1e-3);
0351    unit_test_cmp('inv_solve: 2b',  std(img.elem_data), 1.5e-2, 1e-2);
0352 
0353    img= inv_solve(imdl,vh*ones(1,nd),vi);
0354    unit_test_cmp('inv_solve: 3a', mean(img.elem_data), 3e-3, 1e-3);
0355    unit_test_cmp('inv_solve: 3b',  std(img.elem_data), 1.5e-2, 1e-2);
0356 
0357    vim= eidors_obj('data','nn','meas',vi);
0358    img= inv_solve(imdl,vhm,vim);
0359    unit_test_cmp('inv_solve: 4a', mean(img.elem_data), 3e-3, 1e-3);
0360    unit_test_cmp('inv_solve: 4b',  std(img.elem_data), 1.5e-2, 1e-2);
0361 
0362    vhm= eidors_obj('data','nn','meas',vh*ones(1,nd));
0363    img= inv_solve(imdl,vhm,vi);
0364    unit_test_cmp('inv_solve: 5a', mean(img.elem_data), 3e-3, 1e-3);
0365    unit_test_cmp('inv_solve: 5b',  std(img.elem_data), 1.5e-2, 1e-2);
0366 
0367    vhm(1)= eidors_obj('data','nn','meas',vh*ones(1,2));
0368    vhm(2)= eidors_obj('data','nn','meas',vh*ones(1,nd-2));
0369    img= inv_solve(imdl,vhm,vi);
0370    unit_test_cmp('inv_solve: 6a', mean(img.elem_data), 3e-3, 1e-3);
0371    unit_test_cmp('inv_solve: 6b',  std(img.elem_data), 1.5e-2, 1e-2);
0372 
0373    vim(1)= eidors_obj('data','nn','meas',vi(:,1:3));
0374    vim(2)= eidors_obj('data','nn','meas',vi(:,4:end));
0375    img= inv_solve(imdl,vhm,vim);
0376    unit_test_cmp('inv_solve: 7a', mean(img.elem_data), 3e-3, 1e-3);
0377    unit_test_cmp('inv_solve: 7b',  std(img.elem_data), 1.5e-2, 1e-2);
0378 
0379    im2 = imdl; im2.inv_solve.select_parameters = 1:5;
0380    img= inv_solve(im2,vh,vi);
0381    unit_test_cmp('inv_solve: 10', size(img.elem_data), [5,nd]);
0382 
0383 
0384    im2 = imdl;
0385    im2.inv_solve.scale_solution.offset = 1;
0386    im2.inv_solve.scale_solution.scale = 2;
0387    img= inv_solve(im2,vh,vi); img.show_slices.img_cols = 5;
0388    unit_test_cmp('inv_solve: 20a', mean(img.elem_data), 1.006, 2e-3);
0389    unit_test_cmp('inv_solve: 20b',  std(img.elem_data), 3e-2, 1e-2);
0390    im2.inv_solve.calc_solution_error = 1;
0391    img= inv_solve(im2,vh,vi);
0392    unit_test_cmp('inv_solve: 20e', mean(img.elem_data), 1.006, 2e-3);
0393    
0394    im2.inv_solve.scale_solution.offset = 0;
0395    d = interp_mesh( imdl.fwd_model); d= sqrt(sum(d.^2,2));
0396    im2.inv_solve.scale_solution.scale = spdiags(1-d,0,length(d),length(d));
0397    img= inv_solve(imdl,vh,vi); img.show_slices.img_cols = 5;
0398    k=k+1; subplot(N,1,k); show_slices(img);
0399 
0400    im2 = select_imdl(imdl, {'Nodal GN dif'} );
0401    img= inv_solve(im2,vh,vi);
0402    unit_test_cmp('inv_solve: 30a', mean(img.node_data), 3e-3, 1e-3);
0403    unit_test_cmp('inv_solve: 30b',  std(img.node_data), 1.5e-2, 1e-2);
0404 
0405    im2.inv_solve.scale_solution.offset = 1;
0406    im2.inv_solve.scale_solution.scale = 2;
0407    img= inv_solve(im2,vh,vi); img.show_slices.img_cols = 5;
0408    unit_test_cmp('inv_solve: 31a', mean(img.node_data), 1.006, 2e-3);
0409    unit_test_cmp('inv_solve: 31b',  std(img.node_data), 3e-2, 1.5e-2);
0410 
0411    im2 = select_imdl(imdl, {'Basic GN abs'} );
0412    im2.inv_solve.calc_solution_error = 0; % ensure accepted
0413    img= inv_solve(im2,vi(:,1));
0414    unit_test_cmp('inv_solve: 40a', mean(img.elem_data), 1.004, 1e-3);
0415    unit_test_cmp('inv_solve: 40b',  std(img.elem_data), 1.5e-2, 1e-2);
0416    im2.inv_solve.calc_solution_error = 1;
0417    img= inv_solve(im2,vi(:,1));
0418    unit_test_cmp('inv_solve: 40e', mean(img.elem_data), 1.004, 1e-3);
0419 
0420
inv_solve

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
