inv_solve_gn

function img= inv_solve_gn( inv_model, data1);
 INV_SOLVE_GN
 This function calls INV_SOLVE_CORE to find a Gauss-Newton
 iterative solution.

 img = inv_solve_gn( inv_model, data1 )
   img        => output image data (or vector of images)
   inv_model  => inverse model struct
   data1      => measurements

 Example:
  imdl=mk_common_model('a2c');
  imdl.reconst_type='absolute'; % ***
  imdl.solve=@inv_solve_gn; % ***
  fimg=mk_image(imdl,1);
  fimg.elem_data(5:10)=1.1;
  vi=fwd_solve(fimg);
  img=inv_solve(imdl,vi); % ***
  show_fem(img,1);
 Example#2:
imdl = mk_common_model('b2c2',16);
imdl.jacobian_bkgnd.value = 2;
img=mk_image(imdl);    vh=fwd_solve(img);
img.elem_data([5])=5;   vi=fwd_solve(img);
imgr = inv_solve(imdl,vh,vi);  % regular solve

imdl.solve = @inv_solve_gn;
imdl.inv_solve_gn.max_iterations = 5;
imdl.inv_solve_gn.min_value = 2;
imgr = inv_solve(imdl,vh,vi); % new solve

 See INV_SOLVE_CORE for arguments, options and parameters.

 (C) 2014 Alistair Boyle
 License: GPL version 2 or version 3
 $Id: inv_solve_gn.m 5917 2019-03-05 06:43:37Z aadler $

