jacobian_movement_2p5d_1st_order

PURPOSE ^

JACOBIAN_MOVEMENT_2P5D: J = jacobian_movement_2p5d_1st_order( img )

SYNOPSIS ^

function J = jacobian_movement_2p5d_1st_order( fwd_model, img)

DESCRIPTION ^

 JACOBIAN_MOVEMENT_2P5D: J = jacobian_movement_2p5d_1st_order( img )
 Calculate Jacobian Matrix for current stimulation EIT
 J = Jacobian matrix
 img.fwd_model = forward model
 img.elem_data = background for jacobian calculations

 fwd_model.normalize_measurements if param exists, calculate
                                  a Jacobian for normalized
                                  difference measurements

    img.fwd_solve_2p5d_1st_order.k = [ a .. b ]
        solve, integrating over the range k = a .. b      (default: [0 Inf])
        - provide a single k to get a point solution
        - solve over a reduced range a = 0, b = 3 for a faster solution
    img.fwd_solve_2p5d_1st_order.method = 'name'
        perform numerical integration using the selected method (default: 'quadv')
        'trapz' - trapezoidal integration across the listed points k
        'quadv' - adaptive quadrature (vectorized), from a to b
        'integral' - adaptive quadrature (matlab2012+), from a to b

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function J = jacobian_movement_2p5d_1st_order( fwd_model, img)
0002 % JACOBIAN_MOVEMENT_2P5D: J = jacobian_movement_2p5d_1st_order( img )
0003 % Calculate Jacobian Matrix for current stimulation EIT
0004 % J = Jacobian matrix
0005 % img.fwd_model = forward model
0006 % img.elem_data = background for jacobian calculations
0007 %
0008 % fwd_model.normalize_measurements if param exists, calculate
0009 %                                  a Jacobian for normalized
0010 %                                  difference measurements
0011 %
0012 %    img.fwd_solve_2p5d_1st_order.k = [ a .. b ]
0013 %        solve, integrating over the range k = a .. b      (default: [0 Inf])
0014 %        - provide a single k to get a point solution
0015 %        - solve over a reduced range a = 0, b = 3 for a faster solution
0016 %    img.fwd_solve_2p5d_1st_order.method = 'name'
0017 %        perform numerical integration using the selected method (default: 'quadv')
0018 %        'trapz' - trapezoidal integration across the listed points k
0019 %        'quadv' - adaptive quadrature (vectorized), from a to b
0020 %        'integral' - adaptive quadrature (matlab2012+), from a to b
0021 
0022 
0023 % (C) 2016 A Boyle
0024 % License: GPL version 2 or version 3
0025 
0026 % correct input paralemeters if function was called with only img
0027 if ischar(fwd_model) && strcmp(fwd_model,'UNIT_TEST'); do_unit_test; return; end
0028 
0029 if nargin == 1
0030    img = fwd_model;
0031 elseif  strcmp(getfield(warning('query','EIDORS:DeprecatedInterface'),'state'),'on')
0032    warning('EIDORS:DeprecatedInterface', ...
0033       ['Calling JACOBIAN_ADJOINT with two arguments is deprecated and will cause' ...
0034        ' an error in a future version. First argument ignored.']);
0035 end
0036 
0037 img = data_mapper(img);
0038 if ~ismember(img.current_params, supported_params)
0039     error('EIDORS:PhysicsNotSupported', '%s does not support %s', ...
0040     'JACOBIAN_ADJOINT',img.current_params);
0041 end
0042 
0043 org_params = img.current_params;
0044 % all calcs use conductivity
0045 img = convert_img_units(img, 'conductivity');
0046 
0047 img.elem_data = check_elem_data_and_apply_c2f(img);
0048 % TODO we don't really want to expand out the c2f
0049 
0050 fwd_model = img.fwd_model;
0051 
0052 pp = fwd_model_parameters( fwd_model );
0053 gnd = img.fwd_model.gnd_node;
0054 
0055 img.fwd_model.system_mat = @system_mat_2p5d_1st_order;
0056 gnd = fwd_model.gnd_node;
0057 img.fwd_model.system_mat_2p5d_1st_order.k = 0;
0058 img.fwd_model.system_mat_2p5d_1st_order.factory = 1;
0059 
0060 k = [0 Inf]; try k = img.fwd_model.jacobian_movement_2p5d_1st_order.k; end
0061 method = 'quadv'; try method = img.fwd_model.jacobian_movement_2p5d_1st_order.method; end
0062 
0063 pp.Sf = system_mat_2p5d_1st_order( img ); % returns a function Sf(k) that builds a system matrix for 'k'
0064 pp.CC = connectivity_matrix( pp );
0065 % build sparse DD * CC: conductivity * connectivity
0066 lCC = size(pp.CC,1);
0067 DD = kron(img.elem_data,ones(pp.n_dims,1));
0068 DD(end+1:lCC) = 1; % CEM
0069 pp.DD = spdiags(DD, 0, lCC, lCC); % sparse diagonal matrix
0070 pp.DC = pp.DD * pp.CC;
0071 clear DD;
0072 % shape derivatives
0073 [pp.dSx, pp.dTx] = shape_derivatives(pp, img);
0074 
0075 % numerical integration
0076 if length(k) == 1 % singleton k
0077    J = 2*jacobian_k(k, pp, gnd, img.fwd_model.stimulation);
0078 else
0079    switch method
0080       case 'trapz'
0081          % less accurate: trapz
0082          trace = 0;
0083          if trace; fprintf('%8s %12s %16s %16s %16s\n','fcnt','a','b-a','||Q||','d||Q||'); end
0084          n = 0; kil = k(1); nJl=0;
0085          tol = 1e-8;
0086          k(isinf(k)) = tol^(-1/6); % 1/k^2 ^3
0087          Jf = zeros(pp.n_meas, pp.n_elec * pp.n_dims, length(k)); % voltages under electrodes elec x stim, frequency domain
0088          for ki = k
0089             n = n + 1;
0090             Jf(:,:,n) = jacobian_k(ki, pp, gnd, img.fwd_model.stimulation);
0091             if trace
0092                nJ = norm(Jf(:,:,n));
0093                fprintf('%8d     %12e %16e %16e %16e\n',n,ki,ki-kil,nJ,nJ-nJl);
0094                kil = ki; nJl = nJ;
0095             end
0096          end
0097          J = 2/pi*trapz(k,Jf,3);
0098          if 0 % draw J(k) to check we are integrating over a large enough range
0099             Jff = squeeze(reshape(Jf,pp.n_meas*pp.n_elec*pp.n_dims,1,length(k)));
0100             slope = (Jff(:,3)-Jff(:,2))./max(abs(Jff),[],2);
0101             clf; plot(slope); drawnow; pause(1);
0102             Jff_normalized = bsxfun(@rdivide, Jff, Jff(:,1)); % Jff ./ J_0
0103             clf;semilogx(repmat(k,pp.n_meas*pp.n_elec*pp.n_dims,1)',Jff_normalized');
0104             xlabel('k'); ylabel('J_k'); title('J_k/J_0')
0105             drawnow; pause(1);
0106          end
0107       case 'quadv'
0108          % more accurate: adaptive gaussian quadrature
0109          trace = 0;
0110          reltol = 1e-4;
0111          tol = norm(jacobian_k(0, pp, gnd, img.fwd_model.stimulation))*reltol;
0112          kend = min(tol^(-1/6), max(k)); % don't go too far... k = Inf is a singular matrix, stop adjacent to numeric singularity
0113          % quadv is scheduled to be removed from matlab eventually... but it is
0114          % WAY faster than integral with any tolerance configuration I could identify
0115          % disp('       9     0.0000000000     2.71580000e+15 3286857082845.8784179688');
0116          if trace; fprintf('%8s %12s %16s %16s\n','fcnt','a','b-a','Q(1)'); end
0117          J = 2/pi*quadv(@(kk) jacobian_k(kk, pp, gnd, img.fwd_model.stimulation), k(1), kend, tol, trace);
0118       case 'integral'
0119          reltol = 1e-4;
0120          tol = norm(jacobian_k(0, pp, gnd, img.fwd_model.stimulation))*reltol;
0121          kend = min(tol^(-1/6), max(k)); % don't go too far... k = Inf is a singular matrix, stop adjacent to numeric singularity
0122          opts = {'ArrayValued', true,
0123                  'AbsTol', tol, % default: 1e-10
0124                  'RelTol', reltol}; % default:  1e-6
0125          opts = opts';
0126          % the integral solution is about 10x slower (5.47 seconds vs. 0.60 seconds for UNIT_TEST)
0127          % ... I played with AbsTol and RelTol but wasn't able to affect the outcome
0128          J = 2/pi*integral(@(kk) jacobian_k(kk, pp, gnd, img.fwd_model.stimulation), k(1), kend, opts{:});
0129       otherwise
0130          error(['unrecognized method: ' method]);
0131    end
0132 end
0133 
0134 if ~strcmp(org_params,'conductivity')
0135     J = apply_chain_rule(J, img, org_params);
0136     if isfield(fwd_model, 'coarse2fine') && ...
0137           size(img.elem_data,1) == size(fwd_model.coarse2fine,1)
0138             J = J*fwd_model.coarse2fine;
0139             nparam = size(fwd_model.coarse2fine,2);
0140     end
0141 end
0142 
0143 %restore img to original condition
0144 if ~strcmp(org_params,'conductivity') || isfield(img, org_params)
0145     img = rmfield(img,'elem_data');
0146     img.current_params = [];
0147 end
0148 
0149 % calculate normalized Jacobian
0150 if pp.normalize
0151    data = fwd_solve( img );
0152    J = J ./ (data.meas(:)*ones(1,nparam));
0153 end
0154 
0155 function J = jacobian_k(k, pp, gnd, stim)
0156    SS = pp.Sf(k);
0157    idx = 1:size(SS,1);
0158    idx( gnd ) = [];
0159    
0160    sv = potential_k(SS, pp, gnd);
0161    zi2E = -(pp.N2E(:,idx)/ SS(idx,idx)) * pp.CC(:,idx)';
0162    SE = jacobian_calc(pp, zi2E, pp.dSx, k, pp.dTx, sv);
0163    
0164    idx = 0;
0165    J = zeros( pp.n_meas, pp.n_elec * pp.n_dims );
0166    for j = 1:pp.n_stim
0167       meas_pat = stim(j).meas_pattern;
0168       m = size(meas_pat,1); % new measurements added
0169       SEj = reshape( SE(:,j,:), pp.n_elec, pp.n_elec * pp.n_dims);
0170       J( idx+(1:m),: ) = meas_pat * SEj;
0171       idx = idx + m;
0172    end
0173 
0174 function v = potential_k(S, pp, gnd)
0175    idx = 1:size(S,1);
0176    idx( gnd ) = [];
0177    idx2 = find(any(pp.N2E));
0178    v = zeros(pp.n_node, pp.n_stim); % voltages at all nodes x number of stim, frequency domain
0179    v(idx,:) = left_divide( S(idx,idx), 1/2*pp.QQ(idx,:) );
0180 
0181 % SE_{i,j,k} is dV_i,j / dS_k
0182 %  where V_i is change in voltage on electrode i for
0183 %        stimulation pattern j
0184 %        S_k is k=n*ne+e for change in position on electrode e of ne, for movement axis n=(0,1)
0185 function SE = jacobian_calc(pp, zi2E, dSx, k, dTx, sv)
0186    N = pp.n_elec*pp.n_dims;
0187    SE = zeros(pp.n_elec, pp.n_stim, N);
0188    for n = 1:N;
0189       dSxk2dTx = dSx{n} + k^2 * dTx{n};
0190       SE(:,:,n) = (zi2E * (dSxk2dTx * pp.DC)) * sv;
0191    end
0192 
0193 function J = apply_chain_rule(J, img, org_params)
0194    switch(org_params)
0195       case 'resistivity'
0196          dCond_dPhys = -img.elem_data.^2;
0197       case 'log_resistivity'
0198          dCond_dPhys = -img.elem_data;
0199       case 'log_conductivity'
0200          dCond_dPhys = img.elem_data;
0201       otherwise
0202          error('not implemented yet')
0203    end
0204    J = J.*repmat(dCond_dPhys ,1,size(J,1))';
0205 
0206 function elem_data = check_elem_data_and_apply_c2f(img)
0207    elem_data = img.elem_data;
0208    sz_elem_data = size(elem_data);
0209    if sz_elem_data(2) ~= 1;
0210       error('jacobian_adjoin: can only solve one image (sz_elem_data = %)', ...
0211             sz_elem_data);
0212    end
0213 
0214    if isfield(img.fwd_model, 'coarse2fine');
0215       c2f = img.fwd_model.coarse2fine;
0216       sz_c2f = size(c2f);
0217       switch sz_elem_data(1)
0218          case sz_c2f(1); % Ok
0219          case sz_c2f(2); elem_data = c2f * elem_data;
0220          otherwise; error(['jacobian_adjoint: provided elem_data ' ...
0221               ' (sz = %d) does not match c2f (sz = %d %d)'], sz_elem_data(1), sz_c2f);
0222       end
0223    else
0224       if sz_elem_data(1) ~= num_elems(img.fwd_model)
0225          error(['jacobian_adjoint: provided elem_data (sz = %d) does ' ...
0226             ' not match fwd_model (sz = %d)'], sz_elem_data(1), num_elems(sz_c2f));
0227       end
0228    end
0229 
0230 function str = supported_params
0231     str = {'conductivity'
0232            'resistivity'
0233            'log_conductivity'
0234            'log_resistivity'};
0235 
0236 % Define the element connectivity matrix Ce
0237 function Ce = connectivity_matrix( pp, img );
0238    m = (pp.n_dims+1)*pp.n_elem;
0239    n = size(pp.N2E,2);
0240    ii = 1:m; % row indices
0241    jj = pp.ELEM(:); % col indices
0242    ss = ones(m,1); % values = 1
0243    Ce = sparse(ii, jj, ss, m, n); % m x n sparse matrix
0244 
0245 function [dS, dT] = shape_derivatives( pp, img )
0246    d1 = pp.n_dims+1; % 2d --> always 3
0247    sz = d1 * pp.n_elem; % dS and dT matrices should be 'sz' rows x cols
0248    Ec = cell(pp.n_elem,1);
0249    for elec = 1:length(img.fwd_model.electrode)
0250       ii = []; jj = []; ss = []; tt = [];
0251       for node = img.fwd_model.electrode(elec).nodes(:)';
0252          % find elements touching that node
0253          [ idx, elem ] = find(pp.ELEM == node);
0254          for j = 1:length(elem(:))
0255             e = elem(j);
0256             % build shape matrices
0257             if ~any(size(Ec{e})) % already calculated?
0258                Ec{e} = inv([ ones(d1,1), pp.NODE( :, pp.ELEM(:,e) )' ]); %%% (12) [Boyle2016]
0259             end
0260             E = Ec{e};
0261             E1 = E(2:end,:); % E_/1
0262             absdetE = 0.5/pp.VOLUME(e); % | det E | = 1 / (n_dim! * element volume) for n_dim=2
0263             for d = 1:pp.n_dims
0264                u = zeros(d1,1); v = u; % init
0265                u(idx(j)) = 1; % element node# to perturb
0266                v(d+1) = 1; % direction for perturbation: x,z?
0267                B = (v'*E*u) / absdetE; %%% (13) [Boyle2016]
0268                dE1 = -E*u*v'*E; dE1 = dE1(2:end,:); %%% (14) [Boyle2016]
0269                dSe = (v'*E*u*E1'*E1 + dE1'*E1 + E1'*dE1)*pp.VOLUME(e); %%% eqn (11) [Boyle2016]
0270 %               dTe = -absdetE * (v'*E*u) * (ones(d1)+eye(d1))/12; %%% eqn (20) [Boyle2016]
0271                dTe = (v'*E*u)/(2*absdetE) * (ones(d1)+eye(d1))/12; %%% eqn (20) [Boyle2016]
0272                dSs(:,d) = dSe(:); % square matrix to column
0273                dTs(:,d) = dTe(:);
0274                if 0 % test
0275                   dfact = factorial(pp.n_dims);
0276                   Ae = pp.NODE(:,pp.ELEM(:, e))'; Ae = [ones(pp.n_dims+1,1), Ae]; Ae = inv(Ae);
0277                   Be = Ae(2:pp.n_dims+1,:);
0278                   absdetAe = abs(det(Ae));
0279                   u = zeros(pp.n_dims+1,1); v = u; u(idx(j)) = 1; v(d+1) = 1;
0280                   delBe = -Ae*u*v'*Ae; delBe = delBe(2:pp.n_dims+1,:);
0281                   dSe_old = 1/dfact*(1/absdetAe*v'*Ae*u*Be'*Be + delBe'*Be/absdetAe + Be'*delBe/absdetAe);
0282                   %%%
0283                   Ae = pp.NODE(:,pp.ELEM(:, e))'; Ae = [ones(pp.n_dims+1,1), Ae]; Ae = inv(Ae);
0284                   Be = Ae(2:pp.n_dims+1,:);
0285                   Se = Be'*Be/dfact/abs(det(Ae));
0286                   Te = 1/(2*abs(det(Ae)))*(ones(d1)+eye(d1))/12;
0287                   delta = 1e-8; d_NODE = pp.NODE; d_NODE(d,node) = d_NODE(d,node) + delta;
0288                   Ae = d_NODE(:,pp.ELEM(:, e))'; Ae = [ones(pp.n_dims+1,1), Ae]; Ae = inv(Ae);
0289                   Be = Ae(2:pp.n_dims+1,:);
0290                   Se_delta = Be'*Be/dfact/abs(det(Ae));
0291                   Te_delta = 1/(2*abs(det(Ae)))*(ones(d1)+eye(d1))/12;
0292                   dSe_pert = (Se_delta - Se) / delta;
0293                   dTe_pert = (Te_delta - Te) / delta;
0294 
0295                   if norm(dSe_pert - dSe,1)/norm(dSe_pert) > 5e-6
0296                      eidors_msg(sprintf('@@@ dSe: elec#%d, elem#%d (%d of %d), %s-axis: calc wrong', ...
0297                                         elec, e, j, length(elem(:)), 'w'+d), ...
0298                                 1);
0299                      error('stop');
0300                   end
0301                   if norm(dTe_pert - dTe)/norm(dTe_pert) > 5e-6
0302                      eidors_msg(sprintf('@@@ dTe: elec#%d, elem#%d (%d of %d), %s-axis: calc wrong', ...
0303                                         elec, e, j, length(elem(:)), 'w'+d), ...
0304                                 1);
0305                      error('stop');
0306                   end
0307                end
0308             end
0309             se_idx = (1:d1)+(e-1)*d1;
0310             [xx,yy] = meshgrid(se_idx, se_idx);
0311             ii(end+1:end+d1^2) = xx(:);
0312             jj(end+1:end+d1^2) = yy(:);
0313             ss(end+1:end+d1^2,:) = dSs;
0314             tt(end+1:end+d1^2,:) = dTs;
0315             % note that for any element, we usually have more than one node
0316             % that is perturbed; all ss (and, separately, tt) for that element
0317             % are summed by the sparse() operation so that we have a sparse
0318             % block-wise matrix with 3x3 blocks on the diagonal in 2D
0319          end
0320       end
0321       % return a dS and dT cell array with a cell per dimension
0322       for d = 1:pp.n_dims
0323          idx = (d-1)*pp.n_elec + elec;
0324          dS{idx} = sparse(ii,jj,ss(:,d),sz,sz);
0325          dT{idx} = sparse(ii,jj,tt(:,d),sz,sz);
0326          if 0 % test dS against perturbation
0327             elec_nodes = img.fwd_model.electrode(elec).nodes(:);
0328             dS_old = calc_delVm( elec_nodes, pp, img.fwd_model, img, d, 1, 1);
0329             clf; spy(pp.CC'*dS_old*pp.CC,'go'); hold on; spy(pp.CC'*dS{idx}*pp.CC,'rx'); hold off;
0330             if norm(dS{idx} - dS_old, 1) > 2e-5
0331                eidors_msg(sprintf('@@@ dS{%d}: elec#%d, %s-axis movement: calc wrong', ...
0332                                   idx, ...
0333                                   mod(idx-1,pp.n_elec)+1, ...
0334                                   'x'+floor((idx-1)/pp.n_elec)), ...
0335                           1);
0336                error('stop');
0337             end
0338          end
0339       end
0340    end
0341 
0342 %%% TODO rm start (2d movement Jacobian)
0343 function delVm=  calc_delVm( elec_nodes_array, pp, fwd_model, img, colidx, Re_Ce, cond_Ce_Vc)
0344    I = []; J=[]; S= [];
0345    pp.dfact = factorial(pp.n_dims);
0346    for elec_nodes= elec_nodes_array(:)';
0347       [rowidx, elemidx] = find(pp.ELEM == elec_nodes);
0348       % Define the system sensitivity matrix to movement delSm
0349       sz= (pp.n_dims+1)*pp.n_elem;
0350       %delSm = sparse(sz,sz);
0351       % For each touching element, calculate the perturbation
0352       jcount = 1;
0353       for j = elemidx'
0354           % Extract the coordinates of the element's four nodes
0355           Ae = pp.NODE(:,pp.ELEM(:, j))';
0356           % Define the invertible matrix P: augment Ae by adding a
0357           % column of ones to invert
0358           P = [ones(pp.n_dims+1,1), Ae];
0359           Ae = inv(P);
0360           absdetAe = abs(det(Ae));
0361           % Define Be as the matrix Ae with row 1 deleted
0362           Be = Ae(2:pp.n_dims+1,:);
0363           % For this coordinate, perturb P by [rowidx,colidx], which are
0364           % our paper's perturbation vectors [a,b]
0365           a = zeros(pp.n_dims+1,1);
0366           b = a;
0367           a(rowidx(jcount)) = 1;
0368           jcount = jcount + 1;
0369           b(colidx+1) = 1;
0370           % Calculate the system submatrix subSm for the element j by
0371           % asymmetric perturbation of the electrode node k
0372           deldetAe =   1/absdetAe*b'*Ae*a;
0373           delBe = -Ae*a*b'*Ae;
0374           delBe = delBe(2:pp.n_dims+1,:);
0375           subSm = 1/pp.dfact*(...
0376               deldetAe*Be'*Be + ...
0377               delBe'*Be/absdetAe + ...
0378               Be'*delBe/absdetAe);
0379       
0380           % Embed subSm into delSm such that subSm(1,1) is the
0381           % (4j+1,4j+1) element of delSm
0382           se_idx= (pp.n_dims+1)*j+(-pp.n_dims : 0);
0383           switch pp.n_dims
0384              case 2
0385                 Iidx = vertcat(se_idx,se_idx,se_idx);
0386                 I = [I Iidx(:)];
0387                 J = [J, se_idx,se_idx,se_idx];
0388                 
0389              case 3
0390                 Iidx = vertcat(se_idx,se_idx,se_idx,se_idx);
0391                 I = [I Iidx(:)];
0392                 J = [J, se_idx,se_idx,se_idx,se_idx];
0393           end
0394           S = [S subSm(:)];
0395       end
0396    end
0397    delSm = sparse(I,J,S,sz,sz);
0398    delVm = Re_Ce * delSm * cond_Ce_Vc;
0399 
0400 function do_unit_test()
0401    nn = 16; % number of electrodes
0402    imdl2 = mk_geophysics_model('h22c',nn);
0403    img2 = mk_image(imdl2,1);
0404    nm = size(stim_meas_list(img2.fwd_model.stimulation),1); % number of measurements
0405    ne = size(imdl2.rec_model.elems,1); % number of elements in coarse model
0406    img2.fwd_model.conductivity_jacobian = zeros(nm,ne);
0407 
0408    % for the 3d model, we throw out the rec_model and inject the
0409    % imdl2.rec_model, then recalculate the c2f so that we can compare apples to
0410    % apples when we take ||Jxxx - J3||
0411    imdl3 = mk_geophysics_model('h32c',nn);
0412    for s = {'nodes', 'elems', 'boundary', 'name','electrode'}
0413       imdl3.rec_model.(s{:}) = imdl2.rec_model.(s{:});
0414    end
0415    disp('recalculating 3D c2f, so that coarse meshes agree (2D vs 3D)');
0416    [c2f,bkgnd] = mk_approx_c2f(imdl3.fwd_model, imdl3.rec_model);
0417    imdl3.fwd_model.coarse2fine = c2f;
0418    imdl3.fwd_model.background = bkgnd;
0419    img3 = mk_image(imdl3,1);
0420    img3.fwd_model.conductivity_jacobian = zeros(nm,ne);
0421 
0422    % confirm c2f is sane
0423    ctr = interp_mesh(imdl2.rec_model);
0424    sel = ((ctr(:,1)-100).^2 < 80) & ((ctr(:,2)+50).^2<80);
0425    img2.elem_data = imdl2.fwd_model.coarse2fine*(1 + sel*9) + -10*imdl2.fwd_model.background;
0426    img3.elem_data = imdl3.fwd_model.coarse2fine*(1 + sel*9) + -10*imdl3.fwd_model.background;
0427    clf;
0428    subplot(121);show_fem(img2,1);
0429    subplot(122);show_fem(img3,1);
0430    drawnow;
0431    % set data to homogeneous
0432    img2.elem_data = imdl2.fwd_model.background*0+1;
0433    img3.elem_data = imdl3.fwd_model.background*0+1;
0434 
0435    t = tic;
0436    img2.fwd_model.nodes(1,:) = img2.fwd_model.nodes(1,:) + rand(1,2)*1e-8; % defeat cache
0437    J2p = jacobian_movement_perturb(img2);
0438    J2p(:,1:(end-2*nn)) = [];
0439    fprintf(' 2D perturb = %.2f sec\n', toc(t));
0440    t = tic;
0441    img2.fwd_model.nodes(1,:) = img2.fwd_model.nodes(1,:) + rand(1,2)*1e-8; % defeat cache
0442    J2a = jacobian_movement(img2);
0443    J2a(:,1:(end-2*nn)) = [];
0444    fprintf(' 2D adjoint = %.2f sec\n', toc(t));
0445 
0446    t = tic;
0447    img2.fwd_model.nodes(1,:) = img2.fwd_model.nodes(1,:) + rand(1,2)*1e-8; % defeat cache
0448    img2.fwd_model.jacobian = @jacobian_movement_2p5d_1st_order;
0449    img2.fwd_model.jacobian_movement_2p5d_1st_order.k = 0;
0450    J2p50 = jacobian_movement_2p5d_1st_order(img2);
0451    J2p50(:,1:(end-2*nn)) = [];
0452    fprintf(' 2.5D (k = 0) = %.2f sec\n', toc(t));
0453 
0454    ke = [ 0 0 0 ];
0455    t = tic;
0456    img2.fwd_model.nodes(1,:) = img2.fwd_model.nodes(1,:) + rand(1,2)*1e-8; % defeat cache
0457    img2.fwd_model.jacobian = @jacobian_movement_2p5d_1st_order;
0458    img2.fwd_model.jacobian_movement_2p5d_1st_order.k = [0 logspace(-3,+1,20)];
0459    img2.fwd_model.jacobian_movement_2p5d_1st_order.method = 'trapz';
0460    J2p5kt = jacobian_movement_2p5d_1st_order(img2);
0461    J2p5kt(:,1:(end-2*nn)) = [];
0462    ke(1) = img2.fwd_model.jacobian_movement_2p5d_1st_order.k(end);
0463    fprintf(' 2.5D (k = 0..%.1f, trapz) = %.2f sec\n', ke(1), toc(t));
0464    t = tic;
0465    img2.fwd_model.nodes(1,:) = img2.fwd_model.nodes(1,:) + rand(1,2)*1e-8; % defeat cache
0466    img2.fwd_model.jacobian = @jacobian_movement_2p5d_1st_order;
0467    img2.fwd_model.jacobian_movement_2p5d_1st_order.k = [0 Inf];
0468    img2.fwd_model.jacobian_movement_2p5d_1st_order.method = 'quadv';
0469    J2p5kq = jacobian_movement_2p5d_1st_order(img2);
0470    J2p5kq(:,1:(end-2*nn)) = [];
0471    ke(2) = img2.fwd_model.jacobian_movement_2p5d_1st_order.k(end);
0472    fprintf(' 2.5D (k = 0..%.1f, quadv) = %.2f sec\n', ke(2), toc(t));
0473    t = tic;
0474    img2.fwd_model.nodes(1,:) = img2.fwd_model.nodes(1,:) + rand(1,2)*1e-8; % defeat cache
0475    img2.fwd_model.jacobian = @jacobian_movement_2p5d_1st_order;
0476    img2.fwd_model.jacobian_movement_2p5d_1st_order.k = [0 Inf];
0477    img2.fwd_model.jacobian_movement_2p5d_1st_order.method = 'integral';
0478    J2p5ki = jacobian_movement_2p5d_1st_order(img2);
0479    J2p5ki(:,1:(end-2*nn)) = [];
0480    ke(3) = img2.fwd_model.jacobian_movement_2p5d_1st_order.k(end);
0481    fprintf(' 2.5D (k = 0..%.1f, integral) = %.2f sec\n', ke(3), toc(t));
0482 
0483    t = tic;
0484    img3.fwd_model.nodes(1,:) = img3.fwd_model.nodes(1,:) + rand(1,3)*1e-8; % defeat cache
0485    J3a = jacobian_movement(img3);
0486    J3a(:,1:(end-3*nn)) = []; J3a(:,(1*nn+1):(2*nn)) = []; % delete y-axis
0487    fprintf(' 3D adjoint = %.2f sec\n', toc(t));
0488    if 0 % skip really slow calculations
0489       t = tic;
0490       img3.fwd_model.nodes(1,:) = img3.fwd_model.nodes(1,:) + rand(1,3)*1e-8; % defeat cache
0491       J3p = jacobian_movement_perturb(img3);
0492       J3p(:,1:(end-3*nn)) = []; J3p(:,(1*nn+1):(2*nn)) = []; % delete y-axis
0493       fprintf(' 3D perturb = %.2f sec\n', toc(t));
0494    else
0495       J3p = J3a;
0496       fprintf(' 3D perturb = <SKIP>\n');
0497    end
0498 
0499 
0500    tol = 1e-8;
0501    reltol = 1e2;
0502    fprintf('tol = %0.1e, reltol = %0.1e, ||J_0|| = %0.1e, ||J|| = %0.1e\n', tol, reltol,norm(J2a),norm(J3a));
0503    unit_test_cmp('2D perturb                 vs 3D', max(max(abs(J2p-J3p)))>tol,1);
0504    unit_test_cmp('2D adjoint                 vs 3D', max(max(abs(J2a-J3a)))>tol,1);
0505    unit_test_cmp('2D adjoint vs 2D perturb        ', J2a, J2p,  tol*10);
0506    unit_test_cmp('3D adjoint vs 3D perturb        ', J3a, J3p,  tol);
0507    unit_test_cmp('2.5D (k = 0)                  vs 2D',  J2p50, J2a, tol);
0508    unit_test_cmp('2.5D (k = 0) matches          == 2D',  max(max(abs(J2p50-J2a)))<tol,1);
0509    unit_test_cmp('2.5D (k = 0) does not match   != 3D',  max(max(abs(J2p50-J3a)))>tol,1);
0510    unit_test_cmp(sprintf('2.5D (k = 0..%-4.1f) (trapz)    != 2D',ke(1)), max(max(abs(J2a./J2p5kt-1)))>reltol, 1);
0511    unit_test_cmp(sprintf('2.5D (k = 0..%-4.1f) (trapz)    == 3D',ke(1)), max(max(abs(J3a./J2p5kt-1)))<reltol, 1);
0512    unit_test_cmp(sprintf('2.5D (k = 0..%-4.1f) (trapz)    vs 3D',ke(1)), J3a./J2p5kt, 1, reltol);
0513    unit_test_cmp(sprintf('2.5D (k = 0..%-4.1f) (quadv)    vs 3D',ke(2)), J3a./J2p5kq, 1, reltol);
0514    unit_test_cmp(sprintf('2.5D (k = 0..%-4.1f) (integral) vs 3D',ke(3)), J3a./J2p5ki, 1, reltol);
Generated on Tue 31-Dec-2019 17:03:26 by m2html © 2005