% Example 10.1
myfun = @(x) (6*x-2).^2.*sin(12*x-4);                    % In-line function
xs=[0 0.5 0.68 1]';                                      % Sample locations
%xs=[0 0.246 0.5 0.68 1]';                               % Sample locations
%xs=[0 0.246 0.5 0.625 0.68 1]';                         % Sample locations
%xs=[0 0.246 0.5 0.625 0.68 0.743 1]';                   % Sample locations
%xs=[0 0.246 0.5 0.625 0.68 0.743 0.762 1]';             % Sample locations
%xs=[0 0.246 0.5 0.625 0.68 0.743 0.758 0.762 1]';       % Sample locations
ys=myfun(xs);                                              % output samples
xp=linspace(0,1,501)';                                  % Prediction points
[yp, ysig, nb]=Kriging(xs,ys,xp,0,0);             % Build Kriging surrogate
ny=size(ys,1);                                          % Number of samples
cilb = yp + tinv(0.16,ny-nb)*ysig;   % Lower-bounds of prediction intervals
ciub = yp + tinv(0.84,ny-nb)*ysig;   % Upper-bounds of prediction intervals
ytrue=myfun(xp);                % True function values at prediction points
figure; hold on;                                         % Plotting results
xShaded=[xp;sort(xp,'descend')]; yShaded=[cilb;ciub(end:-1:1)];
patch(xShaded,yShaded,[.24 .73 .89],'EdgeColor','none');
plot(xp,yp,'r-',xs,ys,'k+',xp,ytrue,'k-','Markersize',10,'Linewidth',1.5)
xlabel('x');ylabel('y(x)');
legend('68% PI','Kriging prediction','Samples','True function');
%
[xopt, yopt]=fminsearch(@myfun,0.6);             % Optimum of true function
[yPBS,I]=min(ys);xPBS=xs(I);                          % Present best sample
u=(yPBS-yp)./ysig;                                        % Normalized zPBS
EI=ysig.*(u.*cdf('Normal',u,0,1)+pdf('Normal',u,0,1));%Expected improvement
[EImax, I]=max(EI); xEImax=xp(I);            % Maximum expected improvement
fprintf('     xPBS     yPBS   xEImax    EImax\n')
fprintf('%8.3g %8.3g %8.3g %8.3g\n',xPBS,yPBS,xEImax,EImax);
figure;
plot(xp,EI,'k',xEImax,EImax,'pentagram','Markersize',10,'Linewidth',1.5);
xlabel('x');ylabel('EI(x)');
eRMS=sqrt(sum((ytrue-yp).^2)/501);             % RMS error at N grid points
eav=sum(abs(ytrue-yp))/501;              % Average error at NxN grid points
emax=max(abs(ytrue-yp));                 % Maximum error at NxN grid points
fprintf('RMS error:%8.3g\nAve error:%8.3g\nMax error:%8.3g\n',eRMS,eav,emax);
%-----------------------------------------------------------------------------
%
% Example 10.2
myfun = @(x) (6*x-2).^2.*sin(12*x-4);                    % In-line function
xs=[0 0.5 0.68 1]';                                         % Input samples
%xs=[0 0.5 0.633 0.68 1]';
%xs=[0 0.5 0.633 0.68 0.723 1]';
%xs=[0 0.244 0.5 0.633 0.68 0.723 1]';
%xs=[0 0.12 0.244 0.5 0.633 0.68 0.723 1]';
%xs=[0 0.12 0.244 0.5 0.633 0.68 0.723 0.794 1]';
ys=myfun(xs);                                              % output samples
xp=linspace(0,1,501)';                                  % Prediction points
[yp, ysig, nb]=Kriging(xs,ys,xp,0,0);             % Build Kriging surrogate
ny=size(ys,1);                                          % Number of samples
cilb = yp + tinv(0.16,ny-nb)*ysig;   % Lower-bounds of prediction intervals
ciub = yp + tinv(0.84,ny-nb)*ysig;   % Upper-bounds of prediction intervals
ytrue=myfun(xp);                 % True function values at predictio points
figure; hold on;                                         % Plotting results
xShaded=[xp;sort(xp,'descend')]; yShaded=[cilb;ciub(end:-1:1)];
patch(xShaded,yShaded,[.24 .73 .89],'EdgeColor','none');
plot(xp,yp,'r-',xs,ys,'k+',xp,ytrue,'k-','Markersize',10,'Linewidth',1.5)
xlabel('x');ylabel('y(x)');
legend('68% PI','Kriging prediction','Samples','True function');
%
[yPBS,I]=min(ys);xPBS=xs(I);                          % Present best sample
ytarget=yPBS-0.25*abs(yPBS);                                 % Target value
ztarget=(ytarget-yp)./ysig;                       % Normalized target value
PI=cdf('Normal',ztarget,0,1);                  % PI
[PImax, I]=max(PI); xPImax=xp(I);      % Max. of PI
fprintf('    xPBS     yPBS   xPImax    PImax  ytarget\n')
fprintf('%8.3g %8.3g %8.3g %8.3g %8.3g\n',xPBS,yPBS,xPImax,PImax,ytarget);
figure;
plot(xp,PI,'k',xPImax,PImax,'pentagram','Markersize',10,'Linewidth',1.5);
xlabel('x');ylabel('PI(x)');
eRMS=sqrt(sum((ytrue-yp).^2)/501);             % RMS error at N grid points
eav=sum(abs(ytrue-yp))/501;              % Average error at NxN grid points
emax=max(abs(ytrue-yp));               % Maximum error at NxN grid points
fprintf('RMS error:%8.3g\nAve error:%8.3g\nMax error:%8.3g\n',eRMS,eav,emax);
%-----------------------------------------------------------------------------
%
% Example 10.4
ytarget=yPBS-0.25*abs(yPBS);
while 1
 ztarget=(ytarget-yp)./ysig;
 PI=cdf('Normal',ztarget,0,1);
 [PImax, I]=max(PI); xPImax=xp(I);
 if PImax<0.05, ytarget=(ytarget+yPBS)/2; else, break; end
