Matrices positive definite issues

4 files changed, 13 insertions, 10 deletions

diff --git a/src/debug.cpp b/src/debug.cpp
index 529d603..3a62d2a 100644
--- a/src/debug.cpp
+++ b/src/debug.cpp
@@ -332,11 +332,11 @@ void do_validate_K(const gsl_matrix *K, const char *__pretty_function, const cha
       warning_at_msg(__file,__line,"K is ill conditioned!");
     if (!isMatrixSymmetric(K))
       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) {
+    const bool negative_eigen_values = has_negative_values_but_one(eigenvalues);
+    if (negative_eigen_values) {
       warning_at_msg(__file,__line,"K has more than one negative eigenvalues!");
     }
-    if (count_small>1 && negative_values && !isMatrixPositiveDefinite(K))
+    if (count_small>1 && negative_eigen_values && !isMatrixPositiveDefinite(K))
       warnfail_at_msg(is_strict_mode(),__pretty_function,__file,__line,"K is not positive definite!");
     gsl_vector_free(eigenvalues);
   }
diff --git a/src/gemma.cpp b/src/gemma.cpp
index c3caa96..4e3ae88 100644
--- a/src/gemma.cpp
+++ b/src/gemma.cpp
@@ -1612,6 +1612,7 @@ void GEMMA::Assign(int argc, char **argv, PARAM &cPar) {
     } else if (strcmp(argv[i], "-strict") == 0) {
       // cPar.mode_strict = true;
       debug_set_strict_mode(true);
+      debug_set_debug_mode(true);
     } else if (strcmp(argv[i], "-legacy") == 0) {
       debug_set_legacy_mode(true);
       warning_msg("you are running in legacy mode - support may drop in future versions of gemma");
@@ -3567,5 +3568,7 @@ void GEMMA::WriteLog(int argc, char **argv, PARAM &cPar) {
 
   outfile.close();
   outfile.clear();
+
+  info_msg("Done");
   return;
 }
diff --git a/src/lapack.cpp b/src/lapack.cpp
index a66705b..bf73338 100644
--- a/src/lapack.cpp
+++ b/src/lapack.cpp
@@ -267,7 +267,7 @@ double EigenDecomp_Zeroed(gsl_matrix *G, gsl_matrix *U, gsl_vector *eval,
     // checks
     if (gsl_vector_get(eval,i) == 0.0)
       count_zero_eigenvalues += 1;
-    if (gsl_vector_get(eval,i) < 0.0) // count smaller than -EIGEN_MINVALUE
+    if (gsl_vector_get(eval,i) < -EIGEN_MINVALUE) // count smaller than -EIGEN_MINVALUE
       count_negative_eigenvalues += 1;
     d += gsl_vector_get(eval, i);
   }
@@ -279,11 +279,11 @@ double EigenDecomp_Zeroed(gsl_matrix *G, gsl_matrix *U, gsl_vector *eval,
     msg += " eigenvalues close to zero";
     warning_msg(msg);
   }
-  if (count_negative_eigenvalues > 0) {
+  const bool negative_eigen_values = has_negative_values_but_one(eval);
+  if (negative_eigen_values) {
     write(eval,"eigenvalues");
-    warning_msg("Matrix G has more than one negative eigenvalues!");
+    warning_msg("K has more than one negative eigenvalues!");
   }
-
   return d;
 }
 
diff --git a/src/mathfunc.cpp b/src/mathfunc.cpp
index 542093e..8acbd52 100644
--- a/src/mathfunc.cpp
+++ b/src/mathfunc.cpp
@@ -104,7 +104,7 @@ bool is_float(const std::string & s){
 }
 
 double safe_log(const double d) {
-  if (!is_legacy_mode() && !is_check_mode())
+  if (!is_legacy_mode() && (is_check_mode() || is_debug_mode()))
     enforce_msg(d > 0.0, (std::string("Trying to take the log of ") + std::to_string(d)).c_str());
   return log(d);
 }
@@ -113,7 +113,7 @@ double safe_sqrt(const double d) {
   double d1 = d;
   if (fabs(d < 0.001))
     d1 = fabs(d);
-  if (!is_legacy_mode() && !is_check_mode())
+  if (!is_legacy_mode() && (is_check_mode() || is_debug_mode()))
     enforce_msg(d1 >= 0.0, (std::string("Trying to take the sqrt of ") + std::to_string(d)).c_str());
   if (d1 < 0.0 )
     return nan("");
@@ -354,7 +354,7 @@ tuple<double, double, double> abs_minmax(const gsl_vector *v) {
 bool has_negative_values_but_one(const gsl_vector *v) {
   bool one_skipped = false;
   for (size_t i=0; i<v->size; i++) {
-    if (v->data[i] < 0.0) {
+    if (v->data[i] < -EIGEN_MINVALUE) {
       if (one_skipped)
         return true;
       one_skipped = true;