/*

                                                         ##
                                                         ##
                         #########                       ##
          ##########    ##########    ###        ###     ##
          ##########    ##            ####      ####     ##      ##########
          ##            ##            ##  ##  ##  ##     ##      ##########
          ##            ########      ##  ##  ##  ##     ##      ##
          ######         ########     ##    ##    ##     ##      ##
          ######               ###    ##    ##    ##     ##############
          ##     ###     ##    ##########         ##     ##############
          ##     ###    #################        ########        ##
          ###########   #########     ##         #########       ##
          ###########    ##    ##     ##         ###    ##       ##########
                 ##  ##  ##    ######            ##     ##       ##########
                 ##  ##  ##    ######            ########
                 ##    # ##    ##                #########
                 ##    # ##    ##                ##      ##
                 ##     ###    ##########        ##      ##
                 ##     ###    ##########        ##########
                                                 #########
                                                                                   */

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
**                                                                                 **
**                                                                                 **
**                                 E N S E M B L E                                 **
**                                                                                 **
**                                   Version 2.1                                   **
**                                                                                 **
**                                                                                 **
**                   Principal Investigator: Julie D. Forman-Kay                   **
**                                                                                 **
**                           Author:  Mickaël Krzeminski                           **
**                                                                                 **
**                               Date: November 2012                               **
**                                                                                 **
**                                                                                 **
*************************************************************************************
**                                                                                 **
**                                                                                 **
**   Copyright (C) The Hospital for Sick Children, 2001                            **
**                                                                                 **
**   Distribution of substantively modified versions of this module is prohibited  **
**  without the explicit permission of the copyright holder.                       **
**                                                                                 **
**   Any use of this work or derivative works in whole or in part for any          **
**  commercial purpose or for monetary gain is prohibited.                         **
**                                                                                 **
**                                                                                 **
**                                   NO WARRANTY                                   **
**   This software package is provided 'as is' without warranty of any kind,       **
**  expressed or implied.                                                          **
**                                                                                 **
**                                                                                 **
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
\* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

/*
   COMPILATION
                gcc -O3 -lm -o analyzens analyzens.c
*/

#include "Include/ens.h"


enum toolsPrograms TP = __ANALYZE__;

void _cleanAnalyze(PDBListData *ENSData, residueList *residueNbs, secondaryStructure *statSecStr, confList *confNbsList, float *atomMasses, flaver **distances, char *ensFileName, char *dirAnalyze, char *dirConformers, char *dirCurrent, char *sequence)
{
    counter count;

    if (ENSData != NULL) _cleanPDB(ENSData);

    if (residueNbs != NULL)
    {
        if (distances != NULL)
        {
            for (count=0; count<(*residueNbs).nb; count++)
            {
                free(distances[count]); distances[count] = NULL;
            }
        }

        free((*residueNbs).nbs); (*residueNbs).nbs = NULL;
    }

    if (statSecStr != NULL)
    {
        free((*statSecStr).Helix);    (*statSecStr).Helix     = NULL;
        free((*statSecStr).Pihelix);  (*statSecStr).Pihelix   = NULL;
        free((*statSecStr).Helix310); (*statSecStr).Helix310  = NULL;
        free((*statSecStr).Strand);   (*statSecStr).Strand    = NULL;
        free((*statSecStr).Coil);     (*statSecStr).Coil      = NULL;
        free((*statSecStr).Turn);     (*statSecStr).Turn      = NULL;
        free((*statSecStr).Bridge);   (*statSecStr).Bridge    = NULL;
    }

    free(distances); distances = NULL;
    free(atomMasses); atomMasses = NULL;

    if (confNbsList != NULL)
    {
        free((*confNbsList).confNb);
        (*confNbsList).confNb = NULL;
    }

    ensFileName = NULL;
    free(dirAnalyze); dirAnalyze = NULL;
    free(dirConformers); dirConformers = NULL;
    free(dirCurrent); dirCurrent = NULL;
    free(sequence); sequence = NULL;
}


