function Y=an_kmeans(X,K)
% initialization X: Nxd
beta=0.008; [N d]=size(X);
mx=sum(X,1)/N;
Y = ones(K,1)*mx+rand(K,d)*0.005-0.0025;
hc=0;
while ~hc
    % A calculate cross distance of X and Y
    D=cross_dis(X,Y);
    % B determine overlapping memberships
    % V NxK  pp20
    expD=exp(-beta*D);
    sum_expD=sum(expD');
    V=expD./(sum_expD'*ones(1,K));
    % C update Y
    % Y Kxd  pp18
    for i=1:K
        Y(i,:)=V(:,i)'*X(:,:);
        if sum(V(:,i).^2) > 1/K*0.01
        Y(i,:)=Y(i,:)/sum(V(:,i));
        else
            Y(i,:)=Y(i,:)+rand(1,d)*0.005-0.0025
        end
    end
    beta = beta/0.9985;
    st = mean(sum(V.^2,2));
    if st < 1/K+0.005
        Y = Y+rand(K,d)*0.005-0.0025;
    end
    if st > 0.95
        hc=1;
    end
    fprintf(' %f  %f\n',beta, st);
end