clear all;
close all;
clc;
% f sim GP(m,k) here f is a function this is a Gaussian Process
% the function f is distributrd as a GP woth mean m and covariance 
% function k
% m(x)=1/4*x^2
% k(x,x')=exp(-1/2 / l^2 *(x-x')^2)


% f sim N(mu,Sigma) here f is a vector this is a Gaussian distribution

% --FROM RASMUSSEN
figure;
hold on;
xs = (-5:0.2:5)'; ns = size(xs,1); keps = 1e-9;
n = randperm(ns,10); 
m = inline('0.25*x.^2'); % mu = 1/4*xi.^2;
K = inline('exp(-0.5/ 1.2^2 * ( repmat(p'',size(q)) - repmat(q,size(p'')) ).^2)'); % Sigma = exp(-1/2*pdist2(xi,xi).^2);
mu =  m(xs);
Sigma = K(xs,xs)+keps*eye(ns);
fs = mu + chol(Sigma)' * randn(ns,1);
plot(xs,fs,'--')
fs = mu + chol(Sigma)' * randn(ns,1);
plot(xs,fs,'--')
fs = mu + chol(Sigma)' * randn(ns,1);
plot(xs,fs,'--')
lb = mu - 2*sqrt(diag(Sigma)) ;
ub = mu + 2*sqrt(diag(Sigma)) ;
color = 'b';
edgecolor ='k';
transparency = 0.2;
fill([xs; flipud(xs)], [lb; flipud(ub)], color, 'EdgeColor',edgecolor,'FaceAlpha',transparency,'EdgeAlpha',transparency);

figure;
imagesc(Sigma);
colorbar

for p=1:length(n)
    figure;
    hold on;
    x = xs(n(1:p))+0.01;
    f = fs(n(1:p));
    plot(x,f,'ro')
    nx = size(x,1)
    a = m(xs);
    b = m(x); 
    A = K(xs,xs)+keps*eye(ns);
    B = K(x,xs)+keps*eye(ns,nx);
    C = K(x,x)+keps*eye(nx);
    mu = a+B*inv(C)*(f-b);
    Sigma = A-B*inv(C)*B';
    ff = mu + chol(Sigma)' * randn(ns,1);
    plot(xs,ff,'--')
    ff = mu + chol(Sigma)' * randn(ns,1);
    plot(xs,ff,'--')
    ff = mu + chol(Sigma)' * randn(ns,1);
    plot(xs,ff,'--') 
    lb = mu - 2*sqrt(diag(Sigma)) ;
    ub = mu + 2*sqrt(diag(Sigma)) ;
    color = 'b';
    edgecolor ='k';
    transparency = 0.2;
    fill([xs; flipud(xs)], [lb; flipud(ub)], color, 'EdgeColor',edgecolor,'FaceAlpha',transparency,'EdgeAlpha',transparency);

    figure;
    imagesc(Sigma);
    colorbar
end