calc_move_jacobian

 CALC_JACOBIAN   Computes the Jacobian matrix for conductivity and
 electrode movement variables in 3D EIT.
 Args:     fwd_model - the EIDORS object forward model
            img_bkgd - the image background conductivity
 Returns:          J - the Jacobian matrix [Jc, Jm]

 WARNING: THIS CODE IS EXPERIMENTAL AND GIVES PROBLEMS
 SEE: Camille Gomez-Laberge, Andy Adler
 Direct EIT Jacobian calculations for conductivity change
  and electrode movement,  Physiol. Meas., 29:S89-S99, 2008

0001 function J = calc_move_jacobian(fwd_model, img_bkgd)
0002 % CALC_JACOBIAN   Computes the Jacobian matrix for conductivity and
0003 % electrode movement variables in 3D EIT.
0004 % Args:     fwd_model - the EIDORS object forward model
0005 %            img_bkgd - the image background conductivity
0006 % Returns:          J - the Jacobian matrix [Jc, Jm]
0007 %
0008 % WARNING: THIS CODE IS EXPERIMENTAL AND GIVES PROBLEMS
0009 % SEE: Camille Gomez-Laberge, Andy Adler
0010 % Direct EIT Jacobian calculations for conductivity change
0011 %  and electrode movement,  Physiol. Meas., 29:S89-S99, 2008
0012 
0013 % (C) 2007, Camille Gomez-Laberge and Andy Adler.
0014 %  License: GPL version 2 or version 3
0015 % $Id: calc_move_jacobian.html 2819 2011-09-07 16:43:11Z aadler $
0016 
0017 warning('THIS CODE IS KNOWN TO HAVE BUGS - use with care');
0018 
0019 % System matrix and its parameters
0020 
0021 pp = aa_fwd_parameters( fwd_model );
0022 pp.dfact = factorial(pp.n_dims);
0023 pp.DEBUG = 1;
0024 if pp.DEBUG
0025     pp.ss_mat = calc_unconnected_system_mat( fwd_model, img_bkgd);
0026     pp.fwd_meas =fwd_solve( fwd_model, img_bkgd);
0027 end
0028 
0029 pp.Ce= connectivity_matrix( pp );
0030 s_mat= calc_system_mat( fwd_model, img_bkgd );
0031 [pp.Vc, pp.Re] = Vc_Re_matrices( pp, fwd_model, s_mat.E );
0032 Jc = calc_conductivity_jacobian(pp, fwd_model, img_bkgd);
0033 Jm = calc_movement_jacobian(pp, fwd_model, img_bkgd);
0034 J=[Jc,Jm];
0035 
0036 if pp.normalize
0037     data= fwd_solve( img_bkgd );
0038     J= J ./ (data.meas(:)*ones(1,size(J,2)));
0039 end
0040 
0041 
0042 
0043 % Define the element connectivity matrix Ce
0044 function Ce= connectivity_matrix( pp );
0045 lengthX = (pp.n_dims+1)*pp.n_elem;
0046 lengthY = pp.n_node;
0047 Xidx = pp.ELEM(:);
0048 Yidx = ones(lengthX, 1);
0049 Ce = sparse(1:lengthX, Xidx, Yidx, lengthX, lengthY);
0050 
0051 
0052 
0053 % Calculate fwd_solution and Impedance mapper matrices
0054 function [Vc, Re] = Vc_Re_matrices( pp, fwd_model, s_mat );
0055 % Define the stimulation matrix Vc for each node I, and pattern
0056 % Ground node is never excited; remove it from the index
0057 nodeidx = 1:pp.n_node;
0058 nodeidx( fwd_model.gnd_node ) = [];
0059 
0060 % The stimulation matrix Vc is the voltage at each node (row) for a
0061 % stimulation (column)
0062 Vc = zeros(pp.n_node, pp.n_stim);
0063 Vc(nodeidx, :) = s_mat(nodeidx, nodeidx) \ pp.QQ(nodeidx,:);
0064 
0065 % Define the electrode resistance matrix Re
0066 % Calculate the resistance between electrodes (row) and all nodes (col)
0067 % N2E matrix maps each electrode to its node(s); we exclude GND
0068 Re = zeros(pp.n_elec, pp.n_node);
0069 Re(:, nodeidx) = pp.N2E(:, nodeidx) / s_mat(nodeidx, nodeidx);
0070 
0071 % FIXME: why do we calculate the negative??
0072 Re = -Re;
0073 
0074 
0075 
0076 % Calculate Meas jacobian from derivative on nodes
0077 % Input delVc
0078 % Ouput J
0079 function J= nodes_to_stim_jacobian( delV, fwd_model, pp )
0080 
0081 sz_out= size(delV,3);
0082 % Define the conductivity Jacobian Jc
0083 J = zeros(pp.n_meas, sz_out);
0084 % Calculate the Jacobian columns for each stimulation pattern
0085 rowidx = 0;
0086 for j = 1:pp.n_stim
0087     % Get the measurement pattern for the stimulation pattern j
0088     meas_pattern = fwd_model.stimulation(j).meas_pattern;
0089     n_measj = size(meas_pattern, 1);
0090     % Extract the voltage sensitivity for electrode j
0091     delVcj = reshape( delV(:,j,:), pp.n_elec, sz_out);
0092     % Calculate sensitivity block for measurements during stimulation j
0093     J(rowidx+(1:n_measj), :) = meas_pattern*delVcj;
0094     rowidx = rowidx+n_measj;
0095 end
0096 
0097 
0098 
0099 % CONDUCTIVITY JACOBIAN (Based on Andy Adler's 1996 algorithms)
0100 % Define the voltage sensitivity delVc on electrode I, for stimulation
0101 % pattern J, for a change in conductivity of element K as a 3D array
0102 function Jc = calc_conductivity_jacobian(pp, fwd_model, img_bkgd);
0103 delVc = zeros(pp.n_elec, pp.n_stim, pp.n_elem);
0104 % Calculate the values for the voltage sensitivity for each element
0105 for k = 1:pp.n_elem
0106     if ~mod(k,500)
0107         fprintf('   JC: element # %d\n',k);
0108     end
0109     % Extract the coordinates of the element's four nodes
0110     Ae = pp.NODE(:,pp.ELEM(:,k))';
0111     % Augment Ae by adding a column of ones to invert
0112     Ae = inv([ones(pp.n_dims+1,1), Ae]);
0113     % Define Be as the matrix Ae with row 1 deleted
0114     Be = Ae(2:pp.n_dims+1,:);
0115     % Calculate the system submatrix subSe for the element i
0116     subSe = 2*Be'*Be/pp.dfact/abs(det(Ae));
0117     % Select the same submatrix of Ce
0118     subidx = (pp.n_dims+1)*(k-1)+1 : (pp.n_dims+1)*k;
0119     % The system submatrix is given by the product
0120     if ~pp.DEBUG
0121         delVc(:,:,k) = pp.Re * pp.Ce(subidx,:)' * subSe * pp.Ce(subidx,:)...
0122             * pp.Vc;
0123     else
0124         sz= (pp.n_dims+1)*pp.n_elem;
0125         delSe = sparse(sz,sz);
0126         se_idx= (pp.n_dims+1)*k+(-pp.n_dims : 0);
0127         delSe(se_idx, se_idx) = subSe;
0128 
0129         delVc(:,:,k) = pp.Re * pp.Ce' * delSe * pp.Ce * pp.Vc;
0130 
0131         if mod(k,50) == 0
0132             delta=1e-6;
0133             img_delta = img_bkgd;
0134             img_delta.elem_data(k) = img_delta.elem_data(k) + delta;
0135             ss_mat_delta= calc_unconnected_system_mat( fwd_model, img_delta);
0136             delSe_pert = (ss_mat_delta - pp.ss_mat) / delta;
0137 
0138             if norm(delSe -delSe_pert ,1) > 1e-6
0139                 eidors_msg('delSe calc wrong',1);
0140             end
0141         end
0142     end
0143 
0144 end
0145 Jc= nodes_to_stim_jacobian( delVc, fwd_model, pp );
0146 
0147 
0148 
0149 % MOVEMENT JACOBIAN
0150 function Jm = calc_movement_jacobian(pp, fwd_model, img_bkgd)
0151 % The movement Jacobian is defined for each coordinate Jm = [Jmx Jmy Jmz]
0152 % Define the voltage sensitivity delVm on electrode I, for stimulation
0153 % pattern J, for a movement of electrode K as a 3D array
0154 delVm = zeros(pp.n_elec, pp.n_stim, pp.n_elec*pp.n_dims);
0155 for colidx = 1:pp.n_dims
0156     fprintf('   JM: direction # %d\n',colidx);
0157     % Calculate the values for the voltage sensitivity for each electrode
0158     for k = 1:pp.n_elec
0159         % Find which elements touch this electrode
0160         elec_nodes = fwd_model.electrode(k).nodes;
0161  
0162         %for compound electrodes, average jacobian for each node
0163         delVm_part = zeros(pp.n_elec);
0164         for each_elec_node= elec_nodes(:)';
0165            delVm_part =  delVm_part + ...
0166                 calc_delVm(each_elec_node,pp,fwd_model,img_bkgd,colidx);
0167         end
0168         delVm_part = delVm_part/length(elec_nodes);
0169 
0170         vm_idx= k + pp.n_elec*(colidx-1);
0171         delVm(:,:,vm_idx) = delVm_part;
0172 
0173         if pp.DEBUG
0174             delta=1e-8;
0175             mdl_delta = fwd_model;
0176             mdl_delta.nodes(elec_nodes, colidx) = ...
0177                 mdl_delta.nodes(elec_nodes, colidx) + delta;
0178             S= calc_system_mat( mdl_delta, img_bkgd); S=S.E;
0179 keyboard
0180             [Vc_delta] = Vc_Re_matrices( pp, mdl_delta, S);
0181             delVm_pert = pp.N2E*(Vc_delta - pp.Vc) / delta;
0182             nn = norm(delVm_part - delVm_pert,1 ); % WHY NEGATIVE?
0183 
0184             if nn > 5e-5 ; keyboard; end
0185         end
0186     end
0187 end
0188 Jm= nodes_to_stim_jacobian( delVm, fwd_model, pp );
0189 
0190 
0191 
0192 function delVm=  calc_delVm( elec_nodes, pp, fwd_model, img_bkgd, colidx)
0193 [rowidx, elemidx] = find(pp.ELEM == elec_nodes);
0194 % Define the system sensitivity matrix to movement delSm
0195 sz= (pp.n_dims+1)*pp.n_elem;
0196 delSm = sparse(sz,sz);
0197 % For each touching element, calculate the perturbation
0198 jcount = 1;
0199 for j = elemidx'
0200     % Extract the coordinates of the element's four nodes
0201     Ae = pp.NODE(:,pp.ELEM(:, j))';
0202     % Define the invertible matrix P: augment Ae by adding a
0203     % column of ones to invert
0204     P = [ones(pp.n_dims+1,1), Ae];
0205     Ae = inv(P);
0206     absdetAe = abs(det(Ae));
0207     % Define Be as the matrix Ae with row 1 deleted
0208     Be = Ae(2:pp.n_dims+1,:);
0209     % For this coordinate, perturb P by [rowidx,colidx], which are
0210     % our paper's perturbation vectors [a,b]
0211     a = zeros(pp.n_dims+1,1);
0212     b = a;
0213     a(rowidx(jcount)) = 1;
0214     jcount = jcount + 1;
0215     b(colidx+1) = 1;
0216     % Calculate the system submatrix subSm for the element j by
0217     % asymmetric perturbation of the electrode node k
0218     deldetAe =   1/absdetAe*b'*Ae*a;
0219     delBe = -Ae*a*b'*Ae;
0220     delBe = delBe(2:pp.n_dims+1,:);
0221     subSm = 2/pp.dfact*(...
0222         deldetAe*Be'*Be + ...
0223         delBe'*Be/absdetAe + ...
0224         Be'*delBe/absdetAe);
0225 
0226     if pp.DEBUG
0227         delta=1e-8;
0228         subSe = 2*Be'*Be/pp.dfact/abs(det(Ae));
0229         d_NODE= pp.NODE;
0230         d_NODE(colidx,elec_nodes) =  d_NODE(colidx,elec_nodes) + delta;
0231         Ae = d_NODE(:,pp.ELEM(:, j))';
0232         Ae = inv( [ones(pp.n_dims+1,1), Ae] );
0233         absdetAe_pert = abs(det(Ae));
0234         deldetAe_pert = (absdetAe_pert - absdetAe) / delta;
0235         % Define Be as the matrix Ae with row 1 deleted
0236         Be = Ae(2:pp.n_dims+1,:);
0237         subSe_delta = 2*Be'*Be/pp.dfact/abs(det(Ae));
0238         subSm_pert= (subSe_delta - subSe ) / delta;
0239         if norm(subSm_pert - subSm,1) > 1e-5
0240             eidors_msg('subSm calc wrong',1);
0241             dd= (subSm_pert - 2/pp.dfact/absdetAe *...
0242                 (delBe'*Be + Be'*delBe) )./(Be'*Be);
0243 
0244             fprintf('colidx=%d, j=%d std=%6.4f >',...
0245                 colidx,j, std(dd(:)));
0246             keyboard
0247             subSm= subSm_pert;
0248         end
0249     end
0250 
0251     % Embed subSm into delSm such that subSm(1,1) is the
0252     % (4j+1,4j+1) element of delSm
0253     se_idx= (pp.n_dims+1)*j+(-pp.n_dims : 0);
0254     delSm(se_idx, se_idx) = subSm;
0255 end
0256 
0257 % The system submatrix is given by the product where delSm is
0258 % non-zero only in submatrices corresponding to touching elements
0259 delVm = pp.Re * pp.Ce' * delSm * pp.Ce * pp.Vc;
0260 if pp.DEBUG
0261     delta=1e-8;
0262     mdl_delta = fwd_model;
0263     mdl_delta.nodes(elec_nodes, colidx) = ...
0264         mdl_delta.nodes(elec_nodes, colidx) + delta;
0265     ss_mat_delta= calc_unconnected_system_mat( mdl_delta, img_bkgd );
0266     delSm_pert = (ss_mat_delta - pp.ss_mat) / delta;
0267     % delSe_pert shound be Ce'*delSe*Ce
0268     if norm(delSm -delSm_pert ,1) > 1e-5
0269         eidors_msg('delSm calc wrong',1);
0270         delVm = pp.Re * pp.Ce' * delSm_pert * pp.Ce * pp.Vc;
0271         keyboard
0272     end
0273 end
0274 
0275 
0276 
0277 function SS= calc_unconnected_system_mat( fwd_model, img)
0278 % Calc system matrix for Andy Adler's EIT code
0279 % fwd_model = forward model
0280 % img       = image background for system matrix calc
0281 % s_mat = CC' * SS * conductivites * CC;
0282 % where:
0283 %   SS  = Unconnected system Matrix
0284 %   CC  = Connectivity Matrix
0285 
0286 p= aa_fwd_parameters( fwd_model );
0287 
0288 d= p.n_dims+1;
0289 e= p.n_elem;
0290 n= p.n_node;
0291 
0292 SSiidx= floor([0:d*e-1]'/d)*d*ones(1,d) + ones(d*e,1)*(1:d) ;
0293 SSjidx= [1:d*e]'*ones(1,d);
0294 SSdata= zeros(d*e,d);
0295 dfact= (d-1)*(d-2); % Valid for d<=3
0296 for j=1:e
0297     a=  inv([ ones(d,1), p.NODE( :, p.ELEM(:,j) )' ]);
0298     idx= d*(j-1)+1 : d*j;
0299     SSdata(idx,1:d)= 2*a(2:d,:)'*a(2:d,:)/dfact/abs(det(a));
0300 end %for j=1:ELEMs
0301 idx= 1:e*d;
0302 SS= sparse(SSiidx,SSjidx,SSdata) * ...
0303     sparse(idx,idx, img.elem_data(ceil(idx/d)) );
0304 
0305 return
0306 
0307 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0308 % TEST CODE FOR MATRIX DERIVATIVES
0309 
0310 % TEST dertiv of det(X + t*a*b')
0311 
0312 d= 1e-6;
0313 X= rand(5);
0314 a=zeros(5,1); a(ceil(5*rand))=1;
0315 b=zeros(5,1); b(ceil(5*rand))=1;
0316 dX_p= (det(X + d*a*b') - det(X) )/d;
0317 dX  = b'*inv(X)*a*det(X);
0318 disp([dX_p dX]);
0319 
0320 % TEST dertiv of inv(X + t*a*b')
0321 dX_p= (inv(X + d*a*b') - inv(X) )/d;
0322 dX  = -inv(X)*a*b'*inv(X);
0323 disp(norm([dX_p-dX],1));
0324 
0325 % TEST dertiv of 1/abs(det(X + t*a*b')) = abs(1/det(X+t*a*b'))
0326 for i=1:20
0327     X= rand(5);
0328     a=zeros(5,1); a(ceil(5*rand))=1;
0329     b=zeros(5,1); b(ceil(5*rand))=1;
0330     dX_p= (1/abs(det(X + d*a*b')) - 1/abs(det(X)) )/d;
0331     dX  = - 1/abs(det(X))*b'*inv(X)*a;
0332     disp(norm([dX_p-dX])/norm(dX));
0333 end