From a5e0b3eb45a3d9bbf8760c49c9ac4f54817a2974 Mon Sep 17 00:00:00 2001 From: Pjotr Prins Date: Sun, 27 Aug 2017 07:07:41 +0000 Subject: Removed large sections of commented out code --- src/gemma.cpp | 234 +------------------------------------- src/lapack.cpp | 347 --------------------------------------------------------- 2 files changed, 1 insertion(+), 580 deletions(-) diff --git a/src/gemma.cpp b/src/gemma.cpp index d3401a6..5983461 100644 --- a/src/gemma.cpp +++ b/src/gemma.cpp @@ -2005,11 +2005,6 @@ void GEMMA::BatchRun(PARAM &cPar) { cPar.WriteMatrix(Vq, "Vq"); cPar.WriteVector(q, "q"); cPar.WriteVector(s, "size"); - /* - for (size_t i=0; isize1; (cPar.v_traceG).push_back(d); } - /* - //eigen-decomposition and calculate trace_G - cout<<"Start Eigen-Decomposition..."<size; i++) { - if (gsl_vector_get (eval, i)<1e-10) {gsl_vector_set (eval, i, 0);} - cPar.trace_G+=gsl_vector_get (eval, i); - } - cPar.trace_G/=(double)eval->size; - - cPar.time_eigen=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0); -} else { - ReadFile_eigenU (cPar.file_ku, cPar.error, U); - if (cPar.error==true) {cout<<"error! fail to read the U file. "<size; i++) { - if (gsl_vector_get(eval, i)<1e-10) {gsl_vector_set(eval, i, 0);} - cPar.trace_G+=gsl_vector_get(eval, i); - } - cPar.trace_G/=(double)eval->size; -} -*/ // fit multiple variance components if (cPar.n_ph == 1) { // if (cPar.n_vc==1) { - /* - //calculate UtW and Uty - CalcUtX (U, W, UtW); - CalcUtX (U, Y, UtY); - - gsl_vector_view beta=gsl_matrix_row (B, 0); - gsl_vector_view se_beta=gsl_matrix_row (se_B, 0); - gsl_vector_view UtY_col=gsl_matrix_column (UtY, 0); - - CalcLambda ('L', eval, UtW, &UtY_col.vector, cPar.l_min, cPar.l_max, - cPar.n_region, cPar.l_mle_null, cPar.logl_mle_H0); - CalcLmmVgVeBeta (eval, UtW, &UtY_col.vector, cPar.l_mle_null, - cPar.vg_mle_null, cPar.ve_mle_null, &beta.vector, &se_beta.vector); - - cPar.beta_mle_null.clear(); - cPar.se_beta_mle_null.clear(); - for (size_t i=0; isize2; i++) { - cPar.beta_mle_null.push_back(gsl_matrix_get(B, 0, i) ); - cPar.se_beta_mle_null.push_back(gsl_matrix_get(se_B, 0, i) ); - } - - CalcLambda ('R', eval, UtW, &UtY_col.vector, cPar.l_min, cPar.l_max, - cPar.n_region, cPar.l_remle_null, cPar.logl_remle_H0); - CalcLmmVgVeBeta (eval, UtW, &UtY_col.vector, cPar.l_remle_null, - cPar.vg_remle_null, cPar.ve_remle_null, &beta.vector, &se_beta.vector); - cPar.beta_remle_null.clear(); - cPar.se_beta_remle_null.clear(); - for (size_t i=0; isize2; i++) { - cPar.beta_remle_null.push_back(gsl_matrix_get(B, 0, i) ); - cPar.se_beta_remle_null.push_back(gsl_matrix_get(se_B, 0, i) ); - } - - CalcPve (eval, UtW, &UtY_col.vector, cPar.l_remle_null, cPar.trace_G, - cPar.pve_null, cPar.pve_se_null); - cPar.PrintSummary(); - - //calculate and output residuals - if (cPar.a_mode==5) { - gsl_vector *Utu_hat=gsl_vector_alloc (Y->size1); - gsl_vector *Ute_hat=gsl_vector_alloc (Y->size1); - gsl_vector *u_hat=gsl_vector_alloc (Y->size1); - gsl_vector *e_hat=gsl_vector_alloc (Y->size1); - gsl_vector *y_hat=gsl_vector_alloc (Y->size1); - - //obtain Utu and Ute - gsl_vector_memcpy (y_hat, &UtY_col.vector); - gsl_blas_dgemv (CblasNoTrans, -1.0, UtW, &beta.vector, 1.0, y_hat); - - double d, u, e; - for (size_t i=0; isize; i++) { - d=gsl_vector_get (eval, i); - u=cPar.l_remle_null*d/(cPar.l_remle_null*d+1.0)*gsl_vector_get(y_hat, - i); - e=1.0/(cPar.l_remle_null*d+1.0)*gsl_vector_get(y_hat, i); - gsl_vector_set (Utu_hat, i, u); - gsl_vector_set (Ute_hat, i, e); - } - - //obtain u and e - gsl_blas_dgemv (CblasNoTrans, 1.0, U, Utu_hat, 0.0, u_hat); - gsl_blas_dgemv (CblasNoTrans, 1.0, U, Ute_hat, 0.0, e_hat); - - //output residuals - cPar.WriteVector(u_hat, "residU"); - cPar.WriteVector(e_hat, "residE"); - - gsl_vector_free(u_hat); - gsl_vector_free(e_hat); - gsl_vector_free(y_hat); - } -*/ // } else { gsl_vector_view Y_col = gsl_matrix_column(Y, 0); VC cVc; @@ -2587,15 +2460,6 @@ return;} MFILEXwz(0, cPar.file_mgeno, cPar.d_pace, cPar.indicator_idv, cPar.mindicator_snp, vec_cat, w1, z, Xz); } - /* - cout<<"Xz: "<size, cPar.n_cvt); - gsl_matrix *G=gsl_matrix_alloc (1, 1); - vector > Xt; - - //set covariates matrix W and phenotype vector y - //an intercept is included in W - cPar.CopyCvtPhen (W, y, 0); - - //read in genotype matrix X - cPar.ReadGenotypes (Xt, G, false); - - LDR cLdr; - cLdr.CopyFromParam(cPar); - time_start=clock(); - - cLdr.VB(Xt, W, y); - - cPar.time_opt=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0); - cLdr.CopyToParam(cPar); - - gsl_vector_free (y); - gsl_matrix_free (W); - gsl_matrix_free (G); - } - */ - cPar.time_total = (clock() - time_begin) / (double(CLOCKS_PER_SEC) * 60.0); return; @@ -3514,18 +3324,6 @@ void GEMMA::WriteLog(int argc, char **argv, PARAM &cPar) { } outfile << endl; } - /* - outfile<<"## beta estimate in the null model = "; - for (size_t i=0; isize1, LDA = A->size1, INFO; - char UPLO = 'L'; - - if (N != (int)A->size2) { - cout << "Matrix needs to be symmetric and same dimension in " - << "lapack_cholesky_decomp." << endl; - return; - } - - spotrf_(&UPLO, &N, A->data, &LDA, &INFO); - if (INFO != 0) { - cout << "Cholesky decomposition unsuccessful in " - << "lapack_cholesky_decomp." << endl; - return; - } - - return; -} -*/ - // Cholesky decomposition, A is destroyed. void lapack_cholesky_decomp(gsl_matrix *A) { int N = A->size1, LDA = A->size1, INFO; @@ -105,30 +65,6 @@ void lapack_cholesky_decomp(gsl_matrix *A) { return; } -/* -// Cholesky solve, A is decomposed. -void lapack_float_cholesky_solve(gsl_matrix_float *A, const gsl_vector_float *b, - gsl_vector_float *x) { - int N = A->size1, NRHS = 1, LDA = A->size1, LDB = b->size, INFO; - char UPLO = 'L'; - - if (N != (int)A->size2 || N != LDB) { - cout << "Matrix needs to be symmetric and same dimension in " - << "lapack_cholesky_solve." << endl; - return; - } - - gsl_vector_float_memcpy(x, b); - spotrs_(&UPLO, &N, &NRHS, A->data, &LDA, x->data, &LDB, &INFO); - if (INFO != 0) { - cout << "Cholesky solve unsuccessful in lapack_cholesky_solve." << endl; - return; - } - - return; -} -*/ - // Cholesky solve, A is decomposed. void lapack_cholesky_solve(gsl_matrix *A, const gsl_vector *b, gsl_vector *x) { int N = A->size1, NRHS = 1, LDA = A->size1, LDB = b->size, INFO; @@ -150,61 +86,6 @@ void lapack_cholesky_solve(gsl_matrix *A, const gsl_vector *b, gsl_vector *x) { return; } -/* -void lapack_sgemm(char *TransA, char *TransB, float alpha, - const gsl_matrix_float *A, const gsl_matrix_float *B, - float beta, gsl_matrix_float *C) { - int M, N, K1, K2, LDA = A->size1, LDB = B->size1, LDC = C->size2; - - if (*TransA == 'N' || *TransA == 'n') { - M = A->size1; - K1 = A->size2; - } else if (*TransA == 'T' || *TransA == 't') { - M = A->size2; - K1 = A->size1; - } else { - cout << "need 'N' or 'T' in lapack_sgemm" << endl; - return; - } - - if (*TransB == 'N' || *TransB == 'n') { - N = B->size2; - K2 = B->size1; - } else if (*TransB == 'T' || *TransB == 't') { - N = B->size1; - K2 = B->size2; - } else { - cout << "need 'N' or 'T' in lapack_sgemm" << endl; - return; - } - - if (K1 != K2) { - cout << "A and B not compatible in lapack_sgemm" << endl; - return; - } - if (C->size1 != (size_t)M || C->size2 != (size_t)N) { - cout << "C not compatible in lapack_sgemm" << endl; - return; - } - - gsl_matrix_float *A_t = gsl_matrix_float_alloc(A->size2, A->size1); - gsl_matrix_float_transpose_memcpy(A_t, A); - gsl_matrix_float *B_t = gsl_matrix_float_alloc(B->size2, B->size1); - gsl_matrix_float_transpose_memcpy(B_t, B); - gsl_matrix_float *C_t = gsl_matrix_float_alloc(C->size2, C->size1); - gsl_matrix_float_transpose_memcpy(C_t, C); - - sgemm_(TransA, TransB, &M, &N, &K1, &alpha, A_t->data, &LDA, B_t->data, &LDB, - &beta, C_t->data, &LDC); - gsl_matrix_float_transpose_memcpy(C, C_t); - - gsl_matrix_float_free(A_t); - gsl_matrix_float_free(B_t); - gsl_matrix_float_free(C_t); - return; -} -*/ - void lapack_dgemm(char *TransA, char *TransB, double alpha, const gsl_matrix *A, const gsl_matrix *B, double beta, gsl_matrix *C) { int M, N, K1, K2, LDA = A->size1, LDB = B->size1, LDC = C->size2; @@ -258,91 +139,6 @@ void lapack_dgemm(char *TransA, char *TransB, double alpha, const gsl_matrix *A, return; } -/* -// Eigen value decomposition, matrix A is destroyed, float seems to -// have problem with large matrices (in mac). -void lapack_float_eigen_symmv(gsl_matrix_float *A, gsl_vector_float *eval, - gsl_matrix_float *evec, - const size_t flag_largematrix) { - if (flag_largematrix == 1) { - int N = A->size1, LDA = A->size1, INFO, LWORK = -1; - char JOBZ = 'V', UPLO = 'L'; - - if (N != (int)A->size2 || N != (int)eval->size) { - cout << "Matrix needs to be symmetric and same " - << "dimension in lapack_eigen_symmv." << endl; - return; - } - - LWORK = 3 * N; - float *WORK = new float[LWORK]; - ssyev_(&JOBZ, &UPLO, &N, A->data, &LDA, eval->data, WORK, &LWORK, &INFO); - if (INFO != 0) { - cout << "Eigen decomposition unsuccessful in " - << "lapack_eigen_symmv." << endl; - return; - } - - gsl_matrix_float_view A_sub = gsl_matrix_float_submatrix(A, 0, 0, N, N); - gsl_matrix_float_memcpy(evec, &A_sub.matrix); - gsl_matrix_float_transpose(evec); - - delete[] WORK; - } else { - int N = A->size1, LDA = A->size1, LDZ = A->size1, INFO, LWORK = -1, - LIWORK = -1; - char JOBZ = 'V', UPLO = 'L', RANGE = 'A'; - float ABSTOL = 1.0E-7; - - // VL, VU, IL, IU are not referenced; M equals N if RANGE='A'. - float VL = 0.0, VU = 0.0; - int IL = 0, IU = 0, M; - - if (N != (int)A->size2 || N != (int)eval->size) { - cout << "Matrix needs to be symmetric and same " - << "dimension in lapack_float_eigen_symmv." << endl; - return; - } - - int *ISUPPZ = new int[2 * N]; - - float WORK_temp[1]; - int IWORK_temp[1]; - ssyevr_(&JOBZ, &RANGE, &UPLO, &N, A->data, &LDA, &VL, &VU, &IL, &IU, - &ABSTOL, &M, eval->data, evec->data, &LDZ, ISUPPZ, WORK_temp, - &LWORK, IWORK_temp, &LIWORK, &INFO); - if (INFO != 0) { - cout << "Work space estimate unsuccessful in " - << "lapack_float_eigen_symmv." << endl; - return; - } - LWORK = (int)WORK_temp[0]; - LIWORK = (int)IWORK_temp[0]; - - float *WORK = new float[LWORK]; - int *IWORK = new int[LIWORK]; - - ssyevr_(&JOBZ, &RANGE, &UPLO, &N, A->data, &LDA, &VL, &VU, &IL, &IU, - &ABSTOL, &M, eval->data, evec->data, &LDZ, ISUPPZ, WORK, &LWORK, - IWORK, &LIWORK, &INFO); - if (INFO != 0) { - cout << "Eigen decomposition unsuccessful in " - << "lapack_float_eigen_symmv." << endl; - return; - } - - gsl_matrix_float_transpose(evec); - - delete[] ISUPPZ; - delete[] WORK; - delete[] IWORK; - } - - return; -} - -*/ - // Eigenvalue decomposition, matrix A is destroyed. Returns eigenvalues in // 'eval'. Also returns matrix 'evec' (U). void lapack_eigen_symmv(gsl_matrix *A, gsl_vector *eval, gsl_matrix *evec, @@ -473,23 +269,6 @@ double EigenDecomp_Zeroed(gsl_matrix *G, gsl_matrix *U, gsl_vector *eval, return d; } -/* -// DO NOT set eigen values to be positive. -double EigenDecomp(gsl_matrix_float *G, gsl_matrix_float *U, - gsl_vector_float *eval, const size_t flag_largematrix) { - lapack_float_eigen_symmv(G, eval, U, flag_largematrix); - - // Calculate track_G=mean(diag(G)). - double d = 0.0; - for (size_t i = 0; i < eval->size; ++i) { - d += gsl_vector_float_get(eval, i); - } - d /= (double)eval->size; - - return d; -} -*/ - double CholeskySolve(gsl_matrix *Omega, gsl_vector *Xty, gsl_vector *OiXty) { double logdet_O = 0.0; @@ -503,55 +282,12 @@ double CholeskySolve(gsl_matrix *Omega, gsl_vector *Xty, gsl_vector *OiXty) { return logdet_O; } -/* -double CholeskySolve(gsl_matrix_float *Omega, gsl_vector_float *Xty, - gsl_vector_float *OiXty) { - double logdet_O = 0.0; - - lapack_float_cholesky_decomp(Omega); - for (size_t i = 0; i < Omega->size1; ++i) { - logdet_O += log(gsl_matrix_float_get(Omega, i, i)); - } - logdet_O *= 2.0; - lapack_float_cholesky_solve(Omega, Xty, OiXty); - - return logdet_O; -} -*/ - // LU decomposition. void LUDecomp(gsl_matrix *LU, gsl_permutation *p, int *signum) { enforce_gsl(gsl_linalg_LU_decomp(LU, p, signum)); return; } -/* -void LUDecomp(gsl_matrix_float *LU, gsl_permutation *p, int *signum) { - gsl_matrix *LU_double = gsl_matrix_alloc(LU->size1, LU->size2); - - // Copy float matrix to double. - for (size_t i = 0; i < LU->size1; i++) { - for (size_t j = 0; j < LU->size2; j++) { - gsl_matrix_set(LU_double, i, j, gsl_matrix_float_get(LU, i, j)); - } - } - - // LU decomposition. - gsl_linalg_LU_decomp(LU_double, p, signum); - - // Copy float matrix to double. - for (size_t i = 0; i < LU->size1; i++) { - for (size_t j = 0; j < LU->size2; j++) { - gsl_matrix_float_set(LU, i, j, gsl_matrix_get(LU_double, i, j)); - } - } - - // Free matrix. - gsl_matrix_free(LU_double); - return; -} -*/ - // LU invert. Returns inverse. Note that GSL does not recommend using // this function void LUInvert(const gsl_matrix *LU, const gsl_permutation *p, gsl_matrix *ret_inverse) { @@ -561,37 +297,6 @@ void LUInvert(const gsl_matrix *LU, const gsl_permutation *p, gsl_matrix *ret_in enforce_gsl(gsl_linalg_LU_invert(LU, p, ret_inverse)); } -/* -void LUInvert(const gsl_matrix_float *LU, const gsl_permutation *p, - gsl_matrix_float *inverse) { - gsl_matrix *LU_double = gsl_matrix_alloc(LU->size1, LU->size2); - gsl_matrix *inverse_double = gsl_matrix_alloc(inverse->size1, inverse->size2); - - // Copy float matrix to double. - for (size_t i = 0; i < LU->size1; i++) { - for (size_t j = 0; j < LU->size2; j++) { - gsl_matrix_set(LU_double, i, j, gsl_matrix_float_get(LU, i, j)); - } - } - - // LU decomposition. - gsl_linalg_LU_invert(LU_double, p, inverse_double); - - // Copy float matrix to double. - for (size_t i = 0; i < inverse->size1; i++) { - for (size_t j = 0; j < inverse->size2; j++) { - gsl_matrix_float_set(inverse, i, j, gsl_matrix_get(inverse_double, i, j)); - } - } - - // Free matrix. - gsl_matrix_free(LU_double); - gsl_matrix_free(inverse_double); - return; -} - -*/ - // LU lndet. double LULndet(const gsl_matrix *LU) { return gsl_linalg_LU_lndet((gsl_matrix *)LU); @@ -625,43 +330,6 @@ void LUSolve(const gsl_matrix *LU, const gsl_permutation *p, return; } -/* -void LUSolve(const gsl_matrix_float *LU, const gsl_permutation *p, - const gsl_vector_float *b, gsl_vector_float *x) { - gsl_matrix *LU_double = gsl_matrix_alloc(LU->size1, LU->size2); - gsl_vector *b_double = gsl_vector_alloc(b->size); - gsl_vector *x_double = gsl_vector_alloc(x->size); - - // Copy float matrix to double. - for (size_t i = 0; i < LU->size1; i++) { - for (size_t j = 0; j < LU->size2; j++) { - gsl_matrix_set(LU_double, i, j, gsl_matrix_float_get(LU, i, j)); - } - } - - for (size_t i = 0; i < b->size; i++) { - gsl_vector_set(b_double, i, gsl_vector_float_get(b, i)); - } - - for (size_t i = 0; i < x->size; i++) { - gsl_vector_set(x_double, i, gsl_vector_float_get(x, i)); - } - - // LU decomposition. - gsl_linalg_LU_solve(LU_double, p, b_double, x_double); - - // Copy float matrix to double. - for (size_t i = 0; i < x->size; i++) { - gsl_vector_float_set(x, i, gsl_vector_get(x_double, i)); - } - - // Free matrix. - gsl_matrix_free(LU_double); - gsl_vector_free(b_double); - gsl_vector_free(x_double); - return; -} -*/ bool lapack_ddot(vector &x, vector &y, double &v) { bool flag = false; int incx = 1; @@ -674,18 +342,3 @@ bool lapack_ddot(vector &x, vector &y, double &v) { return flag; } - -/* -bool lapack_sdot(vector &x, vector &y, double &v) { - bool flag = false; - int incx = 1; - int incy = 1; - int n = (int)x.size(); - if (x.size() == y.size()) { - v = sdot_(&n, &x[0], &incx, &y[0], &incy); - flag = true; - } - - return flag; -} -*/ -- cgit v1.2.3