From 08a6b07cf7eb95de94423c84b613b1c4258d2825 Mon Sep 17 00:00:00 2001
From: Minh Bui <m.bui@anu.edu.au>
Date: Fri, 15 Jun 2018 15:26:39 +1000
Subject: [PATCH] Fix bug initFromCatMinusOne() creating negative rates in edge
 cases (reported by Paul Frandsen)

---
 main/phylotesting.cpp |  4 ++++
 model/ratefree.cpp    | 14 ++++++++++----
 2 files changed, 14 insertions(+), 4 deletions(-)

diff --git a/main/phylotesting.cpp b/main/phylotesting.cpp
index e994b268..e80b6a20 100644
--- a/main/phylotesting.cpp
+++ b/main/phylotesting.cpp
@@ -973,6 +973,7 @@ int getModelList(Params &params, Alignment *aln, StrVector &models, bool separat
 	bool test_options_noASC_I_new[] = {true,false,  false, true,  false,    false, true,    false};
 	bool test_options_asc_new[]   ={false,false,   true,false,  false,     true,false,     true};
 	bool test_options_pomo[]      = {true, false,  false, true,   false,   false,false,    false};
+    bool test_options_norate[]    = {true, false,  false, false,   false,   false,false,    false};
     bool *test_options = test_options_default;
 //	bool test_options_codon[] =  {true,false,  false,false,  false,    false};
 	const int noptions = sizeof(rate_options) / sizeof(char*);
@@ -1163,6 +1164,9 @@ int getModelList(Params &params, Alignment *aln, StrVector &models, bool separat
             test_options = test_options_morph;
         else
             test_options = test_options_noASC_I;
+    } else if (aln->frac_invariant_sites >= 1.0) {
+        // 2018-06-12: alignment with only invariant sites, no rate variation added
+        test_options = test_options_norate;
 	} else {
         // normal data, use +I instead
         if (with_new) {
diff --git a/model/ratefree.cpp b/model/ratefree.cpp
index eca95c01..8a57e0a5 100644
--- a/model/ratefree.cpp
+++ b/model/ratefree.cpp
@@ -128,13 +128,13 @@ void RateFree::initFromCatMinusOne() {
     restoreCheckpoint();
     ncategory++;
 
-    int first = 0, second = -1, i;
+    int first = 0, i;
     // get the category k with largest proportion
     for (i = 1; i < ncategory-1; i++)
         if (prop[i] > prop[first]) {
             first = i;
         }
-    second = (first == 0) ? 1 : 0;
+    int second = (first == 0) ? 1 : 0;
     for (i = 0; i < ncategory-1; i++)
         if (prop[i] > prop[second] && second != first)
             second = i;
@@ -143,9 +143,15 @@ void RateFree::initFromCatMinusOne() {
 //    memmove(prop, input->prop, (k+1)*sizeof(double));
 
     // divide highest category into 2 of the same prop
-    rates[ncategory-1] = (-rates[second] + 3*rates[first])/2.0;
+    // 2018-06-12: fix bug negative rates
+    if (-rates[second] + 3*rates[first] > 0.0) {
+        rates[ncategory-1] = (-rates[second] + 3*rates[first])/2.0;
+        rates[first] = (rates[second]+rates[first])/2.0;
+    } else {
+        rates[ncategory-1] = (3*rates[first])/2.0;
+        rates[first] = (rates[first])/2.0;
+    }
     prop[ncategory-1] = prop[first]/2;
-    rates[first] = (rates[second]+rates[first])/2.0;
     prop[first] = prop[first]/2;
 //    if (k < ncategory-2) {
 //        memcpy(&rates[k+2], &input->rates[k+1], (ncategory-2-k)*sizeof(double));