Description of jacobian

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 5608 2017-06-23 15:01:09Z htregidgo $
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 % changed behaviour to only factorise system matrix once. only if system
0051 % matrix is symmetric.
0052 if isequal((s_mat.E(idx,idx)).',s_mat.E(idx,idx))
0053     Einvtemp=left_divide(s_mat.E(idx,idx) , [pp.QQ( idx,: ),(-pp.N2E(:,idx)).']);
0054     
0055     sv= zeros(n, pp.n_stim );
0056     % sv( idx,:) = left_divide(s_mat.E(idx,idx) , pp.QQ( idx,: ));
0057     sv( idx,:) = Einvtemp(:,1:pp.n_stim);
0058     
0059     zi2E= zeros(pp.n_elec, n);
0060     % the minus below used to be missing
0061     % zi2E(:, idx)= -pp.N2E(:,idx)/ s_mat.E(idx,idx) ;
0062     zi2E(:, idx)= (Einvtemp(:,pp.n_stim+1:end)).';
0063 else
0064     sv= zeros(n, pp.n_stim );
0065     sv( idx,:) = left_divide(s_mat.E(idx,idx) , pp.QQ( idx,: ));
0066     
0067     zi2E= zeros(pp.n_elec, n);
0068     % the minus below used to be missing
0069     zi2E(:, idx)= -pp.N2E(:,idx)/ s_mat.E(idx,idx) ;
0070 end
0071 
0072 FC= system_mat_fields( fwd_model );
0073 
0074 if isfield(fwd_model,'coarse2fine') && strcmp(org_params, 'conductivity');
0075    DE = jacobian_calc(pp, zi2E, FC, sv, fwd_model.coarse2fine);
0076    nparam= size(fwd_model.coarse2fine,2);
0077 else
0078    DE = jacobian_calc(pp, zi2E, FC, sv);
0079    nparam= e;
0080 end
0081 
0082 J = zeros( pp.n_meas, nparam );
0083 idx=0;
0084 for j= 1:pp.n_stim
0085    meas_pat= fwd_model.stimulation(j).meas_pattern;
0086    n_meas  = size(meas_pat,1);
0087    DEj = reshape( DE(:,j,:), pp.n_elec, nparam );
0088    J( idx+(1:n_meas),: ) = meas_pat*DEj;
0089    idx= idx+ n_meas;
0090 end
0091 
0092 if ~strcmp(org_params,'conductivity')
0093     J = apply_chain_rule(J, img, org_params);
0094     if isfield(fwd_model, 'coarse2fine') && ...
0095           size(img.elem_data,1)==size(fwd_model.coarse2fine,1)
0096             J=J*fwd_model.coarse2fine;
0097             nparam = size(fwd_model.coarse2fine,2);
0098     end
0099 end
0100 
0101 %restore img to original condition
0102 if ~strcmp(org_params,'conductivity') || isfield(img, org_params)
0103     img = rmfield(img,'elem_data');
0104     img.current_params = [];
0105 end
0106 
0107 % calculate normalized Jacobian
0108 if pp.normalize
0109    data= fwd_solve( img );
0110    J= J ./ (data.meas(:)*ones(1,nparam));
0111    
0112 end
0113 
0114 % This was a correction for the missing minus above
0115 % J= -J;
0116 
0117 % DE_{i,j,k} is dV_i,j / dS_k
0118 %  where V_i is change in voltage on electrode i for
0119 %        stimulation pattern j
0120 %        S_k is change in conductivity on element k
0121 function DE = jacobian_calc(pp, zi2E, FC, sv, c2f)
0122 d= pp.n_dims+1;
0123 dfact= factorial(d);
0124 
0125 do_c2f = ( nargin==5 );
0126 zi2E_FCt = zi2E * FC';
0127 
0128 FC_sv   = FC * sv;
0129 
0130 if ~do_c2f
0131    DE= zeros(pp.n_elec, pp.n_stim, pp.n_elem);
0132    for k= 1:pp.n_elem
0133        idx= (d-1)*(k-1)+1 : (d-1)*k;
0134        dq= zi2E_FCt(:,idx) * FC_sv(idx,:);
0135        DE(:,:,k)= dq;
0136    end
0137 else
0138    DE= zeros(pp.n_elec, pp.n_stim, size(c2f,2) );
0139    if 0 % Code is slower
0140       de= pp.n_elem * (d-1);
0141       for k= 1:size(c2f,2);
0142           chg_col = kron( c2f(:,k), ones(d-1,1));
0143           dDD_dEj = spdiags(chg_col,0, de, de);
0144           dq= zi2E_FCt * dDD_dEj * FC_sv;
0145           DE(:,:,k)= dq;
0146       end
0147    else
0148       de= pp.n_elem * (d-1);
0149       for k= 1:size(c2f,2);
0150           ff = find( c2f(:,k) );
0151           lff= length(ff)*(d-1);
0152           ff1= ones(d-1,1) * ff(:)';
0153           ffd= (d-1)*ff1 + (-(d-2):0)'*ones(1,length(ff));
0154           dDD_dEj = spdiags(c2f(ff1,k), 0, lff, lff);
0155           dq= zi2E_FCt(:,ffd) * dDD_dEj * FC_sv(ffd,:);
0156           DE(:,:,k)= dq;
0157       end
0158    end
0159 end
0160 
0161 function J = apply_chain_rule(J, img, org_params)
0162 
0163 switch(org_params)
0164     case 'resistivity'
0165         dCond_dPhys = -img.elem_data.^2;
0166     case 'log_resistivity'
0167         dCond_dPhys = -img.elem_data;
0168     case 'log_conductivity'
0169         dCond_dPhys = img.elem_data;
0170     otherwise
0171         error('not implemented yet')
0172 end
0173 
0174 J = J.*repmat(dCond_dPhys ,1,size(J,1))';
0175 
0176 function elem_data = check_elem_data(img)
0177    elem_data = img.elem_data; 
0178    sz_elem_data = size(elem_data);
0179    if sz_elem_data(2) ~= 1;
0180       error('jacobian_adjoin: can only solve one image (sz_elem_data=%)', ...
0181             sz_elem_data);
0182    end
0183 
0184    if isfield(img.fwd_model, 'coarse2fine');
0185      c2f = img.fwd_model.coarse2fine;
0186      sz_c2f = size(c2f);
0187      switch sz_elem_data(1)
0188        case sz_c2f(1); % Ok
0189        case sz_c2f(2); elem_data = c2f * elem_data;
0190        otherwise; error(['jacobian_adjoint: provided elem_data ' ...
0191             ' (sz=%d) does not match c2f (sz=%d %d)'], sz_elem_data(1), sz_c2f);
0192      end
0193    else
0194      if sz_elem_data(1) ~= num_elems(img.fwd_model)
0195        error(['jacobian_adjoint: provided elem_data (sz=%d) does ' ...
0196           ' not match fwd_model (sz=%d)'], sz_elem_data(1), num_elems(sz_c2f));
0197      end
0198    end
0199 
0200 
0201 
0202 function str = supported_params
0203     str = {'conductivity'
0204            'resistivity'
0205            'log_conductivity'
0206            'log_resistivity'};
0207       
0208 
0209 function do_unit_test
0210    current = 4; measure=1;
0211    [R,img] = test_2d_resistor(current,measure);
0212    img.fwd_solve.get_all_nodes = 1;
0213    vs = fwd_solve_1st_order( img);
0214    va= measure*current*sum(R); % analytic
0215    unit_test_cmp('2D resistor test', va, vs.meas, 1e-12);
0216 
0217    J = jacobian_adjoint(img);
0218    unit_test_cmp('2D resistor Jacobian', size(J), ...
0219       [length(img.fwd_model.stimulation), size(img.fwd_model.coarse2fine,2)]);
0220    unit_test_cmp('2D resistor Jacobian', std(J),0, 1e-12);
0221 %  unit_test_cmp('2D R voltages', vs.volt(1:3:10)-vs.volt(1), ...
0222 
0223    [R,img] = test_3d_resistor(current,measure);
0224    img.fwd_solve.get_all_nodes = 1;
0225    vs = fwd_solve_1st_order( img);
0226    va= current*sum(R);
0227    unit_test_cmp('3D resistor test', va, vs.meas, 1e-10);
0228    J = jacobian_adjoint(img);
0229    unit_test_cmp('3D resistor Jacobian', size(J), ...
0230       [length(img.fwd_model.stimulation), size(img.fwd_model.coarse2fine,2)]);
0231    unit_test_cmp('3D resistor Jacobian', std(J),0, 1e-12);
0232 
0233 function [R,img] = test_2d_resistor(current,measure)
0234    conduc=  .4 + 2*pi*j*10; % conductivity in Ohm-meters
0235    z_contact= .1; wid = 3; len = 12; 
0236 
0237    fmdl=mk_grid_model([],linspace(0,wid,3), linspace(0,len,4));
0238    fmdl.electrode(1).nodes = find(fmdl.nodes(:,2) ==   0);
0239    fmdl.electrode(2).nodes = find(fmdl.nodes(:,2) == 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 / conduc;
0245    Contact_R = z_contact/wid;
0246    R = [Block_R, 2*Contact_R];
0247 
0248 function [R,img] = test_3d_resistor(current,measure);;
0249    conduc=  .4 + 2*pi*j*10; % conductivity in Ohm-meters
0250    z_contact= .1; wid = 2; len = 5; hig=3; 
0251 
0252    fmdl=mk_grid_model([],0:wid, 0:hig, 0:len);
0253    fmdl.electrode(1).nodes = find(fmdl.nodes(:,3) ==   0);
0254    fmdl.electrode(2).nodes = find(fmdl.nodes(:,3) == len);
0255    [fmdl.electrode(:).z_contact] = deal(z_contact);
0256    fmdl.stimulation = stim_meas_list([1,2,1,2],2,current,measure);
0257    img= mk_image(fmdl,conduc);
0258 
0259    Block_R =  len / wid / hig / conduc;
0260    Contact_R = z_contact/(wid*hig);
0261    R = [Block_R, 2*Contact_R];
jacobian_adjoint

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
