jacobian_adjoint_higher_order

Find the Jacobian associated with an image (and forward model)
Derivative of discretization method
 
 Example (2D):
  imdl = mk_common_model('c2C',16); img=mk_image(imdl.fwd_model,1);
  img.fwd_model.solve = @fwd_solve_higher_order;
  img.fwd_model.system_mat = @system_mat_higher_order;
  img.fwd_model.jacobian = @jacobian_adjoint_higher_order;
  
  vve=[]; JJ4=[];
  for i= 1:3; switch i;
     case 1; img.fwd_model.approx_type = 'tri3'; % linear
     case 2; img.fwd_model.approx_type = 'tri6'; % quadratic
     case 3; img.fwd_model.approx_type = 'tri10'; % cubic;
     end %switch
     vv=fwd_solve(img);      vve(:,i)=vv.meas;
     JJ=calc_jacobian(img);  JJ4(:,i)=JJ(4,:)';
  end

 Example (3D):
  imdl = mk_common_model('b3cr',16);  img=mk_image(imdl.fwd_model,1);
  img.fwd_model.solve = @fwd_solve_higher_order;
  img.fwd_model.system_mat = @system_mat_higher_order;
  img.fwd_model.jacobian = @jacobian_adjoint_higher_order;
  
  vve=[]; JJ4=[];
  for i= 1:2; switch i;
     case 1; img.fwd_model.approx_type = 'tet4'; % linear
     case 2; img.fwd_model.approx_type = 'tet10'; % quadratic
     end %switch
     vv=fwd_solve(img);      vve(:,i)=vv.meas;
     JJ=calc_jacobian(img);  JJ4(:,i)=JJ(4,:)';
  end

