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 // This function is taken from Numerical Recipes in C++ by Press et al., 2nd edition, pg. 479
58 int LinearAlgebra::tred2(vector<vector<double> >& a, vector<double>& d, vector<double>& e){
60 double scale, hh, h, g, f;
67 for(int i=n-1;i>0;i--){
71 for(int k=0;k<l+1;k++){
72 scale += fabs(a[i][k]);
78 for(int k=0;k<l+1;k++){
80 h += a[i][k] * a[i][k];
83 g = (f >= 0.0 ? -sqrt(h) : sqrt(h));
88 for(int j=0;j<l+1;j++){
89 a[j][i] = a[i][j] / h;
91 for(int k=0;k<j+1;k++){
92 g += a[j][k] * a[i][k];
94 for(int k=j+1;k<l+1;k++){
95 g += a[k][j] * a[i][k];
101 for(int j=0;j<l+1;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]);
120 for(int i=0;i<n;i++){
123 for(int j=0;j<l;j++){
125 for(int k=0;k<l;k++){
126 g += a[i][k] * a[k][j];
128 for(int k=0;k<l;k++){
129 a[k][j] -= g * a[k][i];
135 for(int j=0;j<l;j++){
136 a[j][i] = a[i][j] = 0.0;
142 catch(exception& e) {
143 m->errorOut(e, "LinearAlgebra", "tred2");
148 /*********************************************************************************************************************************/
150 double LinearAlgebra::pythag(double a, double b) { return(pow(a*a+b*b,0.5)); }
152 /*********************************************************************************************************************************/
154 // This function is taken from Numerical Recipes in C++ by Press et al., 2nd edition, pg. 479
156 int LinearAlgebra::qtli(vector<double>& d, vector<double>& e, vector<vector<double> >& z) {
159 double s, r, p, g, f, dd, c, b;
162 for(int i=1;i<=n;i++){
167 for(int l=0;l<n;l++){
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;
175 if(iter++ == 3000) cerr << "Too many iterations in tqli\n";
176 g = (d[l+1]-d[l]) / (2.0 * e[l]);
178 g = d[myM] - d[l] + e[l] / (g + SIGN(r,g));
181 for(i=myM-1;i>=l;i--){
184 e[i+1] = (r=pythag(f,g));
193 r = (d[i] - g) * s + 2.0 * c * b;
194 d[i+1] = g + ( p = s * r);
196 for(int k=0;k<n;k++){
198 z[k][i+1] = s * z[k][i] + c * f;
199 z[k][i] = c * z[k][i] - s * f;
202 if(r == 0.00 && i >= l) continue;
211 for(int i=0;i<n;i++){
213 for(int j=i;j<n;j++){
221 for(int j=0;j<n;j++){
231 catch(exception& e) {
232 m->errorOut(e, "LinearAlgebra", "qtli");
236 /*********************************************************************************************************************************/
237 //groups by dimension
238 vector< vector<double> > LinearAlgebra::calculateEuclidianDistance(vector< vector<double> >& axes, int dimensions){
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); }
244 if (dimensions == 1) { //one dimension calc = abs(x-y)
246 for (int i = 0; i < dists.size(); i++) {
248 if (m->control_pressed) { return dists; }
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];
256 }else if (dimensions > 1) { //two dimension calc = sqrt ((x1 - y1)^2 + (x2 - y2)^2)...
258 for (int i = 0; i < dists.size(); i++) {
260 if (m->control_pressed) { return dists; }
262 for (int j = 0; j < i; j++) {
264 for (int k = 0; k < dimensions; k++) {
265 sum += ((axes[i][k] - axes[j][k]) * (axes[i][k] - axes[j][k]));
268 dists[i][j] = sqrt(sum);
269 dists[j][i] = dists[i][j];
277 catch(exception& e) {
278 m->errorOut(e, "LinearAlgebra", "calculateEuclidianDistance");
282 /*********************************************************************************************************************************/
283 //returns groups by dimensions from dimensions by groups
284 vector< vector<double> > LinearAlgebra::calculateEuclidianDistance(vector< vector<double> >& axes){
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); }
290 if (axes.size() == 1) { //one dimension calc = abs(x-y)
292 for (int i = 0; i < dists.size(); i++) {
294 if (m->control_pressed) { return dists; }
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];
302 }else if (axes.size() > 1) { //two dimension calc = sqrt ((x1 - y1)^2 + (x2 - y2)^2)...
304 for (int i = 0; i < dists[0].size(); i++) {
306 if (m->control_pressed) { return dists; }
308 for (int j = 0; j < i; j++) {
310 for (int k = 0; k < axes.size(); k++) {
311 sum += ((axes[k][i] - axes[k][j]) * (axes[k][i] - axes[k][j]));
314 dists[i][j] = sqrt(sum);
315 dists[j][i] = dists[i][j];
323 catch(exception& e) {
324 m->errorOut(e, "LinearAlgebra", "calculateEuclidianDistance");
328 /*********************************************************************************************************************************/
329 //assumes both matrices are square and the same size
330 double LinearAlgebra::calcPearson(vector< vector<double> >& euclidDists, vector< vector<double> >& userDists){
333 //find average for - X
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];
342 averageEuclid = averageEuclid / (float) count;
344 //find average for - Y
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];
353 averageUser = averageUser / (float) count;
355 double numerator = 0.0;
356 double denomTerm1 = 0.0;
357 double denomTerm2 = 0.0;
359 for (int i = 0; i < euclidDists.size(); i++) {
361 for (int k = 0; k < i; k++) { //just lt dists
363 float Yi = userDists[i][k];
364 float Xi = euclidDists[i][k];
366 numerator += ((Xi - averageEuclid) * (Yi - averageUser));
367 denomTerm1 += ((Xi - averageEuclid) * (Xi - averageEuclid));
368 denomTerm2 += ((Yi - averageUser) * (Yi - averageUser));
372 double denom = (sqrt(denomTerm1) * sqrt(denomTerm2));
373 double r = numerator / denom;
375 //divide by zero error
376 if (isnan(r) || isinf(r)) { r = 0.0; }
381 catch(exception& e) {
382 m->errorOut(e, "LinearAlgebra", "calcPearson");
386 /*********************************************************************************************************************************/
387 //assumes both matrices are square and the same size
388 double LinearAlgebra::calcSpearman(vector< vector<double> >& euclidDists, vector< vector<double> >& userDists){
393 map<float, int> tableX;
394 map<float, int>::iterator itTable;
395 vector<spearmanRank> scores;
397 for (int i = 0; i < euclidDists.size(); i++) {
398 for (int j = 0; j < i; j++) {
399 spearmanRank member(toString(scores.size()), euclidDists[i][j]);
400 scores.push_back(member);
402 //count number of repeats
403 itTable = tableX.find(euclidDists[i][j]);
404 if (itTable == tableX.end()) {
405 tableX[euclidDists[i][j]] = 1;
407 tableX[euclidDists[i][j]]++;
413 sort(scores.begin(), scores.end(), compareSpearman);
417 for (itTable = tableX.begin(); itTable != tableX.end(); itTable++) {
418 double tx = (double) itTable->second;
419 Lx += ((pow(tx, 3.0) - tx) / 12.0);
423 map<string, float> rankEuclid;
424 vector<spearmanRank> ties;
426 for (int j = 0; j < scores.size(); j++) {
428 ties.push_back(scores[j]);
430 if (j != (scores.size()-1)) { // you are not the last so you can look ahead
431 if (scores[j].score != scores[j+1].score) { // you are done with ties, rank them and continue
433 for (int k = 0; k < ties.size(); k++) {
434 float thisrank = rankTotal / (float) ties.size();
435 rankEuclid[ties[k].name] = thisrank;
440 }else { // you are the last one
442 for (int k = 0; k < ties.size(); k++) {
443 float thisrank = rankTotal / (float) ties.size();
444 rankEuclid[ties[k].name] = thisrank;
451 map<float, int> tableY;
454 for (int i = 0; i < userDists.size(); i++) {
455 for (int j = 0; j < i; j++) {
456 spearmanRank member(toString(scores.size()), userDists[i][j]);
457 scores.push_back(member);
459 //count number of repeats
460 itTable = tableY.find(userDists[i][j]);
461 if (itTable == tableY.end()) {
462 tableY[userDists[i][j]] = 1;
464 tableY[userDists[i][j]]++;
470 sort(scores.begin(), scores.end(), compareSpearman);
474 for (itTable = tableY.begin(); itTable != tableY.end(); itTable++) {
475 double ty = (double) itTable->second;
476 Ly += ((pow(ty, 3.0) - ty) / 12.0);
480 map<string, float> rankUser;
483 for (int j = 0; j < scores.size(); j++) {
485 ties.push_back(scores[j]);
487 if (j != (scores.size()-1)) { // you are not the last so you can look ahead
488 if (scores[j].score != scores[j+1].score) { // you are done with ties, rank them and continue
490 for (int k = 0; k < ties.size(); k++) {
491 float thisrank = rankTotal / (float) ties.size();
492 rankUser[ties[k].name] = thisrank;
497 }else { // you are the last one
499 for (int k = 0; k < ties.size(); k++) {
500 float thisrank = rankTotal / (float) ties.size();
501 rankUser[ties[k].name] = thisrank;
509 for (int i = 0; i < userDists.size(); i++) {
510 for (int j = 0; j < i; j++) {
512 float xi = rankEuclid[toString(count)];
513 float yi = rankUser[toString(count)];
515 di += ((xi - yi) * (xi - yi));
521 double n = (double) count;
523 double SX2 = ((pow(n, 3.0) - n) / 12.0) - Lx;
524 double SY2 = ((pow(n, 3.0) - n) / 12.0) - Ly;
526 r = (SX2 + SY2 - di) / (2.0 * sqrt((SX2*SY2)));
528 //divide by zero error
529 if (isnan(r) || isinf(r)) { r = 0.0; }
534 catch(exception& e) {
535 m->errorOut(e, "LinearAlgebra", "calcSpearman");
540 /*********************************************************************************************************************************/
541 //assumes both matrices are square and the same size
542 double LinearAlgebra::calcKendall(vector< vector<double> >& euclidDists, vector< vector<double> >& userDists){
547 vector<spearmanRank> scores;
548 for (int i = 0; i < euclidDists.size(); i++) {
549 for (int j = 0; j < i; j++) {
550 spearmanRank member(toString(scores.size()), euclidDists[i][j]);
551 scores.push_back(member);
556 sort(scores.begin(), scores.end(), compareSpearman);
559 map<string, float> rankEuclid;
560 vector<spearmanRank> ties;
562 for (int j = 0; j < scores.size(); j++) {
564 ties.push_back(scores[j]);
566 if (j != (scores.size()-1)) { // you are not the last so you can look ahead
567 if (scores[j].score != scores[j+1].score) { // you are done with ties, rank them and continue
569 for (int k = 0; k < ties.size(); k++) {
570 float thisrank = rankTotal / (float) ties.size();
571 rankEuclid[ties[k].name] = thisrank;
576 }else { // you are the last one
578 for (int k = 0; k < ties.size(); k++) {
579 float thisrank = rankTotal / (float) ties.size();
580 rankEuclid[ties[k].name] = thisrank;
585 vector<spearmanRank> scoresUser;
586 for (int i = 0; i < userDists.size(); i++) {
587 for (int j = 0; j < i; j++) {
588 spearmanRank member(toString(scoresUser.size()), userDists[i][j]);
589 scoresUser.push_back(member);
594 sort(scoresUser.begin(), scoresUser.end(), compareSpearman);
597 map<string, float> rankUser;
600 for (int j = 0; j < scoresUser.size(); j++) {
602 ties.push_back(scoresUser[j]);
604 if (j != (scoresUser.size()-1)) { // you are not the last so you can look ahead
605 if (scoresUser[j].score != scoresUser[j+1].score) { // you are done with ties, rank them and continue
607 for (int k = 0; k < ties.size(); k++) {
608 float thisrank = rankTotal / (float) ties.size();
609 rankUser[ties[k].name] = thisrank;
614 }else { // you are the last one
616 for (int k = 0; k < ties.size(); k++) {
617 float thisrank = rankTotal / (float) ties.size();
618 rankUser[ties[k].name] = thisrank;
627 vector<spearmanRank> user;
628 for (int l = 0; l < scores.size(); l++) {
629 spearmanRank member(scores[l].name, rankUser[scores[l].name]);
630 user.push_back(member);
634 for (int l = 0; l < scores.size(); l++) {
636 int numWithHigherRank = 0;
637 int numWithLowerRank = 0;
638 float thisrank = user[l].score;
640 for (int u = l+1; u < scores.size(); u++) {
641 if (user[u].score > thisrank) { numWithHigherRank++; }
642 else if (user[u].score < thisrank) { numWithLowerRank++; }
646 numCoor += numWithHigherRank;
647 numDisCoor += numWithLowerRank;
650 r = (numCoor - numDisCoor) / (float) count;
652 //divide by zero error
653 if (isnan(r) || isinf(r)) { r = 0.0; }
658 catch(exception& e) {
659 m->errorOut(e, "LinearAlgebra", "calcKendall");
664 /*********************************************************************************************************************************/