jacobian_adjoint

PURPOSE ^

JACOBIAN_ADJOINT: J= jacobian_adjoint( img )

SYNOPSIS ^

function J= jacobian_adjoint( fwd_model, img)

DESCRIPTION ^

 JACOBIAN_ADJOINT: J= jacobian_adjoint( 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function J= jacobian_adjoint( fwd_model, img)
0002 % JACOBIAN_ADJOINT: J= jacobian_adjoint( 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 
0013 % (C) 2005 Andy Adler. License: GPL version 2 or version 3
0014 % $Id: jacobian_adjoint.m 5466 2017-05-10 10:38:25Z aadler $
0015 
0016 if ischar(fwd_model) && strcmp(fwd_model,'UNIT_TEST'); do_unit_test; return; end
0017 
0018 if nargin == 1
0019    img= fwd_model;
0020 elseif  strcmp(getfield(warning('query','EIDORS:DeprecatedInterface'),'state'),'on')
0021    warning('EIDORS:DeprecatedInterface', ...
0022       ['Calling JACOBIAN_ADJOINT with two arguments is deprecated and will cause' ...
0023        ' an error in a future version. First argument ignored.']);
0024 end
0025 
0026 img = data_mapper(img);
0027 if ~ismember(img.current_params, supported_params)
0028     error('EIDORS:PhysicsNotSupported', '%s does not support %s', ...
0029     'JACOBIAN_ADJOINT',img.current_params);
0030 end
0031 
0032 org_params = img.current_params;
0033 % all calcs use conductivity
0034 img = convert_img_units(img, 'conductivity');
0035 
0036 img.elem_data = check_elem_data(img);
0037 
0038 fwd_model= img.fwd_model;
0039 
0040 pp= fwd_model_parameters( fwd_model, 'skip_VOLUME' );
0041 s_mat= calc_system_mat( img );
0042 
0043 d= pp.n_dims+1;
0044 e= pp.n_elem;
0045 n= pp.n_node;
0046 
0047 idx= 1:size(s_mat.E,1);
0048 idx( fwd_model.gnd_node ) = [];
0049 
0050 sv= zeros(n, pp.n_stim );
0051 sv( idx,:) = left_divide(s_mat.E(idx,idx) , pp.QQ( idx,: ));
0052 
0053 zi2E= zeros(pp.n_elec, n);
0054 % the minus below used to be missing
0055 zi2E(:, idx)= -pp.N2E(:,idx)/ s_mat.E(idx,idx) ;
0056 
0057 FC= system_mat_fields( fwd_model );
0058 
0059 if isfield(fwd_model,'coarse2fine') && strcmp(org_params, 'conductivity');
0060    DE = jacobian_calc(pp, zi2E, FC, sv, fwd_model.coarse2fine);
0061    nparam= size(fwd_model.coarse2fine,2);
0062 else
0063    DE = jacobian_calc(pp, zi2E, FC, sv);
0064    nparam= e;
0065 end
0066 
0067 J = zeros( pp.n_meas, nparam );
0068 idx=0;
0069 for j= 1:pp.n_stim
0070    meas_pat= fwd_model.stimulation(j).meas_pattern;
0071    n_meas  = size(meas_pat,1);
0072    DEj = reshape( DE(:,j,:), pp.n_elec, nparam );
0073    J( idx+(1:n_meas),: ) = meas_pat*DEj;
0074    idx= idx+ n_meas;
0075 end
0076 
0077 if ~strcmp(org_params,'conductivity')
0078     J = apply_chain_rule(J, img, org_params);
0079     if isfield(fwd_model, 'coarse2fine') && ...
0080           size(img.elem_data,1)==size(fwd_model.coarse2fine,1)
0081             J=J*fwd_model.coarse2fine;
0082             nparam = size(fwd_model.coarse2fine,2);
0083     end
0084 end
0085 
0086 %restore img to original condition
0087 if ~strcmp(org_params,'conductivity') || isfield(img, org_params)
0088     img = rmfield(img,'elem_data');
0089     img.current_params = [];
0090 end
0091 
0092 % calculate normalized Jacobian
0093 if pp.normalize
0094    data= fwd_solve( img );
0095    J= J ./ (data.meas(:)*ones(1,nparam));
0096    
0097 end
0098 
0099 % This was a correction for the missing minus above
0100 % J= -J;
0101 
0102 % DE_{i,j,k} is dV_i,j / dS_k
0103 %  where V_i is change in voltage on electrode i for
0104 %        stimulation pattern j
0105 %        S_k is change in conductivity on element k
0106 function DE = jacobian_calc(pp, zi2E, FC, sv, c2f)
0107 d= pp.n_dims+1;
0108 dfact= factorial(d);
0109 
0110 do_c2f = ( nargin==5 );
0111 zi2E_FCt = zi2E * FC';
0112 
0113 FC_sv   = FC * sv;
0114 
0115 if ~do_c2f
0116    DE= zeros(pp.n_elec, pp.n_stim, pp.n_elem);
0117    for k= 1:pp.n_elem
0118        idx= (d-1)*(k-1)+1 : (d-1)*k;
0119        dq= zi2E_FCt(:,idx) * FC_sv(idx,:);
0120        DE(:,:,k)= dq;
0121    end
0122 else
0123    DE= zeros(pp.n_elec, pp.n_stim, size(c2f,2) );
0124    if 0 % Code is slower
0125       de= pp.n_elem * (d-1);
0126       for k= 1:size(c2f,2);
0127           chg_col = kron( c2f(:,k), ones(d-1,1));
0128           dDD_dEj = spdiags(chg_col,0, de, de);
0129           dq= zi2E_FCt * dDD_dEj * FC_sv;
0130           DE(:,:,k)= dq;
0131       end
0132    else
0133       de= pp.n_elem * (d-1);
0134       for k= 1:size(c2f,2);
0135           ff = find( c2f(:,k) );
0136           lff= length(ff)*(d-1);
0137           ff1= ones(d-1,1) * ff(:)';
0138           ffd= (d-1)*ff1 + (-(d-2):0)'*ones(1,length(ff));
0139           dDD_dEj = spdiags(c2f(ff1,k), 0, lff, lff);
0140           dq= zi2E_FCt(:,ffd) * dDD_dEj * FC_sv(ffd,:);
0141           DE(:,:,k)= dq;
0142       end
0143    end
0144 end
0145 
0146 function J = apply_chain_rule(J, img, org_params)
0147 
0148 switch(org_params)
0149     case 'resistivity'
0150         dCond_dPhys = -img.elem_data.^2;
0151     case 'log_resistivity'
0152         dCond_dPhys = -img.elem_data;
0153     case 'log_conductivity'
0154         dCond_dPhys = img.elem_data;
0155     otherwise
0156         error('not implemented yet')
0157 end
0158 
0159 J = J.*repmat(dCond_dPhys ,1,size(J,1))';
0160 
0161 function elem_data = check_elem_data(img)
0162    elem_data = img.elem_data; 
0163    sz_elem_data = size(elem_data);
0164    if sz_elem_data(2) ~= 1;
0165       error('jacobian_adjoin: can only solve one image (sz_elem_data=%)', ...
0166             sz_elem_data);
0167    end
0168 
0169    if isfield(img.fwd_model, 'coarse2fine');
0170      c2f = img.fwd_model.coarse2fine;
0171      sz_c2f = size(c2f);
0172      switch sz_elem_data(1)
0173        case sz_c2f(1); % Ok
0174        case sz_c2f(2); elem_data = c2f * elem_data;
0175        otherwise; error(['jacobian_adjoint: provided elem_data ' ...
0176             ' (sz=%d) does not match c2f (sz=%d %d)'], sz_elem_data(1), sz_c2f);
0177      end
0178    else
0179      if sz_elem_data(1) ~= num_elems(img.fwd_model)
0180        error(['jacobian_adjoint: provided elem_data (sz=%d) does ' ...
0181           ' not match fwd_model (sz=%d)'], sz_elem_data(1), num_elems(sz_c2f));
0182      end
0183    end
0184 
0185 
0186 
0187 function str = supported_params
0188     str = {'conductivity'
0189            'resistivity'
0190            'log_conductivity'
0191            'log_resistivity'};
0192       
0193 
0194 function do_unit_test
0195    current = 4; measure=1;
0196    [R,img] = test_2d_resistor(current,measure);
0197    img.fwd_solve.get_all_nodes = 1;
0198    vs = fwd_solve_1st_order( img);
0199    va= measure*current*sum(R); % analytic
0200    unit_test_cmp('2D resistor test', va, vs.meas, 1e-12);
0201 
0202    J = jacobian_adjoint(img);
0203    unit_test_cmp('2D resistor Jacobian', size(J), ...
0204       [length(img.fwd_model.stimulation), size(img.fwd_model.coarse2fine,2)]);
0205    unit_test_cmp('2D resistor Jacobian', std(J),0, 1e-12);
0206 %  unit_test_cmp('2D R voltages', vs.volt(1:3:10)-vs.volt(1), ...
0207 
0208    [R,img] = test_3d_resistor(current,measure);
0209    img.fwd_solve.get_all_nodes = 1;
0210    vs = fwd_solve_1st_order( img);
0211    va= current*sum(R);
0212    unit_test_cmp('3D resistor test', va, vs.meas, 1e-10);
0213    J = jacobian_adjoint(img);
0214    unit_test_cmp('3D resistor Jacobian', size(J), ...
0215       [length(img.fwd_model.stimulation), size(img.fwd_model.coarse2fine,2)]);
0216    unit_test_cmp('3D resistor Jacobian', std(J),0, 1e-12);
0217 
0218 function [R,img] = test_2d_resistor(current,measure)
0219    conduc=  .4 + 2*pi*j*10; % conductivity in Ohm-meters
0220    z_contact= .1; wid = 3; len = 12; 
0221 
0222    fmdl=mk_grid_model([],linspace(0,wid,3), linspace(0,len,4));
0223    fmdl.electrode(1).nodes = find(fmdl.nodes(:,2) ==   0);
0224    fmdl.electrode(2).nodes = find(fmdl.nodes(:,2) == len);
0225    [fmdl.electrode(:).z_contact] = deal(z_contact);
0226    fmdl.stimulation = stim_meas_list([1,2,1,2],2,current,measure);
0227    img= mk_image(fmdl,conduc);
0228 
0229    Block_R = len / wid / conduc;
0230    Contact_R = z_contact/wid;
0231    R = [Block_R, 2*Contact_R];
0232 
0233 function [R,img] = test_3d_resistor(current,measure);;
0234    conduc=  .4 + 2*pi*j*10; % conductivity in Ohm-meters
0235    z_contact= .1; wid = 2; len = 5; hig=3; 
0236 
0237    fmdl=mk_grid_model([],0:wid, 0:hig, 0:len);
0238    fmdl.electrode(1).nodes = find(fmdl.nodes(:,3) ==   0);
0239    fmdl.electrode(2).nodes = find(fmdl.nodes(:,3) == len);
0240    [fmdl.electrode(:).z_contact] = deal(z_contact);
0241    fmdl.stimulation = stim_meas_list([1,2,1,2],2,current,measure);
0242    img= mk_image(fmdl,conduc);
0243 
0244    Block_R =  len / wid / hig / conduc;
0245    Contact_R = z_contact/(wid*hig);
0246    R = [Block_R, 2*Contact_R];

Generated on Wed 21-Jun-2017 09:29:07 by m2html © 2005