/**
 * Author: Mark Larkin
 * 
 * Copyright (c) 2007 Des Higgins, Julie Thompson and Toby Gibson.  
 */
/**
 * The SymMatrix class is used to store the distance matrix. It stores it in a double  
 * array.
 * This class throws an out_of_range exception if we try to access an element outside 
 * of the array bounds. It will throw a bad_alloc if we cannot allocate enough memory for 
 * the distance matrix.
 */
#ifndef SYMMATRIX_H
#define SYMMATRIX_H

#include<iostream>
#include <vector>
#include <iomanip>
#include <cmath>
#include <stdexcept>
#include <new>
#include <stdlib.h>

namespace clustalw
{
using namespace std;

class SymMatrix
{
    public:
        SymMatrix() : elements(0), numSeqs(0), subElements(0), firstSeq(0), numSeqsInSub(0){;}
        SymMatrix(int size)
         : elements(0),
           numSeqs(0),
           subElements(0),
           firstSeq(0),
           numSeqsInSub(0) 
        {
            // Size is the numSeqs + 1
            numSeqs = size - 1;
            sizeElements = ((numSeqs + 1) * (numSeqs + 2)) >> 1;
            try
            {
                elements = new double[sizeElements];
                setAllArrayToVal(elements, sizeElements, 0.0);
            }
            catch(bad_alloc& e)
            {
                cout << "Could not allocate a distance matrix for " << numSeqs 
                     << " seqs.\n";
                throw e;
            }
            
        }
        ~SymMatrix()
        {
            if(elements)
            {
                delete [] elements;
            }
            if(subElements)
            {
                delete [] subElements;
            }
        }
        
        inline void SetAt(int nRow, int nCol, const double& value)
        {
            index = getIndex(nRow, nCol, numSeqs);
            elements[index] = value; 
        }
        
        inline double GetAt(int nRow, int nCol) 
        {
            index = getIndex(nRow, nCol, numSeqs);
            return elements[index];
        }
   
        inline void ResizeRect(int size, double val = 0.0)
        {
            numSeqs = size - 1;
            sizeElements = ((numSeqs + 1) * (numSeqs + 2)) >> 1;
            if(elements)
            {
                delete [] elements;
            }
            try
            {
                elements = new double[sizeElements];
                setAllArrayToVal(elements, sizeElements, val);
            }
            catch(bad_alloc& e)
            {
                cout << "Could not allocate a distance matrix for " << numSeqs 
                     << " seqs. Need to terminate program.\n";
                throw e;
            }            
        }
        
        inline void setAllArrayToVal(double* array, int size, double val)
        {
            for(int i = 0; i < size; i++)
            {
                array[i] = val;
            }
        } 
        
        inline int getSize()
        {
            return numSeqs + 1;
        }

        inline void clearArray()
        {
            sizeElements = 0;
            numSeqs = 0;
            if(elements)
            {
                delete [] elements;
                elements = 0;
            }
            
            numSeqsInSub = 0;
            sizeSubElements = 0;
            
            if(subElements)
            {
                delete [] subElements;
                subElements = 0;                
            }
        }
         
        void printArray()
        {
            printArray(elements, sizeElements);
        }
        
        void printSubArray()
        {
            printArray(subElements, sizeSubElements);
        }
        
        void printArray(double* array, int size)
        {
            int numThisTime = 1;
            int numSoFar = 0;
            for(int i = 0; i < size; i++)
            {
                
                numSoFar++;
                
                if(numSoFar == numThisTime)
                {
                    cout << " " << setw(4) << array[i] << "\n";
                    numSoFar = 0;
                    numThisTime++;
                }
                else
                {
                    cout << " " << setw(4) << array[i];
                }    
            }
            cout << "\n";
        }
        
        inline int getIndex(const int &i, const int &j, const int &nSeqs) const
        {
            if(i == 0 || j == 0)
            {
                return 0;
            }
            
            int _i = i - 1, _j = j - 1;
            if (_i == _j) 
            {
                if ((_i >= nSeqs) || (_i < 0))
                {
                    throw out_of_range("index out of range\n");
                }
                return (_i * (_i + 3)) >> 1;
            }

            if (_i > _j) 
            {
                if ((_i >= nSeqs) || (_j < 0))
                {
                    throw out_of_range("index out of range\n");
                }
                return ((_i * (_i + 1)) >> 1) + _j;
            }

            if ((_j >= nSeqs) || (_i < 0))
            {
                throw out_of_range("index out of range\n");
            }
            return  ((_j * (_j + 1)) >> 1) + _i;            
        }
        
        //inline
        inline double operator() (unsigned row, unsigned col) const
        {
            return elements[getIndex(row, col, numSeqs)];
        }
        
        inline double& operator() (unsigned row, unsigned col)
        {
            return elements[getIndex(row, col, numSeqs)];
        }
        
        inline double* getElements()
        {
            return elements;
        }
        
        inline double* getDistMatrix(int fSeq, int nSeqsInSub)
        {
            if(fSeq == 1 && numSeqs == nSeqsInSub)
            {
                return elements; // We want the whole array.
            }
            else
            {
                // Calculate the subMatrix and return it!!!
                try
                {
                    if(subElements)
                    {
                        delete [] subElements;
                    }
                    sizeSubElements = ((nSeqsInSub + 1) * (nSeqsInSub + 2)) >> 1;
                    numSeqsInSub = nSeqsInSub;                    
                    
                    subElements = new double[sizeSubElements];
                    setAllArrayToVal(subElements, sizeSubElements, 0.0);
                    
                    int currIndex = 0;
                    subElements[0] = 0.0;
                    int lSeq = fSeq + numSeqsInSub - 1; // NOTE this is wrong!!!!! Need to fix

                    for(int i = fSeq; i <= lSeq; i++)
                    {
                        for(int j = i + 1; j <= lSeq; j++)
                        {
                            currIndex = getIndex(i - fSeq + 1, j - fSeq + 1, numSeqsInSub);
                            subElements[currIndex] = elements[getIndex(i, j, numSeqs)];    
                        }
                    }

                    return subElements;
                }
                catch(bad_alloc& e)
                {
                    cout << "Out of Memory!\n";
                    throw e;
                }                
            }
        }
        
        void clearSubArray()
        {
            if(subElements)
            {
                delete [] subElements;
            }
            subElements = 0; 
            sizeSubElements = 0;
            numSeqsInSub = 0;       
        }
        
        void makeSimilarityMatrix()
        {
            double value;
            for (int i = 0; i < numSeqs; i++)
            {
                SetAt(i + 1, i + 1, 0.0);
                for (int j = 0; j < i; j++)
                {
                    value = 100.0 - (GetAt(i + 1, j + 1)) * 100.0;
                    SetAt(i + 1, j + 1, value);
                    //distMat->SetAt(j + 1, i + 1, value);
                }
            }
        }
        
    private:
        double* elements;
        int sizeElements;
        int numSeqs;
        int index;
        double* subElements; // To be used to return a sub matrix.
        int firstSeq, numSeqsInSub;
        int sizeSubElements;
};

}

#endif