CMakeLists.txt

@@ -51,7 +51,7 @@ add_definitions(-DIQ_TREE)
 # The version number.
 set (iqtree_VERSION_MAJOR 1)
 set (iqtree_VERSION_MINOR 6)
-set (iqtree_VERSION_PATCH "3")
+set (iqtree_VERSION_PATCH "5")
 
 set(BUILD_SHARED_LIBS OFF)
 
@@ -71,6 +71,7 @@ endif()
 
 if (CMAKE_GENERATOR MATCHES "Xcode")
     set(CMAKE_XCODE_ATTRIBUTE_DEBUG_INFORMATION_FORMAT "dwarf-with-dsym")
+    set(CMAKE_XCODE_ATTRIBUTE_COMPILER_INDEX_STORE_ENABLE "No")
 endif()
 
 include_directories("${PROJECT_SOURCE_DIR}")
@@ -278,6 +279,7 @@ if (NOT IQTREE_FLAGS MATCHES "single")
   	elseif (CLANG)
 		set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -pthread")
   		set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fopenmp=libomp -pthread")
+        set (CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -fopenmp=libomp")
   	endif()
 else()
 	message("OpenMP        : NONE")

alignment/alignment.cpp

@@ -123,7 +123,7 @@ int Alignment::checkAbsentStates(string msg) {
                 rare_states += ", ";
             rare_states += convertStateBackStr(i);
         }
-    if (count >= num_states-1)
+    if (count >= num_states-1 && Params::getInstance().fixed_branch_length != BRLEN_FIX)
         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;
@@ -138,12 +138,7 @@ void Alignment::checkSeqName() {
     StrVector::iterator it;
     for (it = seq_names.begin(); it != seq_names.end(); it++) {
         string orig_name = (*it);
-        for (string::iterator i = it->begin(); i != it->end(); i++) {
-            if (!isalnum(*i) && (*i) != '_' && (*i) != '-' && (*i) != '.') {
-                (*i) = '_';
-            }
-        }
-        if (orig_name != (*it))
+        if (renameString(*it))
             warn_str << orig_name << " -> " << (*it) << endl;
     }
     if (warn_str.str() != "") {
@@ -394,7 +389,12 @@ void Alignment::checkGappySeq(bool force_error) {
     }
 }
 
-Alignment::Alignment(char *filename, char *sequence_type, InputType &intype) : vector<Pattern>() {
+Alignment::Alignment(char *filename, char *sequence_type, InputType &intype, string model) : vector<Pattern>() {
+    name = "Noname";
+    this->model_name = model;
+    if (sequence_type)
+        this->sequence_type = sequence_type;
+    aln_file = filename;
     num_states = 0;
     frac_const_sites = 0.0;
     frac_invariant_sites = 0.0;
@@ -2636,6 +2636,11 @@ void Alignment::extractSubAlignment(Alignment *aln, IntVector &seq_id, int min_t
         ASSERT(*it >= 0 && *it < aln->getNSeq());
         seq_names.push_back(aln->getSeqName(*it));
     }
+    name = aln->name;
+    model_name = aln->model_name;
+    sequence_type = aln->sequence_type;
+    position_spec = aln->position_spec;
+    aln_file = aln->aln_file;
     num_states = aln->num_states;
     seq_type = aln->seq_type;
     STATE_UNKNOWN = aln->STATE_UNKNOWN;
@@ -2684,6 +2689,11 @@ void Alignment::extractPatterns(Alignment *aln, IntVector &ptn_id) {
     for (i = 0; i < aln->getNSeq(); i++) {
         seq_names.push_back(aln->getSeqName(i));
     }
+    name = aln->name;
+    model_name = aln->model_name;
+    sequence_type = aln->sequence_type;
+    position_spec = aln->position_spec;
+    aln_file = aln->aln_file;
     num_states = aln->num_states;
     seq_type = aln->seq_type;
     STATE_UNKNOWN = aln->STATE_UNKNOWN;
@@ -2720,6 +2730,11 @@ void Alignment::extractPatternFreqs(Alignment *aln, IntVector &ptn_freq) {
     for (i = 0; i < aln->getNSeq(); i++) {
         seq_names.push_back(aln->getSeqName(i));
     }
+    name = aln->name;
+    model_name = aln->model_name;
+    sequence_type = aln->sequence_type;
+    position_spec = aln->position_spec;
+    aln_file = aln->aln_file;
     num_states = aln->num_states;
     seq_type = aln->seq_type;
     genetic_code = aln->genetic_code;
@@ -2756,6 +2771,11 @@ void Alignment::extractSites(Alignment *aln, IntVector &site_id) {
     for (i = 0; i < aln->getNSeq(); i++) {
         seq_names.push_back(aln->getSeqName(i));
     }
+    name = aln->name;
+    model_name = aln->model_name;
+    sequence_type = aln->sequence_type;
+    position_spec = aln->position_spec;
+    aln_file = aln->aln_file;
     num_states = aln->num_states;
     seq_type = aln->seq_type;
     STATE_UNKNOWN = aln->STATE_UNKNOWN;
@@ -2799,6 +2819,11 @@ void Alignment::convertToCodonOrAA(Alignment *aln, char *gene_code_id, bool nt2a
     for (i = 0; i < aln->getNSeq(); i++) {
         seq_names.push_back(aln->getSeqName(i));
     }
+    name = aln->name;
+    model_name = aln->model_name;
+    sequence_type = aln->sequence_type;
+    position_spec = aln->position_spec;
+    aln_file = aln->aln_file;
 //    num_states = aln->num_states;
     seq_type = SEQ_CODON;
     initCodon(gene_code_id);
@@ -2971,6 +2996,11 @@ void Alignment::extractSites(Alignment *aln, const char* spec) {
 
 void Alignment::createBootstrapAlignment(Alignment *aln, IntVector* pattern_freq, const char *spec) {
     if (aln->isSuperAlignment()) outError("Internal error: ", __func__);
+    name = aln->name;
+    model_name = aln->model_name;
+    sequence_type = aln->sequence_type;
+    position_spec = aln->position_spec;
+    aln_file = aln->aln_file;
     int site, nsite = aln->getNSite();
     seq_names.insert(seq_names.begin(), aln->seq_names.begin(), aln->seq_names.end());
     num_states = aln->num_states;
@@ -3235,6 +3265,11 @@ void Alignment::createBootstrapAlignment(int *pattern_freq, const char *spec, in
 void Alignment::buildFromPatternFreq(Alignment & aln, IntVector new_pattern_freqs){
 	int nsite = aln.getNSite();
     seq_names.insert(seq_names.begin(), aln.seq_names.begin(), aln.seq_names.end());
+    name = aln.name;
+    model_name = aln.model_name;
+    sequence_type = aln.sequence_type;
+    position_spec = aln.position_spec;
+    aln_file = aln.aln_file;
     num_states = aln.num_states;
     seq_type = aln.seq_type;
 
@@ -3274,6 +3309,11 @@ void Alignment::createGapMaskedAlignment(Alignment *masked_aln, Alignment *aln)
 
     int site, nsite = aln->getNSite(), nseq = aln->getNSeq();
     seq_names.insert(seq_names.begin(), aln->seq_names.begin(), aln->seq_names.end());
+    name = aln->name;
+    model_name = aln->model_name;
+    sequence_type = aln->sequence_type;
+    position_spec = aln->position_spec;
+    aln_file = aln->aln_file;
     num_states = aln->num_states;
     seq_type = aln->seq_type;
     genetic_code = aln->genetic_code;
@@ -3343,6 +3383,11 @@ void Alignment::concatenateAlignment(Alignment *aln) {
 void Alignment::copyAlignment(Alignment *aln) {
     int site, nsite = aln->getNSite();
     seq_names.insert(seq_names.begin(), aln->seq_names.begin(), aln->seq_names.end());
+    name = aln->name;
+    model_name = aln->model_name;
+    sequence_type = aln->sequence_type;
+    position_spec = aln->position_spec;
+    aln_file = aln->aln_file;
     num_states = aln->num_states;
     seq_type = aln->seq_type;
     genetic_code = aln->genetic_code;

alignment/alignment.h

@@ -65,7 +65,7 @@ Multiple Sequence Alignment. Stored by a vector of site-patterns
 
         @author BUI Quang Minh, Steffen Klaere, Arndt von Haeseler <minh.bui@univie.ac.at>
  */
-class Alignment : public vector<Pattern> {
+class Alignment : public vector<Pattern>, public CharSet {
     friend class SuperAlignment;
 
 public:
@@ -81,7 +81,7 @@ public:
             @param sequence_type type of the sequence, either "BIN", "DNA", "AA", or NULL
             @param intype (OUT) input format of the file
      */
-    Alignment(char *filename, char *sequence_type, InputType &intype);
+    Alignment(char *filename, char *sequence_type, InputType &intype, string model);
 
     /**
             destructor

alignment/maalignment.h

@@ -34,7 +34,7 @@ class MaAlignment : public Alignment
 public:
     MaAlignment() : Alignment() {};
 
-    MaAlignment(char *filename,  char *sequence_type, InputType &intype) : Alignment(filename, sequence_type, intype){};
+    MaAlignment(char *filename,  char *sequence_type, InputType &intype, string model) : Alignment(filename, sequence_type, intype, model){};
 	
 	MaAlignment(Alignment &align) : Alignment(align){};
 

alignment/superalignment.h

@@ -24,12 +24,6 @@
 
 
 struct PartitionInfo {
-	string name; // partition name
-	string model_name; // model name
-	string aln_file; // alignment file associated
-	string sequence_type; // sequence type (DNA/AA/BIN)
-	string position_spec; // position specification, e.g., "1-100\1 1-100\2"
-
 	double cur_score;	// current log-likelihood
 	double part_rate;	// partition heterogeneity rate
 	int    evalNNIs;	// number of evaluated NNIs on subtree
@@ -73,7 +67,7 @@ class SuperAlignment : public Alignment
 {
 public:
 	/** constructor initialize from a supertree */
-    SuperAlignment(PhyloSuperTree *super_tree);
+    SuperAlignment(Params &params);
 
 	/** constructor initialize empty alignment */
     SuperAlignment();
@@ -81,9 +75,27 @@ public:
     /** destructor */
     ~SuperAlignment();
 
+    void init(StrVector *sequence_names = NULL);
+    
     /** return that this is a super-alignment structure */
 	virtual bool isSuperAlignment() { return true; }
 
+    /** read partition model file */
+    void readPartition(Params &params);
+    
+    /** read RAxML-style partition file */
+    void readPartitionRaxml(Params &params);
+    
+    /** read partition model file in NEXUS format into variable info */
+    void readPartitionNexus(Params &params);
+    
+    void printPartition(const char *filename);
+    
+    void printPartitionRaxml(const char *filename);
+    
+    void printBestPartition(const char *filename);
+    void printBestPartitionRaxml(const char *filename);
+
 	/**
 	 * create taxa_index from super-alignment to sub-alignment
 	 * @param part index of sub-alignment
@@ -209,7 +221,7 @@ public:
 	 * print all sub alignments into files with prefix, suffix is the charset name
 	 * @param prefix prefix of output files
 	 */
-	void printSubAlignments(Params &params, vector<PartitionInfo> &part_info);
+	void printSubAlignments(Params &params);
 
 	/**
 		@return unconstrained log-likelihood (without a tree)

debian/changelog

+iqtree (1.6.5+dfsg-1) unstable; urgency=medium
+
+  * New upstream version
+  * Delete RPATH from /usr/bin/iqtree-mpi
+
+ -- Andreas Tille <tille@debian.org>  Tue, 22 May 2018 17:43:03 +0200
+
 iqtree (1.6.3+dfsg-1) unstable; urgency=medium
 
   * New upstream version

debian/control

@@ -11,7 +11,8 @@ Build-Depends: debhelper (>= 11~),
                mpi-default-dev,
                zlib1g-dev,
                help2man,
-               time
+               time,
+               chrpath
 Standards-Version: 4.1.4
 Vcs-Browser: https://salsa.debian.org/med-team/iqtree
 Vcs-Git: https://salsa.debian.org/med-team/iqtree.git

debian/rules

@@ -35,6 +35,10 @@ dh_auto_build_%:
 dh_auto_install_%:
 	dh_auto_install -Bbuild.$(subst dh_auto_install_,,$@)
 
+override_dh_install:
+	dh_install
+	find debian -name iqtree-mpi -exec chrpath --delete \{\} \;
+
 dh_auto_clean_%:
 	dh_auto_clean -Bbuild.$(subst dh_auto_clean_,,$@)
 

main/main.cpp

@@ -1635,7 +1635,7 @@ void compare(Params &params){
 /**MINH ANH: to compute 'guided bootstrap' alignment*/
 void guidedBootstrap(Params &params)
 {
-	MaAlignment inputAlign(params.aln_file,params.sequence_type, params.intype);
+	MaAlignment inputAlign(params.aln_file,params.sequence_type, params.intype, params.model_name);
 	inputAlign.readLogLL(params.siteLL_file);
 
 	string outFre_name = params.out_prefix;
@@ -1673,8 +1673,8 @@ void guidedBootstrap(Params &params)
 /**MINH ANH: to compute the probability of an alignment given the multinomial distribution of patterns frequencies derived from a reference alignment*/
 void computeMulProb(Params &params)
 {
-	Alignment refAlign(params.second_align, params.sequence_type, params.intype);
-	Alignment inputAlign(params.aln_file, params.sequence_type, params.intype);
+	Alignment refAlign(params.second_align, params.sequence_type, params.intype, params.model_name);
+	Alignment inputAlign(params.aln_file, params.sequence_type, params.intype, params.model_name);
 	double prob;
 	inputAlign.multinomialProb(refAlign,prob);
 	//Printing

main/phyloanalysis.h

@@ -38,8 +38,10 @@ class IQTree;
 */
 void runPhyloAnalysis(Params &params, Checkpoint *checkpoint);
 
-void runTreeReconstruction(Params &params, string &original_model,
-		IQTree* &tree, ModelCheckpoint &model_info);
+void startTreeReconstruction(Params &params, IQTree* &iqtree,
+        ModelCheckpoint &model_info);
+
+void runTreeReconstruction(Params &params, IQTree* &tree);
 
 /**
 	take the collection of trees from input_trees, it assign support values to target_tree

main/phylotesting.h

@@ -145,6 +145,14 @@ string criterionName(ModelTestCriterion mtc);
 
 bool checkModelFile(string model_file, bool is_partitioned, ModelCheckpoint &infos);
 
+/**
+ perform ModelFinder to find the best-fit model
+ @param params program parameters
+ @param iqtree phylogenetic tree
+ @param model_info (IN/OUT) information for all models considered
+ */
+void runModelFinder(Params &params, IQTree &iqtree, ModelCheckpoint &model_info);
+
 /**
  testing the best-fit model
  return in params.freq_type and params.rate_type
@@ -155,9 +163,9 @@ bool checkModelFile(string model_file, bool is_partitioned, ModelCheckpoint &inf
  @param print_mem_usage true to print RAM memory used (default: false) 
  @return name of best-fit-model
  */
-string testModel(Params &params, PhyloTree* in_tree, ModelCheckpoint &model_info,
-		ModelsBlock *models_block, int num_threads, int brlen_type,
-        string set_name = "", bool print_mem_usage = false, string in_model_name = "");
+//string testModel(Params &params, PhyloTree* in_tree, ModelCheckpoint &model_info,
+//		ModelsBlock *models_block, int num_threads, int brlen_type,
+//        string set_name = "", bool print_mem_usage = false, string in_model_name = "");
 
 /**
  * print site log likelihoods to a fileExists
@@ -234,9 +242,9 @@ void evaluateTrees(Params &params, IQTree *tree);
     @param seq_type (OUT) sequence type, SEQ_UNKNOWN if is not determined
     @return 1 for parametric model, 2 for empirical model
 */
-int getSeqType(const char *model_name, SeqType &seq_type);
+int detectSeqType(const char *model_name, SeqType &seq_type);
 
-string getSeqType(string model_name);
+string detectSeqTypeName(string model_name);
 
 
 #endif /* PHYLOTESTING_H_ */

model/modelcodon.cpp

@@ -10,7 +10,7 @@
 
 
 const double MIN_OMEGA_KAPPA = 0.001;
-const double MAX_OMEGA_KAPPA = 100.0;
+const double MAX_OMEGA_KAPPA = 50.0;
 
 /* Empirical codon model restricted (Kosiol et al. 2007), source: http://www.ebi.ac.uk/goldman/ECM/ */
 string model_ECMrest1 =

model/modeldna.cpp

@@ -183,6 +183,8 @@ void ModelDNA::init(const char *model_name, string model_params, StateFreqType f
 		    full_name = "Time reversible ("+name+")";
 		} else {
 			readParameters(model_name);
+            name = full_name = model_name;
+            freq = FREQ_USER_DEFINED;
 			//name += " (user-defined)";
 		}
 	}
@@ -225,7 +227,10 @@ void ModelDNA::restoreCheckpoint() {
     string rate_spec = param_spec;
     for (auto i = rate_spec.begin(); i != rate_spec.end(); i++)
         *i = *i + '0';
-    ASSERT(setRateType(rate_spec));
+    
+    if (!rate_spec.empty())
+        if (!setRateType(rate_spec))
+            ASSERT(0 && "Cannot set rate_spec");
 
     decomposeRateMatrix();
     if (phylo_tree)
@@ -239,24 +244,25 @@ void ModelDNA::readRates(string str) throw(const char*) {
 	for (j = 0; j < param_spec.length(); j++)
 		rates[j] = 1.0;
 	num_params = 0;
-	for (i = 0; i < nrates && end_pos < str.length(); i++) {
+	for (i = 0; i <= nrates && end_pos < str.length(); i++) {
 		int new_end_pos;
 		double rate = 0;
+        int id = (i < nrates) ? i+1 : 0;
 		if (str[end_pos] == '?') {
-			param_fixed[i+1] = false;
+			param_fixed[id] = false;
 			end_pos++;
 			rate = 1.0;
 			num_params++;
 		} else {
             if (Params::getInstance().optimize_rate_matrix) {
                 num_params++;
-                param_fixed[i+1] = false;
+                param_fixed[id] = false;
             } else
 	        if (Params::getInstance().optimize_from_given_params) {
                     num_params++;
-                    param_fixed[i+1] = false;
+                    param_fixed[id] = false;
 	        } else {
-	            param_fixed[i+1] = true;
+	            param_fixed[id] = true;
 	        }
 			try {
 				rate = convert_double(str.substr(end_pos).c_str(), new_end_pos);
@@ -267,7 +273,7 @@ void ModelDNA::readRates(string str) throw(const char*) {
 		}
 		if (rate < 0.0)
 			outError("Negative rates found");
-		if (i == nrates-1 && end_pos < str.length())
+		if (i == nrates && end_pos < str.length())
 			outError("More than " + convertIntToString(nrates) + " rate parameters specified in " + str);
 		if (i < nrates-1 && end_pos >= str.length())
 			outError("Unexpected end of string ", str);
@@ -275,7 +281,7 @@ void ModelDNA::readRates(string str) throw(const char*) {
 			outError("Comma to separate rates not found in ", str);
 		end_pos++;
 		for (j = 0; j < param_spec.length(); j++)
-			if (param_spec[j] == i+1)
+			if (param_spec[j] == id)
 				rates[j] = rate;
 	}
 }
@@ -360,7 +366,10 @@ bool ModelDNA::setRateType(string rate_str) {
 	for (i = 0; i <= num_params; i++)
 		avg_rates[i] /= num_rates[i];
 	for (i = 0; i < param_spec.size(); i++) {
+        if (avg_rates[0] > 0.0)
             rates[i] = avg_rates[(int)param_spec[i]] / avg_rates[0];
+        else
+            rates[i] = avg_rates[(int)param_spec[i]];
 	}
 	if (verbose_mode >= VB_DEBUG) {
 		cout << "Initialized rates: ";

model/modelfactory.cpp

@@ -44,17 +44,24 @@
 #include "nclextra/myreader.h"
 #include <sstream>
 
-string::size_type posRateHeterotachy(string &model_name) {
-	string::size_type pos1 = 0, pos2 = 0;
-    do {
-        pos1 = model_name.find("+H", pos1+1);
-        if (pos1 == string::npos) break;
-    } while (pos1 < model_name.length()-2 && isalpha(model_name[pos1+2]));
+string::size_type findSubStr(string &name, string sub1, string sub2) {
+    string::size_type pos1, pos2;
+    for (pos1 = 0; pos1 != string::npos; pos1++) {
+        pos1 = name.find(sub1, pos1);
+        if (pos1 == string::npos)
+            break;
+        if (pos1+2 >= name.length() || !isalpha(name[pos1+2])) {
+            break;
+        }
+    }
     
-    do {
-        pos2 = model_name.find("*H", pos2+1);
-        if (pos2 == string::npos) break;
-    } while (pos2 < model_name.length()-2 && isalpha(model_name[pos2+2]));
+    for (pos2 = 0; pos2 != string::npos; pos2++) {
+        pos2 = name.find(sub2, pos2);
+        if (pos2 == string::npos)
+            break;
+        if (pos2+2 >= name.length() ||!isalpha(name[pos2+2]))
+            break;
+    }
     
     if (pos1 != string::npos && pos2 != string::npos) {
         return min(pos1, pos2);
@@ -62,49 +69,18 @@ string::size_type posRateHeterotachy(string &model_name) {
         return pos1;
     else
         return pos2;
+}
 
+string::size_type posRateHeterotachy(string &model_name) {
+    return findSubStr(model_name, "+H", "*H");
 }
 
 string::size_type posRateFree(string &model_name) {
-	string::size_type pos1 = 0, pos2 = 0;
-    do {
-        pos1 = model_name.find("+R", pos1+1);
-        if (pos1 == string::npos) break;
-    } while (pos1 < model_name.length()-2 && isalpha(model_name[pos1+2]));
-
-    do {
-        pos2 = model_name.find("*R", pos2+1);
-        if (pos2 == string::npos) break;
-    } while (pos2 < model_name.length()-2 && isalpha(model_name[pos2+2]));
-
-    if (pos1 != string::npos && pos2 != string::npos) {
-        return min(pos1, pos2);
-    } else if (pos1 != string::npos)
-        return pos1;
-    else
-        return pos2;
-
+    return findSubStr(model_name, "+R", "*R");
 }
 
 string::size_type posPOMO(string &model_name) {
-	string::size_type pos1 = 0, pos2 = 0;
-    do {
-        pos1 = model_name.find("+P", pos1+1);
-        if (pos1 == string::npos) break;
-    } while (pos1 < model_name.length()-2 && isalpha(model_name[pos1+2]));
-
-    do {
-        pos2 = model_name.find("*P", pos2+1);
-        if (pos2 == string::npos) break;
-    } while (pos2 < model_name.length()-2 && isalpha(model_name[pos2+2]));
-
-    if (pos1 != string::npos && pos2 != string::npos) {
-        return min(pos1, pos2);
-    } else if (pos1 != string::npos)
-        return pos1;
-    else
-        return pos2;
-
+    return findSubStr(model_name, "+P", "*P");
 }
 
 ModelsBlock *readModelsDefinition(Params &params) {
@@ -242,7 +218,8 @@ ModelFactory::ModelFactory(Params &params, string &model_name, PhyloTree *tree,
 				outError("Model name has wrong bracket notation '{...}'");
 				rate_str = model_str.substr(pos+1);
 				model_str = model_str.substr(0, pos+1);
-        } else {
+        }
+    else {
             rate_str = model_str.substr(spec_pos);
             model_str = model_str.substr(0, spec_pos);
         }
@@ -304,7 +281,6 @@ ModelFactory::ModelFactory(Params &params, string &model_name, PhyloTree *tree,
     // models, the heterozygosity can be set separately for each model and the
     // +P{}, +GXX and +I flags should already be inside the model definition.
     if (model_str.substr(0, 3) != "MIX" && pomo) {
-
       // +P{} flag.
       p_pos = posPOMO(rate_str);
       if (p_pos != string::npos) {
@@ -369,7 +345,7 @@ ModelFactory::ModelFactory(Params &params, string &model_name, PhyloTree *tree,
 	if (freq_type == FREQ_UNKNOWN) {
 		switch (tree->aln->seq_type) {
 		case SEQ_BINARY: freq_type = FREQ_ESTIMATE; break; // default for binary: optimized frequencies
-		case SEQ_PROTEIN: freq_type = FREQ_USER_DEFINED; break; // default for protein: frequencies of the empirical AA matrix
+		case SEQ_PROTEIN: break; // let ModelProtein decide by itself
 		case SEQ_MORPH: freq_type = FREQ_EQUAL; break;
 		case SEQ_CODON: freq_type = FREQ_UNKNOWN; break;
             break;
@@ -576,6 +552,12 @@ ModelFactory::ModelFactory(Params &params, string &model_name, PhyloTree *tree,
 	if ((posasc = rate_str.find("+ASC")) != string::npos) {
 		// ascertainment bias correction
 		unobserved_ptns = tree->aln->getUnobservedConstPatterns();
+        
+        // delete rarely observed state
+        for (int i = unobserved_ptns.length()-1; i >= 0; i--)
+            if (model->state_freq[(int)unobserved_ptns[i]] < 1e-8)
+                unobserved_ptns.erase(i);
+                
 		// rebuild the seq_states to contain states of unobserved constant patterns
 		tree->aln->buildSeqStates(true);
 //		if (unobserved_ptns.size() <= 0)

model/modelprotein.cpp

@@ -579,13 +579,16 @@ void ModelProtein::init(const char *model_name, string model_params, StateFreqTy
     } else if (name_upper == "GTR20") {
         outWarning("GTR20 model will estimate 189 substitution rates that might be overfitting!");
         outWarning("Please only use GTR20 with very large data and always test for model fit!");
-        if (freq == FREQ_UNKNOWN || freq == FREQ_USER_DEFINED)
+        if (freq == FREQ_UNKNOWN)
             freq = FREQ_EMPIRICAL;
         // initialize rate matrix with LG
         nxs_model = models_block->findModel("LG");
         ASSERT(nxs_model);
         readParametersString(nxs_model->description);
         rescaleRates(rates, getNumRateEntries());
+
+        // 2018-05-08 bug fix: GTR20 rates are not optimized
+        num_params = getNumRateEntries()-1;
 	} else {
 		// if name does not match, read the user-defined model
 		readParameters(model_name);

model/partitionmodel.cpp

@@ -50,7 +50,7 @@ PartitionModel::PartitionModel(Params &params, PhyloSuperTree *tree, ModelsBlock
     }
     for (it = tree->begin(), part = 0; it != tree->end(); it++, part++) {
         ASSERT(!((*it)->getModelFactory()));
-        string model_name = tree->part_info[part].model_name;
+        string model_name = (*it)->aln->model_name;
         if (model_name == "") // if empty, take model name from command option
         	model_name = params.model_name;
         (*it)->setModelFactory(new ModelFactory(params, model_name, (*it), models_block));
@@ -86,7 +86,7 @@ void PartitionModel::saveCheckpoint() {
     PhyloSuperTree *tree = (PhyloSuperTree*)site_rate->getTree();
     int part = 0;
     for (PhyloSuperTree::iterator it = tree->begin(); it != tree->end(); it++, part++) {
-        checkpoint->startStruct(tree->part_info[part].name);
+        checkpoint->startStruct((*it)->aln->name);
         (*it)->getModelFactory()->saveCheckpoint();
         checkpoint->endStruct();
     }
@@ -103,7 +103,7 @@ void PartitionModel::restoreCheckpoint() {
     PhyloSuperTree *tree = (PhyloSuperTree*)site_rate->getTree();
     int part = 0;
     for (PhyloSuperTree::iterator it = tree->begin(); it != tree->end(); it++, part++) {
-        checkpoint->startStruct(tree->part_info[part].name);
+        checkpoint->startStruct((*it)->aln->name);
         (*it)->getModelFactory()->restoreCheckpoint();
         checkpoint->endStruct();
     }
@@ -165,7 +165,7 @@ double PartitionModel::optimizeParameters(int fixed_len, bool write_info, double
         #endif
         {
     		cout << "Optimizing " << tree->at(part)->getModelName() <<
-        		" parameters for partition " << tree->part_info[part].name <<
+        		" parameters for partition " << tree->at(part)->aln->name <<
         		" (" << tree->at(part)->getModelFactory()->getNParameters(fixed_len) << " free parameters)" << endl;
         }
         tree_lh += tree->at(part)->getModelFactory()->optimizeParameters(fixed_len, write_info && verbose_mode >= VB_MED, 
@@ -199,7 +199,7 @@ double PartitionModel::optimizeParametersGammaInvar(int fixed_len, bool write_in
         #endif
         {
     		cout << "Optimizing " << tree->at(part)->getModelName() <<
-        		" parameters for partition " << tree->part_info[part].name <<
+        		" parameters for partition " << tree->at(part)->aln->name <<
         		" (" << tree->at(part)->getModelFactory()->getNParameters(fixed_len) << " free parameters)" << endl;
         }
         tree_lh += tree->at(part)->getModelFactory()->optimizeParametersGammaInvar(fixed_len, write_info && verbose_mode >= VB_MED, 

nclextra/msetsblock.cpp

@@ -18,6 +18,7 @@
  *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
  ***************************************************************************/
 #include "msetsblock.h"
+#include "utils/tools.h"
 
 MSetsBlock::MSetsBlock()
  : NxsBlock()
@@ -105,7 +106,9 @@ void MSetsBlock::Read(NxsToken &token)
 			token.SetLabileFlagBit(NxsToken::preserveUnderscores);
 			token.GetNextToken();
 			do {
-				myset->taxlist.push_back(token.GetToken());
+                string taxname = token.GetToken();
+                renameString(taxname);
+				myset->taxlist.push_back(taxname);
 				token.SetLabileFlagBit(NxsToken::preserveUnderscores);
 				token.GetNextToken();
 			} while (!token.AtEOF() && !token.Equals(";"));