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 //groups by dimension
 238 vector< vector<double> > LinearAlgebra::calculateEuclidianDistance(vector< vector<double> >& axes, int dimensions){
 239         try {
 240                 //make square matrix
 241                 vector< vector<double> > dists; dists.resize(axes.size());
 242                 for (int i = 0; i < dists.size(); i++) {  dists[i].resize(axes.size(), 0.0); }
 243
 244                 if (dimensions == 1) { //one dimension calc = abs(x-y)
 245
 246                         for (int i = 0; i < dists.size(); i++) {
 247
 248                                 if (m->control_pressed) { return dists; }
 249
 250                                 for (int j = 0; j < i; j++) {
 251                                         dists[i][j] = abs(axes[i][0] - axes[j][0]);
 252                                         dists[j][i] = dists[i][j];
 253                                 }
 254                         }
 255
 256                 }else if (dimensions > 1) { //two dimension calc = sqrt ((x1 - y1)^2 + (x2 - y2)^2)...
 257
 258                         for (int i = 0; i < dists.size(); i++) {
 259
 260                                 if (m->control_pressed) { return dists; }
 261
 262                                 for (int j = 0; j < i; j++) {
 263                                         double sum = 0.0;
 264                                         for (int k = 0; k < dimensions; k++) {
 265                                                 sum += ((axes[i][k] - axes[j][k]) * (axes[i][k] - axes[j][k]));
 266                                         }
 267
 268                                         dists[i][j] = sqrt(sum);
 269                                         dists[j][i] = dists[i][j];
 270                                 }
 271                         }
 272
 273                 }
 274
 275                 return dists;
 276         }
 277         catch(exception& e) {
 278                 m->errorOut(e, "LinearAlgebra", "calculateEuclidianDistance");
 279                 exit(1);
 280         }
 281 }
 282 /*********************************************************************************************************************************/
 283 //returns groups by dimensions from dimensions by groups
 284 vector< vector<double> > LinearAlgebra::calculateEuclidianDistance(vector< vector<double> >& axes){
 285         try {
 286                 //make square matrix
 287                 vector< vector<double> > dists; dists.resize(axes[0].size());
 288                 for (int i = 0; i < dists.size(); i++) {  dists[i].resize(axes[0].size(), 0.0); }
 289
 290                 if (axes.size() == 1) { //one dimension calc = abs(x-y)
 291
 292                         for (int i = 0; i < dists.size(); i++) {
 293
 294                                 if (m->control_pressed) { return dists; }
 295
 296                                 for (int j = 0; j < i; j++) {
 297                                         dists[i][j] = abs(axes[0][i] - axes[0][j]);
 298                                         dists[j][i] = dists[i][j];
 299                                 }
 300                         }
 301
 302                 }else if (axes.size() > 1) { //two dimension calc = sqrt ((x1 - y1)^2 + (x2 - y2)^2)...
 303
 304                         for (int i = 0; i < dists[0].size(); i++) {
 305
 306                                 if (m->control_pressed) { return dists; }
 307
 308                                 for (int j = 0; j < i; j++) {
 309                                         double sum = 0.0;
 310                                         for (int k = 0; k < axes.size(); k++) {
 311                                                 sum += ((axes[k][i] - axes[k][j]) * (axes[k][i] - axes[k][j]));
 312                                         }
 313
 314                                         dists[i][j] = sqrt(sum);
 315                                         dists[j][i] = dists[i][j];
 316                                 }
 317                         }
 318
 319                 }
 320
 321                 return dists;
 322         }
 323         catch(exception& e) {
 324                 m->errorOut(e, "LinearAlgebra", "calculateEuclidianDistance");
 325                 exit(1);
 326         }
 327 }
 328 /*********************************************************************************************************************************/
 329 //assumes both matrices are square and the same size
 330 double LinearAlgebra::calcPearson(vector< vector<double> >& euclidDists, vector< vector<double> >& userDists){
 331         try {
 332
 333                 //find average for - X
 334                 int count = 0;
 335                 float averageEuclid = 0.0;
 336                 for (int i = 0; i < euclidDists.size(); i++) {
 337                         for (int j = 0; j < i; j++) {
 338                                 averageEuclid += euclidDists[i][j];
 339                                 count++;
 340                         }
 341                 }
 342                 averageEuclid = averageEuclid / (float) count;
 343
 344                 //find average for - Y
 345                 count = 0;
 346                 float averageUser = 0.0;
 347                 for (int i = 0; i < userDists.size(); i++) {
 348                         for (int j = 0; j < i; j++) {
 349                                 averageUser += userDists[i][j];
 350                                 count++;
 351                         }
 352                 }
 353                 averageUser = averageUser / (float) count;
 354
 355                 double numerator = 0.0;
 356                 double denomTerm1 = 0.0;
 357                 double denomTerm2 = 0.0;
 358
 359                 for (int i = 0; i < euclidDists.size(); i++) {
 360
 361                         for (int k = 0; k < i; k++) { //just lt dists
 362
 363                                 float Yi = userDists[i][k];
 364                                 float Xi = euclidDists[i][k];
 365
 366                                 numerator += ((Xi - averageEuclid) * (Yi - averageUser));
 367                                 denomTerm1 += ((Xi - averageEuclid) * (Xi - averageEuclid));
 368                                 denomTerm2 += ((Yi - averageUser) * (Yi - averageUser));
 369                         }
 370                 }
 371
 372                 double denom = (sqrt(denomTerm1) * sqrt(denomTerm2));
 373                 double r = numerator / denom;
 374
 375                 //divide by zero error
 376                 if (isnan(r) || isinf(r)) { r = 0.0; }
 377
 378                 return r;
 379
 380         }
 381         catch(exception& e) {
 382                 m->errorOut(e, "LinearAlgebra", "calculateEuclidianDistance");
 383                 exit(1);
 384         }
 385 }
 386 /*********************************************************************************************************************************/
 387
 388