%c1da = [(randn(25,1))-2 .5*randn(25,1)];
%c1db = [(randn(25,1))-5 .5*randn(25,1)-1];
%c1dc = [(randn(25,1))-5 .5*randn(25,1)-2];
%c1dd = [(randn(25,1))-2 .5*randn(25,1)-3];
%c1d =  [c1da;c1db;c1dc;c1dd];
%c2d = randn(100,2)+1.5;

%load /stuff/s1.point;
%load /stuff/s2.stop;
%point = s1;
%stop = s2;


nclust = 4;
choose = 20;

pointsize = size(cpoint);
stopsize = size(cstop);
A = [];
ncprekern = [];
prekern = [];
chosen = [];
whichclust=[];
c1d =[];
c2d =[];
selected = 0;
chosen = zeros(pointsize(1)+1, 1);
while(selected < choose)
	pick = floor(rand(1)*pointsize(1))+1;
	while(chosen(pick) == 1)
		pick = floor(rand(1)*pointsize(1))+1;
	end
	newone = cpoint(pick, :);
	chosen(pick) = 1;
	selected = selected+1;
	c1d = [c1d; newone];
end

selected = 0;
chosen = zeros(stopsize(1)+1, 1);
while(selected < choose)
	pick = floor(rand(1)*stopsize(1))+1;
	while(chosen(pick) == 1)
		pick = floor(rand(1)*stopsize(1))+1;
	end
	newone = cstop(pick, :);
	chosen(pick) = 1;
	selected = selected+1;
	c2d = [c2d; newone];
end


%c1d = point(1:200, 2:23);
%c2d = stop(1:200,2:23);

data = [c1d;c2d];

d1size = size(c1d);
d1len = d1size(1);
d2size = size(c2d);
d2len = d1size(1);
dsize = size(data);
dlen = dsize(1);
path1 = [];
path2 = [];


% gamma(i,j) is mean i projected onto phi(x_j)
% kernel(a, b) is the kernel function
% M(i,j) is the assignment matrix.  It equals 1 if x_i is nearest mean j
% distance metric is as equation (39)
% Update:
% 	zeta = (is the new point in cluster)/(number of points in cluster)
% 	gamma(a,j) = gamma(a,j)*(1-zeta) (except for gamma(max, :) = zeta)

M = zeros(dsize(1), nclust);
gamma = zeros(nclust, dsize(1));
dist = zeros(nclust,dsize(1));

datamean = mean(data);
dm = ones(dsize(1),1)*datamean;
data = data - dm;

disp('Starting kernel calculations...');
for i = 1:dsize(1)
	for j = i:dsize(1)
		tmp = kernel(data(i,:), data(j,:));
		ncprekern(i,j) = tmp;
		ncprekern(j,i) = tmp;
	end
end
disp('Kernel calculations complete, centering...')
fM = ones(dsize(1))*(1./dsize(1));
prekern = ncprekern - fM*ncprekern - ncprekern*fM + fM*ncprekern*fM;
disp('High dimensional centering complete.')
oldM = zeros(dsize(1), nclust);
for t = 1:nclust
		gamma(t,:) = zeros(1,dsize(1));
		gamma(t,t) = 1;
		M(t,t) = 1;
end

for it = 1:50
it
for t = 1:dsize(1)
	tmp = 0;

% Faster, better, matlab friendlier
% The first term of the distance value for cluster c, point t
% A = gamma(c,:)'*gamma(c,:);
% Fist term:  sum(sum(A.*prekern))
% Second term: -2*sum(gamma(c,:).*prekern(t,:))

% The first term is independent of t, except to the extent that
% as t changes, gamma will.  Does this matter?

	A = gamma(1,:)'*gamma(1,:);
	tmp = sum(sum(A.*prekern)) - 2*sum(gamma(1,:).*prekern(t,:));

	best = tmp;
	bestcluster = 1;
	dist(1,t) = tmp;
	for c = 2:nclust
		tmp = 0;
		A = gamma(c,:)'*gamma(c,:);
		tmp = sum(sum(A.*prekern)) - 2*sum(gamma(c,:).*prekern(t,:));
		dist(c,t) = tmp;
		if(tmp < best)
			best = tmp;
			bestcluster = c;
		end
	end
	M(t, :) = zeros(1,nclust);
	M(t, bestcluster) = 1;
	whichclust(t) = bestcluster;
	[t bestcluster]
%	t
%	bestcluster
	zeta = 1./(sum(M(1:t, bestcluster)));

	gammacor = ones(1, dsize(1)).*(1-zeta);
	gamma(bestcluster,:) = gamma(bestcluster,:).*gammacor;
	gamma(bestcluster,t) = gamma(bestcluster,t) + zeta;
end
if(norm(M-oldM) == 0)
	M
	sum(M(1:choose,:))./sum(M)
	sum(M(choose+1:2*choose,:))./sum(M)
%	kmerr
	kmcv
	return
end
M
oldM = M;
end



return;





c1 = [0 0];
c2 = [.5 -.5];

clust = ones(nclust,1)*zeros(1,dsize(2));

for it = 1:8
it
newc = ones(nclust, 1)*zeros(1,dsize(2));
cct = zeros(nclust,1);

for n = 1:nclust
	clust(n,:) = data(n,:);
end

for p = 1:dlen
	d = norm(data(p,:)-clust(1,:));
	dmin = d;
	whichcluster = 1;
	for n = 1:nclust
		d = norm(data(p,:)-clust(n,:));
		if(d < dmin)
			dmin = d;
			whichcluster = n;
		end
	end
	newc(whichcluster,:) = newc(whichcluster,:)+data(p,:);
	cct(whichcluster) = cct(whichcluster)+1;
end

for n = 1:nclust
	clust(n,:) = newc(n,:)./cct(n);
end
end
hold off;
plot(clust(:,1), clust(:,2), 'magentao');
hold on;
plot(c1d(:,1), c1d(:,2), 'red*');
plot(c2d(:,1), c2d(:,2), 'blue*');

c1d = point;
c2d = stop;

data = [c1d;c2d];

dsize = size(data);
dlen = dsize(1);
d1size = size(c1d);
d1len = d1size(1);
d2size = size(c2d);
d2len = d1size(1);
path1 = [];
path2 = [];

for p = 1:dlen
	d = norm(data(p,:)-clust(1,:));
	dmin = d;
	whichcluster = 1;
	for n = 1:nclust
		d = norm(data(p,:)-clust(n,:));
		if(d < dmin)
			dmin = d;
			whichcluster = n;
		end
	end
	categ(p) = whichcluster;
%	text(data(p,1)+.08, data(p,2), num2str(whichcluster));
end


return;

hold off;
plot(0,0,'white.');
hold on;
for p = 1:d2len
	d1 = norm(c2d(p,:)-c1);
	d2 = norm(c2d(p,:)-c2);
	if(d1 < d2)
		plot(c2d(p,1), c2d(p,2), 'black+');
	else
		plot(c2d(p,1), c2d(p,2), 'black*');
	end
end
for p = 1:d1len
	d1 = norm(c1d(p,:)-c1);
	d2 = norm(c1d(p,:)-c2);
	if(d1 < d2)
		plot(c1d(p,1), c1d(p,2), 'red*');
	else
		plot(c1d(p,1), c1d(p,2), 'red+');
	end
end
plot(path1(:,1), path1(:,2), 'red')
plot(path2(:,1), path2(:,2), 'black')
plot(c1(:,1), c1(:,2), 'redo');
plot(c2(:,1), c2(:,2), 'blacko');
c1
c2
