
function [V, U] = OLYMPUS(data, Nclusters, time_points, NP, vec_param)
    
    % This function implements OLYMPUS Clustering Tool for Short Time 
    % Series gene expression data
    
    % ------------------------------------------------------------------- %
    % CopyRight (c) 2012-2013 Aristidis G. Vrahatis
    %
    % Contact:  vrachatis at ceid dot upatras dot gr
    %           aris.vrahatis at gmail dot com
    % ------------------------------------------------------------------- %
    
    % ------------------------------------------------------------------- % 
    % INPUT PARAMETERS :
    % data          : Matrix(Rows: Genes, Cols: expression values of each 
    %                 gene across the time points of the experiment)
    % Nclusters     : User Defined Number of clusters 
    % time_points   : Horizontal Vector with the Time points of experiment
    % NP            : Size of Individuals [Note: Number of individuals must 
    %                 be > 5, the minimun number in the proposed mutation 
    %                 scheme DE\rand\2 and DE\best\2.
    % vec_param     : Vector [7x1] with parameters (real numbers)
    %                 [1] MaxIterOLYMPUS = Number of OLYMPUS iterations
    %                 [2] MaxIterDE      = Number of DE iterations
    %                 [3] MaxIterFSTS   = Number of FSTS Method iterations
    %                 [4] F e [0,1]      = DE constant (default value is 0.5)
    %                 [5] Cr e [0,1]     = DE constant (default value is 0.75)
    %                 [6] m              = fuzziness m>0 (default value is 1.4)
    %                 [7] strategy       = Number  e [1, 4] (Note:Mutation Scheme
    %                                                        [1]: DE-rand-1, 
    %                                                        [2]: DE-best-1, 
    %                                                        [3]: DE-rand-2, 
    %                                                        [4]: DE-best-2 
    % OUTPUT :
    %           V   : All the Cluster Centers in matrix [Nclusters x time_points]
    %           U   : Partition matrix [Nclusters x Number of genes]
    % ------------------------------------------------------------------- % 

% ------------------------------------------------------------------------ %
    % 1. Initialize parameters
    % Print error message if user gives wrong number of arguments
    if nargin < 4 || nargin > 5
       error('Too many or too few input arguments!');
    end
    % The vector vec_param may defined as default 
    if nargin < 5
        vec_param = defaultParameters();
    end
    
    MaxIterOLYMPUS  = vec_param(1);
    MaxIterDE       = vec_param(2);
    MaxIterFSTS    = vec_param(3);
    F               = vec_param(4);
    Cr              = vec_param(5);
    m               = vec_param(6);
    strategy        = vec_param(7);
    
    % ------------------------------------------------------------------- %
    
    % Data scaling in range [-1, 1]. 
    data = data-repmat(data(:,1),1,size(data,2));
    width = max(data') - min(data');
    data = data./repmat(width',1,size(data,2)); %X: xekinaei apo 0, times -1 ws 1
    % ------------------------------------------------------------------- %
     
    data=[zeros(size(data,1), 2), data];
    % Extra parameters
    dim = size(data,2);
    LB = min(min(data));
    UB = max(max(data));
    p_success = 0;
    
    % Initialize population of NP individuals with: 
    % a) Lower Bound: the minimum value of dataset and 
    % b) Upper Bound: the maximum value of dataset
    individuals = initPopulation(NP, dim, Nclusters, LB, UB);
    
    % Set iterative counts
    Iter1 = 0;
    best_idx = -1;
    
    % If Iter1 < MaxIterOLYMPUS continue, otherwise terminate
    while (Iter1 < MaxIterOLYMPUS)    
        % Set inner iterative counts
        Iter2 = 0;
        Iter3 = 0;
        
        % --------------------------------------------------------------- %
        % Differential Evolution (DE) method
        % --------------------------------------------------------------- %
        while (Iter2 < MaxIterDE)
             
            % Update the NP individuals according to equations 8, 11, 12 as  
            % stated in the published manuscript 
            %
            % Output--> particles: The new particles(if exist)
            %           impr: [1] if individual improve its position
            %                 [0] otherwise  
            [individuals impr] = updateIndividuals(individuals, Nclusters, time_points, F, Cr, m, data, LB, UB, strategy, best_idx);
            % Count the number of improvements of particles
            if impr == 1
                p_success = p_success + 1;
            end
            Iter2 = Iter2 + 1;
        end
        % Initialization with extreme values as a first attempt to find the 
        % best individuals that minimize the value of fitness function
        best = Inf;
        best_idx = -1;
        % For each individual calculate the fitness function and select the
        % best one
        for i=1:NP
            f = calculateFitnessSimple(Nclusters, time_points, individuals(i, :), data, m);
            if (f<best)
                best = f;
                best_idx = i;
            end
        end
        
        % --------------------------------------------------------------- %
        % Fuzzy  Clustering Method
        % --------------------------------------------------------------- %
        % Define each individual i as the initial cluster center vector in
        % Fuzzy Clustering Method
        % For each Individual process the Fuzzy Clustering Method for MaxIterFSTS times
        while (Iter3 < MaxIterFSTS)
            % Recalculate each cluster center using eq. 11
            % Calculate Partition matrix using eq. 12
            % Calculate Fitness Function using eq. 8  
            % Output: The position of all Individuals after DE and
            % Fuzzy Clustering Methods
            individuals(best_idx, :) = calculatePartition(Nclusters, time_points, individuals(best_idx, :), data, m);
            Iter3 = Iter3 + 1;
        end
        % Increase OLYMPUS iteration
        Iter1 = Iter1 + 1;
        % End of OLYMPUS iterations
    end
    
    % Print messages
    % Print the individual of population (index and value)
%     fprintf('Best individual index: %d,with fitness function value: %.8f\n', best_idx, best);
    % Print the number of DE succesful improvements during the iterations 
%     fprintf('DE succesful improvements: %d / %d\n', p_success, MaxIterOLYMPUS * MaxIterDE);
    % V : The final cluster center reshaped in a matrix [Nclusters x
    % time_points]     
    V = reshape(individuals(best_idx, :), dim, Nclusters)'; 
    % We delete the two first cols, which were added in the beginning for 
    % solving the equation system
    V = V(:, 3:size(V,2)); 
    % U: The final partition matrix [Nclusters x Ngenes]
    [~, U] = calculatePartition(Nclusters, time_points, individuals, data, m, best_idx);
    
    % End of ALGORITHM Implementation
end
% ------------------------------------------------------------------------ %