int main(int argc, char *argv[])
{
    /***********************************************************\
    |*                                                         *|
    |*                                                         *|
    |*                  Variables declaration                  *|
    |*                                                         *|
    |*                                                         *|
    \***********************************************************/
    counter count, countx, county, confNb, resiNb;

    char *ensFileName = NULL;
    char *dirAnalyze = NULL;
    char *dirConformers = NULL;
    char *dirCurrent = getcwd(NULL, 0);
    char *sequence = NULL;
    char fileType[4] = {"\0"};

    float *atomMasses = NULL;
    confList confNbsList = {0, NULL};

    mode_t processMask;
    residueList residueNbs = {0, NULL};
    PDBListData ENSData = {0, 0, NULL, NULL};
    secondaryStructure statSecStr = {NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0};
    flaver **distances = NULL;
    color ldRGB = {0., 0., 0.}, udRGB = {1., 1., 1.};

    FILE *ensFile, *logFile, *resFile;

    /***********************************************************\
    |*                                                         *|
    |*                                                         *|
    |*                       Preparation                       *|
    |*                                                         *|
    |*                                                         *|
    \***********************************************************/
    if (argc != 2)
    {
        __USAGE__
    }

    printf("Checking arguments and creating the analysis directory... ");

    // Instanciation of the provided ENSEMBLE file name
    ensFileName = &argv[1][0];

    // Opening the provided file in binary mode and checking its existence
    if ((ensFile = fopen(ensFileName, "rb")) == NULL) fatal_error("NO_FILE", ensFileName);

    // Verifying that it is a ENS file
    fread(&fileType, CHAR_SIZE, 3, ensFile);
    fclose(ensFile);

    if (strcmp(fileType, "Ens")) fatal_error("ENS_FILE", ensFileName);

    // Creating the Analysis directory
    dirAnalyze = malloc((strlen(ensFileName)+1)*CHAR_SIZE);
    strcpy(dirAnalyze, ensFileName);

    for (count=strlen(ensFileName)-1; count>=0; count--)
    {
        if (ensFileName[count] == '/') break;
        dirAnalyze[count] = toupper(ensFileName[count]);
    }

    if ((dirAnalyze = realloc(dirAnalyze, (strlen(dirAnalyze)+9)*CHAR_SIZE)) == NULL)
    {
        _cleanAnalyze(&ENSData, &residueNbs, &statSecStr, &confNbsList, atomMasses, distances, ensFileName, dirAnalyze, dirConformers, dirCurrent, sequence);
        fatal_error("The program encountered a critical problem while allocating memory for the analyzing directory name (%s)", dirAnalyze);
    }
    strcat(dirAnalyze, "_ANALYZE");

    processMask = umask(0);
    mkdir(dirAnalyze, S_IRWXU);

    if (errno==EACCES)
    {
        printf ("\n\tCannot create the directory %s (Write permission is denied)!\n\n", dirAnalyze);

        _cleanAnalyze(&ENSData, &residueNbs, &statSecStr, &confNbsList, atomMasses, distances, ensFileName, dirAnalyze, dirConformers, dirCurrent, sequence);
        exit(EACCES);
    }
    else if (errno==EEXIST)
    {
        printf ("\n\tThe directory %s already exists!\n\n", dirAnalyze);

        _cleanAnalyze(&ENSData, &residueNbs, &statSecStr, &confNbsList, atomMasses, distances, ensFileName, dirAnalyze, dirConformers, dirCurrent, sequence);
        exit(EEXIST);
    }
    else if (errno==EROFS)
    {
        printf ("\n\tThe parent directory of %s is on a read-only file system and cannot be modified!\n\n", dirAnalyze);

        _cleanAnalyze(&ENSData, &residueNbs, &statSecStr, &confNbsList, atomMasses, distances, ensFileName, dirAnalyze, dirConformers, dirCurrent, sequence);
        exit(EROFS);
    }

    umask(processMask);

    printf("DONE\n");

    /***********************************************************\
    |*                                                         *|
    |*                                                         *|
    |*                     Recording  Data                     *|
    |*                                                         *|
    |*                                                         *|
    \***********************************************************/
    // Recording the conformer numbers
    printf ("Recording data... ");

    if (getConfListNb(&confNbsList, ensFileName))
    {
        _cleanAnalyze(&ENSData, &residueNbs, &statSecStr, &confNbsList, atomMasses, distances, ensFileName, dirAnalyze, dirConformers, dirCurrent, sequence);
        fatal_error("The analyzing script encountered a problem.", NULL);
    }

    // Getting all info about conformers
    if (getCoords(&ENSData, ensFileName, NULL))
    {
        _cleanAnalyze(&ENSData, &residueNbs, &statSecStr, &confNbsList, atomMasses, distances, ensFileName, dirAnalyze, dirConformers, dirCurrent, sequence);
        fatal_error("The program encountered a critical problem while recording coordinates from the provided list file.", NULL);
    }

    if (getMasses(&atomMasses, &ENSData))
    {
        _cleanAnalyze(&ENSData, &residueNbs, &statSecStr, &confNbsList, atomMasses, distances, ensFileName, dirAnalyze, dirConformers, dirCurrent, sequence);
        fatal_error("The program encountered a critical problem while recording masses from provided list file.", NULL);
    }

    if (getResidueNumbers(&residueNbs, &ENSData))
    {
        _cleanAnalyze(&ENSData, &residueNbs, &statSecStr, &confNbsList, atomMasses, distances, ensFileName, dirAnalyze, dirConformers, dirCurrent, sequence);
        fatal_error("The program encountered a critical problem while recording residue numbers from provided list file.", NULL);
    }
    if (!residueNbs.nb)
    {
        _cleanAnalyze(&ENSData, &residueNbs, &statSecStr, &confNbsList, atomMasses, distances, ensFileName, dirAnalyze, dirConformers, dirCurrent, sequence);
        fatal_error("The program encountered a critical problem while recording residue numbers from atom headers.", NULL);
    }

    printf ("DONE\n");

    /***********************************************************\
    |*                                                         *|
    |*                                                         *|
    |*                        Analysis                         *|
    |*                                                         *|
    |*                                                         *|
    \***********************************************************/
    if (chdir(dirAnalyze))
    {
        fatal_error("Impossible to work in the analyzing directory: %s", dirAnalyze);        // Go to the analyzing directory

        _cleanAnalyze(&ENSData, &residueNbs, &statSecStr, &confNbsList, atomMasses, distances, ensFileName, dirAnalyze, dirConformers, dirCurrent, sequence);
	return 1;
    }

    /***********************************************************\
    |*                       Information                       *|
    \***********************************************************/
    // Writing in the log file all important data
    if ((logFile = fopen("analyze.log", "wt")) == NULL)
    {
        error("OPEN_FILE", strcat(dirAnalyze, "/analyze.log"));
        _cleanAnalyze(&ENSData, &residueNbs, &statSecStr, &confNbsList, atomMasses, distances, ensFileName, dirAnalyze, dirConformers, dirCurrent, sequence);
        return 1;
    }

    fprintf(logFile, "================= General Informations =================\n");

    pdb2seq(&sequence, &ENSData);
    fprintf(logFile, "Sequence: %s\n", sequence);
    fprintf(logFile, "========================================================\n\n\n");

    /***********************************************************\
    |*                  Extracting structures                  *|
    \***********************************************************/
    // Creating the Conformer directory (Where structures are stored)
    printf("Extracting structures... ");
    fprintf(logFile, "================= Conformers extraction ================\n");

    if ((dirConformers = strdup("Conformers")) == NULL)
    {
        fprintf(logFile, "Memory allocation occured while attempting to create the 'Conformers' directory!\nThe analysis has aborted!\n");
        fclose(logFile);

        _cleanAnalyze(&ENSData, &residueNbs, &statSecStr, &confNbsList, atomMasses, distances, ensFileName, dirAnalyze, dirConformers, dirCurrent, sequence);
        fatal_error("MEMORY", "create the 'Conformers' directory");
    }

    processMask = umask(0);
    mkdir(dirConformers, S_IRWXU);        // No need any test, as the parent directory has already been tested
    umask(processMask);

    if (extractCSP(&ENSData, &confNbsList, dirConformers, "conf_", 1))
    {
        fprintf(logFile, "Conformers extraction failed!\n");
        fprintf(logFile, "========================================================\n\n\n");
    }
    else
    {
        fprintf(logFile, "Conformers extracted successfully in the Conformers directory.\n");
    }
    fprintf(logFile, "========================================================\n\n\n");

    printf("DONE\n");

    /***********************************************************\
    |*                   Radius of Gyration                    *|
    \***********************************************************/
    printf("Calculating radia of gyration... ");
    radiaOfGyration rog = getRoG(&ENSData, &atomMasses);

    if ((resFile = fopen("radiusGyration.txt", "wt")) == NULL)
    {
        fprintf(logFile, "Impossible to access the radiusGyration.txt file!\nThe analysis has aboarted!\n");
        fclose(logFile);

        error("OPEN_FILE", strcat(dirAnalyze, "radiusGyration.txt"));
        _cleanAnalyze(&ENSData, &residueNbs, &statSecStr, &confNbsList, atomMasses, distances, ensFileName, dirAnalyze, dirConformers, dirCurrent, sequence);
        return 1;
    }

    for (confNb=0; confNb<rog.nb; confNb++) fprintf(resFile, "Conformer %u: %.2f\n", confNb+1, rog.rg[confNb]);
    fprintf(resFile, "========================================================\n");
    fprintf(resFile, "Average radius of gyration: %.2f +/- %.2f A\n", rog.mean.average, rog.mean.stdv);
    fprintf(resFile, "========================================================\n");

    fclose(resFile);

    fprintf(logFile, "================== Radius of gyration ==================\n");
    fprintf(logFile, "Radius of gyration calculated successfully.\n");
    fprintf(logFile, "See the file radiusGyration.txt for details.\n");
    fprintf(logFile, "========================================================\n\n\n");

    // Releasing memory
    free(rog.rg); rog.rg = NULL;

    printf("DONE\n");

    /***********************************************************\
    |*                  Radius of Hydration                    *|
    \***********************************************************/


    /***********************************************************\
    |*                  Secondary Structures                   *|
    \***********************************************************/
    printf("Calculating secondary structures... ");
    fprintf(logFile, "========== Secondary Structures Distribution ==========\n");

    statSecStr = getSecondaryStructure(&ENSData);
    if (!statSecStr.nbResidues)
    {
        fprintf(logFile, "Impossible to retrieve the secondary structures of conformers.\n");\
        fprintf(logFile, "Secondary Structure calculations failed!\n");
        fprintf(logFile, "========================================================\n");
        fclose(logFile);

        _cleanAnalyze(&ENSData, &residueNbs, &statSecStr, &confNbsList, atomMasses, distances, ensFileName, dirAnalyze, dirConformers, dirCurrent, sequence);
        fatal_error("The program encountered a critical problem while attempting to get the secondary structure distribution of the ensemble!", NULL);
    }

    if ((resFile = fopen("ss_distr.txt", "wt")) == NULL)
    {
        error("OPEN_FILE", strcat(dirAnalyze, "ss_distr.txt"));

        _cleanAnalyze(&ENSData, &residueNbs, &statSecStr, &confNbsList, atomMasses, distances, ensFileName, dirAnalyze, dirConformers, dirCurrent, sequence);
        return 1;        
    }    
    fprintf(resFile, "=========== Secondary Structure Distribution ===========\n");

    fprintf(resFile, "Residue    Helix     Pi_Helix  Helix3_10 Strand    Coil      Turn      Bridge\n");
    for (resiNb=0; resiNb<residueNbs.nb; resiNb++)
    {
        fprintf(resFile, "%-7d    %-6.2f    %-6.2f    %-6.2f    %-6.2f    %-6.2f    %-6.2f    %-6.2f\n", residueNbs.nbs[resiNb],\
                                                                               (float)statSecStr.Helix[resiNb]/ENSData.nbConfs,\
                                                                               (float)statSecStr.Pihelix[resiNb]/ENSData.nbConfs,\
                                                                               (float)statSecStr.Helix310[resiNb]/ENSData.nbConfs,\
                                                                               (float)statSecStr.Strand[resiNb]/ENSData.nbConfs,\
                                                                               (float)statSecStr.Coil[resiNb]/ENSData.nbConfs,\
                                                                               (float)statSecStr.Turn[resiNb]/ENSData.nbConfs,\
                                                                               (float)statSecStr.Bridge[resiNb]/ENSData.nbConfs);
    }

    fprintf (resFile, "\n\n\nResidue        Helix          Strand         Coil&Turn      Bridge\n");
    for (resiNb=0; resiNb<statSecStr.nbResidues; resiNb++)
    {
        fprintf(resFile, "%-7d        %-6.2f         %-6.2f         %-6.2f         %-6.2f\n", residueNbs.nbs[resiNb],\
                                                    (float)(statSecStr.Helix[resiNb]+statSecStr.Pihelix[resiNb]+statSecStr.Helix310[resiNb])/ENSData.nbConfs,\
                                                    (float)statSecStr.Strand[resiNb]/ENSData.nbConfs,\
                                                    (float)(statSecStr.Coil[resiNb]+statSecStr.Turn[resiNb])/ENSData.nbConfs,\
                                                    (float)statSecStr.Bridge[resiNb]/ENSData.nbConfs);
    }

    fclose(resFile);

    fprintf(logFile, "Secondary structures estimated successfully.\n");
    fprintf(logFile, "========================================================\n\n\n");

    printf("DONE\n");

    /***********************************************************\
    |*                 Ca-Ca distance mapping                  *|
    \***********************************************************/
    printf("Building the Ca-Ca distance map... ");

    // This MUST be the last part of the analysis, as only Ca atoms will be kept
    fprintf(logFile, "================ Ca-Ca distance  matrix ================\n");

    // Recording the CA coordinates of all atoms of the PDB list file
    if (getAtomCoords(&ENSData, NULL, "CA"))
    {
        fprintf(logFile, "Impossible to retrieve the Ca atom coordinates.\n");
        fprintf(logFile, "Ca-Ca distance matrix calculations failed!\n");
        fprintf(logFile, "========================================================\n\n\n");
        fclose(logFile);

        _cleanAnalyze(&ENSData, &residueNbs, &statSecStr, &confNbsList, atomMasses, distances, ensFileName, dirAnalyze, dirConformers, dirCurrent, sequence);
        fatal_error("A problem occured while selecting Ca atoms", NULL);
    }

    if (caca_map(&distances, &ENSData, &residueNbs, "./caca_map.eps", "./caca_map_sdtv.eps", 3.5, 50., 50., &ldRGB, &udRGB))
    {
        fprintf(logFile, "Impossible to build the Ca-Ca matrix.\n");
        fprintf(logFile, "Ca-Ca distance matrix calculations failed!\n");
        fprintf(logFile, "========================================================\n\n\n");
        fclose(logFile);

        _cleanAnalyze(&ENSData, &residueNbs, &statSecStr, &confNbsList, atomMasses, distances, ensFileName, dirAnalyze, dirConformers, dirCurrent, sequence);
        fatal_error("A problem occured while selecting Ca atoms", NULL);
    }

    if ((resFile = fopen("caca_map.txt", "wt")) == NULL)
    {
        fprintf(logFile, "Impossible to access the caca_map.txt file.\n");
        fprintf(logFile, "Ca-Ca distance matrix calculations failed!\n");
        fprintf(logFile, "========================================================\n\n\n");
        fclose(logFile);

        error("OPEN_FILE", strcat(dirAnalyze, "caca_map.txt"));
        _cleanAnalyze(&ENSData, &residueNbs, &statSecStr, &confNbsList, atomMasses, distances, ensFileName, dirAnalyze, dirConformers, dirCurrent, sequence);
        return 1;
    }

    fprintf(resFile, "================ Ca-Ca distance  matrix ================\n");
    for (countx=0; countx<residueNbs.nb-1; countx++)
    {
        for (county=countx+1; county<residueNbs.nb; county++)
        {
            fprintf(resFile, "Ca_%d-Ca_%d:\t%6.2f +/- %6.2f\n", residueNbs.nbs[countx], residueNbs.nbs[county], distances[countx][county].average, distances[countx][county].stdv);
        }
    }
    fprintf(resFile, "========================================================\n");
    fclose(resFile);

    fprintf(logFile, "================ Ca-Ca distance  matrix ================\n");
    fprintf(logFile, "Ca-Ca distance  matrix calculated successfully.\n");
    fprintf(logFile, "Minimum and maximum distances are: 3.5 and 10 Angstroms.\n");
    fprintf(logFile, "Graphic display of the matrix can be found in the\n");
    fprintf(logFile, "caca_map.eps and caca_map_sdtv.eps files.\n");
    fprintf(logFile, "Raw data are store in the caca_map.txt file.\n");
    fprintf(logFile, "========================================================\n\n\n");

    printf("DONE\n");

    /***********************************************************\
    |*           Going back to the current directory           *|
    |*                   Cleaning up memory                    *|
    \***********************************************************/
    fclose(logFile);

    chdir(dirCurrent);

    _cleanAnalyze(&ENSData, &residueNbs, &statSecStr, &confNbsList, atomMasses, distances, ensFileName, dirAnalyze, dirConformers, dirCurrent, sequence);


    return 0;
}