#include "main_cu.h"
#include "neighbour_joining.h"
#include "back_up_struct.h"
#include <math.h>

#define PRI 0
#define EXT 1

//struct LibRead
//{
//    int x;   // index within first sequence (smaller one)
//    int y;   // index within second sequence
//    int idx; // index of a cell in library from which the value was read
//};
//
//struct Position
//{
//    int x; // position in original sequence
//    int p; // how many gap elements left (before we get to the next letter)
//};

MultipleSequenceAlignment::MultipleSequenceAlignment()
{
    elb = new ExtendedLibraryBuilder();
    plb = new PrimaryLibraryBuilder();
}

MultipleSequenceAlignment::~MultipleSequenceAlignment()
{
    delete elb;
    delete plb;
}

const char* MultipleSequenceAlignment::getAlgorithmName()
{
    return "msa";
}

char* MultipleSequenceAlignment::defaultSettingsFile()
{
    return "etc/msa-defaults.cnf";
}

void MultipleSequenceAlignment::actuallyArgumentsLoading(ArgumentsManager* am)
{
    // works only on Linux, Solaris and AIX
    this->nCPUs = sysconf(_SC_NPROCESSORS_ONLN);

    if(am->containParam("ncpus", "cpus"))
        this->nCPUs = min(am->getIntParam("ncpus", "cpus"), this->nCPUs);

    if(am->containParam("msagapop", "msago"))
        this->msaGapOp = am->getIntParam("msagapop", "msago");
    if(am->containParam("msagapex", "msage"))
        this->msaGapEx = am->getIntParam("msagapex", "msage");

    if(am->containParam("msaoutclustal", "msaoc"))
        this->msaoutclustal          = am->getParam("msaoutclustal", "msaoc");
    else
        this->msaoutclustal          = NULL;

    if(am->containParam("msaoutfasta", "msaof"))
        this->msaoutfasta            = am->getParam("msaoutfasta", "msaof");
    else
        this->msaoutfasta            = NULL;


    this->libType = EXT;
    if(am->containParam("library", "lib"))
    {
        const char* library = am->getParam("library", "lib");
        if(!strcmp(library, "primary"))
            this->libType = PRI;
        if(!strcmp(library, "extended"))
            this->libType = EXT;
    }

    if(this->libType == EXT)
    {
        elb->actuallyArgumentsLoading(am);
        elb->plb->computeNJ = true; // we need NJ tree to compute MSA
    }
    else if(this->libType == PRI)
    {
        plb->actuallyArgumentsLoading(am);
        plb->computeNJ = true; // we need NJ tree to compute MSA
    }
}

void helloCPU(ThreadManager* This, void* data)
{
    int cpuNo = This->getThreadsInfo().gpuNo;
    printf("CPU core: %d\n",cpuNo);
}

