perturb_jacobian

 PERTURB_JACOBIAN: J= perturb_jacobian( fwd_model, img)
 Calculate Jacobian Matrix, based on small perturbations
   in the forward model. This will tend to be slow, but
   should be best used to 'sanity check' other code

 J         = Jacobian matrix
 fwd_model = forward model
 fwd_model.perturb_jacobian.delta   - delta perturbation to use
 img = image background for jacobian calc

0001 function J= perturb_jacobian( fwd_model, img)
0002 % PERTURB_JACOBIAN: J= perturb_jacobian( fwd_model, img)
0003 % Calculate Jacobian Matrix, based on small perturbations
0004 %   in the forward model. This will tend to be slow, but
0005 %   should be best used to 'sanity check' other code
0006 %
0007 % J         = Jacobian matrix
0008 % fwd_model = forward model
0009 % fwd_model.perturb_jacobian.delta   - delta perturbation to use
0010 % img = image background for jacobian calc
0011 
0012 % (C) 2006 Andy Adler. License: GPL version 2 or version 3
0013 % $Id: perturb_jacobian.html 2819 2011-09-07 16:43:11Z aadler $
0014 
0015 if isstr(fwd_model) && strcmp(fwd_model,'UNIT_TEST'); do_unit_test; return;end
0016 
0017 if isfield(fwd_model,'perturb_jacobian')
0018    delta = fwd_model.perturb_jacobian.delta;
0019 else
0020    delta= 1e-6; % tests indicate this is a good value
0021 end
0022 
0023 n_elem = size(fwd_model.elems,1);
0024 % force image to use provided fwd_model
0025 img.fwd_model= fwd_model;
0026 
0027 % solve one time to get the size
0028 d0= fwd_solve( img );
0029 
0030 if isfield(img.fwd_model,'coarse2fine');
0031    Jcol= perturb_c2f(img, 1, delta, d0);
0032    Jrows= size(img.fwd_model.coarse2fine,2);
0033    J= zeros(length(Jcol), Jrows );
0034    J(:,1)= Jcol;
0035    for i=2:size(img.fwd_model.coarse2fine,2);
0036      J(:,i)= perturb_c2f(img, i, delta, d0);
0037      if rem(i,50)==0; fprintf('+'); end
0038    end
0039 else
0040    Jcol= perturb(img, 1, delta, d0);
0041 
0042    J= zeros(length(Jcol), n_elem);
0043    J(:,1)= Jcol;
0044 
0045    for i=2:n_elem
0046      J(:,i)= perturb(img, i, delta, d0);
0047    end
0048 end
0049 
0050 
0051 function Jcol= perturb( img, i, delta, d0)
0052    img.elem_data(i)= img.elem_data(i) + delta;
0053    di= fwd_solve( img );
0054    Jcol = (1/delta) * (di.meas - d0.meas);
0055 
0056 function Jcol= perturb_c2f( img, i, delta, d0)
0057    img.elem_data= img.elem_data + delta*img.fwd_model.coarse2fine(:,i);
0058    di= fwd_solve( img );
0059    Jcol = (1/delta) * (di.meas - d0.meas);
0060 
0061 
0062 function do_unit_test
0063 % Perturbation Jacobians
0064 % $Id: perturb_jacobian.html 2819 2011-09-07 16:43:11Z aadler $
0065 
0066 % imdl= mk_common_model('c2c2',16);
0067   imdl= mk_common_model('a3cr',16);
0068 % imdl= mk_common_model('n3r2',16);
0069   imdl.fwd_model.nodes = imdl.fwd_model.nodes*.25;
0070   img= calc_jacobian_bkgnd(imdl);
0071 
0072   img.fwd_model.normalize_measurements= 0;
0073 
0074   img.fwd_model.jacobian=   @aa_calc_jacobian;
0075   img.fwd_model.system_mat= @aa_calc_system_mat;
0076   img.fwd_model.solve=      @aa_fwd_solve;
0077   J_aa= calc_jacobian( img ); % 2 for bug in my code
0078 
0079   img.fwd_model.jacobian=   @perturb_jacobian;
0080   img.fwd_model.system_mat= @aa_calc_system_mat;
0081   img.fwd_model.solve=      @aa_fwd_solve;
0082   J_aa_p= calc_jacobian( img ); % 2 for bug in my code
0083 
0084   img.fwd_model.jacobian=   @np_calc_jacobian;
0085   img.fwd_model.system_mat= @np_calc_system_mat;
0086   img.fwd_model.solve=      @np_fwd_solve;
0087   J_np= calc_jacobian( img );
0088 
0089   img.fwd_model.jacobian=   @perturb_jacobian;
0090   img.fwd_model.system_mat= @np_calc_system_mat;
0091   img.fwd_model.solve=      @np_fwd_solve;
0092   J_np_p= calc_jacobian( img );
0093 
0094   tol = 4e-4;
0095   unit_test_cmp('J_aa - J_aa_p', norm(J_aa - J_aa_p,'fro')/norm(J_aa,'fro'), 0, tol);
0096   unit_test_cmp('J_np - J_np_p', norm(J_np - J_np_p,'fro')/norm(J_np,'fro'), 0, tol);
0097   unit_test_cmp('J_np - J_aa_p', norm(J_np - J_aa_p,'fro')/norm(J_np,'fro'), 0, tol);
0098   unit_test_cmp('J_np - J_aa  ', norm(J_np - J_aa  ,'fro')/norm(J_np,'fro'), 0, tol);
0099 
0100 % Demo model with EIDORS3D
0101   imdl= mk_common_model('n3r2',16);
0102   img= calc_jacobian_bkgnd(imdl);
0103   for i=1:16; imdl.fwd_model.electrode(i).z_contact= .01;end
0104 
0105 % Calculate the normal Jacobian
0106   img.fwd_model.jacobian=   @np_calc_jacobian;
0107   img.fwd_model.system_mat= @np_calc_system_mat;
0108   img.fwd_model.solve=      @np_fwd_solve;
0109   J_np= calc_jacobian( img );
0110 
0111 function display_jacobian_deltas
0112 % Calculate the perturbation Jacobian
0113   vh= fwd_solve(img);
0114 
0115   for el= 1:50:size(img.fwd_model.elems,1);
0116      fprintf('\nelem#%03d: ',el);
0117      for delta= [1e-4,1e-6,1e-8] % delta is perturb
0118         img_i = img;
0119         img_i.elem_data(el)= img_i.elem_data(el) + delta;
0120         vi= fwd_solve(img_i);
0121         J_p = (vi.meas - vh.meas)/delta; % perturb Jacobian
0122         fprintf(' %8.6f', norm(J_p - J_np(:,el))/norm(J_np(:,el)));
0123      end
0124   end
0125