calc_jacobian

 CALC_JACOBIAN: calculate jacobian from an inv_model
 
  J = calc_jacobian( img )
      calc Jacobian on img.fwd_model at conductivity given
      in image (fwd_model is for forward and reconstruction)
 
 The actual work is done by the jacobian calculator specified in 
    img.fwd_model.jacobian, unless that field is numeric, in which case
    calc_jacobian returns its contents.

 For reconstructions on dual meshes, the interpolation matrix
    is defined as img.fwd_model.coarse2fine. This takes
    coarse2fine * x_coarse = x_fine

 If the underlying jacobian calculator doesn't understand dual
    meshes, then calc_jacobian will automatically postmultiply
    by fwd_model.coarse2fine.

 img       is an image structure, with 'elem_data' or
           'node_data' parameters

0001 function J = calc_jacobian( fwd_model, img)
0002 % CALC_JACOBIAN: calculate jacobian from an inv_model
0003 %
0004 %  J = calc_jacobian( img )
0005 %      calc Jacobian on img.fwd_model at conductivity given
0006 %      in image (fwd_model is for forward and reconstruction)
0007 %
0008 % The actual work is done by the jacobian calculator specified in
0009 %    img.fwd_model.jacobian, unless that field is numeric, in which case
0010 %    calc_jacobian returns its contents.
0011 %
0012 % For reconstructions on dual meshes, the interpolation matrix
0013 %    is defined as img.fwd_model.coarse2fine. This takes
0014 %    coarse2fine * x_coarse = x_fine
0015 %
0016 % If the underlying jacobian calculator doesn't understand dual
0017 %    meshes, then calc_jacobian will automatically postmultiply
0018 %    by fwd_model.coarse2fine.
0019 %
0020 % img       is an image structure, with 'elem_data' or
0021 %           'node_data' parameters
0022 
0023 % (C) 2005-08 Andy Adler. License: GPL version 2 or version 3
0024 % $Id: calc_jacobian.m 6716 2024-04-02 17:59:16Z aadler $
0025 
0026 if ischar(fwd_model) && strcmp(fwd_model,'UNIT_TEST'); do_unit_test; return; end
0027 
0028 if nargin == 1
0029    img= fwd_model;
0030 else
0031    warning('EIDORS:DeprecatedInterface', ...
0032       ['Calling CALC_JACOBIAN with two arguments is deprecated and will cause' ...
0033        ' an error in a future version. First argument ignored.']);
0034 end
0035 ws = warning('query','EIDORS:DeprecatedInterface');
0036 warning off EIDORS:DeprecatedInterface
0037 
0038 try 
0039    fwd_model= img.fwd_model;
0040 catch
0041    error('CALC_JACOBIAN requires an eidors image structure');
0042 end
0043 
0044 fwd_model_check(fwd_model);
0045 
0046 if isnumeric(fwd_model.jacobian)             % we have the Jacobian matrix
0047    J = fwd_model.jacobian;
0048 else                                         % we need to calculate
0049    
0050    copt.cache_obj= jacobian_cache_params( fwd_model, img );
0051    copt.fstr = 'jacobian';
0052    try
0053        fwd_model.jacobian = str2func(fwd_model.jacobian);
0054    end
0055    J = eidors_cache(fwd_model.jacobian, {fwd_model, img}, copt);
0056    
0057 end
0058 
0059 if isfield(fwd_model,'coarse2fine')
0060    c2f= fwd_model.coarse2fine;
0061    if size(J,2)==size(c2f,1)
0062 %     calc_jacobian did not take into account the coarse2fine
0063       J=J*c2f;
0064    end
0065 end
0066 
0067 warning(ws.state, 'EIDORS:DeprecatedInterface');
0068 
0069         
0070 
0071 
0072 
0073 % Make the Jacobian only depend on
0074 function cache_obj= jacobian_cache_params( fwd_model, img );
0075    img = data_mapper(img);
0076    if isfield(img, 'elem_data')
0077       cache_obj = {fwd_model, img.elem_data, img.current_params};
0078    elseif isfield(img, 'node_data')
0079       cache_obj = {fwd_model, img.node_data, img.current_params};
0080    else
0081       error('calc_jacobian: execting elem_data or node_data in image');
0082    end
0083 
0084 function fwd_model_check(fmdl)
0085    try; if fmdl.jacobian_dont_check
0086       return
0087    end; end
0088    pp = fwd_model_parameters(fmdl); % they cache, so no problem
0089    if pp.n_elec == 0
0090        error('Cannot calculate Jacobian. No electrodes found.');
0091    end
0092    if pp.n_stim == 0
0093        error('Cannot calculate Jacobian. No stimulation found.');
0094    end
0095    if pp.n_meas == 0
0096        error('Cannot calculate Jacobian. No measurements found.');
0097    end
0098 
0099 function do_unit_test
0100     delta = 1e-6;
0101     testvec= [5,20,40,130];
0102     mdl= test_aa_mdl2;
0103     run_jacobian_test( mdl, delta, testvec );
0104     mdl = mdl_normalize(mdl,1); mdl.name = [mdl.name,':normalize'];
0105     run_jacobian_test( mdl, delta, testvec );
0106 
0107     mdl= test_aa_mdl3;
0108     run_jacobian_test( mdl, delta, testvec );
0109     mdl = mdl_normalize(mdl,1); mdl.name = [mdl.name,':normalize'];
0110     run_jacobian_test( mdl, delta, testvec );
0111 
0112     warning('off','EIDORS:deprecated');
0113     mdl= test_np_mdl;
0114     % Need lower tol because np models use a solver with lower tol
0115     run_jacobian_test( mdl, delta, testvec, 1e-3 );
0116     mdl = mdl_normalize(mdl,1); mdl.name = [mdl.name,':normalize'];
0117     run_jacobian_test( mdl, delta, testvec, 1e-3 );
0118 
0119 function run_jacobian_test( mdl, delta, testvec,tol ); 
0120     img= mk_image(mdl,1);
0121     homg_data=fwd_solve( img);
0122 
0123     J= calc_jacobian( img );
0124     if nargin<4;tol= 1e-5;end
0125     tol = tol*std(J(:));
0126 
0127     % J = dF/dx = [F(x+d)  - F(x)]/d
0128     sumdiff= 0;
0129     bkgnd_elem_data= img.elem_data;
0130     for testelem = testvec
0131        img.elem_data= bkgnd_elem_data;
0132        img.elem_data(testelem)= bkgnd_elem_data(testelem)+delta;
0133        inh_data=fwd_solve( img);
0134 
0135        if mdl_normalize(mdl)
0136           simJ= 1/delta* (inh_data.meas ./ homg_data.meas - 1);
0137        else
0138           simJ= 1/delta* (inh_data.meas-homg_data.meas);
0139        end
0140        
0141        plot([J(:,testelem) simJ]);
0142        sumdiff = sumdiff + std( J(:,testelem) - simJ );
0143     end
0144 
0145     dif= sumdiff/length(testvec);
0146     unit_test_cmp(['Jacobian calc: ',mdl.name],sumdiff/length(testvec),0,tol);
0147 
0148 function mdl= test_aa_mdl2;
0149     n_elec= 16;
0150     n_rings= 1;
0151     options = {'no_meas_current','no_rotate_meas'};
0152     params= mk_circ_tank(8, [], n_elec);
0153 
0154     params.stimulation= mk_stim_patterns(n_elec, n_rings, '{ad}','{ad}', ...
0155                                 options, 10);
0156     params.solve=      'fwd_solve_1st_order';
0157     params.system_mat= 'system_mat_1st_order';
0158     params.jacobian=   'jacobian_adjoint';
0159     params.normalize_measurements = 0;
0160     mdl = eidors_obj('fwd_model', params);
0161     mdl.name= 'AA_1996 mdl';
0162 
0163 
0164 
0165 function mdl= test_aa_mdl3;
0166     i_mdl = mk_common_model('b3cz2',16);
0167     mdl= i_mdl.fwd_model;
0168     mdl.name= 'AA_1996 mdl';
0169     mdl.solve=      'fwd_solve_1st_order';
0170     mdl.system_mat= 'system_mat_1st_order';
0171     mdl.jacobian=   'jacobian_adjoint';
0172     
0173     
0174 
0175 function mdl= test_np_mdl;
0176     i_mdl = mk_common_model('n3r2',[16,2]);
0177     mdl= i_mdl.fwd_model;
0178     mdl.name=     'NP_2003 mdl';
0179     mdl.solve=      'np_fwd_solve';
0180     mdl.system_mat= 'np_calc_system_mat';
0181     mdl.jacobian=   'np_calc_jacobian';
0182 
0183 function mdl= test_ms_mdl;
0184     i_mdl = mk_common_model('n3r2',16);
0185     mdl= i_mdl.fwd_model;
0186     mdl.name=       'MS_2005 mdl';
0187     mdl.solve=      'np_fwd_solve';
0188     mdl.system_mat= 'np_calc_system_mat';
0189     mdl.jacobian=   'ms_calc_jacobian';