inv_solve_abs_CG

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

 img = inv_solve_abs_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_abs_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_abs_CG.m 4918 2015-05-07 20:57:25Z alistair_boyle $

0001 function img= inv_solve_abs_CG( inv_model, data1);
0002 %function img= inv_solve_abs_CG( inv_model, data1);
0003 % INV_SOLVE_ABS_CG
0004 % This function calls INV_SOLVE_ABS_CORE to find a Conjugate-Gradient
0005 % iterative solution.
0006 %
0007 % img = inv_solve_abs_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_abs_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_abs_CG.m 4918 2015-05-07 20:57:25Z alistair_boyle $
0027 
0028 %--------------------------
0029 % UNIT_TEST?
0030 if isstr(inv_model) && strcmp(inv_model,'UNIT_TEST'); img = do_unit_test; return; end
0031 
0032 % merge legacy options locations
0033 inv_model = deprecate_imdl_opt(inv_model, 'parameters');
0034 inv_model = deprecate_imdl_opt(inv_model, 'inv_solve');
0035 inv_model = deprecate_imdl_opt(inv_model, 'inv_solve_abs_core');
0036 
0037 % fixed working data... otherwise we wouldn't be calling this function!
0038 if ~isfield(inv_model, 'inv_solve_abs_CG')
0039    inv_model.inv_solve_abs_CG = struct;
0040 end
0041 if ~isfield(inv_model.inv_solve_abs_CG, 'update_func')
0042    inv_model.inv_solve_abs_CG.update_func = @CG_update;
0043 end
0044 if ~isfield(inv_model.inv_solve_abs_CG, 'beta_func')
0045    inv_model.inv_solve_abs_CG.beta_func = @beta_reset_polak_ribiere;
0046 end
0047 
0048 % inv_model.inv_solve_abs_CG -> inv_solve_abs_core
0049 if isfield(inv_model, 'inv_solve_abs_CG')
0050    inv_model.inv_solve_abs_core = inv_model.inv_solve_abs_CG;
0051    inv_model = rmfield(inv_model, 'inv_solve_abs_CG');
0052 end
0053 
0054 img = inv_solve_abs_core(inv_model, data1);
0055 
0056 if isfield(img, 'inv_solve_abs_core')
0057   img.inv_solve_abs_CG = img.inv_solve_abs_core;
0058   img=rmfield(img, 'inv_solve_abs_core');
0059 end
0060 
0061 function dx = CG_update(J, W, hp2RtR, dv, de)
0062 %  % the actual update
0063   dx = (J'*W*J + hp2RtR)\(J'*dv + hp2RtR*de);
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 imdl = deprecate_imdl_opt(imdl,opt)
0119    if ~isfield(imdl, opt)
0120       return;
0121    end
0122    if ~isstruct(imdl.(opt))
0123       error(['unexpected inv_model.' opt ' where ' opt ' is not a struct... i do not know what to do']);
0124    end
0125 
0126    % warn on anything but inv_model.inv_solve.calc_solution_error
0127    Af = fieldnames(imdl.(opt));
0128    if ~strcmp(opt, 'inv_solve') || (length(Af(:)) ~= 1) || ~strcmp(Af(:),'calc_solution_error')
0129       disp(imdl)
0130       disp(imdl.(opt))
0131       warning('EIDORS:deprecatedParameters',['INV_SOLVE inv_model.' opt '.* are deprecated in favor of inv_model.inv_solve_abs_CG.* as of 30-Apr-2014.']);
0132    end
0133 
0134    if ~isfield(imdl, 'inv_solve_abs_CG')
0135       imdl.inv_solve_abs_CG = imdl.(opt);
0136    else % we merge
0137       % merge struct trick from:
0138       %  http://stackoverflow.com/questions/38645
0139       for i = fieldnames(imdl.(opt))'
0140          imdl.inv_solve_abs_CG.(i{1})=imdl.(opt).(i{1});
0141       end
0142    end
0143    imdl = rmfield(imdl, opt);
0144 
0145 function pass = do_unit_test()
0146    pass=1;
0147    imdl=mk_common_model('a2c');
0148    imdl.reconst_type='absolute'; % ***
0149    imdl.solve=@inv_solve_abs_CG; % ***
0150    fimg=mk_image(imdl,1);
0151    fimg.elem_data(5:10)=1.1;
0152    vi=fwd_solve(fimg);
0153    img=inv_solve(imdl,vi); % ***
0154    clf; subplot(121); show_fem(fimg,1); title('forward model');
0155         subplot(122); show_fem(img,1);  title('reconstruction');
0156    try unit_test_cmp('fwd vs. reconst', fimg.elem_data, img.elem_data, 0.08);
0157    catch me; disp(me.message); pass=0; end
0158 %   pass = inv_solve_abs_core('UNIT_TEST', 'inv_solve_abs_CG');