calc_data_prior_test

 Verify dataprior:
   for  difference EIT: dataprior should be 1
   for  normalized EIT: dataprior should be 1 / homogeneous
 $Id: calc_data_prior_test.m 3097 2012-06-08 14:07:14Z bgrychtol $

0001 function ok= calc_data_prior_test
0002 % Verify dataprior:
0003 %   for  difference EIT: dataprior should be 1
0004 %   for  normalized EIT: dataprior should be 1 / homogeneous
0005 % $Id: calc_data_prior_test.m 3097 2012-06-08 14:07:14Z bgrychtol $
0006 
0007 ok= 1;
0008 
0009 % Create simple 16 electrode 2D model
0010 %
0011 n_elec= 16;
0012 n_rings= 1;
0013  options = {'no_meas_current','no_rotate_meas'};
0014 params= mk_circ_tank(8, [], n_elec);
0015 
0016 params.stimulation= mk_stim_patterns(n_elec, n_rings, '{ad}','{ad}', ...
0017                             options, 10);
0018 params.solve=      'fwd_solve_1st_order';
0019 params.system_mat= 'system_mat_1st_order';
0020 params.jacobian=   'jacobian_adjoint';
0021 params.normalize_measurements = 0;
0022 mdl_2d = eidors_obj('fwd_model', params);
0023 
0024 % create homogeneous image + simulate data
0025 mat= ones( size(mdl_2d.elems,1) ,1);
0026 
0027 homg_img= eidors_obj('image', 'homg image', ...
0028                      'elem_data', mat, 'fwd_model', mdl_2d );
0029 
0030 homg_data=fwd_solve( homg_img);
0031 J= calc_jacobian( homg_img );
0032 
0033 sumdiff= 0;
0034 delta = 2e-5;
0035 testvec= [5,20,40,130];
0036 for testelem = testvec
0037     mat= ones( size(mdl_2d.elems,1) ,1);
0038     mat(testelem)= 1+delta;
0039     inh_img= eidors_obj('image', 'inh', ...
0040                          'elem_data', mat, 'fwd_model', mdl_2d );
0041     inh_data=fwd_solve( inh_img);
0042 
0043     simJ= 1/delta* (inh_data.meas - homg_data.meas);
0044     
0045 %   plot([J(:,testelem) simJ]);
0046     sumdiff = sumdiff + std( J(:,testelem) - simJ );
0047 end
0048 
0049 tol= 1e-4*std(J(:));
0050 if sumdiff/length(testvec) > tol
0051    error('Jacobian calculation error');
0052    ok=0;
0053 end
0054 
0055 
0056 % create normalized model
0057 params.normalize_measurements = 1;
0058 mdl_2d_norm = eidors_obj('fwd_model', params);
0059 % create homogeneous image with normalize_measurements
0060 mat= ones( size(mdl_2d.elems,1) ,1);
0061 homg_img= eidors_obj('image', 'homg image', ...
0062                      'elem_data', mat, 'fwd_model', mdl_2d_norm );
0063 
0064 J= calc_jacobian( homg_img );
0065 sumdiff= 0;
0066 delta = 2e-5;
0067 testvec= [5,20,40,130];
0068 for testelem = testvec
0069     mat= ones( size(mdl_2d.elems,1) ,1);
0070     mat(testelem)= 1+delta;
0071     inh_img= eidors_obj('image', 'inh', ...
0072                          'elem_data', mat, 'fwd_model', mdl_2d );
0073     inh_data=fwd_solve( inh_img);
0074 
0075     simJ= 1/delta* (inh_data.meas ./ homg_data.meas - 1);
0076     
0077 %   plot([J(:,testelem) simJ]);
0078     sumdiff = sumdiff + std( J(:,testelem) - simJ );
0079 end
0080 
0081 tol= 1e-4*std(J(:));
0082 if sumdiff/length(testvec) > tol
0083    error('normalize Jacobian calculation error');
0084    ok=0;
0085 end