% Course "Kernel methods"
%
% Homepage : http://www.bioinfo.ensmp.fr/~vert/teaching/2007master
%
% Copyright 2007 Jean-Philippe Vert
%
%
%
% Homework 1
%
% Kernel PCA and Kernel ridge regression
%
% Files required:
%   rbf.m
%   centering.m
%   kernelPCA.m
%   krr.m


% Load the data
load xtrain.txt
load ytrain.txt
load xtest.txt
load ytest.txt

ntrain = size(xtrain,1);
ntest = size(xtest,1);
xall = [xtrain;xtest];

% Create a few kernel matrices
kname = {'Linear','Gaussian (sigma = 10)','Gaussian (sigma = 1)','Gaussian (sigma = 0.1)'}
kmat{1} = xall * xall';
kmat{2} = rbf(10 , xall);
kmat{3} = rbf(1 , xall);
kmat{4} = rbf(0.1 , xall);

nk = length(kmat); % number of kernels

%%%%%%%%%%%%%%%%%%%%%%%%%%
% Exercice 1: kernel PCA %
%%%%%%%%%%%%%%%%%%%%%%%%%%

% Check the code of kernelPCA.m for details

for ik = 1:nk
    % projet all points onto the first 2 eigenvectors computed on the
    % training set only
    proj = kernelPCA( kmat{ik} , [1:ntrain] , 2 );
    % Visualize
    plot(proj(1:ntrain,1),proj(1:ntrain,2),'ob',proj((ntrain+1):end,1),proj((ntrain+1):end,2),'*r')
    legend('Train','Test')
    xlabel('PC1')
    ylabel('PC2')
    title(kname{ik})
    input('Next picture?')
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Exercice 2: kernel ridge regression with offset %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% check krr.m for details

% the lambdas to be tested. It is a good idea to search a regular grid in
% the log because we don't know a priori what would be a good scale for
% lambda.
nlambda = 20;
lambdas = logspace( -6, 6, nlambda );
trainerror = zeros(1,nlambda);
testerror = zeros(1,nlambda);

for ik = 1:nk
    for i = 1:nlambda
        % run kernel ridge regression
        [ alphas , b ] = krr( lambdas(i) , kmat{ik}( 1:ntrain , 1:ntrain ) , ytrain );
        % make prediction on all data
        ypred = kmat{ik}( : , 1:ntrain ) * alphas + b;
        % compute the trainiing and test errors
        trainerror(i) = sum( ( ypred(1:ntrain) - ytrain ).^2 ) / ntrain;
        testerror(i) = sum( ( ypred((ntrain+1):end) - ytest ).^2 ) / ntest;
    end
    semilogx(lambdas,trainerror,'b',lambdas,testerror,'r')
    legend('Training error','Test error')
    title(kname{ik})
    xlabel('Lambda')
    ylabel('Mean square error')
    input('Next picture?')
end
