function Y=annealed_kmeans(X,K)
[N d]=size(X);
mean_x = mean(X);
Y=randn(K,d)*0.1+ones(K,1)*mean_x;
change=1; v=ceil(rand(N,1)*size(Y,1))';
while change > 0
    D=cross_dis(X,Y);
    v_old = v;
    [dd v]=min(D');
    mean(dd)
    for k=1:K
        ind=find(v == k);
        if length(ind) > 0
        Y(k,:) = mean(X(ind,:));
        end
    end
    change=length(find(v~=v_old));
end



function D=cross_dis(X,Y)
K=size(Y,1);N=size(X,1);
A=sum(X.^2,2)*ones(1,K);
C=ones(N,1)*sum(Y.^2,2)';
B=X*Y';
D=sqrt(A-2*B+C);