%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% 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) = [];