aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Makefile2
-rw-r--r--src/debug.cpp33
-rw-r--r--src/gemma.cpp6
-rw-r--r--src/lapack.cpp2
-rw-r--r--src/lmm.cpp48
-rw-r--r--src/mvlmm.cpp5
6 files changed, 74 insertions, 22 deletions
diff --git a/Makefile b/Makefile
index e247a87..83a6845 100644
--- a/Makefile
+++ b/Makefile
@@ -198,7 +198,7 @@ $(OBJS): $(HDR)
./bin/unittests-gemma: contrib/catch-1.9.7/catch.hpp $(TEST_SRC_DIR)/unittests-main.o $(TEST_SRC_DIR)/unittests-math.o $(OBJS)
$(CPP) $(CPPFLAGS) $(TEST_SRC_DIR)/unittests-main.o $(TEST_SRC_DIR)/unittests-math.o $(filter-out src/main.o, $(OBJS)) $(LIBS) -o ./bin/unittests-gemma
-unittests: ./bin/unittests-gemma
+unittests: all ./bin/unittests-gemma
./bin/unittests-gemma
fast-check: all unittests
diff --git a/src/debug.cpp b/src/debug.cpp
index 8df4334..3efcce6 100644
--- a/src/debug.cpp
+++ b/src/debug.cpp
@@ -160,15 +160,22 @@ void write(const gsl_vector *v, const char *msg) {
void write(const gsl_matrix *m, const char *msg) {
if (msg) cout << "// " << msg << endl;
- cout << "// matrix size: " << m->size1 << " cols, " << m->size2 << " rows" << endl;
+ // Matrices are stored in row-major order, meaning that each row of
+ // elements forms a contiguous block in memory. This is the standard
+ // “C-language ordering” of two-dimensional arrays. The number of
+ // rows is size1.
+ auto rows = m->size1; // see https://www.gnu.org/software/gsl/manual/html_node/Matrices.html#Matrices
+ auto cols = m->size2;
+
+ cout << "// matrix size: " << cols << " cols, " << rows << " rows" << endl;
cout << "double " << msg << "[] = {";
- for (size_t j=0; j < m->size2; j++) {
- for (size_t i=0; i < m->size1; i++) {
+ for (size_t row=0; row < rows; row++) {
+ for (size_t col=0; col < cols; col++) {
// cout << "(" << i << "," << j << ")";
- cout << gsl_matrix_safe_get(m,i,j);
+ cout << gsl_matrix_safe_get(m,row,col);
cout << ",";
}
- cout << "// row " << j << endl;
+ cout << "// row " << row << endl;
}
cout << "}" << endl;
}
@@ -254,21 +261,21 @@ gsl_vector *gsl_vector_safe_alloc(size_t n) {
return v;
}
-double do_gsl_matrix_safe_get(const gsl_matrix * m, const size_t i, const size_t j,
+double do_gsl_matrix_safe_get(const gsl_matrix * m, const size_t row, const size_t col,
const char *__pretty_function, const char *__file, int __line) {
enforce(m);
if (!is_legacy_mode() && (is_debug_mode() || is_check_mode() || is_strict_mode())) {
- auto rows = m->size1;
+ auto rows = m->size1; // see above write function
auto cols = m->size2;
- if (i >= cols || j >= rows) {
- std::string msg = "Matrix out of bounds (" + std::to_string(cols) + "," + std::to_string(rows) + ") ";
- msg += std::to_string(i);
- msg += ",";
- msg += std::to_string(j);
+ if (col >= cols || row >= rows) {
+ std::string msg = "Matrix out of bounds (" + std::to_string(rows) + "," + std::to_string(cols) + ") ";
+ msg += std::to_string(row);
+ msg += "r,";
+ msg += std::to_string(col);
fail_at_msg(__file,__line,msg.c_str());
}
}
- return gsl_matrix_get(m,i,j);
+ return gsl_matrix_get(m,row,col);
}
diff --git a/src/gemma.cpp b/src/gemma.cpp
index ddb937e..94d05dc 100644
--- a/src/gemma.cpp
+++ b/src/gemma.cpp
@@ -1759,6 +1759,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
cout << "Start Eigen-Decomposition..." << endl;
time_start = clock();
cPar.trace_G = EigenDecomp_Zeroed(G, U, eval, 0);
+ write(eval,"eval zeroed");
cPar.time_eigen = (clock() - time_start) / (double(CLOCKS_PER_SEC) * 60.0);
// calculate UtW and Uty
@@ -2603,6 +2604,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
} else {
cPar.trace_G = EigenDecomp_Zeroed(G, U, eval, 0);
}
+ write(eval,"eval");
if (!cPar.file_weight.empty()) {
double wi;
@@ -2642,6 +2644,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
}
cPar.trace_G /= (double)eval->size;
}
+ // write(eval,"eval2");
if (cPar.a_mode == 31) {
cPar.WriteMatrix(U, "eigenU");
@@ -2700,6 +2703,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
assert(!std::isnan(cPar.beta_mle_null.front()));
assert(!std::isnan(cPar.se_beta_mle_null.front()));
+ // the following functions do not modify eval
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,
@@ -2798,6 +2802,8 @@ void GEMMA::BatchRun(PARAM &cPar) {
// if (is_check_mode()) disable_segfpe(); // disable fast NaN checking
+ // write(eval,"eval3");
+
if (!cPar.file_bfile.empty()) {
if (cPar.file_gxe.empty()) {
cMvlmm.AnalyzePlink(U, eval, UtW, UtY);
diff --git a/src/lapack.cpp b/src/lapack.cpp
index 6e43bd9..165a82d 100644
--- a/src/lapack.cpp
+++ b/src/lapack.cpp
@@ -239,6 +239,8 @@ void lapack_eigen_symmv(gsl_matrix *A, gsl_vector *eval, gsl_matrix *evec,
double EigenDecomp(gsl_matrix *G, gsl_matrix *U, gsl_vector *eval,
const size_t flag_largematrix) {
lapack_eigen_symmv(G, eval, U, flag_largematrix);
+ assert(!has_nan(eval));
+ // write(eval,"eval");
// Calculate track_G=mean(diag(G)).
double d = 0.0;
diff --git a/src/lmm.cpp b/src/lmm.cpp
index bba2e28..b368c0b 100644
--- a/src/lmm.cpp
+++ b/src/lmm.cpp
@@ -270,6 +270,8 @@ $7 = 3
$8 = 6
Hi_eval [0..ind] x Uab [ind, n_index] x ab [n_index]
+
+Iterating through a dataset Hi_eval differs and Uab (last row)
*/
void CalcPab(const size_t n_cvt, const size_t e_mode, const gsl_vector *Hi_eval,
@@ -290,9 +292,13 @@ void CalcPab(const size_t n_cvt, const size_t e_mode, const gsl_vector *Hi_eval,
double p_ab = 0.0;
+ debug_msg("CalcPab");
write(Hi_eval,"Hi_eval");
write(Uab,"Uab");
write(ab,"ab");
+ assert(!has_nan(Hi_eval));
+ assert(!has_nan(Uab));
+ // assert(!has_nan(ab));
for (size_t p = 0; p <= n_cvt + 1; ++p) { // p walks phenotypes + covariates
for (size_t a = p + 1; a <= n_cvt + 2; ++a) { // a in p+1..rest
@@ -303,9 +309,10 @@ void CalcPab(const size_t n_cvt, const size_t e_mode, const gsl_vector *Hi_eval,
gsl_matrix_const_column(Uab, index_ab);
gsl_blas_ddot(Hi_eval, &Uab_col.vector, &p_ab);
if (e_mode != 0) {
- if (! is_legacy_mode()) assert(false); // disabling to see when it is used
+ if (! is_legacy_mode()) assert(false); // disabling to see when it is used; allow with legacy mode
p_ab = gsl_vector_get(ab, index_ab) - p_ab;
}
+ cout << p << "r," << index_ab << ":" << p_ab << endl;
gsl_matrix_set(Pab, 0, index_ab, p_ab);
} else {
size_t index_aw = GetabIndex(a, p, n_cvt);
@@ -325,12 +332,15 @@ void CalcPab(const size_t n_cvt, const size_t e_mode, const gsl_vector *Hi_eval,
if (ps_ww != 0.0)
p_ab = ps_ab - ps_aw * ps_bw / ps_ww;
+ cout << p << "r," << index_ab << ":" << p_ab << endl;
gsl_matrix_set(Pab, p, index_ab, p_ab);
}
}
}
}
write(Pab,"Pab");
+ if (is_strict_mode() && (has_nan(Uab) || has_nan(Pab) || has_nan(Hi_eval) || has_nan(ab)))
+ exit(2);
return;
}
@@ -341,6 +351,11 @@ void CalcPPab(const size_t n_cvt, const size_t e_mode,
double p2_ab;
double ps2_ab, ps_aw, ps_bw, ps_ww, ps2_aw, ps2_bw, ps2_ww;
+ debug_msg("CalcPPab");
+ write(HiHi_eval,"Hi_eval");
+ write(Uab,"Uab");
+ write(ab,"ab");
+
for (size_t p = 0; p <= n_cvt + 1; ++p) {
for (size_t a = p + 1; a <= n_cvt + 2; ++a) {
for (size_t b = a; b <= n_cvt + 2; ++b) {
@@ -374,6 +389,10 @@ void CalcPPab(const size_t n_cvt, const size_t e_mode,
}
}
}
+ write(Pab,"Pab");
+ write(PPab,"PPab");
+ if (is_strict_mode() && (has_nan(Uab) || has_nan(PPab) || has_nan(HiHi_eval) || has_nan(ab)))
+ exit(2);
return;
}
@@ -385,6 +404,10 @@ void CalcPPPab(const size_t n_cvt, const size_t e_mode,
double p3_ab;
double ps3_ab, ps_aw, ps_bw, ps_ww, ps2_aw, ps2_bw, ps2_ww, ps3_aw, ps3_bw,
ps3_ww;
+ debug_msg("CalcPPPab");
+ write(HiHiHi_eval,"HiHiHi_eval");
+ write(Uab,"Uab");
+ write(ab,"ab");
for (size_t p = 0; p <= n_cvt + 1; ++p) {
for (size_t a = p + 1; a <= n_cvt + 2; ++a) {
@@ -428,6 +451,11 @@ void CalcPPPab(const size_t n_cvt, const size_t e_mode,
}
}
}
+ write(Pab,"Pab");
+ write(PPab,"PPab");
+ write(PPPab,"PPPab");
+ if (is_strict_mode() && (has_nan(Uab) || has_nan(PPPab) || has_nan(HiHiHi_eval) || has_nan(ab)))
+ exit(2);
return;
}
@@ -832,7 +860,8 @@ double LogRL_dev1(double l, void *params) {
gsl_vector *HiHi_eval = gsl_vector_safe_alloc((p->eval)->size);
gsl_vector *v_temp = gsl_vector_safe_alloc((p->eval)->size);
- gsl_vector_safe_memcpy(v_temp, p->eval);
+ // write(p->eval, "p->eval");
+ gsl_vector_safe_memcpy(v_temp, p->eval); // initialize with eval
gsl_vector_scale(v_temp, l);
if (p->e_mode == 0) {
gsl_vector_set_all(Hi_eval, 1.0);
@@ -852,6 +881,8 @@ double LogRL_dev1(double l, void *params) {
trace_Hi = (double)ni_test - trace_Hi;
}
+ write(p->eval, "p->eval2");
+ write(p->ab, "p->ab");
CalcPab(n_cvt, p->e_mode, Hi_eval, p->Uab, p->ab, Pab);
CalcPPab(n_cvt, p->e_mode, HiHi_eval, p->Uab, p->ab, Pab, PPab);
@@ -1216,12 +1247,9 @@ void CalcUab(const gsl_matrix *UtW, const gsl_vector *Uty,
return;
}
-/*
void Calcab(const gsl_matrix *W, const gsl_vector *y, gsl_vector *ab) {
- size_t index_ab;
size_t n_cvt = W->size2;
- double d;
gsl_vector *v_a = gsl_vector_safe_alloc(y->size);
gsl_vector *v_b = gsl_vector_safe_alloc(y->size);
@@ -1242,7 +1270,7 @@ void Calcab(const gsl_matrix *W, const gsl_vector *y, gsl_vector *ab) {
continue;
}
- index_ab = GetabIndex(a, b, n_cvt);
+ auto index_ab = GetabIndex(a, b, n_cvt);
if (b == n_cvt + 2) {
gsl_vector_safe_memcpy(v_b, y);
@@ -1251,6 +1279,7 @@ void Calcab(const gsl_matrix *W, const gsl_vector *y, gsl_vector *ab) {
gsl_vector_safe_memcpy(v_b, &W_col.vector);
}
+ double d;
gsl_blas_ddot(v_a, v_b, &d);
gsl_vector_set(ab, index_ab, d);
}
@@ -1288,7 +1317,6 @@ void Calcab(const gsl_matrix *W, const gsl_vector *y, const gsl_vector *x,
gsl_vector_safe_free(v_b);
return;
}
-*/
void LMM::AnalyzeGene(const gsl_matrix *U, const gsl_vector *eval,
const gsl_matrix *UtW, const gsl_vector *Utx,
@@ -1766,7 +1794,7 @@ void CalcLambda(const char func_name, FUNC_PARAM &params, const double l_min,
lambda_l = l_min * exp(lambda_interval * i);
lambda_h = l_min * exp(lambda_interval * (i + 1.0));
- if (func_name == 'R' || func_name == 'r') {
+ if (func_name == 'R' || func_name == 'r') { // log-restricted likelihood
dev1_l = LogRL_dev1(lambda_l, &params);
dev1_h = LogRL_dev1(lambda_h, &params);
} else {
@@ -1909,6 +1937,9 @@ void CalcLambda(const char func_name, const gsl_vector *eval,
const gsl_matrix *UtW, const gsl_vector *Uty,
const double l_min, const double l_max, const size_t n_region,
double &lambda, double &logl_H0) {
+
+ write(eval,"eval6");
+
if (func_name != 'R' && func_name != 'L' && func_name != 'r' &&
func_name != 'l') {
cout << "func_name only takes 'R' or 'L': 'R' for "
@@ -1924,6 +1955,7 @@ void CalcLambda(const char func_name, const gsl_vector *eval,
gsl_matrix_set_zero(Uab);
CalcUab(UtW, Uty, Uab);
+ Calcab(UtW, Uty, ab);
FUNC_PARAM param0 = {true, ni_test, n_cvt, eval, Uab, ab, 0};
diff --git a/src/mvlmm.cpp b/src/mvlmm.cpp
index 14f8b4a..0a259a7 100644
--- a/src/mvlmm.cpp
+++ b/src/mvlmm.cpp
@@ -2754,6 +2754,7 @@ void MphInitial(const size_t em_iter, const double em_prec,
const size_t n_region, gsl_matrix *V_g, gsl_matrix *V_e,
gsl_matrix *B) {
+ debug_msg("MphInitial");
gsl_matrix_set_zero(V_g);
gsl_matrix_set_zero(V_e);
gsl_matrix_set_zero(B);
@@ -2816,6 +2817,8 @@ void MphInitial(const size_t em_iter, const double em_prec,
gsl_matrix *Ve_sub = gsl_matrix_alloc(2, 2);
gsl_matrix *B_sub = gsl_matrix_alloc(2, c_size);
+ write(eval,"eval5");
+
for (size_t i = 0; i < d_size; i++) {
gsl_vector_view Y_sub1 = gsl_matrix_row(Y_sub, 0);
gsl_vector_const_view Y_1 = gsl_matrix_const_row(Y, i);
@@ -3477,6 +3480,8 @@ void MVLMM::AnalyzePlink(const gsl_matrix *U, const gsl_vector *eval,
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);
+ write(eval,"eval4");
+
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);
generated by cgit v1.2.3 (git 2.44.0) at 2025-01-18 23:50:48 +0000