Home > eidors > solvers > inverse > inv_solve_gn.m

inv_solve_gn

PURPOSE ^

function img= inv_solve_gn( inv_model, data1);

SYNOPSIS ^

function img= inv_solve_gn( inv_model, data1, data2);

DESCRIPTION ^

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);

 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 5613 2017-07-07 18:35:01Z alistair_boyle $

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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 %
0022 % See INV_SOLVE_CORE for arguments, options and parameters.
0023 %
0024 % (C) 2014 Alistair Boyle
0025 % License: GPL version 2 or version 3
0026 % $Id: inv_solve_gn.m 5613 2017-07-07 18:35:01Z alistair_boyle $
0027 
0028 %--------------------------
0029 % UNIT_TEST?
0030 if ischar(inv_model) && strcmp(inv_model,'UNIT_TEST'); do_unit_test; return; end
0031 
0032 
0033 % inv_model.inv_solve_gn -> inv_solve_core
0034 if isfield(inv_model, 'inv_solve_gn')
0035    if isfield(inv_model, 'inv_solve_core')
0036       error('inv_model.inv_solve_gn replaces inv_model.inv_solve_core, parameters will be lost');
0037    end
0038    inv_model.inv_solve_core = inv_model.inv_solve_gn;
0039    inv_model = rmfield(inv_model, 'inv_solve_gn');
0040 end
0041 
0042 if nargin > 2
0043   img = inv_solve_core(inv_model, data1, data2);
0044 else
0045   img = inv_solve_core(inv_model, data1);
0046 end
0047 
0048 if isfield(img, 'inv_solve_core')
0049   img.inv_solve_gn = img.inv_solve_core;
0050   img=rmfield(img, 'inv_solve_core');
0051 end
0052 
0053 function do_unit_test()
0054    imdl=mk_common_model('a2c');
0055    imdl.reconst_type='absolute'; % ***
0056    imdl.solve=@inv_solve_gn; % ***
0057    fimg=mk_image(imdl,1);
0058    fimg.elem_data(5:10)=1.1;
0059    vi=fwd_solve(fimg);
0060    img=inv_solve(imdl,vi); % ***
0061    clf; subplot(121); show_fem(fimg,1); title('forward model');
0062         subplot(122); show_fem(img,1);  title('reconstruction');
0063    unit_test_cmp('fwd vs. reconst', fimg.elem_data, img.elem_data, 0.08);
0064 
0065    do_unit_test_abs_diff();
0066 
0067 function do_unit_test_abs_diff()
0068    imdl = mk_geophysics_model('h2c',32);
0069    imdl.solve = 'inv_solve_gn';
0070    imdl.RtR_prior = 'prior_tikhonov';
0071    imdl.inv_solve_gn.verbose = 0;
0072    imdl.inv_solve_gn.elem_working = 'log_conductivity';
0073 %   imdl.inv_solve_gn.meas_working = 'apparent_resistivity';
0074    % build fwd simulations
0075    ctrs = interp_mesh(imdl.fwd_model);
0076    x = ctrs(:,1);
0077    y = ctrs(:,2);
0078    r1 = sqrt((x-20).^2  + (y+25).^2);
0079    r2 = sqrt((x-125).^2 + (y+45).^2);
0080    r3 = sqrt((x-50).^2  + (y+35).^2);
0081    imga = mk_image(imdl,1);
0082    imga.elem_data(r1<15)= 0.5;
0083    imga.elem_data(r2<20)= 2;
0084    va = fwd_solve(imga);
0085    imgb = imga;
0086    imgb.elem_data(r3<15)= 1.3;
0087    vb = fwd_solve(imgb);
0088    clf; subplot(121); show_fem(imga,1); subplot(122); show_fem(imgb,1); drawnow;
0089    % reconstruct meas va as absolute
0090    disp('*** INVERSE CRIME ALERT **** INVERSE CRIME ALERT ***');
0091    disp(' 1. This reconstruction has no measurement noise.');
0092    disp(' 2. This reconstruction is being performed on the same mesh the');
0093    disp('    measurements were simulated from: no discretization errors have');
0094    disp('    been introduced.');
0095    imdl.reconst_type = 'absolute';
0096    imgaa = inv_solve(imdl,va);
0097    % set background as abs sol'n, then diff solve
0098    imdl.reconst_type = 'difference';
0099    imdl.inv_solve_gn.prior_data = imgaa.elem_data;
0100    imgab = inv_solve(imdl,va,vb);
0101    clf; subplot(221); show_fem(imga,1); title('A'); subplot(222); show_fem(imgb,1); title('B');
0102         subplot(223); show_fem(imgaa,1); title('rec A'); subplot(224); show_fem(imgab,1); title('rec B-A');
0103    unit_test_cmp('abs ', imga.elem_data, imgaa.elem_data, max(imga.elem_data)/2);
0104    imgba = imga; imgba.elem_data = imgb.elem_data - imga.elem_data;
0105    imgba.elem_data = imgba.elem_data/max(imgba.elem_data);
0106    imgab.elem_data = imgab.elem_data/max(imgab.elem_data);
0107    unit_test_cmp('diff', norm(imgba.elem_data - imgab.elem_data), 0, 20);

Generated on Fri 01-Jun-2018 15:59:55 by m2html © 2005