% kernel
%

clear;
N=100;
syms u v;
g=u*sin(2*pi*(u+v)) + v*sin(3*pi*(u-v));

% generate input (data) matrix, X
X=rand(2,N);

w=normrnd(0,sqrt(0.1),1,100);

% generate output vector, y
for i=1:N,
  u=X(1,i);
  v=X(2,i);
  y(i)=eval(g)+w(i);
end
y=y';

% make a test data set, t
% by withholding last 10
for k=1:10,
    for i=1:2,
        t(i,k)=X(i,k+(N-10));
    end
end


%% train kernel

% step 1: define kernels
% a) linear
% b) gaussian
% c) other

K_lin=X'*X;     % linear kernel
              
C_g=1;          % gaussian = exp(-C||x_1 - x_2||^2)
K_gauss= exp(- C_g * (repmat(sum(X.^2,1),N,1) + repmat(sum(X.^2,1),N,1)' - 2*X'*X) );

K

% step 2: compute c = (lambda * identity + K)^-1 * y

% pick lambda between 10^-5 and 1
lambda = 10^-3;

c_lin=((lambda*ones(N,N)+K_lin)^1)*y;
c_gauss=((lambda*ones(N,N)+K_gauss)^1)*y;

%% step 3: given a data set X, find f(x) = sum ( c(i) * K(X(i),X) )

% compute T, which is f(x) = sum...
for j=1:10,    
    for i=1:N-10,
        a=[ X(1,i) X(2,i) ]';
        b=[ t(1,j) t(2,j) ]';
        
        T_partial_lin(i) = c_lin(i) * a'*b;
        T_lin(j) = sum(T_partial_lin);
        
        T_partial_gauss(i) = c_gauss(i) * exp(- C_g * ( a'*a + b'*b - 2*a'*b ) );
        T_gauss(j) = sum(T_partial_gauss);
    end        
end


% how far off is f(x) from y?  should be within variance of w

for j=1:10,
  u=t(1,j);
  v=t(2,j);
  y_t(j)=eval(g);
end

clf;    
xaxis=linspace(1,10,10);
plot(xaxis,T_lin,'go',xaxis,T_gauss,'bd',xaxis,y_t,'rx');
title('kernel');
legend('f(x) linear','f(x) gaussian','y');