end
%-----------------------------------------------------------------------------
%
% Example 10.5
rng('default');                      % Initializing random number generator
xs=[0 0.25 0.5 0.68 1]';ys=EGOfn(xs);                % Initial sample pairs
N=100;xp=linspace(0,1,N)';                              % Prediction points
[xs, ys, yp, ysig, nb]=EGOPRS(xs,ys,xp,20,[0; 1],4);          % Calling EGO
ny=size(ys,1);                                  % Number of updated samples
cilb = yp + tinv(0.025,ny-nb)*ysig;  % Lower-bounds of prediction intervals
ciub = yp + tinv(0.975,ny-nb)*ysig;  % Upper-bounds of prediction intervals
ytrue=EGOfn(xp);                % True function values at prediction points
figure; hold on;                                         % Plotting results
xShaded=[xp;sort(xp,'descend')]; yShaded=[cilb;ciub(end:-1:1)];
patch(xShaded,yShaded,[.24 .73 .89],'EdgeColor','none');
plot(xp,yp,'r-',xs,ys,'k+',xp,ytrue,'k-','Markersize',10,'Linewidth',1.5)
xlabel('x');ylabel('y(x)');
legend('95% PI','Kriging prediction','Samples','True function');
[yPBS,I]=min(ys);
fprintf('Opt at x=');fprintf('%8.3g',xs(I,:));fprintf(', yPBS=%8.3g\n',yPBS);
eRMS=sqrt(sum((ytrue-yp).^2)/N);             % RMS error at NxN grid points
eav=sum(abs(ytrue-yp))/N;                % Average error at NxN grid points
emax=max(abs(ytrue-yp));                 % Maximum error at NxN grid points
fprintf('RMS error:%8.3g\nAve error:%8.3g\nMax error:%8.3g\n',eRMS,eav,emax);
%-----------------------------------------------------------------------------
%
% Example 10.6
function f=EGOfn(x)
 f = sin(x(:,1)).*cos(x(:,2));
end
%
rng('default'); N=21; dN=2*pi/(N-1); ny=20;
xs=2*pi*lhsdesign(ny,2,'iteration',5000);                  % 20 LHS samples
ys=EGOfn(xs);                                  % Samples from true function
[XX,YY]=meshgrid(-0:dN:2*pi, 0:dN:2*pi);          % NxN grid of predictions
xp=[reshape(XX,N^2,1) reshape(YY,N^2,1)];               % Prediction points
figure(1); scatter(xs(:,1),xs(:,2),'filled');hold on% Plot sample locations
axis([0 2*pi 0 2*pi]); xlabel('x_1');ylabel('x_2');
[xs, ys, yp, ysig, nb]=EGOPRS(xs,ys,xp,40,[0 0; 2*pi 2*pi],3);
YP=reshape(yp,N,N);                                   % Predictions at grid
YT=reshape(EGOfn(xp),N,N);                          % True function at grid
scatter(xs(ny+1:end,1),xs(ny+1:end,2),'d');     % Plot new sample locations
axis([0 2*pi 0 2*pi]); xlabel('x_1');ylabel('x_2');
figure(2); hold on; surf(XX,YY,YP);                   % Plot of predictions
plot3(xs(:,1),xs(:,2),ys,'ro');                           % Plot of samples
axis([0 2*pi 0 2*pi]); xlabel('x_1');ylabel('x_2');zlabel('y(x)');
figure(3); hold on; surf(XX,YY,YT)                  % PLot of true function
plot3(xs(:,1),xs(:,2),ys,'ro');                           % Plot of samples
axis([0 2*pi 0 2*pi]); xlabel('x_1');ylabel('x_2');zlabel('y(x)');
[yPBS,I]=min(ys);
fprintf('Opt at x=');fprintf('%8.3g',xs(I,:));fprintf(', yPBS=%8.3g\n',yPBS);
eRMS=sqrt(sum(sum((YT-YP).^2))/N^2);         % RMS error at NxN grid points
eav=sum(sum(abs(YT-YP)))/N^2;            % Average error at NxN grid points
emax=max(max(abs(YT-YP)));               % Maximum error at NxN grid points
fprintf('RMS error:%8.3g\nAve error:%8.3g\nMax error:%8.3g\n',eRMS,eav,emax);
%-----------------------------------------------------------------------------
%
% Exaxmple 10.7
xs=[0 0.5 0.68 1]';ys=EGOfn(xs);                     % Initial sample pairs
[xs, ys, yp, ysig, nb]=EGO(xs,ys,xp,0,20,[0; 1],0);
%-----------------------------------------------------------------------------
%
% Example 10.8
[xs, ys, yp, ysig, nb]=EGO(xs,ys,xp,0,40,[0 0; 2*pi 2*pi],0);
%-----------------------------------------------------------------------------