function [index,c]=TSD_clustering(X,beta);

%TSD Data set clustering using Transitional State Discrimination (TSD) algorithm with two states.
%
%   The TSD algorithm groups elements according to the transition of their consecutive time points. 
%   The number of clusters C is given by C=2^(N_T-1), where N_T is the number of time points.
%
%   [INDEX,C] = TSD_clustering(X,BETA) finds C clusters in the data set X. X is size N_G-by-N_S, where N_G is 
%   the number of genes and N_S is the number of time points. INDEX is a structure (INDEX.CLUSTER) 
%   which contains the indices of X to form each cluster. BETA is a threashold for "change" filtering. A gene is 
%   filtered out if it does not reaches the threashold of change BETA in any of its transitions. If BETA=0 there
%   is no filtering.
%
%   Example
%        % Generate the data set X
%        X=rand(20,3);
%        % Set BETA and obtain the clusters
%        beta=0;
%        [index,c] = TSD_clustering(X,beta);
%        % Plot the resulting clusters in different figures
%        figure;
%        for i=1:c;
%            if ~isempty(index(i).cluster)
%               subplot(c,1,i);
%               plot(X(index(i).cluster,:)');
%            end;
%        end;

[n_g,n_t]=size(X);  %Get number of genes and number of samples

%%%%%%%%% Scoring the data, by states 1 if increases, 0 if decreases
xa=X;
for g=1:n_g;               % For all the genes
    for i=1:(n_t-1);       % For all the time points
        if ((xa(g,(i+1))>xa(g,(i)))&((abs((xa(g,(i+1))-xa(g,i))/xa(g,(i+1))))>beta));
            xt(g,i)=1;
        elseif ((xa(g,i)>xa(g,(i+1)))&((abs((xa(g,i)-xa(g,(i+1)))/xa(g,i)))>beta));
            xt(g,i)=0;
        else xt(g,i)=NaN;
        end;
    end;
end;

c=2^(n_t-1);              % Number of clusters or combinations

%%%%%%%%% Creating the combinations

a=2;                      % A will increase in multiples of 2
num=1;                    % Index of combinations 
while a<=c                % While A is less or equal number of possible combinations
    for i=1:c/a           % Creates a pattern for each A
        com_1(i)=1;
        com_0(i)=0;
    end;
    comb=[];
    for j=1:(a/2);        % Repeats the pattern to complete each state
        comb=[comb com_1 com_0];
    end;
    all_combinations(num,:)=comb(1,:);  % Save the complete series of patterns 
    a=a*2;
    num=num+1;                          % Next combination 
    clear com_1 com_0;
end
all_combinations=all_combinations';

%%%%%%%%% Creating the clusters

index(c).cluster=[];             % Creates the structure for the clusters
for i=1:c                        % For all the combinations
   a=1;                          % Index of a match 
   for j=1:n_g                   % For all the genes
        subtraction=all_combinations(i,:)-xt(j,:);
        if (isempty(find(subtraction==1)))&(isempty(find(subtraction==-1)))&(isempty(find(sum(isnan(subtraction)))));    %if the combination is equal to the gene
           index(i).cluster(a)=j;                % Index of genes for each combination
           a=a+1;                                % Increases for the next match in the same combination
        end;
        clear subtraction;
   end;
end;



