PCA分析EasyPCA-1.0（c++ 0资源分）

/********************************************************************************** * EasyPCA performs a PCA analysis of a binary input file * * Copyright (C) 2008 Javier Cruz Mota * * E-Mail: javier.cruz (AT) epfl.ch * * Postal Address: EPFL / ENAC / INTER / TRANSP-OR * * GC B3 425 (Bâtiment GC), Station 18 * * CH-1015 Lausanne (Switzerland) * * * * This file is part of EasyPCA. * * * * EasyPCA is free software; you can redistribute it and/or modify * * it under the terms of the GNU General Public License as published by * * the Free Software Foundation; either version 2 of the License, or * * (at your option) any later version. * * * * EasyPCA is distributed in the hope that it will be useful, * * but WITHOUT ANY WARRANTY; without even the implied warranty of * * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * * GNU General Public License for more details. * * * * You should have received a copy of the GNU General Public License * * along with EasyPCA; if not, write to the Free Software * * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA * **********************************************************************************/ #include "pca.hpp" void PCA::PCA_LoadMatrix(void) { PCA_N = 4; PCA_M = 5; PCA_X = new double*[PCA_M]; for(int i=0; i<PCA_M; i++) PCA_X[i] = new double[PCA_N]; //Small matrix only for testing purposes PCA_X[0][0] = 4.0; PCA_X[0][1] = 1.0; PCA_X[0][2] = 7.0; PCA_X[0][3] = 1.0; PCA_X[1][0] = 1.0; PCA_X[1][1] = 3.0; PCA_X[1][2] = 1.0; PCA_X[1][3] = 0.0; PCA_X[2][0] = 1.0; PCA_X[2][1] = 0.0; PCA_X[2][2] = -5.0; PCA_X[2][3] = -3.0; PCA_X[3][0] = 0.0; PCA_X[3][1] = 2.0; PCA_X[3][2] = 1.0; PCA_X[3][3] = -7.0; PCA_X[4][0] = 0.0; PCA_X[4][1] = 3.0; PCA_X[4][2] = 1.0; PCA_X[4][3] = 4.0; PCA_IsTransposed = false; } void PCA::PCA_LoadMatrix(const char *PCA_File) { std::ifstream inputStream; inputStream.open(PCA_File); if( !inputStream ) { throw OpenERROR; } inputStream.read( (char *) &PCA_M, sizeof(int) ); inputStream.read( (char *) &PCA_N, sizeof(int) ); PCA_X = new double*[PCA_M]; for(int i=0; i<PCA_M; i++) PCA_X[i] = new double[PCA_N]; for(int i=0; i<PCA_M; i++) inputStream.read( (char *) &PCA_X[i][0], PCA_N*sizeof(double) ); PCA_IsTransposed = false; } void PCA::PCA_PerformPCA(void) { PCA_WriteX("XMatrix.bin"); std::cout << "\tEmpirical Mean..."; std::cout.flush(); PCA_EmpiricalMean(); std::cout << "Done!" << std::endl; std::cout << "\tSubstract Mean..."; std::cout.flush(); PCA_SubstractMean(); std::cout << "Done!" << std::endl; if (PCA_M > PCA_N) { std::cout << "\tTransposing to speed up PCA calculation..."; std::cout.flush(); PCA_Transpose(); std::cout << "Done!" << std::endl; } //PCA_WriteB("BMatrix.bin"); std::cout << "\tCovariance Matrix..."; std::cout.flush(); PCA_CovarianceMatrix(); std::cout << "Done!" << std::endl; //PCA_WriteC("covarMatrix.bin"); std::cout << "\tHouse-Holder Tridiagonalization..."; std::cout.flush(); PCA_HouseHolderTridiag(); //PCA_Print(); std::cout << "Done!" << std::endl; std::cout << "\tQL Decomposition..."; std::cout.flush(); PCA_QL(); PCA_QuickSortEV(); std::cout << "Done!" << std::endl; std::cout << "\tCumulative Energy..."; std::cout.flush(); PCA_CumulativeEnergy(); std::cout << "Done!" << std::endl; if (PCA_IsTransposed) { std::cout << "\tDe-Transposing matrix..."; std::cout.flush(); PCA_DeTranspose(); std::cout << "Done!" << std::endl; } } void PCA::PCA_EmpiricalMean(void) { PCA_u = new double[PCA_M]; for(int i=0; i<PCA_M; i++) PCA_u[i] = 0.0; for(int i=0; i<PCA_M; i++) { for(int j=0; j<PCA_N; j++) PCA_u[i]+=PCA_X[i][j]; PCA_u[i] = PCA_u[i]/PCA_N; } } void PCA::PCA_SubstractMean(void) { PCA_B = new double*[PCA_M]; for(int i=0; i<PCA_M; i++) PCA_B[i] = new double[PCA_N]; for(int i=0; i<PCA_M; i++) for(int j=0; j<PCA_N; j++) PCA_B[i][j] = PCA_X[i][j] - PCA_u[i]; } void PCA::PCA_CovarianceMatrix(void) { //Important: Assumming real matrices!! PCA_C = new double*[PCA_M]; for(int i=0; i<PCA_M; i++) PCA_C[i] = new double[PCA_M]; for(int i=0; i<PCA_M; i++) for(int j=i; j<PCA_M; j++) { PCA_C[i][j] = 0.0; for(int k=0; k<PCA_N; k++) PCA_C[i][j] += PCA_B[i][k]*PCA_B[j][k]; PCA_C[i][j] = PCA_C[i][j]/PCA_N; PCA_C[j][i] = PCA_C[i][j]; } } void PCA::PCA_HouseHolderTridiag(void) { if (DEBUG>0) std::cout << "PCA_HouseHolderTridiag:\tPCA_M = " << PCA_M << "\tPCA_N = " << PCA_N << std::endl; double scale,hh,h,g,f; int l,i,k,j; PCA_d = new double[PCA_M]; PCA_e = new double[PCA_M]; for(i=PCA_M-1; i>=1; i--) { if (DEBUG>0) std::cout << "PCA_HouseHolderTridiag:\tIteration " << PCA_M-i << std::endl; l = i-1; h = scale = 0.0; if (l>0) { for(k=0; k<=l; k++) scale += fabs(PCA_C[i][k]); if (scale == 0.0) PCA_e[i] = PCA_C[i][l]; else { for(k=0; k<=l; k++) { PCA_C[i][k] /= scale; h += PCA_C[i][k]*PCA_C[i][k]; } f = PCA_C[i][l]; g = (f>=0.0?-sqrt(h):sqrt(h)); PCA_e[i] = scale*g; h -= f*g; PCA_C[i][l] = f-g; f = 0.0; for(j=0; j<=l; j++) { PCA_C[j][i] = PCA_C[i][j]/h; g = 0.0; for(k=0; k<=j; k++) g += PCA_C[j][k]*PCA_C[i][k]; for(k=j+1; k<=l; k++) g += PCA_C[k][j]*PCA_C[i][k]; PCA_e[j] = g/h; f += PCA_e[j]*PCA_C[i][j]; } hh = f/(h+h); for(j=0; j<=l; j++) { f = PCA_C[i][j]; PCA_e[j] = g = PCA_e[j]-hh*f; for(k=0;k<=j;k++) PCA_C[j][k] -= (f*PCA_e[k]+g*PCA_C[i][k]); } } } else { PCA_e[i] = PCA_C[i][l]; } PCA_d[i] = h; } PCA_d[0] = 0.0; PCA_e[0] = 0.0; for(i=0; i<PCA_M; i++) { l = i-1; if (PCA_d[i]) { for(j=0; j<=l; j++) { g = 0.0; for(k=0; k<=l; k++) g += PCA_C[i][k]*PCA_C[k][j]; for(k=0; k<=l; k++) PCA_C[k][j] -= g*PCA_C[k][i]; } } PCA_d[i] = PCA_C[i][i]; PCA_C[i][i] = 1.0; for(j=0; j<=l; j++) { PCA_C[j][i] = PCA_C[i][j] = 0.0; } } } void PCA::PCA_QL(void) { int iter, m, i, l, k; double dd, g, r, s, c, p, f, b; for(i=1; i<PCA_M; i++) PCA_e[i-1] = PCA_e[i]; PCA_e[PCA_M-1] = 0.0; for(l=0; l<PCA_M; l++) { if (DEBUG>0) std::cout << "PCA_QL:\tIteration " << l << std::endl; iter=0; do{ for(m=l; m<PCA_M-1; m++) { dd = fabs(PCA_d[m]) + fabs(PCA_d[m+1]); if((double)(fabs(PCA_e[m])+dd) == dd) break; } if (DEBUG>0) std::cout << "PCA_QL:\t\tInner Iteration " << iter << std::endl; if(m!=l) { if(iter++ == 30) //Too many iterations { throw TooManyIterationERROR; } g = (PCA_d[l+1]-PCA_d[l])/(2.0*PCA_e[l]); r = pythag(g,1.0); g = PCA_d[m]-PCA_d[l]+PCA_e[l]/(g+SIGN(r,g)); s = c = 1.0; p = 0.0; for(i=m-1; i>=l; i--) { f = s*PCA_e[i]; b = c*PCA_e[i]; PCA_e[i+1] = (r=pythag(f,g)); if(r == 0.0) { PCA_d[i+1] -= p; PCA_e[m] = 0.0; break; } s = f/r; c = g/r; g = PCA_d[i+1]-p; r = (PCA_d[i]-g)*s+2.0*c*b; PCA_d[i+1] = g+(p=s*r)