jacobian_adjoint_2p5d_1st_order

 JACOBIAN_ADJOINT_2P5D: J= jacobian_adjoint_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

0001 function J= jacobian_adjoint_2p5d_1st_order( fwd_model, img)
0002 % JACOBIAN_ADJOINT_2P5D: J= jacobian_adjoint_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(img);
0048 
0049 fwd_model= img.fwd_model;
0050 
0051 pp= fwd_model_parameters( fwd_model, 'skip_VOLUME' );
0052 gnd = img.fwd_model.gnd_node;
0053 
0054 img.fwd_model.system_mat = @system_mat_2p5d_1st_order;
0055 gnd = fwd_model.gnd_node;
0056 img.fwd_model.system_mat_2p5d_1st_order.k = 0;
0057 img.fwd_model.system_mat_2p5d_1st_order.factory = 1;
0058 
0059 pp.Sf = system_mat_2p5d_1st_order( img ); % returns a function Sf(k) that builds a system matrix for 'k'
0060 pp.FC= system_mat_fields( fwd_model );
0061 pp.FT= system_mat_2p5d_fields( fwd_model );
0062 
0063 c2f = speye(pp.n_elem); try c2f = img.fwd_model.coarse2fine; end % make code uniform by assigning 1:1 c2f if no c2f provided
0064 k = [0 Inf]; try k = img.fwd_model.jacobian_adjoint_2p5d_1st_order.k; end
0065 method = 'quadv'; try method = img.fwd_model.jacobian_adjoint_2p5d_1st_order.method; end
0066 
0067 if length(k) == 1 % singleton k
0068    J = 2*jacobian_k(k, pp, gnd, img.fwd_model.stimulation, c2f);
0069 else
0070    switch method
0071       case 'trapz'
0072          % less accurate: trapz
0073          trace = 0;
0074          if trace; fprintf('%8s %12s %16s %16s\n','fcnt','a','b-a','||Q||'); end
0075          n = 0; kil = k(1);
0076          tol = 1e-8;
0077          k(isinf(k)) = tol^(-1/6);
0078          Jf = zeros(pp.n_meas, size(c2f,2), length(k)); % voltages under electrodes elec x stim, frequency domain
0079          for ki = k
0080             n = n + 1;
0081             Jf(:,:,n) = jacobian_k(ki, pp, gnd, img.fwd_model.stimulation, c2f);
0082             if trace; fprintf('%8d     %12e %16e %16e\n',n,ki,ki-kil,norm(Jf(1,1,n))); kil = ki; end
0083          end
0084          J=2/pi*trapz(k,Jf,3);
0085          % check
0086          Jff = squeeze(reshape(Jf,pp.n_meas*size(c2f,2),1,length(k)));
0087          assert(max(abs(Jff(:,end)) < tol), sprintf('trapz k=%e truncated too early as k->Inf',k(end)));
0088          slope = (Jff(:,3)-Jff(:,2))./Jff(:,2); slope(Jff(:,2) == 0) = nan;
0089          %clf; plot(slope); drawnow;
0090          %[idx]=find(abs(slope) > median(abs(slope))+sqrt(tol)); clf;semilogx(k,Jff(idx,:)); drawnow;
0091          assert(median(abs(slope)) < sqrt(tol)*10, sprintf('trapz k->0 (tol=%e), slope=%e != 0',tol,median(abs(slope))));
0092          if 0 % draw J(k) to check we are integrating over a large enough range
0093             clf;semilogx(repmat(k,pp.n_meas*size(c2f,2),1)',Jff');
0094             xlabel('k'); ylabel('J_k'); title('J_k')
0095             drawnow; pause(1);
0096          end
0097          clear Jff;
0098          clear slope;
0099       case 'quadv'
0100          % more accurate: adaptive gaussian quadrature
0101          trace = 0;
0102          tol = norm(jacobian_k(0, pp, gnd, img.fwd_model.stimulation, 1))*1e-2;
0103          kend = min(tol^(-1/6), max(k)); % don't go too far... k=Inf is a singular matrix, stop adjacent to numeric singularity
0104          % quadv is scheduled to be removed from matlab eventually... but it is
0105          % WAY faster than integral with any tolerance configuration I could identify
0106          if trace; fprintf('%8s %12s %16s %16s\n','fcnt','a','b-a','Q(1)'); end
0107          J=2/pi*quadv(@(kk) jacobian_k(kk, pp, gnd, img.fwd_model.stimulation, c2f), k(1), kend, tol, trace);
0108       case 'integral'
0109          reltol = 1e-8;
0110          tol = norm(jacobian_k(0, pp, gnd, img.fwd_model.stimulation, 1))*reltol;
0111          kend = min(tol^(-1/6), max(k)); % don't go too far... k=Inf is a singular matrix, stop adjacent to numeric singularity
0112          opts = {'ArrayValued', true,
0113                  'AbsTol', tol, % default: 1e-10
0114                  'RelTol',reltol}; % default:  1e-6
0115          opts = opts';
0116          % the integral solution is about 10x slower (5.47 seconds vs. 0.60 seconds for UNIT_TEST)
0117          % ... I played with AbsTol and RelTol but wasn't able to affect the outcome
0118          J=2/pi*integral(@(kk) jacobian_k(kk, pp, gnd, img.fwd_model.stimulation, c2f), k(1), kend, 'ArrayVAlued', true);
0119       otherwise
0120          error(['unrecognized method: ' method]);
0121    end
0122 end
0123 
0124 if ~strcmp(org_params,'conductivity')
0125     J = apply_chain_rule(J, img, org_params);
0126     if isfield(fwd_model, 'coarse2fine') && ...
0127           size(img.elem_data,1)==size(fwd_model.coarse2fine,1)
0128             J=J*fwd_model.coarse2fine;
0129             nparam = size(fwd_model.coarse2fine,2);
0130     end
0131 end
0132 
0133 %restore img to original condition
0134 if ~strcmp(org_params,'conductivity') || isfield(img, org_params)
0135     img = rmfield(img,'elem_data');
0136     img.current_params = [];
0137 end
0138 
0139 % calculate normalized Jacobian
0140 if pp.normalize
0141    data= fwd_solve( img );
0142    J= J ./ (data.meas(:)*ones(1,nparam));
0143 
0144 end
0145 
0146 function J = jacobian_k(k, pp, gnd, stim, c2f)
0147    SS = pp.Sf(k);
0148    idx = 1:size(SS,1);
0149    idx( gnd ) = [];
0150    
0151    sv = potential_k(SS, pp, gnd);
0152    
0153    zi2E= zeros(pp.n_elec, pp.n_node);
0154    zi2E(:, idx)= -pp.N2E(:,idx)/ SS(idx,idx);
0155    
0156    % NOTE: this is k*FT here, not k^2*FT: FT is squared in jacobian_calc (FT*DE*FT)
0157    DE = jacobian_calc(pp, zi2E, pp.FC, k*pp.FT, sv, c2f);
0158    
0159    J = zeros( pp.n_meas, size(c2f,2) );
0160    idx=0;
0161    for j= 1:pp.n_stim
0162       meas_pat= stim(j).meas_pattern;
0163       n_meas  = size(meas_pat,1);
0164       DEj = reshape( DE(:,j,:), pp.n_elec, size(c2f,2) );
0165       J( idx+(1:n_meas),: ) = meas_pat*DEj;
0166       idx= idx+ n_meas;
0167    end
0168 
0169 function v = potential_k(S, pp, gnd)
0170    idx = 1:size(S,1);
0171    idx( gnd ) = [];
0172    idx2 = find(any(pp.N2E));
0173    v = zeros(pp.n_node, pp.n_stim); % voltages at all nodes x number of stim, frequency domain
0174    v(idx,:) = left_divide( S(idx,idx), 1/2*pp.QQ(idx,:) );
0175 
0176 % DE_{i,j,k} is dV_i,j / dS_k
0177 %  where V_i is change in voltage on electrode i for
0178 %        stimulation pattern j
0179 %        S_k is change in conductivity on element k
0180 function DE = jacobian_calc(pp, zi2E, FC, FT, sv, c2f)
0181    d0 = pp.n_dims + 0;
0182    d1 = pp.n_dims + 1;
0183    
0184    zi2E_FCt = zi2E * FC';
0185    zi2E_FTt = zi2E * FT';
0186    
0187    FC_sv = FC * sv;
0188    FT_sv = FT * sv;
0189 
0190    DE= zeros(pp.n_elec, pp.n_stim, size(c2f,2) );
0191    de0= pp.n_elem * d0;
0192    de1= pp.n_elem * d1;
0193    if all(all(speye(size(c2f)) == c2f))
0194       for k= 1:size(c2f,2);
0195          idx = d0*(k-1)+1 : d0*k;
0196          dq1= zi2E_FCt(:,idx) * FC_sv(idx,:);
0197          dq2= zi2E_FTt(:,idx) * FT_sv(idx,:);
0198          DE(:,:,k)= dq1 + dq2;
0199       end
0200    elseif 0 % Code is slower
0201       for k= 1:size(c2f,2);
0202          chg_col = kron( c2f(:,k), ones(d0,1));
0203          dDD_dEj = spdiags(chg_col, 0, de0, de0);
0204          dq1 = zi2E_FCt * dDD_dEj * FC_sv;
0205          chg_col = kron( c2f(:,k), ones(d1,1));
0206          dDD_dEj = spdiags(chg_col, 0, de1, de1);
0207          dq2 = zi2E_FTt * dDD_dEj * FT_sv;
0208          DE(:,:,k)= dq1 + dq2;
0209       end
0210    else
0211       for k= 1:size(c2f,2);
0212          ff = find( c2f(:,k) > 1e-8 );
0213          lff= length(ff)*d0;
0214          ff1= ones(d0,1) * ff(:)';
0215          ffd= d0*ff1 + (-(d0-1):0)'*ones(1,length(ff));
0216          dDD_dEj = spdiags(c2f(ff1,k), 0, lff, lff);
0217          dq1= zi2E_FCt(:,ffd) * dDD_dEj * FC_sv(ffd,:);
0218          lff= length(ff)*d1;
0219          ff1= ones(d1,1) * ff(:)';
0220          ffd= d1*ff1 + (-(d1-1):0)'*ones(1,length(ff));
0221          dDD_dEj = spdiags(c2f(ff1,k), 0, lff, lff);
0222          dq2= zi2E_FTt(:,ffd) * dDD_dEj * FT_sv(ffd,:);
0223          DE(:,:,k)= dq1 + dq2;
0224       end
0225    end
0226 
0227 function J = apply_chain_rule(J, img, org_params)
0228    switch(org_params)
0229       case 'resistivity'
0230          dCond_dPhys = -img.elem_data.^2;
0231       case 'log_resistivity'
0232          dCond_dPhys = -img.elem_data;
0233       case 'log_conductivity'
0234          dCond_dPhys = img.elem_data;
0235       otherwise
0236          error('not implemented yet')
0237    end
0238    J = J.*repmat(dCond_dPhys ,1,size(J,1))';
0239 
0240 function elem_data = check_elem_data(img)
0241    elem_data = img.elem_data;
0242    sz_elem_data = size(elem_data);
0243    if sz_elem_data(2) ~= 1;
0244       error('jacobian_adjoin: can only solve one image (sz_elem_data=%)', ...
0245             sz_elem_data);
0246    end
0247 
0248    if isfield(img.fwd_model, 'coarse2fine');
0249       c2f = img.fwd_model.coarse2fine;
0250       sz_c2f = size(c2f);
0251       switch sz_elem_data(1)
0252          case sz_c2f(1); % Ok
0253          case sz_c2f(2); elem_data = c2f * elem_data;
0254          otherwise; error(['jacobian_adjoint: provided elem_data ' ...
0255               ' (sz=%d) does not match c2f (sz=%d %d)'], sz_elem_data(1), sz_c2f);
0256       end
0257    else
0258       if sz_elem_data(1) ~= num_elems(img.fwd_model)
0259          error(['jacobian_adjoint: provided elem_data (sz=%d) does ' ...
0260             ' not match fwd_model (sz=%d)'], sz_elem_data(1), num_elems(sz_c2f));
0261       end
0262    end
0263 
0264 function str = supported_params
0265     str = {'conductivity'
0266            'resistivity'
0267            'log_conductivity'
0268            'log_resistivity'};
0269 
0270 function do_unit_test()
0271    imdl2 = mk_geophysics_model('h22c',16);
0272    assert(length(imdl2.rec_model.elems) > 15, 'expect sufficient rec_model density');
0273    img2 = mk_image(imdl2,1);
0274    if 0 % check model c2f
0275       clf; subplot(121); show_fem(imdl2.fwd_model); subplot(122); show_fem(imdl2.rec_model);
0276       c2f=imdl2.fwd_model.coarse2fine; img2r = mk_image(imdl2); img2r.fwd_model = imdl2.rec_model; img2r.elem_data = zeros(size(c2f,2),1);
0277       for i=1:size(c2f,2)
0278          subplot(121);
0279          img2.elem_data = c2f(:,i); show_fem(img2); title(sprintf('%d',i));
0280          subplot(122);
0281          img2r.elem_data(:) = 0; img2r.elem_data(i)=1; show_fem(img2r);  title(sprintf('%d',i));
0282          drawnow; pause(0.25);
0283       end
0284    end
0285 
0286    % for the 3d model, we throw out the rec_model and inject the
0287    % imdl2.rec_model, then recalculate the c2f so that we can compare apples to
0288    % apples when we take ||Jxxx - J3||
0289    imdl3 = mk_geophysics_model('h32c',16);
0290    for s = {'nodes', 'elems', 'boundary', 'name','electrode'}
0291       imdl3.rec_model.(s{:}) = imdl2.rec_model.(s{:});
0292    end
0293    disp('recalculating 3D c2f, so that coarse meshes agree (2D vs 3D)');
0294    [c2f,bkgnd] = mk_approx_c2f(imdl3.fwd_model, imdl3.rec_model);
0295    imdl3.fwd_model.coarse2fine = c2f;
0296    imdl3.fwd_model.background = bkgnd;
0297    img3 = mk_image(imdl3,1);
0298    if 0 % check model c2f
0299       clf; subplot(121); show_fem(imdl3.fwd_model); subplot(122); show_fem(imdl3.rec_model);
0300       c2f=imdl3.fwd_model.coarse2fine; img3r = mk_image(imdl3); img3r.fwd_model = imdl3.rec_model; img3r.elem_data = zeros(size(c2f,2),1);
0301       for i=1:size(c2f,2)
0302          subplot(121);
0303          img3.elem_data = c2f(:,i); show_fem(img3); title(sprintf('%d',i));
0304          subplot(122);
0305          img3r.elem_data(:) = 0; img3r.elem_data(i)=1; show_fem(img3r);  title(sprintf('%d',i));
0306          drawnow; pause(0.25);
0307       end
0308    end
0309 
0310    t = tic;
0311    eidors_cache clear_name jacobian
0312    J2 = calc_jacobian(img2);
0313    fprintf(' 2D                         = %.2f sec\n', toc(t));
0314 
0315    t = tic;
0316    eidors_cache clear_name jacobian
0317    img2.fwd_model.jacobian = @jacobian_adjoint_2p5d_1st_order;
0318    img2.fwd_model.jacobian_adjoint_2p5d_1st_order.k = 0;
0319    J2p50 = calc_jacobian(img2);
0320    fprintf(' 2.5D (k=0)                 = %.2f sec\n', toc(t));
0321 
0322    ke = [0 0 0];
0323    t = tic;
0324    eidors_cache clear_name jacobian
0325    img2.fwd_model.jacobian = @jacobian_adjoint_2p5d_1st_order;
0326    img2.fwd_model.jacobian_adjoint_2p5d_1st_order.k = [0 logspace(-4,1,100)]; % capture all the effects to (k^2 T)!
0327    img2.fwd_model.jacobian_adjoint_2p5d_1st_order.method = 'trapz';
0328    J2p5kt = calc_jacobian(img2);
0329    ke(1) = img2.fwd_model.jacobian_adjoint_2p5d_1st_order.k(end);
0330    fprintf(' 2.5D (k=0..%.1f, trapz)    = %.2f sec\n', ke(1), toc(t));
0331    t = tic;
0332    img2.fwd_model.nodes(1,:) = img2.fwd_model.nodes(1,:) + rand(1,2)*1e-8; % defeat cache
0333    img2.fwd_model.jacobian = @jacobian_adjoint_2p5d_1st_order;
0334    img2.fwd_model.jacobian_adjoint_2p5d_1st_order.k = [0 Inf];
0335    img2.fwd_model.jacobian_adjoint_2p5d_1st_order.method = 'quadv';
0336    J2p5kq = calc_jacobian(img2);
0337    ke(2) = img2.fwd_model.jacobian_adjoint_2p5d_1st_order.k(end);
0338    fprintf(' 2.5D (k=0..%.1f, quadv)    = %.2f sec\n', ke(2), toc(t));
0339    t = tic;
0340    eidors_cache clear_name jacobian
0341    img2.fwd_model.jacobian = @jacobian_adjoint_2p5d_1st_order;
0342    img2.fwd_model.jacobian_adjoint_2p5d_1st_order.k = [0 Inf];
0343    img2.fwd_model.jacobian_adjoint_2p5d_1st_order.method = 'integral';
0344    J2p5ki = calc_jacobian(img2);
0345    ke(3) = img2.fwd_model.jacobian_adjoint_2p5d_1st_order.k(end);
0346    fprintf(' 2.5D (k=0..%.1f, integral) = %.2f sec\n', ke(3), toc(t));
0347 
0348    t = tic;
0349    eidors_cache clear_name jacobian
0350    J3 = calc_jacobian(img3);
0351    fprintf(' 3D                         = %.2f sec\n', toc(t));
0352 
0353 
0354    tol = 1e-8;
0355    reltol = norm(J3)/2; %*2e-2;
0356    unit_test_cmp('2D                          vs 3D', J2, J3, -tol);
0357    unit_test_cmp('2.5D (k=0)                  vs 2D', J2p50, J2, tol);
0358    unit_test_cmp(sprintf('2.5D (k=0..%-4.1f) (trapz)    vs 3D',ke(1)), J2p5kt, J3, 2*reltol);
0359    unit_test_cmp(sprintf('2.5D (k=0..%-4.1f) (quadv)    vs 3D',ke(2)), J2p5kq, J3, reltol);
0360    unit_test_cmp(sprintf('2.5D (k=0..%-4.1f) (integral) vs 3D',ke(3)), J2p5ki, J3, reltol);
0361 
0362    imgr = img2;
0363    imgr.fwd_model = imdl2.rec_model;
0364    clf; subplot(221); imgr.elem_data = sens(J3); show_fem(imgr,1); title('3D [log_{10}]');
0365         subplot(222); imgr.elem_data = sens(J2); show_fem(imgr,1); title('2D [log_{10}]');
0366         subplot(223); imgr.elem_data = sens(J2p5kq); show_fem(imgr,1); title('2.5D (k=0..3.0) [log_{10}]');
0367         subplot(224); imgr.elem_data = sens(J2p50); show_fem(imgr,1); title('2.5D (k=0) [log_{10}]');
0368    if 1
0369         subplot(222); imgr.elem_data = sens(J2p5kt); show_fem(imgr,1); title('2.5D (k=0..3.0) trapz [log_{10}]');
0370         subplot(223); imgr.elem_data = sens(J2p5kq); show_fem(imgr,1); title('2.5D (k=0..Inf) quadv [log_{10}]');
0371         subplot(224); imgr.elem_data = sens(J2p5ki); show_fem(imgr,1); title('2.5D (k=0..Inf) integral [log_{10}]');
0372    end
0373 
0374 function S = sens(J)
0375    S = sum(J.^2,1);
0376    S = log10(max(S,1e-10)'); % AVOID -inf