%by ankur kumar
%grand tour animation - http://www.galaxy.gmu.edu/papers/tecrep.ps
%pollen data from: http://lib.stat.cmu.edu/datasets/
%dimensions of pollen data: ridge, nub, crack, weight, density
%first, graph the pollen data

%set up all the video files
fname='C:\Documents and Settings\labuser\Desktop\HiD\gt1.avi';
mov=avifile(fname);
mov.compression='cinepak';%compression is 'none' for unix
mov.quality=60;
fname2='C:\Documents and Settings\labuser\Desktop\HiD\gt2.avi';
mov2=avifile(fname2);
mov2.compression='cinepak';
mov2.quality=60;
%load the dataset
load 'C:\Documents and Settings\labuser\Desktop\data\pollen.csv';
pollen = pollen(:,[1:5]);
ndim = length(pollen(1,:));
kclust = 6; %number of clusters desired
[idx,c,sumd,d] = kmeans(pollen,kclust);
map = jet(kclust);
minx = sumd(1);
mini = 0;

%not sure how the dataset is pruned to reveal the internal structure
%because none of the research papers explain how they did this, and don't
%know how to do it without SGI machine, you can use density plots to do it
%with decent accuracy

%the devised pruning scheme used: find clusters, then find clusters in the found
%clusters.. this goes on 3 times until the data is pruned to about 100
%points and Grand Tour is run on this
for i=1:length(sumd)
    if(sumd(i)<minx)
        minx = sumd(i);
        mini = i;
    end
end
polleny = [];
figure(1)
hold on
%graph the entire dataset with clusters
for i=1:length(pollen)
    if(idx(i)==mini)
        polleny = [polleny;pollen(i,:)]; %store minimum sumd cluster
    end
    plot(pollen(i,:),'color',map(idx(i),:));
end
%the axes for parallel coords
axisx = [1;ndim];
axisy = [0;0];
label = [{'r'},{'n'},{'c'},{'w'},{'d'}];
for i = 1:ndim
    axisx = [axisx,[i;i]];
    axisy = [axisy,[-40;40]];
end
line(axisx, axisy, 'color', 'b');           %draw the axes
set(gca, 'xtick',[1:1:ndim]);
set(gca, 'xticklabel', label);
hold off

%filter out everything to get the internal structure using Brushing with
%Hue and Saturation - can't add pixel color values - no silicon graphics
%machine - need to think of another scheme to cut off noise and be able to
%use the concept of hue-saturation, use kmeans to find clusters - color
%them, then using sumd, get the min(sumd) cluster - possibly internal structure 
%will use density.m, density plots, to show internal structure, but
%reverting back to the points from the density plot is difficult since
%density plot is an image not a graph

%cluster
[idx2,c2,sumd2] = kmeans(polleny,kclust);
map = jet(kclust);
minx2 = sumd2(1);
mini2 = 0;
for i=1:length(sumd2)
    if(sumd2(i)<minx)
        minx2 = sumd2(i);
        mini2 = i;
    end
end
figure(2)
pollenx=[];
hold on
%graph in parallel coordinates
for i=1:length(polleny)
    if(idx2(i)==mini2)
        pollenx = [pollenx;polleny(i,:)]; %store minimum sumd cluster
    end
        plot(polleny(i,:),'color',map(idx2(i),:));
end
sprintf('length of polleny %d',length(polleny))
%the axes for parallel coords
axisx = [1;ndim];
axisy = [0;0];
label = [{'r'},{'n'},{'c'},{'w'},{'d'}];
for i = 1:ndim
    axisx = [axisx,[i;i]];
    axisy = [axisy,[-40;40]];
end
line(axisx, axisy, 'color', 'b');           %draw the axes
set(gca, 'xtick',[1:1:ndim]);
set(gca, 'xticklabel', label);
hold off

%---------change in dataset---------%
load gtdata.mat %-- for retreive.m
pollenx = newmat;
% pollenx = pollenx - repmat(mean(pollenx),length(pollenx),1);
%---------end change in dataset-------%

%further reduce data set (cluster) for easy viewing, and timing issues
[idx3,c3,sumd3] = kmeans(pollenx,kclust);
map = jet(kclust);
minx3 = sumd3(1);
mini3= 0;
for i=1:length(sumd3)
    if(sumd3(i)<minx3)
        minx3 = sumd3(i);
        mini3 = i;
    end
end
figure(3)
hold on
%graph in parallel coordinates
for i=1:length(pollenx)
        plot(pollenx(i,[1:5]),'color',map(idx3(i),:));
end
sprintf('length of pollenx %d',length(pollenx))
%the axes for parallel coords
axisx = [1;ndim];
axisy = [0;0];
label = [{'r'},{'n'},{'c'},{'w'},{'d'}];
for i = 1:ndim
    axisx = [axisx,[i;i]];
    axisy = [axisy,[-40;40]];
end
line(axisx, axisy, 'color', 'b');           %draw the axes
set(gca, 'xtick',[1:1:ndim]);
set(gca, 'xticklabel', label);
hold off



%run grand tour on pollenx
%http://www.galaxy.gmu.edu/papers/tecrep.ps
%Lambda is a sequence of mutually irrational real numbers to fill-up space 
% of a p-torus. This Lamda*t will become theta in the R matrix
%irrational numbers from: http://en.wikipedia.org/wiki/List_of_numbers#Transcendental_numbers
%10 lamdas for 0.5(d^2 - d) permutations for d=5

