demo_real_test2

 Perform tests based on the demo_real function with new structs
 $Id: demo_real_test2.m 3630 2012-11-18 18:29:52Z aadler $

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