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