integrated Rupali's changes - replaced OpenCV matrix math with LAPACK

David Mayerich
1 parent 89676102
Showing 1 changed file with 67 additions and 25 deletions Show diff stats
src/main.cu
@@ -50,6 +50,43 @@ typedef struct {
 }gnome;
 gnome gnom;
+//computing matrix determinant using LU decomposition
+template<typename T>
+T mtxdeterminant(T* M, size_t r, size_t c) {
+	//----DETERMINANT using LU decomposition using LAPACK
+	/* DGETRF computes an LU factorization of a general M-by-N matrix A using partial pivoting with row interchanges.
+	   The factorization has the form
+		A = P * L * U
+	   where P is a permutation matrix, L is lower triangular with unit diagonal elements (lower trapezoidal if m > n), and U is upper  triangular (upper trapezoidal if m < n).*/
+	int* ipiv = (int*)malloc(sizeof(int) * r);
+	memset(ipiv, 0, r * sizeof(int));
+	LAPACKE_sgetrf(LAPACK_COL_MAJOR, (int)r, (int)c, M, (int)r, ipiv);
+
+	//determinant of U
+	T product = 1;
+	for (size_t i = 0; i < r; i++) {
+		for (size_t j = 0; j < r; j++) {
+			if (i == j) {
+				product = product * M[i * r + j];
+			}
+		}
+	}
+
+	//determinant of P
+	int j;
+	double detp = 1.;
+	for (j = 0; j < r; j++) {
+		if (j + 1 != ipiv[j]) {
+			// j+1 : following feedback of ead : ipiv is from Fortran, hence starts at 1.
+			// hey ! This is a transpose !
+			detp = -detp;
+		}
+	}
+	T detSw = product * detp * 1; //det(U)*det(P)*det(L)
+	if (ipiv != NULL) std::free(ipiv);
+	return detSw;
+
+}
 void gpuComputeEignS( size_t g, size_t fea){	
 	//eigen value computation will return  r = (nC-1) eigen vectors so new projected data will have dimension of r rather than f
@@ -87,7 +124,7 @@ void gpuComputeEignS( size_t g, size_t fea){
 	int *IPIV = (int*) malloc(sizeof(int) * f);
 	//computing inverse of matrix Sw
 	memset(IPIV, 0, f * sizeof(int));
-	LAPACKE_sgetrf(LAPACK_COL_MAJOR, (int)f, (int)f, gSw_a, (int)f, IPIV);						
+	LAPACKE_sgetrf(LAPACK_COL_MAJOR, (int)f, (int)f, gSw_a, (int)f, IPIV);		
 	// DGETRI computes the inverse of a matrix using the LU factorization computed by DGETRF.
 	LAPACKE_sgetri(LAPACK_COL_MAJOR, (int)f, gSw_a, (int)f, IPIV);
@@ -146,23 +183,22 @@ void gpuComputeEignS( size_t g, size_t fea){
 	mtxMul(tempgSw, &gnom.lda[g * r * f], &gnom.Sw[g * f * f], r, f, f,f);
 	float* nSw = (float*)calloc(r * r, sizeof(float));
 	mtxMultranspose(nSw, tempgSw, &gnom.lda[g * r * f], r, f, r, f);
-	if(debug){
-		std::cout<<"From Eigen function: projected Sb sn Sw"<<std::endl;
+		
+	//determinant of Sb and Sw
+	float detSw = mtxdeterminant(nSw, r, r);
+	float detSb = mtxdeterminant(nSb, r, r);
+
+	if (debug) {
+		std::cout << "From Eigen function: projected Sb sn Sw" << std::endl;
 		displayS(nSb, r);									//display Sb				
 		displayS(nSw, r);									//display Sw
-		std::cout<<std::endl;
+		std::cout << std::endl;
+		std::cout <<"det(Sb)= "<< detSb<< std::endl;
+		std::cout << "det(Sw)= " << detSw << std::endl;
 	}
-	///DETERMINANT CURRENTLY REQUIRES OPENCV
-	std::cout<<"ERROR: This code requires a fix to remove an OpenCV dependence."<<std::endl;
-	mtxOutputFile("newSw.csv", nSw, r, r);
-	exit(1);
-	//FIX BY REPLACING THE FOLLOWING THREE LINES OF CODE USING LAPACK
-	//cv::Mat newSw = cv::Mat((int)r, (int)r, CV_32FC1, nSw);		//within scatter of gnome g in the population 
-	//cv::Mat newSb = cv::Mat((int)r, (int)r, CV_32FC1, nSb);		//within scatter of gnome g in the population 
-			
 	//fisher's ratio from ratio of projected sb and sw
-	float fisherRatio = 0;// = cv::determinant(newSb) /cv::determinant(newSw);
+	float fisherRatio = detSb /detSw;
 	gnom.S[g] = 1/fisherRatio;
 	if (debug) {
 		std::cout<<"Score["<<g<<"]: "<< gnom.S[g]<<std::endl;
@@ -172,6 +208,7 @@ void gpuComputeEignS( size_t g, size_t fea){
 			std::cout << ga.P[ga.f * g + i] << ", ";
 		std::cout << std::endl;
 	}
+
 	if(IPIV!= NULL) std::free(IPIV);	
@@ -430,27 +467,30 @@ int main(int argc, char* argv[]){
  //fitness values of best gnome is 
-	csv<<"best gnome's fitness value is "<<best_S[ga.gnrtn-1]<<std::endl;									//output fitness value of best gnome in last generation
-	//output gnome i.e. band index of selected featurs
+	csv<<"fitness value for best gnome:  "<<best_S[ga.gnrtn-1]<<std::endl;									//output fitness value of best gnome in last generation
+	/*//output gnome i.e. band index of selected featurs
 	csv<<(bestgnome.at(0));									//output feature subset
 	for(size_t i = 1; i < ga.f; i++)
 		csv<<","<<(bestgnome.at(i));
-	csv<<std::endl;
+	csv<<std::endl;*/
-	//output LDA basis of size r X f, r is nC - 1 as LDA projection is rank limited by number of classes - 1
-	for (size_t i = 0; i < nC-1; i++){						
-		csv<<lda[ldaindx + i * ga.f ];
-		for (size_t j = 1; j < ga.f; j++){
-			csv<<","<<lda[ldaindx + i * ga.f +j];
-		}
-		csv << std::endl;
-	}
 	//output actual wavelenths corresponding to those band indices 
+	csv << "selected wavenumbers:"<<std::endl;
 	csv << (wavelengths[bestgnome.at(0)]);
 	for (size_t i = 1; i < ga.f; i++)
-		csv << "," << (wavelengths[bestgnome.at(i)]);
+		csv << ", " << (wavelengths[bestgnome.at(i)]);
 	csv << std::endl;
+	//output LDA basis of size r X f, r is nC - 1 as LDA projection is rank limited by number of classes - 1
+	csv << "LDA basis:" << std::endl;
+	for (size_t i = 0; i < nC - 1; i++) {
+		csv << lda[ldaindx + i * ga.f];
+		for (size_t j = 1; j < ga.f; j++) {
+			csv << ", " << lda[ldaindx + i * ga.f + j];
+		}
+		csv << std::endl;
+	}
+
 	if (args["trim"].is_set()) {
 		csv << "trim info" << std::endl;
@@ -489,6 +529,8 @@ int main(int argc, char* argv[]){
 		////output actual wavelenths corresponding to those trim indices 
 		////these wavenumber are grouped in pairs, check each pair, if duplicated numbers are there in pair delete one and keep other as band to trim, if 2nd wavnumber is smaller than 1st in a pair ignore that pair
 		////e.g [1, 228, 230, 230, 232, 350,352, 351, 353, 1200] pairas [1(start)-228,230-230, 232-350, 352-351, 353-1200(end)], trimming wavenumbers are [1-228, 230, 233-350, 353-1200]
+		
+		
 		csv << (wavelengths[finaltrim_ind.at(0)]);
 		for (size_t i = 1; i < ga.f * + 2 ; i++)
 			csv << "," << (wavelengths[finaltrim_ind.at(i)]);