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 "5")
+set (iqtree_VERSION_PATCH "7.1")
 
 set(BUILD_SHARED_LIBS OFF)
 
@@ -513,12 +513,21 @@ if (IQTREE_FLAGS MATCHES "mpi")
 	add_library(mympi utils/TreeCollection.cpp utils/ObjectStream.cpp)
 endif()
 
+if (NOT IQTREE_FLAGS MATCHES "single")
+    if (BINARY32)
+        link_directories(${PROJECT_SOURCE_DIR}/lib32)
+    else()
+        link_directories(${PROJECT_SOURCE_DIR}/lib)
+    endif()
+endif()
+
 add_executable(iqtree
 main/main.cpp
 main/phyloanalysis.cpp
 main/phyloanalysis.h
 main/phylotesting.cpp
 main/phylotesting.h
+obsolete/parsmultistate.cpp
 )
 
 if(Backtrace_FOUND)
@@ -612,7 +621,7 @@ endif()
 set_target_properties(iqtree PROPERTIES OUTPUT_NAME "iqtree${EXE_SUFFIX}")
 
 # strip the release build
-if (CMAKE_BUILD_TYPE STREQUAL "Release" AND (GCC OR CLANG) AND NOT APPLE) # strip is not necessary for MSVC
+if (NOT IQTREE_FLAGS MATCHES "nostrip" AND CMAKE_BUILD_TYPE STREQUAL "Release" AND (GCC OR CLANG) AND NOT APPLE) # strip is not necessary for MSVC
 	if (WIN32)
 		ADD_CUSTOM_COMMAND(TARGET iqtree POST_BUILD COMMAND strip $<TARGET_FILE:iqtree>)
 	elseif (NOT APPLE)

alignment/alignment.cpp

@@ -532,18 +532,26 @@ int Alignment::readNexus(char *filename) {
 		return 0;
 	}
 
-    if (char_block->GetNTax() == 0) { char_block = data_block; }
-
-    if (char_block->GetNTax() == 0) {
-        outError("No data is given in the input file");
+    if (data_block->GetNTax() == 0 && char_block->GetNTax() == 0) {
+        outError("No DATA or CHARACTERS blocks found");
         return 0;
     }
-    if (verbose_mode >= VB_DEBUG)
-        char_block->Report(cout);
-
 
+    if (char_block->GetNTax() > 0) {
         extractDataBlock(char_block);
+        if (verbose_mode >= VB_DEBUG)
+            char_block->Report(cout);
+    } else {
+        extractDataBlock(data_block);
+        if (verbose_mode >= VB_DEBUG)
+            data_block->Report(cout);
+    }
 
+    delete trees_block;
+    delete char_block;
+    delete data_block;
+    delete assumptions_block;
+    delete taxa_block;
     return 1;
 }
 
