/**
 * Author: Mark Larkin
 * 
 * Copyright (c) 2007 Des Higgins, Julie Thompson and Toby Gibson.  
 */
#ifdef HAVE_CONFIG_H
    #include "config.h"
#endif
#include <stdio.h>
#include <string>
#include <cstdlib>
#include <exception>
#include <iostream>
#include <sstream>
#include "SubMatrix.h"
#include "matrices.h"
#include "../general/InvalidCombination.cpp"
#include "../general/debuglogObject.h"

namespace clustalw
{
using namespace std;

/**
 * 
 * @param log 
 */
SubMatrix::SubMatrix()
: sizenAAMatrix(276),
  sizeDNAMatrix(153),
  matrixAvgScore(0),
  QTDNAHistMatNum(DNAIUB),
  QTAAHistMatNum(AAHISTGONNETPAM250),
  QTsegmentDNAMatNum(DNAIUB),
  QTsegmentAAMatNum(QTAASEGGONNETPAM250)
{
    userSeries = false;
         
    setUpCrossReferences();

    #if DEBUGFULL 
        if(logObject && DEBUGLOG)
        {
            logObject->logMsg("Creating the SubMatrix object\n");
        }
    #endif
    
    try
    {              
        /* Set up the vectors with the matrices defined in matrices.h 
         * The matrices are intially defined as a short array, as there is no way
         * to intiailize a vector with a {....} list.  
         */
        blosum30mtVec = new Matrix(blosum30mt, blosum30mt + sizenAAMatrix);
        blosum40mtVec = new Matrix(blosum40mt, blosum40mt + sizenAAMatrix);
        blosum45mtVec = new Matrix(blosum45mt, blosum45mt + sizenAAMatrix);
        blosum62mt2Vec = new Matrix(blosum62mt2, blosum62mt2 + sizenAAMatrix);
        blosum80mtVec = new Matrix(blosum80mt, blosum80mt + sizenAAMatrix);
        pam20mtVec = new Matrix(pam20mt, pam20mt + sizenAAMatrix);
        pam60mtVec = new Matrix(pam60mt, pam60mt + sizenAAMatrix);
        pam120mtVec = new Matrix(pam120mt, pam120mt + sizenAAMatrix);
        pam350mtVec = new Matrix(pam350mt, pam350mt + sizenAAMatrix);
        idmatVec = new Matrix(idmat, idmat + sizenAAMatrix);
        gon40mtVec = new Matrix(gon40mt, gon40mt + sizenAAMatrix);
        gon80mtVec = new Matrix(gon80mt, gon80mt + sizenAAMatrix);
        gon120mtVec = new Matrix(gon120mt, gon120mt + sizenAAMatrix);
        gon160mtVec = new Matrix(gon160mt, gon160mt + sizenAAMatrix);
        gon250mtVec = new Matrix(gon250mt, gon250mt + sizenAAMatrix);
        gon350mtVec = new Matrix(gon350mt, gon350mt + sizenAAMatrix);
        clustalvdnamtVec = new Matrix(clustalvdnamt, clustalvdnamt + sizeDNAMatrix);
        swgapdnamtVec = new Matrix(swgapdnamt, swgapdnamt + sizeDNAMatrix);
        
        /*
         * Set up the vectors for user defined types.
         * Probably dont need to do this, as they may not be used. It would be better
         * to initialise their size if thye are used.
         */
        userMat.resize(NUMRES * NUMRES);
        pwUserMat.resize(NUMRES * NUMRES);
        userDNAMat.resize(NUMRES * NUMRES);
        pwUserDNAMat.resize(NUMRES * NUMRES);
        QTscoreUserMatrix.resize(NUMRES * NUMRES);
        QTscoreUserDNAMatrix.resize(NUMRES * NUMRES);
        QTsegmentDNAMatrix.resize(NUMRES * NUMRES);
        QTsegmentAAMatrix.resize(NUMRES * NUMRES);
        
        userMatSeries.resize(MAXMAT); // Maximum number of matrices        
        vector<Matrix>::iterator firstM = userMatSeries.begin();
        vector<Matrix>::iterator lastM = userMatSeries.end();

        while(firstM != lastM)
        {
            firstM->resize(NUMRES * NUMRES);
            firstM++;
        }
        
        // Maybe I should put this in with the other xref intialisations!
        AAXrefseries.resize(MAXMAT);
        vector<Xref>::iterator firstX = AAXrefseries.begin();
        vector<Xref>::iterator lastX = AAXrefseries.end();
        
        while(firstX != lastX)
        {
            firstX->resize(NUMRES + 1);
            firstX++;
        }
        
        // Set the defaults.
        matrixNum = 3;
        matrixName = new string("gonnet");
        DNAMatrixNum = 1;
        DNAMatrixName = new string("iub");
        pwMatrixNum = 3;
        pwMatrixName = new string("gonnet");
        pwDNAMatrixNum = 1;
        pwDNAMatrixName = new string("iub");
    }
    catch(const exception &ex)
    {
        cerr << ex.what() << endl;
        cerr << "Terminating program. Cannot continue\n";
        exit(1);
    }  
}

void SubMatrix::setValuesToDefault()
{
    matrixAvgScore = 0;
    QTDNAHistMatNum = DNAIUB;
    QTAAHistMatNum = AAHISTGONNETPAM250;
    QTsegmentDNAMatNum = DNAIUB;
    QTsegmentAAMatNum = QTAASEGGONNETPAM250;
    userSeries = false;
    matrixNum = 3;
    DNAMatrixNum = 1;
    pwMatrixNum = 3;
    pwDNAMatrixNum = 1;      
}

/**
 * The destructor frees up any dynamically allocated memory.
 */
SubMatrix::~SubMatrix()
{
    delete blosum30mtVec;
    delete blosum40mtVec;
    delete blosum45mtVec;
    delete blosum62mt2Vec;
    delete blosum80mtVec;
    delete pam20mtVec;
    delete pam60mtVec;
    delete pam120mtVec;
    delete pam350mtVec;
    delete idmatVec;
    delete gon40mtVec;
    delete gon80mtVec;
    delete gon120mtVec;
    delete gon160mtVec;
    delete gon250mtVec;
    delete gon350mtVec;
    delete clustalvdnamtVec;
    delete swgapdnamtVec;
    delete matrixName;
    delete DNAMatrixName;
    delete pwMatrixName;
    delete pwDNAMatrixName;
}

/**
 * This function sets up the initial cross references. 
 */
void SubMatrix::setUpCrossReferences()
{
    char c1, c2;
    short i, j, maxRes;  
    defaultAAXref.resize(NUMRES + 1);
    defaultDNAXref.resize(NUMRES + 1);

    string aminoAcidOrder = "ABCDEFGHIKLMNPQRSTVWXYZ";
    string nucleicAcidOrder = "ABCDGHKMNRSTUVWXY"; 
    /* 
     * I also need to resize the user defined xrefs.
     */
    DNAXref.resize(NUMRES + 1);
    AAXref.resize(NUMRES + 1);
    pwAAXref.resize(NUMRES + 1);
    pwDNAXref.resize(NUMRES + 1);
    QTscoreXref.resize(NUMRES + 1);
    QTscoreDNAXref.resize(NUMRES + 1);
    QTsegmentDNAXref.resize(NUMRES + 1);
    QTsegmentAAXref.resize(NUMRES + 1);
    /*
     * set up cross-reference for default matrices hard-coded in matrices.h
     */
    for (i = 0; i < NUMRES; i++)
    {
        defaultAAXref[i] = -1;
    }
    for (i = 0; i < NUMRES; i++)
    {
        defaultDNAXref[i] = -1;
    }

    maxRes = 0;

    for (i = 0; (c1 = aminoAcidOrder[i]); i++)
    {
        for (j = 0; (c2 = userParameters->getAminoAcidCode(j)); j++)
        {
            if (c1 == c2)
            {
                defaultAAXref[i] = j;

                maxRes++;
                break;
            }
        }
        if ((defaultAAXref[i] ==  - 1) && (aminoAcidOrder[i] != '*'))
        {
            utilityObject->error("residue %c in matrices.h is not recognised",
                aminoAcidOrder[i]);
        }
    } 

    maxRes = 0;
    for (i = 0; (c1 = nucleicAcidOrder[i]); i++)
    {
        for (j = 0; (c2 = userParameters->getAminoAcidCode(j)); j++)
        {
            if (c1 == c2)
            {
                defaultDNAXref[i] = j;
                maxRes++;
                break;
            }
        }
        if ((defaultDNAXref[i] ==  - 1) && (nucleicAcidOrder[i] != '*'))
        {
            utilityObject->error("nucleic acid %c in matrices.h is not recognised",
                nucleicAcidOrder[i]);
        }
    }
    
}

/**
 * The function getPairwiseMatrix is called by the user to get the correct sub matrix for the
 * pairwise alignment stage. It calls getMatrix to actually calculate the matrix.
 * This function provides an interface for the user. It allows the user to get the matrix 
 * they wish to use for the pairwise alignment part. 
 * @param matrix[][] 
 * @param scale 
 * @param matAvg 
 * @return 
 */
int SubMatrix::getPairwiseMatrix(int matrix[NUMRES][NUMRES], PairScaleValues& scale,
                                 int& matAvg)
{
    int _maxRes; // Local copy.
    /* Pointers to Matrix and xref to use in calculation */
    Matrix* _matPtr;
    Xref* _matXref;
    
    #if DEBUGFULL 
        if(logObject && DEBUGLOG)
        {
            logObject->logMsg("In the function getPairwiseMatrix: \n");
        }    
    #endif
    
    string matrixPointer;
    string xrefPointer;
        
    #ifdef OS_MAC
        scale.intScale = 10;
    #else
        scale.intScale = 100;
    #endif
    scale.gapOpenScale = scale.gapExtendScale = 1.0;
    
    if (userParameters->getDNAFlag())
    {
        #if DEBUGFULL 
            if(logObject && DEBUGLOG)
            {
                string msg = "    (DNA AND Pairwise) " + *pwDNAMatrixName + "\n";
                logObject->logMsg(msg);
            }
        #endif
        if (pwDNAMatrixName->compare("iub") == 0)
        {
            matrixPointer = "swgapdnamtVec"; xrefPointer = "defaultDNAXref";
            _matPtr = swgapdnamtVec;
            _matXref = &defaultDNAXref;
        }
        else if (pwDNAMatrixName->compare("clustalw") == 0)
        {
            matrixPointer = "clustalvdnamtVec"; xrefPointer = "defaultDNAXref";
            _matPtr = clustalvdnamtVec;
            _matXref = &defaultDNAXref;
            scale.gapOpenScale = 0.6667;
            scale.gapExtendScale = 0.751;
        }
        else
        {
            matrixPointer = "pwUserDNAMat"; xrefPointer = "pwDNAXref";
            _matPtr = &pwUserDNAMat;
            _matXref = &pwDNAXref;
        }
        _maxRes = getMatrix(_matPtr, _matXref, matrix, true, scale.intScale);
        
        if (_maxRes == 0)
        {
            return ((int) -1);
        }
        float _transitionWeight = userParameters->getTransitionWeight();
        // Mark change 17-5-07
        matrix[0][4] = static_cast<int>(_transitionWeight * matrix[0][0]);
        matrix[4][0] = static_cast<int>(_transitionWeight * matrix[0][0]);
        matrix[2][11] = static_cast<int>(_transitionWeight * matrix[0][0]);
        matrix[11][2] = static_cast<int>(_transitionWeight * matrix[0][0]);
        matrix[2][12] = static_cast<int>(_transitionWeight * matrix[0][0]);
        matrix[12][2] = static_cast<int>(_transitionWeight * matrix[0][0]);
        
    }
    else
    {
        #if DEBUGFULL 
            if(logObject && DEBUGLOG)
            {
                string msg = "    (Protein AND Pairwise) " + *pwMatrixName + "\n";
                logObject->logMsg(msg);
            }
        #endif
        
        if (pwMatrixName->compare("blosum") == 0)
        {
            matrixPointer = "blosum30mtVec"; xrefPointer = "defaultAAXref";
            _matPtr = blosum30mtVec;
            _matXref = &defaultAAXref;
        }
        else if (pwMatrixName->compare("pam") == 0)
        {
            matrixPointer = "pam350mtVec"; xrefPointer = "defaultAAXref";
            _matPtr = pam350mtVec;
            _matXref = &defaultAAXref;
        }
        else if (pwMatrixName->compare("gonnet") == 0)
        {
            matrixPointer = "gon250mtVec"; xrefPointer = "defaultAAXref";
            _matPtr = gon250mtVec;
            scale.intScale /= 10;
            _matXref = &defaultAAXref;
        }
        else if (pwMatrixName->compare("id") == 0)
        {
            matrixPointer = "idmatVec"; xrefPointer = "defaultAAXref";
            _matPtr = idmatVec;
            _matXref = &defaultAAXref;
        }
        else
        {
            matrixPointer = "pwUserMat"; xrefPointer = "pwAAXref";
            _matPtr = &pwUserMat;
            _matXref = &pwAAXref;
        }

        _maxRes = getMatrix(_matPtr, _matXref, matrix, true, scale.intScale);
        
        if (_maxRes == 0)
        {
            return ((int) -1);
        }
    }
    
    #if DEBUGFULL 
        if(logObject && DEBUGLOG)
        {
            ostringstream outs;
            outs << "    Called getMatrix with "
                 << matrixPointer << " and " << xrefPointer << ".\n"
                 << "    intScale = " << scale.intScale << ", gapOpenScale = "
                 << scale.gapOpenScale << ", gapExtendScale = " << scale.gapExtendScale 
                 << "\n\n";
            logObject->logMsg(outs.str());
        }
    #endif
    matAvg = matrixAvgScore;  
    return _maxRes; 
}

/**
 * The function getProfileAlignMatrix provides an interface for the user to get
 * the matrix to be used in the profile alignment. This depends on the matrix series
 * that was chosen, and also on the percent identity.
 * @param matrix[][] 
 * @param pcid 
 * @param minLen 
 * @param scaleParam 
 * @param matAvg 
 * @return 
 */
int SubMatrix::getProfileAlignMatrix(int matrix[NUMRES][NUMRES], double pcid, int minLen, 
                                  PrfScaleValues& scaleParam, int& matAvg)
{
    bool found = false;
    bool errorGiven = false;
    bool _negMatrix = userParameters->getUseNegMatrix();
    int i = 0, j = 0;
    int _maxRes = 0;
    scaleParam.intScale = 100;
    string matrixPointer;
    string xrefPointer;
    
    #if DEBUGFULL 
        if(logObject && DEBUGLOG)
        {
            logObject->logMsg("In the function getProfileAlignMatrix: \n");
        }    
    #endif    
    if (userParameters->getDNAFlag())
    {
        #if DEBUGFULL 
            if(logObject && DEBUGLOG)
            {
                ostringstream outs;
                outs << "    (DNA AND Multiple align) "<< DNAMatrixName->c_str() << "\n";
                logObject->logMsg(outs.str());
            }    
        #endif        
        scaleParam.scale = 1.0;
        if (DNAMatrixName->compare("iub") == 0)
        {
            _matPtr = swgapdnamtVec;
            _matXref = &defaultDNAXref;
            matrixPointer = "swgapdnamtVec"; xrefPointer = "defaultDNAXref";
        }
        else if (DNAMatrixName->compare("clustalw") == 0)
        {
            _matPtr = clustalvdnamtVec;
            _matXref = &defaultDNAXref;
            scaleParam.scale = 0.66;
            matrixPointer = "clustalvdnamtVec"; xrefPointer = "defaultDNAXref";
        }
        else
        {
            _matPtr = &userDNAMat;
            _matXref = &DNAXref;
            matrixPointer = "userDNAMat"; xrefPointer = "DNAXref";
        }
        
        _maxRes = getMatrix(_matPtr, _matXref, matrix, _negMatrix,
                            static_cast<int>(scaleParam.intScale)); // Mark change 17-5-07
        if (_maxRes == 0)
        {
            return ((int) - 1);
        }
        
        float _transitionWeight = userParameters->getTransitionWeight();
        // fix suggested by Chanan Rubin at Compugen 
        matrix[(*_matXref)[0]][(*_matXref)[4]] = 
                            static_cast<int>(_transitionWeight * matrix[0][0]);
        matrix[(*_matXref)[4]][(*_matXref)[0]] = 
                            static_cast<int>(_transitionWeight * matrix[0][0]);
        matrix[(*_matXref)[2]][(*_matXref)[11]] = 
                            static_cast<int>(_transitionWeight * matrix[0][0]);
        matrix[(*_matXref)[11]][(*_matXref)[2]] = 
                            static_cast<int>(_transitionWeight * matrix[0][0]);
        matrix[(*_matXref)[2]][(*_matXref)[12]] = 
                            static_cast<int>(_transitionWeight * matrix[0][0]);
        matrix[(*_matXref)[12]][(*_matXref)[2]] = 
                            static_cast<int>(_transitionWeight * matrix[0][0]);

    }
    else // Amino acid alignment!!!!
    {   
        #if DEBUGFULL 
            if(logObject && DEBUGLOG)
            {
                ostringstream outs;
                outs << "    (Protein AND Multiple align) "<< matrixName->c_str() << "\n";
                logObject->logMsg(outs.str());
            }    
        #endif   
             
        scaleParam.scale = 0.75;
        if (matrixName->compare("blosum") == 0)
        {
            scaleParam.scale = 0.75;
            if (_negMatrix || userParameters->getDistanceTree() == false)
            {
                _matPtr = blosum40mtVec;
                matrixPointer = "blosum40mtVec";
            }
            else if (pcid > 80.0)
            {
                _matPtr = blosum80mtVec;
                matrixPointer = "blosum80mtVec";
            }
            else if (pcid > 60.0)
            {
                _matPtr = blosum62mt2Vec;
                matrixPointer = "blosum62mt2Vec";
            }
            else if (pcid > 40.0)
            {
                _matPtr = blosum45mtVec;
                matrixPointer = "blosum45mtVec";
            }
            else if (pcid > 30.0)
            {
                scaleParam.scale = 0.5;
                _matPtr = blosum45mtVec;
                matrixPointer = "blosum45mtVec";
            }
            else if (pcid > 20.0)
            {
                scaleParam.scale = 0.6;
                _matPtr = blosum45mtVec;
                matrixPointer = "blosum45mtVec";
            }
            else
            {
                scaleParam.scale = 0.6;
                _matPtr = blosum30mtVec;
                matrixPointer = "blosum30mtVec";
            }
            _matXref = &defaultAAXref;
            xrefPointer = "defaultAAXref";
        }
        else if (matrixName->compare("pam") == 0)
        {
            scaleParam.scale = 0.75;
            if (_negMatrix || userParameters->getDistanceTree() == false)
            {
                _matPtr = pam120mtVec;
                matrixPointer = "pam120mtVec";
            }
            else if (pcid > 80.0)
            {
                _matPtr = pam20mtVec;
                matrixPointer = "pam20mtVec";
            }
            else if (pcid > 60.0)
            {
                _matPtr = pam60mtVec;
                matrixPointer = "pam60mtVec";
            }
            else if (pcid > 40.0)
            {
                _matPtr = pam120mtVec;
                matrixPointer = "pam120mtVec";
            }
            else
            {
                _matPtr = pam350mtVec;
                matrixPointer = "pam350mtVec";
            }
            _matXref = &defaultAAXref;
            xrefPointer = "defaultAAXref";
        }
        else if (matrixName->compare("gonnet") == 0)
        {
            scaleParam.scale /= 2.0;
            if (_negMatrix || userParameters->getDistanceTree() == false)
            {
                _matPtr = gon250mtVec;
                matrixPointer = "gon250mtVec";
            }
            else if (pcid > 35.0)
            {
                _matPtr = gon80mtVec;
                scaleParam.scale /= 2.0;
                matrixPointer = "gon80mtVec";
            }
            else if (pcid > 25.0)
            {
                if (minLen < 100)
                {
                    _matPtr = gon250mtVec;
                    matrixPointer = "gon250mtVec";
                }
                else
                {
                    _matPtr = gon120mtVec;
                    matrixPointer = "gon120mtVec";
                }
            }
            else
            {
                if (minLen < 100)
                {
                    _matPtr = gon350mtVec;
                    matrixPointer = "gon350mtVec";
                }
                else
                {
                    _matPtr = gon160mtVec;
                    matrixPointer = "gon160mtVec";
                }
            }
            _matXref = &defaultAAXref;
            xrefPointer = "defaultAAXref";
            scaleParam.intScale /= 10;
        }
        else if (matrixName->compare("id") == 0)
        {
            _matPtr = idmatVec;
            _matXref = &defaultAAXref;
            xrefPointer = "defaultAAXref";
            matrixPointer = "idmatVec";
        }
        else if (userSeries)
        {
            _matPtr = NULL;
            found = false;
            for (i = 0; i < matSeries.nmat; i++)
            {
                if (pcid >= matSeries.mat[i].llimit && pcid <=
                    matSeries.mat[i].ulimit)
                {
                    j = i;
                    found = true;
                    break;
                }
            }
            if (found == false)
            {
                if (!errorGiven)
                {
                    utilityObject->warning(
                        "\nSeries matrix not found for sequence percent identity = %d.\n""(Using first matrix in series as a default.)\n""This alignment may not be optimal!\n""SUGGESTION: Check your matrix series input file and try again.", (int)pcid);
                }
                errorGiven = true;
                j = 0;
            }

            _matPtr = matSeries.mat[j].matptr;
            _matXref = matSeries.mat[j].AAXref;
            // this gives a scale of 0.5 for pcid=llimit and 1.0 for pcid=ulimit 
            scaleParam.scale = 0.5 + (pcid - matSeries.mat[j].llimit) / (
                (matSeries.mat[j].ulimit - matSeries.mat[j].llimit) *2.0);
            xrefPointer = "matSeries.mat[j].AAXref";
            matrixPointer = "matSeries.mat[j].matptr";
        }
        else
        {
            _matPtr = &userMat;
            _matXref = &AAXref;
            xrefPointer = "AAXref";
            matrixPointer = "userMat";
        }
        
        _maxRes = getMatrix(_matPtr, _matXref, matrix, _negMatrix,
                            static_cast<int>(scaleParam.intScale));
        if (_maxRes == 0)
        {
            cerr << "Error: matrix " << matrixName << " not found\n";
            return ( - 1);
        }       
    }
    
    #if DEBUGFULL 
        if(logObject && DEBUGLOG)
        {
            ostringstream outs;
            outs << "    Called getMatrix with "
                 << matrixPointer << " and " << xrefPointer << ".\n"
                 << "    intScale = " << scaleParam.intScale << ", scale = "
                 << scaleParam.scale << ", pcid = " << pcid << "\n\n";
            logObject->logMsg(outs.str());
        }    
    #endif    
    
    matAvg = matrixAvgScore;
    return _maxRes;
}

/**
 * The function getMatrix is what the other parts of the code call to get a useable 
 * substitution matrix. This is stored in matrix[NUMRES][NUMRES].
 * @param matptr 
 * @param xref 
 * @param matrix[][] 
 * @param negFlag 
 * @param scale 
 * @return 
 */
int SubMatrix::getMatrix(Matrix* matptr, Xref* xref, int matrix[NUMRES][NUMRES], 
                         bool negFlag, int scale, bool minimise)
{
    int ggScore = 0;
    int grScore = 0;
    int i, j, k, ix = 0;
    int ti, tj;
    int maxRes;
    int av1, av2, av3, min, max;  

    for (i = 0; i < NUMRES; i++)
    {
        for (j = 0; j < NUMRES; j++)
        {
            matrix[i][j] = 0;
        }
    }
    
    ix = 0;
    maxRes = 0;
    for (i = 0; i <= userParameters->getMaxAA(); i++)
    {
        ti = (*xref)[i];
        for (j = 0; j <= i; j++)
        {
            tj = (*xref)[j];
            if ((ti !=  - 1) && (tj !=  - 1))
            {
                k = (*matptr)[ix];
                if (ti == tj)
                {
                    matrix[ti][ti] = k * scale;
                    maxRes++;
                }
                else
                {
                    matrix[ti][tj] = k * scale;
                    matrix[tj][ti] = k * scale;
                }
                ix++;
            }
        }
    }

    --maxRes;

    av1 = av2 = av3 = 0;
    for (i = 0; i <= userParameters->getMaxAA(); i++)
    {
        for (j = 0; j <= i; j++)
        {
            av1 += matrix[i][j];
            if (i == j)
            {
                av2 += matrix[i][j];
            }
            else
            {
                av3 += matrix[i][j];
            }
        }
    }

    av1 /= (maxRes * maxRes) / 2;
    av2 /= maxRes;
    av3 = static_cast<int>(av3 / (((float)(maxRes * maxRes - maxRes)) / 2));
    matrixAvgScore =  - av3;

    min = max = matrix[0][0];
    for (i = 0; i <= userParameters->getMaxAA(); i++)
    for (j = 1; j <= i; j++)
    {
        if (matrix[i][j] < min)
        {
            min = matrix[i][j];
        }
        if (matrix[i][j] > max)
        {
            max = matrix[i][j];
        }
    }
    //cout << "MAX = " << max << "\n";
    if(!minimise)
    {
        /*
            if requested, make a positive matrix - add -(lowest score) to every entry
        */
        if (negFlag == false)
        {
            if (min < 0)
            {
                for (i = 0; i <= userParameters->getMaxAA(); i++)
                {
                    ti = (*xref)[i];
                    if (ti !=  - 1)
                    {
                        for (j = 0; j <= userParameters->getMaxAA(); j++)
                        {
                            tj = (*xref)[j];

                            if (tj !=  - 1)
                            {
                                matrix[ti][tj] -= min;
                            }
                        }
                    }
                }
            }
        }
    
        // local copies of the gap positions
        int _gapPos1 = userParameters->getGapPos1();
        int _gapPos2 = userParameters->getGapPos2();

        for (i = 0; i < userParameters->getGapPos1(); i++)
        {
            matrix[i][_gapPos1] = grScore;
            matrix[_gapPos1][i] = grScore;
            matrix[i][_gapPos2] = grScore;
            matrix[_gapPos2][i] = grScore;
        }
        matrix[_gapPos1][_gapPos1] = ggScore;
        matrix[_gapPos2][_gapPos2] = ggScore;
        matrix[_gapPos2][_gapPos1] = ggScore;
        matrix[_gapPos1][_gapPos2] = ggScore;
    }
    else
    {
        // DO THE SAGA MATRIX
        for (i = 0; i <= userParameters->getMaxAA(); i++)
        {
            for (j = 0; j <= userParameters->getMaxAA(); j++)
            {
                matrix[i][j] = max - matrix[i][j];
            }
        }    
    }
    maxRes += 2;

    return (maxRes);
}

/**
 * The function getUserMatFromFile is used to read in a user defined matrix.
 * @param str 
 * @param alignResidueType 
 * @param alignType 
 * @return 
 */
bool SubMatrix::getUserMatFromFile(char *str, int alignResidueType, int alignType)
{
    int maxRes;

    FILE *infile;
    // Need to check if the values are a valid combination!
    checkResidueAndAlignType(alignResidueType, alignType);
     
    if(userParameters->getMenuFlag())
    {
        utilityObject->getStr(string("Enter name of the matrix file"), line2);
    }
    else
    {
        line2 = string(str);
    }

    if(line2.size() == 0)
    {
        return false;
    }

    if((infile = fopen(line2.c_str(), "r")) == NULL) 
    {
        utilityObject->error("Cannot find matrix file [%s]", line2.c_str());
        return false;
    }

    strcpy(str, line2.c_str());
    // Find out which part of the code we are reading the matrix in for.
    mat = getUserMatAddress(alignResidueType, alignType);
    xref = getUserXrefAddress(alignResidueType, alignType);
    
    if ((alignResidueType == Protein) && (alignType == MultipleAlign))
    {
        // Try read in a matrix series.
        maxRes = readMatrixSeries(str, userMat, AAXref);
    }
    else // Read in a single matrix!!!
    {
        maxRes = readUserMatrix(str, *mat, *xref);
    }
    
    if (maxRes <= 0) return false;

    return true;
}

/**
 * The function compareMatrices is used to compare 2 matrices that have been read in.
 * It will compare them in a given region. It will not compare all of them, as some of it
 * will be random memory.
 * @param mat1[][] 
 * @param mat2[][] 
 */
void SubMatrix::compareMatrices(int mat1[NUMRES][NUMRES], int mat2[NUMRES][NUMRES])
{
    int same = 1;
    for(int row = 0; row < NUMRES; row++)
    {
        for(int col = 0; col < NUMRES; col++)
        {
            if(mat1[row][col] != mat2[row][col])
            {
                same = 0;
                cout << "The row is " << row << ". The column is " << col << endl;
                break; // It is not the same. End the loop.
            }
        }
    }
    
    if(same == 0)
    {
        cout << "It was not the same\n";
    }
    else
    {
        cout << "It is the same\n";
    }
}

/**
 * This function is simply to display the results of the getMatrix function.
 * This is so that I can compare it to the original clustal version of it.
 * @param mat[][] 
 */
void SubMatrix::printGetMatrixResults(int mat[NUMRES][NUMRES])
{
    ofstream outfile("getmatrix.out");
    
    if(!outfile)
        cerr<<"oops failed to open !!!\n";
    
    for(int row = 0; row < NUMRES; row++)
    {
        for(int col = 0; col < NUMRES; col++)
        {
            if((mat[row][col] > 9) || (mat[row][col] < 0))
            {
                outfile <<" "<< mat[row][col]<<",";
            }
            else
            {
                outfile <<"  "<< mat[row][col]<<",";
            }
        }
        outfile<<"\n";
    }
}

/**
 * This function is from the old interface. It simply calls the readMatrixSeries
 * function. It seems to only be called from the commandline part.
 * @param str 
 * @return 
 */
bool SubMatrix::getUserMatSeriesFromFile(char *str)
{
    int maxRes;

    FILE *infile;

    if(userParameters->getMenuFlag()) // Check if we are using menu!
    {
        utilityObject->getStr(string("Enter name of the matrix file"), line2);
    }
    else
    {
        //strcpy(lin2,str);
        line2 = string(str);
    }
    
    if(line2.size() == 0) return false;

    if((infile = fopen(line2.c_str(), "r"))==NULL) 
    {
        utilityObject->error("Cannot find matrix file [%s]",line2.c_str());
        return false;
    }

    strcpy(str, line2.c_str());

    maxRes = readMatrixSeries(str, userMat, AAXref);
    if (maxRes <= 0) return false;

    return true;
}

/*
 * The function readMatrixSeries is used to read in a series of matrices from a file.
 * It calls readUserMatrix to read in the individual matrices. The matrices are stored
 * in userMatSeries. 
 */
int SubMatrix::readMatrixSeries(const char *fileName, Matrix& userMat, Xref& xref)
{
    FILE *fd = NULL;
    char mat_fileName[FILENAMELEN];
    char inline1[1024];
    int maxRes = 0;
    int nmat;
    int n, llimit, ulimit;

    if (fileName[0] == '\0')
    {
        utilityObject->error("comparison matrix not specified");
        return ((int)0);
    }
    if ((fd = fopen(fileName, "r")) == NULL)
    {
        utilityObject->error("cannot open %s", fileName);
        return ((int)0);
    }

    /* check the first line to see if it's a series or a single matrix */
    while (fgets(inline1, 1024, fd) != NULL)
    {
        if (commentline(inline1))
        {
            continue;
        }
        if (utilityObject->lineType(inline1, "CLUSTAL_SERIES"))
        {
            userSeries = true;
        }
        else
        {
            userSeries = false;
        }
        break;
    }

    /* it's a single matrix */
    if (userSeries == false)
    {
        fclose(fd);
        maxRes = readUserMatrix(fileName, userMat, xref);
        return (maxRes);
    }

    /* it's a series of matrices, find the next MATRIX line */
    nmat = 0;
    matSeries.nmat = 0;
    while (fgets(inline1, 1024, fd) != NULL)
    {
        if (commentline(inline1))
        {
            continue;
        }
        if (utilityObject->lineType(inline1, "MATRIX")) // Have found a matrix
        {
            if (sscanf(inline1 + 6, "%d %d %s", &llimit, &ulimit, mat_fileName)
                != 3)
            {
                utilityObject->error("Bad format in file %s\n", fileName);
                fclose(fd);
                return ((int)0);
            }
            if (llimit < 0 || llimit > 100 || ulimit < 0 || ulimit > 100)
            {
                utilityObject->error("Bad format in file %s\n", fileName);
                fclose(fd);
                return ((int)0);
            }
            if (ulimit <= llimit)
            {
                utilityObject->error("in file %s: lower limit is greater than upper (%d-%d)\n",
                    fileName, llimit, ulimit);
                fclose(fd);
                return ((int)0);
            }
            
            n = readUserMatrix(mat_fileName, userMatSeries[nmat],
                AAXrefseries[nmat]);
            
            //cout << "Read in matrix number " << nmat << "\n"; // NOTE Testing!!!!! 
            char nameOfFile[] = "matrix"; 
            printInFormat(userMatSeries[nmat], nameOfFile);
            if (n <= 0)
            {
                utilityObject->error("Bad format in matrix file %s\n", mat_fileName);
                fclose(fd);
                return ((int)0);
            }
            matSeries.mat[nmat].llimit = llimit;
            matSeries.mat[nmat].ulimit = ulimit;
            matSeries.mat[nmat].matptr = &userMatSeries[nmat];
            matSeries.mat[nmat].AAXref = &AAXrefseries[nmat];
            nmat++;
            
            if(nmat >= MAXMAT)
            {
               // We have read in all the matrices that we can read into this vector
               // Write a message to the screen, and break out of loop.
               cerr << "The matrix series file has more entries than allowed in \n"
                    << "a user defined series. The most that are allowed is "
                    << MAXMAT << ".\n"
                    << "The first " << MAXMAT << " have been read in and will be used.\n";
               break; // Get out of the loop! 
            }
        }
    }
    fclose(fd);
    matSeries.nmat = nmat;

    maxRes = n;
    return (maxRes);
}

/*
 * This function is used to read a single user matrix from a file.
 * It can be called repeatedly if there are multiple matrices in the file.
 */
int SubMatrix::readUserMatrix(const char *fileName, Matrix& userMat, Xref& xref)
{
    double f;
    FILE *fd;
    int numargs, farg;
    int i, j, k = 0;
    char codes[NUMRES];
    char inline1[1024];
    char *args[NUMRES + 4];
    char c1, c2;
    int ix1, ix = 0;
    int maxRes = 0;
    float scale;

    if (fileName[0] == '\0')
    {
        utilityObject->error("comparison matrix not specified");
        return ((int)0);
    }

    if ((fd = fopen(fileName, "r")) == NULL)
    {
        utilityObject->error("cannot open %s", fileName);
        return ((int)0);
    }
    maxRes = 0;
    while (fgets(inline1, 1024, fd) != NULL)
    {
        if (commentline(inline1))
        {
            continue;
        }
        if (utilityObject->lineType(inline1, "CLUSTAL_SERIES"))
        {
            utilityObject->error("in %s - single matrix expected.", fileName);
            fclose(fd);
            return ((int)0);
        }
        /*
        read residue characters.
         */
        k = 0;
        for (j = 0; j < (int)strlen(inline1); j++)
        {
            if (isalpha((int)inline1[j]))
            {
                codes[k++] = inline1[j];
            }
            if (k > NUMRES)
            {
                utilityObject->error("too many entries in matrix %s", fileName);
                fclose(fd);
                return ((int)0);
            }
        }
        codes[k] = '\0';
        break;
    }

    if (k == 0)
    {
        utilityObject->error("wrong format in matrix %s", fileName);
        fclose(fd);
        return ((int)0);
    }

    /*
    cross-reference the residues
     */
    for (i = 0; i < NUMRES; i++)
    {
        xref[i] =  - 1;
    }

    maxRes = 0;
    for (i = 0; (c1 = codes[i]); i++)
    {
        for (j = 0; (c2 = userParameters->getAminoAcidCode(j)); j++)
        if (c1 == c2)
        {
            xref[i] = j;
            maxRes++;
            break;
        }
        if ((xref[i] ==  - 1) && (codes[i] != '*'))
        {
            utilityObject->warning("residue %c in matrix %s not recognised", codes[i],
                fileName);
        }
    }


    /*
    get the weights
     */

    ix = ix1 = 0;
    while (fgets(inline1, 1024, fd) != NULL)
    {
        if (inline1[0] == '\n')
        {
            continue;
        }
        if (inline1[0] == '#' || inline1[0] == '!')
        {
            break;
        }
        numargs = getArgs(inline1, args, (int)(k + 1));
        if (numargs < maxRes)
        {
            utilityObject->error("wrong format in matrix %s", fileName);
            fclose(fd);
            return ((int)0);
        }
        if (isalpha(args[0][0]))
        {
            farg = 1;
        }
        else
        {
            farg = 0;
        }

        /* decide whether the matrix values are float or decimal */
        scale = 1.0;
        for (i = 0; i < (int)strlen(args[farg]); i++)
        if (args[farg][i] == '.')
        {
            /* we've found a float value */
            scale = 10.0;
            break;
        }

        for (i = 0; i <= ix; i++)
        {
            if (xref[i] !=  - 1)
            {
                f = atof(args[i + farg]);
                userMat[ix1++] = (short)(f *scale);
            }
        }
        ix++;
    }
    if (ix != k + 1)
    {
        utilityObject->error("wrong format in matrix %s", fileName);
        fclose(fd);
        return ((int)0);
    }

    userMat.resize(ix1 + 1);
    
    maxRes += 2;
    fclose(fd);

    return (maxRes);
}

/**
 * 
 * @param inline1 
 * @param args[] 
 * @param max 
 * @return 
 */
int SubMatrix::getArgs(char *inline1,char *args[],int max)
{
    char *inptr;
    int i;

    inptr = inline1;
    for (i = 0; i <= max; i++)
    {
        if ((args[i] = strtok(inptr, " \t\n")) == NULL)
        {
            break;
        }
        inptr = NULL;
    }

    return (i);
}


/**
 * 
 * @return 
 */
int SubMatrix::getMatrixNum()
{
    return matrixNum;
}

/**
 * 
 * @return 
 */
int SubMatrix::getDNAMatrixNum()
{
    return DNAMatrixNum;
}

/**
 * 
 * @return 
 */
int SubMatrix::getPWMatrixNum()
{
    return pwMatrixNum;
}

/**
 * 
 * @return 
 */
int SubMatrix::getPWDNAMatrixNum()
{
    return pwDNAMatrixNum;
}

/**
 * The function setCurrentNameAndNum is used to select a matrix series.
 * This will then be used for the alignment. The matrices will change, but the
 * series remains the same. We can set the series for pairwise/full for both 
 * protein and DNA. NOTE: The default can be set to user defined matrices.
 * @param _matrixName 
 * @param _matrixNum 
 * @param alignResidueType 
 * @param alignType 
 */
void SubMatrix::setCurrentNameAndNum(string _matrixName, int _matrixNum, 
                                          int alignResidueType,int alignType)
{       
    // Check if the values are valid.
    checkResidueAndAlignType(alignResidueType, alignType);
    
    string residue;
    string align;
    if((alignResidueType == Protein) && (alignType == Pairwise))
    {
        residue = "Protein"; align = "Pairwise";
        pwMatrixNum = _matrixNum;
        pwMatrixName = new string(_matrixName);
    }
    else if((alignResidueType == Protein) && (alignType == MultipleAlign))
    {
        residue = "Protein"; align = "MultipleAlign";
        matrixNum = _matrixNum;
        matrixName = new string(_matrixName); 
    }
    else if((alignResidueType == DNA) && (alignType == Pairwise))
    {
        residue = "DNA"; align = "Pairwise";
        pwDNAMatrixNum = _matrixNum;
        pwDNAMatrixName = new string(_matrixName);     
    }
    else if((alignResidueType == DNA) && (alignType == MultipleAlign))
    {
        residue = "DNA"; align = "MultipleAlign";    
        DNAMatrixNum = _matrixNum;
        DNAMatrixName = new string(_matrixName);
            
    }
    
    #if DEBUGFULL 
        if(logObject && DEBUGLOG)
        {
            ostringstream outs;
            outs << "The matrix/matrix series has been changed for "
                 << "(" << residue << " AND " << align << ")."
                 << " New value: " << _matrixName << "\n\n";           
            logObject->logMsg(outs.str());
        }
    #endif    
}

/**
 * 
 * @param alignResidueType 
 * @param alignType 
 * @return 
 */
int SubMatrix::getMatrixNumForMenu(int alignResidueType, int alignType)
{
    checkResidueAndAlignType(alignResidueType, alignType);
    
    if((alignResidueType == Protein) && (alignType == Pairwise))
    {
        return pwMatrixNum;
    }
    else if((alignResidueType == Protein) && (alignType == MultipleAlign))
    {
        return matrixNum;  
    }
    else if((alignResidueType == DNA) && (alignType == Pairwise))
    {
        return pwDNAMatrixNum;   
    }
    else if((alignResidueType == DNA) && (alignType == MultipleAlign))
    {
        return DNAMatrixNum;   
    }
    else
        return -100; // NOTE NONE of these. I need to put in better error checking
}

/**
 * 
 * @param line 
 * @return 
 */
bool SubMatrix::commentline(char* line)
{
    int i;

    if (line[0] == '#')
    {
        return true;
    }
    for (i = 0; line[i] != '\n' && line[i] != EOS; i++)
    {
        if (!isspace(line[i]))
        {
            return false;
        }
    }
    return true;
}

/**
 * This function prints out the vector to the file specified by name. 
 * The vector is printed out in a triangular format, the same as the way
 * the arrays are displayed in matrices.h. This function is for testing purposes.
 * @param temp 
 * @param name 
 */
void SubMatrix::printInFormat(vector<short>& temp, char* name)
{
    char nameOfFile[30];
    strcpy(nameOfFile, name);
    strcat(nameOfFile, ".out");

    ofstream outfile(nameOfFile);
    
    if(!outfile)
        cerr<<"oops failed to open !!!\n";
    
    outfile<<"short "<<name<<"[]{\n";
            
    int numOnCurrentLine = 1;
    int soFar = 0;
    for(int i = 0; i < (int)temp.size(); i++)
    {
        if(soFar == numOnCurrentLine)
        {
            outfile<<"\n";
            soFar = 0;
            numOnCurrentLine++;            
        }
        if((temp[i] > 9) || (temp[i] < 0))
        {
            outfile <<" "<< temp[i]<<",";
        }
        else
        {
            outfile <<"  "<< temp[i]<<",";
        }
        
        // Now increment so far
        soFar++;
        // Check to see if the next element is the last element
        if((i + 1) == (int)temp.size() - 1)
        {
            // Print out the last element. Then a curly brace, and break.
            if((temp[i+1] > 9) || (temp[i+1] < 0))
            {
                outfile <<" "<< temp[i + 1]<<"};\n";
            }
            else
            {
                outfile <<"  "<< temp[i + 1]<<"};\n";
            }
            break;
        }  
    }
    
    ofstream outfile2("temp.out");
    for(int i = 0; i < (int)temp.size(); i++)
    {
        outfile2 << temp[i] << " ";
    }   
}


/**
 * 
 * @param temp 
 * @param name 
 */
void SubMatrix::printVectorToFile(vector<short>& temp, char* name)
{
    char nameOfFile[30];
    strcpy(nameOfFile, name);
    strcat(nameOfFile, ".out");

    ofstream outfile(nameOfFile);
    
    if(!outfile)
        cerr<<"oops failed to open !!!\n";
    
    for(int i = 0; i < (int)temp.size(); i++)
    {
        if((temp[i] > 9) || (temp[i] < 0))
        {
            outfile <<" "<< temp[i]<<",";
        }
        else
        {
            outfile <<"  "<< temp[i]<<",";
        }
    }
    outfile.close();
}


/**
 * 
 * @param alignResidueType 
 * @param alignType 
 * @return 
 */
Matrix* SubMatrix::getUserMatAddress(int alignResidueType, int alignType)
{
    if((alignResidueType == Protein) && (alignType == Pairwise))
    {
        return &pwUserMat;
    }
    else if((alignResidueType == Protein) && (alignType == MultipleAlign))
    {
        return &userMat; 
    }
    else if((alignResidueType == DNA) && (alignType == Pairwise))
    {
        return &pwUserDNAMat; 
    }
    else if((alignResidueType == DNA) && (alignType == MultipleAlign))
    {
        return &userDNAMat;  
    }
    return NULL;
}


/**
 * 
 * @param alignResidueType 
 * @param alignType 
 * @return 
 */
Xref* SubMatrix::getUserXrefAddress(int alignResidueType, int alignType)
{
    if((alignResidueType == Protein) && (alignType == Pairwise))
    {
        return &pwAAXref;
    }
    else if((alignResidueType == Protein) && (alignType == MultipleAlign))
    {
        return &AAXref; 
    }
    else if((alignResidueType == DNA) && (alignType == Pairwise))
    {
        return &pwDNAXref;  
    }
    else if((alignResidueType == DNA) && (alignType == MultipleAlign))
    {
        return &DNAXref;   
    }
    return NULL;
}

/**
 * This is an error handling routine. If an incorrect combination of values
 * is given, it will terminate the program.
 * @param alignResidueType 
 * @param alignType 
 */
void SubMatrix::checkResidueAndAlignType(int alignResidueType, int alignType)
{
    if(((alignResidueType != 0) && (alignResidueType != 1))
           || ((alignType != 0) && (alignType != 1)))
    {
        InvalidCombination ex(alignResidueType, alignType);
        ex.whatHappened();
        exit(1);
    }
}

/**
 * The function tempInterface is used to call the SubMatrix in the way it is
 * supposed to be used. It is for testing purposes.
 * @param alignResidueType 
 * @param alignType 
 */
void SubMatrix::tempInterface(int alignResidueType, int alignType)
{
    char userFile[FILENAMELEN + 1];

    userParameters->setDNAFlag(true);
    strcpy(userFile, "mat1");
    userParameters->setMenuFlag(false);   
    getUserMatFromFile(userFile, DNA, Pairwise);
    setCurrentNameAndNum(userFile, 4, 3, Pairwise);
    
    setCurrentNameAndNum("gonnet", 4, Protein, Pairwise);
}

/**
 * A single matrix is used in scoring the alignment. This is Blosum45. This is the 
 * function to get it.
 * @param matrix[][] 
 * @return 
 */
int SubMatrix::getAlnScoreMatrix(int matrix[NUMRES][NUMRES])
{
    int _maxNumRes;
    /* 
       //_maxNumRes = getMatrix(blosum45mtVec, &defaultAAXref, matrix, true, 100);
       _maxNumRes = getMatrix(blosum45mtVec, &defaultAAXref, matrix, false, 1, true);
       //_maxNumRes = getMatrix(blosum62mt2Vec, &defaultAAXref, matrix, true, 100);
    */
    // 1.83 style
    _maxNumRes = getMatrix(blosum45mtVec, &defaultAAXref, matrix, true, 100);

    return _maxNumRes;
}

/**
 * This function is used to get the matrix that will be used for the calculation of 
 * the histogram. The histogram values are used in ClustalQt. 
 * @param matrix[][] 
 * @param matNum 
 * @param dnaMatNum 
 */
void SubMatrix::getQTMatrixForHistogram(int matrix[NUMRES][NUMRES])
{
    Matrix* _matPtrLocal;
    Xref* _matXrefLocal;
    int maxRes;
    if(userParameters->getDNAFlag())
    {
        if (QTDNAHistMatNum == DNAUSERDEFINED)
        {
            _matPtrLocal = &QTscoreUserDNAMatrix;
            _matXrefLocal = &QTscoreDNAXref;
        }
        else if (QTDNAHistMatNum == DNACLUSTALW)
        {
            _matPtrLocal = clustalvdnamtVec;
            _matXrefLocal = &defaultDNAXref;
        }
        else
        {
            _matPtrLocal = swgapdnamtVec;
            _matXrefLocal = &defaultDNAXref;            
        }
    }
    else 
    {
        if (QTAAHistMatNum == AAHISTIDENTITY)
        {
            _matPtrLocal = idmatVec;
            _matXrefLocal = &defaultAAXref;
        }
        else if (QTAAHistMatNum == AAHISTGONNETPAM80)
        {
            _matPtrLocal = gon80mtVec;
            _matXrefLocal = &defaultAAXref;
        }
        else if (QTAAHistMatNum == AAHISTGONNETPAM120)
        {
            _matPtrLocal = gon120mtVec;
            _matXrefLocal = &defaultAAXref;
        }
        else if (QTAAHistMatNum == AAHISTUSER)
        {
            _matPtrLocal = &QTscoreUserMatrix;
            _matXrefLocal = &QTscoreXref;
        }
        else if (QTAAHistMatNum == AAHISTGONNETPAM350)
        {
            _matPtrLocal = gon350mtVec;
            _matXrefLocal = &defaultAAXref;
        }
        else // Default
        {
            _matPtrLocal = gon250mtVec;
            _matXrefLocal = &defaultAAXref;            
        }
    }
    maxRes = getMatrix(_matPtrLocal, _matXrefLocal, matrix, false, 100);

}

void SubMatrix::getQTMatrixForLowScoreSeg(int matrix[NUMRES][NUMRES])
{
    Matrix* _matPtrLocal;
    Xref* _matXrefLocal;
    int maxRes;
    int _maxAA = userParameters->getMaxAA();
    int max = 0;
    int offset;
        
    if(userParameters->getDNAFlag())
    {
        if (QTsegmentDNAMatNum == DNAUSERDEFINED)
        {
            _matPtrLocal = &QTsegmentDNAMatrix;
            _matXrefLocal = &QTsegmentDNAXref;            
        }
        else if (QTsegmentDNAMatNum == DNACLUSTALW)
        {
            _matPtrLocal = clustalvdnamtVec;
            _matXrefLocal = &defaultDNAXref;
        }
        else
        {
            _matPtrLocal = swgapdnamtVec;
            _matXrefLocal = &defaultDNAXref;
        }
        /* get a positive matrix - then adjust it according to scale */
        maxRes = getMatrix(_matPtrLocal, _matXrefLocal, matrix, false, 100);
        /* find the maximum value */
        for(int i = 0; i <= _maxAA; i++)
        {
            for(int j = 0; j <= _maxAA; j++)
            {
                if(matrix[i][j] > max)
                { 
                    max = matrix[i][j];
                }
            }
        }
        /* subtract max * scale / 2 from each matrix value */
        offset = static_cast<int>(static_cast<float>(max *
                                  userParameters->getQTlowScoreDNAMarkingScale()) / 20.0);

        for(int i = 0; i <= _maxAA; i++)
        {
            for(int j = 0; j <= _maxAA; j++)
            {
                matrix[i][j] -= offset;
            }
        }
    }
    else
    {
        if (QTsegmentAAMatNum == QTAASEGGONNETPAM80)
        {
            _matPtrLocal = gon80mtVec;
            _matXrefLocal = &defaultAAXref;
        }
        else if (QTsegmentAAMatNum == QTAASEGGONNETPAM120)
        {
            _matPtrLocal = gon120mtVec;
            _matXrefLocal = &defaultAAXref;
        }
        else if (QTsegmentAAMatNum == QTAASEGUSER)
        {
            _matPtrLocal = &QTsegmentAAMatrix;
            _matXrefLocal = &QTsegmentAAXref;
        }
        else if (QTsegmentAAMatNum == QTAASEGGONNETPAM350)
        {
            _matPtrLocal = gon350mtVec;
            _matXrefLocal = &defaultAAXref;
        }
        else
        {
            _matPtrLocal = gon250mtVec;
            _matXrefLocal = &defaultAAXref;            
        }
        /* get a negative matrix */
        maxRes = getMatrix(_matPtrLocal, _matXrefLocal, matrix, true, 100);
    }
}

bool SubMatrix::getQTLowScoreMatFromFile(char* fileName, bool dna)
{
    int maxRes;

    FILE *infile;
     
    line2 = string(fileName);

    if(line2.size() == 0)
    {
        return false;
    }

    if((infile = fopen(line2.c_str(), "r")) == NULL) 
    {
        utilityObject->error("Cannot find matrix file [%s]", line2.c_str());
        return false;
    }

    strcpy(fileName, line2.c_str());
    
    if(dna)
    {
        maxRes = readUserMatrix(fileName, QTsegmentDNAMatrix, QTsegmentDNAXref);
    }
    else
    {
        maxRes = readUserMatrix(fileName, QTsegmentAAMatrix, QTsegmentAAXref);
    }
    
    if (maxRes <= 0) 
    {
        return false;
    }

    return true;
}

bool SubMatrix::getAAScoreMatFromFile(char *str)
{
    int maxRes;

    FILE *infile;
     
    line2 = string(str);

    if(line2.size() == 0)
    {
        return false;
    }

    if((infile = fopen(line2.c_str(), "r")) == NULL) 
    {
        utilityObject->error("Cannot find matrix file [%s]", line2.c_str());
        return false;
    }

    strcpy(str, line2.c_str());
    
    maxRes = readUserMatrix(str, QTscoreUserMatrix, QTscoreXref);
    
    if (maxRes <= 0) 
    {
        return false;
    }

    return true;
}

bool SubMatrix::getDNAScoreMatFromFile(char *str)
{
    int maxRes;

    FILE *infile;
     
    line2 = string(str);

    if(line2.size() == 0)
    {
        return false;
    }

    if((infile = fopen(line2.c_str(), "r")) == NULL) 
    {
        utilityObject->error("Cannot find matrix file [%s]", line2.c_str());
        return false;
    }

    strcpy(str, line2.c_str());
    
    maxRes = readUserMatrix(str, QTscoreUserDNAMatrix, QTscoreDNAXref);
    
    if (maxRes <= 0) 
    {
        return false;
    }

    return true;
}

}