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
 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
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 5559 2017-06-18 17:50:49Z 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 5559 2017-06-18 17:50:49Z aadler $
0080 
0081   models ={'c2c2','a3cr','n3r2'};
0082   for i = 1:length(models)
0083      imdl= mk_common_model(models{i},16);
0084      imdl.fwd_model.nodes = imdl.fwd_model.nodes*.25;
0085      img= calc_jacobian_bkgnd(imdl);
0086 
0087      img.fwd_model = mdl_normalize(img.fwd_model, 0);
0088 
0089      img.fwd_model.jacobian=   @jacobian_adjoint;
0090      img.fwd_model.system_mat= @system_mat_1st_order;
0091      img.fwd_model.solve=      @fwd_solve_1st_order;
0092      J_aa= calc_jacobian( img ); % 2 for bug in my code
0093 
0094      img.fwd_model.jacobian=   @jacobian_perturb;
0095      J_aa_p= calc_jacobian( img ); % 2 for bug in my code
0096 
0097      unit_test_cmp('J_ - J_p', norm(J_aa - J_aa_p,'fro')/norm(J_aa,'fro'), 0, 1e-3);
0098    end
0099 
0100 
0101 function display_jacobian_deltas
0102 % Calculate the perturbation Jacobian
0103   vh= fwd_solve(img);
0104 
0105   for el= 1:50:size(img.fwd_model.elems,1);
0106      fprintf('\nelem#%03d: ',el);
0107      for delta= [1e-4,1e-6,1e-8] % delta is perturb
0108         img_i = img;
0109         img_i.elem_data(el)= img_i.elem_data(el) + delta;
0110         vi= fwd_solve(img_i);
0111         J_p = (vi.meas - vh.meas)/delta; % perturb Jacobian
0112         fprintf(' %8.6f', norm(J_p - J_np(:,el))/norm(J_np(:,el)));
0113      end
0114   end
0115