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