void MultipleSequenceAlignment::saveToFile()
{

    if(msaoutclustal)
    {
        printf("Saving MSA to file %s ", msaoutclustal);
        HiResTimer timer;
        timer.start();

        int maxSeqNameLength = 0;
        for (int i = 0; i < sequenceNumber; i++)
        {
            maxSeqNameLength = MAX(maxSeqNameLength, strlen(seqs->getSeqName(i)));
        }
        char* seqFormat = new char[16]; //"%-23s\t"
        sprintf(seqFormat, "%%-%ds  ", maxSeqNameLength);


        FILE* file = fopen(msaoutclustal, "w");

        fprintf(file, "CLUSTAL FORMAT for GPU-T-COFFEE, SCORE=%d, Nseq=%d, Len=%d\n\n",
                     (int)scoreMSA, sequenceNumber, totalLengthMSA);

        int* carret = new int[sequenceNumber];
        int* letter = new int[sequenceNumber];
        for(int i=0; i<sequenceNumber; i++)
        {
            carret[i] = 0;
            letter[i] = 0;
        }

        int LIR = 50; // LETTERS_IN_ROW

        // CONSERVATION INFORMATION
        char* alignment = new char[LIR*sequenceNumber];
        const char* strongGroups[] = {"STA", "NEQK", "NHQK", "NDEQ", "QHRK", "MILV", "MILF", "HY", "FYW"};
        const char* weakGroups[] = {"CSA", "ATV", "SAG", "STNK", "STPA", "SGND", "SNDEQK", "NDEQHK", "NEQHRK", "FVLIM", "HFY"};
        int strongGroupCount = 9;
        int weakGroupCount = 11;


        for(int offset=0; offset<totalLengthMSA; offset+=LIR)
        {
            for(int i=0; i<sequenceNumber; i++)
            {
                fprintf(file, seqFormat, seqs->getSeqName(i));

                int start  = seqs->getStarts() [i];
                int length = seqs->getLengths()[i];
                for (;carret[i] < offset + LIR; carret[i]++)
                {
                    if (carret[i] >= totalLengthMSA)
                        break;

                    if (letter[i] >= length)
                    {
                        alignment[i*LIR + carret[i]%LIR] = '-';
                        fprintf(file, "-");
                        continue;
                    }

                    if (carret[i]==MSA[letter[i] + start])
                    {
                        alignment[i*LIR + carret[i]%LIR] =
                                seqs->getSubtitutionMatrix()->revConvert(
                                seqs->getSequences()[start+length - letter[i] - 1]);


                        fprintf(file, "%c", alignment[i*LIR + carret[i]%LIR]);
                        letter[i]++;
                    }
                    else
                    {
                        alignment[i*LIR + carret[i]%LIR] = '-';
                        fprintf(file, "-");
                    }
                }
                fprintf(file, " %d\n", letter[i]);
            }

            // write conservation line
            fprintf(file, seqFormat, "");
            for(int i=0; i<LIR; i++)
            {
                if(offset+i >= totalLengthMSA)
                    break;

                int groupCount;
                int identCount = 1;
                bool done = false;

                for(int j=1; j<sequenceNumber; j++)
                {
                    if(alignment[j*LIR + i] == alignment[i])
                        identCount++;
                }
                if(identCount == sequenceNumber)
                {
                    fprintf(file, "*"); // all residues in this column are identical
                    continue;
                }
                for(int gr=0; gr<strongGroupCount; gr++)
                {
                    groupCount = 0;
                    for (int let=0; strongGroups[gr][let]; let++) // let - letter in group
                    {
                        for(int j=0; j<sequenceNumber; j++)
                        {
                            if(alignment[j*LIR + i] == strongGroups[gr][let])
                                groupCount++;
                        }
                    }
                    if(groupCount == sequenceNumber)
                    {
                        done = true;
                        fprintf(file, ":"); // all residues in this column are in strong group
                        break;
                    }
                }
                if(done)
                    continue;

                for(int gr=0; gr<weakGroupCount; gr++)
                {
                    groupCount = 0;
                    for (int let=0; weakGroups[gr][let]; let++) // let - letter in group
                    {
                        for(int j=0; j<sequenceNumber; j++)
                        {
                            if(alignment[j*LIR + i] == weakGroups[gr][let])
                                groupCount++;
                        }
                    }
                    if(groupCount == sequenceNumber)
                    {
                        done = true;
                        fprintf(file, "."); // all residues in this column are in weak group
                        break;
                    }
                }
                if(done)
                    continue;

                fprintf(file, " ");

            }

            fprintf(file, "\n\n");
        }



        fclose(file);

        timer.stop();
        printf("(%dms)\n", (int)timer.getElapsedTime());
    }

    if(msaoutfasta)
    {
        printf("Saving MSA to file %s ", msaoutfasta);
        HiResTimer timer;
        timer.start();

        FILE* file = fopen(msaoutfasta, "w");
        
        int LIR = 50; // LETTERS_IN_ROW

        for(int i=0; i<sequenceNumber; i++)
        {
            int carret = 0;
            int start  = seqs->getStarts() [i];
            int length = seqs->getLengths()[i];

            fprintf(file, ">%s\n", seqs->getSeqName(i));

            for(int j=0; j<length; j++)
            {
                for(int k=carret; k<MSA[start + j]; k++)
                {
                    fprintf(file, "-");
                    carret++;
                    if(!(carret%LIR))
                        fprintf(file, "\n");
                }
                fprintf(file, "%c",
                       seqs->getSubtitutionMatrix()->revConvert(seqs->getSequences()[start+length - j - 1])
                        );
                carret ++;
                if(!(carret%LIR))
                    fprintf(file, "\n");
            }
            for(int j=carret; j<totalLengthMSA; j++)
            {
                fprintf(file, "-");
                carret++;
                if(!(carret%LIR))
                    fprintf(file, "\n");
            }
            if(carret%LIR)
                fprintf(file, "\n");
        }

        fclose(file);

        timer.stop();
        printf("(%dms)\n", (int)timer.getElapsedTime());
    }




//    // printing alignment (without conservation)
//    for(int i=0; i<sequenceNumber; i++)
//    {
//        int carret = 0;
//        int start  = seqs->getStarts() [i];
//        int length = seqs->getLengths()[i];
//
//        printf("seq %3d ", i);
//
//        for(int j=0; j<length; j++)
//        {
//            for(int k=carret; k<MSA[start + j]; k++)
//            {
//                printf("-");
//            }
//            printf("%c",
//                   seqs->getSubtitutionMatrix()->revConvert(seqs->getSequences()[start+length - j - 1])
//                    );
//            carret = MSA[start + j] + 1;
//        }
//        for(int j=carret; j<totalLengthMSA; j++)
//            printf("-");
//        printf("\n");
//    }

}

