linearalgebra.cpp

   1 /*
   2  *  linearalgebra.cpp
   3  *  mothur
   4  *
   5  *  Created by westcott on 1/7/11.
   6  *  Copyright 2011 Schloss Lab. All rights reserved.
   7  *
   8  */
   9
  10 #include "linearalgebra.h"
  11
  12 /*********************************************************************************************************************************/
  13
  14 inline double SIGN(const double a, const double b)
  15 {
  16     return b>=0 ? (a>=0 ? a:-a) : (a>=0 ? -a:a);
  17 }
  18 /*********************************************************************************************************************************/
  19
  20 vector<vector<double> > LinearAlgebra::matrix_mult(vector<vector<double> > first, vector<vector<double> > second){
  21         try {
  22                 vector<vector<double> > product;
  23
  24                 int first_rows = first.size();
  25                 int first_cols = first[0].size();
  26                 int second_cols = second[0].size();
  27
  28                 product.resize(first_rows);
  29                 for(int i=0;i<first_rows;i++){
  30                         product[i].resize(second_cols);
  31                 }
  32
  33                 for(int i=0;i<first_rows;i++){
  34                         for(int j=0;j<second_cols;j++){
  35
  36                                 if (m->control_pressed) { return product; }
  37
  38                                 product[i][j] = 0.0;
  39                                 for(int k=0;k<first_cols;k++){
  40                                         product[i][j] += first[i][k] * second[k][j];
  41                                 }
  42                         }
  43                 }
  44
  45                 return product;
  46         }
  47         catch(exception& e) {
  48                 m->errorOut(e, "LinearAlgebra", "matrix_mult");
  49                 exit(1);
  50         }
  51
  52 }
  53
  54 /*********************************************************************************************************************************/
  55
  56 //  This function is taken from Numerical Recipes in C++ by Press et al., 2nd edition, pg. 479
  57
  58 int LinearAlgebra::tred2(vector<vector<double> >& a, vector<double>& d, vector<double>& e){
  59         try {
  60                 double scale, hh, h, g, f;
  61
  62                 int n = a.size();
  63
  64                 d.resize(n);
  65                 e.resize(n);
  66
  67                 for(int i=n-1;i>0;i--){
  68                         int l=i-1;
  69                         h = scale = 0.0000;
  70                         if(l>0){
  71                                 for(int k=0;k<l+1;k++){
  72                                         scale += fabs(a[i][k]);
  73                                 }
  74                                 if(scale == 0.0){
  75                                         e[i] = a[i][l];
  76                                 }
  77                                 else{
  78                                         for(int k=0;k<l+1;k++){
  79                                                 a[i][k] /= scale;
  80                                                 h += a[i][k] * a[i][k];
  81                                         }
  82                                         f = a[i][l];
  83                                         g = (f >= 0.0 ? -sqrt(h) : sqrt(h));
  84                                         e[i] = scale * g;
  85                                         h -= f * g;
  86                                         a[i][l] = f - g;
  87                                         f = 0.0;
  88                                         for(int j=0;j<l+1;j++){
  89                                                 a[j][i] = a[i][j] / h;
  90                                                 g = 0.0;
  91                                                 for(int k=0;k<j+1;k++){
  92                                                         g += a[j][k] * a[i][k];
  93                                                 }
  94                                                 for(int k=j+1;k<l+1;k++){
  95                                                         g += a[k][j] * a[i][k];
  96                                                 }
  97                                                 e[j] = g / h;
  98                                                 f += e[j] * a[i][j];
  99                                         }
 100                                         hh = f / (h + h);
 101                                         for(int j=0;j<l+1;j++){
 102                                                 f = a[i][j];
 103                                                 e[j] = g = e[j] - hh * f;
 104                                                 for(int k=0;k<j+1;k++){
 105                                                         a[j][k] -= (f * e[k] + g * a[i][k]);
 106                                                 }
 107                                         }
 108                                 }
 109                         }
 110                         else{
 111                                 e[i] = a[i][l];
 112                         }
 113
 114                         d[i] = h;
 115                 }
 116
 117                 d[0] = 0.0000;
 118                 e[0] = 0.0000;
 119
 120                 for(int i=0;i<n;i++){
 121                         int l = i;
 122                         if(d[i] != 0.0){
 123                                 for(int j=0;j<l;j++){
 124                                         g = 0.0000;
 125                                         for(int k=0;k<l;k++){
 126                                                 g += a[i][k] * a[k][j];
 127                                         }
 128                                         for(int k=0;k<l;k++){
 129                                                 a[k][j] -= g * a[k][i];
 130                                         }
 131                                 }
 132                         }
 133                         d[i] = a[i][i];
 134                         a[i][i] = 1.0000;
 135                         for(int j=0;j<l;j++){
 136                                 a[j][i] = a[i][j] = 0.0;
 137                         }
 138                 }
 139
 140                 return 0;
 141         }
 142         catch(exception& e) {
 143                 m->errorOut(e, "LinearAlgebra", "tred2");
 144                 exit(1);
 145         }
 146
 147 }
 148 /*********************************************************************************************************************************/
 149
 150 double LinearAlgebra::pythag(double a, double b)        {       return(pow(a*a+b*b,0.5));       }
 151
 152 /*********************************************************************************************************************************/
 153
 154 //  This function is taken from Numerical Recipes in C++ by Press et al., 2nd edition, pg. 479
 155
 156 int LinearAlgebra::qtli(vector<double>& d, vector<double>& e, vector<vector<double> >& z) {
 157         try {
 158                 int myM, i, iter;
 159                 double s, r, p, g, f, dd, c, b;
 160
 161                 int n = d.size();
 162                 for(int i=1;i<=n;i++){
 163                         e[i-1] = e[i];
 164                 }
 165                 e[n-1] = 0.0000;
 166
 167                 for(int l=0;l<n;l++){
 168                         iter = 0;
 169                         do {
 170                                 for(myM=l;myM<n-1;myM++){
 171                                         dd = fabs(d[myM]) + fabs(d[myM+1]);
 172                                         if(fabs(e[myM])+dd == dd) break;
 173                                 }
 174                                 if(myM != l){
 175                                         if(iter++ == 30) cerr << "Too many iterations in tqli\n";
 176                                         g = (d[l+1]-d[l]) / (2.0 * e[l]);
 177                                         r = pythag(g, 1.0);
 178                                         g = d[myM] - d[l] + e[l] / (g + SIGN(r,g));
 179                                         s = c = 1.0;
 180                                         p = 0.0000;
 181                                         for(i=myM-1;i>=l;i--){
 182                                                 f = s * e[i];
 183                                                 b = c * e[i];
 184                                                 e[i+1] = (r=pythag(f,g));
 185                                                 if(r==0.0){
 186                                                         d[i+1] -= p;
 187                                                         e[myM] = 0.0000;
 188                                                         break;
 189                                                 }
 190                                                 s = f / r;
 191                                                 c = g / r;
 192                                                 g = d[i+1] - p;
 193                                                 r = (d[i] - g) * s + 2.0 * c * b;
 194                                                 d[i+1] = g + ( p = s * r);
 195                                                 g = c * r - b;
 196                                                 for(int k=0;k<n;k++){
 197                                                         f = z[k][i+1];
 198                                                         z[k][i+1] = s * z[k][i] + c * f;
 199                                                         z[k][i] = c * z[k][i] - s * f;
 200                                                 }
 201                                         }
 202                                         if(r == 0.00 && i >= l) continue;
 203                                         d[l] -= p;
 204                                         e[l] = g;
 205                                         e[myM] = 0.0;
 206                                 }
 207                         } while (myM != l);
 208                 }
 209
 210                 int k;
 211                 for(int i=0;i<n;i++){
 212                         p=d[k=i];
 213                         for(int j=i;j<n;j++){
 214                                 if(d[j] >= p){
 215                                         p=d[k=j];
 216                                 }
 217                         }
 218                         if(k!=i){
 219                                 d[k]=d[i];
 220                                 d[i]=p;
 221                                 for(int j=0;j<n;j++){
 222                                         p=z[j][i];
 223                                         z[j][i] = z[j][k];
 224                                         z[j][k] = p;
 225                                 }
 226                         }
 227                 }
 228
 229                 return 0;
 230         }
 231         catch(exception& e) {
 232                 m->errorOut(e, "LinearAlgebra", "qtli");
 233                 exit(1);
 234         }
 235 }
 236 /*********************************************************************************************************************************/
 237
 238