% 071120 gabriel@pads.ufrj.br

function [Yout,p,D,k] = XYltoYpDk(X,Yin,l,lambda);

p = zeros(size(Yin,2),1); Y = Yin; Ynew = zeros(size(Y)); D = 0;
k = zeros(1,size(X,2));
for n=1:size(X,2), 
	% Partition
	j = sum((repmat(X(:,n),1,size(Y,2)) - Y).^2,1) + lambda*l; 
	k(n) = min(find(j==min(j))); 
	D = D + sum((X(:,n)-Y(:,k(n))).^2);
	% Centroid
	Ynew(:,k(n)) = Ynew(:,k(n)) + X(:,n);
	p(k(n)) = p(k(n)) + 1;
end; 
D = D/size(X,2);
for j=1:size(Yin,2), 
	if p(j)~=0, Yout(:,j) = Ynew(:,j)/p(j); end; 
end;
