%                       Example 18.1.m and 18.4.m 
%
%   Code for Examples 18.1 and 18.4
%

  clear all
%  close all

%
% Set the coefficients and training points used to construct the Gaussian covariance function
% 
%        c(q,q') = exp[-(q-q')^2/2 \ell^2]
%
%

  sigma0 = 1e-14;
  gamma = 2;
  ell = 1;
  xlim = 4;
  dx = 0.05;
  x = -xlim:dx:xlim;
  n = length(x);

  xvals = [-3.0 -1.5 -0.2 1.3 2.9]';
  yvals = [-1.6 0.5 1.2 -1.0 -0.4]';
  k = length(yvals);

  xp = 2;

  for i=1:n
    for j=1:n
      C(i,j) = exp(-((x(i)-x(j))^gamma)/(2*ell^2));
    end
  end
  C = C + sigma0*eye(n,n); 

%
%  Construct the matrices used for training and predictions
% 

  for i=1:k
    for j=1:k
      C_t(i,j) = exp(-((xvals(i)-xvals(j))^gamma)/(2*ell^2));
    end
  end
  C_tinv = inv(C_t);

  for i=1:k
    for j=1:n
      C_s(i,j) = exp(-((xvals(i)-x(j))^gamma)/(2*ell^2));    
    end
  end 

  for i=1:n
    for j=1:n
      C_ss(i,j) = exp(-((x(i)-x(j))^gamma)/(2*ell^2));
    end
  end
  
  mean = C_s'*C_tinv*yvals;
  C_cond = C_ss - C_s'*C_tinv*C_s; 
  C_cond = (C_cond+C_cond.')/2 + sigma0*eye(n,n);
  V = diag(C_cond);
  SD2 = 2*sqrt(V);


  for i=1:k
    c_s(i,1) = exp(-((xvals(i)-xp)^gamma)/(2*ell^2));
  end
  c_ss = 1;
  c_cond = c_ss-c_s'*C_tinv*c_s;
  meanp = c_s'*C_tinv*yvals;


%
%  Compute the Cholesky decompositions C = R R^T, where R is lower triangular
%

  R = chol(C,'lower');
  R_cond = chol(C_cond,'lower');
  r_cond = sqrt(c_cond);

%
%  Sample n random variables z~N(0,1) and compute the prior functions f1, f2 and f3 and
%  posterior functions F1, F2, F3
%

  for j=1:3
    z = randn(n,1);
    f(:,j) = R*z;
  end

  for j=1:3
    z = randn(n,1);
    F(:,j) = mean + R_cond*z;
  end 

  for j=1:1000
    z = randn(1,1);
    g(j) = meanp + r_cond*z;
  end

%
%  Plot the prior functions 
%

  y1 = -2*ones(1,n);
  y2 = 2*ones(1,n);
  X = [x,fliplr(x)];
  Y = [y1,fliplr(y2)];

  figure(1)
  dimc = [0.85 0.85 0.85];
  fill(X,Y,dimc)
  hold on
  plot(x,f(:,1),'k-',x,f(:,2),'r--',x,f(:,3),'b-.','linewidth',2)
  hold off
  axis([-4 4 -2.5 2.5])
  xlabel('Input q')
  ylabel('Responses')
  set(gca,'Fontsize',[22]);
%  print -depsc GP_Prior2

%
%  Plot the posterior functions which are priors conditioned on the data
%   

  y3 = (mean+SD2)';
  y4 = (mean-SD2)';
  Y2 = [y3,fliplr(y4)];

  figure(2)
  fill(X,Y2,dimc,'HandleVisibility','off')
  hold on
  plot(x,mean,'r-',x,F(:,1),'k--',x,F(:,2),'k--',x,F(:,3),'k--','linewidth',2)
  plot(xvals,yvals,'+k','linewidth',5)
  plot([x(121),x(121)],[mean(121)-SD2(121),mean(121)+SD2(121)],'k-','linewidth',3) 
  hold off 
  axis([-4 4 -3.0 3.0])
  xlabel('Input q')
  ylabel('Responses')
  set(gca,'Fontsize',[22]);
  legend('Mean','Realizations','Location','South')
  %print -depsc GP_Posterior
   

  [bandwidth,density,mesh,cdf]=kde(g);
  figure(3)
  plot(mesh,density,'k-','linewidth',3)
  set(gca,'Fontsize',[22]);
  xlabel('f(q*)')
  ylabel('KDE of Responses')
  %print -depsc GP_Density