function [newCenters retU] = calculatePartition(Nc, times, centers, data, m, single_center_index)
    
    % This function calculate and returns the partition matrix U(i, j)  
 
    % ------------------------------------------------------------------- %
    % Input  :  Nc      = Number of Clusters
    %           times   = Time points of experiment     
    %           centers = Cluster centers (V matrix)
    %           data    = Dataset
    %           m       = fuzziness of CLustering
    % Output :  Partition matrix U [Nc X Ngenes]
    %           (Each element of the partition matrix is setting the 
    %           membership degree of each gene to each cluster) 
    % ------------------------------------------------------------------- %
    
    
    % ------------------------------------------------------------------- %
    % Fill the distance matrix of c,data(cluster centers,dataset) with STS 
    % metric
    function [nc u] = new_center(c)
        dist = distance_matrix(c, times, data);
        
        function u12 = equation_12(i,j)
        % This function calculate the partition matrix U(i, j) as mentioned 
        % in Eq.12 of pseudocode
            s=0;
            exponent=1/(m-1);% 1/w-2
            %exponent=1;
            for q=1:Nc
                s=s+(dist(i,j)/dist(q,j))^exponent;
            end
            u12=1/s;
        end
        
        function mik_ret = equation_10(i,k)
        % This function calculate m(i, k) as mentioned in
        % Eq.10 of pseudocode
            global GLOBAL_d;
            global GLOBAL_e;
            global GLOBAL_f;
            if k == 1
                mik_ret=0;
                return
            end

            coeff=[GLOBAL_d(k) ; GLOBAL_e(k); GLOBAL_f(k)];   %3*1
            tempX=[data(:,k-1) , data(:,k), data(:,k+1)];    %Ng*3
            tempU=u(i, :);   %1*ng
            tempU=tempU.^m;
            mik_ret=tempU*(tempX*coeff);   % 1*ng x ng*1 = 1
            mik_ret=-mik_ret/sum(tempU);
        end
        
        function V = equation_11(i,n)
        % This function calculate the prototype V(i, n) as mentioned in
        % Eq.11 of pseudocode    
            global GLOBAL_a;
            global GLOBAL_c;            
            V=0;

            for r=2:(n-3)
                S=0;
                for pp=(r+3):n
                    S=S+prod(GLOBAL_c((pp-1):(n-1)))*prod(GLOBAL_a((r+1):(pp-2)));
                end

                V=V+ equation_10(i,r)*prod(GLOBAL_c(1:(r-1)))*(prod(GLOBAL_a((r+1):(n-1)))+prod(GLOBAL_c((r+1):(n-1)))+ S);
            end
            V=V+ equation_10(i,n-1)*prod(GLOBAL_c(1:(n-2)))+equation_10(i,n-2)*(GLOBAL_a(n-1)+GLOBAL_c(n-1))*prod(GLOBAL_c(1:(n-3)));
            V=V/prod(GLOBAL_c(2:(n-1)));
        end
        
        u = zeros(Nc, N);
        for i=1:Nc
            for j=1:N
                u(i,j)=equation_12(i, j);
            end
        end
        for j=1:N
            temp=sum(u(:,j));
            u(:,j)=u(:,j)./temp;
        end        
        V=zeros(size(c));
        for i=1:Nc
            for n=3:dim
                V(i,n)=equation_11(i,n);
            end
        end
        nc = V;
    end

    Nparticles = size(centers,1);
    N = size(data, 1);
    dim = size(data, 2);
    % Reshape the vector to matrix
    if nargin == 5
        for p=1:Nparticles
            centers(p, :) = reshape(        ...
                new_center (                ...
                    reshape(centers(p, :),  ...
                        dim,                ...
                        Nc)'         ...
                )',                         ...
                dim*Nc, 1)';
        end
        newCenters = centers;
        retU = 0;
    else
        [V U] = new_center (                ...
                    reshape(centers(single_center_index, :),  ...
                        dim,                ...
                        Nc)'         ...
                );
        V = reshape(V', dim*Nc, 1)';
        newCenters = V;
        retU = U;
    end   
    % ------------------------------------------------------------------- %

end