adjoint_spin

function [JTb] = adjoint_spin(vtx,simp,elec,x,gnd_ind,zc,I,no_pl,Vmes);

The function calculates the product J'*b, i.e. Jacobian transpose times 
a (measurements Vmes) vector b, using the adjoint sources formulation.
 


JTb     = J'*Vmes
vtx     = The vertices matrix
simp    = The simplices matrix
elec    = The electrodes matrix
x       = The solution estimate based on which J is supposed to be caclulated
gnd_ind = The grounf index
zc      = The electrode contact impedance
I       = The current patterns
no_pl   = Number of electrode planes
Vmes    = The measurements vector

0001 function [JTb] = adjoint_spin(vtx,simp,elec,x,gnd_ind,zc,I,no_pl,Vmes);
0002 %function [JTb] = adjoint_spin(vtx,simp,elec,x,gnd_ind,zc,I,no_pl,Vmes);
0003 %
0004 %The function calculates the product J'*b, i.e. Jacobian transpose times
0005 %a (measurements Vmes) vector b, using the adjoint sources formulation.
0006 %
0007 %
0008 %
0009 %JTb     = J'*Vmes
0010 %vtx     = The vertices matrix
0011 %simp    = The simplices matrix
0012 %elec    = The electrodes matrix
0013 %x       = The solution estimate based on which J is supposed to be caclulated
0014 %gnd_ind = The grounf index
0015 %zc      = The electrode contact impedance
0016 %I       = The current patterns
0017 %no_pl   = Number of electrode planes
0018 %Vmes    = The measurements vector
0019 
0020 
0021 spin = 1;%size(I,2);
0022 [vr] = size(vtx,1);
0023 ptr = 0;
0024     
0025      %Update model based on last x
0026      [E,D,Ela,pp] = fem_master_full(vtx,simp,x,gnd_ind,elec,zc,'{n}');
0027      EB = kron(eye(spin),inv(E));    
0028      
0029           
0030      for j=1:spin:size(I,2); %For each group of current patterns
0031      
0032      Ib = I(:,j:j+spin-1);
0033      IB = reshape(Ib,size(Ib,1)*size(Ib,2),1);
0034      
0035      VB = EB*IB;
0036          
0037      
0038      %The current fields for I(: ,spin) only
0039      VB1 = reshape(VB,size(vtx,1)+size(elec,1),spin);
0040      VB2 = VB1(1:size(vtx,1),:); 
0041      %These are the standard current fields used in the calculation of the Jacobian
0042      Vjcf = reshape(VB2,size(vtx,1)*spin,1);
0043      %Electrode potentials removed
0044      
0045      volts = []; %Some simulated data based on x
0046      ind = [];
0047      dfh = [];
0048      
0049      [voltageH,voltageV,indH,indV,df] = get_3d_meas(elec,vtx,VB1,... 
0050                                                        I(size(vtx,1)+1:end,j:j+spin-1),no_pl);
0051      volts = [volts;voltageH];
0052      ind = [ind;indH];
0053      dfh = [dfh;df(1:2:end)];
0054        
0055   
0056           
0057          %The measurement fields for the measurements during I(:,spin) are active
0058          VmI = Vmes(ptr+1:ptr+sum(dfh)); %Part of the measurements
0059          ptr  = ptr+sum(dfh);
0060      
0061                 
0062          %Set up the measurement patterns for all I(:,spin)
0063          Imb = [];
0064          MC = [];
0065          kap = 1;
0066 
0067              for ij=1:size(ind,1)
0068                  m_n = zeros(size(elec,1),1);
0069                  m_n(ind(ij,1)) = 1;
0070                  m_n(ind(ij,2)) = -1;
0071                  MC = [MC,m_n];
0072                  if ij == sum(dfh(1:kap))
0073                      if ij ~= size(ind,1)
0074                      kap = kap+1;
0075                      end
0076                      Imb = sparse(blkdiag(Imb,MC(:,1+ij-dfh(kap):end)));
0077                  end
0078             end
0079              
0080                  
0081          %These are measurement residuals of the j'th current pattern
0082          b = VmI;
0083               
0084          Im = Imb*b; 
0085          
0086          %Reshape it in columns.
0087          Im_col = reshape(Im,size(elec,1),spin);
0088          mul = 1;
0089          
0090                     I_s = [];
0091                     for t=1:spin
0092                         I_s = [I_s; zeros(size(vtx,1),1); Im_col(:,t)*mul];
0093                     end
0094      
0095                     Vmf = EB*I_s;
0096                     Vmf1 = reshape(Vmf,size(vtx,1)+size(elec,1),spin);
0097                     Vmf2 = Vmf1(1:size(vtx,1),:);
0098                     Vjmf = reshape(Vmf2,size(vtx,1)*spin,1); %
0099                     
0100                     DM = kron(eye(spin),D); % Gradient operator
0101                     Vjmf_n = -DM * Vjmf;
0102                     Vjcf_n = -DM * Vjcf;
0103                     JTb_p = Vjmf_n .* Vjcf_n ;  
0104                     J_m = reshape(JTb_p,size(simp,1)*3,spin);
0105                     if spin > 1
0106                     JTb_x3 = sum(J_m.');
0107                 else
0108                     JTb_x3 = J_m.';
0109                 end
0110                     JTb_u = JTb_x3(1:3:end) + JTb_x3(2:3:end) + JTb_x3(3:3:end);
0111         
0112                     JTb = -JTb_u.' .* diag(Ela(1:3:end,1:3:end));
0113         
0114    end
0115    
0116    
0117 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0118 % This is part of the EIDORS suite.
0119 % Copyright (c) N. Polydorides 2003
0120 % Copying permitted under terms of GNU GPL
0121 % See enclosed file gpl.html for details.
0122 % EIDORS 3D version 2.0
0123 % MATLAB version 5.3 R11
0124 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0125