/*******************************************************************
*  -/_|:|_|_\-
*
*  File:           crush.c
*
*  Function:       CRUSH: Maximum likelihood alignment and
*                  superpositioning of
*                  multiple macromolecular structures
*
*  Author(s):      Douglas L. Theobald
*                  Biochemistry Department
*                  Brandeis University
*                  MS 009
*                  415 South St
*                  Waltham, MA  02454-9110
*
*                  dtheobald@gmail.com
*                  dtheobald@brandeis.edu
*
*  Copyright:      Copyright (c) 2014-2015 Douglas L. Theobald
*
*  CRUSH is free software; you can redistribute it and/or modify
*  it under the terms of the GNU General Public License as
*  published by the Free Software Foundation; either version 3 of
*  the License, or (at your option) any later version.
*
*  CRUSH is distributed in the hope that it will be useful, but
*  WITHOUT ANY WARRANTY; without even the implied warranty of
*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
*  GNU General Public License for more details.
*
*  You should have received a copy of the GNU General Public
*  License along with CRUSH in the file 'COPYING'; if not, write
*  to the:
*
*  Free Software Foundation, Inc.,
*  59 Temple Place, Suite 330,
*  Boston, MA  02111-1307  USA
*
*  -/_|:|_|_\-
*
*******************************************************************/
#include "theseus.h"
#include "theseuslib.h"

#include "crush.h"
#include "StructAlign.h"
#include "NWfill.h"

static MSA*
AlignPairsThenMultiple(PDBCdsArray *, char **);

void
leave(int sig);

/* global declarations (necessary for leave(), I think) */
const gsl_rng_type     *T = NULL;
gsl_rng                *r2 = NULL;

CdsArray       *cdsA = NULL; /* main array of selected pdb cds, never modified */
PDBCdsArray    *pdbA = NULL; /* pdb file coordinate info, much of it unused */
Algorithm      *algo = NULL;
NWalgo         *nw_algo = NULL;
Statistics     *stats = NULL;

NWtable        *nw_table = NULL;
PDBCdsArray    *alignA = NULL;

int             crush_again;
double          ave_gapscore;
double         *egap = NULL;


