inv_solve

 INV_SOLVE: calculate imag from an inv_model and data
 
 inv_solve can be called as
     img= inv_solve( inv_model, data1, data2)
   if inv_model.reconst_type = 'difference'
 or
     img= inv_solve( inv_model, data )
   if inv_model.reconst_type = 'static'

   if inv_model.reconst_to = 'nodes' then output
      img.node_data has output data
   else   reconst_to = 'elems' (DEFAULT) then output to
      img.elem_data has output data

 in each case it will call the inv_model.solve

 data      is a measurement data structure
 inv_model is a inv_model structure
 img       is an image structure
           or a vector of images is data1 or data2 are vectors

 For difference EIT:
 data1      => difference data at earlier time (ie homogeneous)
 data2      => difference data at later time   (ie inhomogeneous)

 data can be:
   - an EIDORS data object

   - an M x S matrix, where M is the total number
         of measurements expected by inv_model

   - an M x S matrix, where M is n_elec^2
        when not using data from current injection
        electrodes, it is common to be given a full
        measurement set.  For example, 16 electrodes give
        208 measures, but 256 measure sets are common.
        Data will be selected based on fwd_model.meas_select.

 If S > 1 for both data1 and data2 then the matrix sizes must be equal

 Parameters:
   inv_model.inv_solve.select_parameters: indices of parameters to return
                         DEFAULT: return all paramteres
  Scale solution (to correct for amplitude or other defects)
   inv_model.inv_solve.scale_solution.offset
   inv_model.inv_solve.scale_solution.scale
  Disable solution error calculations
   inv_model.inv_solve.calc_solution_error = 0

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