clear
load fisheriris
X = meas;
m  = size(X,2);
 % load machine
%  load census
% % X = meas;
% X=X(1:2000,:);
d = pdist(X,'euclidean');

Z = linkage(d,'ward');%Create a hierarchical binary cluster tree using linkage
% [Hr,Tr]=dendrogram(Z, 0);  %generates a dendrogram plot of the hierarchical binary cluster tree
% h_gca = gca;
% h_gca.TickDir = 'out';
% h_gca.TickLength = [.002 0];
% h_gca.XTickLabel = [];
% c = cluster

val0=eval(vpa(Z(:,3),5));

t=1; n(1)=0; val=val0(1); f(1)=val;
for i=1:length(val0)
    if val0(i)==val
        n(t)=n(t)+1;
        continue
    else
        t=t+1;
        val=val0(i);
        n(t)=1;
        f(t)=val;
        continue
    end
end
figure(2)
fsum=cumsum(n);
x= length(val0)-fsum;
plot(x, f, 'o-','LineWidth', 2)
xlabel('Number of Clusters')
ylabel('Merge Distance')
title('A sample evalution graph')

% for i = 1: length(x)
%     y= cluster(Z, x(i));
%     B = Amean(X, y)
% end

x =x(end:-1:1);
f=f(end:-1:1);

i=0;
flag = 1;  b =  length(x); minerror_c_old = b-1;
while flag ~=0 & b>=20
    RMSE = [];
    for c = 2: (b-2)
        x1 = x(1:c); y1 = f(1 : c);
        x2 = x(c:b); y2 = f(c : b);
        [~, s1] = polyfit(x1, y1,1);
        [~, s2] = polyfit(x2, y2,1);
        RMSE(c-1) = ((c-1)/(b-1))*s1.normr + ( (b-c)/(b-1))*s2.normr;
    end
     [~, minerror_c ]= min(RMSE);
     minerror_c= minerror_c+1;
     if  minerror_c >= minerror_c_old
         break
     else
         minerror_c_old = minerror_c;
%           b = b -1;
          b = 2*minerror_c ;
     end
     i=i+1
end

hold on 
plot(x( minerror_c), f( minerror_c),'ro')
c =x( minerror_c)

labels = cluster(Z, c);

labs = unique(labels);
CORR = zeros();
for i = 1:length(labs)
    corri = cov(X(labels == i, :)); 
    CORR = CORR +  corri; 
end