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 5176 2016-02-04 16:28:29Z 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 5176 2016-02-04 16:28:29Z alistair_boyle $
0027 
0028 %--------------------------
0029 % UNIT_TEST?
0030 if ischar(inv_model) && strcmp(inv_model,'UNIT_TEST'); img = 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 = @CG_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 function dx = CG_update(J, W, hp2RtR, dv, de, opt)
0061 %  % the actual update
0062   dx = -(J'*W*J + hp2RtR)\(J'*dv + hp2RtR*de);
0063 %   dx = J'*dv;
0064 
0065 % Fletcher-Reeves
0066 % cite: http://en.wikipedia.org/wiki/Nonlinear_conjugate_gradient_method
0067 %            dx_k^T  dx_k
0068 % beta = ---------------------
0069 %         dx_{k-1}^T dx_{k-1}
0070 function beta = beta_fletcher_reeves(dx_k, dx_km1, sx_km1)
0071    beta = (dx_k' * dx_k)/(dx_km1' * dx_km1);
0072    if isinf(beta) || isnan(beta)
0073       beta = 0;
0074    end
0075 
0076 % Polak-Ribiere
0077 % cite: http://en.wikipedia.org/wiki/Nonlinear_conjugate_gradient_method
0078 %         dx_k^T ( dx_k - dx_{k-1} )
0079 % beta = ----------------------------
0080 %             dx_{k-1}^T dx_{k-1}
0081 function beta = beta_polak_ribiere(dx_k, dx_km1, sx_km1)
0082    ddx = dx_k - dx_km1;
0083    beta = (dx_k' * ddx)/(dx_km1' * dx_km1);
0084    if isinf(beta) || isnan(beta)
0085       beta = 0;
0086    end
0087 
0088 % Hestenes-Stiefel
0089 % cite: http://en.wikipedia.org/wiki/Nonlinear_conjugate_gradient_method
0090 %          dx_k^T    ( dx_k - dx_{k-1} )
0091 % beta = - -------------------------------
0092 %          s_{k-1}^T ( dx_k - dx_{k-1} )
0093 function beta = beta_hestenes_stiefel(dx_k, dx_km1, sx_km1)
0094    ddx = dx_k - dx_km1;
0095    beta = -(dx_k' * ddx)/(sx_km1' * ddx);
0096    if isinf(beta) || isnan(beta)
0097       beta = 0;
0098    end
0099 
0100 % Polak-Ribiere with reset
0101 % cite: http://en.wikipedia.org/wiki/Nonlinear_conjugate_gradient_method
0102 % beta = max(0, beta_pr)
0103 function beta = beta_reset_polak_ribiere(dx_k, dx_km1, sx_km1)
0104    beta = beta_polak_ribiere(dx_k, dx_km1, sx_km1);
0105    if beta < 0
0106       beta = 0;
0107    end
0108 
0109 % Dai-Yuan
0110 % cite: http://en.wikipedia.org/wiki/Nonlinear_conjugate_gradient_method
0111 %                dx_k^T    dx_k
0112 % beta = - -------------------------------
0113 %          s_{k-1}^T ( dx_k - dx_{k-1} )
0114 function beta = beta_dai_yuan(dx_k, dx_km1, sx_km1)
0115    beta = -dx_k'*dx_k/(sx_km1'*(dx_k-dx_km1));
0116 
0117 function pass = do_unit_test()
0118    pass=1;
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    try unit_test_cmp('fwd vs. reconst', fimg.elem_data, img.elem_data, 0.08);
0129    catch me; disp(me.message); pass=0; end
0130 %   pass = inv_solve_core('UNIT_TEST', 'inv_solve_cg');

Generated on Wed 21-Jun-2017 09:29:07 by m2html © 2005