function GPopt_it2d;

% Setting up grids
tv=[0.02:0.02:1]';
n1=size(tv,1);
sv=[0.02:0.02:1]';
n2=size(sv,1);
t=tv*ones(1,n2);
s=ones(n1,1)*sv';
n=n1*n2;

% True function 
x=1.5*cos(15*s)+0.5*sin(10*t)-cos(15*s.*t);
figure(4); clf; 
axis xy;
imagesc(x);
impixelinfo;

% Vectorized values
tA=reshape(t,n,1);
sA=reshape(s,n,1);
% evaluations stored in vector
xA=reshape(x,n,1);


% ------------ Iterative evaluation -------


% Iterate
iter=0;
m0=20;

ITERMAX=10;
while (iter<ITERMAX)

    % 
    % Design new evaluation points 
    if (iter==0)
        %
        % Sample random locations to initialize
        %
        ii=randsample(n,m0);
        tB=tA(ii,1); 
        sB=sA(ii,1);
        % Evaluate
        xB=xA(ii,1);
    else
        % New data points:
        [eindi,indi]=sort(-eiAN);
        ii=indi(1:m0,1);
        tB=[tB;tA(ii,1)]; 
        sB=[sB;sA(ii,1)];
        xB=[xB;xA(ii,1)];
    end
     

    %
    % GP pred
    [mpost,Spost]=gp_update([sA tA],[sB tB],xB);
    sv=sqrt(diag(Spost));
  
    figure(5);
    clf;
    imagesc(reshape(mpost,n1,n2));
    axis xy;
    colorbar;
    title('Prediction');
    figure(6);
    clf;
    imagesc(reshape(sv,n1,n2));
    axis xy;
    colorbar;
    title('Std dev');

    vstar=max(xB);
    vstarVEC=vstar*ones(n,1);
    
    % Expected improvement
    eiAN=expectedimprovement(mpost,sv,vstarVEC);
    figure(7);
    clf;
    imagesc(reshape(eiAN,n1,n2));
    axis xy;
    colorbar;
    title('Expected improvement');
    

    % Monitor EI at each iteration
    iter=iter+1;
    emmax(iter)=max(eiAN);
    emmean(iter)=mean(eiAN);
    disp(sprintf('Iteration %d. max EI is %0.3g. mean EI is %0.3g.',iter,max(eiAN),mean(eiAN)));
    pause(1);

end;
figure(8);
clf;
plot([1:iter],emmax,'k-x',[1:iter],emmean,'r');
legend('EImax','EImean');
xlabel('Iteration'); ylabel('Expected improvement');


%%---------------------

function [mpost,Spost]=gp_update(xa,xb,yb);

TTa=size(xa,1);
ttua=0.1;
mm=mean(yb);
sig=std(yb);
TTb=size(xb,1);
ma=mm*ones(TTa,1);
mb=mm*ones(TTb,1);
phimm=10;
% b set
Hb=formdist(xb,xb);
Sb=sig^2*(1+phimm*Hb).*exp(-phimm*Hb)+ttua^2*eye(TTb); % Added observation noise
Hab=formdist(xa,xb);
Sab=sig^2*(1+phimm*Hab).*exp(-phimm*Hab);
Ha=formdist(xa,xa);
Sa=sig^2*(1+phimm*Ha).*exp(-phimm*Ha);
mpost=ma+Sab*(Sb\(yb-mb));
Spost=Sa-Sab*(Sb\Sab');

function H=formdist(xa,xb);

Ta=size(xa,1);
Tb=size(xb,1);
Hb1=abs(xa(:,1)*ones(1,Tb)-ones(Ta,1)*xb(:,1)');
Hb2=abs(xa(:,2)*ones(1,Tb)-ones(Ta,1)*xb(:,2)');
H=sqrt(Hb1.^2+Hb2.^2);

function eiAN=expectedimprovement(mux,sigx,at);

eiAN=(mux-at).*normcdf((mux-at)./sigx)+sigx.*normpdf((mux-at)./sigx);

