D = 2; %% Dimensions
N_PTS = [100 100 100 100]; N_GRPS = length(N_PTS); %% How many points per cluster, how many clusters
RADIUS = 0.8;
%DENSITY = 'uniform'; 
DENSITY = 'gaussian';

clrs = ['r'; 'g'; 'b'; 'k'; 'c']; nclrs = length(clrs);
    
%% Generate data and show them as well as original densities
X = [];
switch DENSITY,
case 'uniform',
    X = rand(sum(N_PTS),D);
    figure(1); plot(X(:,1),X(:,2),'.'); set(1,'DoubleBuffer','on');
    title('Ground truth is uniform');
case 'gaussian',
    MIX_TRUE = gmm(D, N_GRPS, 'full');
    %MIX_TRUE.centres = MIX_TRUE.centres * RADIUS;
    figure(1); clf;  set(1,'DoubleBuffer','on');
    for g=1:N_GRPS,
        clr = clrs(mod(g-1,nclrs)+1);
        SS = RADIUS*(rand(D,D)-0.5); MIX_TRUE.covars(:,:,g) = SS*SS';
        X1 = gsamp(MIX_TRUE.centres(g,:),MIX_TRUE.covars(:,:,g),N_PTS(g));
        X = [X ; X1];
        draw_ellipse(MIX_TRUE.centres(g,:),MIX_TRUE.covars(:,:,g),clr, MIX_TRUE.priors(g)*30);
        hold on;plot(X1(:,1),X1(:,2),['o' clr]);hold off;
    end;
title(['Ground truth is a mixture model with ' num2str(N_GRPS) ' components']);
end;
pause(4);



%% Starting from data try and estimate densities using gmmem
SHOW_ELLIPSES = 5;
model_types = {'full','diag', 'spherical'};
P_RATIO = [6 9 12 16];  %Ratio between n. of training examples and number of model parameters
model_params = [ D+D*(D+1)/2; 2*D; D+1];
N_GROUPS_TEST = ceil( 0.5*sum(N_PTS)./(model_params * P_RATIO));

OPTIONS = zeros(1,14); OPTIONS(1) = 0; OPTIONS(3) = 1; OPTIONS(5)=1; OPTIONS(14)=SHOW_ELLIPSES;
N_ITE = 100;

for m=1:length(model_types),
    for NG=N_GROUPS_TEST(m,:),
        i_g = find(NG==N_GROUPS_TEST(m,:));
        
        MIX = gmm(D, NG, model_types{m});
        MIX = gmminit(MIX, X(1:2:end,:), OPTIONS);
        loger = [];
        
        ERR_MIN = 0;
        for it=1:ceil(N_ITE/SHOW_ELLIPSES),
            figure(2); clf
            [MIX, OPTIONS, ERRLOG] = gmmem(MIX,X(1:2:end,:),OPTIONS);
            loger(it) = -OPTIONS(8);
            if it>1, if abs(loger(it)-loger(it-1))<OPTIONS(3), ERR_MIN=1; end; end;
            subplot(1,2,1);
            for g=1:NG,
                clr = clrs(mod(g-1,nclrs)+1);
                if MIX.covar_type(1:4)=='full', SS = MIX.covars(1:2,1:2,g); elseif MIX.covar_type(1:4)=='diag',SS = diag(MIX.covars(g,1:2)) ;else, SS = diag([1 1])*MIX.covars(1,g); end;
                draw_ellipse(MIX.centres(g,1:2),SS,clr, MIX.priors(g)*30);
            end;
            hold on;plot(X(:,1),X(:,2),'ko');hold off;
            title([num2str(NG) ' components with ' model_types{m} ' covariance']);
            subplot(1,2,2); plot(SHOW_ELLIPSES*[1:it],loger);
            title('Log likelihood');
            pause(0.1);
        end;
        pause(1);
        
        loglik_train(m,i_g) = sum(log(gmmprob(MIX,X(1:2:end,:))));
        loglik_test(m,i_g) = sum(log(gmmprob(MIX,X(2:2:end,:))));
    end;
end;

figure(100); clf; hold on;
ls = {'-', '--', '-.'};
for m=1:length(model_types),
N_PARAM = N_GROUPS_TEST(m,:)*(D+(D.*(D+1))/2);
P_RATIO1 = 0.5*sum(N_PTS)./N_PARAM;
semilogx(P_RATIO1,loglik_test(m,:),[ 'o' clrs(m) '-'],'LineWidth',2); plot_legend{2*m-1} = [model_types{m} ' test'];
semilogx(P_RATIO1,loglik_train(m,:),[ '*' clrs(m) '-.'],'LineWidth',2); plot_legend{2*m} = [model_types{m} ' train'];
end;
hold off;
legend(plot_legend);
title('Log likelihood as a function of model size');
xlabel('Ratio between n. training samples and n model params');
ylabel('Log likelihood');