function [xk, f, xl] = trust_region (obj, x0, eta, reg, reg_max, epsilon)
%TRUST_REGION   Trust regin Newton CG method (ACM 113, Problem 3.5)
%
%   [xmin, fmin, xlist] = trust_region(obj, x0, eta, reg, reg_max, epsilon)
%
%       obj:         Objective function
%                    (obj(x, n) is the n-th deriv of obj at x)
%       x0:          Initial guess of x_min
%
%       eta:         Sufficient trust parameter (0<eta<1/4)
%       reg:         Initial region size
%       reg_max:     Maximal region size
%       epsilon:     epsilon for outer Newton iteration.
%
%   By default, eta = 0.2; reg = 1; reg_max = 5; epsilon = 1e-12.
%
%   See also: TRUNC_NEWTON.

%   Ling Li, Caltech
%   05/05/2001

if nargin < 3 | eta < 0 | eta > 1/4,
    eta = 0.2;
end
if nargin < 4 | reg < 0,
    reg = 1;
end
if nargin < 5,
    reg_max = 5;
end
if nargin < 6,
    epsilon = 1e-12;
end

% Trust region control
%% Initialize
if reg > reg_max,  %% reg in (0, reg_max)
    reg = reg_max;
end

out_iter = 0;
xk = x0; xl = xk;
f = eval([obj '(xk, 0)']); df = eval([obj '(xk, 1)']);

while norm(df) > epsilon * (1+abs(f)),
    A = eval([obj '(xk, 2)']);
    A = (A + A') / 2;

    %% Inner CG:     1  T        T
    %%    min m(p) = - p A p + df p, where A is symmetric
    %%               2
    %%    subject to |p| <= reg.

    %% xi: CG solution; pi: CG search direction
    xi = zeros(size(df));
    ri = df; %% A*xi - df;
    pi = - ri;

    ndf = norm(df);
    ndf = min(0.5, sqrt(ndf)) * ndf;
    while norm(ri) > ndf,
        Api = A * pi;
        piApi = pi' * Api;
        if piApi <= 0,
            %% Find tau that minimize m(xi + tao*pi)
            tau = bordertau(xi, pi, reg);
            xi = xi + tau * pi;
            break;
        end

        ri2 = ri' * ri;
        ai = ri2 / piApi;
        xin = xi + ai * pi;
        if norm(xin) >= reg,
            %% Find tau s.t. |xi + tao*pi| == reg
            %% Note xi is orthogonal to pi
            tau = bordertau(xi, pi, reg);
            xi = xi + tau * pi;
            break;
        end
        xi = xin;
        ri = ri + ai .* Api;
        bi = ri' * ri / ri2;
        pi = bi * pi - ri;
    end
    %% Inner linear CG: END

    pk = xi;

    fn = eval([obj '(xk+pk, 0)']);
    est = df'*pk + pk'*A*pk/2;
    if est <= 0,
        rho = 1;
    else
        rho = (fn - f) / est;
    end
    if rho < 1/4,
        reg = norm(pk) / 4;
    elseif rho > 3/4 & norm(pk) >= 0.9999*reg,
        reg = min(2*reg, reg_max);
    end

    if rho > eta,
        xk = xk + pk;
        f = fn;
        df = eval([obj '(xk, 1)']);
    end

    out_iter = out_iter + 1;
    xl = [xl, xk];
end

disp('Result: ');
disp(['  x_min = [' num2str(xk') ']'' with norm = ' num2str(norm(xk)) '.']);
disp(['  f_min = ' num2str(f)]);
disp(['  iterations = ' int2str(out_iter)]);


function tau = bordertau (p, d, reg)
pd = p' * d;
dd = d' * d;
pp = p' * p;

tau = (sqrt(pd*pd + dd*(reg*reg - pp)) - pd) / dd;
