diff -Nru fasttree-2.1.9/changelog fasttree-2.1.10/changelog --- fasttree-2.1.9/changelog 2016-04-20 20:13:57.000000000 +0000 +++ fasttree-2.1.10/changelog 2017-08-18 11:30:32.000000000 +0000 @@ -1,3 +1,14 @@ +Version 2.1.10: April 11, 2017 + + Fix a bug when using GTR models with huge alignments with over 2 + billion non-gap characters. The frequencies of A, C, G, and T were + being computed with a 32-bit (signed) counter. This led to + negative frequencies and eventually to a crash with + "FastTree.c:9769: tqli: Assertion `iter < 30' failed.". SetMLGtr() + now uses 64-bit counters. Also, more information about the + optimization of the GTR model is savd in the log file (if using + the -log option). To support this, gtr_opt_t now includes fpLog. + Version 2.1.9: March 29, 2016 Add the -lg option to use Le/Gascuel model for amino acid diff -Nru fasttree-2.1.9/debian/changelog fasttree-2.1.10/debian/changelog --- fasttree-2.1.9/debian/changelog 2017-06-29 11:46:37.000000000 +0000 +++ fasttree-2.1.10/debian/changelog 2017-09-09 17:00:00.000000000 +0000 @@ -1,3 +1,9 @@ +fasttree (2.1.10-1) unstable; urgency=medium + + * New upstream version + + -- Thorsten Alteholz Sat, 09 Sep 2017 19:00:00 +0200 + fasttree (2.1.9-2) unstable; urgency=medium * Team upload diff -Nru fasttree-2.1.9/fasttree.c fasttree-2.1.10/fasttree.c --- fasttree-2.1.9/fasttree.c 2016-04-20 20:13:54.000000000 +0000 +++ fasttree-2.1.10/fasttree.c 2017-08-18 11:30:30.000000000 +0000 @@ -343,7 +343,7 @@ #endif /* USE_SSE3 */ -#define FT_VERSION "2.1.9" +#define FT_VERSION "2.1.10" char *usage = " FastTree protein_alignment > tree\n" @@ -1380,6 +1380,7 @@ double freq[4]; double rates[6]; int iRate; /* which rate to set x from */ + FILE *fpLog; /* OPTIONAL WRITE */ } gtr_opt_t; /* Returns -log_likelihood for the tree with the given rates @@ -5848,11 +5849,14 @@ assert(nCodes==4); gtr_opt_t gtr; gtr.NJ = NJ; + gtr.fpLog = fpLog; if (freq_in != NULL) { for (i=0; i<4; i++) gtr.freq[i]=freq_in[i]; } else { - int n[4] = {1,1,1,1}; /* pseudocounts */ + /* n[] and sum were int in FastTree 2.1.9 and earlier -- this + caused gtr analyses to fail on analyses with >2e9 positions */ + long n[4] = {1,1,1,1}; /* pseudocounts */ for (i=0; inSeq; i++) { unsigned char *codes = NJ->profiles[i]->codes; int iPos; @@ -5860,7 +5864,7 @@ if (codes[iPos] < 4) n[codes[iPos]]++; } - int sum = n[0]+n[1]+n[2]+n[3]; + long sum = n[0]+n[1]+n[2]+n[3]; for (i=0; i<4; i++) gtr.freq[i] = n[i]/(double)sum; } @@ -5903,6 +5907,7 @@ } double GTRNegLogLk(double x, void *data) { + gtr_opt_t *gtr = (gtr_opt_t*)data; assert(nCodes == 4); assert(gtr->NJ != NULL); @@ -5915,6 +5920,12 @@ rates[i] = gtr->rates[i]; rates[gtr->iRate] = x; + FILE *fpLog = gtr->fpLog; + if (fpLog) + fprintf(fpLog, "GTR_Opt\tfreq %.5f %.5f %.5f %.5f rates %.5f %.5f %.5f %.5f %.5f %.5f\n", + gtr->freq[0], gtr->freq[1], gtr->freq[2], gtr->freq[3], + rates[0], rates[1], rates[2], rates[3], rates[4], rates[5]); + gtr->NJ->transmat = CreateGTR(rates, gtr->freq); RecomputeMLProfiles(/*IN/OUT*/gtr->NJ); double loglk = TreeLogLk(gtr->NJ, /*site_loglk*/NULL); @@ -5923,7 +5934,10 @@ /* Do not recompute profiles -- assume the caller will do that */ if (verbose > 2) fprintf(stderr, "GTR LogLk(%.5f %.5f %.5f %.5f %.5f %.5f) = %f\n", - rates[0], rates[1], rates[2], rates[3], rates[4], rates[5], loglk); + rates[0], rates[1], rates[2], rates[3], rates[4], rates[5], loglk); + if (fpLog) + fprintf(fpLog, "GTR_Opt\tGTR LogLk(%.5f %.5f %.5f %.5f %.5f %.5f) = %f\n", + rates[0], rates[1], rates[2], rates[3], rates[4], rates[5], loglk); return(-loglk); }