version 1.5.6

Cibiv · Dec 4, 2017 · 0b47ae8 · 0b47ae8
1 parent b4f0006
commit 0b47ae8
Show file tree

Hide file tree

Showing 2 changed files with 57 additions and 58 deletions.
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -50,7 +50,7 @@ add_definitions(-DIQ_TREE)
 # The version number.
 set (iqtree_VERSION_MAJOR 1)
 set (iqtree_VERSION_MINOR 5)
-set (iqtree_VERSION_PATCH "5")
+set (iqtree_VERSION_PATCH "6")
 
 set(BUILD_SHARED_LIBS OFF)
 

diff --git a/alignment.cpp b/alignment.cpp
@@ -119,6 +119,8 @@ int Alignment::checkAbsentStates(string msg) {
                 rare_states += ", ";
             rare_states += convertStateBackStr(i);
         }
+    if (absent_states.size() >= num_states-1)
+        outError("Only one state is 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;
     if (!rare_states.empty())
@@ -165,67 +167,64 @@ void Alignment::checkSeqName() {
     computeStateFreq(state_freq);
 //    computeStateFreqPerSequence(freq_per_sequence);
     countStatePerSequence(count_per_seq);
+
+    int i, df = -1;
+    for (i = 0; i < num_states; i++)
+        if (state_freq[i] > 0.0) df++;
 
-    /*if (verbose_mode >= VB_MIN)*/ {
-        int max_len = getMaxSeqNameLength()+1;
-//        cout << "  ID  ";
-//        cout <<  "  Sequence";
-        cout.width(max_len+14);
-        cout << right << "Gap/Ambiguity" << "  Composition  p-value"<< endl;
-        int num_problem_seq = 0;
-        int total_gaps = 0;
-        cout.precision(2);
-        int num_failed = 0;
-        for (int i = 0; i < seq_names.size(); i++) {
-            int j;
-            int num_gaps = getNSite() - countProperChar(i);
-            total_gaps += num_gaps;
-            double percent_gaps = ((double)num_gaps / getNSite())*100.0;
-			cout.width(4);
-			cout << right << i+1 << "  ";
-            cout.width(max_len);
-            cout << left << seq_names[i] << " ";
-			cout.width(6);
-//			cout << num_gaps << " (" << percent_gaps << "%)";
-            cout << right << percent_gaps << "%";
-            if (percent_gaps > 50) {
-//				cout << " !!!";
-				num_problem_seq++;
-			}
-//            cout << "\t" << seq_states[i].size();
+    int max_len = getMaxSeqNameLength()+1;
+    cout.width(max_len+14);
+    cout << right << "Gap/Ambiguity" << "  Composition  p-value"<< endl;
+    int num_problem_seq = 0;
+    int total_gaps = 0;
+    cout.precision(2);
+    int num_failed = 0;
+    for (i = 0; i < seq_names.size(); i++) {
+        int j;
+        int num_gaps = getNSite() - countProperChar(i);
+        total_gaps += num_gaps;
+        double percent_gaps = ((double)num_gaps / getNSite())*100.0;
+        cout.width(4);
+        cout << right << i+1 << "  ";
+        cout.width(max_len);
+        cout << left << seq_names[i] << " ";
+        cout.width(6);
+        cout << right << percent_gaps << "%";
+        if (percent_gaps > 50) {
+            num_problem_seq++;
+        }
 
-            double chi2 = 0.0;
-            unsigned sum_count = 0;
-            for (j = 0; j < num_states; j++)
-                sum_count += count_per_seq[i*num_states+j];
-            double sum_inv = 1.0/sum_count;
-            for (j = 0; j < num_states; j++)
-                freq_per_sequence[j] = count_per_seq[i*num_states+j]*sum_inv;
-            for (j = 0; j < num_states; j++)
+        double chi2 = 0.0;
+        unsigned sum_count = 0;
+        for (j = 0; j < num_states; j++)
+            sum_count += count_per_seq[i*num_states+j];
+        double sum_inv = 1.0/sum_count;
+        for (j = 0; j < num_states; j++)
+            freq_per_sequence[j] = count_per_seq[i*num_states+j]*sum_inv;
+        for (j = 0; j < num_states; j++)
+            if (state_freq[j] > 0.0)
                 chi2 += (state_freq[j] - freq_per_sequence[j]) * (state_freq[j] - freq_per_sequence[j]) / state_freq[j];
-
-//            chi2 *= getNSite();
-            chi2 *= sum_count;
-            double pvalue = chi2prob(num_states-1, chi2);
-            if (pvalue < 0.05) {
-                cout << "    failed ";
-                num_failed++;
-            } else
-                cout << "    passed ";
-            cout.width(9);
-            cout << right << pvalue*100 << "%";
-//            cout << "  " << chi2;
-			cout << endl;
-        }
-        if (num_problem_seq) cout << "WARNING: " << num_problem_seq << " sequences contain more than 50% gaps/ambiguity" << endl;
-        cout << "**** ";
-        cout.width(max_len+2);
-        cout << left << " TOTAL  ";
-        cout.width(6);
-        cout << right << ((double)total_gaps/getNSite())/getNSeq()*100 << "% ";
-        cout << " " << num_failed << " sequences failed composition chi2 test (p-value<5%; df=" << num_states-1 << ")" << endl;
-        cout.precision(3);
+
+        chi2 *= sum_count;
+        double pvalue = chi2prob(df, chi2);
+        if (pvalue < 0.05) {
+            cout << "    failed ";
+            num_failed++;
+        } else
+            cout << "    passed ";
+        cout.width(9);
+        cout << right << pvalue*100 << "%";
+        cout << endl;
     }
+    if (num_problem_seq) cout << "WARNING: " << num_problem_seq << " sequences contain more than 50% gaps/ambiguity" << endl;
+    cout << "**** ";
+    cout.width(max_len+2);
+    cout << left << " TOTAL  ";
+    cout.width(6);
+    cout << right << ((double)total_gaps/getNSite())/getNSeq()*100 << "% ";
+    cout << " " << num_failed << " sequences failed composition chi2 test (p-value<5%; df=" << df << ")" << endl;
+    cout.precision(3);
+
     delete [] count_per_seq;
     delete [] freq_per_sequence;
     delete [] state_freq;