Merge pull request #73 from StefanFlaumberg/bugfix

Allow for single-state alignments and remove misleading warnings

Author: bqminh

alignment/alignment.cpp

@@ -123,33 +123,37 @@ double chi2prob (int deg, double chi2)
 } /* chi2prob */
 
 
-int Alignment::checkAbsentStates(string msg) {
-    double *state_freq = new double[num_states];
-    computeStateFreq(state_freq);
-    string absent_states, rare_states;
-    int count = 0;
-    // Skip check for PoMo.
+void Alignment::checkAbsentStates(string msg) {
+    // skip checking for PoMo
     if (seq_type == SEQ_POMO)
-      return 0;
-    for (int i = 0; i < num_states; i++)
-        if (state_freq[i] == 0.0) {
+        return;
+    string absent_states, rare_states;
+    int absent_cnt = 0;
+    double *state_freqs = new double[num_states];
+    computeStateFreq(state_freqs);
+    for (int x = 0; x < num_states; ++x) {
+        if (state_freqs[x] == 0.0) {
             if (!absent_states.empty())
                 absent_states += ", ";
-            absent_states += convertStateBackStr(i);
-            count++;
-        } else if (state_freq[i] <= Params::getInstance().min_state_freq) {
+            absent_states += convertStateBackStr(x);
+            absent_cnt++;
+        } else if (state_freqs[x] <= Params::getInstance().min_state_freq) {
             if (!rare_states.empty())
                 rare_states += ", ";
-            rare_states += convertStateBackStr(i);
+            rare_states += convertStateBackStr(x);
         }
-    if (count >= num_states-1 && Params::getInstance().fixed_branch_length != BRLEN_FIX)
-        outError("Only one state is observed in " + msg);
+    }
+    delete [] state_freqs;
+    if (absent_cnt == num_states)
+        outError("Only gaps observed in " + msg);
+    if (absent_cnt == num_states - 1)
+        outWarning("Only one state observed in " + msg);
+    if (absent_cnt > 0)
+        outWarning(convertIntToString(absent_cnt) + " states (see below) not observed in " + msg);
     if (!absent_states.empty())
-        cout << "NOTE: State(s) " << absent_states << " not present in " << msg << " and thus removed from Markov process to prevent numerical problems" << endl;
+        outWarning("State(s) " + absent_states + " not present in " + msg + " and may cause numerical problems");
     if (!rare_states.empty())
-        cout << "WARNING: States(s) " << rare_states << " rarely appear in " << msg << " and may cause numerical problems" << endl;
-    delete[] state_freq;
-    return count;
+        outWarning("State(s) " + rare_states + " rarely appear in " + msg + " and may cause numerical problems");
 }
 
 void Alignment::checkSeqName() {

alignment/alignment.h

@@ -501,11 +501,10 @@ class Alignment : public vector<Pattern>, public CharSet, public StateSpace {
     int getMaxSeqNameLength();
 
     /*
-        check if some state is absent, which may cause numerical issues
+        check if some states are absent, which may cause numerical issues
         @param msg additional message into the warning
-        @return number of absent states in the alignment
     */
-    virtual int checkAbsentStates(string msg);
+    virtual void checkAbsentStates(string msg);
 
     /**
             check proper and undupplicated sequence names

alignment/superalignment.cpp

@@ -1148,11 +1148,9 @@ Alignment *SuperAlignment::removeIdenticalSeq(string not_remove, bool keep_two,
 	return aln;
 }
 
-int SuperAlignment::checkAbsentStates(string msg) {
-    int count = 0;
-    for (auto it = partitions.begin(); it != partitions.end(); ++it)
-        count += (*it)->checkAbsentStates("partition " + convertIntToString((it-partitions.begin())+1));
-    return count;
+void SuperAlignment::checkAbsentStates(string msg) {
+    for (vector<Alignment*>::iterator it = partitions.begin(); it != partitions.end(); ++it)
+        (*it)->checkAbsentStates("partition " + convertIntToString((it-partitions.begin())+1));
 }
 
 /*

alignment/superalignment.h

@@ -171,11 +171,10 @@ class SuperAlignment : public Alignment
 
 
     /*
-        check if some state is absent, which may cause numerical issues
+        check if some states are absent, which may cause numerical issues
         @param msg additional message into the warning
-        @return number of absent states in the alignment
     */
-    virtual int checkAbsentStates(string msg);
+    virtual void checkAbsentStates(string msg);
 
 	/**
 		Quit if some sequences contain only gaps or missing data

main/phyloanalysis.cpp

@@ -3420,22 +3420,8 @@ bool isTreeMixture(Params& params) {
 
 void runTreeReconstruction(Params &params, IQTree* &iqtree) {
 
-    //    string dist_file;
-    // params.startCPUTime = getCPUTime();
-    // params.start_real_time = getRealTime();
-    
-    int absent_states = 0;
-    if (iqtree->isSuperTree()) {
-        PhyloSuperTree *stree = (PhyloSuperTree*)iqtree;
-        for (auto i = stree->begin(); i != stree->end(); i++)
-            absent_states += (*i)->aln->checkAbsentStates("partition " + (*i)->aln->name);
-    } else {
-        absent_states = iqtree->aln->checkAbsentStates("alignment");
-    }
-    if (absent_states > 0) {
-        cout << "NOTE: " << absent_states << " states (see above) are not present and thus removed from Markov process to prevent numerical problems" << endl;
-    }
-    
+    iqtree->aln->checkAbsentStates("alignment");
+
     // Make sure that no partial likelihood of IQ-TREE is initialized when PLL is used to save memory
     if (params.pll) {
         iqtree->deleteAllPartialLh();