aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/debug.cpp55
-rw-r--r--src/debug.h19
-rw-r--r--src/gemma.cpp188
-rw-r--r--src/lmm.cpp288
-rw-r--r--src/mathfunc.cpp69
-rw-r--r--src/mathfunc.h4
6 files changed, 357 insertions, 266 deletions
diff --git a/src/debug.cpp b/src/debug.cpp
index d133344..a00190d 100644
--- a/src/debug.cpp
+++ b/src/debug.cpp
@@ -56,6 +56,55 @@ gsl_matrix *gsl_matrix_safe_alloc(size_t rows,size_t cols) {
return m;
}
+int gsl_matrix_safe_memcpy (gsl_matrix *dest, const gsl_matrix *src) {
+ enforce(dest->size1 == src->size1);
+ enforce(dest->size2 == src->size2);
+ return gsl_matrix_memcpy(dest,src);
+}
+
+void do_gsl_matrix_safe_free (gsl_matrix *m, const char *__pretty_function, const char *__file, int __line) {
+ enforce(m);
+ if (is_check_mode() && is_debug_mode()) {
+ bool has_NaN = has_nan(m);
+ bool has_Inf = has_inf(m);
+ if (has_NaN || has_Inf) {
+ std::string msg = "Matrix (size ";
+ msg += std::to_string(m->size1);
+ msg += "x";
+ msg += std::to_string(m->size2);
+ msg += ")";
+ if (has_Inf)
+ warnfail_at_msg(is_strict_mode(),__pretty_function,__file,__line,(msg+" contains Infinite on free!").c_str());
+ if (has_NaN)
+ warnfail_at_msg(is_strict_mode(),__pretty_function,__file,__line,(msg+" contains NaN on free!").c_str());
+ }
+ }
+ return gsl_matrix_free(m);
+}
+
+int gsl_vector_safe_memcpy (gsl_vector *dest, const gsl_vector *src) {
+ enforce(dest->size == src->size);
+ return gsl_vector_memcpy(dest,src);
+}
+
+void do_gsl_vector_safe_free (gsl_vector *v, const char *__pretty_function, const char *__file, int __line) {
+ enforce(v);
+ if (is_check_mode() && is_debug_mode()) {
+ bool has_NaN = has_nan(v);
+ bool has_Inf = has_inf(v);
+ if (has_NaN || has_Inf) {
+ std::string msg = "Vector (size ";
+ msg += std::to_string(v->size);
+ msg += ")";
+ if (has_Inf)
+ warnfail_at_msg(is_strict_mode(),__pretty_function,__file,__line,(msg+" contains Infinite on free!").c_str());
+ if (has_NaN)
+ warnfail_at_msg(is_strict_mode(),__pretty_function,__file,__line,(msg+" contains NaN on free!").c_str());
+ }
+ }
+ return gsl_vector_free(v);
+}
+
/*
Helper function to make sure gsl allocations do their job because
gsl_vector_alloc does not initiatize values (behaviour that changed
@@ -81,7 +130,7 @@ char *do_strtok_safe(char *tokenize, const char *delimiters, const char *__prett
// Helper function called by macro validate_K(K, check). K is validated
// unless -no-check option is used.
-void do_validate_K(const gsl_matrix *K, const char *__file, int __line) {
+void do_validate_K(const gsl_matrix *K, const char *__pretty_function, const char *__file, int __line) {
if (is_check_mode()) {
// debug_msg("Validating K");
auto eigenvalues = getEigenValues(K);
@@ -95,13 +144,13 @@ void do_validate_K(const gsl_matrix *K, const char *__file, int __line) {
if (!isMatrixIllConditioned(eigenvalues))
warning_at_msg(__file,__line,"K is ill conditioned!");
if (!isMatrixSymmetric(K))
- warnfail_at_msg(is_strict_mode(),__file,__line,"K is not symmetric!" );
+ warnfail_at_msg(is_strict_mode(),__pretty_function,__file,__line,"K is not symmetric!" );
const bool negative_values = has_negative_values_but_one(eigenvalues);
if (negative_values) {
warning_at_msg(__file,__line,"K has more than one negative eigenvalues!");
}
if (count_small>1 && negative_values && !isMatrixPositiveDefinite(K))
- warnfail_at_msg(is_strict_mode(),__file,__line,"K is not positive definite!");
+ warnfail_at_msg(is_strict_mode(),__pretty_function,__file,__line,"K is not positive definite!");
gsl_vector_free(eigenvalues);
}
}
diff --git a/src/debug.h b/src/debug.h
index 6b82b3b..69b0a7c 100644
--- a/src/debug.h
+++ b/src/debug.h
@@ -26,25 +26,32 @@ bool is_issue(uint issue);
bool is_legacy_mode();
gsl_matrix *gsl_matrix_safe_alloc(size_t rows,size_t cols);
+int gsl_matrix_safe_memcpy (gsl_matrix *dest, const gsl_matrix *src);
+void gsl_matrix_safe_free (gsl_matrix *v);
+void do_gsl_matrix_safe_free (gsl_matrix *m, const char *__pretty_function, const char *__file, int __line);
+
gsl_vector *gsl_vector_safe_alloc(size_t n);
+int gsl_vector_safe_memcpy (gsl_vector *dest, const gsl_vector *src);
+void gsl_vector_safe_free (gsl_vector *v);
+void do_gsl_vector_safe_free (gsl_vector *v, const char *__pretty_function, const char *__file, int __line);
char *do_strtok_safe(char *tokenize, const char *delimiters, const char *__pretty_function, const char *__file, int __line);
#define strtok_safe(string,delimiters) do_strtok_safe(string,delimiters,__PRETTY_FUNCTION__,__FILE__,__LINE__)
// Validation routines
-void do_validate_K(const gsl_matrix *K, const char *__file, int __line);
+void do_validate_K(const gsl_matrix *K, const char*__pretty_func, const char *__file, int __line);
#define ROUND(f) round(f * 10000.)/10000
-#define validate_K(K) do_validate_K(K,__FILE__,__LINE__)
+#define validate_K(K) do_validate_K(K,__PRETTY_FUNCTION__,__FILE__,__LINE__)
#define warning_at_msg(__file,__line,msg) cerr << "**** WARNING: " << msg << " in " << __file << " at line " << __line << endl;
-inline void warnfail_at_msg(bool strict, const char *__file, int __line, const char *msg) {
+inline void warnfail_at_msg(bool strict, const char *__function, const char *__file, int __line, const char *msg) {
if (strict)
std::cerr << "**** STRICT FAIL: ";
else
std::cerr << "**** WARNING: ";
- std::cerr << msg << " in " << __file << " at line " << __line << std::endl;
+ std::cerr << msg << " in " << __file << " at line " << __line << " in " << __function << std::endl;
if (strict)
exit(1);
}
@@ -127,5 +134,9 @@ inline void __enforce_fail(const char *__assertion, const char *__file,
enforce_msg(stat(fn.c_str(), &fileInfo) == 0, \
((std::string(__STRING(fn)) + " " + fn + ": " + msg).c_str()));
+#define gsl_matrix_safe_free(m) \
+ do_gsl_matrix_safe_free(m,__ASSERT_FUNCTION,__FILE__,__LINE__);
+#define gsl_vector_safe_free(v) \
+ do_gsl_vector_safe_free(v,__ASSERT_FUNCTION,__FILE__,__LINE__);
#endif
diff --git a/src/gemma.cpp b/src/gemma.cpp
index 2c54672..8977255 100644
--- a/src/gemma.cpp
+++ b/src/gemma.cpp
@@ -1673,8 +1673,8 @@ void GEMMA::BatchRun(PARAM &cPar) {
// read u
cPRDT.AddBV(G, u_hat, y_prdt);
- gsl_matrix_free(G);
- gsl_vector_free(u_hat);
+ gsl_matrix_safe_free(G);
+ gsl_vector_safe_free(u_hat);
}
// add beta
@@ -1701,7 +1701,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
cPRDT.WriteFiles(y_prdt);
- gsl_vector_free(y_prdt);
+ gsl_vector_safe_free(y_prdt);
}
// Prediction with kinship matrix only; for one or more phenotypes
@@ -1790,8 +1790,8 @@ void GEMMA::BatchRun(PARAM &cPar) {
gsl_matrix_add(H_full, G_full);
// free matrices
- gsl_vector_free(beta);
- gsl_vector_free(se_beta);
+ gsl_vector_safe_free(beta);
+ gsl_vector_safe_free(se_beta);
} else {
gsl_matrix *Vg = gsl_matrix_safe_alloc(cPar.n_ph, cPar.n_ph);
gsl_matrix *Ve = gsl_matrix_safe_alloc(cPar.n_ph, cPar.n_ph);
@@ -1838,10 +1838,10 @@ void GEMMA::BatchRun(PARAM &cPar) {
}
// free matrices
- gsl_matrix_free(Vg);
- gsl_matrix_free(Ve);
- gsl_matrix_free(B);
- gsl_matrix_free(se_B);
+ gsl_matrix_safe_free(Vg);
+ gsl_matrix_safe_free(Ve);
+ gsl_matrix_safe_free(B);
+ gsl_matrix_safe_free(se_B);
}
PRDT cPRDT;
@@ -1855,19 +1855,19 @@ void GEMMA::BatchRun(PARAM &cPar) {
cPRDT.WriteFiles(Y_full);
- gsl_matrix_free(Y);
- gsl_matrix_free(W);
- gsl_matrix_free(G);
- gsl_matrix_free(U);
- gsl_matrix_free(UtW);
- gsl_matrix_free(UtY);
- gsl_vector_free(eval);
-
- gsl_matrix_free(Y_full);
- gsl_matrix_free(Y_hat);
- gsl_matrix_free(W_full);
- gsl_matrix_free(G_full);
- gsl_matrix_free(H_full);
+ gsl_matrix_safe_free(Y);
+ gsl_matrix_safe_free(W);
+ gsl_matrix_safe_free(G);
+ gsl_matrix_safe_free(U);
+ gsl_matrix_safe_free(UtW);
+ gsl_matrix_safe_free(UtY);
+ gsl_vector_safe_free(eval);
+
+ gsl_matrix_safe_free(Y_full);
+ gsl_matrix_safe_free(Y_hat);
+ gsl_matrix_safe_free(W_full);
+ gsl_matrix_safe_free(G_full);
+ gsl_matrix_safe_free(H_full);
}
// Generate Kinship matrix (optionally using LOCO)
@@ -1895,7 +1895,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
cPar.WriteMatrix(G, "sXX");
}
- gsl_matrix_free(G);
+ gsl_matrix_safe_free(G);
}
// Compute the LDSC weights (not implemented yet)
@@ -1959,14 +1959,14 @@ void GEMMA::BatchRun(PARAM &cPar) {
cPar.WriteVector(ns, "size");
cPar.WriteVar("snps");
- gsl_matrix_free(S);
- gsl_vector_free(ns);
+ gsl_matrix_safe_free(S);
+ gsl_vector_safe_free(ns);
- gsl_matrix_free(A);
- gsl_matrix_free(K);
+ gsl_matrix_safe_free(A);
+ gsl_matrix_safe_free(K);
- gsl_vector_free(y);
- gsl_matrix_free(K);
+ gsl_vector_safe_free(y);
+ gsl_matrix_safe_free(K);
}
// Compute the q vector, that is used for variance component estimation using
@@ -2008,9 +2008,9 @@ void GEMMA::BatchRun(PARAM &cPar) {
cPar.WriteMatrix(Vq, "Vq");
cPar.WriteVector(q, "q");
cPar.WriteVector(s, "size");
- gsl_matrix_free(Vq);
- gsl_vector_free(q);
- gsl_vector_free(s);
+ gsl_matrix_safe_free(Vq);
+ gsl_vector_safe_free(q);
+ gsl_vector_safe_free(s);
}
// Calculate SNP covariance.
@@ -2057,8 +2057,8 @@ void GEMMA::BatchRun(PARAM &cPar) {
cLm.CopyToParam(cPar);
}
// release all matrices and vectors
- gsl_matrix_free(Y);
- gsl_matrix_free(W);
+ gsl_matrix_safe_free(Y);
+ gsl_matrix_safe_free(W);
}
// VC estimation with one or multiple kinship matrices
@@ -2186,15 +2186,15 @@ void GEMMA::BatchRun(PARAM &cPar) {
cPar.WriteVector(q, "q");
cPar.WriteVector(s, "size");
- gsl_matrix_free(S);
- gsl_matrix_free(Vq);
- gsl_vector_free(q);
- gsl_vector_free(s);
+ gsl_matrix_safe_free(S);
+ gsl_matrix_safe_free(Vq);
+ gsl_vector_safe_free(q);
+ gsl_vector_safe_free(s);
- gsl_matrix_free(A);
- gsl_matrix_free(K);
- gsl_vector_free(y);
- gsl_matrix_free(W);
+ gsl_matrix_safe_free(A);
+ gsl_matrix_safe_free(K);
+ gsl_vector_safe_free(y);
+ gsl_matrix_safe_free(W);
} else if (!cPar.file_study.empty() || !cPar.file_mstudy.empty()) {
if (!cPar.file_study.empty()) {
string sfile = cPar.file_study + ".size.txt";
@@ -2270,7 +2270,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
assert(!has_nan(cPar.v_se_pve));
gsl_vector_view s_sub = gsl_vector_subvector(s, 0, cPar.n_vc);
- gsl_vector_memcpy(&s_sub.vector, s_ref);
+ gsl_vector_safe_memcpy(&s_sub.vector, s_ref);
gsl_vector_set(s, cPar.n_vc, cPar.ni_ref);
cPar.WriteMatrix(S, "S");
@@ -2278,14 +2278,14 @@ void GEMMA::BatchRun(PARAM &cPar) {
cPar.WriteVector(q, "q");
cPar.WriteVector(s, "size");
- gsl_matrix_free(S);
- gsl_matrix_free(Vq);
- // gsl_matrix_free (V);
- // gsl_matrix_free (Vslope);
- gsl_vector_free(q);
- gsl_vector_free(s_study);
- gsl_vector_free(s_ref);
- gsl_vector_free(s);
+ gsl_matrix_safe_free(S);
+ gsl_matrix_safe_free(Vq);
+ // gsl_matrix_safe_free (V);
+ // gsl_matrix_safe_free (Vslope);
+ gsl_vector_safe_free(q);
+ gsl_vector_safe_free(s_study);
+ gsl_vector_safe_free(s_ref);
+ gsl_vector_safe_free(s);
} else {
gsl_matrix *Y = gsl_matrix_safe_alloc(cPar.ni_test, cPar.n_ph);
gsl_matrix *W = gsl_matrix_safe_alloc(Y->size1, cPar.n_cvt);
@@ -2462,7 +2462,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
cPar.mindicator_snp, vec_cat, w1, z, Xz);
}
if (cPar.a_mode == 66) {
- gsl_matrix_memcpy(XWz, Xz);
+ gsl_matrix_safe_memcpy(XWz, Xz);
} else if (cPar.a_mode == 67) {
cout << "Calculating XWz ... " << endl;
@@ -2507,17 +2507,17 @@ void GEMMA::BatchRun(PARAM &cPar) {
cPar.v_se_sigma2, cPar.v_enrich, cPar.v_se_enrich);
assert(!has_nan(cPar.v_se_pve));
- gsl_matrix_free(S);
- gsl_matrix_free(Svar);
- gsl_vector_free(s_ref);
-
- gsl_matrix_free(Xz);
- gsl_matrix_free(XWz);
- gsl_matrix_free(XtXWz);
- gsl_vector_free(w);
- gsl_vector_free(w1);
- gsl_vector_free(z);
- gsl_vector_free(s_vec);
+ gsl_matrix_safe_free(S);
+ gsl_matrix_safe_free(Svar);
+ gsl_vector_safe_free(s_ref);
+
+ gsl_matrix_safe_free(Xz);
+ gsl_matrix_safe_free(XWz);
+ gsl_matrix_safe_free(XtXWz);
+ gsl_vector_safe_free(w);
+ gsl_vector_safe_free(w1);
+ gsl_vector_safe_free(z);
+ gsl_vector_safe_free(s_vec);
}
// LMM or mvLMM or Eigen-Decomposition
@@ -2716,7 +2716,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
gsl_vector *y_hat = gsl_vector_safe_alloc(Y->size1);
// obtain Utu and Ute
- gsl_vector_memcpy(y_hat, &UtY_col.vector);
+ gsl_vector_safe_memcpy(y_hat, &UtY_col.vector);
gsl_blas_dgemv(CblasNoTrans, -1.0, UtW, &beta.vector, 1.0, y_hat);
double d, u, e;
@@ -2737,9 +2737,9 @@ void GEMMA::BatchRun(PARAM &cPar) {
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);
+ gsl_vector_safe_free(u_hat);
+ gsl_vector_safe_free(e_hat);
+ gsl_vector_safe_free(y_hat);
}
}
@@ -2802,16 +2802,16 @@ void GEMMA::BatchRun(PARAM &cPar) {
}
// release all matrices and vectors
- gsl_matrix_free(Y);
- gsl_matrix_free(W);
- gsl_matrix_free(B);
- gsl_matrix_free(se_B);
- gsl_matrix_free(G);
- gsl_matrix_free(U);
- gsl_matrix_free(UtW);
- gsl_matrix_free(UtY);
- gsl_vector_free(eval);
- gsl_vector_free(env);
+ gsl_matrix_safe_free(Y);
+ gsl_matrix_safe_free(W);
+ gsl_matrix_safe_free(B);
+ gsl_matrix_safe_free(se_B);
+ gsl_matrix_safe_free(G);
+ gsl_matrix_safe_free(U);
+ gsl_matrix_safe_free(UtW);
+ gsl_matrix_safe_free(UtY);
+ gsl_vector_safe_free(eval);
+ gsl_vector_safe_free(env);
}
// BSLMM
@@ -2912,15 +2912,15 @@ void GEMMA::BatchRun(PARAM &cPar) {
}
// release all matrices and vectors
- gsl_matrix_free(G);
- gsl_matrix_free(U);
- gsl_matrix_free(UtW);
- gsl_vector_free(eval);
- gsl_vector_free(Uty);
+ gsl_matrix_safe_free(G);
+ gsl_matrix_safe_free(U);
+ gsl_matrix_safe_free(UtW);
+ gsl_vector_safe_free(eval);
+ gsl_vector_safe_free(Uty);
}
- gsl_matrix_free(W);
- gsl_vector_free(y);
- gsl_matrix_free(UtX);
+ gsl_matrix_safe_free(W);
+ gsl_vector_safe_free(y);
+ gsl_matrix_safe_free(UtX);
}
// BSLMM-DAP
@@ -3016,16 +3016,16 @@ void GEMMA::BatchRun(PARAM &cPar) {
cBslmmDap.CopyToParam(cPar);
// release all matrices and vectors
- gsl_matrix_free(G);
- gsl_matrix_free(U);
- gsl_matrix_free(UtW);
- gsl_vector_free(eval);
- gsl_vector_free(Uty);
+ gsl_matrix_safe_free(G);
+ gsl_matrix_safe_free(U);
+ gsl_matrix_safe_free(UtW);
+ gsl_vector_safe_free(eval);
+ gsl_vector_safe_free(Uty);
}
- gsl_matrix_free(W);
- gsl_vector_free(y);
- gsl_matrix_free(UtX);
+ gsl_matrix_safe_free(W);
+ gsl_vector_safe_free(y);
+ gsl_matrix_safe_free(UtX);
} else if (cPar.a_mode == 15) {
// perform EM algorithm and estimate parameters
vector<string> vec_rs;
@@ -3074,7 +3074,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
cPar.time_opt = (clock() - time_start) / (double(CLOCKS_PER_SEC) * 60.0);
cBslmmDap.CopyToParam(cPar);
- gsl_matrix_free(Ac);
+ gsl_matrix_safe_free(Ac);
gsl_matrix_int_free(Ad);
gsl_vector_int_free(dlevel);
} else {
diff --git a/src/lmm.cpp b/src/lmm.cpp
index 683858d..1a2614b 100644
--- a/src/lmm.cpp
+++ b/src/lmm.cpp
@@ -367,12 +367,12 @@ double LogL_f(double l, void *params) {
gsl_vector *Hi_eval = gsl_vector_safe_alloc((p->eval)->size);
gsl_vector *v_temp = gsl_vector_safe_alloc((p->eval)->size);
- gsl_vector_memcpy(v_temp, p->eval);
+ gsl_vector_safe_memcpy(v_temp, p->eval);
gsl_vector_scale(v_temp, l);
if (p->e_mode == 0) {
gsl_vector_set_all(Hi_eval, 1.0);
} else {
- gsl_vector_memcpy(Hi_eval, v_temp);
+ gsl_vector_safe_memcpy(Hi_eval, v_temp);
}
gsl_vector_add_constant(v_temp, 1.0);
gsl_vector_div(Hi_eval, v_temp);
@@ -391,9 +391,9 @@ double LogL_f(double l, void *params) {
double P_yy = gsl_matrix_get(Pab, nc_total, index_yy);
f = c - 0.5 * logdet_h - 0.5 * (double)ni_test * log(P_yy);
- gsl_matrix_free(Pab);
- gsl_vector_free(Hi_eval);
- gsl_vector_free(v_temp);
+ gsl_matrix_safe_free(Pab); // FIXME
+ gsl_vector_safe_free(Hi_eval);
+ gsl_vector_safe_free(v_temp);
return f;
}
@@ -419,17 +419,17 @@ double LogL_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_memcpy(v_temp, p->eval);
+ gsl_vector_safe_memcpy(v_temp, p->eval);
gsl_vector_scale(v_temp, l);
if (p->e_mode == 0) {
gsl_vector_set_all(Hi_eval, 1.0);
} else {
- gsl_vector_memcpy(Hi_eval, v_temp);
+ gsl_vector_safe_memcpy(Hi_eval, v_temp);
}
gsl_vector_add_constant(v_temp, 1.0);
gsl_vector_div(Hi_eval, v_temp);
- gsl_vector_memcpy(HiHi_eval, Hi_eval);
+ gsl_vector_safe_memcpy(HiHi_eval, Hi_eval);
gsl_vector_mul(HiHi_eval, Hi_eval);
gsl_vector_set_all(v_temp, 1.0);
@@ -451,11 +451,11 @@ double LogL_dev1(double l, void *params) {
double yPKPy = (P_yy - PP_yy) / l;
dev1 = -0.5 * trace_HiK + 0.5 * (double)ni_test * yPKPy / P_yy;
- gsl_matrix_free(Pab);
- gsl_matrix_free(PPab);
- gsl_vector_free(Hi_eval);
- gsl_vector_free(HiHi_eval);
- gsl_vector_free(v_temp);
+ gsl_matrix_safe_free(Pab); // FIXME: may contain NaN
+ gsl_matrix_safe_free(PPab); // FIXME: may contain NaN
+ gsl_vector_safe_free(Hi_eval);
+ gsl_vector_safe_free(HiHi_eval);
+ gsl_vector_safe_free(v_temp);
return dev1;
}
@@ -484,19 +484,19 @@ double LogL_dev2(double l, void *params) {
gsl_vector *HiHiHi_eval = gsl_vector_safe_alloc((p->eval)->size);
gsl_vector *v_temp = gsl_vector_safe_alloc((p->eval)->size);
- gsl_vector_memcpy(v_temp, p->eval);
+ gsl_vector_safe_memcpy(v_temp, p->eval);
gsl_vector_scale(v_temp, l);
if (p->e_mode == 0) {
gsl_vector_set_all(Hi_eval, 1.0);
} else {
- gsl_vector_memcpy(Hi_eval, v_temp);
+ gsl_vector_safe_memcpy(Hi_eval, v_temp);
}
gsl_vector_add_constant(v_temp, 1.0);
gsl_vector_div(Hi_eval, v_temp);
- gsl_vector_memcpy(HiHi_eval, Hi_eval);
+ gsl_vector_safe_memcpy(HiHi_eval, Hi_eval);
gsl_vector_mul(HiHi_eval, Hi_eval);
- gsl_vector_memcpy(HiHiHi_eval, HiHi_eval);
+ gsl_vector_safe_memcpy(HiHiHi_eval, HiHi_eval);
gsl_vector_mul(HiHiHi_eval, Hi_eval);
gsl_vector_set_all(v_temp, 1.0);
@@ -526,13 +526,13 @@ double LogL_dev2(double l, void *params) {
0.5 * (double)ni_test * (2.0 * yPKPKPy * P_yy - yPKPy * yPKPy) /
(P_yy * P_yy);
- gsl_matrix_free(Pab);
- gsl_matrix_free(PPab);
- gsl_matrix_free(PPPab);
- gsl_vector_free(Hi_eval);
- gsl_vector_free(HiHi_eval);
- gsl_vector_free(HiHiHi_eval);
- gsl_vector_free(v_temp);
+ gsl_matrix_safe_free(Pab); // FIXME
+ gsl_matrix_safe_free(PPab);
+ gsl_matrix_safe_free(PPPab);
+ gsl_vector_safe_free(Hi_eval);
+ gsl_vector_safe_free(HiHi_eval);
+ gsl_vector_safe_free(HiHiHi_eval);
+ gsl_vector_safe_free(v_temp);
return dev2;
}
@@ -561,19 +561,19 @@ void LogL_dev12(double l, void *params, double *dev1, double *dev2) {
gsl_vector *HiHiHi_eval = gsl_vector_safe_alloc((p->eval)->size);
gsl_vector *v_temp = gsl_vector_safe_alloc((p->eval)->size);
- gsl_vector_memcpy(v_temp, p->eval);
+ gsl_vector_safe_memcpy(v_temp, p->eval);
gsl_vector_scale(v_temp, l);
if (p->e_mode == 0) {
gsl_vector_set_all(Hi_eval, 1.0);
} else {
- gsl_vector_memcpy(Hi_eval, v_temp);
+ gsl_vector_safe_memcpy(Hi_eval, v_temp);
}
gsl_vector_add_constant(v_temp, 1.0);
gsl_vector_div(Hi_eval, v_temp);
- gsl_vector_memcpy(HiHi_eval, Hi_eval);
+ gsl_vector_safe_memcpy(HiHi_eval, Hi_eval);
gsl_vector_mul(HiHi_eval, Hi_eval);
- gsl_vector_memcpy(HiHiHi_eval, HiHi_eval);
+ gsl_vector_safe_memcpy(HiHiHi_eval, HiHi_eval);
gsl_vector_mul(HiHiHi_eval, Hi_eval);
gsl_vector_set_all(v_temp, 1.0);
@@ -606,13 +606,13 @@ void LogL_dev12(double l, void *params, double *dev1, double *dev2) {
0.5 * (double)ni_test * (2.0 * yPKPKPy * P_yy - yPKPy * yPKPy) /
(P_yy * P_yy);
- gsl_matrix_free(Pab);
- gsl_matrix_free(PPab);
- gsl_matrix_free(PPPab);
- gsl_vector_free(Hi_eval);
- gsl_vector_free(HiHi_eval);
- gsl_vector_free(HiHiHi_eval);
- gsl_vector_free(v_temp);
+ gsl_matrix_safe_free(Pab); // FIXME: may contain NaN
+ gsl_matrix_safe_free(PPab); // FIXME: may contain NaN
+ gsl_matrix_safe_free(PPPab); // FIXME: may contain NaN
+ gsl_vector_safe_free(Hi_eval);
+ gsl_vector_safe_free(HiHi_eval);
+ gsl_vector_safe_free(HiHiHi_eval);
+ gsl_vector_safe_free(v_temp);
return;
}
@@ -641,12 +641,12 @@ double LogRL_f(double l, void *params) {
gsl_vector *Hi_eval = gsl_vector_safe_alloc((p->eval)->size);
gsl_vector *v_temp = gsl_vector_safe_alloc((p->eval)->size);
- gsl_vector_memcpy(v_temp, p->eval);
+ gsl_vector_safe_memcpy(v_temp, p->eval);
gsl_vector_scale(v_temp, l);
if (p->e_mode == 0) {
gsl_vector_set_all(Hi_eval, 1.0);
} else {
- gsl_vector_memcpy(Hi_eval, v_temp);
+ gsl_vector_safe_memcpy(Hi_eval, v_temp);
}
gsl_vector_add_constant(v_temp, 1.0);
gsl_vector_div(Hi_eval, v_temp);
@@ -675,10 +675,10 @@ double LogRL_f(double l, void *params) {
double c = 0.5 * df * (log(df) - log(2 * M_PI) - 1.0);
f = c - 0.5 * logdet_h - 0.5 * logdet_hiw - 0.5 * df * log(P_yy);
- gsl_matrix_free(Pab);
- gsl_matrix_free(Iab);
- gsl_vector_free(Hi_eval);
- gsl_vector_free(v_temp);
+ gsl_matrix_safe_free(Pab);
+ gsl_matrix_safe_free(Iab);
+ gsl_vector_safe_free(Hi_eval);
+ gsl_vector_safe_free(v_temp);
return f;
}
@@ -707,17 +707,17 @@ 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_memcpy(v_temp, p->eval);
+ gsl_vector_safe_memcpy(v_temp, p->eval);
gsl_vector_scale(v_temp, l);
if (p->e_mode == 0) {
gsl_vector_set_all(Hi_eval, 1.0);
} else {
- gsl_vector_memcpy(Hi_eval, v_temp);
+ gsl_vector_safe_memcpy(Hi_eval, v_temp);
}
gsl_vector_add_constant(v_temp, 1.0);
gsl_vector_div(Hi_eval, v_temp);
- gsl_vector_memcpy(HiHi_eval, Hi_eval);
+ gsl_vector_safe_memcpy(HiHi_eval, Hi_eval);
gsl_vector_mul(HiHi_eval, Hi_eval);
gsl_vector_set_all(v_temp, 1.0);
@@ -749,11 +749,11 @@ double LogRL_dev1(double l, void *params) {
dev1 = -0.5 * trace_PK + 0.5 * df * yPKPy / P_yy;
- gsl_matrix_free(Pab);
- gsl_matrix_free(PPab);
- gsl_vector_free(Hi_eval);
- gsl_vector_free(HiHi_eval);
- gsl_vector_free(v_temp);
+ gsl_matrix_safe_free(Pab); // FIXME: may contain NaN
+ gsl_matrix_safe_free(PPab); // FIXME: may contain NaN
+ gsl_vector_safe_free(Hi_eval);
+ gsl_vector_safe_free(HiHi_eval);
+ gsl_vector_safe_free(v_temp);
return dev1;
}
@@ -785,19 +785,19 @@ double LogRL_dev2(double l, void *params) {
gsl_vector *HiHiHi_eval = gsl_vector_safe_alloc((p->eval)->size);
gsl_vector *v_temp = gsl_vector_safe_alloc((p->eval)->size);
- gsl_vector_memcpy(v_temp, p->eval);
+ gsl_vector_safe_memcpy(v_temp, p->eval);
gsl_vector_scale(v_temp, l);
if (p->e_mode == 0) {
gsl_vector_set_all(Hi_eval, 1.0);
} else {
- gsl_vector_memcpy(Hi_eval, v_temp);
+ gsl_vector_safe_memcpy(Hi_eval, v_temp);
}
gsl_vector_add_constant(v_temp, 1.0);
gsl_vector_div(Hi_eval, v_temp);
- gsl_vector_memcpy(HiHi_eval, Hi_eval);
+ gsl_vector_safe_memcpy(HiHi_eval, Hi_eval);
gsl_vector_mul(HiHi_eval, Hi_eval);
- gsl_vector_memcpy(HiHiHi_eval, HiHi_eval);
+ gsl_vector_safe_memcpy(HiHiHi_eval, HiHi_eval);
gsl_vector_mul(HiHiHi_eval, Hi_eval);
gsl_vector_set_all(v_temp, 1.0);
@@ -837,13 +837,13 @@ double LogRL_dev2(double l, void *params) {
dev2 = 0.5 * trace_PKPK -
0.5 * df * (2.0 * yPKPKPy * P_yy - yPKPy * yPKPy) / (P_yy * P_yy);
- gsl_matrix_free(Pab);
- gsl_matrix_free(PPab);
- gsl_matrix_free(PPPab);
- gsl_vector_free(Hi_eval);
- gsl_vector_free(HiHi_eval);
- gsl_vector_free(HiHiHi_eval);
- gsl_vector_free(v_temp);
+ gsl_matrix_safe_free(Pab); // FIXME
+ gsl_matrix_safe_free(PPab);
+ gsl_matrix_safe_free(PPPab);
+ gsl_vector_safe_free(Hi_eval);
+ gsl_vector_safe_free(HiHi_eval);
+ gsl_vector_safe_free(HiHiHi_eval);
+ gsl_vector_safe_free(v_temp);
return dev2;
}
@@ -875,19 +875,19 @@ void LogRL_dev12(double l, void *params, double *dev1, double *dev2) {
gsl_vector *HiHiHi_eval = gsl_vector_safe_alloc((p->eval)->size);
gsl_vector *v_temp = gsl_vector_safe_alloc((p->eval)->size);
- gsl_vector_memcpy(v_temp, p->eval);
+ gsl_vector_safe_memcpy(v_temp, p->eval);
gsl_vector_scale(v_temp, l);
if (p->e_mode == 0) {
gsl_vector_set_all(Hi_eval, 1.0);
} else {
- gsl_vector_memcpy(Hi_eval, v_temp);
+ gsl_vector_safe_memcpy(Hi_eval, v_temp);
}
gsl_vector_add_constant(v_temp, 1.0);
gsl_vector_div(Hi_eval, v_temp);
- gsl_vector_memcpy(HiHi_eval, Hi_eval);
+ gsl_vector_safe_memcpy(HiHi_eval, Hi_eval);
gsl_vector_mul(HiHi_eval, Hi_eval);
- gsl_vector_memcpy(HiHiHi_eval, HiHi_eval);
+ gsl_vector_safe_memcpy(HiHiHi_eval, HiHi_eval);
gsl_vector_mul(HiHiHi_eval, Hi_eval);
gsl_vector_set_all(v_temp, 1.0);
@@ -929,13 +929,13 @@ void LogRL_dev12(double l, void *params, double *dev1, double *dev2) {
*dev2 = 0.5 * trace_PKPK -
0.5 * df * (2.0 * yPKPKPy * P_yy - yPKPy * yPKPy) / (P_yy * P_yy);
- gsl_matrix_free(Pab);
- gsl_matrix_free(PPab);
- gsl_matrix_free(PPPab);
- gsl_vector_free(Hi_eval);
- gsl_vector_free(HiHi_eval);
- gsl_vector_free(HiHiHi_eval);
- gsl_vector_free(v_temp);
+ gsl_matrix_safe_free(Pab); // FIXME
+ gsl_matrix_safe_free(PPab);
+ gsl_matrix_safe_free(PPPab);
+ gsl_vector_safe_free(Hi_eval);
+ gsl_vector_safe_free(HiHi_eval);
+ gsl_vector_safe_free(HiHiHi_eval);
+ gsl_vector_safe_free(v_temp);
return;
}
@@ -951,12 +951,12 @@ void LMM::CalcRLWald(const double &l, const FUNC_PARAM &params, double &beta,
gsl_vector *Hi_eval = gsl_vector_safe_alloc(params.eval->size);
gsl_vector *v_temp = gsl_vector_safe_alloc(params.eval->size);
- gsl_vector_memcpy(v_temp, params.eval);
+ gsl_vector_safe_memcpy(v_temp, params.eval);
gsl_vector_scale(v_temp, l);
if (params.e_mode == 0) {
gsl_vector_set_all(Hi_eval, 1.0);
} else {
- gsl_vector_memcpy(Hi_eval, v_temp);
+ gsl_vector_safe_memcpy(Hi_eval, v_temp);
}
gsl_vector_add_constant(v_temp, 1.0);
gsl_vector_div(Hi_eval, v_temp);
@@ -976,9 +976,9 @@ void LMM::CalcRLWald(const double &l, const FUNC_PARAM &params, double &beta,
se = sqrt(1.0 / (tau * P_xx));
p_wald = gsl_cdf_fdist_Q((P_yy - Px_yy) * tau, 1.0, df);
- gsl_matrix_free(Pab);
- gsl_vector_free(Hi_eval);
- gsl_vector_free(v_temp);
+ gsl_matrix_safe_free(Pab);
+ gsl_vector_safe_free(Hi_eval);
+ gsl_vector_safe_free(v_temp);
return;
}
@@ -993,12 +993,12 @@ void LMM::CalcRLScore(const double &l, const FUNC_PARAM &params, double &beta,
gsl_vector *Hi_eval = gsl_vector_safe_alloc(params.eval->size);
gsl_vector *v_temp = gsl_vector_safe_alloc(params.eval->size);
- gsl_vector_memcpy(v_temp, params.eval);
+ gsl_vector_safe_memcpy(v_temp, params.eval);
gsl_vector_scale(v_temp, l);
if (params.e_mode == 0) {
gsl_vector_set_all(Hi_eval, 1.0);
} else {
- gsl_vector_memcpy(Hi_eval, v_temp);
+ gsl_vector_safe_memcpy(Hi_eval, v_temp);
}
gsl_vector_add_constant(v_temp, 1.0);
gsl_vector_div(Hi_eval, v_temp);
@@ -1020,9 +1020,9 @@ void LMM::CalcRLScore(const double &l, const FUNC_PARAM &params, double &beta,
p_score =
gsl_cdf_fdist_Q((double)ni_test * P_xy * P_xy / (P_yy * P_xx), 1.0, df);
- gsl_matrix_free(Pab);
- gsl_vector_free(Hi_eval);
- gsl_vector_free(v_temp);
+ gsl_matrix_safe_free(Pab);
+ gsl_vector_safe_free(Hi_eval);
+ gsl_vector_safe_free(v_temp);
return;
}
@@ -1038,10 +1038,10 @@ void CalcUab(const gsl_matrix *UtW, const gsl_vector *Uty, gsl_matrix *Uab) {
}
if (a == n_cvt + 2) {
- gsl_vector_memcpy(u_a, Uty);
+ gsl_vector_safe_memcpy(u_a, Uty);
} else {
gsl_vector_const_view UtW_col = gsl_matrix_const_column(UtW, a - 1);
- gsl_vector_memcpy(u_a, &UtW_col.vector);
+ gsl_vector_safe_memcpy(u_a, &UtW_col.vector);
}
for (size_t b = a; b >= 1; --b) {
@@ -1053,17 +1053,17 @@ void CalcUab(const gsl_matrix *UtW, const gsl_vector *Uty, gsl_matrix *Uab) {
gsl_vector_view Uab_col = gsl_matrix_column(Uab, index_ab);
if (b == n_cvt + 2) {
- gsl_vector_memcpy(&Uab_col.vector, Uty);
+ gsl_vector_safe_memcpy(&Uab_col.vector, Uty);
} else {
gsl_vector_const_view UtW_col = gsl_matrix_const_column(UtW, b - 1);
- gsl_vector_memcpy(&Uab_col.vector, &UtW_col.vector);
+ gsl_vector_safe_memcpy(&Uab_col.vector, &UtW_col.vector);
}
gsl_vector_mul(&Uab_col.vector, u_a);
}
}
- gsl_vector_free(u_a);
+ gsl_vector_safe_free(u_a);
return;
}
@@ -1077,12 +1077,12 @@ void CalcUab(const gsl_matrix *UtW, const gsl_vector *Uty,
gsl_vector_view Uab_col = gsl_matrix_column(Uab, index_ab);
if (b == n_cvt + 2) {
- gsl_vector_memcpy(&Uab_col.vector, Uty);
+ gsl_vector_safe_memcpy(&Uab_col.vector, Uty);
} else if (b == n_cvt + 1) {
- gsl_vector_memcpy(&Uab_col.vector, Utx);
+ gsl_vector_safe_memcpy(&Uab_col.vector, Utx);
} else {
gsl_vector_const_view UtW_col = gsl_matrix_const_column(UtW, b - 1);
- gsl_vector_memcpy(&Uab_col.vector, &UtW_col.vector);
+ gsl_vector_safe_memcpy(&Uab_col.vector, &UtW_col.vector);
}
gsl_vector_mul(&Uab_col.vector, Utx);
@@ -1105,10 +1105,10 @@ void Calcab(const gsl_matrix *W, const gsl_vector *y, gsl_vector *ab) {
}
if (a == n_cvt + 2) {
- gsl_vector_memcpy(v_a, y);
+ gsl_vector_safe_memcpy(v_a, y);
} else {
gsl_vector_const_view W_col = gsl_matrix_const_column(W, a - 1);
- gsl_vector_memcpy(v_a, &W_col.vector);
+ gsl_vector_safe_memcpy(v_a, &W_col.vector);
}
for (size_t b = a; b >= 1; --b) {
@@ -1119,10 +1119,10 @@ void Calcab(const gsl_matrix *W, const gsl_vector *y, gsl_vector *ab) {
index_ab = GetabIndex(a, b, n_cvt);
if (b == n_cvt + 2) {
- gsl_vector_memcpy(v_b, y);
+ gsl_vector_safe_memcpy(v_b, y);
} else {
gsl_vector_const_view W_col = gsl_matrix_const_column(W, b - 1);
- gsl_vector_memcpy(v_b, &W_col.vector);
+ gsl_vector_safe_memcpy(v_b, &W_col.vector);
}
gsl_blas_ddot(v_a, v_b, &d);
@@ -1130,8 +1130,8 @@ void Calcab(const gsl_matrix *W, const gsl_vector *y, gsl_vector *ab) {
}
}
- gsl_vector_free(v_a);
- gsl_vector_free(v_b);
+ gsl_vector_safe_free(v_a);
+ gsl_vector_safe_free(v_b);
return;
}
@@ -1147,19 +1147,19 @@ void Calcab(const gsl_matrix *W, const gsl_vector *y, const gsl_vector *x,
index_ab = GetabIndex(n_cvt + 1, b, n_cvt);
if (b == n_cvt + 2) {
- gsl_vector_memcpy(v_b, y);
+ gsl_vector_safe_memcpy(v_b, y);
} else if (b == n_cvt + 1) {
- gsl_vector_memcpy(v_b, x);
+ gsl_vector_safe_memcpy(v_b, x);
} else {
gsl_vector_const_view W_col = gsl_matrix_const_column(W, b - 1);
- gsl_vector_memcpy(v_b, &W_col.vector);
+ gsl_vector_safe_memcpy(v_b, &W_col.vector);
}
gsl_blas_ddot(x, v_b, &d);
gsl_vector_set(ab, index_ab, d);
}
- gsl_vector_free(v_b);
+ gsl_vector_safe_free(v_b);
return;
}
@@ -1260,10 +1260,10 @@ void LMM::AnalyzeGene(const gsl_matrix *U, const gsl_vector *eval,
}
cout << endl;
- gsl_vector_free(y);
- gsl_vector_free(Uty);
- gsl_matrix_free(Uab);
- gsl_vector_free(ab);
+ gsl_vector_safe_free(y);
+ gsl_vector_safe_free(Uty);
+ gsl_matrix_safe_free(Uab);
+ gsl_vector_safe_free(ab);
infile.close();
infile.clear();
@@ -1324,7 +1324,7 @@ void LMM::Analyze(std::function< SnpNameValues(size_t) >& fetch_snp,
for (size_t i = 0; i < l; i++) {
// for every batch...
gsl_vector_view UtXlarge_col = gsl_matrix_column(UtXlarge, i);
- gsl_vector_memcpy(Utx, &UtXlarge_col.vector);
+ gsl_vector_safe_memcpy(Utx, &UtXlarge_col.vector);
CalcUab(UtW, Uty, Utx, Uab);
@@ -1425,7 +1425,7 @@ void LMM::Analyze(std::function< SnpNameValues(size_t) >& fetch_snp,
// copy genotype values for SNP into Xlarge cache
gsl_vector_view Xlarge_col = gsl_matrix_column(Xlarge, c % msize);
- gsl_vector_memcpy(&Xlarge_col.vector, x);
+ gsl_vector_safe_memcpy(&Xlarge_col.vector, x);
c++; // count SNPs going in
if (c % msize == 0)
@@ -1436,14 +1436,14 @@ void LMM::Analyze(std::function< SnpNameValues(size_t) >& fetch_snp,
// cout << "Counted SNPs " << c << " sumStat " << sumStat.size() << endl;
cout << endl;
- gsl_vector_free(x);
- gsl_vector_free(x_miss);
- gsl_vector_free(Utx);
- gsl_matrix_free(Uab);
- gsl_vector_free(ab);
+ gsl_vector_safe_free(x);
+ gsl_vector_safe_free(x_miss);
+ gsl_vector_safe_free(Utx);
+ gsl_matrix_safe_free(Uab);
+ gsl_vector_safe_free(ab);
- gsl_matrix_free(Xlarge);
- gsl_matrix_free(UtXlarge);
+ gsl_matrix_safe_free(Xlarge);
+ gsl_matrix_safe_free(UtXlarge);
}
@@ -1598,10 +1598,10 @@ void MatrixCalcLR(const gsl_matrix *U, const gsl_matrix *UtX,
pos_loglr.push_back(make_pair(i, log_lr));
}
- gsl_vector_free(w);
- gsl_matrix_free(Utw);
- gsl_matrix_free(Uab);
- gsl_vector_free(ab);
+ gsl_vector_safe_free(w);
+ gsl_matrix_safe_free(Utw);
+ gsl_matrix_safe_free(Uab);
+ gsl_vector_safe_free(ab);
return;
}
@@ -1790,8 +1790,8 @@ void CalcLambda(const char func_name, const gsl_vector *eval,
CalcLambda(func_name, param0, l_min, l_max, n_region, lambda, logl_H0);
- gsl_matrix_free(Uab);
- gsl_vector_free(ab);
+ gsl_matrix_safe_free(Uab);
+ gsl_vector_safe_free(ab);
return;
}
@@ -1816,8 +1816,8 @@ void CalcPve(const gsl_vector *eval, const gsl_matrix *UtW,
pve = trace_G * lambda / (trace_G * lambda + 1.0);
pve_se = trace_G / ((trace_G * lambda + 1.0) * (trace_G * lambda + 1.0)) * se;
- gsl_matrix_free(Uab);
- gsl_vector_free(ab);
+ gsl_matrix_safe_free(Uab);
+ gsl_vector_safe_free(ab);
return;
}
@@ -1843,14 +1843,14 @@ void CalcLmmVgVeBeta(const gsl_vector *eval, const gsl_matrix *UtW,
gsl_matrix_set_zero(Uab);
CalcUab(UtW, Uty, Uab);
- gsl_vector_memcpy(v_temp, eval);
+ gsl_vector_safe_memcpy(v_temp, eval);
gsl_vector_scale(v_temp, lambda);
gsl_vector_set_all(Hi_eval, 1.0);
gsl_vector_add_constant(v_temp, 1.0);
gsl_vector_div(Hi_eval, v_temp);
// Calculate beta.
- gsl_matrix_memcpy(HiW, UtW);
+ gsl_matrix_safe_memcpy(HiW, UtW);
for (size_t i = 0; i < UtW->size2; i++) {
gsl_vector_view HiW_col = gsl_matrix_column(HiW, i);
gsl_vector_mul(&HiW_col.vector, Hi_eval);
@@ -1881,15 +1881,15 @@ void CalcLmmVgVeBeta(const gsl_vector *eval, const gsl_matrix *UtW,
gsl_vector_set(se_beta, i, sqrt(gsl_matrix_get(Vbeta, i, i)));
}
- gsl_matrix_free(Uab);
- gsl_matrix_free(Pab);
- gsl_vector_free(ab);
- gsl_vector_free(Hi_eval);
- gsl_vector_free(v_temp);
- gsl_matrix_free(HiW);
- gsl_matrix_free(WHiW);
- gsl_vector_free(WHiy);
- gsl_matrix_free(Vbeta);
+ gsl_matrix_safe_free(Uab);
+ gsl_matrix_safe_free(Pab);
+ gsl_vector_safe_free(ab);
+ gsl_vector_safe_free(Hi_eval);
+ gsl_vector_safe_free(v_temp);
+ gsl_matrix_safe_free(HiW);
+ gsl_matrix_safe_free(WHiW);
+ gsl_vector_safe_free(WHiy);
+ gsl_matrix_safe_free(Vbeta);
gsl_permutation_free(pmt);
return;
@@ -1929,7 +1929,7 @@ void LMM::AnalyzeBimbamGXE(const gsl_matrix *U, const gsl_vector *eval,
gsl_matrix *UtW_expand = gsl_matrix_safe_alloc(U->size1, UtW->size2 + 2);
gsl_matrix_view UtW_expand_mat =
gsl_matrix_submatrix(UtW_expand, 0, 0, U->size1, UtW->size2);
- gsl_matrix_memcpy(&UtW_expand_mat.matrix, UtW);
+ gsl_matrix_safe_memcpy(&UtW_expand_mat.matrix, UtW);
gsl_vector_view UtW_expand_env = gsl_matrix_column(UtW_expand, UtW->size2);
gsl_blas_dgemv(CblasTrans, 1.0, U, env, 0.0, &UtW_expand_env.vector);
gsl_vector_view UtW_expand_x = gsl_matrix_column(UtW_expand, UtW->size2 + 1);
@@ -2030,13 +2030,13 @@ void LMM::AnalyzeBimbamGXE(const gsl_matrix *U, const gsl_vector *eval,
}
cout << endl;
- gsl_vector_free(x);
- gsl_vector_free(x_miss);
- gsl_vector_free(Utx);
- gsl_matrix_free(Uab);
- gsl_vector_free(ab);
+ gsl_vector_safe_free(x);
+ gsl_vector_safe_free(x_miss);
+ gsl_vector_safe_free(Utx);
+ gsl_matrix_safe_free(Uab);
+ gsl_vector_safe_free(ab);
- gsl_matrix_free(UtW_expand);
+ gsl_matrix_safe_free(UtW_expand);
infile.close();
infile.clear();
@@ -2078,7 +2078,7 @@ void LMM::AnalyzePlinkGXE(const gsl_matrix *U, const gsl_vector *eval,
gsl_matrix *UtW_expand = gsl_matrix_safe_alloc(U->size1, UtW->size2 + 2);
gsl_matrix_view UtW_expand_mat =
gsl_matrix_submatrix(UtW_expand, 0, 0, U->size1, UtW->size2);
- gsl_matrix_memcpy(&UtW_expand_mat.matrix, UtW);
+ gsl_matrix_safe_memcpy(&UtW_expand_mat.matrix, UtW);
gsl_vector_view UtW_expand_env = gsl_matrix_column(UtW_expand, UtW->size2);
gsl_blas_dgemv(CblasTrans, 1.0, U, env, 0.0, &UtW_expand_env.vector);
gsl_vector_view UtW_expand_x = gsl_matrix_column(UtW_expand, UtW->size2 + 1);
@@ -2208,12 +2208,12 @@ void LMM::AnalyzePlinkGXE(const gsl_matrix *U, const gsl_vector *eval,
}
cout << endl;
- gsl_vector_free(x);
- gsl_vector_free(Utx);
- gsl_matrix_free(Uab);
- gsl_vector_free(ab);
+ gsl_vector_safe_free(x);
+ gsl_vector_safe_free(Utx);
+ gsl_matrix_safe_free(Uab);
+ gsl_vector_safe_free(ab);
- gsl_matrix_free(UtW_expand);
+ gsl_matrix_safe_free(UtW_expand);
infile.close();
infile.clear();
diff --git a/src/mathfunc.cpp b/src/mathfunc.cpp
index ba71b64..908a52d 100644
--- a/src/mathfunc.cpp
+++ b/src/mathfunc.cpp
@@ -64,6 +64,33 @@ bool has_nan(const vector<double> v) {
return false;
}
+bool has_nan(const gsl_vector *v) {
+ for (size_t i = 0; i < v->size; ++i)
+ if (gsl_isnan(gsl_vector_get(v,i))) return true;
+ return false;
+}
+bool has_inf(const gsl_vector *v) {
+ for (size_t i = 0; i < v->size; ++i) {
+ auto value = gsl_vector_get(v,i);
+ if (gsl_isinf(value) != 0) return true;
+ }
+ return false;
+}
+bool has_nan(const gsl_matrix *m) {
+ for (size_t i = 0; i < m->size1; ++i)
+ for (size_t j = 0; j < m->size2; ++j)
+ if (gsl_isnan(gsl_matrix_get(m,i,j))) return true;
+ return false;
+}
+bool has_inf(const gsl_matrix *m) {
+ for (size_t i = 0; i < m->size1; ++i)
+ for (size_t j = 0; j < m->size2; ++j) {
+ auto value = gsl_matrix_get(m,i,j);
+ if (gsl_isinf(value) != 0) return true;
+ }
+ return false;
+}
+
// calculate variance of a vector
double VectorVar(const gsl_vector *v) {
double d, m = 0.0, m2 = 0.0;
@@ -96,8 +123,8 @@ void CenterMatrix(gsl_matrix *G) {
}
}
- gsl_vector_free(w);
- gsl_vector_free(Gw);
+ gsl_vector_safe_free(w);
+ gsl_vector_safe_free(Gw);
return;
}
@@ -120,7 +147,7 @@ void CenterMatrix(gsl_matrix *G, const gsl_vector *w) {
}
}
- gsl_vector_free(Gw);
+ gsl_vector_safe_free(Gw);
return;
}
@@ -156,12 +183,12 @@ void CenterMatrix(gsl_matrix *G, const gsl_matrix *W) {
gsl_matrix_add(G, Gtmp);
- gsl_matrix_free(WtW);
- gsl_matrix_free(WtWi);
- gsl_matrix_free(WtWiWt);
- gsl_matrix_free(GW);
- gsl_matrix_free(WtGW);
- gsl_matrix_free(Gtmp);
+ gsl_matrix_safe_free(WtW);
+ gsl_matrix_safe_free(WtWi);
+ gsl_matrix_safe_free(WtWiWt);
+ gsl_matrix_safe_free(GW);
+ gsl_matrix_safe_free(WtGW);
+ gsl_matrix_safe_free(Gtmp);
return;
}
@@ -228,7 +255,7 @@ bool isMatrixSymmetric(const gsl_matrix *G) {
bool isMatrixPositiveDefinite(const gsl_matrix *G) {
enforce(G->size1 == G->size2);
auto G2 = gsl_matrix_safe_alloc(G->size1, G->size2);
- enforce_gsl(gsl_matrix_memcpy(G2,G));
+ enforce_gsl(gsl_matrix_safe_memcpy(G2,G));
auto handler = gsl_set_error_handler_off();
#if GSL_MAJOR_VERSION >= 2 && GSL_MINOR_VERSION >= 3
auto s = gsl_linalg_cholesky_decomp1(G2);
@@ -236,20 +263,20 @@ bool isMatrixPositiveDefinite(const gsl_matrix *G) {
auto s = gsl_linalg_cholesky_decomp(G2);
#endif
gsl_set_error_handler(handler);
- gsl_matrix_free(G2);
+ gsl_matrix_safe_free(G2);
return (s == GSL_SUCCESS);
}
gsl_vector *getEigenValues(const gsl_matrix *G) {
enforce(G->size1 == G->size2);
auto G2 = gsl_matrix_safe_alloc(G->size1, G->size2);
- enforce_gsl(gsl_matrix_memcpy(G2,G));
+ enforce_gsl(gsl_matrix_safe_memcpy(G2,G));
auto eworkspace = gsl_eigen_symm_alloc(G->size1);
enforce(eworkspace);
gsl_vector *eigenvalues = gsl_vector_safe_alloc(G->size1);
enforce_gsl(gsl_eigen_symm(G2, eigenvalues, eworkspace));
gsl_eigen_symm_free(eworkspace);
- gsl_matrix_free(G2);
+ gsl_matrix_safe_free(G2);
return eigenvalues;
}
@@ -359,9 +386,9 @@ void CenterVector(gsl_vector *y, const gsl_matrix *W) {
gsl_blas_dgemv(CblasNoTrans, -1.0, W, WtWiWty, 1.0, y);
- gsl_matrix_free(WtW);
- gsl_vector_free(Wty);
- gsl_vector_free(WtWiWty);
+ gsl_matrix_safe_free(WtW);
+ gsl_vector_safe_free(Wty);
+ gsl_vector_safe_free(WtWiWty);
return;
}
@@ -388,9 +415,9 @@ void StandardizeVector(gsl_vector *y) {
// Calculate UtX.
void CalcUtX(const gsl_matrix *U, gsl_matrix *UtX) {
gsl_matrix *X = gsl_matrix_safe_alloc(UtX->size1, UtX->size2);
- gsl_matrix_memcpy(X, UtX);
+ gsl_matrix_safe_memcpy(X, UtX);
fast_dgemm("T", "N", 1.0, U, X, 0.0, UtX);
- gsl_matrix_free(X);
+ gsl_matrix_safe_free(X);
}
void CalcUtX(const gsl_matrix *U, const gsl_matrix *X, gsl_matrix *UtX) {
@@ -407,7 +434,7 @@ void Kronecker(const gsl_matrix *K, const gsl_matrix *V, gsl_matrix *H) {
for (size_t j = 0; j < K->size2; j++) {
gsl_matrix_view H_sub = gsl_matrix_submatrix(
H, i * V->size1, j * V->size2, V->size1, V->size2);
- gsl_matrix_memcpy(&H_sub.matrix, V);
+ gsl_matrix_safe_memcpy(&H_sub.matrix, V);
gsl_matrix_scale(&H_sub.matrix, gsl_matrix_get(K, i, j));
}
}
@@ -420,13 +447,13 @@ void KroneckerSym(const gsl_matrix *K, const gsl_matrix *V, gsl_matrix *H) {
for (size_t j = i; j < K->size2; j++) {
gsl_matrix_view H_sub = gsl_matrix_submatrix(
H, i * V->size1, j * V->size2, V->size1, V->size2);
- gsl_matrix_memcpy(&H_sub.matrix, V);
+ gsl_matrix_safe_memcpy(&H_sub.matrix, V);
gsl_matrix_scale(&H_sub.matrix, gsl_matrix_get(K, i, j));
if (i != j) {
gsl_matrix_view H_sub_sym = gsl_matrix_submatrix(
H, j * V->size1, i * V->size2, V->size1, V->size2);
- gsl_matrix_memcpy(&H_sub_sym.matrix, &H_sub.matrix);
+ gsl_matrix_safe_memcpy(&H_sub_sym.matrix, &H_sub.matrix);
}
}
}
diff --git a/src/mathfunc.h b/src/mathfunc.h
index 1319a64..42e76e5 100644
--- a/src/mathfunc.h
+++ b/src/mathfunc.h
@@ -30,6 +30,10 @@ using namespace std;
bool has_nan(const vector<double> v);
+bool has_nan(const gsl_vector *v);
+bool has_inf(const gsl_vector *v);
+bool has_nan(const gsl_matrix *m);
+bool has_inf(const gsl_matrix *m);
double VectorVar(const gsl_vector *v);
void CenterMatrix(gsl_matrix *G);
generated by cgit v1.2.3 (git 2.44.0) at 2025-01-18 19:21:44 +0000