% Generates some random points and then fits the smallest hypersphere that
% contains them all

% Iain Murray 2003

% Useful excercise in understanding mathematical programming and duality

% This relies on a quadprog solver being available, this might involve adding
% something like:
% addpath('/opt/matlab6.5/toolbox/optim')
% to this code

% params
N = 20;
d = 2; % number of dimensions (d==2 assumed in some of the plotting stuff)
D = randn(d,N);

% quadratic program on dual problem
Lambda=quadprog(2*D'*D,-sum(D.*D,1),[ones(1,N);-eye(N)],[1;zeros(N,1)]);

% solving for answer in primal problem
c=(D*Lambda)/sum(Lambda);
r=sqrt(max(sum((D-repmat(c,1,N)).*(D-repmat(c,1,N)),1)));

% All that follows is just plotting stuff
allangles=[0:0.01:2*pi];
line(1,:)=c(1)+r*cos(allangles);
line(2,:)=c(2)+r*sin(allangles);
clf; hold on;
plot(D(1,:),D(2,:),'.');
plot(line(1,:),line(2,:),'g');
axis square;