Add extant sequence reconstruction (-esr) and gap-aware ASR/ESR (-gap-asr, -gap-esr)

alignment/alignment.cpp

@@ -179,6 +179,27 @@ void Alignment::checkAbsentStates(string msg) {
     }
 }
 
+bool Alignment::containSingleStateOnly(const int& state)
+{
+    // the alignment contains only one single state if:
+    // 1. it contains only one pattern
+    // 2. that pattern contains only one state
+    if (size() == 1
+        && at(0).num_chars == 1)
+    {
+        // if a specific state is specified (not -1), check if that state is the only one present in the alignment
+        if (state >= 0)
+        {
+            ASSERT(state < num_states);
+            return at(0).freqs[state] == at(0).size();
+        }
+        // otherwise, don't need to check that constraint
+        return true;
+    }
+
+    return false;
+}
+
 void Alignment::checkSeqName() {
     ostringstream warn_str;
     StrVector::iterator it;
@@ -3986,6 +4007,54 @@ Alignment *Alignment::convertCodonToDNA() const {
     return aln;
 }
 
+void Alignment::convertToBin(Alignment* res, const string& new_model_name)
+{
+    // metadata copy, matching the convention used by initAlignmentCopy()
+    // (see convertToCodonOrAA/convertCodonToAA/convertCodonToDNA)
+    res->seq_names = seq_names;
+    res->name = name;
+    res->position_spec = position_spec;
+    res->model_name = new_model_name; //res->model_name = model_name;
+    res->aln_file = aln_file;
+    res->sequence_type = "BIN"; //res->sequence_type = sequence_type;
+    res->char_partition = char_partition;
+    res->tree_len = tree_len;
+    res->seq_type = SEQ_BINARY;
+    res->num_states = 2;
+
+    res->computeUnknownState();
+
+    res->clear();
+    res->pattern_index.clear();
+
+    VerboseMode save_mode = verbose_mode;
+    verbose_mode = min(verbose_mode, VB_MIN); // to avoid printing gappy sites in addPattern
+    size_t nsite = getNSite();
+    size_t nseq = getNSeq();
+    Pattern pat;
+    pat.resize(nseq);
+
+    for (size_t site = 0; site < nsite; ++site)
+    {
+        for (size_t seq = 0; seq < nseq; ++seq)
+        {
+            StateType state = at(getPatternID(site))[seq];
+            pat[seq] = state == STATE_UNKNOWN ? 0 : 1;
+        }
+        res->addPattern(pat);
+    }
+    verbose_mode = save_mode;
+    res->updateConstPatterns();
+    res->countConstSites();
+}
+
+Alignment* Alignment::convertToBin(const string& new_model_name)
+{
+    Alignment *res = new Alignment;
+    convertToBin(res, new_model_name);
+    return res;
+}
+
 static const char *convert_site_range(const char *str, int &lower, int &upper,
                                       int &step) {
     char *endptr;

alignment/alignment.h

@@ -503,6 +503,13 @@ class Alignment : public vector<Pattern>, public CharSet, public StateSpace {
     */
     virtual void checkAbsentStates(string msg);
 
+    /*
+        check if the alignment contains only one single state
+        @param[in] state the state to check, if not specified, there is no constraint on the only single state
+        @return TRUE if the alignment contains only one single state; otherwise return FALSE
+    */
+    bool containSingleStateOnly(const int& state = -1);
+
     /**
             check proper and undupplicated sequence names
      */
@@ -695,6 +702,20 @@ class Alignment : public vector<Pattern>, public CharSet, public StateSpace {
      */
     Alignment *convertCodonToDNA() const;
 
+    /**
+        convert an alignment into binary (gap/non-gap) alignment
+        @param[in] model_name name of model for the new alignment
+        @return a pointer to a new alignment
+    */
+    virtual Alignment* convertToBin(const string& model_name = "GTR2");
+
+    /**
+        convert an alignment into binary (gap/non-gap) alignment
+        @param[in] model_name name of model for the new alignment
+        @param[out] output_aln a pointer to a new alignment
+    */
+    void convertToBin(Alignment* output_aln, const string& model_name);
+
     /**
      @param quartet ID of four taxa
      @param[out] support number of sites supporting 12|34, 13|24 and 14|23

alignment/superalignment.cpp

@@ -1891,3 +1891,28 @@ void SuperAlignment::orderPatternByNumChars(int pat_type) {
     // TODO compute pars_lower_bound (lower bound of pars score for remaining patterns)
     delete [] pars_lower_bound;
 }
+
+SuperAlignment* SuperAlignment::convertToBin(const string& new_model_name)
+{
+    SuperAlignment* new_aln = new SuperAlignment;
+    // convert the base alignment
+    Alignment::convertToBin(new_aln, new_model_name);
+
+    // convert alignment members one by one
+    for (vector<Alignment*>::iterator it = partitions.begin(); it != partitions.end(); it++) {
+        new_aln->partitions.push_back((*it)->convertToBin(new_model_name));
+    }
+
+    // clone SuperAlignment-specific variables
+    new_aln->max_num_states = 2;
+    new_aln->taxa_index = taxa_index;
+    // discard the binary gap/non-gap patterns written by Alignment::convertToBin() above:
+    // SuperAlignment::buildPattern() (called via init() below) builds its own
+    // partition-presence pattern instead, and asserts the pattern list starts empty
+    new_aln->clear();
+    new_aln->pattern_index.clear();
+    new_aln->site_pattern.clear();
+    new_aln->init();
+
+    return new_aln;
+}

alignment/superalignment.h

@@ -296,6 +296,13 @@ class SuperAlignment : public Alignment
     /** order pattern by number of character states and return in ptn_order
     */
     virtual void orderPatternByNumChars(int pat_type);
+
+    /**
+        convert an alignment into binary (gap/non-gap) alignment
+        @param[in] model_name name of model for the new alignment
+        @return a pointer to a new alignment
+    */
+    virtual SuperAlignment* convertToBin(const string& model_name = "GTR2");
 
 	/**
 		actual partition alignments

main/main.cpp

@@ -2259,43 +2259,9 @@ int main(int argc, char *argv[]) {
     parseArg(argc, argv, Params::getInstance());
 
     // 2015-12-05
-    Checkpoint *checkpoint = new Checkpoint;
     string filename = (string)Params::getInstance().out_prefix +".ckp.gz";
-    checkpoint->setFileName(filename);
-
     bool append_log = false;
-
-    if (!Params::getInstance().ignore_checkpoint && fileExists(filename)) {
-        checkpoint->load();
-        if (checkpoint->hasKey("finished")) {
-            if (checkpoint->getBool("finished")) {
-                if (Params::getInstance().force_unfinished) {
-                    if (MPIHelper::getInstance().isMaster())
-                        cout << "NOTE: Continue analysis although a previous run already finished" << endl;
-                } else {
-                    delete checkpoint;
-                    if (MPIHelper::getInstance().isMaster())
-                        outError("Checkpoint (" + filename + ") indicates that a previous run successfully finished\n" +
-                            "Use `-redo` option if you really want to redo the analysis and overwrite all output files.\n" +
-                            "Use `--redo-tree` option if you want to restore ModelFinder and only redo tree search.\n" +
-                            "Use `--undo` option if you want to continue previous run when changing/adding options."
-                        );
-                    else
-                        exit(EXIT_SUCCESS);
-                    exit(EXIT_FAILURE);
-                } 
-            } else {
-                append_log = true;
-            }
-        } else {
-            if (MPIHelper::getInstance().isMaster())
-                outWarning("Ignore invalid checkpoint file " + filename);
-            checkpoint->clear();
-        }
-    }
-
-    if (MPIHelper::getInstance().isWorker())
-        checkpoint->setFileName("");
+    Checkpoint *checkpoint = initAndCheckCheckpoint(filename, Params::getInstance(), append_log);
 
     _log_file = Params::getInstance().out_prefix;
     _log_file += ".log";