0001 function J = jacobian_adjoint_higher_order(fwd_model,img)
0002 %Find the Jacobian associated with an image (and forward model)
0003 %Derivative of discretization method
0004 %
0005 % Example (2D):
0006 %  imdl = mk_common_model('c2C',16); img=mk_image(imdl.fwd_model,1);
0007 %  img.fwd_model.solve = @fwd_solve_higher_order;
0008 %  img.fwd_model.system_mat = @system_mat_higher_order;
0009 %  img.fwd_model.jacobian = @jacobian_adjoint_higher_order;
0010 %
0011 %  vve=[]; JJ4=[];
0012 %  for i= 1:3; switch i;
0013 %     case 1; img.fwd_model.approx_type = 'tri3'; % linear
0014 %     case 2; img.fwd_model.approx_type = 'tri6'; % quadratic
0015 %     case 3; img.fwd_model.approx_type = 'tri10'; % cubic;
0016 %     end %switch
0017 %     vv=fwd_solve(img);      vve(:,i)=vv.meas;
0018 %     JJ=calc_jacobian(img);  JJ4(:,i)=JJ(4,:)';
0019 %  end
0020 %
0021 % Example (3D):
0022 %  imdl = mk_common_model('b3cr',16);  img=mk_image(imdl.fwd_model,1);
0023 %  img.fwd_model.solve = @fwd_solve_higher_order;
0024 %  img.fwd_model.system_mat = @system_mat_higher_order;
0025 %  img.fwd_model.jacobian = @jacobian_adjoint_higher_order;
0026 %
0027 %  vve=[]; JJ4=[];
0028 %  for i= 1:2; switch i;
0029 %     case 1; img.fwd_model.approx_type = 'tet4'; % linear
0030 %     case 2; img.fwd_model.approx_type = 'tet10'; % quadratic
0031 %     end %switch
0032 %     vv=fwd_solve(img);      vve(:,i)=vv.meas;
0033 %     JJ=calc_jacobian(img);  JJ4(:,i)=JJ(4,:)';
0034 %  end
0035 
0036 if ischar(fwd_model) && strcmp(fwd_model,'UNIT_TEST'); do_unit_test; return; end
0037 
0038 if nargin == 1
0039    img= fwd_model;
0040 elseif  strcmp(getfield(warning('query','EIDORS:DeprecatedInterface'),'state'),'on')
0041    warning('EIDORS:DeprecatedInterface', ...
0042       ['Calling JACOBIAN_ADJOINT_HIGHER_ORDER with two arguments is deprecated and will cause' ...
0043        ' an error in a future version. First argument ignored.']);
0044 end
0045 fwd_model= img.fwd_model;
0046 
0047 if mdl_normalize(fwd_model)
0048      fwd_solve_data= fwd_solve( img );   
0049 end
0050 
0051 %Modify the forward model to be of my type
0052 if ~isfield(fwd_model,'approx_type')    || ...
0053    strcmp(fwd_model.approx_type,'tri3') || ...
0054    strcmp(fwd_model.approx_type,'tet4')
0055     %Do nothing
0056 else
0057     [bound,elem,nodes] = fem_1st_to_higher_order(fwd_model);
0058     fwd_model.boundary=bound; fwd_model.elems=elem; fwd_model.nodes=nodes;
0059     img.fwd_model=fwd_model; %CHANGE THIS
0060 end
0061 
0062 %Calculate the total stiffness matrix and elemental stiffness matrices
0063 s_mat = calc_system_mat(img); At=s_mat.E; elemstiff=s_mat.elemstiff;
0064  
0065 %Find electrode stucture and no.of electrodes
0066 %Find stim strucutre and no. stimulations
0067 %Find node structure and find no.nodes
0068 %Find element structure and create vector of length no. elements
0069 elecstruc=fwd_model.electrode; nelecs=size(elecstruc,2);
0070 stimstruc=fwd_model.stimulation; nstims=size(stimstruc,2); 
0071 nodestruc=fwd_model.nodes; nnodes=size(nodestruc,1); 
0072 elemstruc=fwd_model.elems; nelems=size(elemstruc,1); 
0073 
0074 %Find total number of measurements
0075 nmeass=0;
0076 for k=1:nstims
0077     stimkmeasmatrix = stimstruc(k).meas_pattern;
0078     nmeass=nmeass+size(stimkmeasmatrix,1);
0079 end
0080 
0081 %Complete or Point? - Check first electrode and change index vector of
0082 %'node' number corresponding to electrode
0083 elecnode=zeros(1,nelecs);
0084 if(size(elecstruc(1).nodes,2)==1 && size(elecstruc(1).nodes,1)==1) %POINT
0085     %Initialise node to electrode matrix
0086     Node2Elec=sparse(nelecs,nnodes);
0087     for i=1:nelecs
0088         %Assign electrode index at correct node
0089         elecnode(i)=elecstruc(i).nodes;
0090         Node2Elec(i,elecnode(i))=1;
0091     end
0092     %Assign a matrix for derivative of FEM w.r.t conduc
0093     dA_zero=sparse(nnodes,nnodes);
0094     
0095     %Assign correct size unknowns and right hand side matrix (forward)
0096     datafwd=zeros(nnodes,nstims); 
0097     nodeunknownsfwd=zeros(nnodes,nstims); 
0098 else
0099     %Initialise node to electrode matrix
0100     Node2Elec=sparse(nelecs,nnodes+nelecs);
0101     for i=1:nelecs
0102         %Assign electrode at bottom of list
0103         elecnode(i)=nnodes+i;
0104         Node2Elec(i,elecnode(i))=1;
0105     end
0106     
0107     %Assign a matrix for derivative of FEM w.r.t conduc
0108     dA_zero=sparse(nnodes+nelecs,nnodes+nelecs);
0109         
0110     %Assign correct size unknowns and right hand side matrix (forward)
0111     datafwd=zeros(nnodes+nelecs,nstims); 
0112     nodeunknownsfwd=zeros(nnodes+nelecs,nstims); 
0113 end
0114 
0115 %Loop over stimulations and assign current matrix
0116 %CHANGE THIS BY USING NODE2ELEC MATRIX!!!!
0117 for ii=1:nstims
0118     %The vector of current values for stimulation
0119     curnode=stimstruc(ii).stim_pattern;
0120     for i=1:nelecs
0121         datafwd(elecnode(i),ii)=curnode(i);
0122     end
0123 end
0124 
0125 %Create index vector and eliminate ground node equation from index
0126 groundnode=fwd_model.gnd_node; idx=1:size(At,1); idx(groundnode)=[];
0127 
0128 %Solve the simulated linear system with index
0129 nodeunknownsfwd(idx,:)=left_divide(At(idx,idx),datafwd(idx,:));
0130 
0131 %Calculate Jacobian tensor - DE_{i,j,k} == dV_i,j / dS_k
0132 %V_i,j - voltage change on electrode i for stim j
0133 %S_k - conductivity change on element k
0134 DE= zeros(nelecs,nstims,nelems);
0135 
0136 %First step, we only want to pick off the ith electrode
0137 zi2E(:,idx) = Node2Elec(:,idx)/At(idx,idx);
0138 
0139 %SPEED UP HERE
0140 %Factorise A = C'*S*C  - S diagonal conduc (C=system_mat_fields)
0141 %We don't need extra multiplication in loop below
0142 %only for piecewise linear FEM??
0143 %
0144 %zi2E= zeros(nelecs, nnodes+nelecs);
0145 %zi2E(:,idx) = Node2Elec(:,idx)/At(idx,idx);
0146 %zi2E=zi2E*FC'; sv=Fc*sv;
0147 
0148 %Calculate the partial derivative matrix for kth change
0149 for k=1:nelems    
0150     %kth element stiffness matrix, global nodes and index vector
0151     stiffk=elemstiff(k).elemstiff; nodesk=elemstruc(k,:); idx2=1:size(nodesk,2);
0152         
0153     %Create the FEM derivative matrix
0154     dA_dSk=dA_zero; dA_dSk(nodesk(idx2),nodesk(idx2))=stiffk(idx2,idx2);
0155 
0156     %Now form product with solution
0157     DE(:,:,k) = zi2E(:,idx)*dA_dSk(idx,idx)*nodeunknownsfwd(idx,:);
0158 end
0159 
0160 %Calculate Jacobian matrix (measurement patterns specified here)
0161 cntjac=0; J=zeros(nmeass,nelems);
0162 for j=1:nstims   
0163    meas_pat= fwd_model.stimulation(j).meas_pattern;
0164    n_meas  = size(meas_pat,1);
0165    DEj = reshape( DE(:,j,:), nelecs, nelems);
0166    J( cntjac+(1:n_meas),: ) = meas_pat*DEj;
0167    cntjac = cntjac + n_meas;
0168 end; 
0169 
0170 %Get the Jacobian and normalize measurements (if field exists)
0171 if mdl_normalize(fwd_model)
0172      data= fwd_solve_data; % must calculate first, because fwd_model is changed
0173      J= J ./ (data.meas(:)*ones(1,nelems));
0174 end
0175 
0176 %Negative Jacobian for injected currents??
0177 J= -J;  
0178 end
0179 
0180 function do_unit_test
0181    tol = 1e-14;
0182    JJ = do_unit_test_2D(0); % not normalized
0183    JJ_ref = -1e-4*[
0184    2.440415272063380   2.705754096983918   2.721135010947898
0185    2.578623854199123   2.327923064064513   2.342086727271722
0186    1.438743206711758   1.333580385504260   1.337599904647092
0187    1.300534624576032   1.650702059478922   1.659896278693538];
0188 
0189    unit_test_cmp('2D: 1st order',JJ(1:4,1),JJ_ref(1:4,1),tol);
0190    unit_test_cmp('2D: 2nd order',JJ(1:4,2),JJ_ref(1:4,2),tol);
0191    unit_test_cmp('2D: 3rd order',JJ(1:4,3),JJ_ref(1:4,3),tol);
0192 
0193    [JJ1,vve]= do_unit_test_2D(1); for i=1:3
0194    unit_test_cmp('2D: (normalize)',JJ1(:,i),JJ(:,i)/vve(4,i),tol);
0195    end
0196 
0197    JJ = do_unit_test_3D(0);
0198    JJ_ref = -1e-5*[
0199    1.246064580179371   1.585061706092707
0200    1.332632578853691   1.354929239220232
0201    0.712061825721561   0.443297935900921
0202    0.625493827047241   0.604174950085724];
0203    [JJ1,vve]= do_unit_test_3D(1); for i=1:2
0204    unit_test_cmp('3D: (normalize)',JJ1(:,i),JJ(:,i)/vve(4,i),tol);
0205    end
0206 
0207    unit_test_cmp('3D: 1st order',JJ(1:4,1),JJ_ref(1:4,1),tol);
0208    unit_test_cmp('3D: 2nd order',JJ(1:4,2),JJ_ref(1:4,2),tol);
0209 
0210 end
0211 function [JJ4,vve]=do_unit_test_2D(normalize_flag)
0212    imdl = mk_common_model('c2C',16); img = mk_image(imdl.fwd_model,1);
0213    img.fwd_model.normalize_measurements = normalize_flag;
0214    vv=fwd_solve(img);      v0e=vv.meas;
0215    JJ=calc_jacobian(img);  J04=JJ(4,:)';
0216 
0217    %High-order EIDORS solver %Change default eidors solvers
0218    img.fwd_model.solve = @fwd_solve_higher_order;
0219    img.fwd_model.system_mat = @system_mat_higher_order;
0220    img.fwd_model.jacobian = @jacobian_adjoint_higher_order;
0221 
0222    vve=[]; JJ4=[];
0223    for i= 1:3; switch i;
0224       case 1; img.fwd_model.approx_type = 'tri3'; % linear
0225       case 2; img.fwd_model.approx_type = 'tri6'; % quadratic
0226       case 3; img.fwd_model.approx_type = 'tri10'; % cubic;
0227       end %switch
0228       vv=fwd_solve(img);      vve(:,i)=vv.meas;
0229       JJ=calc_jacobian(img);  JJ4(:,i)=JJ(4,:)';
0230    end
0231 
0232    subplot(321);
0233    plot([v0e,vve,(v0e*[1,1,1]-vve)*10]);
0234    legend('Default','linear','quadratic','cubic','(1-0)x10','(2-0)x10','(3-0)x10');
0235    xlim([1,100]);
0236 
0237    imgJJ=img; imgJJ.elem_data = JJ4;
0238    imgJJ.show_slices.img_cols = 3;
0239 
0240    subplot(323); show_slices(imgJJ); eidors_colourbar(imgJJ);
0241 
0242    imgJJ.elem_data = JJ4 - J04*[1,1,1];
0243    subplot(325); show_slices(imgJJ); eidors_colourbar(imgJJ);
0244 end
0245 function [JJ4,vve]=do_unit_test_3D(normalize_flag)
0246    imdl = mk_common_model('b3cr',16); img = mk_image(imdl.fwd_model,1);
0247    img.fwd_model.normalize_measurements = normalize_flag;
0248    vv=fwd_solve(img);      v0e=vv.meas;
0249    JJ=calc_jacobian(img);  J04=JJ(4,:)';
0250 
0251    %High-order EIDORS solver %Change default eidors solvers
0252    img.fwd_model.solve = @fwd_solve_higher_order;
0253    img.fwd_model.system_mat = @system_mat_higher_order;
0254    img.fwd_model.jacobian = @jacobian_adjoint_higher_order;
0255 
0256    vve=[]; JJ4=[];
0257    for i= 1:2; switch i;
0258       case 1; img.fwd_model.approx_type = 'tet4'; % linear
0259       case 2; img.fwd_model.approx_type = 'tet10'; % quadratic
0260       end %switch
0261       vv=fwd_solve(img);      vve(:,i)=vv.meas;
0262       JJ=calc_jacobian(img);  JJ4(:,i)=JJ(4,:)';
0263    end
0264 
0265    subplot(322);
0266    plot([v0e,vve,(v0e*[1,1]-vve)*10]);
0267    legend('Default','linear','quadratic','(1-0)x10','(2-0)x10');
0268    xlim([1,100]);
0269 
0270    imgJJ=img; imgJJ.elem_data = JJ4;
0271    imgJJ.show_slices.img_cols = 2;
0272 
0273    level = [inf,inf,0.3];
0274    subplot(324); show_slices(imgJJ,level); eidors_colourbar(imgJJ);
0275 
0276    imgJJ.elem_data = JJ4 - J04*[1,1];
0277    subplot(326); show_slices(imgJJ,level); eidors_colourbar(imgJJ);
0278 end