void MultipleSequenceAlignment::run(ThreadManager* tm)
{
    if(this->libType == EXT)
    {
        elb->run(tm);
        library = elb->plb->library; // should be casted to (ExtendedLibrary*) when needed
        seqs    = elb->plb->seqs;
        dm      = elb->plb->dm;
        nj      = elb->plb->nj;
    }
    else if(this->libType == PRI)
    {
        plb->run(tm);
        library = plb->library;
        seqs    = plb->seqs;
        dm      = plb->dm;
        nj      = plb->nj;
    }


    int cores = sysconf(_SC_NPROCESSORS_ONLN);
    if(cores > nCPUs)
        printf("Cores: %d, using: %d\n", cores, nCPUs);
    ThreadManager* tmCPU;
    int* cpus = new int[nCPUs];
    for(int i=0; i<nCPUs; i++)
        cpus[i] = i;
    tmCPU = new ThreadManager(cpus,nCPUs);
    for(int i=0; i<nCPUs; i++)
        tmCPU->request(helloCPU, NULL, i);
    tmCPU->wait();

    
    //HERE TASKS ARE DEFINED

    HiResTimer timer;
    timer.start();


    resPen = new double[32];
    for(int i=0; i<32; i++)
        resPen[i] = 1;
    resPen[ seqs->getSubtitutionMatrix()->convert('A') ] = 1.13;
    resPen[ seqs->getSubtitutionMatrix()->convert('C') ] = 1.13;
    resPen[ seqs->getSubtitutionMatrix()->convert('D') ] = 0.96;
    resPen[ seqs->getSubtitutionMatrix()->convert('E') ] = 1.31;
    resPen[ seqs->getSubtitutionMatrix()->convert('F') ] = 1.20;
    resPen[ seqs->getSubtitutionMatrix()->convert('G') ] = 0.61;
    resPen[ seqs->getSubtitutionMatrix()->convert('H') ] = 1.00;
    resPen[ seqs->getSubtitutionMatrix()->convert('I') ] = 1.32;
    resPen[ seqs->getSubtitutionMatrix()->convert('K') ] = 0.96;
    resPen[ seqs->getSubtitutionMatrix()->convert('L') ] = 1.21;
    resPen[ seqs->getSubtitutionMatrix()->convert('M') ] = 1.29;
    resPen[ seqs->getSubtitutionMatrix()->convert('N') ] = 0.63;
    resPen[ seqs->getSubtitutionMatrix()->convert('P') ] = 0.74;
    resPen[ seqs->getSubtitutionMatrix()->convert('Q') ] = 1.07;
    resPen[ seqs->getSubtitutionMatrix()->convert('R') ] = 0.72;
    resPen[ seqs->getSubtitutionMatrix()->convert('S') ] = 0.76;
    resPen[ seqs->getSubtitutionMatrix()->convert('T') ] = 0.89;
    resPen[ seqs->getSubtitutionMatrix()->convert('V') ] = 1.25;
    resPen[ seqs->getSubtitutionMatrix()->convert('Y') ] = 1.00;
    resPen[ seqs->getSubtitutionMatrix()->convert('W') ] = 1.23;

//    for(int i=0; i<32; i++)
//        printf("%f\n", resPen[i]);


    sequenceNumber = seqs->getSequenceNumber();
    
    // Results of MSA will be held in MSA array.
    // MSA array holds informations about
    // real position in final alignment
    // (including gaps) of each residue
    // for every input sequence.
    // MSA array aggregates then gaps
    // between consecutive residues.
    int seqsEntireLength = seqs->getEntireLength();
    MSA = new short[seqsEntireLength];
    for(int i=0; i<sequenceNumber; i++)
    {
        int startPos = seqs->getStarts()[i];
        int elements = seqs->getLengths()[i];
        for(int j=0; j<elements; j++)
            MSA[startPos + j] = j;
    }


    int partId = 1;
    MSAThreadParams* params;

    int maxTreeHeight = nj->root->maxHeight();
    for(int i=maxTreeHeight; i>=1; i--)
    {
        INJTreeNodeList* nodesCollection = nj->root->collectNodes(i);
        INJTreeNodeList* iterator = nodesCollection;
        while (iterator != NULL)
        {
            params = new MSAThreadParams();
            params->partId = partId++;
            params->node = (NJInternalNode*)iterator->current;

            RunnableWithParams* invokingParams = new RunnableWithParams();
            invokingParams->algorithm = this;
            invokingParams->params = params;
            tmCPU->request(invoker, (void*)invokingParams, -1);

            iterator = iterator->next;
        }

        tmCPU->wait();
    }


    timer.stop();
    printf("%s total: %dms\n", getAlgorithmName(), (int) timer.getElapsedTime());

    totalLengthMSA = params->node->length;
    saveToFile();

}

