%cluster

%n = 11944;
%bs = 10000;

%replica = zeros(n,bs);

%disp('generate random numbers')
%for i=1:bs,
%    vrand = random('unid',n,n,1);
%    replica(:,i) = cook(vrand);
%end

disp('simple clustering');
for l= 1:7,
    vind = find(cook == l);
    vn = length(vind);
    
    vdist = zeros(vn,1);
    vdist(1) = vind(1)-1 + n-vind(vn);
    for i=2:vn,
        vdist(i) = vind(i)-vind(i-1);
    end

    cookdist(l) = sum(vdist)/vn;
 
    for j=1:bs,
        vind = find(replica(:,j) == l);
        vn = length(vind);
        
        if vn == 0,
            repldist(j,l) = n;
        else
            vdist = zeros(vn,1);
            vdist(1) = vind(1)-1 + n-vind(vn);
            for i=2:vn,
                vdist(i) = vind(i)-vind(i-1);
            end
            repldist(j,l) = sum(vdist)/vn;
        end
    end

    meandist(l) = sum(repldist(:,l))/bs;
    p05dist(l) = prctile(repldist(:,l),5);
    p95dist(l) = prctile(repldist(:,l),95);
end

disp('composite clustering')
vind = find((cook == 1) | (cook == 2));
vn = length(vind);
    
vdist = zeros(vn,1);
vdist(1) = vind(1)-1 + n-vind(vn);
for i=2:vn,
    vdist(i) = vind(i)-vind(i-1);
end

cookdist(8) = sum(vdist)/vn;

vind = find((cook == 6) | (cook == 7));
vn = length(vind);
    
vdist = zeros(vn,1);
vdist(1) = vind(1)-1 + n-vind(vn);
for i=2:vn,
    vdist(i) = vind(i)-vind(i-1);
end

cookdist(9) = sum(vdist)/vn;

vind = find((cook == 1) | (cook == 2) | (cook == 3));
vn = length(vind);
    
vdist = zeros(vn,1);
vdist(1) = vind(1)-1 + n-vind(vn);
for i=2:vn,
    vdist(i) = vind(i)-vind(i-1);
end

cookdist(10) = sum(vdist)/vn;

vind = find((cook == 6) | (cook == 7) | (cook == 5));
vn = length(vind);
    
vdist = zeros(vn,1);
vdist(1) = vind(1)-1 + n-vind(vn);
for i=2:vn,
    vdist(i) = vind(i)-vind(i-1);
end

cookdist(11) = sum(vdist)/vn;

vind = find((cook == 1) | (cook == 7));
vn = length(vind);
    
vdist = zeros(vn,1);
vdist(1) = vind(1)-1 + n-vind(vn);
for i=2:vn,
    vdist(i) = vind(i)-vind(i-1);
end

cookdist(12) = sum(vdist)/vn;

for j=1:bs,
    vind = find(replica(:,j) == 1 | replica(:,j) == 2);
    vn = length(vind);
        
    vdist = zeros(vn,1);
    vdist(1) = vind(1)-1 + n-vind(vn);
    for i=2:vn,
        vdist(i) = vind(i)-vind(i-1);
    end
    repldist(j,8) = sum(vdist)/vn;

    vind = find(replica(:,j) == 6 | replica(:,j) == 7);
    vn = length(vind);
        
    vdist = zeros(vn,1);
    vdist(1) = vind(1)-1 + n-vind(vn);
    for i=2:vn,
        vdist(i) = vind(i)-vind(i-1);
    end
    repldist(j,9) = sum(vdist)/vn;

    vind = find(replica(:,j) == 1 | replica(:,j) == 2 | replica(:,j) == 3);
    vn = length(vind);
        
    vdist = zeros(vn,1);
    vdist(1) = vind(1)-1 + n-vind(vn);
    for i=2:vn,
        vdist(i) = vind(i)-vind(i-1);
    end
    repldist(j,10) = sum(vdist)/vn;

    vind = find(replica(:,j) == 6 | replica(:,j) == 7 | replica(:,j) == 5);
    vn = length(vind);
        
    vdist = zeros(vn,1);
    vdist(1) = vind(1)-1 + n-vind(vn);
    for i=2:vn,
        vdist(i) = vind(i)-vind(i-1);
    end
    repldist(j,11) = sum(vdist)/vn;
    
    vind = find(replica(:,j) == 1 | replica(:,j) == 7);
    vn = length(vind);
        
    vdist = zeros(vn,1);
    vdist(1) = vind(1)-1 + n-vind(vn);
    for i=2:vn,
        vdist(i) = vind(i)-vind(i-1);
    end
    repldist(j,12) = sum(vdist)/vn;
      
end

for l=8:12,
    meandist(l) = sum(repldist(:,l))/bs;
    p05dist(l) = prctile(repldist(:,l),5);
    p95dist(l) = prctile(repldist(:,l),95);
end

clear v*