%author: ankur kumar
%function: perform bisecting k-means clustering on a dataset
%CSV files in format: Date, open, high, low, close, volume, adj close
%index of format:      1     2      3     4     5     6       7
%data from yahoo finance: oct-16-2006 to sept-14-2007
%use Lehman Brothers Data


function main()
clear all;
%import the dataset
dat = importfile('data-leh.csv');
tdow = [];
tdow = dat.data(:,[1:6]);
ndim = length(tdow(1,:));
map = jet(length(tdow));
figure(3);
%normalize all the prices for graphing purposes
tdow(:,1) = length(dat.data):-1:1;
scl = 10^6; %scaling the volume because of big number
tdow(:,2) = tdow(:,2);
tdow(:,3) = tdow(:,3);
tdow(:,4) = tdow(:,4);
tdow(:,5) = tdow(:,5);
tdow(:,6) = tdow(:,6)./scl; %volume
hold on
for i=1:length(tdow)
    plot(tdow(i,:),'color',map(i,:));
end
%graph the orginial dataset
%the axes for parallel coords
axisx = [1;ndim];
axisy = [0;0];
label = [{'t'},{'O'},{'H'},{'L'},{'C'},{'V'}];
for i = 1:ndim
    axisx = [axisx,[i;i]];
    axisy = [axisy,[0;100]];
end
line(axisx, axisy, 'color', 'b');           %draw the axes
set(gca, 'xtick',[1:1:ndim]);
set(gca, 'xticklabel', label);
hold off

%start cluster analysis - use bisecting kmeans find centroid
cluster.data = tdow;%will be updated
cluster.flag = zeros(1,length(cluster.data));%will be used to flagging clusters in the array
list.size = 0;
trials = 30;%number of iterations to find the clusters
targetk=4; %number of clusters desired
k = 2;%used in kmeans for bisection
numclust = 1;
sse=[];
for i=1:trials
    sprintf('iteration %d',i)
    [d1,d2]=calckmeans(cluster.data,cluster.flag,k);%calculate k-means
    size(d1)
    mycluster = bisect(cluster.data,cluster.flag,d1); %have all the clusters
    %find sse for all clusters
    length(d2(:,1))
    for j=1:length(d2(:,1))
        sse = [sse;calcSSEclusters(mycluster.data{j}.val,d2(j,:))]; %give arr and centroid into SSE
    end
    %find minimum sse
    [mn,v]=min(sse); %store index of sse, for the cluster index
    %remove the low-sse cluster from the data under examination by flagging
    cluster.flag = flagclust(cluster.flag,mycluster.data{v}.index);
    %add cluster with low sse to list of clusters
    list.clusters{numclust}.val = mycluster.data{v}.val; %store low-sse cluster in a list
    list.size = list.size+1;
    list.clusters{numclust}.sse = sse(v); %store the low-sse for that cluster
    numclust=numclust+1;
    %find if sse of leftover data is > sse of max in list of clusters
    %recompute centroid
    tempcent = centroid(cluster.data,cluster.flag);
    tempsse = calcSSE(cluster.data,cluster.flag,tempcent);
    %break after target # of clusters is reached
    if(list.size >= targetk)
        break;
    end
    [mx,w]=max(sse);
    %if a sse of a cluster in list of clusters > sse of points not in cluster, clustering is
    %not optimal, remove cluster from the list of clusters, recompute
    if(mx>tempsse)
        %merge data
        cluster.flag = unflagclust(cluster.flag,mycluster.data{w}.index);
        %remove from list of clusters, use sse
        for p=1:list.size
            if(list.clusters{p}.sse==mx)
                list.clusters{p}.val = [];%remove cluster from list of clusters
                list.clusters{p}.sse = 0;
                list.size = list.size - 1;
                break;
            end
        end
    end
    %reset variables
    sse=[];
end
%graph all the clusters in different colors
figure(4)
hold on
ii=1;
numpts = 0;
numclst=0;
map= jet(list.size);
for i=1:numclust-1
    if(list.clusters{i}.sse~=0)
        numclst = numclst+1;
        for j=1:length(list.clusters{i}.val(:,1))
            numpts=numpts+1;
            %parallel coordinates plot
            plot(list.clusters{i}.val(j,:),'color',map(numclst,:));
        end
        sprintf('cluster number %d',numclst)
        list.clusters{i}.val
    end
end
numpts
%the axes for parallel coords
axisx = [1;ndim];
axisy = [0;0];
label = [{'t'},{'O'},{'H'},{'L'},{'C'},{'V'}];
for i = 1:ndim
    axisx = [axisx,[i;i]];
    axisy = [axisy,[0;100]];
end
line(axisx, axisy, 'color', 'b');           %draw the axes
set(gca, 'xtick',[1:1:ndim]);
set(gca, 'xticklabel', label);
hold off
    

%calculate centroid without time base
function cent=centroid(arr,flags) %the matrix, and ndim
sum=0;
leng = 0;
for i=1:length(arr)
    if(flags(i)==0) %data is not flagged as cluster
        sum = sum+arr(i,:); %vector in ndim
        leng=leng+1; %number of unflagged points
    end
end
cent = sum./leng; %centroid

%calculate SSE for clusters without flags, this is for list of clusters
function ss=calcSSEclusters(barr,cntrd) %the matrix, and centroid
ss=0;
for i=1:length(barr(:,1))
    %calculate SSE
    ss = ss + sum((cntrd-barr(i,:)).^2);
end

%calculate SSE for nonclusters with flags
function ss=calcSSE(farr,flags,cntrd) %the matrix, and centroid
ss=0;
for i=1:length(farr)
    if(flags(i)==0) %data is not flagged as cluster
        ss = ss + sum((cntrd-farr(i,:)).^2);
    end
    %calculated SSE
end

%calculate kmeans clustering with bisection, k=2
%flags points as well
function [idx,idc]= calckmeans(arr,flags,k)
%create array of points for finding kmeans - noncluster points
temparr=[];
for i=1:length(arr)
    if(flags(i)==0)
        temparr=[temparr;arr(i,:)];
    end
end
[idxc,idc] = kmeans(temparr,k);
size(idxc)
idx = zeros(1,length(arr));
j=1;
%add in 0s for all cluster points based on flags in idx to simplify the other functions
for i=1:length(arr)
    if(flags(i)==0)
        %insert a 0 in the idx vector
        idx(i) = idxc(j);
        j=j+1;
    end
end
size(idx)


%bisect based on vector q
function clust = bisect(arr,flags,q)
for i=1:max(q)
    clust.data{i}.val = [];
    clust.data{i}.index = [];
end
for i=1:length(arr)
    if(flags(i)==0)
        clust.data{q(i)}.val=[clust.data{q(i)}.val; arr(i,:)];
        clust.data{q(i)}.index=[clust.data{q(i)}.index;i];
    end
end

%flag data based on indices in rarr and put flag in flags
function flags = flagclust(flags,rarr)
for i=1:length(rarr)
    flags(rarr(i)) = 1;
end

%unflag data based on indices in rarr and put flag in flags
function flags = unflagclust(flags,rarr)
for i=1:length(rarr)
    flags(rarr(i)) = 0;
end
