function prob2_6(method, x0, C, epsilon)
%   ACM 113, Problem 2.6
%
%       prob2_6(method, x0, C, epsilon)
%
%   the default value of C is 0.5; of epsilon is 1e-8
%
%   See also: LINESCH, OBJECTIVE

%   Ling Li, Caltech
%   04/23/2001

if nargin < 4,
    epsilon = 1e-8;
end
if nargin < 3,
    C = 0.5;
end

x = x0; f = objective(x, 0); ite = 0;
df = objective(x, 1);
if strcmp(method, 'newton'),
    % Newton
    disp('Newton method');
    while norm(df) > epsilon * (1+abs(f)),
        d2f = objective(x, 2);
        p = - inv(d2f) * df;
        [ak, x, f] = linesch(x, p, 1, C, f, df);
        disp(['  alpha_' int2str(ite) ' = ' num2str(ak)]);
        ite = ite + 1;
        df = objective(x, 1);
    end

elseif strcmp(method, 'conjugate'),
    % Conjugate
    disp('Conjugate gradient');
    df2 = df' * df; p = -df; aguess = 1;
    while norm(df) > epsilon * (1+abs(f)),
        [ak, x, f_new, it] = linesch(x, p, aguess, C, f, df);
        disp(['guess = ' num2str(aguess) ...
              ', alpha_' int2str(ite) ' = ' num2str(ak)]);
        ite = ite + 1; % + it

        df_new = objective(x, 1);
        df2_new = df_new' * df_new;
        % Polak-Ribiore method
        %b = (df2_new - df_new' * df) / df2;
        if (df_new' * df) > 0.1 * sqrt(df2_new * df2),
            % restart
            p_new = -df_new;
            aguess = 1;
        else
            b = df2_new / df2;
            p_new = b * p - df_new;
            aguess = (f_new - f) / (p_new' * df_new) / C;
            aguess = min(5, max(0.2, aguess));
        end
        df = df_new; df2 = df2_new; f = f_new; p = p_new;
    end

elseif strcmp(method, 'steepest'),
    % Steepest descent
    disp('Steepest descent');
    aguess = 1;
    while norm(df) > epsilon*(1 + abs(f)),
        p = - df;
        [ak, x, f_new, it] = linesch(x, p, aguess, C, f, df);
        disp(['guess = ' num2str(aguess) ...
              ', alpha_' int2str(ite) ' = ' num2str(ak)]);
        ite = ite + 1; % + it;
        df = objective(x, 1);
        aguess = (f - f_new) / (df' * df) / C;
        aguess = min(5, max(0.2, aguess));
        f = f_new;
    end
    ite
else
    error('Unknow method. Better be newton, conjugate, or steepest');
end

x
f
