"""This program aims to cluster Flow Cytometry data using the previous work of A.Malkassian(2010-2012)
    This file computes the main Execution of the program.
    It is linked with the Controller in the Interface.

Marjorie Roscian CCi 2017-2018

Libraries needed :
- os
- matplotlib
- numpy
- pandas
- sklearn.linear_model
- scipy.interpolate
- math

Files needed :
- FunctionUtils.py
- AnthoRPhytoClus.py
- Your CSV file exported from the cyto software (English format with .)
"""
from Interface import Controller
from Analysis.AnthoRPhytoClus import *
from Analysis.FunctionUtils import *
from Analysis.DonneesTestBIS import *
from sklearn.decomposition import PCA

import os
from os.path import basename
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
rplot = robjects.r('plot')


def main_exe_4prof(filename):
    """ This function is the main action of the program.
    Modify it to change what Start Analysis button does.

    This version is for 4 profils : FWS,SWS,FLO,FLR"""

    # Set the working directory
    os.chdir(os.path.dirname(str(filename)))

    # If test needed, uncomment to get the test file:
    # data_matrix = donnees()

    # Open CSV file
    data = pd.DataFrame(pd.read_csv(basename(filename), sep=","))
    # Reshape file
    data_grouped = data.groupby("ID")
    data_transform = data_grouped.apply(np.array)
    data_transform = np.asarray(data_transform)

    data_matrix = []
    # Beug de la fonction np.copy si je modifie data_matrix raw_data change
    for cell in np.arange(0, len(data_transform)):
        data_temp = data_transform[cell][:, [1, 2, 4, 5]]
        data_matrix.append(data_temp)
    data_matrix = np.asarray(data_matrix)
    # raw_data = np.copy(data_matrix) This command does not work I don't know why. It is not a copy,
    #  it is a link to memomry place

    print("Number of cells : {}".format(str(len(data_matrix))))
    print("DONE open")

    # Global variable
    nb_col = 4

    # Scaling data
    for cell in np.arange(0, data_matrix.size):
        for profil in np.arange(0, nb_col):
            data_matrix[cell][:, profil] = data_matrix[cell][:, profil] / np.max(data_matrix[cell][:, profil])
    print("DONE scale")

    # Data to function
    nb_bases = 41
    final_data_coeffs = np.array(data_to_function(data_matrix, nb_bases, nb_col))

    # Split each profil
    new_profils = np.array(sort_each_profil(final_data_coeffs, nb_bases, nb_col))
    print('DONE split profils')

    # Sum of the square pair of bases/ Parseval's theorem (?)
    mod_coeffs = np.zeros((len(final_data_coeffs), nb_col * (int((nb_bases / 2)) + 1)))
    for cell in np.arange(0, len(final_data_coeffs)):
        nmod = 1
        for profil in np.arange(0, nb_col):
            for basis in np.arange(1, nb_bases, 2):
                mod_coeffs[cell, 0] = np.square(final_data_coeffs[cell][0, profil])
                mod_coeffs[cell, nmod] = np.square(new_profils[profil][cell, basis]) + np.square(
                    new_profils[profil][cell, basis + 1])
                nmod += 1
    print('DONE parseval')

    # Centering data
    for col in np.arange(0, nb_bases):
        mod_coeffs[:, col] = (mod_coeffs[:, col] - np.mean(mod_coeffs[:, col], axis=0))

    print('DONE center')

    # Scaling coefficients by first eigen value of the PCA
    pca = PCA()
    ncoef = 0

    while ncoef < nb_col * int(nb_bases / 2) + 1:
        var_temp = mod_coeffs[:, ncoef:ncoef + int(nb_bases / 2) + 1]
        pca.fit(var_temp)
        mod_coeffs[:, ncoef:ncoef + int(nb_bases / 2) + 1] = np.array(var_temp / np.sqrt(pca.singular_values_[0]))
        ncoef += 21
    print('DONE PCA scale')

    # Functional Principal Component Analysis
    cov_tot = np.cov(np.transpose(mod_coeffs))  # Caution: cov is row = variable, col = observations

    # Eigen solution
    e_values, e_vectors = np.linalg.eig(cov_tot)

    # Explained variance
    tot = sum(e_values)
    var_exp = [(i / tot) * 100 for i in sorted(e_values, reverse=True)]
    cum_var_exp = np.cumsum(var_exp)
    print(cum_var_exp)

    # Principal component (princ_comp), lambda (lamb)
    eig_vect = e_vectors
    lamb = np.diag(1 / np.sqrt(e_values))
    comp = mod_coeffs.dot(eig_vect)
    princ_comp = comp.dot(lamb)
    # Delete NaN
    princ_comp = np.real(np.nan_to_num(princ_comp))

    # Three first's principal components
    principal_component_array = np.column_stack((princ_comp[:, 0], princ_comp[:, 1], princ_comp[:, 2]))

    # Model-based clustering on 3 principal components
    utils.install_packages("mclust")
    mclust = importr("mclust")
    data_clust = mclust.mclustBIC(principal_component_array, G=np.arange(1, 16))
    rplot(data_clust)
    print("How many groups do you want ?")
    # Graphic version :
    # s = Controller.ControlMain().getint()
    # Console version :
    s = input()
    print("You typed", s)
    data_clust2 = mclust.Mclust(principal_component_array, G=s)
    data_clust2 = np.asarray(data_clust2.rx2('classification'))
    data_clust2 = np.transpose(np.atleast_2d(data_clust2))
    print('DONE MBC')

    # Re open data for raw_data
    data = pd.DataFrame(pd.read_csv(basename(filename), sep=","))
    # Reshape file
    data_grouped = data.groupby("ID")
    data_transform = data_grouped.apply(np.array)
    data_transform = np.asarray(data_transform)

    raw_data = []
    # Beug de la fonction np.copy si je modifie data_matrix raw_data change
    for cell in np.arange(0, len(data_transform)):
        data_temp = data_transform[cell][:, [1, 2, 4, 5]]
        raw_data.append(data_temp)
    raw_data = np.asarray(raw_data)

    # Cytogram plot with mclust's results

    tot_flr = np.zeros((raw_data.size, 1))
    tot_flo = np.zeros((raw_data.size, 1))
    tot_fws = np.zeros((raw_data.size, 1))
    max_sws = np.zeros((raw_data.size, 1))
    max_fws = np.zeros((raw_data.size, 1))
    for cell in np.arange(0, raw_data.size):
        tot_flo[cell] = area(raw_data[cell][:, 2], 500)
        tot_flr[cell] = area(raw_data[cell][:, 3], 500)
        tot_fws[cell] = area(raw_data[cell][:, 0], 500)
        max_fws[cell] = max(raw_data[cell][:, 0])
        max_sws[cell] = max(raw_data[cell][:, 1])

    # data_clust2 = np.ndarray.flatten(data_clust2)

    cytogram = np.column_stack((tot_fws, tot_flr, tot_flo, max_fws, max_sws, data_clust2, tot_flo / tot_flr))
    cytogram = pd.DataFrame(cytogram, columns=["fws", "flr", "flo", "maxfws", "maxsws", "group", "flo/flr"])

    sns.lmplot(x="flo", y="flr", data=cytogram, hue="group", fit_reg=False, palette="Set1", scatter_kws={"s": 1},
               legend=False).set(
        xscale="log", yscale="log")
    plt.title("Cytogram 1")
    plt.xlabel("tot_FLO")
    plt.ylabel("tot_FLR")

    sns.lmplot(x="fws", y="flo/flr", data=cytogram, hue="group", fit_reg=False, palette="Set1", scatter_kws={"s": 1},
               legend=False).set(
        xscale="log", yscale="log")
    plt.title("Cytogram 2")
    plt.xlabel("tot_FWS")
    plt.ylabel("tot_FLO/tot_FLR")

    sns.lmplot(x="maxfws", y="flr", data=cytogram, hue="group", fit_reg=False, palette="Set1", scatter_kws={"s": 1},
               legend=False).set(xscale="log", yscale="log")
    plt.title("Cytogram 3")
    plt.xlabel("max_FWS")
    plt.ylabel("tot_FLR")

    sns.lmplot(x="maxsws", y="flr", data=cytogram, hue="group", fit_reg=False, palette="Set1", scatter_kws={"s": 1},
               legend=False).set(xscale="log", yscale="log")
    plt.title("Cytogram 4")
    plt.xlabel("max_SWS")
    plt.ylabel("tot_FLR")

    plt.show()

    # Summary stats
    # Number of cells per group
    number_of_cell = cytogram.groupby("group").count()

    # Mean per group per profil
    mean_per_group = cytogram.groupby("group").mean()

    # STD per group per profil
    std_per_group = cytogram.groupby("group").std()

    # Final table to plot
    final_table = np.empty((int(s) + 1, 12), dtype=np.object)
    final_table[0, 0] = 'File name'
    final_table[0, 1] = 'Nb cells'
    final_table[0, 2] = 'N° group '
    final_table[0, 3] = 'Nb cells/group'
    final_table[0, 4] = 'Mean FWS'
    final_table[0, 5] = 'STD FWS '
    final_table[0, 6] = 'Mean SWS'
    final_table[0, 7] = 'SD SWS '
    final_table[0, 8] = 'Mean FLO'
    final_table[0, 9] = 'SD FLO'
    final_table[0, 10] = 'Mean FLR '
    final_table[0, 11] = 'SD FLR'
    final_table[1:, 0] = basename(filename)
    final_table[1:, 1] = str(raw_data.size)
    final_table[1:, 2] = np.arange(1, int(s) + 1)

    for group in np.arange(1, int(s) + 1):
        final_table[group, 3] = number_of_cell.iloc[group - 1, 0]

    for group in np.arange(1, int(s) + 1):
        profil = 0
        col = 4
        while col < 12 and profil < 7:
            final_table[group, col] = mean_per_group.iloc[group - 1, profil]
            final_table[group, col + 1] = std_per_group.iloc[group - 1, profil]
            profil += 1
            col += 2

    controller_table = Controller.ControlTable(final_table)
    controller_table.show()

    print(final_table)


