% constructs gaussian processes w/training
%

clear
clc
clf

rand(2);

%%%%%%%%%%%%%%%
% NEED TO ADD MATERN KERNEL
%%%%%%%%
%k = @(x,y) 1*x'*y;
k= @(x,y) exp(-((x-y)'*(x-y))/(2.*0.1*0.1));

% test points
Xtest= (0:.001:1);
n = length(Xtest);
N = 10;   %number of GP samples

% noiseless training data

Xtrain= [0.2 ; 0.49 ; 0.59];  % location
ytrain= [0.1 ; -0.4 ; -0.33];  %  sample value

% construct covariance matrix

for i=1:length(Xtrain)
    for j=1:length(Xtrain)
        K11(i,j)= k(Xtrain(i),Xtrain(j));
    end
end

for i=1:length(Xtrain)
    for j=1:length(Xtest)
        K12(i,j)= k(Xtrain(i), Xtest(j));
    end
end

K21= K12';

for i=1:length(Xtest)
    for j=1:length(Xtest)
        K22(i,j)= k(Xtest(i), Xtest(j));
    end
end

K = [K11 K12; K21 K22];

% mean function
m= K21*inv(K11)*ytrain;

C = zeros(n,n);
for i=1:n
    for j=1:n
        C(i,j)= k(Xtest(i),Xtest(j));
    end
end

C= C - K21*inv(K11)*K12;

u = randn(n,1);
[A, S, ~]= svd(C);
z = m + A*sqrt(S)*u;  % z = A S^0.5 u

%plot results
plot(Xtest,z,'.-')   % first sample
axis([0, 1, -3, 3])

hold on

for i=2:N % remaining samples
    u = randn(n,1);
    z = m + A*sqrt(S)*u;  % z = A S^0.5 u
    plot(Xtest,z,'.-')
end

% plot observations
for i=1:length(Xtrain)
    plot(Xtrain(i),ytrain(i),'o','MarkerFaceColor','k', 'MarkerSize', 13);
end


