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++ == 3000) 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 vector< vector<double> > LinearAlgebra::calculateEuclidianDistance(vector< vector<double> >& axes, int dimensions){
 238         try {
 239                 //make square matrix
 240                 vector< vector<double> > dists; dists.resize(axes.size());
 241                 for (int i = 0; i < dists.size(); i++) {  dists[i].resize(axes.size(), 0.0); }
 242
 243                 if (dimensions == 1) { //one dimension calc = abs(x-y)
 244
 245                         for (int i = 0; i < dists.size(); i++) {
 246
 247                                 if (m->control_pressed) { return dists; }
 248
 249                                 for (int j = 0; j < i; j++) {
 250                                         dists[i][j] = abs(axes[i][0] - axes[j][0]);
 251                                         dists[j][i] = dists[i][j];
 252                                 }
 253                         }
 254
 255                 }else if (dimensions == 2) { //two dimension calc = sqrt ((x1 - y1)^2 + (x2 - y2)^2)
 256
 257                         for (int i = 0; i < dists.size(); i++) {
 258
 259                                 if (m->control_pressed) { return dists; }
 260
 261                                 for (int j = 0; j < i; j++) {
 262                                         double firstDim = ((axes[i][0] - axes[j][0]) * (axes[i][0] - axes[j][0]));
 263                                         double secondDim = ((axes[i][1] - axes[j][1]) * (axes[i][1] - axes[j][1]));
 264
 265                                         dists[i][j] = sqrt((firstDim + secondDim));
 266                                         dists[j][i] = dists[i][j];
 267                                 }
 268                         }
 269
 270                 }else if (dimensions == 3) { //two dimension calc = sqrt ((x1 - y1)^2 + (x2 - y2)^2 + (x3 - y3)^2)
 271
 272                         for (int i = 0; i < dists.size(); i++) {
 273
 274                                 if (m->control_pressed) { return dists; }
 275
 276                                 for (int j = 0; j < i; j++) {
 277                                         double firstDim = ((axes[i][0] - axes[j][0]) * (axes[i][0] - axes[j][0]));
 278                                         double secondDim = ((axes[i][1] - axes[j][1]) * (axes[i][1] - axes[j][1]));
 279                                         double thirdDim = ((axes[i][2] - axes[j][2]) * (axes[i][2] - axes[j][2]));
 280
 281                                         dists[i][j] = sqrt((firstDim + secondDim + thirdDim));
 282                                         dists[j][i] = dists[i][j];
 283                                 }
 284                         }
 285
 286                 }else { m->mothurOut("[ERROR]: too many dimensions, aborting."); m->mothurOutEndLine(); m->control_pressed = true; }
 287
 288                 return dists;
 289         }
 290         catch(exception& e) {
 291                 m->errorOut(e, "LinearAlgebra", "calculateEuclidianDistance");
 292                 exit(1);
 293         }
 294 }
 295 /*********************************************************************************************************************************/
 296 double LinearAlgebra::calcPearson(vector< vector<double> >& euclidDists, vector< vector<double> >& userDists){
 297         try {
 298
 299                 //find average for - X
 300                 vector<float> averageEuclid; averageEuclid.resize(euclidDists.size(), 0.0);
 301                 for (int i = 0; i < euclidDists.size(); i++) {
 302                         for (int j = 0; j < euclidDists[i].size(); j++) {
 303                                 averageEuclid[i] += euclidDists[i][j];
 304                         }
 305                 }
 306                 for (int i = 0; i < averageEuclid.size(); i++) {  averageEuclid[i] = averageEuclid[i] / (float) euclidDists.size();   }
 307
 308                 //find average for - Y
 309                 vector<float> averageUser; averageUser.resize(userDists.size(), 0.0);
 310                 for (int i = 0; i < userDists.size(); i++) {
 311                         for (int j = 0; j < userDists[i].size(); j++) {
 312                                 averageUser[i] += userDists[i][j];
 313                         }
 314                 }
 315                 for (int i = 0; i < averageUser.size(); i++) {  averageUser[i] = averageUser[i] / (float) userDists.size();  }
 316
 317                 double numerator = 0.0;
 318                 double denomTerm1 = 0.0;
 319                 double denomTerm2 = 0.0;
 320
 321                 for (int i = 0; i < euclidDists.size(); i++) {
 322
 323                         for (int k = 0; k < i; k++) {
 324
 325                                 float Yi = userDists[i][k];
 326                                 float Xi = euclidDists[i][k];
 327
 328                                 numerator += ((Xi - averageEuclid[k]) * (Yi - averageUser[k]));
 329                                 denomTerm1 += ((Xi - averageEuclid[k]) * (Xi - averageEuclid[k]));
 330                                 denomTerm2 += ((Yi - averageUser[k]) * (Yi - averageUser[k]));
 331                         }
 332                 }
 333
 334                 double denom = (sqrt(denomTerm1) * sqrt(denomTerm2));
 335                 double r = numerator / denom;
 336
 337                 return r;
 338         }
 339         catch(exception& e) {
 340                 m->errorOut(e, "LinearAlgebra", "calculateEuclidianDistance");
 341                 exit(1);
 342         }
 343 }
 344 /*********************************************************************************************************************************/
 345
 346