Fix bug initFromCatMinusOne() creating negative rates in edge cases (…

…reported by Paul Frandsen)

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) {

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));