% Active learning, bayesian optimization and experimental design
% 18th Machine Learning Summer School
%
% Author: Ruben Martinez-Cantin <rmcantin@unizar.es>
%
clear all, close all
seed = 197; randn('seed',seed), rand('seed',seed)


p = @(x) 1./(1+exp(-5*sin(x*10)));                  % "true" underlying probability
criteria = @simple_criterion;
bounds = [0, 1];

% Set initial set of points using latin hypercube sampling
sn = 0.005;
nInit = 20;    nStart = 5;  nIter = 20;
xtr = lhsu(bounds(1),bounds(2),nInit);   
ytr = 2*(p(xtr)>rand(nInit,1))-1;                          % draw labels +1/-1
i = randperm(nInit); nout = 2;                               % add outliers
ytr(i(1:nout)) = -ytr(i(1:nout)); 

% setup the GP
cov = 
mean = 
lik = 
inf = 
hyp0 =


% Once we have learned the hyperparameters we reset the training points 
% to make the problem more interesting. Never do this in a real problem!
xtr = lhsu(bounds(1),bounds(2),nStart);   
ytr = 2*(p(xtr)>rand(nStart,1))-1;                          % draw labels +1/-1

Ncg = 50;
hyp = minimize(hyp0,'gp', -Ncg, inf, mean, cov, lik, xtr, ytr); % opt hypers


% set up DIRECT
opts.maxevals = 1000;
opts.maxits = 200;
opts.showits = 0;
Problem.f = criteria;

for i = 1:nIter
    
    fprintf(1,'- Iter %i of %i \n', i,nIter);
    
    % Sometimes, we can anneal the parameters. Thus, we include the
    % iteration number as a temporal index for annealing.
    params.iter = i;
    
    [minval, newX, hist] = Direct(Problem, bounds, opts, hyp, inf, mean, cov, lik, xtr, ytr, params);
    newY = 2*(p(newX)>rand(1,1))-1;                       
    
    xtr = [xtr; newX];    ytr = [ytr; newY];
    
    fprintf(1,'New evaluation point. x = %2.4f, y = %2.4f \n', newX, newY);

    
    hyp = minimize(hyp0,'gp', -Ncg, inf, mean, cov, lik, xtr, ytr); % opt hypers

    % =====================================================================
    % PLOT RESULTS    =====================================================
    % =====================================================================
    x = linspace(0,1,1e4)';
    y = 2*(p(x)>rand(1e4,1))-1;
    [crit, ymu, ys] = feval(criteria, x, hyp, inf, mean, cov, lik, xtr, ytr, params );
    
    col = {'k',[.8,0,0],[0,.5,0],'b',[0,.75,.75],[.7,0,.5]};                % colors
    lw = 2;
    h = figure(1);
    clf;hold on;
    
    subplot(3,1,[1 2]);   hold on;
    title('Gaussian Process');
    plot(xtr,ytr,'ko','LineWidth',lw);
    plot(newX,newY,'ro','LineWidth',lw);
    
    plot(x,2*p(x)-1, 'Color', col{2},'LineWidth',lw);
    plot(x, ymu, 'Color', col{1},'LineWidth',lw);
    
    % Plot std values
    sdscale = 0.5;
    ysd = sdscale*ys;
    fill([x;flipud(x)],[ymu+ysd;flipud(ymu-ysd)],...
      col{1},'EdgeColor',col{1},'FaceAlpha',0.1,'EdgeAlpha',0.3);

    avalue = axis;   axis([0 1 avalue(3) avalue(4)]);
    
    subplot(3,1,3); hold on;
    title('Criteria');
    bar(x,-crit,'k');
    avalue = axis;   axis([0 1 avalue(3) avalue(4)]);
    
    pause(0.2);
end