def main_exe_5prof(filename):
    """ This function is the main action of the program.
    Modify it to change what Start Analysis button does.

    This version is for 5 profils : FWS,SWS,FLY,FLR,FLO"""

    # Set the working directory
    os.chdir(os.path.dirname(str(filename)))

    # If test needed, uncomment :
    # data_matrix = donnees()

    # Open CSV file
    data = pd.DataFrame(pd.read_csv(basename(filename), sep=","))

    # Reshape file
    data_grouped = data.groupby("ID")
    data_transform = data_grouped.apply(np.array)
    data_transform = np.asarray(data_transform)

    data_matrix = []
    for cell in np.arange(0, len(data_transform)):
        data_temp = data_transform[cell][:, [1, 2, 4, 5, 3]]
        data_matrix.append(data_temp)
    data_matrix = np.asarray(data_matrix)
    # raw_data = np.copy(data_matrix)
    print("Number of cells : {}".format(str(len(data_matrix))))
    print("DONE open")

    # Global variable
    nb_col = 5

    # Scaling data
    for cell in np.arange(0, data_matrix.size):
        for profil in np.arange(0, nb_col):
            data_matrix[cell][:, profil] = data_matrix[cell][:, profil] / np.max(data_matrix[cell][:, profil])
    print("DONE scale")

    # Data to function
    nb_bases = 41
    final_data_coeffs = np.array(data_to_function(data_matrix, nb_bases, nb_col))

    # Split each profil
    new_profils = np.array(sort_each_profil(final_data_coeffs, nb_bases, nb_col))
    print('DONE split profils')

    # Sum of the square pair of bases/ Parseval's theorem (?)
    mod_coeffs = np.zeros((len(final_data_coeffs), nb_col * (int((nb_bases / 2)) + 1)))
    for cell in np.arange(0, len(final_data_coeffs)):
        nmod = 1
        for profil in np.arange(0, nb_col):
            for basis in np.arange(1, nb_bases, 2):
                mod_coeffs[cell, 0] = np.square(final_data_coeffs[cell][0, profil])
                mod_coeffs[cell, nmod] = np.square(new_profils[profil][cell, basis]) + np.square(
                    new_profils[profil][cell, basis + 1])
                nmod += 1
    print('DONE parseval')

    # Centering data
    for col in np.arange(0, nb_bases):
        mod_coeffs[:, col] = (mod_coeffs[:, col] - np.mean(mod_coeffs[:, col], axis=0))
    print('DONE center')

    # Scaling coefficients by first eigen value of the PCA
    pca = PCA()
    ncoef = 0

    while ncoef < nb_col * int(nb_bases / 2) + 1:
        var_temp = mod_coeffs[:, ncoef:ncoef + int(nb_bases / 2) + 1]
        pca.fit(var_temp)
        mod_coeffs[:, ncoef:ncoef + int(nb_bases / 2) + 1] = np.array(var_temp / np.sqrt(pca.singular_values_[0]))
        ncoef += 21
    print('DONE PCA scale')

    # Functional Principal Component Analysis
    cov_tot = np.cov(np.transpose(mod_coeffs))  # Caution: cov is row = variable, col = observations

    # Eigen solution
    e_values, e_vectors = np.linalg.eig(cov_tot)

    # Explained variance
    tot = sum(e_values)
    var_exp = [(i / tot) * 100 for i in sorted(e_values, reverse=True)]
    cum_var_exp = np.cumsum(var_exp)
    print(cum_var_exp)

    # Principal component (princ_comp), lambda (lamb)
    eig_vect = e_vectors
    lamb = np.diag(1 / np.sqrt(e_values))
    comp = mod_coeffs.dot(eig_vect)
    princ_comp = comp.dot(lamb)
    # Delete NaN
    princ_comp = np.real(np.nan_to_num(princ_comp))

    # Three first's principal components
    principal_component_array = np.column_stack((princ_comp[:, 0], princ_comp[:, 1], princ_comp[:, 2]))

    # Model-based clustering on 3 principal components
    utils.install_packages("mclust")
    mclust = importr("mclust")
    data_clust = mclust.mclustBIC(principal_component_array, G=np.arange(1, 16))
    rplot(data_clust)
    print("How many groups do you want ?")
    s = Controller.ControlMain().getint()
    print("You typed", s)
    data_clust2 = mclust.Mclust(principal_component_array, G=s)
    data_clust2 = np.asarray(data_clust2.rx2('classification'))
    data_clust2 = np.transpose(np.atleast_2d(data_clust2))

    print('DONE MBC')

    # Re open data for raw_data
    data = pd.DataFrame(pd.read_csv(basename(filename), sep=","))
    # Reshape file
    data_grouped = data.groupby("ID")
    data_transform = data_grouped.apply(np.array)
    data_transform = np.asarray(data_transform)

    raw_data = []
    # Beug de la fonction np.copy si je modifie data_matrix raw_data change
    for cell in np.arange(0, len(data_transform)):
        data_temp = data_transform[cell][:, [1, 2, 4, 5]]
        raw_data.append(data_temp)
    raw_data = np.asarray(raw_data)

    # Cytogram plot with mclust's results

    tot_flr = np.zeros((raw_data.size, 1))
    tot_flo = np.zeros((raw_data.size, 1))
    tot_fws = np.zeros((raw_data.size, 1))
    max_sws = np.zeros((raw_data.size, 1))
    max_fws = np.zeros((raw_data.size, 1))
    for cell in np.arange(0, raw_data.size):
        tot_flo[cell] = area(raw_data[cell][:, 2], 500)
        tot_flr[cell] = area(raw_data[cell][:, 3], 500)
        tot_fws[cell] = area(raw_data[cell][:, 0], 500)
        max_fws[cell] = max(raw_data[cell][:, 0])
        max_sws[cell] = max(raw_data[cell][:, 1])

    # data_clust2 = np.ndarray.flatten(data_clust2)

    cytogram = np.column_stack((tot_fws, tot_flr, tot_flo, max_fws, max_sws, data_clust2, tot_flo / tot_flr))
    cytogram = pd.DataFrame(cytogram, columns=["fws", "flr", "flo", "maxfws", "maxsws", "group", "flo/flr"])

    sns.lmplot(x="flo", y="flr", data=cytogram, hue="group", fit_reg=False, palette="Set1", scatter_kws={"s": 1},
               legend=False).set(
        xscale="log", yscale="log")
    plt.title("Cytogram 1")
    plt.xlabel("tot_FLO")
    plt.ylabel("tot_FLR")

    sns.lmplot(x="fws", y="flo/flr", data=cytogram, hue="group", fit_reg=False, palette="Set1", scatter_kws={"s": 1},
               legend=False).set(
        xscale="log", yscale="log")
    plt.title("Cytogram 2")
    plt.xlabel("tot_FWS")
    plt.ylabel("tot_FLO/tot_FLR")

    sns.lmplot(x="maxfws", y="flr", data=cytogram, hue="group", fit_reg=False, palette="Set1", scatter_kws={"s": 1},
               legend=False).set(xscale="log", yscale="log")
    plt.title("Cytogram 3")
    plt.xlabel("max_FWS")
    plt.ylabel("tot_FLR")

    sns.lmplot(x="maxsws", y="flr", data=cytogram, hue="group", fit_reg=False, palette="Set1", scatter_kws={"s": 1},
               legend=False).set(xscale="log", yscale="log")
    plt.title("Cytogram 4")
    plt.xlabel("max_SWS")
    plt.ylabel("tot_FLR")

    plt.show()

    # Summary stats
    # Number of cells per group
    number_of_cell = cytogram.groupby("group").count()

    # Mean per group per profil
    mean_per_group = cytogram.groupby("group").mean()

    # STD per group per profil
    std_per_group = cytogram.groupby("group").std()

    # Final table to plot
    final_table = np.empty((int(s) + 1, 12), dtype=np.object)
    final_table[0, 0] = 'File name'
    final_table[0, 1] = 'Nb cells'
    final_table[0, 2] = 'N° group '
    final_table[0, 3] = 'Nb cells/group'
    final_table[0, 4] = 'Mean FWS'
    final_table[0, 5] = 'STD FWS '
    final_table[0, 6] = 'Mean SWS'
    final_table[0, 7] = 'SD SWS '
    final_table[0, 8] = 'Mean FLO'
    final_table[0, 9] = 'SD FLO'
    final_table[0, 10] = 'Mean FLR '
    final_table[0, 11] = 'SD FLR'
    final_table[1:, 0] = basename(filename)
    final_table[1:, 1] = str(raw_data.size)
    final_table[1:, 2] = np.arange(1, int(s) + 1)

    for group in np.arange(1, int(s)+1):
        final_table[group, 3] = number_of_cell.iloc[group-1, 0]

    for group in np.arange(1, int(s)+1):
        profil = 0
        col = 4
        while col < 12 and profil < 7:
            final_table[group, col] = str(mean_per_group.iloc[group-1, profil])
            final_table[group, col+1] = str(std_per_group.iloc[group-1, profil])
            profil += 1
            col += 2

    controller_table = Controller.ControlTable(final_table)
    controller_table.show()

    print(final_table)
