function [alpha,b,sv] = ihlf_svr_rfntrain_fast(samples,targets,param)
% recursive finite Newton algorithm for support vector regression
% precompute the entire kernel matrix and the space complexity is O(n*n)
%
% input parameters
% samples: n*d matrix vector
% targets: n vector
% param:   aus parameters
%   param.kernel:       kernel type
%   param.kernelparam:  kernel parameter
%   param.lamda:        regular constant
%   param.epsilon:      insensitive parameter
%   param.delta:        Huber parameter
%
% output parameters
% alpha:  weight vector
% b:      threshold 
% sv:     support vector
%
% written by Liefeng Bo on 12/12/2005
% version 1.0

st = cputime;
itmax = 100;
[n,d] = size(samples);

% rearrange the samples
perm = randperm(n);
samples = samples(perm,:);
targets = targets(perm,1);

% recursively partition the samples
tempn = n;
nvec(1) = n;
i = 1;
while tempn > 1000
    i = i+1;
    tempn = round(tempn/2);
    nvec(i) = tempn;
end
nvec = fliplr(nvec);

% initialize the parameters
ksamples = evalkernel(samples,samples,param.kernel,param.kernelparam);
oldsv = 1:nvec(1);
oldsv1 = oldsv;
oldsv2 = oldsv;
oldalpha = (ksamples(oldsv,oldsv) + param.lamda*eye(length(oldsv)))\targets(1:nvec(1));
fprintf('\n');

% train the sub-problem
for i = 1:length(nvec)
    
    y = ksamples(1:nvec(i),oldsv)*oldalpha(oldsv);
    it = 1;
    
    while 1
        
        % find the support vectors
        r = y - targets(1:nvec(i));
        sv1 = find(abs(r) >= param.epsilon & abs(r) < param.delta);
        sv2 = find(abs(r) >= param.delta);
        nonsv = find(abs(r) < param.epsilon);
        sv = [sv1;sv2];
        
        % print the related information
        obj(it) = param.lamda*oldalpha'*y + sum((abs(r(sv1))-param.epsilon).^2)...
            + 2*(param.delta - param.epsilon)*sum(abs(r(sv2))) - length(sv2)*(param.delta.^2 - param.epsilon.^2);
        fprintf('n = %d, iter = %d, obj = %f, nb of sv = %d     \r',[nvec(i) it obj(it) length(sv)]);
        
        % stop
        if isempty(setxor(sv1,oldsv1)) & isempty(setxor(sv2,oldsv2)) | it >itmax
            b = 0;
            break;
        end
        
        % update the kernel matrix and targets
        H = ksamples(sv1,sv1) + param.lamda*eye(length(sv1));
        term = (param.delta - param.epsilon)*sign(r(sv2))/param.lamda;
        ttargets = targets(sv1) + param.epsilon*sign(r(sv1)) + ksamples(sv1,sv2)*term;
        
        % compute the step
        tempalpha = zeros(nvec(i),1);
        tempalpha(sv) = [H\ttargets; -term];
        tempy = ksamples(1:nvec(i),sv)*tempalpha(sv);
        h = tempalpha - oldalpha;
        
        % line search
        u = tempy - y;
        coff(1,1) = 2*(u(sv1)'*r(sv1) - param.epsilon*u(sv1)'*sign(r(sv1))...
            + (param.delta - param.epsilon)*u(sv2)'*sign(r(sv2)) + param.lamda*h'*y);
        coff(1,2) = 2*(u(sv1)'*u(sv1) + param.lamda*h'*u);
        t = linesearch(coff,r,u,param,nonsv,sv1,sv2);
        alpha = oldalpha + t*h;
        
        % update y and support vector
        y = t*tempy + (1 - t)*y;
        oldsv = sv;
        oldsv1 = sv1;
        oldsv2 = sv2;
        oldalpha = alpha;
        it = it + 1;
        
    end
    
    if i < length(nvec)
        oldalpha = [oldalpha;zeros(nvec(i+1) - nvec(i),1)];
    end
    fprintf('\n');
    
end

%output the weights and support vectors
alpha = alpha(sv);
sv = perm(sv);
[sv,index] = sort(sv);
alpha = alpha(index);
fprintf('Execution time : %4.1f seconds\n',cputime - st);

function t = linesearch(coff,r,u,param,nonsv,sv1,sv2)

% from zero loss to the quadratic loss, the linear loss
kv = [(sign(u(nonsv))*param.epsilon - r(nonsv))./u(nonsv) nonsv]; % enter the quadratic loss
kv = [kv; [(sign(u(nonsv))*param.delta - r(nonsv))./u(nonsv) -nonsv]]; % leave the quadratic loss

% from quadratic loss to the other loss
index1 = find(r(sv1).* u(sv1) > 0);
index1 = sv1(index1);
kv = [kv; [(sign(u(index1))*param.delta - r(index1))./u(index1) -index1]]; % enter the linear loss
index2 = find(r(sv1).* u(sv1) < 0);
index2 = sv1(index2);
kv = [kv; [(-sign(u(index2))*param.epsilon - r(index2))./u(index2) -index2]]; % enter the zero loss
kv = [kv; [(sign(u(index2))*param.epsilon - r(index2))./u(index2) index2]]; % enter the quadratic loss
kv = [kv; [(sign(u(index2))*param.delta - r(index2))./u(index2) -index2]]; % enter the linear loss

% from linear loss to the other loss
index = find(r(sv2).* u(sv2) < 0);
index = sv2(index);
kv = [kv; [(-sign(u(index))*param.delta - r(index))./u(index) index]]; % enters the quadratic loss
kv = [kv; [(-sign(u(index))*param.epsilon - r(index))./u(index) -index]]; % enters the zero loss
kv = [kv; [(sign(u(index))*param.epsilon - r(index))./u(index) index]]; % again enters the quadratic loss
kv = [kv; [(sign(u(index))*param.delta - r(index))./u(index) -index]]; % enters the linear loss

% sort the kink values
[temp ind] = sort(kv(:,1));
kv = kv(ind,:);

% compute the step length
kv = [kv; [Inf,0]];
leng = size(kv,1);
for i = 1:leng
    if  coff(1,1) + kv(i,1)*coff(1,2) > 0
        t = -coff(1,1)/coff(1,2); 
        break;
    else  % update coefficients
        if kv(i,2) > 0
            ind = abs(kv(i,2));
            coff(1,1) = coff(1,1) - 2*kv(i,1)*u(ind)^2;
            coff(1,2) = coff(1,2) + 2*u(ind)^2;
        else
            ind = abs(kv(i,2));
            coff(1,1) = coff(1,1) + 2*kv(i,1)*u(ind)^2;;
            coff(1,2) = coff(1,2) - 2*u(ind)^2;
        end
    end
end