jacobian_perturb

 JACOBIAN_PERTURB: J= jacobian_perturb( 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.jacobian_perturb.delta   - delta perturbation to use (default 1e-6)
 img = image background for jacobian calc

0001 function J= jacobian_perturb( fwd_model, img)
0002 % JACOBIAN_PERTURB: J= jacobian_perturb( 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.jacobian_perturb.delta   - delta perturbation to use (default 1e-6)
0010 % img = image background for jacobian calc
0011 
0012 % (C) 2006 Andy Adler. License: GPL version 2 or version 3
0013 % $Id: jacobian_perturb.m 6721 2024-04-02 18:37:40Z aadler $
0014 
0015 if ischar(fwd_model) && strcmp(fwd_model,'UNIT_TEST'); do_unit_test; return;end
0016 
0017 if isfield(fwd_model,'jacobian_perturb')
0018    delta = fwd_model.jacobian_perturb.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 jnk = data_mapper(img);
0028 curprms = jnk.current_params;
0029 
0030 if strcmp(curprms, 'conductivity') && ~isfield(img,'conductivity')
0031    curprms = '';
0032 end
0033 % solve one time to get the size
0034 d0= fwd_solve( img );
0035 
0036 if isfield(img.fwd_model,'coarse2fine');
0037    Jcol= perturb_c2f(img, 1, delta, d0, curprms);
0038    Jrows= size(img.fwd_model.coarse2fine,2);
0039    J= zeros(length(Jcol), Jrows );
0040    J(:,1)= Jcol;
0041    for i=2:size(img.fwd_model.coarse2fine,2);
0042      J(:,i)= perturb_c2f(img, i, delta, d0, curprms);
0043      if rem(i,50)==0; fprintf('+'); end
0044    end
0045 else
0046    Jcol= perturb(img, 1, delta, d0, curprms);
0047 
0048    J= zeros(length(Jcol), n_elem);
0049    J(:,1)= Jcol;
0050 
0051    for i=2:n_elem
0052      J(:,i)= perturb(img, i, delta, d0, curprms);
0053    end
0054 end
0055 
0056 
0057 function Jcol= perturb( img, i, delta, d0, curprms)
0058    if ~isempty(curprms)
0059       img.(curprms).elem_data(i)= img.(curprms).elem_data(i) + delta;
0060    else
0061       img.elem_data(i)= img.elem_data(i) + delta;
0062    end
0063    di= fwd_solve( img );
0064    Jcol = (1/delta) * (di.meas - d0.meas);
0065 
0066 function Jcol= perturb_c2f( img, i, delta, d0, curprms)
0067    if ~isempty(curprms)
0068       img.(curprms).elem_data= img.(curprms).elem_data + delta*img.fwd_model.coarse2fine(:,i);
0069    else
0070        img.elem_data(i)= img.elem_data(i) + delta; %*img.fwd_model.coarse2fine(:,i);
0071 %       img.elem_data= img.fwd_model.coarse2fine*img.elem_data + delta*img.fwd_model.coarse2fine(:,i);
0072    end
0073    di= fwd_solve( img );
0074    Jcol = (1/delta) * (di.meas - d0.meas);
0075 
0076 
0077 function do_unit_test
0078 % Perturbation Jacobians
0079 % $Id: jacobian_perturb.m 6721 2024-04-02 18:37:40Z aadler $
0080 
0081   models ={'c2c2','a3cr','n3r2'};
0082   for delta = [1e-7,1e-6,1e-5,1e-4];
0083      fprintf('Using delta = %6.3g\n',delta);
0084   for i = 1:length(models);
0085      Nel = 16; if i==3; Nel=[16,2]; end
0086      imdl= mk_common_model(models{i},Nel);
0087      imdl.fwd_model.nodes = imdl.fwd_model.nodes*.25;
0088      img= calc_jacobian_bkgnd(imdl);
0089 
0090      img.fwd_model = mdl_normalize(img.fwd_model, 0);
0091 
0092      img.fwd_model.jacobian=   @jacobian_adjoint;
0093      img.fwd_model.system_mat= @system_mat_1st_order;
0094      img.fwd_model.solve=      @fwd_solve_1st_order;
0095      J_aa= calc_jacobian( img );
0096 %    display_jacobian_deltas(img, J_aa) % show perturbation size
0097 
0098      img.fwd_model.jacobian=   @jacobian_perturb;
0099      img.fwd_model.jacobian_perturb.delta = delta;
0100      J_aa_p= calc_jacobian( img );
0101 
0102      unit_test_cmp(['J_ - J_p:', img.fwd_model.name] , J_aa, J_aa_p, 1e-4);
0103    end
0104    end
0105 
0106 
0107 function display_jacobian_deltas(img, J);
0108 % Calculate the perturbation Jacobian
0109   vh= fwd_solve(img);
0110 
0111   for el= 1:50:size(img.fwd_model.elems,1);
0112      fprintf('\nelem#%03d: ',el);
0113      for delta= [1e-4,1e-6,1e-8] % delta is perturb
0114         img_i = img;
0115         img_i.elem_data(el)= img_i.elem_data(el) + delta;
0116         vi= fwd_solve(img_i);
0117         J_p = (vi.meas - vh.meas)/delta; % perturb Jacobian
0118         fprintf(' %8.6f', norm(J_p - J(:,el))/norm(J(:,el)));
0119      end
0120   end
0121