demo_real_test3

 Perform tests based on the demo_real function with new structs

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