binaries/src/clustalw/src/general/SymMatrix.h

   1 /**
   2  * Author: Mark Larkin
   3  *
   4  * Copyright (c) 2007 Des Higgins, Julie Thompson and Toby Gibson.
   5  */
   6 /**
   7  * The SymMatrix class is used to store the distance matrix. It stores it in a double
   8  * array.
   9  * This class throws an out_of_range exception if we try to access an element outside
  10  * of the array bounds. It will throw a bad_alloc if we cannot allocate enough memory for
  11  * the distance matrix.
  12  */
  13 #ifndef SYMMATRIX_H
  14 #define SYMMATRIX_H
  15
  16 #include<iostream>
  17 #include <vector>
  18 #include <iomanip>
  19 #include <cmath>
  20 #include <stdexcept>
  21 #include <new>
  22 #include <cstdlib>
  23
  24 namespace clustalw
  25 {
  26 using namespace std;
  27
  28 class SymMatrix
  29 {
  30     public:
  31         SymMatrix() : elements(0), numSeqs(0), subElements(0), firstSeq(0), numSeqsInSub(0){;}
  32         SymMatrix(int size)
  33          : elements(0),
  34            numSeqs(0),
  35            subElements(0),
  36            firstSeq(0),
  37            numSeqsInSub(0)
  38         {
  39             // Size is the numSeqs + 1
  40             numSeqs = size - 1;
  41             sizeElements = ((numSeqs + 1) * (numSeqs + 2)) >> 1;
  42             try
  43             {
  44                 elements = new double[sizeElements];
  45                 setAllArrayToVal(elements, sizeElements, 0.0);
  46             }
  47             catch(bad_alloc& e)
  48             {
  49                 cout << "Could not allocate a distance matrix for " << numSeqs
  50                      << " seqs.\n";
  51                 throw e;
  52             }
  53
  54         }
  55         ~SymMatrix()
  56         {
  57             if(elements)
  58             {
  59                 delete [] elements;
  60             }
  61             if(subElements)
  62             {
  63                 delete [] subElements;
  64             }
  65         }
  66
  67         inline void SetAt(int nRow, int nCol, const double& value)
  68         {
  69             index = getIndex(nRow, nCol, numSeqs);
  70             elements[index] = value;
  71         }
  72
  73         inline double GetAt(int nRow, int nCol)
  74         {
  75             index = getIndex(nRow, nCol, numSeqs);
  76             return elements[index];
  77         }
  78
  79         inline void ResizeRect(int size, double val = 0.0)
  80         {
  81             numSeqs = size - 1;
  82             sizeElements = ((numSeqs + 1) * (numSeqs + 2)) >> 1;
  83             if(elements)
  84             {
  85                 delete [] elements;
  86             }
  87             try
  88             {
  89                 elements = new double[sizeElements];
  90                 setAllArrayToVal(elements, sizeElements, val);
  91             }
  92             catch(bad_alloc& e)
  93             {
  94                 cout << "Could not allocate a distance matrix for " << numSeqs
  95                      << " seqs. Need to terminate program.\n";
  96                 throw e;
  97             }
  98         }
  99
 100         inline void setAllArrayToVal(double* array, int size, double val)
 101         {
 102             for(int i = 0; i < size; i++)
 103             {
 104                 array[i] = val;
 105             }
 106         }
 107
 108         inline int getSize()
 109         {
 110             return numSeqs + 1;
 111         }
 112
 113         inline void clearArray()
 114         {
 115             sizeElements = 0;
 116             numSeqs = 0;
 117             if(elements)
 118             {
 119                 delete [] elements;
 120                 elements = 0;
 121             }
 122
 123             numSeqsInSub = 0;
 124             sizeSubElements = 0;
 125
 126             if(subElements)
 127             {
 128                 delete [] subElements;
 129                 subElements = 0;
 130             }
 131         }
 132
 133         void printArray()
 134         {
 135             printArray(elements, sizeElements);
 136         }
 137
 138         void printSubArray()
 139         {
 140             printArray(subElements, sizeSubElements);
 141         }
 142
 143         void printArray(double* array, int size)
 144         {
 145             int numThisTime = 1;
 146             int numSoFar = 0;
 147             for(int i = 0; i < size; i++)
 148             {
 149
 150                 numSoFar++;
 151
 152                 if(numSoFar == numThisTime)
 153                 {
 154                     cout << " " << setw(4) << array[i] << "\n";
 155                     numSoFar = 0;
 156                     numThisTime++;
 157                 }
 158                 else
 159                 {
 160                     cout << " " << setw(4) << array[i];
 161                 }
 162             }
 163             cout << "\n";
 164         }
 165
 166         inline int getIndex(const int &i, const int &j, const int &nSeqs) const
 167         {
 168             if(i == 0 || j == 0)
 169             {
 170                 return 0;
 171             }
 172
 173             int _i = i - 1, _j = j - 1;
 174             if (_i == _j)
 175             {
 176                 if ((_i >= nSeqs) || (_i < 0))
 177                 {
 178                     throw out_of_range("index out of range\n");
 179                 }
 180                 return (_i * (_i + 3)) >> 1;
 181             }
 182
 183             if (_i > _j)
 184             {
 185                 if ((_i >= nSeqs) || (_j < 0))
 186                 {
 187                     throw out_of_range("index out of range\n");
 188                 }
 189                 return ((_i * (_i + 1)) >> 1) + _j;
 190             }
 191
 192             if ((_j >= nSeqs) || (_i < 0))
 193             {
 194                 throw out_of_range("index out of range\n");
 195             }
 196             return  ((_j * (_j + 1)) >> 1) + _i;
 197         }
 198
 199         //inline
 200         inline double operator() (unsigned row, unsigned col) const
 201         {
 202             return elements[getIndex(row, col, numSeqs)];
 203         }
 204
 205         inline double& operator() (unsigned row, unsigned col)
 206         {
 207             return elements[getIndex(row, col, numSeqs)];
 208         }
 209
 210         inline double* getElements()
 211         {
 212             return elements;
 213         }
 214
 215         inline double* getDistMatrix(int fSeq, int nSeqsInSub)
 216         {
 217             if(fSeq == 1 && numSeqs == nSeqsInSub)
 218             {
 219                 return elements; // We want the whole array.
 220             }
 221             else
 222             {
 223                 // Calculate the subMatrix and return it!!!
 224                 try
 225                 {
 226                     if(subElements)
 227                     {
 228                         delete [] subElements;
 229                     }
 230                     sizeSubElements = ((nSeqsInSub + 1) * (nSeqsInSub + 2)) >> 1;
 231                     numSeqsInSub = nSeqsInSub;
 232
 233                     subElements = new double[sizeSubElements];
 234                     setAllArrayToVal(subElements, sizeSubElements, 0.0);
 235
 236                     int currIndex = 0;
 237                     subElements[0] = 0.0;
 238                     int lSeq = fSeq + numSeqsInSub - 1; // NOTE this is wrong!!!!! Need to fix
 239
 240                     for(int i = fSeq; i <= lSeq; i++)
 241                     {
 242                         for(int j = i + 1; j <= lSeq; j++)
 243                         {
 244                             currIndex = getIndex(i - fSeq + 1, j - fSeq + 1, numSeqsInSub);
 245                             subElements[currIndex] = elements[getIndex(i, j, numSeqs)];
 246                         }
 247                     }
 248
 249                     return subElements;
 250                 }
 251                 catch(bad_alloc& e)
 252                 {
 253                     cout << "Out of Memory!\n";
 254                     throw e;
 255                 }
 256             }
 257         }
 258
 259         void clearSubArray()
 260         {
 261             if(subElements)
 262             {
 263                 delete [] subElements;
 264             }
 265             subElements = 0;
 266             sizeSubElements = 0;
 267             numSeqsInSub = 0;
 268         }
 269
 270         void makeSimilarityMatrix()
 271         {
 272             double value;
 273             for (int i = 0; i < numSeqs; i++)
 274             {
 275                 SetAt(i + 1, i + 1, 0.0);
 276                 for (int j = 0; j < i; j++)
 277                 {
 278                     value = 100.0 - (GetAt(i + 1, j + 1)) * 100.0;
 279                     SetAt(i + 1, j + 1, value);
 280                     //distMat->SetAt(j + 1, i + 1, value);
 281                 }
 282             }
 283         }
 284
 285     private:
 286         double* elements;
 287         int sizeElements;
 288         int numSeqs;
 289         int index;
 290         double* subElements; // To be used to return a sub matrix.
 291         int firstSeq, numSeqsInSub;
 292         int sizeSubElements;
 293 };
 294
 295 }
 296
 297 #endif