5 * Created by westcott on 1/7/11.
6 * Copyright 2011 Schloss Lab. All rights reserved.
10 #include "linearalgebra.h"
12 /*********************************************************************************************************************************/
14 inline double SIGN(const double a, const double b)
16 return b>=0 ? (a>=0 ? a:-a) : (a>=0 ? -a:a);
18 /*********************************************************************************************************************************/
20 vector<vector<double> > LinearAlgebra::matrix_mult(vector<vector<double> > first, vector<vector<double> > second){
22 vector<vector<double> > product;
24 int first_rows = first.size();
25 int first_cols = first[0].size();
26 int second_cols = second[0].size();
28 product.resize(first_rows);
29 for(int i=0;i<first_rows;i++){
30 product[i].resize(second_cols);
33 for(int i=0;i<first_rows;i++){
34 for(int j=0;j<second_cols;j++){
36 if (m->control_pressed) { return product; }
39 for(int k=0;k<first_cols;k++){
40 product[i][j] += first[i][k] * second[k][j];
48 m->errorOut(e, "LinearAlgebra", "matrix_mult");
54 /*********************************************************************************************************************************/
56 void LinearAlgebra::recenter(double offset, vector<vector<double> > D, vector<vector<double> >& G){
60 vector<vector<double> > A(rank);
61 vector<vector<double> > C(rank);
62 for(int i=0;i<rank;i++){
67 double scale = -1.0000 / (double) rank;
69 for(int i=0;i<rank;i++){
71 C[i][i] = 1.0000 + scale;
72 for(int j=i+1;j<rank;j++){
73 A[i][j] = A[j][i] = -0.5 * D[i][j] * D[i][j] + offset;
74 C[i][j] = C[j][i] = scale;
82 m->errorOut(e, "LinearAlgebra", "recenter");
87 /*********************************************************************************************************************************/
89 // This function is taken from Numerical Recipes in C++ by Press et al., 2nd edition, pg. 479
91 int LinearAlgebra::tred2(vector<vector<double> >& a, vector<double>& d, vector<double>& e){
93 double scale, hh, h, g, f;
100 for(int i=n-1;i>0;i--){
104 for(int k=0;k<l+1;k++){
105 scale += fabs(a[i][k]);
111 for(int k=0;k<l+1;k++){
113 h += a[i][k] * a[i][k];
116 g = (f >= 0.0 ? -sqrt(h) : sqrt(h));
121 for(int j=0;j<l+1;j++){
122 a[j][i] = a[i][j] / h;
124 for(int k=0;k<j+1;k++){
125 g += a[j][k] * a[i][k];
127 for(int k=j+1;k<l+1;k++){
128 g += a[k][j] * a[i][k];
134 for(int j=0;j<l+1;j++){
136 e[j] = g = e[j] - hh * f;
137 for(int k=0;k<j+1;k++){
138 a[j][k] -= (f * e[k] + g * a[i][k]);
153 for(int i=0;i<n;i++){
156 for(int j=0;j<l;j++){
158 for(int k=0;k<l;k++){
159 g += a[i][k] * a[k][j];
161 for(int k=0;k<l;k++){
162 a[k][j] -= g * a[k][i];
168 for(int j=0;j<l;j++){
169 a[j][i] = a[i][j] = 0.0;
175 catch(exception& e) {
176 m->errorOut(e, "LinearAlgebra", "tred2");
181 /*********************************************************************************************************************************/
183 double LinearAlgebra::pythag(double a, double b) { return(pow(a*a+b*b,0.5)); }
185 /*********************************************************************************************************************************/
187 // This function is taken from Numerical Recipes in C++ by Press et al., 2nd edition, pg. 479
189 int LinearAlgebra::qtli(vector<double>& d, vector<double>& e, vector<vector<double> >& z) {
192 double s, r, p, g, f, dd, c, b;
195 for(int i=1;i<=n;i++){
200 for(int l=0;l<n;l++){
203 for(myM=l;myM<n-1;myM++){
204 dd = fabs(d[myM]) + fabs(d[myM+1]);
205 if(fabs(e[myM])+dd == dd) break;
208 if(iter++ == 3000) cerr << "Too many iterations in tqli\n";
209 g = (d[l+1]-d[l]) / (2.0 * e[l]);
211 g = d[myM] - d[l] + e[l] / (g + SIGN(r,g));
214 for(i=myM-1;i>=l;i--){
217 e[i+1] = (r=pythag(f,g));
226 r = (d[i] - g) * s + 2.0 * c * b;
227 d[i+1] = g + ( p = s * r);
229 for(int k=0;k<n;k++){
231 z[k][i+1] = s * z[k][i] + c * f;
232 z[k][i] = c * z[k][i] - s * f;
235 if(r == 0.00 && i >= l) continue;
244 for(int i=0;i<n;i++){
246 for(int j=i;j<n;j++){
254 for(int j=0;j<n;j++){
264 catch(exception& e) {
265 m->errorOut(e, "LinearAlgebra", "qtli");
269 /*********************************************************************************************************************************/
270 //groups by dimension
271 vector< vector<double> > LinearAlgebra::calculateEuclidianDistance(vector< vector<double> >& axes, int dimensions){
274 vector< vector<double> > dists; dists.resize(axes.size());
275 for (int i = 0; i < dists.size(); i++) { dists[i].resize(axes.size(), 0.0); }
277 if (dimensions == 1) { //one dimension calc = abs(x-y)
279 for (int i = 0; i < dists.size(); i++) {
281 if (m->control_pressed) { return dists; }
283 for (int j = 0; j < i; j++) {
284 dists[i][j] = abs(axes[i][0] - axes[j][0]);
285 dists[j][i] = dists[i][j];
289 }else if (dimensions > 1) { //two dimension calc = sqrt ((x1 - y1)^2 + (x2 - y2)^2)...
291 for (int i = 0; i < dists.size(); i++) {
293 if (m->control_pressed) { return dists; }
295 for (int j = 0; j < i; j++) {
297 for (int k = 0; k < dimensions; k++) {
298 sum += ((axes[i][k] - axes[j][k]) * (axes[i][k] - axes[j][k]));
301 dists[i][j] = sqrt(sum);
302 dists[j][i] = dists[i][j];
310 catch(exception& e) {
311 m->errorOut(e, "LinearAlgebra", "calculateEuclidianDistance");
315 /*********************************************************************************************************************************/
316 //returns groups by dimensions from dimensions by groups
317 vector< vector<double> > LinearAlgebra::calculateEuclidianDistance(vector< vector<double> >& axes){
320 vector< vector<double> > dists; dists.resize(axes[0].size());
321 for (int i = 0; i < dists.size(); i++) { dists[i].resize(axes[0].size(), 0.0); }
323 if (axes.size() == 1) { //one dimension calc = abs(x-y)
325 for (int i = 0; i < dists.size(); i++) {
327 if (m->control_pressed) { return dists; }
329 for (int j = 0; j < i; j++) {
330 dists[i][j] = abs(axes[0][i] - axes[0][j]);
331 dists[j][i] = dists[i][j];
335 }else if (axes.size() > 1) { //two dimension calc = sqrt ((x1 - y1)^2 + (x2 - y2)^2)...
337 for (int i = 0; i < dists[0].size(); i++) {
339 if (m->control_pressed) { return dists; }
341 for (int j = 0; j < i; j++) {
343 for (int k = 0; k < axes.size(); k++) {
344 sum += ((axes[k][i] - axes[k][j]) * (axes[k][i] - axes[k][j]));
347 dists[i][j] = sqrt(sum);
348 dists[j][i] = dists[i][j];
356 catch(exception& e) {
357 m->errorOut(e, "LinearAlgebra", "calculateEuclidianDistance");
361 /*********************************************************************************************************************************/
362 //assumes both matrices are square and the same size
363 double LinearAlgebra::calcPearson(vector< vector<double> >& euclidDists, vector< vector<double> >& userDists){
366 //find average for - X
368 float averageEuclid = 0.0;
369 for (int i = 0; i < euclidDists.size(); i++) {
370 for (int j = 0; j < i; j++) {
371 averageEuclid += euclidDists[i][j];
375 averageEuclid = averageEuclid / (float) count;
377 //find average for - Y
379 float averageUser = 0.0;
380 for (int i = 0; i < userDists.size(); i++) {
381 for (int j = 0; j < i; j++) {
382 averageUser += userDists[i][j];
386 averageUser = averageUser / (float) count;
388 double numerator = 0.0;
389 double denomTerm1 = 0.0;
390 double denomTerm2 = 0.0;
392 for (int i = 0; i < euclidDists.size(); i++) {
394 for (int k = 0; k < i; k++) { //just lt dists
396 float Yi = userDists[i][k];
397 float Xi = euclidDists[i][k];
399 numerator += ((Xi - averageEuclid) * (Yi - averageUser));
400 denomTerm1 += ((Xi - averageEuclid) * (Xi - averageEuclid));
401 denomTerm2 += ((Yi - averageUser) * (Yi - averageUser));
405 double denom = (sqrt(denomTerm1) * sqrt(denomTerm2));
406 double r = numerator / denom;
408 //divide by zero error
409 if (isnan(r) || isinf(r)) { r = 0.0; }
414 catch(exception& e) {
415 m->errorOut(e, "LinearAlgebra", "calcPearson");
419 /*********************************************************************************************************************************/
420 //assumes both matrices are square and the same size
421 double LinearAlgebra::calcSpearman(vector< vector<double> >& euclidDists, vector< vector<double> >& userDists){
426 map<float, int> tableX;
427 map<float, int>::iterator itTable;
428 vector<spearmanRank> scores;
430 for (int i = 0; i < euclidDists.size(); i++) {
431 for (int j = 0; j < i; j++) {
432 spearmanRank member(toString(scores.size()), euclidDists[i][j]);
433 scores.push_back(member);
435 //count number of repeats
436 itTable = tableX.find(euclidDists[i][j]);
437 if (itTable == tableX.end()) {
438 tableX[euclidDists[i][j]] = 1;
440 tableX[euclidDists[i][j]]++;
446 sort(scores.begin(), scores.end(), compareSpearman);
450 for (itTable = tableX.begin(); itTable != tableX.end(); itTable++) {
451 double tx = (double) itTable->second;
452 Lx += ((pow(tx, 3.0) - tx) / 12.0);
456 map<string, float> rankEuclid;
457 vector<spearmanRank> ties;
459 for (int j = 0; j < scores.size(); j++) {
461 ties.push_back(scores[j]);
463 if (j != (scores.size()-1)) { // you are not the last so you can look ahead
464 if (scores[j].score != scores[j+1].score) { // you are done with ties, rank them and continue
466 for (int k = 0; k < ties.size(); k++) {
467 float thisrank = rankTotal / (float) ties.size();
468 rankEuclid[ties[k].name] = thisrank;
473 }else { // you are the last one
475 for (int k = 0; k < ties.size(); k++) {
476 float thisrank = rankTotal / (float) ties.size();
477 rankEuclid[ties[k].name] = thisrank;
484 map<float, int> tableY;
487 for (int i = 0; i < userDists.size(); i++) {
488 for (int j = 0; j < i; j++) {
489 spearmanRank member(toString(scores.size()), userDists[i][j]);
490 scores.push_back(member);
492 //count number of repeats
493 itTable = tableY.find(userDists[i][j]);
494 if (itTable == tableY.end()) {
495 tableY[userDists[i][j]] = 1;
497 tableY[userDists[i][j]]++;
503 sort(scores.begin(), scores.end(), compareSpearman);
507 for (itTable = tableY.begin(); itTable != tableY.end(); itTable++) {
508 double ty = (double) itTable->second;
509 Ly += ((pow(ty, 3.0) - ty) / 12.0);
513 map<string, float> rankUser;
516 for (int j = 0; j < scores.size(); j++) {
518 ties.push_back(scores[j]);
520 if (j != (scores.size()-1)) { // you are not the last so you can look ahead
521 if (scores[j].score != scores[j+1].score) { // you are done with ties, rank them and continue
523 for (int k = 0; k < ties.size(); k++) {
524 float thisrank = rankTotal / (float) ties.size();
525 rankUser[ties[k].name] = thisrank;
530 }else { // you are the last one
532 for (int k = 0; k < ties.size(); k++) {
533 float thisrank = rankTotal / (float) ties.size();
534 rankUser[ties[k].name] = thisrank;
542 for (int i = 0; i < userDists.size(); i++) {
543 for (int j = 0; j < i; j++) {
545 float xi = rankEuclid[toString(count)];
546 float yi = rankUser[toString(count)];
548 di += ((xi - yi) * (xi - yi));
554 double n = (double) count;
556 double SX2 = ((pow(n, 3.0) - n) / 12.0) - Lx;
557 double SY2 = ((pow(n, 3.0) - n) / 12.0) - Ly;
559 r = (SX2 + SY2 - di) / (2.0 * sqrt((SX2*SY2)));
561 //divide by zero error
562 if (isnan(r) || isinf(r)) { r = 0.0; }
567 catch(exception& e) {
568 m->errorOut(e, "LinearAlgebra", "calcSpearman");
573 /*********************************************************************************************************************************/
574 //assumes both matrices are square and the same size
575 double LinearAlgebra::calcKendall(vector< vector<double> >& euclidDists, vector< vector<double> >& userDists){
580 vector<spearmanRank> scores;
581 for (int i = 0; i < euclidDists.size(); i++) {
582 for (int j = 0; j < i; j++) {
583 spearmanRank member(toString(scores.size()), euclidDists[i][j]);
584 scores.push_back(member);
589 sort(scores.begin(), scores.end(), compareSpearman);
592 map<string, float> rankEuclid;
593 vector<spearmanRank> ties;
595 for (int j = 0; j < scores.size(); j++) {
597 ties.push_back(scores[j]);
599 if (j != (scores.size()-1)) { // you are not the last so you can look ahead
600 if (scores[j].score != scores[j+1].score) { // you are done with ties, rank them and continue
602 for (int k = 0; k < ties.size(); k++) {
603 float thisrank = rankTotal / (float) ties.size();
604 rankEuclid[ties[k].name] = thisrank;
609 }else { // you are the last one
611 for (int k = 0; k < ties.size(); k++) {
612 float thisrank = rankTotal / (float) ties.size();
613 rankEuclid[ties[k].name] = thisrank;
618 vector<spearmanRank> scoresUser;
619 for (int i = 0; i < userDists.size(); i++) {
620 for (int j = 0; j < i; j++) {
621 spearmanRank member(toString(scoresUser.size()), userDists[i][j]);
622 scoresUser.push_back(member);
627 sort(scoresUser.begin(), scoresUser.end(), compareSpearman);
630 map<string, float> rankUser;
633 for (int j = 0; j < scoresUser.size(); j++) {
635 ties.push_back(scoresUser[j]);
637 if (j != (scoresUser.size()-1)) { // you are not the last so you can look ahead
638 if (scoresUser[j].score != scoresUser[j+1].score) { // you are done with ties, rank them and continue
640 for (int k = 0; k < ties.size(); k++) {
641 float thisrank = rankTotal / (float) ties.size();
642 rankUser[ties[k].name] = thisrank;
647 }else { // you are the last one
649 for (int k = 0; k < ties.size(); k++) {
650 float thisrank = rankTotal / (float) ties.size();
651 rankUser[ties[k].name] = thisrank;
660 vector<spearmanRank> user;
661 for (int l = 0; l < scores.size(); l++) {
662 spearmanRank member(scores[l].name, rankUser[scores[l].name]);
663 user.push_back(member);
667 for (int l = 0; l < scores.size(); l++) {
669 int numWithHigherRank = 0;
670 int numWithLowerRank = 0;
671 float thisrank = user[l].score;
673 for (int u = l+1; u < scores.size(); u++) {
674 if (user[u].score > thisrank) { numWithHigherRank++; }
675 else if (user[u].score < thisrank) { numWithLowerRank++; }
679 numCoor += numWithHigherRank;
680 numDisCoor += numWithLowerRank;
683 r = (numCoor - numDisCoor) / (float) count;
685 //divide by zero error
686 if (isnan(r) || isinf(r)) { r = 0.0; }
691 catch(exception& e) {
692 m->errorOut(e, "LinearAlgebra", "calcKendall");
697 /*********************************************************************************************************************************/
699 vector<vector<double> > LinearAlgebra::getObservedEuclideanDistance(vector<vector<double> >& relAbundData){
701 int numSamples = relAbundData.size();
702 int numOTUs = relAbundData[0].size();
704 vector<vector<double> > dMatrix(numSamples);
705 for(int i=0;i<numSamples;i++){
706 dMatrix[i].resize(numSamples);
709 for(int i=0;i<numSamples;i++){
710 for(int j=0;j<numSamples;j++){
713 for(int k=0;k<numOTUs;k++){
714 d += pow((relAbundData[i][k] - relAbundData[j][k]), 2.0000);
716 dMatrix[i][j] = pow(d, 0.50000);
717 dMatrix[j][i] = dMatrix[i][j];
725 /*********************************************************************************************************************************/