void MultipleSequenceAlignment::actualInvokedMethod(void* params)
{
    MSAThreadParams* param = (MSAThreadParams*)params;

    INJTreeNode* nodeLeft  = param->node->childNodeLeft;
    INJTreeNode* nodeUp    = param->node->childNodeRight;
    
    INJSequencesNoList* seqsLeftList = nodeLeft->collectSequenceNumbers();
    int  seqsLeftCount = seqsLeftList->count;
    int* seqsLeft      = seqsLeftList->toArray();

    INJSequencesNoList* seqsUpList = nodeUp->collectSequenceNumbers();
    int  seqsUpCount = seqsUpList->count;
    int* seqsUp      = seqsUpList->toArray();

    delete seqsLeftList;
    delete seqsUpList;



//    // printing sequences names
//    printf("Up sequences: \n");
//    for (int seqU = 0; seqU < seqsUpCount; seqU++)
//    {
//        printf("%s\n", seqs->getSeqName(seqsUp[seqU]));
//    }
//    printf("Left sequences: \n");
//    for (int seqL = 0; seqL < seqsLeftCount; seqL++)
//    {
//        printf("%s\n", seqs->getSeqName(seqsLeft[seqL]));
//    }
//    printf("\n");
//    printf("\n");

//    //printing sequences number
//    printf("Left(%d):\n", seqsLeftCount);
//    for (int i = 0; i < seqsLeftCount; i++)
//        printf("%4d", seqsLeft[i]);
//    printf("\nRight(%d):\n", seqsRightCount);
//    for (int i = 0; i < seqsRightCount; i++)
//        printf("%4d", seqsRight[i]);
//    printf("\n\n");


    //printf("%fMB\n", (nodeLeft->length * nodeUp->length * 4)/(1024.0*1024.0));
    float* SM = new float[nodeLeft->length * nodeUp->length];

    for (int i = 0; i < nodeLeft->length * nodeUp->length; i++)
        SM[i] = 0;

    unsigned int start;
    unsigned int end;
    unsigned int tmp;
    unsigned int posUp; // pos - position
    unsigned int posLeft;
    unsigned int val;

    double arms = 0; // average residue mismatch score (see art. about ClustalW)
    unsigned int armsUpdates = 0;
    double avgIdentityPercent = 0;

    for (int seqL = 0; seqL < seqsLeftCount; seqL++)
    {
        for (int seqU = 0; seqU < seqsUpCount; seqU++)
        {
            avgIdentityPercent += (1.0 -dm->getElement(seqsLeft[seqL], seqsUp[seqU]))
                                  / (seqsLeftCount * seqsUpCount);

            if (seqsUp[seqU] > seqsLeft[seqL])
            {
                start = library->starts[seqsLeft[seqL] + seqsUp[seqU] * library->seqCount];
                end   = library->starts[seqsLeft[seqL] + seqsUp[seqU] * library->seqCount + 1];
            }
            else // (seqsUp[seqU] < seqsLeft[seqL])
            {
                start = library->starts[seqsUp[seqU] + seqsLeft[seqL] * library->seqCount];
                end   = library->starts[seqsUp[seqU] + seqsLeft[seqL] * library->seqCount + 1];
            }



            // i - position inside the library
            for (int i = start; i < end; i++)
            {
                tmp = library->pl[i];
                if (tmp == 0xFFFFFFFF)
                    break;

                if (seqsUp[seqU] > seqsLeft[seqL])
                {
                    posUp   = (tmp >> 8)&0xFFF;
                    posLeft = tmp >> 20;
                }
                else
                {
                    posUp   = tmp >> 20;
                    posLeft = (tmp >> 8)&0xFFF;
                }

                if(this->libType == EXT)
                {
                    val = ((ExtendedLibrary*)library)->el[i];
                }
                else if(this->libType == PRI)
                {
                    val = (tmp & 0xFF);
                }
                if(posUp != posLeft)
                {
                    arms += val;
                    armsUpdates++;
                }

                // to only for tests (should be commented later)
                if(posLeft >= seqs->getLengths()[seqsLeft[seqL]])
                    printf("Error: index exceeded: max=%d, val=%d\n", seqs->getLengths()[seqsLeft[seqL]]-1, posLeft);
                if(posUp >= seqs->getLengths()[seqsUp[seqU]])
                    printf("Error: index exceeded: max=%d, val=%d\n", seqs->getLengths()[seqsUp[seqU]]-1, posUp);

                posLeft = MSA[ seqs->getStarts()[seqsLeft[seqL]] + posLeft ];
                posUp   = MSA[ seqs->getStarts()[seqsUp  [seqU]] + posUp   ];

                SM[posUp + posLeft * nodeUp->length] += (float)
                   (((double)val)/(seqsLeftCount * seqsUpCount)
                     * nj->seqsWeights[seqsLeft[seqL]]       // sequences are weighted
                     * nj->seqsWeights[seqsUp  [seqU]]);
            }
        }
    }

//    // Print identity percent and sequences
//    printf("avgIdentityPercent: %f\n", avgIdentityPercent);
//    if(avgIdentityPercent == 1.0)
//    {
//        for (int seqU = 0; seqU < seqsUpCount; seqU++)
//        {
//            int start  = seqs->getStarts() [ seqsUp[seqU] ];
//            int length = seqs->getLengths()[ seqsUp[seqU] ];
//            for(int i = 0; i < length; i++)
//            {
//                printf("%c", seqs->getSubtitutionMatrix()->revConvert(seqs->getSequences()[start + i]));
//            }
//            printf("\n");
//        }
//        printf("\n");
//        for (int seqL = 0; seqL < seqsLeftCount; seqL++)
//        {
//            int start  = seqs->getStarts() [ seqsLeft[seqL] ];
//            int length = seqs->getLengths()[ seqsLeft[seqL] ];
//            for(int i = 0; i < length; i++)
//            {
//                printf("%c", seqs->getSubtitutionMatrix()->revConvert(seqs->getSequences()[start + i]));
//            }
//            printf("\n");
//        }
//    }

    //arms /= (nodeLeft->length*nodeUp->length);
    arms /= armsUpdates;
    arms /= (nodeLeft->length+nodeUp->length)/16.0;
    //printf("ARMS: %f\n", arms);


    // Modyfication of gap penalties:
    
    double gapOp = ( msaGapOp + log(min(nodeLeft->length, nodeUp->length)) ) * 2 * avgIdentityPercent;
    gapOp *= arms;
    double gapEx = msaGapEx * (1.0f + abs(log((float)nodeLeft->length/(float)nodeUp->length)) );
    //printf("New values: gapOpen=%5.3f, gapExt=%5.3f\n", gapOp, gapEx);

    double* gapOpsLeft = new double[nodeLeft->length];
    double* gapExsLeft = new double[nodeLeft->length];
    double* gapOpsUp   = new double[nodeUp->length];
    double* gapExsUp   = new double[nodeUp->length];
    for(int i=0; i<nodeLeft->length; i++)
    {
        gapOpsLeft[i] = 0;
        gapExsLeft[i] = 0;
    }
    for(int i=0; i<nodeUp->length; i++)
    {
        gapOpsUp[i] = 0;
        gapExsUp[i] = 0;
    }

    // LOWERED GAP PENALTIES AT EXISTING GAP
    for (int seqU = 0; seqU < seqsUpCount; seqU++)
    {
        int start  = seqs->getStarts() [ seqsUp[seqU] ];
        int length = seqs->getLengths()[ seqsUp[seqU] ];

        int pos = 0;

        for(int i = 0; i < length; i++)
        {
            // 0 1 2 5 7
            // 0 1 2 _ _ 5 _ 7
            for (int j = pos; j < MSA[start + i]; j++)
            {
                gapOpsUp[j]++;
            }
            pos = MSA[start + i] + 1;
        }

        for (int i = pos; i<nodeUp->length; i++)
        {
            gapOpsUp[i]++;
        }
    }

    for (int seqL = 0; seqL < seqsLeftCount; seqL++)
    {
        int start  = seqs->getStarts() [ seqsLeft[seqL] ];
        int length = seqs->getLengths()[ seqsLeft[seqL] ];

        int pos = 0;

        for(int i = 0; i < length; i++)
        {
            // 0 1 2 5 7
            // 0 1 2 _ _ 5 _ 7
            for (int j = pos; j < MSA[start + i]; j++)
            {
                gapOpsLeft[j]++;
            }
            pos = MSA[start + i] + 1;
        }

        for (int i = pos; i<nodeLeft->length; i++)
        {
            gapOpsLeft[i]++;
        }
    }


    // INCREASED GAP PENALTIES NEAR EXISTING GAPS
    for (int i = 0; i<nodeUp->length; i++)
    {
        if (gapOpsUp[i] > 0)
        {
            int start = max(i-7, 0);
            int stop  = min(i+8, nodeUp->length);
            for (int j = start; j < stop; j++)
            {
                if (gapOpsUp[j] == 0)
                {
                    gapExsUp[j] = max(gapExsUp[j], 8.0 - abs(i-j));
                }
            }
        }
    }
    for (int i = 0; i<nodeLeft->length; i++)
    {
        if (gapOpsLeft[i] > 0)
        {
            int start = max(i-7, 0);
            int stop  = min(i+8, nodeLeft->length);
            for (int j = start; j < stop; j++)
            {
                if (gapOpsLeft[j] == 0)
                {
                    gapExsLeft[j] = max(gapExsLeft[j], 8.0 - abs(i-j));
                }
            }
        }
    }

    // NO (FOR NOW) REDUCED GAP PENALTIES IN HYDROPHILIC STRECHES

    // RESIDUE-SPECIFIC PENALTIES
    double* resPenUp   = new double[ nodeUp->length ];
    double* resPenLeft = new double[nodeLeft->length];

    for (int i = 0; i < nodeUp->length; i++)
        resPenUp[i] = 0;
    for (int i = 0; i < nodeLeft->length; i++)
        resPenLeft[i] = 0;

    for (int seqU = 0; seqU < seqsUpCount; seqU++)
    {
        int start  = seqs->getStarts() [ seqsUp[seqU] ];
        int length = seqs->getLengths()[ seqsUp[seqU] ];

        for (int i = 0; i < length; i++)
        {
            if (gapOpsUp[ MSA[start + i] ] == 0)
            {
                resPenUp[ MSA[start + i] ] += resPen[ seqs->getSequences()[start + i] ] / seqsUpCount;
            }
        }
    }
    for (int seqL = 0; seqL < seqsLeftCount; seqL++)
    {
        int start  = seqs->getStarts() [ seqsLeft[seqL] ];
        int length = seqs->getLengths()[ seqsLeft[seqL] ];

        for (int i = 0; i < length; i++)
        {
            if (gapOpsLeft[ MSA[start + i] ] == 0)
            {
                resPenLeft[ MSA[start + i] ] += resPen[ seqs->getSequences()[start + i] ] / seqsLeftCount;
            }
        }
    }
    

    // FINAL COMPUTATION OF PENALTIES
    for (int i = 0; i<nodeUp->length; i++)
    {
        if (gapOpsUp[i] > 0)
        {
            gapOpsUp[i] = gapOp * 0.3 * ((seqsUpCount - gapOpsUp[i]) / seqsUpCount);
            gapExsUp[i] = gapEx * 0.5;
        }
        else // no gap at this position
        {
            gapOpsUp[i] = gapOp * resPenUp[i]; // RESIDUE-SPECIFIC PENALTIES
            if(resPenUp[i] == 0)
                printf("Error: something was wrong with residue specific penalties.\n");

            if(gapExsUp[i] > 0) // INCREASED GAP PENALTIES NEAR EXISTING GAPS
            {
                gapOpsUp[i] = gapOp * (2 + gapExsUp[i]/4.0);
                gapOpsUp[i] *= resPenUp[i];
            }
            
            gapExsUp[i] = gapEx;

        }
    }
    for (int i = 0; i<nodeLeft->length; i++)
    {
        if (gapOpsLeft[i] > 0)
        {
            gapOpsLeft[i] = gapOp * 0.3 * ((seqsLeftCount - gapOpsLeft[i]) / seqsLeftCount);
            gapExsLeft[i] = gapEx * 0.5;
        }
        else // no gap at this position
        {
            gapOpsLeft[i] = gapOp * resPenLeft[i]; // RESIDUE-SPECIFIC PENALTIES
            if(resPenLeft[i] == 0)
                printf("Error: something was wrong with residue specific penalties.\n");

            if(gapExsLeft[i] > 0) // INCREASED GAP PENALTIES NEAR EXISTING GAPS
            {
                gapOpsLeft[i] = gapOp * (2 + gapExsLeft[i]/4.0);
                gapOpsLeft[i] *= resPenLeft[i];
            }

            gapExsLeft[i] = gapEx;
        }
    }

    // END OF PENALTIES CHANGES
    //==========================================================================
    // DYNAMIC PROGRAMMING

    float *A;
    float *E; //left matrix
    float *F; //up matrix
    BackUpStruct *B;

    A = new float[nodeUp->length + 1];
    E = new float[nodeUp->length + 1]; //left matrix
    F = new float[nodeUp->length + 1]; //up matrix
    B = new BackUpStruct[(nodeUp->length + 1)*(nodeLeft->length + 1)];

    // Initialization of the arrays

    A[0] = 0;
    E[0] = INT_MIN;
    F[0] = INT_MIN;
    B[0].backDirection = stop;
    int LEN = nodeUp->length + 1;
    for (int i = 1; i < nodeUp->length + 1; i++)
    {
        A[i] = gapOpsUp[0] - i*gapExsUp[i-1];
        E[i] = INT_MIN;
        F[i] = INT_MIN;
        B[i].backDirection = ::left;
    }
    for (int j = 1; j < nodeLeft->length + 1; j++)
    {
        B[j*LEN].backDirection = up;
    }

    // Filling of the DTM matrix
    
    /*
     *   s e q U p
     * s
     * e
     * q
     * L
     * e
     * f
     * t
     */
    //E - responsible for left direction
    //F - responsible for up   direction
    float bufA;
    float resF;
    float diagVal;

    for (int j = 1; j <= nodeLeft->length; j++) // Y
    {
        A[0] = -j*gapEx;
        E[0] = INT_MIN;
        F[0] = INT_MIN;
        bufA = -(j - 1)*gapEx;

        //if (param->partId == 1)
        //        printf("%5.1f ", A[0]);

        for (int i = 1; i <= nodeUp->length; i++) // X
        {

            E[i] = MAX(E[i-1] - gapExsUp[i-1], A[i-1] - gapOpsUp[i-1] - gapExsUp[i-1]);
            B[j*LEN + i - 1].continueLeft = (E[i] == E[i-1] - gapExsUp[i-1]);
            resF = MAX(F[i] - gapExsLeft[j-1], A[i] - gapOpsLeft[j-1] - gapExsLeft[j-1]);
            B[(j-1)*LEN + i].continueUp = (resF == F[i] - gapExsLeft[j-1]);
            F[i] = resF;
//            E[i] = MAX(E[i-1] - gapEx, A[i-1] - gapOp - gapEx);
//            B[j*LEN + i - 1].continueLeft = (E[i] == E[i-1] - gapEx);
//            resF = MAX(F[i] - gapEx, A[i] - gapOp - gapEx);
//            B[(j-1)*LEN + i].continueUp = (resF == F[i] - gapEx);
//            F[i] = resF;

//            //no penalties
//            E[i] = MAX(E[i-1], A[i-1]);
//            B[j*LEN + i - 1].continueLeft = (E[i] == E[i-1]);
//            resF = MAX(F[i], A[i]);
//            B[(j-1)*LEN + i].continueUp = (resF == F[i]);
//            F[i] = resF;

            diagVal = bufA + SM[ (j-1)*nodeUp->length + (i-1) ];
            bufA = A[i];

            A[i] = MAX3(E[i], F[i], diagVal);


            if(A[i] == diagVal)
                B[j*LEN+i].backDirection = crosswise;
            else if(A[i] == E[i])
               B[j*LEN+i].backDirection = ::left;
            else //if(A[i] == F[i])
                B[j*LEN+i].backDirection = ::up;

            //if (param->partId == 1)
            //    printf("%5.1f ", A[i]);

        }
        //if (param->partId == 1)
        //    printf("\n");
    }
    //printf("\n");
    //printf("\n");

    scoreMSA = A[nodeUp->length];

    //BACKWARD MOVING

    int* gapsUp   = new int[nodeUp->length+1];
    int* gapsLeft = new int[nodeLeft->length+1];
    
    for(int i=0; i<=nodeUp->length; i++)
        gapsUp[i] = 0;
    for(int i=0; i<=nodeLeft->length; i++)
        gapsLeft[i] = 0;

    int y = nodeLeft->length;
    int x = nodeUp->length;
    int carret = 0;
    BackDirection prev = crosswise;
    while(B[y*LEN + x].backDirection != stop)
    {
        if (prev == up && B[y*LEN + x].continueUp) //CONTINUE GOING UP
        {                                          //GAP EXTENSION
            gapsUp[x]++;
            carret++;
            y--;
        }
        else if (prev == ::left && B[y*LEN + x].continueLeft) //CONTINUE GOING LEFT
        {                                         //GAP EXTENSION
            gapsLeft[y]++;
            carret++;
            x--;
        }
        else
        {
            prev = B[y*LEN + x].backDirection;
            if(prev == up)
            {
                gapsUp[x]++;
                carret++;
                y--;
            }
            else if(prev == ::left)
            {
                gapsLeft[y]++;
                carret++;
                x--;
            }
            else //prev == crosswise
            {
                carret++;
                x--;
                y--;
            }
        }
    }

//    if(param->partId == 1)
//    {
//        printf("New values: gapOpen=%5.3f, gapExt=%5.3f\n", gapOp, gapEx);
//
//        for(int i = 0; i < nodeUp->length + 1; i++)
//        {
//            printf("GO=%5.3f, GE=%5.3f, GAPS=%d\n", gapOpsUp[i], gapExsUp[i], gapsUp[i]);
//        }
//    }

    // applying results to MSA array

    for (int seqU = 0; seqU < seqsUpCount; seqU++)
    {
        int gapNum = 0; // cumulative number of gaps
        int carret2 = 0;
        int start  = seqs->getStarts() [ seqsUp[seqU] ];
        int length = seqs->getLengths()[ seqsUp[seqU] ];

        for(int i=0; i<length; i++)
        {
            for(int j=carret2; j<MSA[start + i] + 1; j++)
            {
                gapNum += gapsUp[j];
            }
            carret2 = MSA[start + i] + 1;
            MSA[start + i] += gapNum;
        }
//        if(param->partId == 1)
//        {
//            for(int i = 0; i < nodeUp->length + 1; i++)
//            {
//                printf("%d ",gapsUp[i]);
//            }
//            printf("\n");
//            for(int i=0; i<length; i++)
//            {
//                printf("%d ", MSA[start + i]);
//            }
//            printf("\n");
//        }
    }

    for (int seqL = 0; seqL < seqsLeftCount; seqL++)
    {
        int gapNum = 0; // cumulative number of gaps
        int carret2 = 0;
        int start  = seqs->getStarts() [ seqsLeft[seqL] ];
        int length = seqs->getLengths()[ seqsLeft[seqL] ];

        for(int i=0; i<length; i++)
        {
            for(int j=carret2; j<MSA[start + i] + 1; j++)
            {
                gapNum += gapsLeft[j];
            }
            carret2 = MSA[start + i] + 1;
            MSA[start + i] += gapNum;
        }
//        if(param->partId == 1)
//        {
//            for(int i = 0; i < nodeLeft->length + 1; i++)
//            {
//                printf("%d ",gapsLeft[i]);
//            }
//            printf("\n");
//            for(int i=0; i<length; i++)
//            {
//                printf("%d ", MSA[start + i]);
//            }
//            printf("\n");
//        }
    }

//    // printing partial alignment
//    if(param->partId ==1)
//    {
//        for(int i=0; i<sequenceNumber; i++)
//        {
//            int carret = 0;
//            int start  = seqs->getStarts() [i];
//            int length = seqs->getLengths()[i];
//
//            printf("seq %3d ", i);
//
//            for(int j=0; j<length; j++)
//            {
//                for(int k=carret; k<MSA[start + j]; k++)
//                {
//                    printf("-");
//                }
//                printf("%c",
//                       seqs->getSubtitutionMatrix()->revConvert(seqs->getSequences()[start+length - j - 1])
//                        );
//                carret = MSA[start + j] + 1;
//            }
//
//            printf("\n");
//        }
//    }

    

    //printf("L=%d, R=%d, NEW=%d\n", nodeLeft->length, nodeUp->length, carret);        
    //param->node->length = max(nodeLeft->length, nodeUp->length);
    param->node->length = carret;
  

    //printf("%d\n", param->partId);
    delete[] SM;
    delete[] A;
    delete[] E;
    delete[] F;
    delete[] B;
    delete[] seqsLeft;
    delete[] seqsUp;
}

#undef PRI
#undef EXT