static void
ParseCmdLine(int argc, char *argv[])
{
    int            option;
    extern char   *optarg;
    extern int     optopt;
    int            option_index = 0;
    int            i;
    size_t         cmdlinelen, argvlen;
    char           space[] = " ";

    /* save the command line, both in parsed form and as a whole string */
    algo->argc = argc;
    algo->argv = malloc(argc * sizeof(char *));

    algo->cmdline[0] = '\0';

    cmdlinelen = 0;
    for (i = 0; i < argc; ++i)
    {
        algo->argv[i] = NULL;
        argvlen = strlen(argv[i]);
        cmdlinelen += argvlen + 1;
        if (cmdlinelen >= 1024)
            break;
        strncat(algo->cmdline, argv[i], argvlen);
        strncat(algo->cmdline, space, 1);
        algo->argv[i] = malloc((argvlen+1) * sizeof(char));
        strncpy(algo->argv[i], argv[i], argvlen);
    }

    char short_options[] = "a:Ab:BCdD:e:Ef:Fg:GhHi:Ij:Jk:K:lLMnNo:Op:PqQ:r:R:s:ST:uUvVWxXyYz:Z"; /* Colon indicates that the option takes an argument. */

    struct option long_options[] =
    {
        {"alignment",    required_argument, 0, 'A'},
        {"amber",        no_argument,       0,  0 },
        {"bayes",        required_argument, 0, 'b'},
        {"covariance",   no_argument,       0, 'c'}, // aliased to the short_option 'c'
        {"fasta",        no_argument,       0, 'F'},
        {"help",         no_argument,       0, 'h'},
        {"info",         no_argument,       0, 'I'},
        {"invselection", required_argument, 0, 'S'},
        {"iterations",   required_argument, 0, 'i'},
        {"leastsquares", no_argument,       0, 'l'},
        {"mapfile",      required_argument, 0, 'M'},
        {"morphfile",    no_argument,       0, 'q'},
        {"noave"  ,      no_argument,       0, 'y'},
        {"nomp"   ,      no_argument,       0,  0 },
        {"notrans",      no_argument,       0, '0'},
        {"norot",        no_argument,       0, '1'},
        {"nohierarch",   no_argument,       0,  0 },
        {"nocovars",     no_argument,       0,  0 },
        {"orient",       required_argument, 0, 'o'},
        {"pca",          required_argument, 0, 'P'},
        {"precision",    required_argument, 0, 'p'},
        {"randgibbs",    no_argument,       0,  0 },
        {"rootname",     required_argument, 0, 'r'},
        {"scaleanchor",  required_argument, 0,  0 },
        {"seed",         required_argument, 0, 'X'},
        {"selection",    required_argument, 0, 's'},
        {"verbose",      optional_argument, 0, 'W'},
        {"version",      no_argument,       0, 'V'},
        {0,              0,                 0,  0 } // required zero array
    };

    while ((option = getopt_long(argc, argv, short_options, long_options, &option_index)) != -1)
    {
        switch (option)
        {
            case 0:
                if (strcmp(long_options[option_index].name, "notrans") == 0)
                {
                    algo->dotrans = 0;
                }
                else if (strcmp(long_options[option_index].name, "nohierarch") == 0)
                {
                    algo->dohierarch = 0;
                }
                else if (strcmp(long_options[option_index].name, "nocovars") == 0)
                {
                    algo->docovars = 0;
                }
                else if (strcmp(long_options[option_index].name, "nomp") == 0)
                {
                    algo->domp = 0;
                }
                else if (strcmp(long_options[option_index].name, "amber") == 0)
                {
                    algo->amber = 1; /* for dealing with fugly AMBER8 PDB format, with atom in wrong column */
                }
                else if (strcmp(long_options[option_index].name, "bayes") == 0)
                {
                    algo->bayes = (int) strtol(optarg, NULL, 10);
                }
                else if (strcmp(long_options[option_index].name, "scaleanchor") == 0)
                {
                    algo->scaleanchor = (int) strtol(optarg, NULL, 10);
                }
                else if (strcmp(long_options[option_index].name, "randgibbs") == 0)
                {
                    algo->randgibbs = 1;
                }
/*
                else
                {
                    printf("\n  Bad option '--%s' \n", long_options[option_index].name);
                    Usage(0);
                    exit(EXIT_FAILURE);
                }
 */
                break;
            case 'a':
                if (algo->alignment)
                {
                    printf("\n\n    -> Only alpha carbons can be selected <-");
                    printf("\n    -> when superimposing to an alignment <-\n");
                    algo->atoms = 0;
                }

                if (algo->pca > 0)
                {
                    printf("\n\n    -> Only alpha carbons can be selected <-");
                    printf("\n    -> when calculating principal components <-\n");
                    algo->atoms = 0;
                }

                if (isdigit(optarg[0]))
                {
                    algo->atoms = (int) strtol(optarg, NULL, 10);
                    if (algo->atoms > 4 || algo->atoms < 0)
                    {
                        printf("\n  Bad -a string \"%s\" \n", optarg);
                        Usage(0);
                        exit(EXIT_FAILURE);
                    }
                }
                else
                {
                    strtoupper(optarg);
                    algo->atomslxn = (char *) malloc((strlen(optarg) + 1) * sizeof(char));
                    mystrncpy(algo->atomslxn, optarg, FILENAME_MAX - 1);
                    algo->atoms = 10;
                }
                break;
            case 'A':
                nw_algo->angles = 1;
                break;
            case 'b':
                if (isalpha(*optarg))
                    nw_algo->blosum = (int) *optarg;
                else
                    nw_algo->blosum = (int) strtol(optarg, NULL, 10);
                break;
            case 'B':
                //nw_algo->bayespost = 1;
                nw_algo->bayes = 1;
                break;
            case 'C':
                nw_algo->center_ca = 1;
                break;
            case 'd':
                nw_algo->distance = 1;
                break;
            case 'D':
                nw_algo->align_dist = (double) strtod(optarg, NULL);
                break;
            case 'e': /* this is in terms of the fraction of the open penalty (see 'o' below) */
                nw_algo->extend = (double) strtod(optarg, NULL);
                nw_algo->gap = 1;
                break;
            case 'E': /* Electrostatic factors */
                nw_algo->raw_charge = 1;
                break;
            case 'f':
                nw_algo->fraglen = (int) strtol(optarg, NULL, 10);
                if (nw_algo->fraglen % 2 == 0)
                    nw_algo->fraglen++;
                break;
            case 'F':
                mystrncpy(nw_table->outfile_name, optarg, FILENAME_MAX-1);
                break;
            case 'g':
                if (isdigit(optarg[0]))
                    algo->hierarch = (int) strtol(optarg, NULL, 10);
                break;
            case 'G':
                nw_algo->gap = 1; /* use a gap penalty */
                break;
            case 'h':
                Usage(0);
                exit(EXIT_SUCCESS);
                break;
            case 'i':
                algo->iterations = (int) strtol(optarg, NULL, 10);
                break;
            case 'K':
                nw_algo->KA = 1;
                break;
            case 'l':
                nw_algo->logodds = 0;
                nw_algo->bayes = 0;
                break;
            case 'L':
                nw_algo->global = 0;
                break;
            case 'm':
                nw_table->msa_filename = (char *) malloc((strlen(optarg) + 2) * sizeof(char));
                mystrncpy(nw_table->msa_filename, optarg, strlen(optarg) + 1);

                cdsA->msafile_name = (char *) malloc((strlen(optarg) + 2) * sizeof(char));
                mystrncpy(cdsA->msafile_name, optarg, strlen(optarg) + 1);

                if (algo->atoms > 0)
                {
                    printf("\n\n    -> Only alpha carbons can be selected <-");
                    printf("\n    -> when superimposing to an alignment <-\n");
                    algo->atoms = 0;
                }

                algo->alignment = 1;
                algo->fmodel = 1;
                break;
                break;
            case 'M':
                nw_algo->logodds = 2;
                nw_algo->bayes = 0;
                break;
            case 'o':
                nw_algo->open = (double) strtod(optarg, NULL);
                nw_algo->gap = 1;
                break;
            case 'O':
                nw_algo->corder = 1;
                break;
            case 'p':
                algo->precision = (double) strtod(optarg, NULL);
                break;
            case 'P':
                nw_algo->procrustes = 1;
                break;
            case 's':
                algo->selection = (char *) malloc((strlen(optarg) + 1) * sizeof(char));
                mystrncpy(algo->selection, optarg, FILENAME_MAX - 1);
                break;
            case 'S':
                nw_algo->dsspflag = 1;
                break;
            case 'v':
                algo->varweight = 1;
                algo->covweight = 0;
                algo->leastsquares = 0;
                break;
            case 'V':
                Version();
                exit(EXIT_SUCCESS);
                break;
            case '?':
                printf("\n  Bad option '-%c' \n", optopt);
                Usage(0);
                exit(EXIT_FAILURE);
                break;
            default:
                abort();
        }
    }
}


