% 081014 gabriel@pads.ufrj.br

lambda = [0 1e-3 2e-3 5e-3 1e-2 2e-2 5e-2 1e-1 2e-1 5e-1 1 2 3];
% lambda=0;
HDJ = zeros(length(lambda),3);
for n=1:length(lambda),
    randn('state',0);
    X = randn(2,1000);
    D0 = mean(sum((X-repmat(mean(X,2),1,size(X,2))).^2,1));
    C = randn(2,32);
    fim = 0; jf=50; j=0;
    D=zeros(jf,1); H=zeros(jf,1);
    l=ones(1,size(C,2))/size(C,2);
    while (fim==0),
        j=j+1; C_new = zeros(size(C)); p = zeros(size(C,2),1); D(j)=0;
        for k=1:size(X,2),
            d = sum((repmat(X(:,k),1,size(C,2))-C).^2,1);
            jk = d + lambda(n)*l;
            i = min(find(jk==min(jk)));
            D(j) = D(j) + d(i);
            C_new(:,i) = C_new(:,i)+X(:,k);
            p(i) = p(i)+1;
        end;
        C=[]; l=[]; i=0; D(j) = D(j)/size(X,2); H(j) = 0;
        for k=1:size(C_new,2),
            if p(k)~=0,
                i=i+1; C(:,i) = C_new(:,k)/p(k); p(k) = p(k)/size(X,2);
                l(i) = -log2(p(k));
                H(j) = H(j) - p(k)*log2(p(k));
            end;
        end;
        if j==jf,
            fim=1;
        end;
    end;    
    HDJ(n,:) = [H(jf) D(jf) D(jf)+lambda(n)*H(jf)];    
    [n H(jf) D(jf) D(jf)+lambda(n)*H(jf)]
end;
plot(HDJ(:,1),HDJ(:,2)/D0,'k+'); hold on; grid on; axis([0 5 0 1]);