aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/debug.h6
-rw-r--r--src/gemma.cpp37
-rw-r--r--src/main.cpp10
-rw-r--r--src/mathfunc.cpp8
-rw-r--r--src/mathfunc.h1
-rw-r--r--src/param.h1
6 files changed, 60 insertions, 3 deletions
diff --git a/src/debug.h b/src/debug.h
index fafebcf..3fbe9e0 100644
--- a/src/debug.h
+++ b/src/debug.h
@@ -7,10 +7,15 @@
// enforce works like assert but also when NDEBUG is set (i.e., it
// always works). enforce_msg prints message instead of expr
+#define ROUND(f) round(f * 10000.)/10000
#if defined NDEBUG
#define debug_msg(msg)
+#define assert_issue(is_issue, expr)
#else
#define debug_msg(msg) cout << "**** DEBUG: " << msg << endl;
+#define assert_issue(is_issue, expr) \
+ ((is_issue) ? enforce_msg(expr,"FAIL: ISSUE assert") : __ASSERT_VOID_CAST(0))
+
#endif
/* This prints an "Assertion failed" message and aborts. */
@@ -49,4 +54,5 @@ inline void __enforce_fail(const char *__assertion, const char *__file,
: __enforce_fail(gsl_strerror(COMBINE(res, __LINE__)), __FILE__, \
__LINE__, __ASSERT_FUNCTION))
+
#endif
diff --git a/src/gemma.cpp b/src/gemma.cpp
index a678b8c..3197043 100644
--- a/src/gemma.cpp
+++ b/src/gemma.cpp
@@ -697,6 +697,7 @@ void GEMMA::PrintHelp(size_t option) {
cout << " -silence silent terminal display" << endl;
cout << " -debug debug output" << endl;
cout << " -nind [num] read up to num individuals" << endl;
+ cout << " -issue [num] enable tests relevant to issue tracker" << endl;
cout << endl;
}
@@ -1347,6 +1348,16 @@ void GEMMA::Assign(int argc, char **argv, PARAM &cPar) {
str.clear();
str.assign(argv[i]);
cPar.ni_max = atoi(str.c_str()); // for testing purposes
+ enforce(cPar.ni_max > 0);
+ } else if (strcmp(argv[i], "-issue") == 0) {
+ if (argv[i + 1] == NULL || argv[i + 1][0] == '-') {
+ continue;
+ }
+ ++i;
+ str.clear();
+ str.assign(argv[i]);
+ cPar.issue = atoi(str.c_str()); // for testing purposes
+ enforce(cPar.issue > 0);
} else if (strcmp(argv[i], "-emp") == 0) {
if (argv[i + 1] == NULL || argv[i + 1][0] == '-') {
continue;
@@ -2656,11 +2667,12 @@ return;}
gsl_matrix *se_B = gsl_matrix_alloc(Y->size2, W->size2);
gsl_matrix *G = gsl_matrix_alloc(Y->size1, Y->size1);
gsl_matrix *U = gsl_matrix_alloc(Y->size1, Y->size1);
- gsl_matrix *UtW = gsl_matrix_alloc(Y->size1, W->size2);
- gsl_matrix *UtY = gsl_matrix_alloc(Y->size1, Y->size2);
- gsl_vector *eval = gsl_vector_alloc(Y->size1);
+ gsl_matrix *UtW = gsl_matrix_calloc(Y->size1, W->size2);
+ gsl_matrix *UtY = gsl_matrix_calloc(Y->size1, Y->size2);
+ gsl_vector *eval = gsl_vector_calloc(Y->size1);
gsl_vector *env = gsl_vector_alloc(Y->size1);
gsl_vector *weight = gsl_vector_alloc(Y->size1);
+ assert_issue(cPar.issue == 26, UtY->data[0] == 0.0);
// set covariates matrix W and phenotype matrix Y
// an intercept should be included in W,
@@ -2769,6 +2781,8 @@ return;}
CalcUtX(U, W, UtW);
CalcUtX(U, Y, UtY);
+ assert_issue(cPar.issue == 26, ROUND(UtY->data[0]) == -16.6143);
+
LMM cLmm;
cLmm.CopyFromParam(cPar);
@@ -2784,6 +2798,7 @@ return;}
// calculate UtW and Uty
CalcUtX(U, W, UtW);
CalcUtX(U, Y, UtY);
+ assert_issue(cPar.issue == 26, ROUND(UtY->data[0]) == -16.6143);
// calculate REMLE/MLE estimate and pve for univariate model
if (cPar.n_ph == 1) { // one phenotype
@@ -2791,24 +2806,40 @@ return;}
gsl_vector_view se_beta = gsl_matrix_row(se_B, 0);
gsl_vector_view UtY_col = gsl_matrix_column(UtY, 0);
+ assert_issue(cPar.issue == 26, ROUND(UtY->data[0]) == -16.6143);
+
CalcLambda('L', eval, UtW, &UtY_col.vector, cPar.l_min, cPar.l_max,
cPar.n_region, cPar.l_mle_null, cPar.logl_mle_H0);
+ assert(!std::isnan(UtY->data[0]));
+ assert(!std::isnan(B->data[0]));
+ assert(!std::isnan(se_B->data[0]));
+
CalcLmmVgVeBeta(eval, UtW, &UtY_col.vector, cPar.l_mle_null,
cPar.vg_mle_null, cPar.ve_mle_null, &beta.vector,
&se_beta.vector);
+ assert(!std::isnan(UtY->data[0]));
+ assert(!std::isnan(B->data[0]));
+ assert(!std::isnan(se_B->data[0]));
+
cPar.beta_mle_null.clear();
cPar.se_beta_mle_null.clear();
for (size_t i = 0; i < B->size2; 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));
}
+ assert(!std::isnan(UtY->data[0]));
+ assert(!std::isnan(B->data[0]));
+ assert(!std::isnan(se_B->data[0]));
+ assert(!std::isnan(cPar.beta_mle_null.front()));
+ assert(!std::isnan(cPar.se_beta_mle_null.front()));
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; i < B->size2; i++) {
diff --git a/src/main.cpp b/src/main.cpp
index 833136c..e37b154 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -25,10 +25,20 @@
using namespace std;
+void gemma_gsl_error_handler (const char * reason,
+ const char * file,
+ int line, int gsl_errno) {
+ cerr << "GSL ERROR: " << reason << " in " << file
+ << " at line " << line << " errno " << gsl_errno <<endl;
+ exit(22);
+}
+
int main(int argc, char *argv[]) {
GEMMA cGemma;
PARAM cPar;
+ gsl_set_error_handler (&gemma_gsl_error_handler);
+
if (argc <= 1) {
cGemma.PrintHeader();
return EXIT_SUCCESS;
diff --git a/src/mathfunc.cpp b/src/mathfunc.cpp
index 9e19bf1..9a4bb8b 100644
--- a/src/mathfunc.cpp
+++ b/src/mathfunc.cpp
@@ -187,6 +187,14 @@ double ScaleMatrix(gsl_matrix *G) {
return d;
}
+double SumVector(const gsl_vector *v) {
+ double sum = 0;
+ for (int i = 0; i < v->size; i++ ) {
+ sum += gsl_vector_get(v, i);
+ }
+ return( sum );
+}
+
// Center the vector y.
double CenterVector(gsl_vector *y) {
double d = 0.0;
diff --git a/src/mathfunc.h b/src/mathfunc.h
index 29eafe4..b9f3c73 100644
--- a/src/mathfunc.h
+++ b/src/mathfunc.h
@@ -32,6 +32,7 @@ void CenterMatrix(gsl_matrix *G, const gsl_vector *w);
void CenterMatrix(gsl_matrix *G, const gsl_matrix *W);
void StandardizeMatrix(gsl_matrix *G);
double ScaleMatrix(gsl_matrix *G);
+double SumVector(const gsl_vector *v);
double CenterVector(gsl_vector *y);
void CenterVector(gsl_vector *y, const gsl_matrix *W);
void StandardizeVector(gsl_vector *y);
diff --git a/src/param.h b/src/param.h
index 45d8c0f..f4d649f 100644
--- a/src/param.h
+++ b/src/param.h
@@ -114,6 +114,7 @@ public:
// IO-related parameters.
bool mode_silence;
bool mode_debug = false;
+ uint issue; // enable tests for issue on github tracker
int a_mode; // Analysis mode, 1/2/3/4 for Frequentist tests
int k_mode; // Kinship read mode: 1: n by n matrix, 2: id/id/k_value;
vector<size_t> p_column; // Which phenotype column needs analysis.
generated by cgit v1.2.3 (git 2.46.2) at 2025-01-31 19:41:50 +0000