X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=linearalgebra.cpp;h=27d35ac2cd7c40d4502aa469ffa62f2d01ff31ec;hb=45b880733289694b2a23c9dc8a774aa626458f81;hp=e015b5fa1472d9e587a8f24fe008f6af4f049c84;hpb=605ab6fa594317a38f0df7bb6797740c735b2348;p=mothur.git

diff --git a/linearalgebra.cpp b/linearalgebra.cpp
index e015b5f..27d35ac 100644
--- a/linearalgebra.cpp
+++ b/linearalgebra.cpp
@@ -172,7 +172,7 @@ int LinearAlgebra::qtli(vector<double>& d, vector<double>& e, vector<vector<doub
 					if(fabs(e[myM])+dd == dd) break;
 				}
 				if(myM != l){
-					if(iter++ == 30) cerr << "Too many iterations in tqli\n";
+					if(iter++ == 3000) cerr << "Too many iterations in tqli\n";
 					g = (d[l+1]-d[l]) / (2.0 * e[l]);
 					r = pythag(g, 1.0);
 					g = d[myM] - d[l] + e[l] / (g + SIGN(r,g));
@@ -234,5 +234,152 @@ int LinearAlgebra::qtli(vector<double>& d, vector<double>& e, vector<vector<doub
 	}
 }
 /*********************************************************************************************************************************/
+//groups by dimension
+vector< vector<double> > LinearAlgebra::calculateEuclidianDistance(vector< vector<double> >& axes, int dimensions){
+	try {
+		//make square matrix
+		vector< vector<double> > dists; dists.resize(axes.size());
+		for (int i = 0; i < dists.size(); i++) {  dists[i].resize(axes.size(), 0.0); }
+		
+		if (dimensions == 1) { //one dimension calc = abs(x-y)
+			
+			for (int i = 0; i < dists.size(); i++) {
+				
+				if (m->control_pressed) { return dists; }
+				
+				for (int j = 0; j < i; j++) {
+					dists[i][j] = abs(axes[i][0] - axes[j][0]);
+					dists[j][i] = dists[i][j];
+				}
+			}
+			
+		}else if (dimensions > 1) { //two dimension calc = sqrt ((x1 - y1)^2 + (x2 - y2)^2)...
+			
+			for (int i = 0; i < dists.size(); i++) {
+				
+				if (m->control_pressed) { return dists; }
+				
+				for (int j = 0; j < i; j++) {
+					double sum = 0.0;
+					for (int k = 0; k < dimensions; k++) {
+						sum += ((axes[i][k] - axes[j][k]) * (axes[i][k] - axes[j][k]));
+					}
+					
+					dists[i][j] = sqrt(sum);
+					dists[j][i] = dists[i][j];
+				}
+			}
+			
+		}
+		
+		return dists;
+	}
+	catch(exception& e) {
+		m->errorOut(e, "LinearAlgebra", "calculateEuclidianDistance");
+		exit(1);
+	}
+}
+/*********************************************************************************************************************************/
+//returns groups by dimensions from dimensions by groups
+vector< vector<double> > LinearAlgebra::calculateEuclidianDistance(vector< vector<double> >& axes){
+	try {
+		//make square matrix
+		vector< vector<double> > dists; dists.resize(axes[0].size());
+		for (int i = 0; i < dists.size(); i++) {  dists[i].resize(axes[0].size(), 0.0); }
+		
+		if (axes.size() == 1) { //one dimension calc = abs(x-y)
+			
+			for (int i = 0; i < dists.size(); i++) {
+				
+				if (m->control_pressed) { return dists; }
+				
+				for (int j = 0; j < i; j++) {
+					dists[i][j] = abs(axes[0][i] - axes[0][j]);
+					dists[j][i] = dists[i][j];
+				}
+			}
+			
+		}else if (axes.size() > 1) { //two dimension calc = sqrt ((x1 - y1)^2 + (x2 - y2)^2)...
+			
+			for (int i = 0; i < dists[0].size(); i++) {
+				
+				if (m->control_pressed) { return dists; }
+				
+				for (int j = 0; j < i; j++) {
+					double sum = 0.0;
+					for (int k = 0; k < axes.size(); k++) {
+						sum += ((axes[k][i] - axes[k][j]) * (axes[k][i] - axes[k][j]));
+					}
+					
+					dists[i][j] = sqrt(sum);
+					dists[j][i] = dists[i][j];
+				}
+			}
+			
+		}
+		
+		return dists;
+	}
+	catch(exception& e) {
+		m->errorOut(e, "LinearAlgebra", "calculateEuclidianDistance");
+		exit(1);
+	}
+}
+/*********************************************************************************************************************************/
+//assumes both matrices are square and the same size
+double LinearAlgebra::calcPearson(vector< vector<double> >& euclidDists, vector< vector<double> >& userDists){
+	try {
+		
+		//find average for - X
+		int count = 0;
+		float averageEuclid = 0.0; 
+		for (int i = 0; i < euclidDists.size(); i++) {
+			for (int j = 0; j < i; j++) {
+				averageEuclid += euclidDists[i][j];  
+				count++;
+			}
+		}
+		averageEuclid = averageEuclid / (float) count;   
+			
+		//find average for - Y
+		count = 0;
+		float averageUser = 0.0; 
+		for (int i = 0; i < userDists.size(); i++) {
+			for (int j = 0; j < i; j++) {
+				averageUser += userDists[i][j]; 
+				count++;
+			}
+		}
+		averageUser = averageUser / (float) count;  
+
+		double numerator = 0.0;
+		double denomTerm1 = 0.0;
+		double denomTerm2 = 0.0;
+		
+		for (int i = 0; i < euclidDists.size(); i++) {
+			
+			for (int k = 0; k < i; k++) { //just lt dists
+				
+				float Yi = userDists[i][k];
+				float Xi = euclidDists[i][k];
+				
+				numerator += ((Xi - averageEuclid) * (Yi - averageUser));
+				denomTerm1 += ((Xi - averageEuclid) * (Xi - averageEuclid));
+				denomTerm2 += ((Yi - averageUser) * (Yi - averageUser));
+			}
+		}
+		
+		double denom = (sqrt(denomTerm1) * sqrt(denomTerm2));
+		double r = numerator / denom;
+		
+		return r;
+		
+	}
+	catch(exception& e) {
+		m->errorOut(e, "LinearAlgebra", "calculateEuclidianDistance");
+		exit(1);
+	}
+}
+/*********************************************************************************************************************************/