function w = conjgrad_jvs(w0, fun, grd, data, nprnt, nfun, tau, Smax)

% function w = conjgrad(w0, fun, grd, data, [nprnt = -1, nfun = 100*length(w), tau = 1e-10])
%  % ALL VECS ARE COL VECS.
% Implementation of P. Williams conjugate gradient minimisation
% with adaptive quadratic (inexact) line-search
%
% w0    = initial value (column vector)
% fun   = cost function name (string) (returns scalar)
% grd   = cost gradient name (string) (returns column vector)
% data  = extra data for fun and grd
% nprnt = print info every nprnt evaluations (default no output )
% nfun  = max number of function evaluations (termination if n > nfun 
% tau   = gradient tolerance                  or |g| < tau)
%
% w     = parameter vector at minimum

global num_linesearches;

% ANNOUNCE
fprintf('THIS IS CONJGRAD ...\n');

% USE PORRILL's MODIFICATIONS?
PORRILL = 1; % See rostest.m for evidence.
% PORRILL causes lambda to explode with ICA.

w = w0;
if nargin < 5, nprnt = -1; end
if nargin < 6, nfun = 100*length(w); end
if nargin < 7, tau  = 1e-6; end

% epsilon = 1e-3;
epsilon = sqrt(eps); % eps is smallest machine increment to 1
pie = 0.05;

lambda = 1; lambdamin = eps; lambdamax = 1/eps;

f = feval(fun, w, data);
g = feval(grd, w, data);
s = -g;
S = 0; n = 0;
success = 1;
Smax = 2*length(w);
jnostep=0;

% KEEP RECORD OF FUNCTION VALS.
fs=zeros(1,nfun);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
while (n < nfun & norm(g)/norm(w) > tau)
    n = n+1;
    num_linesearches = n;
	
    %%%%%%%%%%%%%%%%%%%%%
    % STEP 1.
    %%%%%%%%%%%%%%%%%%%%%
    if success == 1
        mu = s'*g;
        if mu >= 0
            s = -g;
            mu = s'*g;
            S = 0;
        end
        kappa = s'*s;
		% jsize(s,'s'); % 1x24
        sigma = epsilon/sqrt(kappa);
        gs = feval(grd, w+sigma*s, data);
        gamma = s'*(gs-g)/sigma;
    end
    
    %%%%%%%%%%%%%%%%%%%%%
    % STEP 2.
    %%%%%%%%%%%%%%%%%%%%%
    delta = gamma+lambda*kappa;

    %%%%%%%%%%%%%%%%%%%%%
    % STEP 3.
    %%%%%%%%%%%%%%%%%%%%%
    if delta <= 0
        delta = lambda*kappa;
        lambda = lambda-gamma/kappa;
    end

    %%%%%%%%%%%%%%%%%%%%%
    % STEP 4.
    %%%%%%%%%%%%%%%%%%%%%
    alpha = -mu/delta;
    epsilon = epsilon*(alpha/sigma)^pie;
    
    % f = feval(fun, w, data);
    fa = feval(fun, w+alpha*s, data);

	step_len = norm(alpha*s);

	%plot_cross_section(fun,w,data,s,alpha);

    %%%%%%%%%%%%%%%%%%%%%
    % STEP 5. rho is the ratio of actual/estimated reduction in f.
    %%%%%%%%%%%%%%%%%%%%%
    if PORRILL
        rho = (fa-f)/(alpha*mu+0.5*alpha^2*gamma);
        success = (fa < f+eps);
    else
    	rho = 2*(fa-f)/(alpha*mu); % original
        success = (rho >= 0);      % original
    end;

    %%%%%%%%%%%%%%%%%%%%%
    % STEP 6.
    %%%%%%%%%%%%%%%%%%%%%
    if rho < 0.25
        lambda = min(lambda+delta*(1-rho)/kappa, lambdamax);
    elseif rho > 0.75
        lambda = max(lambda/2, lambdamin);
    end
    
    %%%%%%%%%%%%%%%%%%%%%
    % STEP 7.
    %%%%%%%%%%%%%%%%%%%%%
    if success == 1
        w = w+alpha*s;
        gold = g;
        g = feval(grd, w, data);
		f = fa;
        S = S+1;
    end
    
    %%%%%%%%%%%%%%%%%%%%%
    % STEP 8.
    %%%%%%%%%%%%%%%%%%%%%
    if S == Smax
        s = -g;
        success = 1;
        S = 0;
    elseif success == 1
        beta = (gold-g)'*g/mu;
        s = -g+beta*s;
    end

    if nprnt > 0 & (nprnt == 1 | rem(n, nprnt) == 1)
    fprintf('n=%6d |rho=%12.6f lambda=%12.6f | f=%12.6g | |g|=%12.6g step_len= %12.6g\n', ...
			n,	  rho,		 lambda,		  f,		norm(g)/norm(w),step_len);
	fs(n)=f;
	plot_cross_section(fun,w,data,s,alpha);
	% jfig(20); clf; plot(fs(2:n)); title('Function values');
    end

%fprintf('success=%d jnostep=%d f=%6.12f\n',success,jnostep,f);

    % STOP ?
    if success==0 jnostep=jnostep+1;
    else jnostep=0; end;
    if jnostep>9 break; end;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
