% This Script runs the OLYMPUS Clustering Tool 
    

% ----------------------------------------------------------------------- % 
% INPUT PARAMETERS :
% data          : Matrix(Rows: Genes, Cols: expression values of each 
%                 gene across the time points of the experiment). In this
%                 example we get the dataset from our synthetic data as
%                 their referred in the paper)          
% Nclusters     : User Defined Optimal 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.
% rest          : Vector [1x7] 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 ... (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 :  Print a message with the optimal number of clusters and after
%           running the OLYMPUS with this 'optimal' number get V,U matrices
%           V : All the Cluster Centers in matrix [Nclusters x time_points]
%           U : Partition matrix [Nclusters x Number of genes]
% ----------------------------------------------------------------------- % 
% ----------------------------------------------------------------------- %
% CopyRight (c) 2012-2013 Aristidis G. Vrahatis
%
% Contact:  vrachatis at ceid dot upatras dot gr
%           aris.vrahatis at gmail dot com
% ----------------------------------------------------------------------- %
% ----------------------------------------------------------------------- % 
% ------------------------- EXAMPLE SYNTHETIC DATASET ------------------- %    
% ------------------------- REPLACE WITH YOUR DATA----------------------- % 
% Load Data Set
    % Synthetic Data construction as referred in the paper (In this example
    % we test 1000 expression profiles across 15 time points (dimension size=10)
    % and the number of optimal clusters is 5)
    data = test_sts_data(200, 5, 15);
    % Set the range of number of cluster. The minimum number is 2 and the maximum is 
    % square(X), X is the number of rows. 
    Nc_min  = 2;
    Nc_step = 1;
    Nc_max  = floor(sqrt(size(data,1)));
    % Initialize of sampling step after OLYMPUS approach a canditate
    % cluster number solution (global minimum of fittness function)
    Nc_step_sampling = 5;
    % ------------------------- SET PARAMETERS--------------------------- %     
    % Input Parameters
    % Number of individuals must be > 5, the minimun number when we use mutation scheme
    % DE\rand\2 and DE\best\2.
    NP = 5;
    % Time-points of the experiment(Dimensions) 
    ts = 0:size(data,2)-1;
    % Initialize of 'rest' vector     
    % The vector 'rest' take values : 
    % a) OLYMPUS Iterations                 (suggested is 10) 
    % b) DE Iterations                      (suggested is 5)
    % c) Fuzzy Method Iterations            (suggested is 5)
    % b) F constant take values 0.1 up 0.9  (suggested is 0.5)
    % c) Cr take values 0.1 up 0.9          (suggested is 0.75)
    % d) fuzziness take values 0.1 up INF   (suggested is 1.4)
    % e) Mutation Scheme suggested is [1] (DE\rand\1\bin)
    rest = [1;20;40;0.5;0.75;1.4;1];
% ----------------------------------------------------------------------- % 
% -------------------- INITIALIZATION OF OLYMPUS ------------------------ % 

    fprintf('-----------------------------------------------------------------------------\n');
    fprintf('---------------------- OLYMPUS Clustering Tool ------------------------------\n');
    fprintf('-----------------------------------------------------------------------------\n');
    % Add values -1 and 0 in the first and second column of ts respectively.
    % These values used for solving the equation system
    ts= [-1, 0, ts+1];
    % Data Dimensions
    [N Nt] = size(data);

    % ----------------------------------------------------------------------- % 
    % calculate a,b,c,d,e,f  as mentioned in equations 8, 9, 10, 11 
    global GLOBAL_a;
    GLOBAL_a = zeros(Nt+1,1);
    global GLOBAL_b;
    GLOBAL_b = zeros(Nt+1,1);
    global GLOBAL_c;
    GLOBAL_c = zeros(Nt+1,1);

    GLOBAL_a(1)=-(ts(2)-ts(1))^2;
    % Default Value
    GLOBAL_c(1)=1;

    for k=2:Nt+1
        GLOBAL_a(k)=-(ts(k+1)-ts(k))^2;
        GLOBAL_c(k)=-(ts(k)-ts(k-1))^2;
    end
    GLOBAL_b = -(GLOBAL_a + GLOBAL_c);
    global GLOBAL_d;
    GLOBAL_d=-GLOBAL_a;
    global GLOBAL_f;
    GLOBAL_f=-GLOBAL_c;
    global GLOBAL_e;
    GLOBAL_e=GLOBAL_a+GLOBAL_c;
    % ---------------------------------------------------------------------
    best_Ncluster = Nc_min;
    best_BIC = Inf;
    no_improve_BIC = 0;
    % If BIC score not improved after 'no_improve_BIC_max' clusters, stop the 
    % programme. 
    no_improve_BIC_max = 3;
    no_improve_BIC_sample = 1;

    % Call the main function of our Clustering Tool
    for Nclusters = Nc_min:Nc_step:Nc_max,
        fprintf('-----------------------------------------------------\n');
        fprintf('Calculating partitioning with %d clusters.\n', Nclusters);
        [V U] = OLYMPUS(data, Nclusters, ts, NP, rest);
        BIC = BIC_count(V,U,ts,data,N,Nclusters);
        if (BIC<best_BIC)
            best_Ncluster = Nclusters;
            best_BIC = BIC;
            no_improve_BIC = 0;
        else
        no_improve_BIC = no_improve_BIC + 1;
        if (no_improve_BIC >= no_improve_BIC_max)
            fprintf('Identified a canditate solution and sampling to validate our result\n');
            for i_sample = Nclusters+Nc_step_sampling:Nc_step_sampling:Nc_max
                fprintf('-----------------------------------------------------\n');
                fprintf('Calculating partitioning with %d clusters.\n', i_sample);
                 [V U] = OLYMPUS(data, i_sample, ts, NP, rest);
                 BIC = BIC_count(V,U,ts,data,N,i_sample);
                 sample_improved = 0;
                 if (BIC<best_BIC)
                    best_Ncluster = i_sample;
                    best_BIC = BIC;
                    sample_improved = sample_improved+1;
                 end
            end
          break;
        end
        end
    end
    fprintf('Best clustering achieved with: %d clusters, with BIC score: %d \n', best_Ncluster, best_BIC);

    % Run OLYMPUS with the optimal number of cluster
    [V U] = OLYMPUS(data, best_Ncluster, ts, NP, rest);
    % Plotting
    Plot_CLusters
    % ----------------------------------------------------------------------- % 
