demo_real_test2

 Perform tests based on the demo_real function with new structs
 $Id: demo_real_test2.m 6689 2024-03-19 16:04:23Z bgrychtol $

0001 function ok= demo_real_test2
0002 % Perform tests based on the demo_real function with new structs
0003 % $Id: demo_real_test2.m 6689 2024-03-19 16:04:23Z bgrychtol $
0004 
0005 isOctave= exist('OCTAVE_VERSION');
0006 
0007 datareal= 'datareal.mat';
0008 datacom=  'datacom.mat';
0009 drt=      'demo_real_test.mat';
0010 
0011 % create FEM model structure
0012 
0013 load(datareal,'vtx','simp');
0014 
0015 demo_mdl.name = 'demo real model';
0016 demo_mdl.nodes= vtx;
0017 demo_mdl.elems= simp;
0018 demo_mdl.boundary= find_boundary( simp );
0019 demo_mdl.solve=      'np_fwd_solve';
0020 demo_mdl.jacobian=   'np_calc_jacobian';
0021 demo_mdl.system_mat= 'np_calc_system_mat';
0022 
0023 clear vtx simp
0024 
0025 % create FEM model electrodes definitions
0026 
0027 load(datareal,'gnd_ind','elec','zc','protocol','no_pl');
0028 perm_sym= '{y}';
0029 
0030 demo_mdl.gnd_node= gnd_ind;
0031 for i=1:length(zc)
0032     demo_mdl.electrode(i).z_contact= zc(i);
0033     demo_mdl.electrode(i).nodes=     elec(i,:);
0034 end
0035 
0036 % TODO: generalize the way that protocol sym no_pl are managed
0037 demo_mdl.np_fwd_solve.perm_sym     = perm_sym;
0038 
0039 % create FEM model stimulation and measurement patterns
0040 
0041 % get the current stimulation patterns
0042 [I,Ib] = set_3d_currents(protocol, ...
0043                          elec, ...
0044                          demo_mdl.nodes, ...
0045                          demo_mdl.gnd_node, ...
0046                          no_pl);
0047 % get the measurement patterns, only indH is used in this model
0048 %   here we only want to get the meas pattern from 'get_3d_meas',
0049 %   not the voltages, so we enter zeros
0050 [jnk,jnk,indH,indV,jnk] = get_3d_meas( ...
0051                   elec, demo_mdl.nodes, ...
0052                   zeros(size(I)), ... % Vfwd
0053                   Ib, no_pl );
0054 n_elec= size(elec,1);
0055 n_meas= size(indH,1) / size(Ib,2);
0056 for i=1:size(Ib,2)
0057     demo_mdl.stimulation(i).stimulation= 'Amp';
0058     demo_mdl.stimulation(i).stim_pattern= Ib(:,i);
0059     idx= ( 1+ (i-1)*n_meas ):( i*n_meas );
0060     meas_pat = sparse( (1:n_meas)'*[1,1], ...
0061                        indH( idx, : ), ...
0062                        ones(n_meas,2)*[1,0;0,-1], ...
0063                        n_meas, n_elec );
0064     demo_mdl.stimulation(i).meas_pattern= meas_pat;
0065 end
0066 
0067 clear gnd_ind elec zc protocol no_pl I Ib
0068 clear indH indV indH_sz meas_pat idx jnk
0069 
0070 demo_mdl= eidors_obj('fwd_model', demo_mdl);
0071 
0072 % simulate data for homogeneous medium
0073 
0074 homg_img.name = 'homogeneous image';
0075 homg_img.elem_data= ones( size(demo_mdl.elems,1) ,1);
0076 homg_img.fwd_model= demo_mdl;
0077 
0078 homg_img = eidors_obj('image', homg_img);
0079 
0080 homg_data=fwd_solve( demo_mdl, homg_img);
0081 
0082 % simulate data for inhomogeneous medium
0083 
0084 mat= ones( size(demo_mdl.elems,1) ,1);
0085 load( datacom ,'A','B') %Indices of the elements to represent the inhomogeneity
0086 mat(A)= mat(A)+0.15;
0087 mat(B)= mat(B)-0.20;
0088 
0089 inhomg_img.name = 'inhomogeneous image';
0090 inhomg_img.elem_data= mat;
0091 inhomg_img.fwd_model= demo_mdl;
0092 clear A B mat
0093 inhomg_img = eidors_obj('image', inhomg_img );
0094 
0095 inhomg_data=fwd_solve( demo_mdl, inhomg_img);
0096 
0097 % create inverse model
0098 
0099 % create an inv_model structure of name 'demo_inv'
0100 demo_inv.name= 'Nick Polydorides EIT inverse';
0101 demo_inv.solve=       'np_inv_solve';
0102 demo_inv.hyperparameter.value= 1e-4;
0103 demo_inv.R_prior= 'np_calc_image_prior';
0104 demo_inv.np_calc_image_prior.parameters= [3 1]; % see iso_f_smooth: deg=1, w=1
0105 demo_inv.jacobian_bkgnd.value= 1;
0106 demo_inv.reconst_type= 'difference';
0107 demo_inv.fwd_model= demo_mdl;
0108 demo_inv= eidors_obj('inv_model', demo_inv);
0109 
0110 % solve inverse model
0111 
0112 demo_img= inv_solve( demo_inv, homg_data, inhomg_data);
0113 
0114 % verifications
0115 
0116 load(drt);
0117 
0118 compare_tol( drt.voltageH, inhomg_data.meas, 'voltageH' )
0119 compare_tol( drt.sol, demo_img.elem_data, 'sol' )
0120 
0121 J= calc_jacobian( demo_mdl, homg_img );
0122 Jcolsby100=J(:,1:100:size(J,2));
0123 compare_tol( drt.Jcolsby100, Jcolsby100, 'Jcolsby100' )
0124 
0125 %Diag_Reg_012= [diag(Reg,0),[diag(Reg,1);0],[diag(Reg,2);0;0]];
0126 %compare_tol( drt.Diag_Reg_012, Diag_Reg_012, 'Diag_Reg_012' )
0127 
0128 ok=1;
0129 
0130 
0131 function compare_tol( cmp1, cmp2, errtext )
0132 % compare matrices and give error if not equal
0133 fprintf(2,'testing parameter: %s ...\n',errtext);
0134 
0135 tol= 2e-4;
0136 
0137 vd= mean(mean( abs(cmp1 - cmp2) ));
0138 vs= mean(mean( abs(cmp1) + abs(cmp2) ));
0139 if vd/vs > tol
0140    eidors_msg( ...
0141      'parameter %s exceeds tolerance %g (=%g)', errtext, tol, vd/vs, 1 );
0142 end
0143