lambdaseq = zeros(5,5);
lambdaseq(1,2) = pi;
lambdaseq(1,3) = sqrt(2);
lambdaseq(1,4) = sqrt(3);
lambdaseq(1,5) = sqrt(5);
% lamdbaseq(1,6) = (1+sqrt(5))/2;
lambdaseq(2,3) = sqrt(12);
lambdaseq(2,4) = sqrt(6);
lambdaseq(2,5) = sqrt(7);
% lamdbaseq(2,6) = (sqrt(5)-1)/2;
lambdaseq(3,4) = sqrt(8);
lambdaseq(3,5) = sqrt(10);
% lamdbaseq(3,6) = sqrt(2) + 1;
lambdaseq(4,5) = sqrt(11);
% lamdbaseq(4,6) = nthroot(2,3);
% lamdbaseq(5,6) = nthroot(2,12);

syms t u %create a symbol for time
%construct the Q matrix 
r = [u 0 0 0 0;0 u 0 0 0;0 0 u 0 0;0 0 0 u 0;0 0 0 0 u];
%Q(lamda_1*t, lamda_2*t..) = R(lamda_1*t)xR(lamda_2*t)..
q = [u 0 0 0 0;0 u 0 0 0;0 0 u 0 0;0 0 0 u 0;0 0 0 0 u];%matrix Q is initialized
for i=1:ndim-1
    for j=i+1:ndim
            temp = lambdaseq(i,j)*t;
            r(i,i) = cos(temp);
            r(i,j) = -sin(temp);
            r(j,i) = sin(temp);
            r(j,j) = cos(temp);
            q = q*r;
            q=subs(q,u,1);
            r = [u 0 0 0 0;0 u 0 0 0;0 0 u 0 0;0 0 0 u 0;0 0 0 0 u]; %reset matrix R
    end
end
%now, e_j x Q(t) = a_j
e = eye(ndim);
time = linspace(0,2*pi,length(pollenx)); %time
y=zeros(length(pollenx),ndim);
qq=q;
size(qq)
size(pollenx(1,:))
fig4 = figure(4) %parallel coordinates
fig5 = figure(5) %dimension vs dimension, scatterplot matrix
% fig6 = figure(6) %surf plot
hold on
for k=1:length(time)
    %one timeframe t

    qq=subs(q,t,time(k))'; %q(t) transposed
    qq = mod(qq,2*pi); %modulo 2*pi
   
    for i=1:length(pollenx)
        x = pollenx(i,[1:5]);
        %aj = e(j,:)*qq;%aj(t)
        y(i,:)=x*qq;%inner product to give y_i(t)
    end
    k
    %y(t) has been calculated for all the data for some t
    %add frame to movie
    %plot y(t) in parallel coordinates
    %see the clusters of pollenx evolve in time
    clf(fig4);
    fig4;
    hold on
    for i=1:length(y)
            plot(y(i,[1:5]),'color',map(idx3(i),:));
    end
    %the axes for parallel coords
    axisx = [1;ndim];
    axisy = [0;0];
    label = [{'r'},{'n'},{'c'},{'w'},{'d'}];
    for i = 1:ndim
        axisx = [axisx,[i;i]];
        axisy = [axisy,[-40;40]];
    end
    line(axisx, axisy, 'color', 'b');           %draw the axes
    set(gca, 'xtick',[1:1:ndim]);
    set(gca, 'xticklabel', label);
    hold off
    %for graphing x1 vs x3 dimension for visualization
%     scatter(y(:,1),y(:,3));
%     set(get(gca,'XLabel'),'String','x1');
%     set(get(gca,'YLabel'),'String','x3');
    %movie
    drawnow;
    F = getframe(gcf);
    mov = addframe(mov,F); %the movie
    %scatterplot matrix - dimension vs dimension
    fig5;
    clf(fig5);
    hold on
    %plot all the scatter plots
    subplot(5,2,1)
    scatter(y(:,1),y(:,2));
    set(get(gca,'XLabel'),'String','x1');
    set(get(gca,'YLabel'),'String','x2');
    subplot(5,2,2)
    scatter(y(:,1),y(:,3));
    set(get(gca,'XLabel'),'String','x1');
    set(get(gca,'YLabel'),'String','x3');
    subplot(5,2,3)
    scatter(y(:,1),y(:,4));
    set(get(gca,'XLabel'),'String','x1');
    set(get(gca,'YLabel'),'String','x4');
    subplot(5,2,4)
    scatter(y(:,1),y(:,5));
    set(get(gca,'XLabel'),'String','x1');
    set(get(gca,'YLabel'),'String','x5');
    subplot(5,2,5)
    scatter(y(:,2),y(:,3));
    set(get(gca,'XLabel'),'String','x2');
    set(get(gca,'YLabel'),'String','x3');
    subplot(5,2,6)
    scatter(y(:,2),y(:,4));
    set(get(gca,'XLabel'),'String','x2');
    set(get(gca,'YLabel'),'String','x4');
    subplot(5,2,7)
    scatter(y(:,2),y(:,5));
    set(get(gca,'XLabel'),'String','x2');
    set(get(gca,'YLabel'),'String','x5');
    subplot(5,2,8)
    scatter(y(:,3),y(:,4));
    set(get(gca,'XLabel'),'String','x3');
    set(get(gca,'YLabel'),'String','x4');
    subplot(5,2,9)
    scatter(y(:,3),y(:,5));
    set(get(gca,'XLabel'),'String','x3');
    set(get(gca,'YLabel'),'String','x5');
    subplot(5,2,10)
    scatter(y(:,4),y(:,5));
    set(get(gca,'XLabel'),'String','x4');
    set(get(gca,'YLabel'),'String','x5');
    hold off
    %surf plot
%     fig6;
%     surf(y(:,[1:5]))
    %movie
    drawnow;
    F2 = getframe(gcf);
    mov2 = addframe(mov2,F2); %the movie
end
%close the movie files
mov = close(mov);
mov2 = close(mov2);