@@ -567,8 +575,8 @@ int getDataBlockMorphStates(NxsCharactersBlock *data_block) {
     char ch;
     int nstates = 0;
 
-    for (site = 0; site < nsite; site++)
-        for (seq = 0; seq < nseq; seq++) {
+    for (seq = 0; seq < nseq; seq++)
+        for (site = 0; site < nsite; site++) {
             int nstate = data_block->GetNumStates(seq, site);
             if (nstate == 0)
                 continue;
@@ -580,23 +588,11 @@ int getDataBlockMorphStates(NxsCharactersBlock *data_block) {
                 else if (ch >= 'A' && ch <= 'Z')
                     ch = ch - 'A' + 11;
                 else
-                    outError(data_block->GetTaxonLabel(seq) + " has invalid state at site " + convertIntToString(site));
+                    outError(data_block->GetTaxonLabel(seq) + " has invalid single state " + ch + " at site " + convertIntToString(site+1));
                 if (ch > nstates) nstates = ch;
                 continue;
             }
-            for (int state = 0; state < nstate; state++) {
-                ch = data_block->GetState(seq, site, state);
-                if (!isalnum(ch)) continue;
-                if (ch >= '0' && ch <= '9') ch = ch - '0' + 1;
-                if (ch >= 'A' && ch <= 'Z') ch = ch - 'A' + 11;
-                if (ch >= '0' && ch <= '9')
-                    ch = ch - '0' + 1;
-                else if (ch >= 'A' && ch <= 'Z')
-                    ch = ch - 'A' + 11;
-                else
-                    outError(data_block->GetTaxonLabel(seq) + " has invalid state at site " + convertIntToString(site));
-                if (ch > nstates) nstates = ch;
-            }
+            //cout << "NOTE: " << data_block->GetTaxonLabel(seq) << " has ambiguous state at site " << site+1 << " which is treated as unknown" << endl;
         }
     return nstates;
 }
@@ -609,6 +605,9 @@ void Alignment::extractDataBlock(NxsCharactersBlock *data_block) {
     char char_to_state[NUM_CHAR];
     char state_to_char[NUM_CHAR];
     
+    if (!data_block->GetMatrix())
+        outError("MATRIX command undeclared or invalid");
+    
     NxsCharactersBlock::DataTypesEnum data_type = (NxsCharactersBlock::DataTypesEnum)data_block->GetDataType();
     if (data_type == NxsCharactersBlock::continuous) {
         outError("Continuous characters not supported");
@@ -674,14 +673,24 @@ void Alignment::extractDataBlock(NxsCharactersBlock *data_block) {
                 pat.push_back(STATE_UNKNOWN);
             else if (nstate == 1) {
                 pat.push_back(char_to_state[(int)data_block->GetState(seq, site, 0)]);
-            } else {
-                ASSERT(data_type != NxsCharactersBlock::dna || data_type != NxsCharactersBlock::rna || data_type != NxsCharactersBlock::nucleotide);
+            } else if (data_type == NxsCharactersBlock::dna || data_type == NxsCharactersBlock::rna || data_type == NxsCharactersBlock::nucleotide) {
+                // 2018-06-07: correctly interpret ambiguous nucleotide
                 char pat_ch = 0;
                 for (int state = 0; state < nstate; state++) {
                     pat_ch |= (1 << char_to_state[(int)data_block->GetState(seq, site, state)]);
                 }
                 pat_ch += 3;
                 pat.push_back(pat_ch);
+            } else {
+                // other ambiguous characters are treated as unknown
+                stringstream str;
+                str << "Sequence " << seq_names[seq] << " site " << site+1 << ": {";
+                for (int state = 0; state < nstate; state++) {
+                    str << data_block->GetState(seq, site, state);
+                }
+                str << "} treated as unknown character";
+                outWarning(str.str());
+                pat.push_back(STATE_UNKNOWN);
             }
         }
         num_gaps_only += addPattern(pat, site);
@@ -1103,6 +1112,7 @@ void Alignment::buildStateMap(char *map, SeqType seq_type) {
         map[(unsigned char)'J'] = 22; // I or L
         map[(unsigned char)'*'] = STATE_UNKNOWN; // stop codon
         map[(unsigned char)'U'] = STATE_UNKNOWN; // 21st amino acid
+        map[(unsigned char)'O'] = STATE_UNKNOWN; // 22nd amino acid
 
         return;
     case SEQ_MULTISTATE:
@@ -1191,6 +1201,7 @@ char Alignment::convertState(char state, SeqType seq_type) {
 		if (state == 'J') return 22;
         if (state == '*') return STATE_UNKNOWN; // stop codon
         if (state == 'U') return STATE_UNKNOWN; // 21st amino-acid
+        if (state == 'O') return STATE_UNKNOWN; // 22nd amino-acid
         loc = strchr(symbols_protein, state);
 
         if (!loc) return STATE_INVALID; // unrecognize character
@@ -1287,13 +1298,13 @@ char Alignment::convertStateBack(char state) {
     }
 }
 
-string Alignment::convertStateBackStr(char state) {
+string Alignment::convertStateBackStr(StateType state) {
 	string str;
     if (seq_type == SEQ_POMO)
         return string("POMO")+convertIntToString(state);
-	if (seq_type != SEQ_CODON) {
-		str = convertStateBack(state);
-	} else {
+    if (seq_type == SEQ_MULTISTATE)
+        return " " + convertIntToString(state);
+	if (seq_type == SEQ_CODON) {
         // codon data
         if (state >= num_states) return "???";
         assert(codon_table);
@@ -1301,7 +1312,10 @@ string Alignment::convertStateBackStr(char state) {
         str = symbols_dna[state/16];
         str += symbols_dna[(state%16)/4];
         str += symbols_dna[state%4];
+        return str;
 	}
+    // all other data types
+    str = convertStateBack(state);
 	return str;
 }
 
@@ -1395,9 +1409,9 @@ SeqType Alignment::getSeqType(const char *sequence_type) {
         user_seq_type = SEQ_PROTEIN;
     } else if (strncmp(sequence_type, "NT2AA", 5) == 0) {
         user_seq_type = SEQ_PROTEIN;
-    } else if (strcmp(sequence_type, "NUM") == 0 || strcmp(sequence_type, "MORPH") == 0 || strcmp(sequence_type, "MULTI") == 0) {
+    } else if (strcmp(sequence_type, "NUM") == 0 || strcmp(sequence_type, "MORPH") == 0) {
         user_seq_type = SEQ_MORPH;
-    } else if (strcmp(sequence_type, "TINA") == 0) {
+    } else if (strcmp(sequence_type, "TINA") == 0 || strcmp(sequence_type, "MULTI") == 0) {
         user_seq_type = SEQ_MULTISTATE;
     } else if (strncmp(sequence_type, "CODON", 5) == 0) {
         user_seq_type = SEQ_CODON;
@@ -1492,11 +1506,11 @@ int Alignment::buildPattern(StrVector &sequences, char *sequence_type, int nseq,
             num_states = 20;
             nt2aa = true;
             cout << "Translating to amino-acid sequences with genetic code " << &sequence_type[5] << " ..." << endl;
-        } else if (strcmp(sequence_type, "NUM") == 0 || strcmp(sequence_type, "MORPH") == 0 || strcmp(sequence_type, "MULTI") == 0) {
+        } else if (strcmp(sequence_type, "NUM") == 0 || strcmp(sequence_type, "MORPH") == 0) {
             num_states = getMorphStates(sequences);
             if (num_states < 2 || num_states > 32) throw "Invalid number of states";
             user_seq_type = SEQ_MORPH;
-        } else if (strcmp(sequence_type, "TINA") == 0) {
+        } else if (strcmp(sequence_type, "TINA") == 0 || strcmp(sequence_type, "MULTI") == 0) {
             cout << "Multi-state data with " << num_states << " alphabets" << endl;
             user_seq_type = SEQ_MULTISTATE;
         } else if (strncmp(sequence_type, "CODON", 5) == 0) {
@@ -1623,7 +1637,7 @@ int Alignment::readPhylip(char *filename, char *sequence_type) {
     string line;
     // remove the failbit
     in.exceptions(ios::badbit);
-    bool tina_state = (sequence_type && strcmp(sequence_type,"TINA") == 0);
+    bool tina_state = (sequence_type && (strcmp(sequence_type,"TINA") == 0 || strcmp(sequence_type,"MULTI") == 0));
     num_states = 0;
 
     for (; !in.eof(); line_num++) {
@@ -2924,6 +2938,8 @@ void convert_range(const char *str, int &lower, int &upper, int &step_size, char
     //int d_save = d;
     upper = d;
     step_size = 1;
+    // skip blank chars
+    for (; *endptr == ' '; endptr++) {}
     if (*endptr != '-') return;
 
     // parse the upper bound of the range
@@ -2938,6 +2954,9 @@ void convert_range(const char *str, int &lower, int &upper, int &step_size, char
 
     //lower = d_save;
     upper = d;
+    // skip blank chars
+    for (; *endptr == ' '; endptr++) {}
+
     if (*endptr != '\\') return;
 
     // parse the step size of the range
@@ -2977,7 +2996,7 @@ void extractSiteID(Alignment *aln, const char* spec, IntVector &site_id) {
             for (i = lower; i <= upper; i+=step)
                 site_id.push_back(i);
             if (*str == ',' || *str == ' ') str++;
-            else break;
+            //else break;
         }
         if (aln->seq_type == SEQ_CODON && nchars % 3 != 0)
             throw (string)"Range " + spec + " length is not multiple of 3 (necessary for codon data)";
@@ -3044,18 +3063,11 @@ void Alignment::createBootstrapAlignment(Alignment *aln, IntVector* pattern_freq
     if (!spec) {
 		// standard bootstrap
         int added_sites = 0;
-		for (site = 0; site < nsite; site++) {
-            int site_id;
-            if (Params::getInstance().jackknife_prop == 0.0) {
-                // bootstrap sampling with replacement
-                site_id = random_int(nsite);
-            } else {
-                // jacknife without replacement
-                if (random_double() < Params::getInstance().jackknife_prop)
-                    continue;
-                site_id = site;
-            }
-			int ptn_id = aln->getPatternID(site_id);
+        IntVector sample;
+        random_resampling(nsite, sample);
+		for (site = 0; site < nsite; site++)
+        for (int rep = 0; rep < sample[site]; rep++) {
+			int ptn_id = aln->getPatternID(site);
 			Pattern pat = aln->at(ptn_id);
             int nptn = getNPattern();
 			addPattern(pat, added_sites);
@@ -3165,26 +3177,16 @@ void Alignment::createBootstrapAlignment(int *pattern_freq, const char *spec, in
     if (Params::getInstance().jackknife_prop > 0.0 && spec)
         outError((string)"Unsupported jackknife with " + spec);
 
-    if (!spec ||  strncmp(spec, "SCALE=", 6) == 0) {
+    if (!spec) {
 
-        if (spec) {
-            double scale = convert_double(spec+6);
-            nsite = (int)round(scale * nsite);
-        }
         int nptn = getNPattern();
 
         if (nsite/8 < nptn || Params::getInstance().jackknife_prop > 0.0) {
-            int orig_nsite = getNSite();
-            for (site = 0; site < nsite; site++) {
-                int site_id;
-                if (Params::getInstance().jackknife_prop == 0.0)
-                    site_id = random_int(orig_nsite, rstream);
-                else {
-                    if (random_double() < Params::getInstance().jackknife_prop)
-                        continue;
-                    site_id = site;
-                }
-                int ptn_id = getPatternID(site_id);
+            IntVector sample;
+            random_resampling(nsite, sample, rstream);
+            for (site = 0; site < nsite; site++)
+            for (int rep = 0; rep < sample[site]; rep++) {
+                int ptn_id = getPatternID(site);
                 pattern_freq[ptn_id]++;
             }
         } else {
@@ -3243,7 +3245,11 @@ void Alignment::createBootstrapAlignment(int *pattern_freq, const char *spec, in
 		}
 	} else {
 		// resampling sites within genes
+        try {
             convert_int_vec(spec, site_vec);
+        } catch (...) {
+            outError("-bsam not allowed for non-partition model");
+        }
 		if (site_vec.size() % 2 != 0)
 			outError("Bootstrap specification length is not divisible by 2");
 		int part, begin_site = 0, out_site = 0;

alignment/alignment.h

@@ -258,7 +258,7 @@ public:
 	 * @param state internal state code
 	 * @return user-readable state string
 	 */
-	string convertStateBackStr(char state);
+	string convertStateBackStr(StateType state);
 
 	/**
             get alignment site range from the residue range relative to a sequence

alignment/superalignment.cpp

@@ -799,6 +799,7 @@ void SuperAlignment::printSiteInfo(const char* filename) {
     }
 }
 
+/*
 void SuperAlignment::createBootstrapAlignment(Alignment *aln, IntVector* pattern_freq, const char *spec) {
 	ASSERT(aln->isSuperAlignment());
 	Alignment::copyAlignment(aln);
@@ -881,6 +882,71 @@ void SuperAlignment::createBootstrapAlignment(Alignment *aln, IntVector* pattern
 	taxa_index = super_aln->taxa_index;
     countConstSite();
 }
+*/
+
+void SuperAlignment::createBootstrapAlignment(Alignment *aln, IntVector* pattern_freq, const char *spec) {
+    ASSERT(aln->isSuperAlignment());
+    SuperAlignment *super_aln = (SuperAlignment*) aln;
+    ASSERT(partitions.empty());
+    name = aln->name;
+    model_name = aln->model_name;
+    sequence_type = aln->sequence_type;
+    position_spec = aln->position_spec;
+    aln_file = aln->aln_file;
+    
+    if (!spec) {
+        // resampling sites within genes
+        Alignment::copyAlignment(aln);
+        partitions.reserve(super_aln->partitions.size());
+        for (vector<Alignment*>::iterator it = super_aln->partitions.begin(); it != super_aln->partitions.end(); it++) {
+            Alignment *boot_aln = new Alignment;
+            if (pattern_freq) {
+                IntVector part_pattern_freq;
+                boot_aln->createBootstrapAlignment(*it, &part_pattern_freq);
+                pattern_freq->insert(pattern_freq->end(), part_pattern_freq.begin(), part_pattern_freq.end());
+            } else {
+                boot_aln->createBootstrapAlignment(*it);
+            }
+            partitions.push_back(boot_aln);
+        }
+        taxa_index = super_aln->taxa_index;
+        countConstSite();
+    } else if (strcmp(spec, "GENE") == 0) {
+        ASSERT(!pattern_freq);
+        // resampling whole genes
+        IntVector gene_freq;
+        random_resampling(super_aln->partitions.size(), gene_freq);
+        for (int i = 0; i < gene_freq.size(); i++)
+            if (gene_freq[i] > 0) {
+                Alignment *boot_aln = new Alignment;
+                boot_aln->copyAlignment(super_aln->partitions[i]);
+                if (gene_freq[i] > 1) {
+                    for (auto it = boot_aln->begin(); it != boot_aln->end(); it++)
+                        it->frequency *= gene_freq[i];
+                    auto site_pattern = boot_aln->site_pattern;
+                    for (int j = 1; j < gene_freq[i]; j++)
+                        boot_aln->site_pattern.insert(boot_aln->site_pattern.end(), site_pattern.begin(), site_pattern.end());
+                }
+                partitions.push_back(boot_aln);
+            }
+        init();
+    } else if (strcmp(spec, "GENESITE") == 0) {
+        ASSERT(!pattern_freq);
+        // resampling whole genes then sites within resampled genes
+        IntVector gene_freq;
+        random_resampling(super_aln->partitions.size(), gene_freq);
+        for (int i = 0; i < gene_freq.size(); i++)
+            for (int rep = 0; rep < gene_freq[i]; rep++) {
+            Alignment *boot_aln = new Alignment;
+            boot_aln->createBootstrapAlignment(super_aln->partitions[i]);
+            boot_aln->name = boot_aln->name + "." + convertIntToString(rep);
+            partitions.push_back(boot_aln);
+        }
+        init();
+    } else {
+        outError("Wrong -bsam, either -bsam GENE or -bsam GENESITE");
+    }
+}
 
 void SuperAlignment::createBootstrapAlignment(IntVector &pattern_freq, const char *spec) {
 	ASSERT(isSuperAlignment());
@@ -910,28 +976,29 @@ void SuperAlignment::createBootstrapAlignment(int *pattern_freq, const char *spe
 			nptn += (*it)->getNPattern();
 		}
 		memset(pattern_freq, 0, nptn * sizeof(int));
-		for (int i = 0; i < partitions.size(); i++) {
-			int part = random_int(partitions.size(), rstream);
+        IntVector gene_freq;
+        random_resampling(partitions.size(), gene_freq, rstream);
+		for (int part = 0; part < partitions.size(); part++)
+        for (int rep = 0; rep < gene_freq[part]; rep++){
 			Alignment *aln = partitions[part];
 			if (strncmp(spec,"GENESITE",8) == 0) {
 				// then resampling sites in resampled gene
-				for (int j = 0; j < aln->getNSite(); j++) {
-					int ptn_id = aln->getPatternID(random_int(aln->getNPattern(), rstream));
+                IntVector sample;
+                random_resampling(aln->getNSite(), sample, rstream);
+				for (int site = 0; site < sample.size(); site++)
+                for (int rep2 = 0; rep2 < sample[site]; rep2++) {
+					int ptn_id = aln->getPatternID(site);
 					pattern_freq[ptn_id + part_pos[part]]++;
 				}
-
 			} else {
-				for (int j = 0; j < aln->getNPattern(); j++)
-					pattern_freq[j + part_pos[part]] += aln->at(j).frequency;
+				for (int ptn = 0; ptn < aln->getNPattern(); ptn++)
+					pattern_freq[ptn + part_pos[part]] += aln->at(ptn).frequency;
 			}
 		}
 	} else {
 		// resampling sites within genes
 		int offset = 0;
 		for (vector<Alignment*>::iterator it = partitions.begin(); it != partitions.end(); it++) {
-            if (spec && strncmp(spec, "SCALE=", 6) == 0)
-                (*it)->createBootstrapAlignment(pattern_freq + offset, spec, rstream);
-            else
             (*it)->createBootstrapAlignment(pattern_freq + offset, NULL, rstream);
 			offset += (*it)->getNPattern();
 		}
@@ -1148,21 +1215,25 @@ Alignment *SuperAlignment::concatenateAlignments(set<int> &ids) {
     int site = 0;
     for (it = ids.begin(); it != ids.end(); it++) {
     	int id = *it;
-		string taxa_set;
-        Pattern taxa_pat = getPattern(id);
-        taxa_set.insert(taxa_set.begin(), taxa_pat.begin(), taxa_pat.end());
+        // 2018-08-23: important bugfix in v1.6: taxa_set has wrong correspondance
+		//string taxa_set;
+        //Pattern taxa_pat = getPattern(id);
+        //taxa_set.insert(taxa_set.begin(), taxa_pat.begin(), taxa_pat.end());
     	for (Alignment::iterator it = partitions[id]->begin(); it != partitions[id]->end(); it++) {
     		Pattern pat;
-    		int part_seq = 0;
+    		//int part_seq = 0;
     		for (int seq = 0; seq < union_taxa.size(); seq++)
     			if (union_taxa[seq] == 1) {
     				char ch = aln->STATE_UNKNOWN;
-    				if (taxa_set[seq] == 1) {
-    					ch = (*it)[part_seq++];
-    				}
+                    int seq_part = taxa_index[seq][id];
+                    if (seq_part >= 0)
+                        ch = (*it)[seq_part];
+                    //if (taxa_set[seq] == 1) {
+                    //    ch = (*it)[part_seq++];
+                    //}
     				pat.push_back(ch);
     			}
-    		ASSERT(part_seq == partitions[id]->getNSeq());
+    		//ASSERT(part_seq == partitions[id]->getNSeq());
     		aln->addPattern(pat, site, (*it).frequency);
     		// IMPORTANT BUG FIX FOLLOW
     		int ptnindex = aln->pattern_index[pat];

debian/changelog

+iqtree (1.6.7.1+dfsg-1) unstable; urgency=medium
+
+  * New upstream version
+  * Standards-Version: 4.2.1
+  * Delete RPATH from /usr/bin/iqtree-omp
+  * Respect DEB_BUILD_OPTIONS in override_dh_auto_test
+
+ -- Andreas Tille <tille@debian.org>  Wed, 26 Sep 2018 19:59:01 +0200
+
 iqtree (1.6.5+dfsg-1) unstable; urgency=medium
 
   * New upstream version

debian/control

@@ -13,7 +13,7 @@ Build-Depends: debhelper (>= 11~),
                help2man,
                time,
                chrpath
-Standards-Version: 4.1.4
+Standards-Version: 4.2.1
 Vcs-Browser: https://salsa.debian.org/med-team/iqtree
 Vcs-Git: https://salsa.debian.org/med-team/iqtree.git
 Homepage: http://www.cibiv.at/software/iqtree/

debian/rules

@@ -38,6 +38,7 @@ dh_auto_install_%:
 override_dh_install:
 	dh_install
 	find debian -name iqtree-mpi -exec chrpath --delete \{\} \;
+	find debian -name iqtree-omp -exec chrpath --delete \{\} \;
 
 dh_auto_clean_%:
 	dh_auto_clean -Bbuild.$(subst dh_auto_clean_,,$@)
@@ -63,6 +64,7 @@ override_dh_installman:
 	    --version-string="$(DEB_VERSION_UPSTREAM)" $(CURDIR)/debian/$(DEB_SOURCE)/usr/bin/iqtree-mpi > $(mandir)/iqtree-mpi.1
 
 override_dh_auto_test:
+ifeq (,$(filter nocheck,$(DEB_BUILD_OPTIONS)))
 	# use only the first example for build time test to save time on autobuilders
 #	if [ "`find $(CURDIR) -name iqtree -type f -executable`" = "" ] ; then \
 #		iqtreeomp=`find $(CURDIR) -name iqtree-omp -type f -executable` ; \
@@ -74,3 +76,4 @@ ifneq ($(shell nproc), 1)
 	time sh example.short
 	rm example.short
 endif
+endif

main/main.cpp

@@ -54,7 +54,7 @@
 #include "phyloanalysis.h"
 #include "tree/matree.h"
 //#include "ngs.h"
-//#include "parsmultistate.h"
+#include "obsolete/parsmultistate.h"
 //#include "gss.h"
 #include "alignment/maalignment.h" //added by MA
 #include "tree/ncbitree.h"
@@ -2540,8 +2540,8 @@ int main(int argc, char *argv[]) {
 	// call the main function
 	if (Params::getInstance().tree_gen != NONE) {
 		generateRandomTree(Params::getInstance());
-//	} else if (Params::getInstance().do_pars_multistate) {
-//		doParsMultiState(Params::getInstance());
+    } else if (Params::getInstance().do_pars_multistate) {
+        doParsMultiState(Params::getInstance());
 	} else if (Params::getInstance().rf_dist_mode != 0) {
 		computeRFDist(Params::getInstance());
 	} else if (Params::getInstance().test_input != TEST_NONE) {

main/phyloanalysis.cpp

@@ -634,12 +634,12 @@ void reportTree(ofstream &out, Params &params, PhyloTree &tree, double tree_lh,
         if (params.num_bootstrap_samples && params.compute_ml_tree) {
             if (params.aLRT_replicates > 0 || params.aLRT_test || params.aBayes_test)
                 out << " /";
-            out << " standard " << ((params.jackknife_prop == 0.0) ? "bootstrap" : "jackknife") << " support (%)";
+            out << " standard " << RESAMPLE_NAME << " support (%)";
         }
         if (params.gbo_replicates) {
             if (params.aLRT_replicates > 0 || params.aLRT_test || params.aBayes_test)
                 out << " /";
-            out << " ultrafast " << ((params.jackknife_prop==0.0) ? "bootstrap" : "jackknife") << " support (%)";
+            out << " ultrafast " << RESAMPLE_NAME << " support (%)";
         }
         out << endl;
     }
@@ -798,7 +798,7 @@ void printOutfilesInfo(Params &params, IQTree &tree) {
     }
 
 	if (params.gbo_replicates) {
-        cout << endl << "Ultrafast " << ((params.jackknife_prop==0.0) ? "bootstrap" : "jackknife") << " approximation results written to:" << endl
+        cout << endl << "Ultrafast " << RESAMPLE_NAME << " approximation results written to:" << endl
 			 << "  Split support values:          " << params.out_prefix << ".splits.nex" << endl
 			 << "  Consensus tree:                " << params.out_prefix << ".contree" << endl;
 		if (params.print_ufboot_trees)
@@ -880,11 +880,11 @@ void reportPhyloAnalysis(Params &params, IQTree &tree, ModelCheckpoint &model_in
         if (params.num_bootstrap_samples > 0) {
             if (params.compute_ml_tree)
                 out << " + ";
-            out << "non-parametric " << ((params.jackknife_prop==0.0) ? "bootstrap" : "jackknife") << " (" << params.num_bootstrap_samples
+            out << "non-parametric " << RESAMPLE_NAME << " (" << params.num_bootstrap_samples
                     << " replicates)";
         }
         if (params.gbo_replicates > 0) {
-            out << " + ultrafast " << ((params.jackknife_prop==0.0) ? "bootstrap" : "jackknife") << " (" << params.gbo_replicates << " replicates)";
+            out << " + ultrafast " << RESAMPLE_NAME << " (" << params.gbo_replicates << " replicates)";
         }
 		out << endl;
 		out << "Random seed number: " << params.ran_seed << endl << endl;
@@ -919,7 +919,7 @@ void reportPhyloAnalysis(Params &params, IQTree &tree, ModelCheckpoint &model_in
 				case SEQ_CODON: out << "CODON"; break;
 				case SEQ_DNA: out << "DNA"; break;
 				case SEQ_MORPH: out << "MORPH"; break;
-				case SEQ_MULTISTATE: out << "TINA"; break;
+				case SEQ_MULTISTATE: out << "MULTI"; break;
 				case SEQ_PROTEIN: out << "AA"; break;
 				case SEQ_POMO: out << "POMO"; break;
 				case SEQ_UNKNOWN: out << "???"; break;
@@ -1118,7 +1118,7 @@ void reportPhyloAnalysis(Params &params, IQTree &tree, ModelCheckpoint &model_in
 			out << "CONSENSUS TREE" << endl << "--------------" << endl << endl;
 			out << "Consensus tree is constructed from "
 					<< (params.num_bootstrap_samples ? params.num_bootstrap_samples : params.gbo_replicates)
-					<< ((params.jackknife_prop == 0.0) ? " bootstrap" : " jackknife") << " trees";
+					<< RESAMPLE_NAME << " trees";
             if (params.gbo_replicates || params.num_bootstrap_samples) {
                 out << endl << "Log-likelihood of consensus tree: " << fixed << tree.boot_consense_logl;
             }
@@ -1142,7 +1142,7 @@ void reportPhyloAnalysis(Params &params, IQTree &tree, ModelCheckpoint &model_in
 				out << " (strict consensus)";
 
 			out << endl << "Branch lengths are optimized by maximum likelihood on original alignment" << endl;
-			out << "Numbers in parentheses are " << ((params.jackknife_prop == 0.0) ? "bootstrap" : "jackknife") << " supports (%)" << endl << endl;
+			out << "Numbers in parentheses are " << RESAMPLE_NAME << " supports (%)" << endl << endl;
 
 			bool rooted = false;
 			MTree contree;
@@ -2431,7 +2431,7 @@ void runTreeReconstruction(Params &params, IQTree* &iqtree) {
         if (params.print_ufboot_trees)
             iqtree->writeUFBootTrees(params);
         
-        cout << endl << "Computing " << ((params.jackknife_prop == 0.0) ? "bootstrap" : "jackknife") << " consensus tree..." << endl;
+        cout << endl << "Computing " << RESAMPLE_NAME << " consensus tree..." << endl;
         string splitsfile = params.out_prefix;
         splitsfile += ".splits.nex";
         double weight_threshold = (params.split_threshold<1) ? params.split_threshold : (params.gbo_replicates-1.0)/params.gbo_replicates;
@@ -2616,7 +2616,7 @@ void runMultipleTreeReconstruction(Params &params, Alignment *alignment, IQTree
     // initialize tree and model strucgture
     ModelsBlock *models_block = readModelsDefinition(params);
     tree->setParams(&params);
-    tree->num_threads = params.num_threads;
+    tree->setNumThreads(params.num_threads);
     if (!tree->getModelFactory()) {
         tree->initializeModel(params, tree->aln->model_name, models_block);
     }
@@ -2908,9 +2908,12 @@ void runStandardBootstrap(Params &params, Alignment *alignment, IQTree *tree) {
 
     startTreeReconstruction(params, tree, *model_info);
     
+    // 2018-06-21: bug fix: alignment might be changed by -m ...MERGE
+    alignment = tree->aln;
+    
 	// do bootstrap analysis
 	for (int sample = bootSample; sample < params.num_bootstrap_samples; sample++) {
-		cout << endl << "===> START " << ((params.jackknife_prop==0.0) ? "BOOTSTRAP" : "JACKKNIFE") << " REPLICATE NUMBER "
+		cout << endl << "===> START " << RESAMPLE_NAME_UPPER << " REPLICATE NUMBER "
 				<< sample + 1 << endl << endl;
 
         // 2015-12-17: initialize random stream for creating bootstrap samples
@@ -2919,10 +2922,7 @@ void runStandardBootstrap(Params &params, Alignment *alignment, IQTree *tree) {
         init_random(params.ran_seed + sample);
 
 		Alignment* bootstrap_alignment;
-        if (Params::getInstance().jackknife_prop == 0.0)
-            cout << "Creating bootstrap alignment (seed: " << params.ran_seed+sample << ")..." << endl;
-        else
-            cout << "Creating jackknife alignment (seed: " << params.ran_seed+sample << ")..." << endl;
+        cout << "Creating " << RESAMPLE_NAME << " alignment (seed: " << params.ran_seed+sample << ")..." << endl;
 
 		if (alignment->isSuperAlignment())
 			bootstrap_alignment = new SuperAlignment;
@@ -3034,7 +3034,7 @@ void runStandardBootstrap(Params &params, Alignment *alignment, IQTree *tree) {
 	if (params.consensus_type == CT_CONSENSUS_TREE && MPIHelper::getInstance().isMaster()) {
 
 		cout << endl << "===> COMPUTE CONSENSUS TREE FROM " << params.num_bootstrap_samples
-            << ((params.jackknife_prop==0.0) ? " BOOTSTRAP" : " JACKKNIFE") << " TREES" << endl << endl;
+            << RESAMPLE_NAME_UPPER << " TREES" << endl << endl;
         string root_name = (params.root) ? params.root : alignment->getSeqName(0);
         const char* saved_root = params.root;
         params.root = root_name.c_str();
@@ -3062,7 +3062,7 @@ void runStandardBootstrap(Params &params, Alignment *alignment, IQTree *tree) {
                 optimizeConTree(params, tree);
             }
 
-            cout << endl << "===> ASSIGN " << ((params.jackknife_prop==0.0) ? " BOOTSTRAP" : " JACKKNIFE")
+            cout << endl << "===> ASSIGN " << RESAMPLE_NAME_UPPER
                 << " SUPPORTS TO THE TREE FROM ORIGINAL ALIGNMENT" << endl << endl;
             MExtTree ext_tree;
             assignBootstrapSupport(boottrees_name.c_str(), 0, 1e6,
@@ -3086,12 +3086,12 @@ void runStandardBootstrap(Params &params, Alignment *alignment, IQTree *tree) {
 		cout << endl;
 
     if (MPIHelper::getInstance().isMaster()) {
-        cout << "Total CPU time for " << ((params.jackknife_prop == 0.0) ? "bootstrap" : "jackknife") << ": " << (getCPUTime() - start_time) << " seconds." << endl;
-	cout << "Total wall-clock time for " << ((params.jackknife_prop == 0.0) ? "bootstrap" : "jackknife") << ": " << (getRealTime() - start_real_time) << " seconds." << endl << endl;
-	cout << "Non-parametric " << ((params.jackknife_prop == 0.0) ? "bootstrap" : "jackknife") << " results written to:" << endl;
+        cout << "Total CPU time for " << RESAMPLE_NAME << ": " << (getCPUTime() - start_time) << " seconds." << endl;
+	cout << "Total wall-clock time for " << RESAMPLE_NAME << ": " << (getRealTime() - start_real_time) << " seconds." << endl << endl;
+	cout << "Non-parametric " << RESAMPLE_NAME << " results written to:" << endl;
 	if (params.print_bootaln)
-		cout << ((params.jackknife_prop == 0.0) ? "  Bootstrap" : "  Jackknife") << " alignments:     " << params.out_prefix << ".bootaln" << endl;
-	cout << ((params.jackknife_prop == 0.0) ? "  Bootstrap" : "  Jackknife") << " trees:          " << params.out_prefix << ".boottrees" << endl;
+		cout << RESAMPLE_NAME_I << " alignments:     " << params.out_prefix << ".bootaln" << endl;
+	cout << RESAMPLE_NAME_I << " trees:          " << params.out_prefix << ".boottrees" << endl;
 	if (params.consensus_type == CT_CONSENSUS_TREE)
 		cout << "  Consensus tree:           " << params.out_prefix << ".contree" << endl;
 	cout << endl;
@@ -3104,10 +3104,7 @@ void convertAlignment(Params &params, IQTree *iqtree) {
 	if (params.num_bootstrap_samples || params.print_bootaln) {
 		// create bootstrap alignment
 		Alignment* bootstrap_alignment;
-        if (Params::getInstance().jackknife_prop == 0.0)
-            cout << "Creating bootstrap alignment..." << endl;
-        else
-            cout << "Creating jackknife alignment..." << endl;
+        cout << "Creating " << RESAMPLE_NAME << " alignment..." << endl;
 		if (alignment->isSuperAlignment())
 			bootstrap_alignment = new SuperAlignment;
 		else
@@ -3158,7 +3155,8 @@ void computeSiteFrequencyModel(Params &params, Alignment *alignment) {
     delete models_block;
     tree->setModel(tree->getModelFactory()->model);
     tree->setRate(tree->getModelFactory()->site_rate);
-    tree->setLikelihoodKernel(params.SSE, params.num_threads);
+    tree->setLikelihoodKernel(params.SSE);
+    tree->setNumThreads(params.num_threads);
     
     if (!tree->getModel()->isMixture())
         outError("No mixture model was specified!");
@@ -3526,7 +3524,7 @@ void assignBootstrapSupport(const char *input_trees, int burnin, int max_count,
 	// compute the percentage of appearance
 	//	printSplitSet(sg, hash_ss);
 	//sg.report(cout);
-    cout << "Creating " << ((params->jackknife_prop == 0.0) ? "bootstrap" : "jackknife") << " support values..." << endl;
+    cout << "Creating " << RESAMPLE_NAME << " support values..." << endl;
 	mytree.createBootstrapSupport(taxname, boot_trees, sg, hash_ss, params->support_tag);
 	//mytree.scaleLength(100.0/boot_trees.size(), true);
 	string out_file;

main/phylotesting.cpp

@@ -49,7 +49,9 @@ const char* bin_model_names[] = { "JC2", "GTR2" };
 
 
 /******* Morphological model set ******/
-const char* morph_model_names[] = {"MK", "ORDERED"};
+// 2018-08-20: don't test ORDERED model due to lots of numerical issues
+//const char* morph_model_names[] = {"MK", "ORDERED"};
+const char* morph_model_names[] = {"MK"};
 
 
 /******* DNA model set ******/
@@ -973,6 +975,7 @@ int getModelList(Params &params, Alignment *aln, StrVector &models, bool separat
 	bool test_options_noASC_I_new[] = {true,false,  false, true,  false,    false, true,    false};
 	bool test_options_asc_new[]   ={false,false,   true,false,  false,     true,false,     true};
 	bool test_options_pomo[]      = {true, false,  false, true,   false,   false,false,    false};
+    bool test_options_norate[]    = {true, false,  false, false,   false,   false,false,    false};
     bool *test_options = test_options_default;
 //	bool test_options_codon[] =  {true,false,  false,false,  false,    false};
 	const int noptions = sizeof(rate_options) / sizeof(char*);
@@ -1163,6 +1166,9 @@ int getModelList(Params &params, Alignment *aln, StrVector &models, bool separat
             test_options = test_options_morph;
         else
             test_options = test_options_noASC_I;
+    } else if (aln->frac_invariant_sites >= 1.0) {
+        // 2018-06-12: alignment with only invariant sites, no rate variation added
+        test_options = test_options_norate;
 	} else {
         // normal data, use +I instead
         if (with_new) {
@@ -1764,7 +1770,8 @@ void testPartitionModel(Params &params, PhyloSuperTree* in_tree, ModelCheckpoint
 		best_model.name = testModel(params, this_tree, part_model_info, models_block,
             (parallel_over_partitions ? 1 : num_threads), params.partition_type, in_tree->at(i)->aln->name, false, part_model_name);
 
-        ASSERT(best_model.restoreCheckpoint(&part_model_info));
+        bool check = (best_model.restoreCheckpoint(&part_model_info));
+        ASSERT(check);
 
 		double score = best_model.computeICScore(this_tree->getAlnNSite());
 		in_tree->at(i)->aln->model_name = best_model.name;
@@ -3135,6 +3142,10 @@ void evaluateTrees(Params &params, IQTree *tree, vector<TreeInfo> &info, IntVect
 		}
 		tree->freeNode();
 		tree->readTree(in, tree->rooted);
+        if (!tree->findNodeName(tree->aln->getSeqName(0))) {
+            outError("Taxon " + tree->aln->getSeqName(0) + " not found in tree");
+        }
+
         if (tree->rooted && tree->getModel()->isReversible()) {
             if (tree->leafNum != tree->aln->getNSeq()+1)
                 outError("Tree does not have same number of taxa as alignment");

model/modelfactory.cpp

@@ -1176,7 +1176,17 @@ double ModelFactory::optimizeParameters(int fixed_len, bool write_info,
     cur_lh = tree->computeLikelihood();
     tree->setCurScore(cur_lh);
     if (verbose_mode >= VB_MED || write_info) {
+	int p = -1;
+
+	// SET precision to 17 (temporarily)
+	if (verbose_mode >= VB_DEBUG) p = cout.precision(17);
+
+	// PRINT Log-Likelihood
 	cout << "1. Initial log-likelihood: " << cur_lh << endl;
+
+	// RESTORE previous precision
+	if (verbose_mode >= VB_DEBUG) cout.precision(p);
+
         if (verbose_mode >= VB_MAX) {
             tree->printTree(cout);
             cout << endl;

model/modelliemarkov.cpp

@@ -375,12 +375,13 @@ void ModelLieMarkov::init(const char *model_name, string model_params, StateFreq
 // we could make virtual setRates in ModelMarkov and
 // perhaps move this code all into there. - MDW
 void ModelLieMarkov::startCheckpoint() {
-    checkpoint->startStruct("ModelLieMarkov");
+    checkpoint->startStruct("ModelLieMarkov" + name);
 }
 
 void ModelLieMarkov::saveCheckpoint() {
     // saves model_parameters
     startCheckpoint();
+    if (num_params > 0)
         CKP_ARRAY_SAVE(num_params, model_parameters);
     endCheckpoint();
     ModelMarkov::saveCheckpoint();
@@ -390,6 +391,7 @@ void ModelLieMarkov::restoreCheckpoint() {
     ModelMarkov::restoreCheckpoint();
     // restores model_parameters
     startCheckpoint();
+    if (num_params > 0)
         CKP_ARRAY_RESTORE(num_params, model_parameters);
     endCheckpoint();
     setRates();                        // updates rate matrix

model/modelmixture.cpp

@@ -1179,6 +1179,8 @@ void ModelMixture::initMixture(string orig_model_name, string model_name, string
 	} else {
 		if (freq_params != "")
 			readStateFreq(freq_params);
+        if (freq == FREQ_UNKNOWN)
+            freq = FREQ_USER_DEFINED;
 		ModelMarkov::init(freq);
 	}
 
@@ -1402,14 +1404,16 @@ void ModelMixture::setCheckpoint(Checkpoint *checkpoint) {
 }
 
 void ModelMixture::startCheckpoint() {
-    checkpoint->startStruct("ModelMixture");
+    checkpoint->startStruct("ModelMixture" + convertIntToString(getNMixtures()));
 }
 
 void ModelMixture::saveCheckpoint() {
     startCheckpoint();
 //    CKP_SAVE(fix_prop);
+    if (!fix_prop) {
         int nmix = getNMixtures();
         CKP_ARRAY_SAVE(nmix, prop);
+    }
     int part = 1;
     for (iterator it = begin(); it != end(); it++, part++) {
         checkpoint->startStruct("Component" + convertIntToString(part));
@@ -1426,8 +1430,10 @@ void ModelMixture::restoreCheckpoint() {
 
     startCheckpoint();
 //    CKP_RESTORE(fix_prop);
+    if (!fix_prop) {
         int nmix = getNMixtures();
         CKP_ARRAY_RESTORE(nmix, prop);
+    }
     int part = 1;
     for (iterator it = begin(); it != end(); it++, part++) {
         checkpoint->startStruct("Component" + convertIntToString(part));
@@ -1627,7 +1633,9 @@ double ModelMixture::optimizeWithEM(double gradient_epsilon) {
     tree->copyPhyloTree(phylo_tree);
     tree->optimize_by_newton = phylo_tree->optimize_by_newton;
     tree->setParams(phylo_tree->params);
-    tree->setLikelihoodKernel(phylo_tree->sse, phylo_tree->num_threads);
+    tree->setLikelihoodKernel(phylo_tree->sse);
+    tree->setNumThreads(phylo_tree->num_threads);
+    
     // initialize model
     ModelFactory *model_fac = new ModelFactory();
     model_fac->joint_optimize = phylo_tree->params->optimize_model_rate_joint;

model/ratefree.cpp

@@ -128,13 +128,13 @@ void RateFree::initFromCatMinusOne() {
     restoreCheckpoint();
     ncategory++;
 
-    int first = 0, second = -1, i;
+    int first = 0, i;
     // get the category k with largest proportion
     for (i = 1; i < ncategory-1; i++)
         if (prop[i] > prop[first]) {
             first = i;
         }
-    second = (first == 0) ? 1 : 0;
+    int second = (first == 0) ? 1 : 0;
     for (i = 0; i < ncategory-1; i++)
         if (prop[i] > prop[second] && second != first)
             second = i;
@@ -143,9 +143,15 @@ void RateFree::initFromCatMinusOne() {
 //    memmove(prop, input->prop, (k+1)*sizeof(double));
 
     // divide highest category into 2 of the same prop
+    // 2018-06-12: fix bug negative rates
+    if (-rates[second] + 3*rates[first] > 0.0) {
         rates[ncategory-1] = (-rates[second] + 3*rates[first])/2.0;
-    prop[ncategory-1] = prop[first]/2;
         rates[first] = (rates[second]+rates[first])/2.0;
+    } else {
+        rates[ncategory-1] = (3*rates[first])/2.0;
+        rates[first] = (rates[first])/2.0;
+    }
+    prop[ncategory-1] = prop[first]/2;
     prop[first] = prop[first]/2;
 //    if (k < ncategory-2) {
 //        memcpy(&rates[k+2], &input->rates[k+1], (ncategory-2-k)*sizeof(double));
@@ -496,7 +502,9 @@ double RateFree::optimizeWithEM() {
     tree->copyPhyloTree(phylo_tree);
     tree->optimize_by_newton = phylo_tree->optimize_by_newton;
     tree->setParams(phylo_tree->params);
-    tree->setLikelihoodKernel(phylo_tree->sse, phylo_tree->num_threads);
+    tree->setLikelihoodKernel(phylo_tree->sse);
+    tree->setNumThreads(phylo_tree->num_threads);
+
     // initialize model
     ModelFactory *model_fac = new ModelFactory();
     model_fac->joint_optimize = phylo_tree->params->optimize_model_rate_joint;
@@ -521,8 +529,15 @@ double RateFree::optimizeWithEM() {
         }
         ASSERT(score < 0);
         
-        if (step > 0)
+        if (step > 0) {
+            if (score <= old_score-0.1) {
+                phylo_tree->printTree(cout, WT_BR_LEN+WT_NEWLINE);
+                writeInfo(cout);
+                cout << "Partition " << phylo_tree->aln->name << endl;
+                cout << "score: " << score << "  old_score: " << old_score << endl;
+            }
             ASSERT(score > old_score-0.1);
+        }
             
         old_score = score;
         
@@ -604,7 +619,7 @@ double RateFree::optimizeWithEM() {
             subst_model->setTree(tree);
             model_fac->model = subst_model;
             if (subst_model->isMixture() || subst_model->isSiteSpecificModel() || !subst_model->isReversible())
-                tree->setLikelihoodKernel(phylo_tree->sse, phylo_tree->num_threads);
+                tree->setLikelihoodKernel(phylo_tree->sse);
 
                         
             // initialize likelihood

ncl/nxscharactersblock.cpp

@@ -2681,12 +2681,6 @@ void NxsCharactersBlock::Reset()
 	missing				= '?';
 	gap					= '\0';
 	matchchar			= '\0';
-	matrix				= NULL;
-	charPos				= NULL;
-	taxonPos			= NULL;
-	activeTaxon			= NULL;
-	activeChar			= NULL;
-	symbols				= NULL;
 
 	ResetSymbols();
 

ncl/nxscharactersblock.h

@@ -215,6 +215,7 @@ class NxsCharactersBlock
 		char					GetGapSymbol();
 		char					GetMatchcharSymbol();
 		char					GetMissingSymbol();
+        NxsDiscreteMatrix       *GetMatrix() { return matrix; }
 		bool					IsGapState(unsigned i, unsigned j);
 		bool					IsInterleave();
 		bool					IsLabels();

nclextra/msetsblock.cpp

@@ -161,6 +161,7 @@ void MSetsBlock::Read(NxsToken &token)
 				myset->aln_file = myset->position_spec.substr(0, pos);
 				myset->position_spec = myset->position_spec.substr(pos+1);
 			}
+            trimString(myset->position_spec);
             if ((pos=myset->position_spec.find(',')) != string::npos && isalpha(myset->position_spec[0])) {
                 myset->sequence_type = myset->position_spec.substr(0, pos);
                 myset->position_spec = myset->position_spec.substr(pos+1);
@@ -174,7 +175,11 @@ void MSetsBlock::Read(NxsToken &token)
 				errormsg += " instead";
 				throw NxsException(errormsg, token.GetFilePosition(), token.GetFileLine(), token.GetFileColumn());
 			}
-
+            trimString(myset->aln_file);
+            trimString(myset->char_partition);
+            trimString(myset->model_name);
+            trimString(myset->position_spec);
+            trimString(myset->sequence_type);
 		} // if (token.Equals("CHARSET"))
 		else if (token.Equals("CHARPARTITION"))
 		{

obsolete/parsmultistate.cpp

@@ -22,15 +22,25 @@
 #include "tree/tinatree.h"
 #include "parsmultistate.h"
 #include "alignment/alignment.h"
+#include "tree/parstree.h"
 
 void doParsMultiState(Params &params) {
-	cout << "Here\n";
-    Alignment alignment(params.aln_file, params.sequence_type, params.intype);
-    TinaTree tree;
-    tree.readTree(params.user_file, params.is_rooted);
-	tree.setAlignment(&alignment);
-	tree.drawTree(cout);
-	cout << "Parsimony score is: " << tree.computeParsimonyScore() << endl;
-	cout << "Parsimony score ver2 is: " << tree.computeParsimony() << endl;
+    Alignment alignment(params.aln_file, params.sequence_type, params.intype, "");
+    alignment.orderPatternByNumChars(PAT_VARIANT);
+    ParsTree pars_tree;
+    pars_tree.readTree(params.user_file, params.is_rooted);
+    if (pars_tree.rooted)
+        pars_tree.convertToUnrooted();
+    pars_tree.setAlignment(&alignment);
+    pars_tree.initCostMatrix(CM_LINEAR);
+    int total_length = pars_tree.computeParsimony();
+    cout << "total length: " << total_length << endl;
+    pars_tree.initCostMatrix(CM_UNIFORM);
+    int pars_score = pars_tree.computeParsimony();
+    cout.unsetf(ios::fixed);
+    cout.precision(6);
+    cout << "mean length: " << double(total_length)/pars_score << endl;
+    cout << "Parsimony score is: " << pars_score << endl;
+	//cout << "Parsimony score ver2 is: " << tree.computeParsimony() << endl;
 	//tree.printParsimonyStates();
 }

obsolete/parsmultistate.h

@@ -21,7 +21,7 @@
 #ifndef PARSMULTISTATE_H
 #define PARSMULTISTATE_H
 
-#include "tools.h"
+#include "utils/tools.h"
 
 void doParsMultiState(Params &params);
 

tree/CMakeLists.txt

@@ -39,6 +39,8 @@ supernode.cpp
 supernode.h
 tinatree.cpp
 tinatree.h
+parstree.cpp
+parstree.h
 )
 
 target_link_libraries(tree pll model alignment)