Allow frequency mixture for binary model (thanks Edward Braun)

Cibiv · Jun 5, 2020 · 192cad7 · 192cad7
1 parent b55b280
commit 192cad7
Show file tree

Hide file tree

Showing 3 changed files with 54 additions and 13 deletions.
diff --git a/model/modelbin.cpp b/model/modelbin.cpp
@@ -32,16 +32,47 @@ void ModelBIN::init(const char *model_name, string model_params, StateFreqType f
 	name = model_name;
 	full_name = model_name;
 	if (name == "JC2") {
-		freq = FREQ_EQUAL;
+		def_freq = FREQ_EQUAL;
 	} else if (name == "GTR2") {
-		freq = FREQ_ESTIMATE;
+		def_freq = FREQ_ESTIMATE;
 	} else {
 		readParameters(model_name);
 	}
+    if (freq_params != "") {
+        readStateFreq(freq_params);
+    }
+    if (model_params != "") {
+      readRates(model_params);
+    }
 	if (freq == FREQ_UNKNOWN || def_freq == FREQ_EQUAL) freq = def_freq;
 	ModelMarkov::init(freq);
 }
 
 void ModelBIN::startCheckpoint() {
     checkpoint->startStruct("ModelBIN");
 }
+
+string ModelBIN::getNameParams() {
+    //if (num_params == 0) return name;
+    ostringstream retname;
+    retname << name;
+//    if (!fixed_parameters) {
+//        retname << '{';
+//        int nrates = getNumRateEntries();
+//        for (int i = 0; i < nrates; i++) {
+//            if (i>0) retname << ',';
+//            retname << rates[i];
+//        }
+//        retname << '}';
+//    }
+//    getNameParamsFreq(retname);
+    retname << freqTypeString(freq_type, phylo_tree->aln->seq_type, true);
+    if (freq_type == FREQ_EMPIRICAL || freq_type == FREQ_ESTIMATE ||
+        (freq_type == FREQ_USER_DEFINED)) {
+        retname << "{" << state_freq[0];
+        for (int i = 1; i < num_states; i++)
+            retname << "," << state_freq[i];
+        retname << "}";
+    }
+    return retname.str();
+}
diff --git a/model/modelbin.h b/model/modelbin.h
@@ -48,7 +48,7 @@ class ModelBIN : public ModelMarkov
 	/**
 	 * @return model name with parameters in form of e.g. GTR{a,b,c,d,e,f}
 	 */
-	virtual string getNameParams() { return name; }
+    virtual string getNameParams();
 
     /**
         start structure for checkpointing

diff --git a/model/modelmarkov.cpp b/model/modelmarkov.cpp
@@ -354,6 +354,8 @@ void ModelMarkov::writeInfo(ostream &out) {
         report_rates(out, "Rate parameters", rates);
         report_state_freqs(out);
 		//if (freq_type != FREQ_ESTIMATE) return;
+    } else if (is_reversible && num_states == 2) {
+        report_state_freqs(out);
 	} else if (!is_reversible) {
         // non-reversible
 //        int i;
@@ -401,16 +403,24 @@ void ModelMarkov::report_rates(ostream& out, string title, double *r) {
 }
 
 void ModelMarkov::report_state_freqs(ostream& out, double *custom_state_freq) {
-  double *f;
-  if (custom_state_freq) f = custom_state_freq;
-  else f = state_freq;
-  out << setprecision(3);
-  out << "Base frequencies:";
-  out << "  A: " << f[0];
-  out << "  C: " << f[1];
-  out << "  G: " << f[2];
-  out << "  T: " << f[3];
-  out << endl;
+    double *f;
+    if (custom_state_freq) f = custom_state_freq;
+    else f = state_freq;
+    if (num_states == 4) {
+        out << setprecision(3);
+        out << "Base frequencies:";
+        out << "  A: " << f[0];
+        out << "  C: " << f[1];
+        out << "  G: " << f[2];
+        out << "  T: " << f[3];
+        out << endl;
+    } else if (num_states == 2) {
+        out << setprecision(3);
+        out << "State frequencies:";
+        out << "  0: " << f[0];
+        out << "  1: " << f[1];
+        out << endl;
+    }
 }
 
 void ModelMarkov::computeTransMatrixNonrev(double time, double *trans_matrix, int mixture) {