Copying over ("merging") recent changes from iqtree2's dev branch:

1. Added the files hashrow.h and fancyrapidnj.h.
2. Supporting changes, in other files, for the FancyRapidNJ class:
   (a) added ClusterTree::appendToLastCluster
   (b) advertising of BROKEN and BROKEN-V (temporary names
       for Fancy Rapid NJ)
3. Tweaked the omp parallel directives in BoundingMatrix,
   to squeeze a bit more speed out of BoundingMatrix (which
   implements NJ-R and BIONJ-R).
4. Added VectorizedBoundingMatrix
5. Supporting changes, in other files:
   (a) advertising of NJ-R-V (VectorizedBoundingMatrix)
6. Moved calculation of RMS of T (the distance matrix implied
   by the constructed tree) minus D (the input distance matrix)
   out of upgma.h, and into hashrow.h (the hashing!), and
   clustertree.h (calculateDistancesToLeaves and
   calculateRMSOfTMinusD).
7. Added DoubleArgument to utils/argument.h and utils/argument.cpp.
8. Fixed up definition of the "-dist" StringArgument (so it
   once again works), and added a "-max-dist" DoubleArgument
   (it's no longer "hard-coded", only defaulted as 10.0).
   Plus, supporting changes in:
   (a) sequence.h and sequence.cpp (extra parameters to
       correctedDistance and uncorrectedDistance,
       to indicate max distance), likewise to
       SequenceLoader'S constructor.
8. If compression is requested, via "-c" and "-gz",
   distance output format has ".gz" appended to it (so it will
   be compressed).  Previously, newick tree files honoured "-c"
   and "-gz", but distance matrix files did not.
9. FlatMatrix::writeToDistanceFile now reports how much of
   the file it has written (via a progress_display instance).
10.Fix ups in pigzstream::open().  It now has a much better
   failure path (previously it claimed it had succeeded,
   in what it wrote to standard output!).  In utils/gzstream.cpp.
   This works via a new member function, markAsFailed(),
   in progress_display (in utils/progress.h and .cpp).
   If it's been called, "done" is replaced with "failed" in
   what is written to standard output when done() is called.
11.MAJOR bug fix in utils/hammingdistance.h
   (the vectorized version of countBitsSetInEitherTemplate()
    didn't work properly at all) (W parameter is the size of
    V in uint64_t instances, not its size in bytes, so the
    j loop is (j=0; j<W; ++j)... not (j=0; j<W/8; ++j)...)

Author: James Barbetti

auctionmatrix.h

@@ -30,10 +30,10 @@ namespace StartTree {
         using super::entryToCluster;
         using super::entriesSorted;
         using super::purgeRow;
-        virtual std::string getAlgorithmName() const {
+        virtual std::string getAlgorithmName() const override {
             return "Auction" + M::getAlgorithmName();
         }
-        StartTree::Position<T> getRowMinimum(size_t row, T maxTot, T qBest) const {
+        virtual StartTree::Position<T> getRowMinimum(intptr_t row, T maxTot, T qBest) const override {
             StartTree::Position<T> pos(row, 0, infiniteDistance, 0);
             const T*   rowData   = entriesSorted.rows[row];
             const int* toCluster = entryToCluster.rows[row];

bionj2.cpp

@@ -55,6 +55,7 @@
 #include "nj.h"
 #include "rapidnj.h"
 #include "auctionmatrix.h"
+#include "fancyrapidnj.h"
 
 namespace StartTree
 {
@@ -64,10 +65,22 @@ namespace StartTree
 
 void addBioNJ2020TreeBuilders(Factory& f) {
     const char* defaultName = "RapidNJ";
-    ADVERTISE(NJMatrix<NJFloat>,    "NJ",      "Neighbour Joining (Saitou, Nei [1987])");
-    ADVERTISE(UNJMatrix<NJFloat>,   "UNJ",     "Unweighted Neighbour Joining (Gascel [1997])");
-    ADVERTISE(RapidNJ,              "NJ-R",    "Rapid Neighbour Joining"
-                                               " (Simonsen, Mailund, Pedersen [2011])");
+    ADVERTISE(NJMatrix<NJFloat>,     "NJ",      "Neighbour Joining (Saitou, Nei [1987])");
+    ADVERTISE(UNJMatrix<NJFloat>,    "UNJ",     "Unweighted Neighbour Joining (Gascel [1997])");
+    ADVERTISE(RapidNJ,               "NJ-R",    "Rapid Neighbour Joining"
+                                                " (Simonsen, Mailund, Pedersen [2011])");
+    ADVERTISE(Vectorized_RapidNJ,    "NJ-R-V",  "Rapid Neighbour Joining (Vectorized)"
+                                                " (Simonsen, Mailund, Pedersen [2011])");
+
+    ADVERTISE(RapidNJ,            defaultName,  "Rapid Neighbour Joining"
+                                                " (Simonsen, Mailund, Pedersen [2011])");
+
+    ADVERTISE(FancyNJMatrix<NJFloat>,"BROKEN",  "Rapid Neighbour Joining (Broken Version)"
+                                                " (Simonsen, Mailund, Pedersen [2011])");
+
+    ADVERTISE(VectorizedFancyNJMatrix<NJFloat>, "BROKEN-V",  "Rapid Neighbour Joining (Broken Version)"
+                                                " (Simonsen, Mailund, Pedersen [2011]) (Vectorized)");
+
 #ifdef USE_VECTORCLASS_LIBRARY
     ADVERTISE(VectorNJ,             "NJ-V",    "Vectorized Neighbour Joining (Saitou, Nei [1987])");
 #endif
@@ -82,7 +95,6 @@ void addBioNJ2020TreeBuilders(Factory& f) {
     ADVERTISE(VectorizedUPGMA_Matrix<NJFloat>, "UPGMA-V", "Vectorized UPGMA (Sokal, Michener [1958])");
 #endif
     ADVERTISE(BoundingMatrix<double>,"NJ-R-D", "Double precision Rapid Neighbour Joining");
-    ADVERTISE(RapidNJ,           defaultName,  "Rapid Neighbour Joining (Simonsen, Mailund, Pedersen [2011]) (default)");
     f.setNameOfDefaultTreeBuilder(defaultName);
     ADVERTISE(DistanceAuctionMatrix, "AUCTION",    "Auction Joining");
 }

clustertree.h

@@ -33,6 +33,7 @@
 #include <iostream>                  //for std::istream
 #include <sstream>                   //for std::stringstream
 #include <utils/progress.h>          //for progress_display
+#include <utils/parallel_mergesort.h>//for MergeSorter
 
 template <class T=double> struct Link {
     //
@@ -114,6 +115,108 @@ template <class T> class ClusterTree: public std::vector<Cluster<T>>
         //<< " is (" << a << ", " << b << ", " << c << ")" << std::endl;
         return cluster;
     }
+    Cluster<T>& appendToLastCluster
+    ( size_t c, T length) {
+        back().links.emplace_back(c, length);
+        back().countOfExteriorNodes += at(c).countOfExteriorNodes;
+        return back();
+    }
+    typedef std::pair<intptr_t /*taxon*/, T /*distance*/ >     LeafDistance;
+    typedef std::vector<LeafDistance> LeafDistanceVector;
+    void calculateDistancesToLeaves(intptr_t top, double branch_length,
+                                    LeafDistanceVector& ldv) {
+        //Running time proportional to subtree size
+        std::vector<LeafDistance> stack;
+        stack.emplace_back(top, branch_length);
+        while (!stack.empty()) {
+            LeafDistance b = stack.back();
+            stack.pop_back();
+            Cluster<T>& node = at(b.first);
+            if (node.links.empty()) {
+                ldv.emplace_back(b);
+                continue;
+            } 
+            for (Link<T>& link : node.links ) {
+                stack.emplace_back(link.clusterIndex, 
+                                   b.second + link.linkDistance );
+            }
+        }
+        //
+        //Hack: Reordering makes memory accesses "cache-friendlier"
+        //      for large clusters, in calculateRMSOfTMinusD().
+        //Note: Single-threaded sorting is used, here, because 
+        //      calculateRMSOfTMinusD() is already executing 
+        //      calls to this function in parallel.
+        //
+        MergeSorter<LeafDistance> sorter;
+        sorter.single_thread_sort(ldv.data(), ldv.size());
+    }
+    bool calculateRMSOfTMinusD(const double* matrix, intptr_t rank, double& rms) {
+        //Assumes: rank is at least 3.
+        //
+        //Total running time: proportional to rank*(rank-1)/2.
+        //(that's the number of additions to sum_of_squares).
+        //
+        //Total memory consumption:  depends on thread count, X.
+        //In theory, on the order of (X+1)*rank*(sizeof(intptr_t)+sizeof(T)).
+        //2 * X * (maximum sum of sizes of leaf distance vectors).
+        //The 2 is because of the use of a MergeSorter in 
+        //calculateDistancesToLeaves().  An in-place sort would 
+        //lower that to a 1.
+        //
+        intptr_t cluster_count  = size();
+        double   sum_of_squares = 0.0;
+        #ifdef _OPENMP
+        #pragma omp parallel for reduction(+:sum_of_squares)
+        #endif
+        for (intptr_t h=rank; h<cluster_count; ++h) {
+            //For each (non-leaf) cluster...
+            //1. calculate distances to all the 
+            //   leaves, for each contributing cluster.
+            Cluster<T>& node = at(h);
+            std::vector<LeafDistanceVector> subtrees;
+            for (Link<T>& link : node.links ) {
+                LeafDistanceVector ldv;
+                calculateDistancesToLeaves(link.clusterIndex, 
+                                           link.linkDistance, ldv);
+                subtrees.emplace_back(ldv);
+            }
+            //2. for each pair of LeafDistanceVectors, A, B,
+            //   (for separate contributing clusters)...
+            size_t subtree_count = subtrees.size();
+            for (size_t i=0; i+1<subtree_count; ++i) {
+                for (size_t j=i+1; j<subtree_count; ++j) {
+                    //2b. for each leaf (a.first) from A
+                    for (LeafDistance& a: subtrees[i]) {
+                        auto row = matrix + a.first * rank;
+                        //2c. for each leaf (b.first) from B
+                        //    calculate error; difference between
+                        //    distance(a.first) + distance(b.first)
+                        //    and D[a.first * rank + b.first].
+                        for (LeafDistance& b: subtrees[j]) {
+                            double diff     = a.second + b.second
+                                            - row[b.first];
+                            sum_of_squares += (diff * diff);
+                            #if (0)
+                            std::cout << a.first << " " << b.first << " "
+                                      << a.second + b.second << " "
+                                      << row[b.first] << " "
+                                      << diff << "\n";
+                            #endif
+                        }
+                    }
+                }
+            }
+            //std::cout << "\n";
+        }
+        double double_rank = static_cast<double>(rank);
+        rms = sqrt( sum_of_squares * 2.0 / double_rank / (double_rank - 1.0) );
+        //std::cout << "rank " << (double_rank) << std::endl;
+        //std::cout << "sum " << sum_of_squares << ","
+        //          << " divisor " << double_rank*(double_rank-1)*0.5 << std::endl;
+        return true;
+    }
+
     template <class F> 
     bool writeTreeToFile(int precision, 
                          const std::string &treeFilePath, 

decenttree.cpp

@@ -53,18 +53,19 @@ namespace {
         return answer;
     }
 
-    bool correcting_distances      = true;
-    bool is_DNA                    = true;
-    bool numbered_names            = false;
-    bool filter_problem_sequences  = false;
-    char unknown_char              = 'N';
-    int  precision                 = 8;
-    int  compression_level         = 9;
+    bool   correcting_distances     = true;
+    bool   is_DNA                   = true;
+    bool   numbered_names           = false;
+    bool   filter_problem_sequences = false;
+    char   unknown_char             = 'N';
+    double max_distance             = 10.0;
+    int    precision                = 8;
+    int    compression_level        = 9;
     std::string msaOutputPath; //write .msa formatted version of .fasta input here
     std::string alphabet;      //defaults to ACGT
     std::string unknown_chars; //defaults to .~_-?N
-    std::string format             = "square.interleaved";
-    bool        interleaved_format = true;
+    std::string format              = "square.interleaved";
+    bool        interleaved_format  = true;
 
     std::string stripName;              //characters to strip from names
     std::string nameReplace("_");       //characters to replace stripepd chars with, in names
@@ -116,7 +117,8 @@ void showUsage() {
 bool loadSequenceDistancesIntoMatrix(Sequences& sequences,
                                      const std::vector<char>&   is_site_variant,
                                      bool report_progress, FlatMatrix& m) {
-    SequenceLoader loader(unknown_char, is_DNA, sequences, correcting_distances,
+    SequenceLoader loader(unknown_char, is_DNA, max_distance,
+                          sequences, correcting_distances,
                           precision, compression_level, format, 
                           is_site_variant, report_progress);
     bool success = loader.loadSequenceDistances(m);
@@ -247,15 +249,17 @@ bool prepInput(const std::string& fastaFilePath,
             if (loadMatrix) {
                 useNumberedNamesIfAskedTo(numbered_names, m);
                 return m.writeToDistanceFile(format, precision,
-                                 compression_level,
+                                 compression_level, reportProgress,
                                  distanceOutputFilePath );
             }
             else {
-                SequenceLoader loader(unknown_char, is_DNA, sequences, 
+                SequenceLoader loader(unknown_char, is_DNA, 
+                                      max_distance, sequences, 
                                       correcting_distances, precision, 
                                       compression_level, format,
                                       is_site_variant, reportProgress);
-                bool success = loader.writeDistanceMatrixToFile(numbered_names, distanceOutputFilePath);
+                bool success = loader.writeDistanceMatrixToFile(numbered_names, 
+                                                                distanceOutputFilePath);
                 return success;
             }
         }
@@ -265,7 +269,7 @@ bool prepInput(const std::string& fastaFilePath,
             fixUpSequenceNames(truncateName, stripName, nameReplace, m);
             if (!distanceOutputFilePath.empty()) {
                 return m.writeToDistanceFile(format, precision,
-                                             compression_level,
+                                             compression_level, reportProgress,
                                              distanceOutputFilePath );
             }
             return true;
@@ -327,7 +331,7 @@ class DecentTreeOptions {
         arg_map << new StringArgument("-fasta", "fasta file path",             fastaFilePath);
         arg_map << new StringArgument("-phylip", "phylip alignment file path", phylipFilePath);
         arg_map << new StringArgument("-in",    "distance matrix file path",   inputFilePath);
-        arg_map << new StringArgument("-dist",  "distance matrix file path",   inputFilePath);
+        arg_map << new StringArgument("-dist",  "distance matrix file path",   distanceOutputFilePath);
         arg_map << new IntArgument   ("-c",     "compression level between 1 and 9", 
                                     compression_level);
         arg_map << new IntArgument   ("-f",     "precision level between 4 and 15",
@@ -349,6 +353,7 @@ class DecentTreeOptions {
         arg_map << new SwitchArgument("-gz",          isOutputZipped,           true);
         arg_map << new SwitchArgument("-no-banner",   isBannerSuppressed,       true);
         arg_map << new SwitchArgument("-uncorrected", correcting_distances,     false);
+        arg_map << new DoubleArgument("-max-dist",    "maximum distance", max_distance);
         arg_map << new IntArgument   ("-nt", "thread count", threadCount);
         arg_map << new SwitchArgument("-q",           beSilent,                 true);
         arg_map << new SwitchArgument("-filter",      filter_problem_sequences, true);
@@ -385,9 +390,17 @@ class DecentTreeOptions {
                 break;
             }
         }
+        if (max_distance<=0) {
+            PROBLEM("Maximum distance (" << max_distance << ") too low");
+        }
         range_restrict(0, 9,  compression_level );
         range_restrict(1, 15, precision );
         format = string_to_lower(format);
+        if (isOutputZipped && 
+            format.find(".gz") == std::string::npos) {
+            //Ensure that distance file will be compressed
+            format += ".gz"; 
+        }
         if (isOutputToStandardOutput) {
             outputFilePath = "STDOUT";
         }
@@ -502,7 +515,7 @@ int obeyCommandLineOptions(DecentTreeOptions& options) {
     else if (succeeded && options.isMatrixToBeLoaded) {
         succeeded = algorithm->constructTreeInMemory(m.getSequenceNames(),
                                                      m.getDistanceMatrix(),
-                                                     options.outputFilePath);    
+                                                     options.outputFilePath);
     }
     else if (!options.inputFilePath.empty()) {
         succeeded = algorithm->constructTree(options.inputFilePath, 

decenttree_update.sh

@@ -14,6 +14,7 @@ cd ../iqtree2/distancematrixtree
 cp decenttree.cpp flatmatrix.cpp flatmatrix.h starttree.cpp starttree.h ../../decenttree
 cp bionj.cpp bionj2.cpp bionj2.h upgma.h nj.h rapidnj.h auctionmatrix.h ../../decenttree
 cp clustertree.h distancematrix.h stitchup.cpp sequence.cpp sequence.h  ../../decenttree
+cp hashrow.h fancyrapidnj.h ../../decenttree
 cd ../utils
 cp argument.cpp argument.h               ../../decenttree/utils
 cp gzstream.cpp gzstream.h heapsort.h    ../../decenttree/utils