function  [subspaceInfo, trainInfo] = dualGeneralFeaturesTrain(trainData, params)
%Train dual sparse general features - deflate on left side and all of
%kernel matrix 

if (nargin ~= 2)
    fprintf('%s\n', help(sprintf('%s', mfilename)));
    error('Incorrect number of inputs - see above usage instructions.');
end

%First, figure out which variables to use in the data struct 
[nameX, nameY] = getSpaceNames(params);
useSparse = issparse(getDataFieldValue(trainData, nameX)); 

[numTrainExamples, numFeatures] = getDataFieldSize(trainData, nameX); 
alpha = 10^-8; %Added to the diagonal of matrices to make them non singular
tol = 10^-6; 

%Store all the parameters 
dualFeatureDirection = char(params.dualFeatureDirection); 
T = min(params.iterations, numTrainExamples); 
kernelFunction = char(params.X.kernel.name); 
kernelParams = params.X.kernel.params; 

if isfield(params, 'normalise')
    normaliseFeatures = params.normalise;
else 
    normaliseFeatures = 1; 
end

%doubleDeflation = 0 for left sided, 1 for double 
if isfield(params, 'doubleDeflation') 
    doubleDeflation = params.doubleDeflation; 
else 
    doubleDeflation = 0; 
end 

%Select best dual direction from a subset of the kernel matrix columns 
if isfield(params.X, 'kernelCols') 
    kernelCols = params.X.kernelCols; 
else 
    kernelCols = numTrainExamples;
end 


trainK = feval(kernelFunction, getDataFieldValue(trainData, nameX), getDataFieldValue(trainData, nameX), kernelParams); 

Kj = trainK; 

if ~strcmp(nameY, '')
    Yj = getDataFieldValue(trainData, nameY);
else
    Yj = rand(numTrainExamples, 1); 
end

%Our dual projection directions 
b = zeros(numTrainExamples, T); 
tau = zeros(numTrainExamples, T); 
s = zeros(numTrainExamples, T); 
r = zeros(T, 1); 
normSqTau = zeros(T, 1); 

if useSparse
    nonZeros = nnz(trainK);
    KbbK = spalloc(numTrainExamples, numTrainExamples, nonZeros);
    KKbbK = spalloc(numTrainExamples, numTrainExamples, nonZeros);
end 

fprintf('Iterating ... \n');
tic;

%Compute the projection directions
for i=1:T
    displayCount(i, T, 30);
    
    if ~strcmp(nameY, '')
        b(:, i) = feval(dualFeatureDirection, trainK, Kj, getDataFieldValue(trainData, nameY), Yj);
    else
        %Bit of a hack, but if this is the case, then use subset of columns
        if kernelCols ~= numTrainExamples;
            b(:, i) = feval(dualFeatureDirection, trainK, Kj, Yj, Yj, kernelCols);
        else
            b(:, i) = feval(dualFeatureDirection, trainK, Kj, Yj, Yj);
        end
    end
    
    tau(:, i) = Kj*b(:, i);
    normSqTau(i) = tau(:, i)'*tau(:, i);
    
    if useSparse  
        s(:, i) = sparse(trainK'*tau(:, i))/normSqTau(i);
        KbbK = sparse(tau(:, i)*tau(:, i)');
        KbbKK = sparse(KbbK*Kj/normSqTau(i));
        projMatrix = sparse(KbbK)/normSqTau(i);
    else
        s(:, i) = trainK'*tau(:, i)/normSqTau(i); 
        KbbK = tau(:, i)*tau(:, i)';
        KbbKK = KbbK*Kj/normSqTau(i);
        projMatrix = KbbK/normSqTau(i);
    end
    
    r(i) = s(:, i)'*tau(:, i)/normSqTau(i); 
      
    %Sometimes norm of tau(:, i) becomes too small, so no point in carrying on 
    if abs(normSqTau(i)) < tol
        fprintf('\nFinished early at projections %d, because norm(tau) is close to zero.\n', i);
        b = b(:, 1:i);
        s = s(:, 1:i);
        tau = tau(:, 1:i);
        normSqTau = normSqTau(1:i);
        T = i; 
        break;
    end
   
    %If we deflate as follows we get a slight difference (strange) 
    %Kj = Kj - (tau(:, i)*tau(:, i)')*Kj/(tau(:, i)'*tau(:, i)); 
    %Yj = Yj - (tau(:, i)*tau(:, i)')*Yj/(tau(:, i)'*tau(:, i));
    if doubleDeflation == 0
        Kj = Kj - KbbKK;
    else
        Kj = Kj - KbbKK  - KbbKK' + KbbKK*KbbK/(normSqTau(i));
    end
    
    Yj = Yj - projMatrix*Yj;
   
    %Clear up the large temporary variables
    clear KbbK KKbbK projMatrix; 
end

trainTime = toc; 
fprintf('Completed in %f seconds\n', trainTime); 

clear Kj; 

if normaliseFeatures == 1
    normMatrix = diag(1./sqrt(normSqTau));
else
    normMatrix = eye(T);
end

%If we have double deflated, we need to deflate b too 
if doubleDeflation ~= 0
    for i=1:T
        for j=1:i-1
            b(:, i) = b(:, i) - tau(:, j)'*b(:, i)*tau(:, j)/normSqTau(j);       
        end 
    end 
end 
    
%Compute new features on training and test data 
Q = tau'*trainK*b + alpha*eye(T);
Q = b/(diag(1./normSqTau)*Q ); 
Q = Q*normMatrix; 

%Now save output variables 
trainInfo = struct; 
trainInfo.data = data; 
trainInfo.data = addDataField(trainInfo.data, 'X', tau*normMatrix, 'examples');  

subspaceInfo = struct; 
subspaceInfo.trainTime = trainTime; 
subspaceInfo.X.numFeatures = T; 
subspaceInfo.X.tau = tau; 
subspaceInfo.X.b = b; 
subspaceInfo.X.r = r; 
subspaceInfo.X.s = s; 
subspaceInfo.X.Q = Q; 
subspaceInfo.X.normSqTau = normSqTau;