properly handle zero state frequencies via rate matrix reduction trick

amkozlov · Mar 4, 2019 · 30c6100 · 30c6100
1 parent 6ec2c7d
commit 30c6100
Show file tree

Hide file tree

Showing 5 changed files with 27 additions and 12 deletions.
diff --git a/libs/pll-modules b/libs/pll-modules
diff --git a/src/Model.cpp b/src/Model.cpp
@@ -328,7 +328,7 @@ void Model::set_user_freqs(doubleVector& freqs)
   bool invalid = false;
   for (auto v: freqs)
   {
-    invalid |= (v <= 0. || v >= 1.);
+    invalid |= (v < 0. || v >= 1.);
   }
 
   if (invalid)

diff --git a/src/log.hpp b/src/log.hpp
@@ -113,6 +113,8 @@ Logging& logger();
 #define FMT_MOD(p) setprecision(logger().precision(LogElement::model)) << p
 #define FMT_BL(bl) setprecision(logger().precision(LogElement::brlen)) << bl
 #define FMT_PREC3(val) setprecision(3) << val
+#define FMT_PREC6(val) setprecision(6) << val
+#define FMT_PREC9(val) setprecision(9) << val
 
 template <class T>
 LogStream& operator<<(LogStream& logstream, const T& object)

diff --git a/src/main.cpp b/src/main.cpp
@@ -474,6 +474,7 @@ size_t total_free_params(const RaxmlInstance& instance)
 void check_models(const RaxmlInstance& instance)
 {
   const auto& opts = instance.opts;
+  bool zero_freqs = false;
 
   for (const auto& pinfo: instance.parted_msa->part_list())
   {
@@ -502,17 +503,24 @@ void check_models(const RaxmlInstance& instance)
         const auto& freqs = stats.emp_base_freqs;
         for (unsigned int i = 0; i < freqs.size(); ++i)
         {
-          if (!(freqs[i] > 0.))
+          if (freqs[i] < PLL_EIGEN_MINFREQ)
           {
-            LOG_ERROR << "\nBase frequencies: ";
+            if (!zero_freqs)
+            {
+              LOG_WARN << endl;
+              zero_freqs = true;
+            }
+
+            LOG_WARN << "WARNING: State " << to_string(i) <<
+                (instance.parted_msa->part_count() > 1 ? " in partition " + pinfo.name() : "") <<
+                " has very low frequency (" << FMT_PREC9(freqs[i]) << ")!" << endl;
+
+            LOG_VERB << "Base frequencies: ";
             for (unsigned int j = 0; j < freqs.size(); ++j)
-              LOG_ERROR << freqs[j] <<  " ";
-            LOG_ERROR << endl;
+              LOG_VERB << FMT_PREC9(freqs[j]) <<  " ";
+            LOG_VERB << endl;
+
 
-            throw runtime_error("Frequency of state " + to_string(i) +
-                                " in partition " + pinfo.name() + " is 0!\n"
-                                "Please either change your partitioning scheme or "
-                                "use model state frequencies for this partition!");
           }
         }
       }
@@ -532,7 +540,7 @@ void check_models(const RaxmlInstance& instance)
         {
           LOG_ERROR << "\nBase frequencies: ";
           for (unsigned int j = 0; j < freqs.size(); ++j)
-            LOG_ERROR << freqs[j] <<  " ";
+            LOG_ERROR << FMT_PREC9(freqs[j]) <<  " ";
           LOG_ERROR << endl;
 
           throw runtime_error("User-specified stationary base frequencies"
@@ -564,6 +572,11 @@ void check_models(const RaxmlInstance& instance)
     }
   }
 
+  if (zero_freqs)
+  {
+    LOG_WARN << endl << "WARNING: Some states have very low frequencies, "
+        "which might lead to numerical issues!" << endl;
+  }
 
   /* Check for extreme cases of overfitting (K >= n) */
   if (opts.safety_checks.isset(SafetyCheck::model_overfit))

diff --git a/src/version.h b/src/version.h
@@ -1,2 +1,2 @@
 #define RAXML_VERSION "0.8.0git BETA"
-#define RAXML_DATE "28.02.2019"
+#define RAXML_DATE "04.03.2019"