%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                    %
%     Fuzzy Short Time-Series (FSTS) Clustering      %
%                                                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%FSTS data set clustering.
% 
%   [U, PROTOTYPE, OBJ_FCN, PART_MATRIX] = sts_algorithm(X, T, N_C, W, N_I,E)
%   finds N_C number of clusters in the data set X. X is size N_G-by-N_T, where N_G is the number of genes 
%   and N_T is the number of time points for each gene, T is the time vector. 
%   The membership function matrix U contains the grade of membership of each time-series in each cluster. 
%   The values 0 and 1 indicate no membership and full membership respectively. Grades between 0 and 1 indicate 
%   that the time-series has partial membership in a cluster. The coordinates for each cluster prototype are 
%   returned in the rows of the matrix PROTOTYPE. At each iteration, an objective function is minimized to find the best location for 
%   the clusters and the objective function vector is returned in OBJ_FCN. 
%   The iteration finishes when the number of predefined iterations N_I is reached or when the partition matrix PART_MATRIX has a minimum change 
%   of E from one iteration to the next one. 
%   NOTE: This function utilizes the function "initfcm.m" available from
%   the fuzzy matlab toolbox. This function generates a random partition of
%   the data. Alternatively, the function "fixed_inital_centres.m" can be
%   used to obtain fixed centres (using hierarchical clustering results).
% 
%   Example
%   % Generate the data set X
%       X=[3.5 3 4 4; 2 2 3 3; 1 1 2 1.5];
%       t=[1 2 3 10];
%       % Scale X and t
%       [X, t, X_original, t_original] = sts_scale(X,t);
%       % Set the clustering paramenters and obtain partition matrix U
%       n_c=2; w=1.6; n_i=250; e=0.001; 
%       [prototype, U, obj_fcn, part_matrix] = sts_algorithm(X, t, n_c, w,n_i,e);
%       % Make the n_c clusters using a threashold ALPHA of membership degree.
%       alpha=0.4;
%       for i=1:n_c,
%           indices = find(U(i,:) > alpha);
%           if ~isempty(indices),
%              index(i).cluster=indices;
%           end; 
%       end;
%       % Plot the resulting clusters with different colours
%       figure;
%       plot(t_original,X_original(index(1).cluster,:),'r');
%       hold on;
%       plot(t_original,X_original(index(2).cluster,:),'g');

function [prototype, U, obj_fcn, part_matrix] = sts_algorithm(X, t, n_c,w,n_i,e)

n_g = size(X, 1);
n_t = size(X, 2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% The constants
for k=2:n_t-1
    a(k)=-(t(k+1)-t(k)).^2;
    c(k)=-(t(k)-t(k-1)).^2;
    d(k)=+(t(k+1)-t(k)).^2;
    e(k)=-((t(k+1)-t(k)).^2)-((t(k)-t(k-1)).^2);
    f(k)=+(t(k)-t(k-1)).^2;
end;
a(1)=1;
c(1)=1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% The algorithm 

max_iter = n_i;		% Max. iteration
min_impro = e;		% Min. improvement

obj_fcn = zeros(max_iter, 1);	    % Initialization of objective function array
Uo = initfcm(n_c,n_g);              % Initial fuzzy partition
U = Uo;

% Main loop
for loop = 1:max_iter,  
    part_matrix(loop)= max(max(U-Uo));           % Change of the partition matrix
    prototype =sts_prototype(U,w,X,a,c,d,e,f);      % New prototype
    dist = sts_distance(prototype,X,t);             % Calculate the distance matrix
    % Calculate new memebership degrees
    Uo=U; 
    for i=1:n_c;
        for j=1:n_g;
            if  dist(:,j)>0
                ratio_ds=(ones(n_c,1)*dist(i,j))./dist(:,j);
                exp_ratio=ratio_ds.^(1/(w-1));
                U(i,j)=1./sum(exp_ratio);
            end;
        end;
    end;
    % Calculate objective function
    mf = U.^w;      
    obj_fcn(loop) = sum(sum((dist.^1).*mf));  
    % Check termination condition
    if loop > 1,
         if  abs(obj_fcn(loop) - obj_fcn(loop-1)) < min_impro, break; end;        
    end 
end
% Final objective function vector
part_matrix(1)=[];
iter_n = loop;	
obj_fcn(iter_n+1:max_iter) = [];