% 070412 gabriel@pads.ufrj.br. Source is Example3_3_2_4.m.
% Comparacao: GLA, SA e DA para K = 16 classes.

close all;

% [1] GLA

clear all; K = 16; P = 100;
randn('state',0); a = randn(1,1); rand('state',0);
X0 = randn(2,K);
X=[];
for k=1:size(X0,2),   
    X = [X repmat(X0(:,k),1,P)+randn(2,P)*0.1];    
end;

figure; plot(X(1,:),X(2,:),'k.',X0(1,:),X0(2,:),'b.'); grid on; hold on; voronoi(X0(1,:),X0(2,:)); [Cn,p,Dmin] = ExemploFuncaoJ(X,X0); Dmin = Dmin/size(X,2)

C0 = mean(X,2);
D0 = mean(sum((X-repmat(C0,1,size(X,2))).^2,1));
C = C0; Di = D0; e = 1e-3; delta = 1e-5; 
n = 0; fim = 0; split = 1; 

while not(fim),
    n = n + 1; BKD(n,:) = [size(C,2) Di];
    if split,
        Phi = e*randn(size(C)); C = [C-Phi C+Phi];
    end;
    split = 0;
    [Cn,p,D] = ExemploFuncaoJ(X,C);
    for k = 1:length(p),
        C(:,k) = Cn(:,k)/p(k); p(k) = p(k)/size(X,2);
    end;
    D = D/size(X,2);
    if abs((D-Di)/Di) < delta,
        if size(C,2) < K, split = 1; else fim = 1; end;
    end;
    Di = D;
end;

figure; plot(1:n,BKD(:,2)); grid on; 

figure; plot(X(1,:),X(2,:),'k.',C(1,:),C(2,:),'b.'); grid on; hold on; voronoi(C(1,:),C(2,:)); plot(C0(1,:),C0(2,:),'r.'); D

% [2] SA

clear all; K = 16; P = 100;
randn('state',0); a = randn(1,1); rand('state',0);
X0 = randn(2,K);
X=[];
for k=1:size(X0,2),   
    X = [X repmat(X0(:,k),1,P)+randn(2,P)*0.1];    
end;

C0 = 4*(rand(2,K)-0.5); [Cn,p,D0] = ExemploFuncaoJ(X,C0); Ci = C0; Di = D0/size(X,2);
e = 5e-2; t = 1; T0 = 5e-3; T = T0/log2(1+t); N = 10000; L = 20;
n = 0; fim = 0; Df = Di; Cf = Ci; BKD = zeros(N*(L-1),2);

while not(fim),
    n = n + 1; BKD(n,:) = [T Di];
    C = Ci + e*randn(size(Ci));
    [Cn,p,D] = ExemploFuncaoJ(X,C); D = D/size(X,2);
    if D<Df, Df = D; Cf = C; end;
    if D<Di, 
        Ci = C; Di = D;
    else
        a = rand(1,1); 
        if a < exp((Di-D)/T), Ci = C; Di = D; end;
    end;
    if (rem(n,N)==0),
        [T Df]
        t = t+1;
        T = T0/log2(1+t);
        if t==L, fim = 1; end;
    end;    
end;
C = Cf; 

figure; plot(1:n,BKD(:,2)); grid on; figure; plot(1:n,BKD(:,1)); grid on;

figure; plot(X(1,:),X(2,:),'k.',C(1,:),C(2,:),'b.'); grid on; hold on; voronoi(C(1,:),C(2,:)); plot(C0(1,:),C0(2,:),'r.'); axis([-2 2.5 -2 2]); Df

% [3] DA

clear all; K = 16; P = 100;
randn('state',0); a = randn(1,1); rand('state',0);
X0 = randn(2,K);
X=[];
for k=1:size(X0,2),
    X = [X repmat(X0(:,k),1,P)+randn(2,P)*0.1];
end;
C0 = randn(2,K); T0 = 0.2; Delta = 0.00001; epsilon = 1e-6; alpha = 0.98;

C = C0; T = T0; fim = 0; i=0; n=0;
while not(fim),

    i=i+1;

    % Partition Condition
    for k=1:size(X,2),
        for l=1:size(C,2),
            d = sum((X(:,k)-C(:,l)).^2);
            p_ygivenx(l,k) = exp(-d/T);
            d_aux(l,k) = d;
        end;
    end;
    Zx = sum(p_ygivenx,1);
    p_ygivenx = p_ygivenx./repmat(Zx,size(C,2),1);

    % Cost Functions
    J(i) = (-T/size(X,2))*sum(log(Zx));
    D(i) = mean(sum(p_ygivenx.*d_aux,1));
    LocalT(i) = T;

    % Centroid Condition
    C_new = zeros(size(C));
    for k=1:size(C,2),
        y = zeros(size(C,1),1); w = 0;
        for l=1:size(X,2),
            y = y + p_ygivenx(k,l)*X(:,l);
            w = w + p_ygivenx(k,l);
        end;
        C_new(:,k) = y/w;
    end;

    C = C_new;
    if (i>2),
        if (((abs(J(i)-J(i-1))/abs(J(i-1)))*100)<Delta), 
            T = alpha*T
            n=n+1; BKC(n).C = C; BKC(n).T = T;
            C = C.*(1+epsilon*randn(size(C)));
        end;
    end;
    if T<0.1, fim = 1; end;

end;

figure; plot(X(1,:),X(2,:),'k.'); hold on; plot(X0(1,:),X0(2,:),'c.'); plot(C0(1,:),C0(2,:),'r.'); voronoi(C(1,:),C(2,:),'b.'); grid on;
title('DA - 16 classes');
figure; plot(D,'b.'); hold on; plot(LocalT,'k.'); grid on; axis([0 length(D) 0 0.25]); xlabel('i'); ylabel('D(i) em azul. T(i) em preto');
for k=1:(length(BKC)-1),
    V(k) = mean(sum((BKC(k+1).C-BKC(k).C).^2,1)); T(k) = BKC(k).T;
end;
figure; plot(T,V,'k-'); hold on; plot(T,V,'ko'); axis([0.1 0.2 0 0.01]); grid on; xlabel('T'); ylabel('\Delta C(T)');

[a,b,c] = ExemploFuncaoJ(X,C); D_SA = c/size(X,2)