function [r,c,coreset,info] = meb(P,epsilon);

% Input: P is a matrix whose columns are the coordinates of each of
% the points. P is a (d) by n matrix.
% Assumes number of points to be greater than 2

% Output: Center c and the radius r of the MEB.

% Read input from file
mytime = cputime;
[d,n] = size(P);

% Approximate the diameter
max = 0.0;
index = 0;
v = 0.0;

for i = 2:n
   v = norm(P(:,i) - P(:,1));
   if v > max
      max = v;
      index = i;
   end
end

fprintf( 'Max = %f', max)

for i = 1:d
   temp = P(i,1);
   P(i,1) = P(i,index);
   P(i,index) = temp;
end

max = 0.0;
for i = 2:n
   v = norm(P(:,i) - P(:,1));
   if v > max
      max = v;
      index = i;
   end
end

% Now index , 1 are a diameter approximation
% Q is the intial core set of size 2
Q = [P(:,1) P(:,index)]

k = 1 / epsilon;
ittime = 0;

% Setup the constant for the random sampling function
constant_f = sqrt(d);
check_all  = 0; 	% Right now doesnt check everything
						% to find a violator

% Begin the main iteration loop

for j = 1:k
   
%	[r,c,info] = meb1(Q,epsilon);
	[r,c] = mebsdpt(Q,epsilon);

	itinit = cputime;
	maxnorm = 0.0;
	maxindex = -1;
   
   % Checking output
   
   while (maxindex == -1) & (check_all == 0)
                        	% Check a subset of points to find a violator
      counter = 1;
      start_from = int32(double(vpa(rand,10)) * n) 
   	steps_todo = int32(constant_f * sqrt(n))
   
   	while counter < steps_todo
      	if (counter+double(start_from)) >= (n+0.5)
         	start_from = 0;
     	 	end
      
      	v = norm(P(:,double(counter+double(start_from))) - c);
      	if v > maxnorm
      		maxnorm = v;
      		maxindex = double(counter+double(start_from));
      	end
      	counter = counter + 1;
      end   
      
      if maxindex == -1
         constant_f = constant_f * 2;
      else
         constant_f = constant_f + 1;
         fprintf('violator found\n');
      end
      
      if (uint32(constant_f*sqrt(n))) > (n-1)
            check_all = 1;
      end

   end 	% while

   if check_all == 1   	% Check all points for finding a violator
		for i = 1:n
   		v = norm(P(:,i) - c);
   		if v > maxnorm
      		maxnorm = v;
      		maxindex = i;
   		end
      end
   end
   
   fprintf('Norm , r = %f %f\n',maxnorm,r);

	if (maxnorm < (r + r * epsilon)) & (check_all == 1)
   	fprintf('Done\nTotal Time :');
   	fprintf('%d sec\n', cputime - mytime);
	   fprintf('Iteration Time : %d sec\n',ittime);
   	fprintf('Number of Iterations : %d\n',j);
	   break;
   elseif (maxnorm < (r + r * epsilon)) 
   	check_all = 1;   
   end

	Q = [Q P(:,maxindex)];
	ittime = ittime + (cputime - itinit);
   coreset = Q;

end