0001 function img= inv_solve_gn( inv_model, data1, data2);
0002 %function img= inv_solve_gn( inv_model, data1);
0003 % INV_SOLVE_GN
0004 % This function calls INV_SOLVE_CORE to find a Gauss-Newton
0005 % iterative solution.
0006 %
0007 % img = inv_solve_gn( inv_model, data1 )
0008 %   img        => output image data (or vector of images)
0009 %   inv_model  => inverse model struct
0010 %   data1      => measurements
0011 %
0012 % Example:
0013 %  imdl=mk_common_model('a2c');
0014 %  imdl.reconst_type='absolute'; % ***
0015 %  imdl.solve=@inv_solve_gn; % ***
0016 %  fimg=mk_image(imdl,1);
0017 %  fimg.elem_data(5:10)=1.1;
0018 %  vi=fwd_solve(fimg);
0019 %  img=inv_solve(imdl,vi); % ***
0020 %  show_fem(img,1);
0021 % Example#2:
0022 %imdl = mk_common_model('b2c2',16);
0023 %imdl.jacobian_bkgnd.value = 2;
0024 %img=mk_image(imdl);    vh=fwd_solve(img);
0025 %img.elem_data([5])=5;   vi=fwd_solve(img);
0026 %imgr = inv_solve(imdl,vh,vi);  % regular solve
0027 %
0028 %imdl.solve = @inv_solve_gn;
0029 %imdl.inv_solve_gn.max_iterations = 5;
0030 %imdl.inv_solve_gn.min_value = 2;
0031 %imgr = inv_solve(imdl,vh,vi); % new solve
0032 %
0033 % See INV_SOLVE_CORE for arguments, options and parameters.
0034 %
0035 % (C) 2014 Alistair Boyle
0036 % License: GPL version 2 or version 3
0037 % $Id: inv_solve_gn.m 5917 2019-03-05 06:43:37Z aadler $
0038 
0039 %--------------------------
0040 % UNIT_TEST?
0041 if ischar(inv_model) && strcmp(inv_model,'UNIT_TEST'); do_unit_test; return; end
0042 
0043 
0044 % inv_model.inv_solve_gn -> inv_solve_core
0045 if isfield(inv_model, 'inv_solve_gn')
0046    if isfield(inv_model, 'inv_solve_core')
0047       error('inv_model.inv_solve_gn replaces inv_model.inv_solve_core, parameters will be lost');
0048    end
0049    inv_model.inv_solve_core = inv_model.inv_solve_gn;
0050    inv_model = rmfield(inv_model, 'inv_solve_gn');
0051 end
0052 
0053 if nargin > 2
0054   img = inv_solve_core(inv_model, data1, data2);
0055 else
0056   img = inv_solve_core(inv_model, data1);
0057 end
0058 
0059 if isfield(img, 'inv_solve_core')
0060   img.inv_solve_gn = img.inv_solve_core;
0061   img=rmfield(img, 'inv_solve_core');
0062 end
0063 
0064 function do_unit_test()
0065    imdl=mk_common_model('a2c');
0066    imdl.reconst_type='absolute'; % ***
0067    imdl.solve=@inv_solve_gn; % ***
0068    fimg=mk_image(imdl,1);
0069    fimg.elem_data(5:10)=1.1;
0070    vi=fwd_solve(fimg);
0071    img=inv_solve(imdl,vi); % ***
0072    clf; subplot(121); show_fem(fimg,1); title('forward model');
0073         subplot(122); show_fem(img,1);  title('reconstruction');
0074    unit_test_cmp('fwd vs. reconst', fimg.elem_data, img.elem_data, 0.08);
0075 
0076    do_unit_test_abs_diff();
0077 
0078 function do_unit_test_abs_diff()
0079    imdl = mk_geophysics_model('h2c',32);
0080    imdl.solve = 'inv_solve_gn';
0081    imdl.RtR_prior = 'prior_tikhonov';
0082    imdl.inv_solve_gn.verbose = 0;
0083    imdl.inv_solve_gn.elem_working = 'log_conductivity';
0084 %   imdl.inv_solve_gn.meas_working = 'apparent_resistivity';
0085    % build fwd simulations
0086    ctrs = interp_mesh(imdl.fwd_model);
0087    x = ctrs(:,1);
0088    y = ctrs(:,2);
0089    r1 = sqrt((x-20).^2  + (y+25).^2);
0090    r2 = sqrt((x-125).^2 + (y+45).^2);
0091    r3 = sqrt((x-50).^2  + (y+35).^2);
0092    imga = mk_image(imdl,1);
0093    imga.elem_data(r1<15)= 0.5;
0094    imga.elem_data(r2<20)= 2;
0095    va = fwd_solve(imga);
0096    imgb = imga;
0097    imgb.elem_data(r3<15)= 1.3;
0098    vb = fwd_solve(imgb);
0099    clf; subplot(121); show_fem(imga,1); subplot(122); show_fem(imgb,1); drawnow;
0100    % reconstruct meas va as absolute
0101    disp('*** INVERSE CRIME ALERT **** INVERSE CRIME ALERT ***');
0102    disp(' 1. This reconstruction has no measurement noise.');
0103    disp(' 2. This reconstruction is being performed on the same mesh the');
0104    disp('    measurements were simulated from: no discretization errors have');
0105    disp('    been introduced.');
0106    imdl.reconst_type = 'absolute';
0107    imgaa = inv_solve(imdl,va);
0108    % set background as abs sol'n, then diff solve
0109    imdl.reconst_type = 'difference';
0110    imdl.inv_solve_gn.prior_data = imgaa.elem_data;
0111    imgab = inv_solve(imdl,va,vb);
0112    clf; subplot(221); show_fem(imga,1); title('A'); subplot(222); show_fem(imgb,1); title('B');
0113         subplot(223); show_fem(imgaa,1); title('rec A'); subplot(224); show_fem(imgab,1); title('rec B-A');
0114    unit_test_cmp('abs ', imga.elem_data, imgaa.elem_data, max(imga.elem_data)/2);
0115    imgba = imga; imgba.elem_data = imgb.elem_data - imga.elem_data;
0116    imgba.elem_data = imgba.elem_data/max(imgba.elem_data);
0117    imgab.elem_data = imgab.elem_data/max(imgab.elem_data);
0118    unit_test_cmp('diff', norm(imgba.elem_data - imgab.elem_data), 0, 20);