calc_solution_error

CALC_SOLUTION_ERROR Calculate residuals for a solution
 E = CALC_SOLUTION_ERROR(SOL, IMDL, DATA) calculates the residual error
 where 
     SOL   -  EIDORS 'image' structure containing the solution
     IMDL  -  EIDORS 'inv_model' structure used to calculate the solution
     DATA  -  the data to be fitted (either a vector, or EIDORS 'data'
              struct)
 and 
     E = norm(DATA - SOL_MEAS) / norm( DATA ),
 where SOL_MEAS is the simulated data obtained from SOL.

 E = CALC_SOLUTION_ERROR(SOL, IMDL, VH, VI) allows specifying two data 
 inputs for difference solvers, such that DATA = VI - VH

 [E RES] = CALC_SOLUTION_ERROR(...) also returns the vector of residuals RES

 See also INV_SOLVE, CALC_DIFFERENCE_DATA

0001 function [e res] = calc_solution_error(imgc, imdl, vh, vi)
0002 %CALC_SOLUTION_ERROR Calculate residuals for a solution
0003 % E = CALC_SOLUTION_ERROR(SOL, IMDL, DATA) calculates the residual error
0004 % where
0005 %     SOL   -  EIDORS 'image' structure containing the solution
0006 %     IMDL  -  EIDORS 'inv_model' structure used to calculate the solution
0007 %     DATA  -  the data to be fitted (either a vector, or EIDORS 'data'
0008 %              struct)
0009 % and
0010 %     E = norm(DATA - SOL_MEAS) / norm( DATA ),
0011 % where SOL_MEAS is the simulated data obtained from SOL.
0012 %
0013 % E = CALC_SOLUTION_ERROR(SOL, IMDL, VH, VI) allows specifying two data
0014 % inputs for difference solvers, such that DATA = VI - VH
0015 %
0016 % [E RES] = CALC_SOLUTION_ERROR(...) also returns the vector of residuals RES
0017 %
0018 % See also INV_SOLVE, CALC_DIFFERENCE_DATA
0019 
0020 % (C) 2013 Bartlomiej Grychtol, License: GPL version 2 or 3.
0021 % $Id: calc_solution_error.m 4733 2015-03-23 13:29:22Z aadler $
0022 
0023 % TODO: Generalize coarse2fine
0024 
0025 
0026 switch imdl.reconst_type
0027    case 'difference'
0028       if nargin==4
0029          data = calc_difference_data(vh,vi,imdl.fwd_model);
0030       else
0031          if isstruct(vh)
0032             data = vh.meas;   % eidors object
0033          else
0034             data = vh;        % vector of measurements
0035          end
0036       end
0037       res = calc_diff_residual(imgc,imdl,data);
0038    case {'absolute', 'static'}
0039       if isstruct(vh)
0040          data = vh.meas;
0041       else
0042          data = vh;
0043       end
0044       res = calc_abs_residual(imgc,imdl,data);
0045    otherwise
0046       error('reconst_type not recognized');
0047 end
0048 
0049 
0050 % avarage error
0051 e = norm(res)/norm(data);
0052 
0053 
0054 function res = calc_abs_residual(imgc,imdl,data)
0055 fmdl = imdl.fwd_model;
0056 
0057 % make sure to have elem_data irrespective of parametrization
0058 img = calc_jacobian_bkgnd(imdl); 
0059 img = data_mapper(img);
0060 imgc = data_mapper(imgc);
0061 
0062 % account for coarse2fine
0063 if size(img.elem_data,1) == size(imgc.elem_data,1)
0064    img.elem_data = imgc.elem_data;
0065 else
0066    img.elem_data = fmdl.coarse2fine*imgc.elem_data;
0067 end
0068 img = data_mapper(img);
0069 
0070 % simualate data from solution
0071 sim = fwd_solve(img);
0072 
0073 % residuals
0074 res = sim.meas - data;
0075 
0076 
0077 
0078 function res = calc_diff_residual(imgc,imdl,data)
0079 
0080 fmdl = imdl.fwd_model;
0081 img = calc_jacobian_bkgnd(imdl);
0082 simh = fwd_solve(img);
0083 
0084 % map parametrization to elem_data
0085 img = data_mapper(img);
0086 imgc = data_mapper(imgc);
0087 
0088 if ~isfield(imgc,'elem_data') && isfield(imgc,'node_data')
0089    eidors_msg('Solution error calculation for nodal solvers not supported (yet).',2);
0090    res = NaN;
0091    return
0092 end
0093    
0094 % add solution to jacobian background
0095 e_data = repmat(img.elem_data,1,size(imgc.elem_data,2));
0096 has_c2f = isfield(imgc.fwd_model,'coarse2fine');
0097 if ~has_c2f
0098    n_elems = num_elems(imgc.fwd_model);
0099    img.elem_data = e_data + imgc.elem_data(1:n_elems,:);
0100 else
0101    n_elems = size(imgc.fwd_model.coarse2fine,2);
0102    img.elem_data = e_data + fmdl.coarse2fine*imgc.elem_data(1:n_elems,:);
0103 end
0104 img = data_mapper(img,1);
0105 
0106 % simulate data from solution
0107 jnk = img;
0108 for i = fliplr(1:size(imgc.elem_data,2));
0109    jnk.elem_data = img.elem_data(:,i);
0110    tmp = fwd_solve(jnk);
0111    simi(:,i) = tmp.meas;
0112 end
0113 sim = calc_difference_data(simh,simi,fmdl);
0114 
0115 % residuals
0116 res = data - sim;