function [U, part_matrix, iter, center, eigvect, obj_fcn] = fuzzy_cv_clustering(X, n_c,w,n_i,e);

% FCV Data set clustering using fuzzy c-variaties clustering of linear prototypes.
% 
%   [U, PART_MATRIX, ITER, CENTER, EIGVECT, OBJ_FCN] = FUZZY_CV_CLUSTERING(X, N_C,W,N_I,E) 
%   finds N_C number of 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 samples for each gene. The coordinates for each cluster center are 
%   returned in the rows of the matrix CENTER. The membership function matrix U contains the grade of 
%   membership of each X point in each cluster. The values 0 and 1 indicate no membership and full 
%   membership respectively. Grades between 0 and 1 indicate that the data point has partial membership 
%   in a cluster. At each iteration, an objective function is minimized to find the best location for 
%   the clusters and its values are 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.
% 
%   Example
%       % Generate the data set X
%       x=1:0.25:5;
%       y=x*2-1;
%       x2=1:0.25:5;
%       y2=x2*1+1;
%       X=[x y;x2 y2]';
%       scatter(X(:,1),X(:,2));  
%       % Set the clustering paramenters and obtain partition matrix U
%       n_c=2; w=2.5; n_i=250; e=0.001;
%       [U,center,obj_fcn] = fuzzy_cv_clustering(X,n_c,w,n_i,e);
%       % Make the n_c clusters using a threashold ALPHA of membership degree.
%       alpha=0.8;
%       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;
%       scatter(X(index(1).cluster,1),X(index(1).cluster,2),'r');
%       hold on;
%       scatter(X(index(2).cluster,1),X(index(2).cluster,2),'g');
          

[n_g,n_s] = size(X);                   %Get number of genes and number of samples
iter = 0;                              %Iteration counter

U = initfcm(n_c, n_g);			       %Initial random fuzzy partition
Uo=zeros(n_c,n_g);

center=zeros(n_c,n_s);                 %Initialization of variables
D=zeros(n_c,n_g);
xv=zeros(n_g,n_s);
p=zeros(n_s,n_s);

%%%% Iterate
while  max(max(U-Uo)) > e    
    iter = iter + 1;
    part_matrix(iter)= max(max(U-Uo));  
    
    %%%% Calculate the centers
    mf = U;                                           %MF matrix before exponential modification 
    center = mf*X./((ones(n_s, 1)*sum(mf'))');        %New center
    
    %%%% Calculate the covariance matrices, eigenvalues and distances
    mf = U.^w;                                        %MF matrix after exponential modification
    sumf = sum(mf');
    
    %%%% Calculate the distance
    for j = 1 : n_c,                                    % For all clusters
        xv = X -ones(n_g,1)*center(j,:);                % All rows minus each centre
        p = ones(n_s,1)*mf(j,:).*xv'*xv/sumf(j);        % Calculate covariance matrix
        [ev,ed] = eig(p); ed = diag(ed)';               % Calculate eigen values and eigen vectors
        eigvect(:,j) = ev(:,ed == max(ed));             % Eigenvector corresponding to the r largest eigen value. r=1  
        D(j,:)=((diag((X- ones(n_g,1)*center(j,:))*(X - ones(n_g,1)*center(j,:))'))-((dot((X- ones(n_g,1)*center(j,:)),(ones(n_g,1)*eigvect(:,j)'),2)).^2))';       
     end;
    
    %%%% Calculate new U
    tmp = D  .^(-1/(w-1));  
    Uo=U;
    for i=1:n_c;
        for k=1:n_g;
            if  D(:,k)>0
                ratio_ds=(ones(n_c,1)*D(i,k))./D(:,k);
                exp_ratio=ratio_ds.^(2/(w-1));
                U(i,k)=1./sum(exp_ratio);
            end;
        end;
    end;
    
    %%%% Objective function
    obj_fcn(iter) = sum(sum((D.^2).*mf));
    
    if iter>n_i break; end;
end;   



