function ws = run_scg_fast( ws )

% function ws = run_scg_fast( ws )
%
% Implementation of the Scaled Conjugate Gradient (SCG) algorithm.
% Some redundancies present in "run_scg" have been eliminated.
%
% "ws" stands for "workspace".
%
% The "ws" structure is used for storing all input, transient and output variables.
% The stopping criterion is:
%    max_k( | ws.r_k | ) < ws.cv_thr
%
% For example, "w_kp1" means in latex notation "w_{k + 1}".
%
% "ws.w_1" is the user-defined initiatlization.
% "ws.w_kp1" is the result.
%
% "ws.dE_E_fun" is a pointer to function (e.g. dE_E_fun = @my_function)
% that can return both the cost "E" and its derivative "dE":
% Typical uses: 
%      dE     = feval( ws.dE_E_fun, w )
%     [dE, E] = feval( ws.dE_E_fun, w )
%
% You may want to optimize "dE_E_fun", e.g. when only "dE" is required.
%
%
% For other variables, refer to the article: notations are consistent
% with those of the NN 1993, vol. 6 article from M. Moller.
%
%
% By Guillaume Lathoud, 2005
% lathoud@idiap.ch
  
  check_param( { 'w_1', 'dE_E_fun' }, fieldnames( ws ) );

  % Default parameters

  % According to the article, they don't need much tuning,
  % as long as  0 < sigma <= 1e-4   and    0 < lambda_1 <= 1e-6
  % 
  % Lower a bit "sigma" and/or "lambda_1" if SCG does not behave appropriately

  % Step "1." in the article

  ws_default.sigma    = 1e-5;
  ws_default.lambda_1 = 1e-7;

  % By default, no minimum nb of iter
  ws_default.min_iter = -Inf;
  
  % By default we never stop before convergence
  ws_default.max_iter = +Inf;
  ws_default.max_iter_nosuccess = +Inf;

  % If verbose = 1,
  % then at each iteration "ws" is dumped, and a few messages may be displayed.
  
  ws_default.verbose  = 0;

  % In case you want some special action on "ws" at the beginning of each iteration
  ws_default.change_fun = [];
  ws_default.change_fun_par = {};

  % In case you want to use your own convergence test,
  % ( default is infinite norm of ws.r_k < ws.cv_thr
  ws_default.cv_test_fun = [];
  ws_default.cv_thr      = 1e-10;

  ws = fill_default( ws, ws_default );
  
  if ws.verbose
    disp( ' ' );
    disp( [ mfilename ' starting value of "ws":' ] );
    disp( ws );
    disp( [ mfilename ' starting value of "ws.w_1":' ] );
    disp( ws.w_1.' );
  end

  % Dimensionality
  ws.N = numel( ws.w_1 );
  
  % Do the SCG descent
  
  ws.k = 1;
  
  ws.w_k          = ws.w_1;
  ws.lambda_k     = ws.lambda_1;
  
  ws.lambda_bar_k = 0;
  
  [ ws.dE_w_k,  ws.E_w_k ] = feval( ws.dE_E_fun, ws.w_k );
  ws.p_k = -ws.dE_w_k;
  ws.r_k = -ws.dE_w_k;

  ws.success = 1;
  ws.iter_success   = 0;
  ws.iter_nosuccess = 0;
  ws.iter_scale_increase = 0;
  
  while 1

    % Step "2." in the article
    if ws.success
    
      ws.norm_p_k_square = sum( ws.p_k .^ 2 );
      ws.norm_p_k        = sqrt( ws.norm_p_k_square );
      ws.sigma_k = ws.sigma / ws.norm_p_k;
      
      dE_sigma = feval( ws.dE_E_fun, ws.w_k + ws.sigma_k * ws.p_k );
      
      ws.s_k = ( dE_sigma - ws.dE_w_k ) / ws.sigma_k;
      
      ws.delta_k = ws.p_k.' * ws.s_k;

    end

    % Step "3." in the article
    ws.delta_k = ws.delta_k + ( ws.lambda_k - ws.lambda_bar_k ) * ws.norm_p_k_square;

    % Step "4." in the article
    if ws.delta_k <= 0
      ws.lambda_bar_k = 2 * ( ws.lambda_k - ws.delta_k / ws.norm_p_k_square );
      ws.delta_k      = - ws.delta_k + ws.lambda_k * ws.norm_p_k_square;
      ws.lambda_k     = ws.lambda_bar_k;
    end
    
    % Step "5." in the article
    ws.mu_k    = ws.p_k.' * ws.r_k;
    ws.alpha_k = ws.mu_k / ws.delta_k;

    % Step "6." in the article

    % Tentative w_kp1 (will be kept as "w_kp1" only if success in step 7.)

    tmp_w_kp1 = ws.w_k + ws.alpha_k * ws.p_k;

    [ dE_alpha, E_alpha ] = feval( ws.dE_E_fun, tmp_w_kp1 );
    
    ws.DELTA_k = 2 * ws.delta_k * ( ws.E_w_k - E_alpha ) / ( ws.mu_k ^ 2 );


    % Step "7." in the article
    if ws.DELTA_k >= 0
      
      % Success

      if ws.verbose
	disp( [ mfilename ': success' ] );
      end

      % DEBUG
      disp( sprintf( 'k:%d  E: %.10g', ws.k+1, E_alpha ) );
      
      ws.w_kp1     = tmp_w_kp1;
      ws.E_w_kp1   = E_alpha;
      ws.dE_w_kp1  = dE_alpha;

      ws.r_kp1        = -ws.dE_w_kp1;
      ws.lambda_bar_k = 0;
      
      ws.success = 1;
      ws.iter_success   = ws.iter_success + 1;
      ws.iter_nosuccess = 0;
      
      % Restart or continue
      
      if mod( ws.k, ws.N ) == 0
	
	if ws.verbose
	  disp( [ mfilename ': restarting the algorithm.' ] );
	end
	
	ws.p_kp1 = ws.r_kp1;

      else

	ws.beta_k = ws.r_kp1.' * ( ws.r_kp1 - ws.r_k ) / ws.mu_k;
	ws.p_kp1  = ws.r_kp1 + ws.beta_k * ws.p_k;
	
      end

      % Reduce the scale parameter if necessary

      if ws.DELTA_k >= 0.75
	
	if ws.verbose
	  disp( [ mfilename ': reducing the scale parameter.' ] );
	end

	ws.lambda_k = ws.lambda_k / 4;
      end

    else

      % Failure
      
      if ws.verbose
	disp( [ mfilename ': failure' ] );
      end

      ws.lambda_bar_k = ws.lambda_k;

      ws.success = 0;
      ws.iter_success   = 0;
      ws.iter_nosuccess = ws.iter_nosuccess + 1;

    end

    % Step "8." in the algorithm
    % Increase the scale parameter if necessary

    if ws.DELTA_k < 0.25
      
      if ws.verbose
	disp( [ mfilename ': increasing the scale parameter.' ] );
      end

      ws.lambda_k = ws.lambda_k + ws.delta_k * ( 1 - ws.DELTA_k ) / ws.norm_p_k_square;
      
      ws.iter_scale_increase = ws.iter_scale_increase + 1;
      
      if 0

	% This is a hack, let's not include it
	% In case this is getting boring
	if ws.iter_scale_increase > 5
	  
	  ws.lambda_k = ws.lambda_k * ws.iter_scale_increase;
	  
	end

      end

    else

      ws.iter_scale_increase = 0;

    end
    
    if ws.verbose > 1
      disp( [ mfilename ': ws after step 8. :' ] );
      disp( ws );
    end
    
    % Step "9." in the algorithm
    % Termination test
    
    if isempty( ws.cv_test_fun )

      % Default convergence test (on norm_Inf)
      ws.has_converged = ( max( abs( ws.r_k ) ) < ws.cv_thr );
      
      if 0
	% This is a hack, let's not include it
	% Also if for some reason we are stuck in a weird position
	if isfield( ws, 'w_kp1' )
	  ws.has_converged = ws.has_converged | ( ( ws.iter_success > 1 ) & ( max( abs( ws.w_k - ws.w_kp1 ) ) < ws.cv_thr ) );
	end
      end
      
    else
      
      % User-defined convergence test
      ws.has_converged = feval( ws.cv_test_fun, ws );
      
    end
    
    do_continue = ( ~ws.has_converged )  &  ( ws.k < ws.max_iter )  &  ( ws.iter_nosuccess < ws.max_iter_nosuccess );
    
    do_continue = do_continue | ( ws.k < ws.min_iter );
    
    if do_continue

      % Not finished -> prepare for the next iteration

      ws.k = ws.k + 1;

      if ws.success
	
	% Update
	ws.w_k    = ws.w_kp1;
	ws.E_w_k  = ws.E_w_kp1;
	ws.dE_w_k = ws.dE_w_kp1;
	ws.r_k    = ws.r_kp1;
	ws.p_k    = ws.p_kp1;
	
	% Optional user-defined action on "w"
	if ~isempty( ws.change_fun )
	  
	  ws = feval( ws.change_fun, ws, ws.change_fun_par{:} );
	
	end
        
      end
      
    else
      
      % Finished
      if ws.verbose
	
	if ws.has_converged
	  disp( [ mfilename ': finished converging.' ] );
	else
	  disp( [ mfilename ': stopped before reaching convergence (max iter reached).' ] );
	end
	
      end
      
      break;
      
    end
    
  end % of main loop (iterations)
  