Home > eidors > solvers > inv_solve.m

inv_solve

PURPOSE ^

INV_SOLVE: calculate imag from an inv_model and data

SYNOPSIS ^

function img = inv_solve( inv_model, data1, data2)

DESCRIPTION ^

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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

Generated on Fri 01-Jun-2018 15:59:55 by m2html © 2005