From 92ce6ce0f50f54c7cf4ba5488a423dcc0ded23fb Mon Sep 17 00:00:00 2001 From: Peter Carbonetto Date: Wed, 19 Jul 2017 15:02:25 -0500 Subject: Removed some commented out code from mvlmm.cpp. --- src/mvlmm.cpp | 2777 +++++++++++++++++++++++++++++---------------------------- 1 file changed, 1437 insertions(+), 1340 deletions(-) (limited to 'src/mvlmm.cpp') diff --git a/src/mvlmm.cpp b/src/mvlmm.cpp index 1c65dd1..5eb7e38 100644 --- a/src/mvlmm.cpp +++ b/src/mvlmm.cpp @@ -27,12 +27,12 @@ #include #include #include +#include #include "gsl/gsl_vector.h" #include "gsl/gsl_matrix.h" #include "gsl/gsl_linalg.h" #include "gsl/gsl_blas.h" - #include "gsl/gsl_cdf.h" #include "gsl/gsl_roots.h" #include "gsl/gsl_min.h" @@ -953,8 +953,9 @@ void CalcHiQi (const gsl_vector *eval, const gsl_matrix *X, mat_dd, 0.0, &Hi_k.matrix); } - // Calculate Qi, and multiply I\otimes UtVeh on both side - // and calculate logdet_Q, don't forget to substract c_size*logdet_Ve. + // Calculate Qi, and multiply I\o times UtVeh on both side and + // calculate logdet_Q, don't forget to substract + // c_size*logdet_Ve. logdet_Q=CalcQi (eval, D_l, X, Qi)-(double)c_size*logdet_Ve; for (size_t i=0; isize, d_size=Hiy->size1; - double delta, d_Hiy_i, d_Hiy_j, d_Hi_i1i2, d_Hi_i1j2, d_Hi_j1i2, d_Hi_j1j2; + double delta, d_Hiy_i, d_Hiy_j, d_Hi_i1i2, d_Hi_i1j2; + double d_Hi_j1i2, d_Hi_j1j2; for (size_t k=0; ksize1, d_size=Hi->size1; size_t v=GetIndex(i, j, d_size); double d; - //calculate the first part: trace(HiD) + // Calculate the first part: trace(HiD). Calc_traceHiD (eval, Hi, i, j, tPD_g, tPD_e); - //calculate the second part: -trace(HixQixHiD) + // Calculate the second part: -trace(HixQixHiD). for (size_t k=0; ksize1, d_size=Hi->size1; size_t v_size=d_size*(d_size+1)/2; size_t v1=GetIndex(i1, j1, d_size), v2=GetIndex(i2, j2, d_size); double d; - //calculate the first part: trace(HiDHiD) - Calc_traceHiDHiD (eval, Hi, i1, j1, i2, j2, tPDPD_gg, tPDPD_ee, tPDPD_ge); + // Calculate the first part: trace(HiDHiD). + Calc_traceHiDHiD (eval, Hi, i1, j1, i2, j2, tPDPD_gg, tPDPD_ee, + tPDPD_ge); - //calculate the second and third parts: -trace(HixQixHiDHiD)-trace(HiDHixQixHiD) + // Calculate the second and third parts: + // -trace(HixQixHiDHiD) - trace(HiDHixQixHiD) for (size_t i=0; isize2/eval->size, dc_size=xHi->size1; - size_t v; - - for (size_t i=0; isize2/eval->size, dc_size=xHi->size1; + size_t v; + + for (size_t i=0; isize; - double delta, d_Hi_ij; - - gsl_matrix *mat_dcdc=gsl_matrix_alloc (dc_size, dc_size); - gsl_matrix *mat_dcdc_temp=gsl_matrix_alloc (dc_size, dc_size); - - for (size_t k=0; ksize1; size_t v_size=xHiDHix_all_g->size2/dc_size; for (size_t i=0; isize2; i++) { - gsl_vector_const_view vec_g=gsl_matrix_const_column (vec_all_g, i); - gsl_vector_const_view vec_e=gsl_matrix_const_column (vec_all_e, i); - - gsl_vector_view Qivec_g=gsl_matrix_column (Qivec_all_g, i); - gsl_vector_view Qivec_e=gsl_matrix_column (Qivec_all_e, i); - - gsl_blas_dgemv (CblasNoTrans, 1.0, Qi, &vec_g.vector, 0.0, &Qivec_g.vector); - gsl_blas_dgemv (CblasNoTrans, 1.0, Qi, &vec_e.vector, 0.0, &Qivec_e.vector); + gsl_vector_const_view vec_g=gsl_matrix_const_column (vec_all_g, i); + gsl_vector_const_view vec_e=gsl_matrix_const_column (vec_all_e, i); + + gsl_vector_view Qivec_g=gsl_matrix_column (Qivec_all_g, i); + gsl_vector_view Qivec_e=gsl_matrix_column (Qivec_all_e, i); + + gsl_blas_dgemv(CblasNoTrans,1.0,Qi,&vec_g.vector,0.0, + &Qivec_g.vector); + gsl_blas_dgemv(CblasNoTrans,1.0,Qi,&vec_e.vector,0.0, + &Qivec_e.vector); } return; } - -//calculate Qi(xHiDHix) for each pair of i j (i<=j) -void Calc_QiMat_all (const gsl_matrix *Qi, const gsl_matrix *mat_all_g, const gsl_matrix *mat_all_e, gsl_matrix *Qimat_all_g, gsl_matrix *Qimat_all_e) -{ +// Calculate Qi(xHiDHix) for each pair of i,j (i<=j). +void Calc_QiMat_all (const gsl_matrix *Qi, const gsl_matrix *mat_all_g, + const gsl_matrix *mat_all_e, gsl_matrix *Qimat_all_g, + gsl_matrix *Qimat_all_e) { size_t dc_size=Qi->size1; size_t v_size=mat_all_g->size2/mat_all_g->size1; for (size_t i=0; isize1; size_t v=GetIndex(i, j, d_size); double d; - //first part: ytHiDHiy + // First part: ytHiDHiy. Calc_yHiDHiy (eval, Hiy, i, j, yPDPy_g, yPDPy_e); - //second and third parts: -(yHix)Qi(xHiDHiy)-(yHiDHix)Qi(xHiy) - gsl_vector_const_view xHiDHiy_g=gsl_matrix_const_column (xHiDHiy_all_g, v); - gsl_vector_const_view xHiDHiy_e=gsl_matrix_const_column (xHiDHiy_all_e, v); + // Second and third parts: -(yHix)Qi(xHiDHiy)-(yHiDHix)Qi(xHiy) + gsl_vector_const_view xHiDHiy_g = + gsl_matrix_const_column (xHiDHiy_all_g, v); + gsl_vector_const_view xHiDHiy_e = + gsl_matrix_const_column (xHiDHiy_all_e, v); gsl_blas_ddot(QixHiy, &xHiDHiy_g.vector, &d); yPDPy_g-=d*2.0; gsl_blas_ddot(QixHiy, &xHiDHiy_e.vector, &d); yPDPy_e-=d*2.0; - //fourth part: +(yHix)Qi(xHiDHix)Qi(xHiy) - gsl_vector_const_view xHiDHixQixHiy_g=gsl_matrix_const_column (xHiDHixQixHiy_all_g, v); - gsl_vector_const_view xHiDHixQixHiy_e=gsl_matrix_const_column (xHiDHixQixHiy_all_e, v); + // Fourth part: +(yHix)Qi(xHiDHix)Qi(xHiy). + gsl_vector_const_view xHiDHixQixHiy_g = + gsl_matrix_const_column (xHiDHixQixHiy_all_g, v); + gsl_vector_const_view xHiDHixQixHiy_e = + gsl_matrix_const_column (xHiDHixQixHiy_all_e, v); gsl_blas_ddot(QixHiy, &xHiDHixQixHiy_g.vector, &d); yPDPy_g+=d; @@ -1870,15 +1904,33 @@ void Calc_yPDPy (const gsl_vector *eval, const gsl_matrix *Hiy, const gsl_vector return; } -//calculate yPDPDPy=y(Hi-HixQixHi)D(Hi-HixQixHi)D(Hi-HixQixHi)y -//yPDPDPy=yHiDHiDHiy -//-(yHix)Qi(xHiDHiDHiy)-(yHiDHiDHix)Qi(xHiy) -//-(yHiDHix)Qi(xHiDHiy) -//+(yHix)Qi(xHiDHix)Qi(xHiDHiy)+(yHiDHix)Qi(xHiDHix)Qi(xHiy) -//+(yHix)Qi(xHiDHiDHix)Qi(xHiy) -//-(yHix)Qi(xHiDHix)Qi(xHiDHix)Qi(xHiy) -void Calc_yPDPDPy (const gsl_vector *eval, const gsl_matrix *Hi, const gsl_matrix *xHi, const gsl_matrix *Hiy, const gsl_vector *QixHiy, const gsl_matrix *xHiDHiy_all_g, const gsl_matrix *xHiDHiy_all_e, const gsl_matrix *QixHiDHiy_all_g, const gsl_matrix *QixHiDHiy_all_e, const gsl_matrix *xHiDHixQixHiy_all_g, const gsl_matrix *xHiDHixQixHiy_all_e, const gsl_matrix *QixHiDHixQixHiy_all_g, const gsl_matrix *QixHiDHixQixHiy_all_e, const gsl_matrix *xHiDHiDHiy_all_gg, const gsl_matrix *xHiDHiDHiy_all_ee, const gsl_matrix *xHiDHiDHiy_all_ge, const gsl_matrix *xHiDHiDHix_all_gg, const gsl_matrix *xHiDHiDHix_all_ee, const gsl_matrix *xHiDHiDHix_all_ge, const size_t i1, const size_t j1, const size_t i2, const size_t j2, double &yPDPDPy_gg, double &yPDPDPy_ee, double &yPDPDPy_ge) -{ +// calculate yPDPDPy = y(Hi-HixQixHi)D(Hi-HixQixHi)D(Hi-HixQixHi)y +// yPDPDPy = yHiDHiDHiy +// - (yHix)Qi(xHiDHiDHiy)-(yHiDHiDHix)Qi(xHiy) +// - (yHiDHix)Qi(xHiDHiy) +// + (yHix)Qi(xHiDHix)Qi(xHiDHiy) +// + (yHiDHix)Qi(xHiDHix)Qi(xHiy) +// + (yHix)Qi(xHiDHiDHix)Qi(xHiy) +// - (yHix)Qi(xHiDHix)Qi(xHiDHix)Qi(xHiy) +void Calc_yPDPDPy (const gsl_vector *eval, const gsl_matrix *Hi, + const gsl_matrix *xHi, const gsl_matrix *Hiy, + const gsl_vector *QixHiy, const gsl_matrix *xHiDHiy_all_g, + const gsl_matrix *xHiDHiy_all_e, + const gsl_matrix *QixHiDHiy_all_g, + const gsl_matrix *QixHiDHiy_all_e, + const gsl_matrix *xHiDHixQixHiy_all_g, + const gsl_matrix *xHiDHixQixHiy_all_e, + const gsl_matrix *QixHiDHixQixHiy_all_g, + const gsl_matrix *QixHiDHixQixHiy_all_e, + const gsl_matrix *xHiDHiDHiy_all_gg, + const gsl_matrix *xHiDHiDHiy_all_ee, + const gsl_matrix *xHiDHiDHiy_all_ge, + const gsl_matrix *xHiDHiDHix_all_gg, + const gsl_matrix *xHiDHiDHix_all_ee, + const gsl_matrix *xHiDHiDHix_all_ge, + const size_t i1, const size_t j1, const size_t i2, + const size_t j2, double &yPDPDPy_gg, double &yPDPDPy_ee, + double &yPDPDPy_ge) { size_t d_size=Hi->size1, dc_size=xHi->size1; size_t v1=GetIndex(i1, j1, d_size), v2=GetIndex(i2, j2, d_size); size_t v_size=d_size*(d_size+1)/2; @@ -1887,17 +1939,25 @@ void Calc_yPDPDPy (const gsl_vector *eval, const gsl_matrix *Hi, const gsl_matri gsl_vector *xHiDHiDHixQixHiy=gsl_vector_alloc (dc_size); - //first part: yHiDHiDHiy - Calc_yHiDHiDHiy (eval, Hi, Hiy, i1, j1, i2, j2, yPDPDPy_gg, yPDPDPy_ee, yPDPDPy_ge); - - //second and third parts: -(yHix)Qi(xHiDHiDHiy)-(yHiDHiDHix)Qi(xHiy) - gsl_vector_const_view xHiDHiDHiy_gg1=gsl_matrix_const_column (xHiDHiDHiy_all_gg, v1*v_size+v2); - gsl_vector_const_view xHiDHiDHiy_ee1=gsl_matrix_const_column (xHiDHiDHiy_all_ee, v1*v_size+v2); - gsl_vector_const_view xHiDHiDHiy_ge1=gsl_matrix_const_column (xHiDHiDHiy_all_ge, v1*v_size+v2); - - gsl_vector_const_view xHiDHiDHiy_gg2=gsl_matrix_const_column (xHiDHiDHiy_all_gg, v2*v_size+v1); - gsl_vector_const_view xHiDHiDHiy_ee2=gsl_matrix_const_column (xHiDHiDHiy_all_ee, v2*v_size+v1); - gsl_vector_const_view xHiDHiDHiy_ge2=gsl_matrix_const_column (xHiDHiDHiy_all_ge, v2*v_size+v1); + // First part: yHiDHiDHiy. + Calc_yHiDHiDHiy (eval, Hi, Hiy, i1, j1, i2, j2, yPDPDPy_gg, + yPDPDPy_ee, yPDPDPy_ge); + + // Second and third parts: + // -(yHix)Qi(xHiDHiDHiy) - (yHiDHiDHix)Qi(xHiy). + gsl_vector_const_view xHiDHiDHiy_gg1 = + gsl_matrix_const_column (xHiDHiDHiy_all_gg, v1*v_size+v2); + gsl_vector_const_view xHiDHiDHiy_ee1 = + gsl_matrix_const_column (xHiDHiDHiy_all_ee, v1*v_size+v2); + gsl_vector_const_view xHiDHiDHiy_ge1 = + gsl_matrix_const_column (xHiDHiDHiy_all_ge, v1*v_size+v2); + + gsl_vector_const_view xHiDHiDHiy_gg2 = + gsl_matrix_const_column (xHiDHiDHiy_all_gg, v2*v_size+v1); + gsl_vector_const_view xHiDHiDHiy_ee2 = + gsl_matrix_const_column (xHiDHiDHiy_all_ee, v2*v_size+v1); + gsl_vector_const_view xHiDHiDHiy_ge2 = + gsl_matrix_const_column (xHiDHiDHiy_all_ge, v2*v_size+v1); gsl_blas_ddot(QixHiy, &xHiDHiDHiy_gg1.vector, &d); yPDPDPy_gg-=d; @@ -1913,11 +1973,15 @@ void Calc_yPDPDPy (const gsl_vector *eval, const gsl_matrix *Hi, const gsl_matri gsl_blas_ddot(QixHiy, &xHiDHiDHiy_ge2.vector, &d); yPDPDPy_ge-=d; - //fourth part: -(yHiDHix)Qi(xHiDHiy) - gsl_vector_const_view xHiDHiy_g1=gsl_matrix_const_column (xHiDHiy_all_g, v1); - gsl_vector_const_view xHiDHiy_e1=gsl_matrix_const_column (xHiDHiy_all_e, v1); - gsl_vector_const_view QixHiDHiy_g2=gsl_matrix_const_column (QixHiDHiy_all_g, v2); - gsl_vector_const_view QixHiDHiy_e2=gsl_matrix_const_column (QixHiDHiy_all_e, v2); + // Fourth part: - (yHiDHix)Qi(xHiDHiy). + gsl_vector_const_view xHiDHiy_g1 = + gsl_matrix_const_column (xHiDHiy_all_g, v1); + gsl_vector_const_view xHiDHiy_e1 = + gsl_matrix_const_column (xHiDHiy_all_e, v1); + gsl_vector_const_view QixHiDHiy_g2 = + gsl_matrix_const_column (QixHiDHiy_all_g, v2); + gsl_vector_const_view QixHiDHiy_e2 = + gsl_matrix_const_column (QixHiDHiy_all_e, v2); gsl_blas_ddot(&xHiDHiy_g1.vector, &QixHiDHiy_g2.vector, &d); yPDPDPy_gg-=d; @@ -1926,14 +1990,22 @@ void Calc_yPDPDPy (const gsl_vector *eval, const gsl_matrix *Hi, const gsl_matri gsl_blas_ddot(&xHiDHiy_g1.vector, &QixHiDHiy_e2.vector, &d); yPDPDPy_ge-=d; - //fifth and sixth parts: +(yHix)Qi(xHiDHix)Qi(xHiDHiy)+(yHiDHix)Qi(xHiDHix)Qi(xHiy) - gsl_vector_const_view QixHiDHiy_g1=gsl_matrix_const_column (QixHiDHiy_all_g, v1); - gsl_vector_const_view QixHiDHiy_e1=gsl_matrix_const_column (QixHiDHiy_all_e, v1); - - gsl_vector_const_view xHiDHixQixHiy_g1=gsl_matrix_const_column (xHiDHixQixHiy_all_g, v1); - gsl_vector_const_view xHiDHixQixHiy_e1=gsl_matrix_const_column (xHiDHixQixHiy_all_e, v1); - gsl_vector_const_view xHiDHixQixHiy_g2=gsl_matrix_const_column (xHiDHixQixHiy_all_g, v2); - gsl_vector_const_view xHiDHixQixHiy_e2=gsl_matrix_const_column (xHiDHixQixHiy_all_e, v2); + // Fifth and sixth parts: + // + (yHix)Qi(xHiDHix)Qi(xHiDHiy) + + // (yHiDHix)Qi(xHiDHix)Qi(xHiy) + gsl_vector_const_view QixHiDHiy_g1 = + gsl_matrix_const_column (QixHiDHiy_all_g, v1); + gsl_vector_const_view QixHiDHiy_e1 = + gsl_matrix_const_column (QixHiDHiy_all_e, v1); + + gsl_vector_const_view xHiDHixQixHiy_g1 = + gsl_matrix_const_column (xHiDHixQixHiy_all_g, v1); + gsl_vector_const_view xHiDHixQixHiy_e1 = + gsl_matrix_const_column (xHiDHixQixHiy_all_e, v1); + gsl_vector_const_view xHiDHixQixHiy_g2 = + gsl_matrix_const_column (xHiDHixQixHiy_all_g, v2); + gsl_vector_const_view xHiDHixQixHiy_e2 = + gsl_matrix_const_column (xHiDHixQixHiy_all_e, v2); gsl_blas_ddot(&xHiDHixQixHiy_g1.vector, &QixHiDHiy_g2.vector, &d); yPDPDPy_gg+=d; @@ -1950,50 +2022,70 @@ void Calc_yPDPDPy (const gsl_vector *eval, const gsl_matrix *Hi, const gsl_matri gsl_blas_ddot(&xHiDHixQixHiy_e2.vector, &QixHiDHiy_g1.vector, &d); yPDPDPy_ge+=d; - //seventh part: +(yHix)Qi(xHiDHiDHix)Qi(xHiy) - gsl_matrix_const_view xHiDHiDHix_gg=gsl_matrix_const_submatrix (xHiDHiDHix_all_gg, 0, (v1*v_size+v2)*dc_size, dc_size, dc_size); - gsl_matrix_const_view xHiDHiDHix_ee=gsl_matrix_const_submatrix (xHiDHiDHix_all_ee, 0, (v1*v_size+v2)*dc_size, dc_size, dc_size); - gsl_matrix_const_view xHiDHiDHix_ge=gsl_matrix_const_submatrix (xHiDHiDHix_all_ge, 0, (v1*v_size+v2)*dc_size, dc_size, dc_size); - - gsl_blas_dgemv (CblasNoTrans, 1.0, &xHiDHiDHix_gg.matrix, QixHiy, 0.0, xHiDHiDHixQixHiy); + // Seventh part: + (yHix)Qi(xHiDHiDHix)Qi(xHiy) + gsl_matrix_const_view xHiDHiDHix_gg = + gsl_matrix_const_submatrix (xHiDHiDHix_all_gg, 0, + (v1*v_size+v2)*dc_size, + dc_size, dc_size); + gsl_matrix_const_view xHiDHiDHix_ee = + gsl_matrix_const_submatrix (xHiDHiDHix_all_ee, 0, + (v1*v_size+v2)*dc_size, + dc_size, dc_size); + gsl_matrix_const_view xHiDHiDHix_ge = + gsl_matrix_const_submatrix (xHiDHiDHix_all_ge, 0, + (v1*v_size+v2)*dc_size, + dc_size, dc_size); + + gsl_blas_dgemv (CblasNoTrans, 1.0, &xHiDHiDHix_gg.matrix, + QixHiy, 0.0, xHiDHiDHixQixHiy); gsl_blas_ddot(xHiDHiDHixQixHiy, QixHiy, &d); yPDPDPy_gg+=d; - gsl_blas_dgemv (CblasNoTrans, 1.0, &xHiDHiDHix_ee.matrix, QixHiy, 0.0, xHiDHiDHixQixHiy); + gsl_blas_dgemv (CblasNoTrans, 1.0, &xHiDHiDHix_ee.matrix, + QixHiy, 0.0, xHiDHiDHixQixHiy); gsl_blas_ddot(xHiDHiDHixQixHiy, QixHiy, &d); yPDPDPy_ee+=d; - gsl_blas_dgemv (CblasNoTrans, 1.0, &xHiDHiDHix_ge.matrix, QixHiy, 0.0, xHiDHiDHixQixHiy); + gsl_blas_dgemv (CblasNoTrans, 1.0, &xHiDHiDHix_ge.matrix, + QixHiy, 0.0, xHiDHiDHixQixHiy); gsl_blas_ddot(xHiDHiDHixQixHiy, QixHiy, &d); yPDPDPy_ge+=d; - //eighth part: -(yHix)Qi(xHiDHix)Qi(xHiDHix)Qi(xHiy) - gsl_vector_const_view QixHiDHixQixHiy_g1=gsl_matrix_const_column (QixHiDHixQixHiy_all_g, v1); - gsl_vector_const_view QixHiDHixQixHiy_e1=gsl_matrix_const_column (QixHiDHixQixHiy_all_e, v1); + // Eighth part: - (yHix)Qi(xHiDHix)Qi(xHiDHix)Qi(xHiy). + gsl_vector_const_view QixHiDHixQixHiy_g1 = + gsl_matrix_const_column (QixHiDHixQixHiy_all_g, v1); + gsl_vector_const_view QixHiDHixQixHiy_e1 = + gsl_matrix_const_column (QixHiDHixQixHiy_all_e, v1); - gsl_blas_ddot(&QixHiDHixQixHiy_g1.vector, &xHiDHixQixHiy_g2.vector, &d); + gsl_blas_ddot(&QixHiDHixQixHiy_g1.vector,&xHiDHixQixHiy_g2.vector,&d); yPDPDPy_gg-=d; - gsl_blas_ddot(&QixHiDHixQixHiy_e1.vector, &xHiDHixQixHiy_e2.vector, &d); + gsl_blas_ddot(&QixHiDHixQixHiy_e1.vector,&xHiDHixQixHiy_e2.vector,&d); yPDPDPy_ee-=d; - gsl_blas_ddot(&QixHiDHixQixHiy_g1.vector, &xHiDHixQixHiy_e2.vector, &d); + gsl_blas_ddot(&QixHiDHixQixHiy_g1.vector,&xHiDHixQixHiy_e2.vector,&d); yPDPDPy_ge-=d; - //free memory + // Free memory. gsl_vector_free(xHiDHiDHixQixHiy); return; } - -//calculate Edgeworth correctation factors for small samples -//notation and method follows Thomas J. Rothenberg, Econometirca 1984; 52 (4) -//M=xHiDHix -void CalcCRT (const gsl_matrix *Hessian_inv, const gsl_matrix *Qi, const gsl_matrix *QixHiDHix_all_g, const gsl_matrix *QixHiDHix_all_e, const gsl_matrix *xHiDHiDHix_all_gg, const gsl_matrix *xHiDHiDHix_all_ee, const gsl_matrix *xHiDHiDHix_all_ge, const size_t d_size, double &crt_a, double &crt_b, double &crt_c) -{ +// Calculate Edgeworth correctation factors for small samples notation +// and method follows Thomas J. Rothenberg, Econometirca 1984; 52 (4) +// M=xHiDHix +void CalcCRT (const gsl_matrix *Hessian_inv, const gsl_matrix *Qi, + const gsl_matrix *QixHiDHix_all_g, + const gsl_matrix *QixHiDHix_all_e, + const gsl_matrix *xHiDHiDHix_all_gg, + const gsl_matrix *xHiDHiDHix_all_ee, + const gsl_matrix *xHiDHiDHix_all_ge, + const size_t d_size, double &crt_a, + double &crt_b, double &crt_c) { crt_a=0.0; crt_b=0.0; crt_c=0.0; size_t dc_size=Qi->size1, v_size=Hessian_inv->size1/2; size_t c_size=dc_size/d_size; double h_gg, h_ge, h_ee, d, B=0.0, C=0.0, D=0.0; - double trCg1, trCe1, trCg2, trCe2, trB_gg, trB_ge, trB_ee, trCC_gg, trCC_ge, trCC_ee, trD_gg=0.0, trD_ge=0.0, trD_ee=0.0; + double trCg1, trCe1, trCg2, trCe2, trB_gg, trB_ge, trB_ee; + double trCC_gg, trCC_ge, trCC_ee, trD_gg=0.0, trD_ge=0.0, trD_ee=0.0; gsl_matrix *QiMQi_g1=gsl_matrix_alloc (dc_size, dc_size); gsl_matrix *QiMQi_e1=gsl_matrix_alloc (dc_size, dc_size); @@ -2018,10 +2110,12 @@ void CalcCRT (const gsl_matrix *Hessian_inv, const gsl_matrix *Qi, const gsl_mat gsl_matrix *M_dd=gsl_matrix_alloc (d_size, d_size); gsl_matrix *M_dcdc=gsl_matrix_alloc (dc_size, dc_size); - //invert Qi_sub to Qi_si + // Invert Qi_sub to Qi_si. gsl_matrix *Qi_sub=gsl_matrix_alloc (d_size, d_size); - gsl_matrix_const_view Qi_s=gsl_matrix_const_submatrix (Qi, (c_size-1)*d_size, (c_size-1)*d_size, d_size, d_size); + gsl_matrix_const_view Qi_s = + gsl_matrix_const_submatrix (Qi, (c_size-1)*d_size, + (c_size-1)*d_size, d_size, d_size); int sig; gsl_permutation * pmt=gsl_permutation_alloc (d_size); @@ -2033,195 +2127,244 @@ void CalcCRT (const gsl_matrix *Hessian_inv, const gsl_matrix *Qi, const gsl_mat gsl_permutation_free(pmt); gsl_matrix_free(Qi_sub); - //calculate correctation factors + // Calculate correction factors. for (size_t v1=0; v1size1, d_size=Hi->size1; size_t c_size=dc_size/d_size; @@ -2275,137 +2423,140 @@ void CalcDev (const char func_name, const gsl_vector *eval, const gsl_matrix *Qi gsl_matrix *QixHiDHixQixHiy_all_g=gsl_matrix_alloc (dc_size, v_size); gsl_matrix *QixHiDHixQixHiy_all_e=gsl_matrix_alloc (dc_size, v_size); - gsl_matrix *xHiDHiDHiy_all_gg=gsl_matrix_alloc (dc_size, v_size*v_size); - gsl_matrix *xHiDHiDHiy_all_ee=gsl_matrix_alloc (dc_size, v_size*v_size); - gsl_matrix *xHiDHiDHiy_all_ge=gsl_matrix_alloc (dc_size, v_size*v_size); - gsl_matrix *xHiDHiDHix_all_gg=gsl_matrix_alloc (dc_size, v_size*v_size*dc_size); - gsl_matrix *xHiDHiDHix_all_ee=gsl_matrix_alloc (dc_size, v_size*v_size*dc_size); - gsl_matrix *xHiDHiDHix_all_ge=gsl_matrix_alloc (dc_size, v_size*v_size*dc_size); - - //calculate xHiDHiy_all, xHiDHix_all and xHiDHixQixHiy_all + gsl_matrix *xHiDHiDHiy_all_gg = + gsl_matrix_alloc (dc_size, v_size*v_size); + gsl_matrix *xHiDHiDHiy_all_ee = + gsl_matrix_alloc (dc_size, v_size*v_size); + gsl_matrix *xHiDHiDHiy_all_ge = + gsl_matrix_alloc (dc_size, v_size*v_size); + gsl_matrix *xHiDHiDHix_all_gg = + gsl_matrix_alloc (dc_size, v_size*v_size*dc_size); + gsl_matrix *xHiDHiDHix_all_ee = + gsl_matrix_alloc (dc_size, v_size*v_size*dc_size); + gsl_matrix *xHiDHiDHix_all_ge = + gsl_matrix_alloc (dc_size, v_size*v_size*dc_size); + + // Calculate xHiDHiy_all, xHiDHix_all and xHiDHixQixHiy_all. Calc_xHiDHiy_all (eval, xHi, Hiy, xHiDHiy_all_g, xHiDHiy_all_e); Calc_xHiDHix_all (eval, xHi, xHiDHix_all_g, xHiDHix_all_e); - Calc_xHiDHixQixHiy_all (xHiDHix_all_g, xHiDHix_all_e, QixHiy, xHiDHixQixHiy_all_g, xHiDHixQixHiy_all_e); - - Calc_xHiDHiDHiy_all (v_size, eval, Hi, xHi, Hiy, xHiDHiDHiy_all_gg, xHiDHiDHiy_all_ee, xHiDHiDHiy_all_ge); - Calc_xHiDHiDHix_all (v_size, eval, Hi, xHi, xHiDHiDHix_all_gg, xHiDHiDHix_all_ee, xHiDHiDHix_all_ge); - - //calculate QixHiDHiy_all, QixHiDHix_all and QixHiDHixQixHiy_all - Calc_QiVec_all (Qi, xHiDHiy_all_g, xHiDHiy_all_e, QixHiDHiy_all_g, QixHiDHiy_all_e); - Calc_QiVec_all (Qi, xHiDHixQixHiy_all_g, xHiDHixQixHiy_all_e, QixHiDHixQixHiy_all_g, QixHiDHixQixHiy_all_e); - Calc_QiMat_all (Qi, xHiDHix_all_g, xHiDHix_all_e, QixHiDHix_all_g, QixHiDHix_all_e); - - double tHiD_g, tHiD_e, tPD_g, tPD_e, tHiDHiD_gg, tHiDHiD_ee, tHiDHiD_ge, tPDPD_gg, tPDPD_ee, tPDPD_ge; + Calc_xHiDHixQixHiy_all (xHiDHix_all_g, xHiDHix_all_e, QixHiy, + xHiDHixQixHiy_all_g, xHiDHixQixHiy_all_e); + + Calc_xHiDHiDHiy_all (v_size, eval, Hi, xHi, Hiy, xHiDHiDHiy_all_gg, + xHiDHiDHiy_all_ee, xHiDHiDHiy_all_ge); + Calc_xHiDHiDHix_all (v_size, eval, Hi, xHi, xHiDHiDHix_all_gg, + xHiDHiDHix_all_ee, xHiDHiDHix_all_ge); + + // Calculate QixHiDHiy_all, QixHiDHix_all and QixHiDHixQixHiy_all. + Calc_QiVec_all (Qi, xHiDHiy_all_g, xHiDHiy_all_e, QixHiDHiy_all_g, + QixHiDHiy_all_e); + Calc_QiVec_all (Qi, xHiDHixQixHiy_all_g, xHiDHixQixHiy_all_e, + QixHiDHixQixHiy_all_g, QixHiDHixQixHiy_all_e); + Calc_QiMat_all (Qi, xHiDHix_all_g, xHiDHix_all_e, QixHiDHix_all_g, + QixHiDHix_all_e); + + double tHiD_g, tHiD_e, tPD_g, tPD_e, tHiDHiD_gg, tHiDHiD_ee; + double tHiDHiD_ge, tPDPD_gg, tPDPD_ee, tPDPD_ge; double yPDPy_g, yPDPy_e, yPDPDPy_gg, yPDPDPy_ee, yPDPDPy_ge; - //calculate gradient and Hessian for Vg + // Calculate gradient and Hessian for Vg. for (size_t i1=0; i11) { - CalcCRT (Hessian_inv, Qi, QixHiDHix_all_g, QixHiDHix_all_e, xHiDHiDHix_all_gg, xHiDHiDHix_all_ee, xHiDHiDHix_all_ge, d_size, crt_a, crt_b, crt_c); + CalcCRT(Hessian_inv, Qi, QixHiDHix_all_g, QixHiDHix_all_e, + xHiDHiDHix_all_gg, xHiDHiDHix_all_ee, xHiDHiDHix_all_ge, + d_size, crt_a, crt_b, crt_c); } else { - crt_a=0.0; crt_b=0.0; crt_c=0.0; + crt_a=0.0; crt_b=0.0; crt_c=0.0; } - + gsl_matrix_free(xHiDHiy_all_g); gsl_matrix_free(xHiDHiy_all_e); gsl_matrix_free(xHiDHix_all_g); @@ -2430,10 +2581,9 @@ void CalcDev (const char func_name, const gsl_vector *eval, const gsl_matrix *Qi return; } - -//update Vg, Ve -void UpdateVgVe (const gsl_matrix *Hessian_inv, const gsl_vector *gradient, const double step_scale, gsl_matrix *V_g, gsl_matrix *V_e) -{ +// Update Vg, Ve. +void UpdateVgVe (const gsl_matrix *Hessian_inv, const gsl_vector *gradient, + const double step_scale, gsl_matrix *V_g, gsl_matrix *V_e) { size_t v_size=gradient->size/2, d_size=V_g->size1; size_t v; @@ -2441,7 +2591,7 @@ void UpdateVgVe (const gsl_matrix *Hessian_inv, const gsl_vector *gradient, cons double d; - //vectorize Vg and Ve + // Vectorize Vg and Ve. for (size_t i=0; isize, c_size=X->size1, d_size=Y->size1; size_t dc_size=d_size*c_size; size_t v_size=d_size*(d_size+1)/2; - double logdet_H, logdet_Q, yPy, logl_const, logl_old=0.0, logl_new=0.0, step_scale; + double logdet_H, logdet_Q, yPy, logl_const; + double logl_old=0.0, logl_new=0.0, step_scale; int sig; size_t step_iter, flag_pd; @@ -2506,80 +2662,84 @@ double MphNR (const char func_name, const size_t max_iter, const double max_prec gsl_vector *gradient=gsl_vector_alloc (v_size*2); - //calculate |XXt| and (XXt)^{-1} + // Calculate |XXt| and (XXt)^{-1}. gsl_blas_dsyrk (CblasUpper, CblasNoTrans, 1.0, X, 0.0, XXt); for (size_t i=0; i10 ) && step_iter<10 && t!=0); + gsl_matrix_memcpy (V_g, Vg_save); + gsl_matrix_memcpy (V_e, Ve_save); + + // Update Vg, Ve, and invert Hessian. + if (t!=0) { + UpdateVgVe (Hessian_inv, gradient, step_scale, V_g, V_e); + } + + // Check if both Vg and Ve are positive definite. + flag_pd=1; + gsl_matrix_memcpy (V_temp, V_e); + EigenDecomp(V_temp, U_temp, D_temp, 0); + for (size_t i=0; i10 ) && step_iter<10 && t!=0); - //terminate if change is small + // Terminate if change is small. if (t!=0) { if (logl_newsize1; i++) { - for (size_t j=0; jsize2; j++) { - cout<4) { - //first obtain good initial values - //large matrices for EM - gsl_matrix *U_hat=gsl_matrix_alloc (2, n_size); - gsl_matrix *E_hat=gsl_matrix_alloc (2, n_size); - gsl_matrix *OmegaU=gsl_matrix_alloc (2, n_size); - gsl_matrix *OmegaE=gsl_matrix_alloc (2, n_size); - gsl_matrix *UltVehiY=gsl_matrix_alloc (2, n_size); - gsl_matrix *UltVehiBX=gsl_matrix_alloc (2, n_size); - gsl_matrix *UltVehiU=gsl_matrix_alloc (2, n_size); - gsl_matrix *UltVehiE=gsl_matrix_alloc (2, n_size); - - //large matrices for NR - gsl_matrix *Hi_all=gsl_matrix_alloc (2, 2*n_size); //each dxd block is H_k^{-1} - gsl_matrix *Hiy_all=gsl_matrix_alloc (2, n_size); //each column is H_k^{-1}y_k - gsl_matrix *xHi_all=gsl_matrix_alloc (2*c_size, 2*n_size); //each dcxdc block is x_k\otimes H_k^{-1} - gsl_matrix *Hessian=gsl_matrix_alloc (6, 6); - - //2 by n matrix of Y - gsl_matrix *Y_sub=gsl_matrix_alloc (2, n_size); - gsl_matrix *Vg_sub=gsl_matrix_alloc (2, 2); - gsl_matrix *Ve_sub=gsl_matrix_alloc (2, 2); - gsl_matrix *B_sub=gsl_matrix_alloc (2, c_size); - - for (size_t i=0; i -void MVLMM::Analyzebgen (const gsl_matrix *U, const gsl_vector *eval, const gsl_matrix *UtW, const gsl_matrix *UtY) -{ +// WJA added. +void MVLMM::Analyzebgen (const gsl_matrix *U, const gsl_vector *eval, + const gsl_matrix *UtW, const gsl_matrix *UtY) { string file_bgen=file_oxford+".bgen"; ifstream infile (file_bgen.c_str(), ios::binary); - if (!infile) {cout<<"error reading bgen file:"<size1, msize); gsl_matrix *UtXlarge=gsl_matrix_alloc (U->size1, msize); gsl_matrix_set_zero(Xlarge); - // double lambda_mle=0, lambda_remle=0, beta=0, se=0, ; double logl_H0=0.0, logl_H1=0.0, p_wald=0, p_lrt=0, p_score=0; double crt_a, crt_b, crt_c; int n_miss, c_phen; double geno, x_mean; size_t c=0; - // double s=0.0; size_t n_size=UtY->size1, d_size=UtY->size2, c_size=UtW->size2; size_t dc_size=d_size*(c_size+1), v_size=d_size*(d_size+1)/2; - //large matrices for EM + // Large matrices for EM. gsl_matrix *U_hat=gsl_matrix_alloc (d_size, n_size); gsl_matrix *E_hat=gsl_matrix_alloc (d_size, n_size); gsl_matrix *OmegaU=gsl_matrix_alloc (d_size, n_size); @@ -2950,12 +3024,15 @@ void MVLMM::Analyzebgen (const gsl_matrix *U, const gsl_vector *eval, const gsl_ gsl_matrix *UltVehiU=gsl_matrix_alloc (d_size, n_size); gsl_matrix *UltVehiE=gsl_matrix_alloc (d_size, n_size); - //large matrices for NR - gsl_matrix *Hi_all=gsl_matrix_alloc (d_size, d_size*n_size); //each dxd block is H_k^{-1} - gsl_matrix *Hiy_all=gsl_matrix_alloc (d_size, n_size); //each column is H_k^{-1}y_k - gsl_matrix *xHi_all=gsl_matrix_alloc (dc_size, d_size*n_size); //each dcxdc block is x_k\otimes H_k^{-1} - gsl_matrix *Hessian=gsl_matrix_alloc (v_size*2, v_size*2); + // Large matrices for NR. Each dxd block is H_k^{-1}. + gsl_matrix *Hi_all=gsl_matrix_alloc (d_size, d_size*n_size); + // Each column is H_k^{-1}y_k. + gsl_matrix *Hiy_all=gsl_matrix_alloc (d_size, n_size); + + // Each dcxdc block is x_k\otimes H_k^{-1}. + gsl_matrix *xHi_all=gsl_matrix_alloc (dc_size, d_size*n_size); + gsl_matrix *Hessian=gsl_matrix_alloc (v_size*2, v_size*2); gsl_vector *x=gsl_vector_alloc (n_size); gsl_vector *x_miss=gsl_vector_alloc (n_size); @@ -2967,7 +3044,7 @@ void MVLMM::Analyzebgen (const gsl_matrix *U, const gsl_vector *eval, const gsl_ gsl_vector *beta=gsl_vector_alloc (d_size); gsl_matrix *Vbeta=gsl_matrix_alloc (d_size, d_size); - //null estimates for initial values + // Null estimates for initial values. gsl_matrix *V_g_null=gsl_matrix_alloc (d_size, d_size); gsl_matrix *V_e_null=gsl_matrix_alloc (d_size, d_size); gsl_matrix *B_null=gsl_matrix_alloc (d_size, c_size+1); @@ -2975,7 +3052,8 @@ void MVLMM::Analyzebgen (const gsl_matrix *U, const gsl_vector *eval, const gsl_ gsl_matrix_view X_sub=gsl_matrix_submatrix (X, 0, 0, c_size, n_size); gsl_matrix_view B_sub=gsl_matrix_submatrix (B, 0, 0, d_size, c_size); - gsl_matrix_view xHi_all_sub=gsl_matrix_submatrix (xHi_all, 0, 0, d_size*c_size, d_size*n_size); + gsl_matrix_view xHi_all_sub = + gsl_matrix_submatrix (xHi_all, 0, 0, d_size*c_size, d_size*n_size); gsl_matrix_transpose_memcpy (Y, UtY); @@ -2986,30 +3064,37 @@ void MVLMM::Analyzebgen (const gsl_matrix *U, const gsl_vector *eval, const gsl_ gsl_vector_view B_col=gsl_matrix_column(B, c_size); gsl_vector_set_zero(&B_col.vector); - MphInitial(em_iter, em_prec, nr_iter, nr_prec, eval, &X_sub.matrix, Y, l_min, l_max, n_region, V_g, V_e, &B_sub.matrix); - logl_H0=MphEM ('R', em_iter, em_prec, eval, &X_sub.matrix, Y, U_hat, E_hat, OmegaU, OmegaE, UltVehiY, UltVehiBX, UltVehiU, UltVehiE, V_g, V_e, &B_sub.matrix); - logl_H0=MphNR ('R', nr_iter, nr_prec, eval, &X_sub.matrix, Y, Hi_all, &xHi_all_sub.matrix, Hiy_all, V_g, V_e, Hessian, crt_a, crt_b, crt_c); - MphCalcBeta (eval, &X_sub.matrix, Y, V_g, V_e, UltVehiY, &B_sub.matrix, se_B_null); + MphInitial(em_iter, em_prec, nr_iter, nr_prec, eval, &X_sub.matrix, + Y, l_min, l_max, n_region, V_g, V_e, &B_sub.matrix); + logl_H0=MphEM ('R', em_iter, em_prec, eval, &X_sub.matrix, Y, U_hat, + E_hat, OmegaU, OmegaE, UltVehiY, UltVehiBX, UltVehiU, + UltVehiE, V_g, V_e, &B_sub.matrix); + logl_H0=MphNR ('R', nr_iter, nr_prec, eval, &X_sub.matrix, Y, + Hi_all, &xHi_all_sub.matrix, Hiy_all, V_g, V_e, + Hessian, crt_a, crt_b, crt_c); + MphCalcBeta (eval, &X_sub.matrix, Y, V_g, V_e, UltVehiY, + &B_sub.matrix, se_B_null); c=0; Vg_remle_null.clear(); Ve_remle_null.clear(); for (size_t i=0; isize1; i++) { - for (size_t j=0; jsize2; j++) { - beta_remle_null.push_back(gsl_matrix_get(B, i, j) ); - se_beta_remle_null.push_back(gsl_matrix_get(se_B_null, i, j) ); - } + for (size_t j=0; jsize2; j++) { + beta_remle_null.push_back(gsl_matrix_get(B, i, j) ); + se_beta_remle_null.push_back(gsl_matrix_get(se_B_null, i, j) ); + } } logl_remle_H0=logl_H0; @@ -3040,91 +3125,96 @@ void MVLMM::Analyzebgen (const gsl_matrix *U, const gsl_vector *eval, const gsl_ } cout<<"se(Ve): "<(&bgen_N),4); + infile.read(reinterpret_cast(&bgen_LS),2); + + id.resize(bgen_LS); + infile.read(&id[0], bgen_LS); + + infile.read(reinterpret_cast(&bgen_LR),2); + rs.resize(bgen_LR); + infile.read(&rs[0], bgen_LR); + + infile.read(reinterpret_cast(&bgen_LC),2); + chr.resize(bgen_LC); + infile.read(&chr[0], bgen_LC); + + infile.read(reinterpret_cast(&bgen_SNP_pos),4); + + infile.read(reinterpret_cast(&bgen_LA),4); + bgen_A_allele.resize(bgen_LA); + infile.read(&bgen_A_allele[0], bgen_LA); + + infile.read(reinterpret_cast(&bgen_LB),4); + bgen_B_allele.resize(bgen_LB); + infile.read(&bgen_B_allele[0], bgen_LB); + + uint16_t unzipped_data[3*bgen_N]; + + if (indicator_snp[t]==0) { + if(CompressedSNPBlocks) + infile.read(reinterpret_cast(&bgen_P),4); + else + bgen_P=6*bgen_N; + + infile.ignore(static_cast(bgen_P)); + + continue; + } + if(CompressedSNPBlocks) { -// if (t>1) {break;} - if (t%d_pace==0 || t==(ns_total-1)) {ProgressBar ("Reading SNPs ", t, ns_total-1);} - if (indicator_snp[t]==0) {continue;} - // read SNP header - id.clear(); - rs.clear(); - chr.clear(); - bgen_A_allele.clear(); - bgen_B_allele.clear(); - - infile.read(reinterpret_cast(&bgen_N),4); - infile.read(reinterpret_cast(&bgen_LS),2); - - id.resize(bgen_LS); - infile.read(&id[0], bgen_LS); - - infile.read(reinterpret_cast(&bgen_LR),2); - rs.resize(bgen_LR); - infile.read(&rs[0], bgen_LR); - - infile.read(reinterpret_cast(&bgen_LC),2); - chr.resize(bgen_LC); - infile.read(&chr[0], bgen_LC); - - infile.read(reinterpret_cast(&bgen_SNP_pos),4); - - infile.read(reinterpret_cast(&bgen_LA),4); - bgen_A_allele.resize(bgen_LA); - infile.read(&bgen_A_allele[0], bgen_LA); - - - infile.read(reinterpret_cast(&bgen_LB),4); - bgen_B_allele.resize(bgen_LB); - infile.read(&bgen_B_allele[0], bgen_LB); - - - - - uint16_t unzipped_data[3*bgen_N]; - - if (indicator_snp[t]==0) { - if(CompressedSNPBlocks) - infile.read(reinterpret_cast(&bgen_P),4); - else - bgen_P=6*bgen_N; - - infile.ignore(static_cast(bgen_P)); - - continue; - } - - - if(CompressedSNPBlocks) - { - - - infile.read(reinterpret_cast(&bgen_P),4); - uint8_t zipped_data[bgen_P]; - - unzipped_data_size=6*bgen_N; - - infile.read(reinterpret_cast(zipped_data),bgen_P); - - int result=uncompress(reinterpret_cast(unzipped_data), reinterpret_cast(&unzipped_data_size), reinterpret_cast(zipped_data), static_cast (bgen_P)); - assert(result == Z_OK); - - } - else - { - - bgen_P=6*bgen_N; - infile.read(reinterpret_cast(unzipped_data),bgen_P); - } - - x_mean=0.0; c_phen=0; n_miss=0; - gsl_vector_set_zero(x_miss); - for (size_t i=0; i(unzipped_data[i*3])/32768.0; - bgen_geno_prob_AB=static_cast(unzipped_data[i*3+1])/32768.0; - bgen_geno_prob_BB=static_cast(unzipped_data[i*3+2])/32768.0; - // WJA - bgen_geno_prob_non_miss=bgen_geno_prob_AA+bgen_geno_prob_AB+bgen_geno_prob_BB; - if (bgen_geno_prob_non_miss<0.9) {gsl_vector_set(x_miss, c_phen, 0.0); n_miss++;} - else { - - bgen_geno_prob_AA/=bgen_geno_prob_non_miss; - bgen_geno_prob_AB/=bgen_geno_prob_non_miss; - bgen_geno_prob_BB/=bgen_geno_prob_non_miss; + infile.read(reinterpret_cast(&bgen_P),4); + uint8_t zipped_data[bgen_P]; + + unzipped_data_size=6*bgen_N; + + infile.read(reinterpret_cast(zipped_data),bgen_P); + + int result=uncompress(reinterpret_cast(unzipped_data), + reinterpret_cast(&unzipped_data_size), + reinterpret_cast(zipped_data), + static_cast (bgen_P)); + assert(result == Z_OK); + + } else { + + bgen_P=6*bgen_N; + infile.read(reinterpret_cast(unzipped_data),bgen_P); + } + + x_mean=0.0; c_phen=0; n_miss=0; + gsl_vector_set_zero(x_miss); + for (size_t i=0; i(unzipped_data[i*3])/32768.0; + bgen_geno_prob_AB = + static_cast(unzipped_data[i*3+1])/32768.0; + bgen_geno_prob_BB = + static_cast(unzipped_data[i*3+2])/32768.0; + + // WJA. + bgen_geno_prob_non_miss=bgen_geno_prob_AA + + bgen_geno_prob_AB+bgen_geno_prob_BB; + if (bgen_geno_prob_non_miss<0.9) { + gsl_vector_set(x_miss, c_phen, 0.0); + n_miss++; + } + else { + + bgen_geno_prob_AA/=bgen_geno_prob_non_miss; + bgen_geno_prob_AB/=bgen_geno_prob_non_miss; + bgen_geno_prob_BB/=bgen_geno_prob_non_miss; + + geno=2.0*bgen_geno_prob_BB+bgen_geno_prob_AB; + + gsl_vector_set(x, c_phen, geno); + gsl_vector_set(x_miss, c_phen, 1.0); + x_mean+=geno; + } + c_phen++; + } + + x_mean/=static_cast(ni_test-n_miss); + + for (size_t i=0; isize1, l); + gsl_matrix_view UtXlarge_sub = + gsl_matrix_submatrix(UtXlarge, 0, 0, UtXlarge->size1, l); + + time_start=clock(); + eigenlib_dgemm ("T", "N", 1.0, U, &Xlarge_sub.matrix, 0.0, + &UtXlarge_sub.matrix); + time_UtX+=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0); + + gsl_matrix_set_zero (Xlarge); + + for (size_t i=0; i(ni_test-n_miss); - - for (size_t i=0; i1) { - //gsl_vector_set(x, i, 2-geno); - //} + } + + if (a_mode==1 || a_mode==4) { + logl_H1=MphEM ('R', em_iter/10, em_prec*10, eval, X, Y, U_hat, + E_hat, OmegaU, OmegaE, UltVehiY, UltVehiBX, + UltVehiU, UltVehiE, V_g, V_e, B); + p_wald=MphCalcP (eval, &X_row.vector, &X_sub.matrix, Y, V_g, + V_e, UltVehiY, beta, Vbeta); + + if (p_waldsize1, l); - gsl_matrix_view UtXlarge_sub=gsl_matrix_submatrix(UtXlarge, 0, 0, UtXlarge->size1, l); - - time_start=clock(); - eigenlib_dgemm ("T", "N", 1.0, U, &Xlarge_sub.matrix, 0.0, &UtXlarge_sub.matrix); - time_UtX+=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0); - - gsl_matrix_set_zero (Xlarge); - - for (size_t i=0; i1) {gsl_vector_scale(beta, -1.0);} - - time_opt+=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0); - - //store summary data - //SUMSTAT SNPs={snpInfo[t].get_chr(), snpInfo[t].get_rs(), snpInfo[t].get_pos(), n_miss, beta, se, lambda_remle, lambda_mle, p_wald, p_lrt, p_score}; - for (size_t i=0; isize1, d_size=UtY->size2, c_size=UtW->size2; size_t dc_size=d_size*(c_size+1), v_size=d_size*(d_size+1)/2; - //create a large matrix + // Create a large matrix. size_t msize=10000; gsl_matrix *Xlarge=gsl_matrix_alloc (U->size1, msize); gsl_matrix *UtXlarge=gsl_matrix_alloc (U->size1, msize); gsl_matrix_set_zero(Xlarge); - //large matrices for EM + // Large matrices for EM. gsl_matrix *U_hat=gsl_matrix_alloc (d_size, n_size); gsl_matrix *E_hat=gsl_matrix_alloc (d_size, n_size); gsl_matrix *OmegaU=gsl_matrix_alloc (d_size, n_size); -- cgit v1.2.3