% EXPERIMENTAL DESIGN FOR LINEAR MODEL USING A AND D OPTIMALITY.
% Written by Ruben Cantin and Nando de Freitas.
clear; echo off;

% PARAMETERS
% ==========
theta = [2 1]';           % parameters: intercept and slope
[d,tmp] = size(theta);    % d=number of features.
n = 40;                   % number of points.
sigma = 2.0;              % Noise variance.
R = [1 0;0 1];            % Prior information matrix.
theta0 = [2 0]';          % Prior mean.
nqueries = 10;            % Number of queries.
A = [1 0;0 1];            % Matrix for A-Designs.

% GENERATE DATA
% =============
x=zeros(n,1);
y=zeros(n,1);
for i=1:n,
 x(i)=i-n/2;
 y(i)= [1 x(i)]*theta + sigma*randn;
end;
%save data x y;
%load data;

% PLOT DATA AND GENERATING LINE
% =============================
figure(1)
clf;
line1X = [x(1):.05:x(n)]';
line1Y = [ones(length(line1X),1) line1X]*theta;
plot(x,y,'bo','linewidth',2);
hold on;
plot(line1X,line1Y,'k.');
xlabel('x','fontsize',15); ylabel('y','fontsize',15);
figure(3)
clf;
[xx,yy] = meshgrid(-1.4:.2:4,-1.4:.2:4);
[r,c] = size(xx);
x3 = reshape(xx,r*r,1);
y3 = reshape(yy,r*r,1);
xy = [x3 y3];
prior = mvnpdf(xy,theta0',inv(R));
prior = reshape(prior,r,r);
contour(xx,yy,prior,1);
xlabel('Intercept \theta_1');
ylabel('slope \theta_2');
hold on;


% EXPERIMENTS
% ===========
nleft = n;         % Number of points left to query.
e = ceil(n*rand);  % first point to query.
xin = x(e);        % x points that have been queried already.
yin = y(e);
figure(1)
plot(xin,y(e),'r+','linewidth',2);
legend('Unknown data','Generating line','Experiment','location','southeast');
pause;
xleft = [x(1:e-1);x(e+1:nleft)]; % Points that still must be queried.
yleft = [y(1:e-1);y(e+1:nleft)];
nleft = nleft -1;
while nleft>n-nqueries, % While there's points to be queried.
 u=zeros(nleft,1);      % Utility of experiment.
 for e=1:nleft,         % For all possible experiments.
   X = [ones(n-nleft,1) xin; 1 xleft(e)];
   XX = X'*X;
   u(e) = log(prod(eigs(XX+R))*(sigma^(-d))); % D-Design.
   %u(e) = - trace(inv(A*(XX+R)))*(sigma^(-d)); % A-Design.
 end;
 
 figure(2);
 plot(xleft,u,'.');
 xlabel('input x=e','fontsize',15); ylabel('u(e)','fontsize',15);
 axislimits = axis; axislimits(1) = x(1); axislimits(2) = x(n);
 axis(axislimits);
 pause;
 
 [uebest,ebest] = max(u);    % Pick best design (x) point.
 xin = [xin; xleft(ebest)];
 yin = [yin; yleft(ebest)];
 xleft = [xleft(1:ebest-1); xleft(ebest+1:nleft)];
 yleft = [yleft(1:ebest-1); yleft(ebest+1:nleft)];
 nleft = nleft-1;
 X = [ones(n-nleft,1) xin];
 XX = X'*X;
 sigmaP = inv(XX + R);
 varP = (sigma^2)*sigmaP;            % Posterior variance.
 meanP = sigmaP*(X'*yin + R*theta0); % Posterior mean.
 
 figure(1)
 plot(xin,yin,'r+','linewidth',2);
 
 figure(3)
 posterior = mvnpdf(xy,meanP',varP);
 posterior = reshape(posterior,r,r);
 contour(xx,yy,posterior,1);
 hold on;
 pause;
end;
figure(1)
line2X = line1X;
line2Y = [ones(length(line1X),1) line2X]*meanP;
plot(line2X,line2Y,'g.');