inv_solve_cg

PURPOSE ^

function img= inv_solve_cg( inv_model, data1);

SYNOPSIS ^

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

DESCRIPTION ^

function img= inv_solve_cg( inv_model, data1);
 INV_SOLVE_CG
 This function calls INV_SOLVE_ABS_CORE to find a Conjugate-Gradient
 iterative solution.

 img = inv_solve_cg( 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_cg; % ***
  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_ABS_CORE for arguments, options and parameters.

 (C) 2014 Alistair Boyle
 License: GPL version 2 or version 3
 $Id: inv_solve_cg.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_cg( inv_model, data1, data2);
0002 %function img= inv_solve_cg( inv_model, data1);
0003 % INV_SOLVE_CG
0004 % This function calls INV_SOLVE_ABS_CORE to find a Conjugate-Gradient
0005 % iterative solution.
0006 %
0007 % img = inv_solve_cg( 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_cg; % ***
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_ABS_CORE for arguments, options and parameters.
0023 %
0024 % (C) 2014 Alistair Boyle
0025 % License: GPL version 2 or version 3
0026 % $Id: inv_solve_cg.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 % fixed working data... otherwise we wouldn't be calling this function!
0033 if ~isfield(inv_model, 'inv_solve_cg')
0034    inv_model.inv_solve_cg = struct;
0035 end
0036 if ~isfield(inv_model.inv_solve_cg, 'update_func')
0037    inv_model.inv_solve_cg.update_func = 'GN_update';
0038 end
0039 if ~isfield(inv_model.inv_solve_cg, 'beta_func')
0040    inv_model.inv_solve_cg.beta_func = @beta_reset_polak_ribiere;
0041 end
0042 
0043 % inv_model.inv_solve_cg -> inv_solve_core
0044 if isfield(inv_model, 'inv_solve_cg')
0045    inv_model.inv_solve_core = inv_model.inv_solve_cg;
0046    inv_model = rmfield(inv_model, 'inv_solve_cg');
0047 end
0048 
0049 if nargin > 2
0050    img = inv_solve_core(inv_model, data1, data2);
0051 else
0052    img = inv_solve_core(inv_model, data1);
0053 end
0054 
0055 if isfield(img, 'inv_solve_core')
0056   img.inv_solve_cg = img.inv_solve_core;
0057   img=rmfield(img, 'inv_solve_core');
0058 end
0059 
0060 % Note: we prefer to still be able to apply regularization, so use the
0061 % GN_update by default, but you can try this if you are determined. It
0062 % generally doesn't give good reconstructions for most EIT inverse problems.
0063 function dx = CG_update(J, W, hp2RtR, dv, de, opt)
0064    dx = J'*dv;
0065 
0066 % Fletcher-Reeves
0067 % cite: http://en.wikipedia.org/wiki/Nonlinear_conjugate_gradient_method
0068 %            dx_k^T  dx_k
0069 % beta = ---------------------
0070 %         dx_{k-1}^T dx_{k-1}
0071 function beta = beta_fletcher_reeves(dx_k, dx_km1, sx_km1)
0072    beta = (dx_k' * dx_k)/(dx_km1' * dx_km1);
0073    if isinf(beta) || isnan(beta)
0074       beta = 0;
0075    end
0076 
0077 % Polak-Ribiere
0078 % cite: http://en.wikipedia.org/wiki/Nonlinear_conjugate_gradient_method
0079 %         dx_k^T ( dx_k - dx_{k-1} )
0080 % beta = ----------------------------
0081 %             dx_{k-1}^T dx_{k-1}
0082 function beta = beta_polak_ribiere(dx_k, dx_km1, sx_km1)
0083    ddx = dx_k - dx_km1;
0084    beta = (dx_k' * ddx)/(dx_km1' * dx_km1);
0085    if isinf(beta) || isnan(beta)
0086       beta = 0;
0087    end
0088 
0089 % Hestenes-Stiefel
0090 % cite: http://en.wikipedia.org/wiki/Nonlinear_conjugate_gradient_method
0091 %          dx_k^T    ( dx_k - dx_{k-1} )
0092 % beta = - -------------------------------
0093 %          s_{k-1}^T ( dx_k - dx_{k-1} )
0094 function beta = beta_hestenes_stiefel(dx_k, dx_km1, sx_km1)
0095    ddx = dx_k - dx_km1;
0096    beta = -(dx_k' * ddx)/(sx_km1' * ddx);
0097    if isinf(beta) || isnan(beta)
0098       beta = 0;
0099    end
0100 
0101 % Polak-Ribiere with reset
0102 % cite: http://en.wikipedia.org/wiki/Nonlinear_conjugate_gradient_method
0103 % beta = max(0, beta_pr)
0104 function beta = beta_reset_polak_ribiere(dx_k, dx_km1, sx_km1)
0105    beta = beta_polak_ribiere(dx_k, dx_km1, sx_km1);
0106    if beta < 0
0107       beta = 0;
0108    end
0109 
0110 % Dai-Yuan
0111 % cite: http://en.wikipedia.org/wiki/Nonlinear_conjugate_gradient_method
0112 %                dx_k^T    dx_k
0113 % beta = - -------------------------------
0114 %          s_{k-1}^T ( dx_k - dx_{k-1} )
0115 function beta = beta_dai_yuan(dx_k, dx_km1, sx_km1)
0116    beta = -dx_k'*dx_k/(sx_km1'*(dx_k-dx_km1));
0117 
0118 function do_unit_test()
0119    imdl=mk_common_model('a2c');
0120    imdl.reconst_type='absolute'; % ***
0121    imdl.solve=@inv_solve_cg; % ***
0122    fimg=mk_image(imdl,1);
0123    fimg.elem_data(5:10)=1.1;
0124    vi=fwd_solve(fimg);
0125    img=inv_solve(imdl,vi); % ***
0126    clf; subplot(121); show_fem(fimg,1); title('forward model');
0127         subplot(122); show_fem(img,1);  title('reconstruction');
0128    unit_test_cmp('fwd vs. reconst', fimg.elem_data, img.elem_data, 0.08);
0129 %  inv_solve_core('UNIT_TEST', 'inv_solve_cg');

Generated on Fri 30-Dec-2022 19:44:54 by m2html © 2005