static MSA*
AlignPairsThenMultiple(PDBCdsArray *myPDBs, char *argv[])
{
    typedef struct PairParser
    {
        int     slen; // String length (including dashes, both the common and paired sequences have this same length).
        int     posn; // Position that we're currently looking at (next character to parse).
        char   *commonSeq; // The actual common sequence (with dashes, if any) of this pair
        char   *pairedSeq; // The actual paired sequence (with dashes, if any)
    } PairParser;

    NWtable        *nwTable = NULL;
    PairParser     *pairParsers = NULL;
    const int       cnum = myPDBs->cnum;
    const int       numPairs = cnum - 1;
    int             ii, jj, jk, tmpLen, unmerged;
    int             iCommon, msaCommonLen, msaTotLen;
    char            currCh, altChar;

    pairParsers = calloc(sizeof(PairParser), numPairs);
    if (pairParsers == NULL)
    {
        // Punch out, Goose!!!
        perror("\n ERROR");
        printf("\n ERROR: Failed to allocate storage (%d pair parser structs) in AlignPairsThenMultiple().\n", cnum);
        exit(EXIT_FAILURE);
    }

    cdsA = GetNWCdsSel(myPDBs);

    // Store the original coordinates (Cds's), and find the longest.
    msaCommonLen = msaTotLen = 0L;
    iCommon = 0;
    for (ii = 0; ii < cnum; ii++)
    {
        tmpLen = cdsA->cds[ii]->vlen;
        msaTotLen += tmpLen;
        // Choose the longest struct, which all pairwise alignments will have in common.
        if (tmpLen > msaCommonLen)
        {
            msaCommonLen = tmpLen;
            iCommon = ii;
        }
    }
    msaTotLen = msaCommonLen;

    printf("    Aligning coordinates pairwise for superposition ... \n");
    fflush(NULL);

    nw_table = NWinit();
    nw_table->msa = MSAinit();
    MSAalloc(nw_table->msa, cnum, msaTotLen, SEQ_NAME_LENGTH);

    StructAlign(myPDBs, iCommon, cdsA, &nwTable); // Actually needs an array of numPairs nwTables here, not just a single nwTable!!!!!

#if (DEBUG == 1)
    for (ii = 0; ii < numPairs; ii++)
    {
        jj = pairParsers[ii].slen;
        if (jj > 40)
            jj = 40;

        currCh = pairParsers[ii].commonSeq[jj];
        pairParsers[ii].commonSeq[jj] = '\0'; // Premature stop. ;-)
        printf("\n%s", pairParsers[ii].commonSeq);
        pairParsers[ii].commonSeq[jj] = currCh;

        currCh = pairParsers[ii].pairedSeq[jj];
        pairParsers[ii].pairedSeq[jj] = '\0';
        printf("\n%s\n", pairParsers[ii].pairedSeq);
        pairParsers[ii].pairedSeq[jj] = currCh;
    }
#endif

    printf("\n    Creating multiple alignment for superposition ... \n");

    // Combine all of the pairwise alignments into one multiple sequence alignment w/a max length <= msaTotLen
    unmerged = numPairs; // A sequence is being worked on if the end position is > start and < pairParsers[jj].slen
    jk = 0;
    msaTotLen = 0;
    while (unmerged)
    {
        // First deal with any "pair" sequences that start before their "common" partner.
        tmpLen = 0;
        for (ii = 0; ii < numPairs; ii++)
        {
            int endLen = pairParsers[ii].slen;
            jj = pairParsers[ii].posn;
            while ((jj < endLen) && (pairParsers[ii].commonSeq[jj] == '-'))
            {
                // In this pair, the current char in the common sequence is a dash: append the character from the partnerSeq to the MSA.
                for (jk = 0; jk < numPairs; jk++)
                {
                    if (jk == ii)
                        nw_table->msa->seq[jk][msaTotLen] = pairParsers[ii].pairedSeq[jj];
                    else
                        nw_table->msa->seq[jk][msaTotLen] = '-';
                }
                nw_table->msa->seq[numPairs][msaTotLen] = '-'; // The common sequence goes last.
                msaTotLen++;
                tmpLen++;
                jj++;
            }
            pairParsers[ii].posn = jj;
            if (jj >= endLen)
                unmerged--;
        }

        if (unmerged == 0)
            continue;

        // Now update the MSA with the "common" sequences up to (but not including) the next dash.
        do {
            // Find a common sequence character (if any non-dash one(s) exist) to append to the MSA
            currCh = 0;
            for (ii = 0; ii < numPairs; ii++)
            {
                jj = pairParsers[ii].posn;
                if (jj < pairParsers[ii].slen)
                {
                    altChar = pairParsers[ii].commonSeq[jj];
                    if (currCh)
                    {
                        if (currCh != altChar)
                            currCh = '-'; // If they aren't the same, one must be a dash (and any dash is a non-starter).
                    }
                    else
                        currCh = altChar;
                    if (currCh == '-')
                        break; // Nothing left to look at here; skip ahead
                }
            }

            if (currCh != '-')
            {
                // We have a non-dash character in the current position of the common sequence of each active pair!
                nw_table->msa->seq[numPairs][msaTotLen] = currCh; // The common sequence.

                for (ii = 0; ii < numPairs; ii++)
                {
                    jj = pairParsers[ii].posn;
                    if (jj < pairParsers[ii].slen)
                    {
                        currCh = pairParsers[ii].pairedSeq[jj];
                        pairParsers[ii].posn++;
                        if (pairParsers[ii].posn >= pairParsers[ii].slen)
                            unmerged--;
                    }
                    else
                        currCh = '-';
                    nw_table->msa->seq[ii][msaTotLen] = currCh;
                }
                msaTotLen++;
            }
        } while (currCh != '-');
    }
    nw_table->msa->seqlen = msaTotLen;

    for (ii = 0; ii < numPairs; ii++)
    {
        free(pairParsers[ii].pairedSeq);
        free(pairParsers[ii].commonSeq);
    }
    free(pairParsers);

    return nw_table->msa;
}


int
main(int argc, char *argv[])
{
    char           *sup_name = NULL, *ave_name = NULL, *transf_name = NULL,
                   *rand_transf_name = NULL, *mean_ip_name = NULL, *sup_var_name = NULL;
    extern int      optind;
    const int       kMaxCrushableStructs = 25;
    int             cnum;

    signal(SIGINT, leave);
    signal(SIGABRT, leave);

//    init_genrand(seed);

    cdsA = CdsArrayInit();
    algo = AlgorithmInit();
    stats = StatsInit();
    nw_algo = NWalgoInit();

    crush_again = 0;
    ave_gapscore = 0.0;

    ParseCmdLine(argc, argv);
    algo->infiles = &argv[optind];
    algo->filenum = argc - optind;

    if (algo->filenum < 1)
    {
        printf("\n  -> No pdb files specified. <- \n");
        Usage(0);
        exit(EXIT_FAILURE);
    }
    else if (algo->filenum > 2)
    {
        printf("\n  More than two pdb files were specified-- a multiple sequence alignment will be attempted.\n");
    }

    PrintTheseusPre();

    if (algo->filenum == 1)
        printf("    Reading pdb file ... \n");
    else
        printf("    Reading %d pdb files ... \n", algo->filenum);
    fflush(NULL);

    pdbA = GetPDBCds(algo->infiles, algo->filenum, algo->fmodel, algo->amber, algo->atom_names);
    /* PrintPDBCds(stdout, pdbA->cds[0]); */

    cnum = pdbA->cnum;

    if (cnum < 2)
    {
        printf("\n  -> Found less than two PDB cds. Could not do superposition. <- \n");
        Usage(0);
        exit(EXIT_FAILURE);
    }

    printf("    Successfully read %d models and/or structures \n", cnum);

    egap = calloc(cnum, sizeof(double));

    if (cnum == 2) // Exactly two structures.
    {
        if (algo->alignment)
        {
            printf("    Reading multiple sequence alignment ... \n");
            fflush(NULL);

            Align2MSA(pdbA, cdsA, cdsA->msafile_name, cdsA->mapfile_name);
            PDBCdsArrayAllocLen(pdbA, cdsA->vlen);
        }
        else
        {   // Exactly two structures to align.
            printf("    Aligning structures for superposition ... \n");
            fflush(NULL);

            nw_table = NWinit();
            cdsA = GetNWCdsSel(pdbA);
            StructAlign(pdbA, 0, cdsA, &nw_table);

            strncpy(nw_table->msa->name[0], argv[0 + optind], strlen(argv[0 + optind]) + 1);
            strncpy(nw_table->msa->name[1], argv[1 + optind], strlen(argv[1 + optind]) + 1);
            //strcpy(nw_table->msa->name[2], " ");
            nw_table->msa->seqnum = 2;
            writealn(nw_table->msa, 0, nw_table->msa->seqlen, "theseus");
            writea2m(nw_table->msa, 0, nw_table->msa->seqlen, "theseus");
        }
    }
    else if (cnum > kMaxCrushableStructs)
    {
        printf("\n  -> Crush is currently limited to %d or fewer structures. <- \n", kMaxCrushableStructs);
        exit(EXIT_FAILURE);
    }
    else
    {   // Aligning multiple structures.
        MSA *       mergemsa = AlignPairsThenMultiple(pdbA, argv);
#if (DEBUG == 1)
        int         ii;
        char        wasChar;

        for (ii = 0; ii < cnum; ii++)
        {
            wasChar = mergemsa->seq[ii][30];
            mergemsa->seq[ii][30] = '\0';
            printf("\n%s", mergemsa->seq[ii]);
            mergemsa->seq[ii][30] = wasChar;
        }
        printf("\n");
#endif
        writea2m(mergemsa, 0, mergemsa->seqlen, "CrushMSA");

        nw_table->alignment = malloc(sizeof(int *));  // The destroyer wants to see these.
        nw_table->alignment[0] = malloc(sizeof(int)); // The destroyer wants to see these.
        NWdestroy(&nw_table);
    }

    printf("\n    Done.\n");
    printf("I-===-===-===-===-===-===-===-===-===-===-===-===-===-===-===-===-===I\n");
    printf("                          < END CRUSH %s >                          \n", VERSION);
    fflush(NULL);

    //AlgorithmDestroy(&algo);
    //free(stats);
    //stats = NULL;
    //CdsArrayDestroy(&cdsA);
    //PDBCdsArrayDestroy(&pdbA);

    return (EXIT_SUCCESS);

    /* //////////////////////////////////////////////////////////////////////////////// */
    /* //////////////////////////////////////////////////////////////////////////////// */

    PrintTheseusTag();

    AlgorithmDestroy(&algo);
    free(stats);
    stats = NULL;

    CdsArrayDestroy(&cdsA);
    PDBCdsArrayDestroy(&pdbA);

    if (sup_name)
        free(sup_name);
    if (sup_var_name)
        free(sup_var_name);
    if (ave_name)
        free(ave_name);
    if (transf_name)
        free(transf_name);
    if (rand_transf_name)
        free(rand_transf_name);
    if (mean_ip_name)
        free(mean_ip_name);

    exit(EXIT_SUCCESS);
}



void
leave(int sig)
{
    if (algo->rounds > 0)
    {
        printf("    Aborting at iteration %d ....\n", algo->rounds+1);
        fflush(NULL);

        algo->abort = 1;
        signal(sig, SIG_IGN);
    }
    else
    {
        fflush(NULL);
        signal(sig, SIG_DFL);
    }
}