fastcluster_dm.cxx

Go to the documentation of this file.

/*
  fastcluster: Fast hierarchical clustering routines for R and Python
 
  Copyright © 2011 Daniel Müllner
  <http://danifold.net>
 
  This library implements various fast algorithms for hierarchical,
  agglomerative clustering methods:
 
  (1) Algorithms for the "stored matrix approach": the input is the array of
      pairwise dissimilarities.
 
      MST_linkage_core: single linkage clustering with the "minimum spanning
      tree algorithm (Rohlfs)
 
      NN_chain_core: nearest-neighbor-chain algorithm, suitable for single,
      complete, average, weighted and Ward linkage (Murtagh)
 
      generic_linkage: generic algorithm, suitable for all distance update
      formulas (Müllner)
 
  (2) Algorithms for the "stored data approach": the input are points in a
      vector space.
 
      MST_linkage_core_vector: single linkage clustering for vector data
 
      generic_linkage_vector: generic algorithm for vector data, suitable for
      the Ward, centroid and median methods.
 
      generic_linkage_vector_alternative: alternative scheme for updating the
      nearest neighbors. This method seems faster than "generic_linkage_vector"
      for the centroid and median methods but slower for the Ward method.
 
  All these implementation treat infinity values correctly. They throw an
  exception if a NaN distance value occurs.
*/
 
// Older versions of Microsoft Visual Studio do not have the fenv header.
#ifdef _MSC_VER
#if (_MSC_VER == 1500 || _MSC_VER == 1600)
#define NO_INCLUDE_FENV
#endif
#endif
// NaN detection via fenv might not work on systems with software
// floating-point emulation (bug report for Debian armel).
#ifdef __SOFTFP__
#define NO_INCLUDE_FENV
#endif
#ifdef NO_INCLUDE_FENV
#pragma message("Do not use fenv header.")
#else
#pragma message("Use fenv header. If there is a warning about unknown #pragma STDC FENV_ACCESS, this can be ignored.")
#pragma STDC FENV_ACCESS on
#include <fenv.h>
#endif
 
#include <algorithm> // for std::fill_n
#include <cfloat>    // also for DBL_MAX, DBL_MIN
#include <cmath>     // for std::pow, std::sqrt
#include <cstddef>   // for std::ptrdiff_t
#include <limits>    // for std::numeric_limits<...>::infinity()
#include <stdexcept> // for std::runtime_error
#include <string>    // for std::string
#ifndef DBL_MANT_DIG
#error The constant DBL_MANT_DIG could not be defined.
#endif
#define T_FLOAt_MANT_DIG DBL_MANT_DIG
 
#ifndef LONG_MAX
#include <climits>
#endif
#ifndef LONG_MAX
#error The constant LONG_MAX could not be defined.
#endif
#ifndef INT_MAX
#error The constant INT_MAX could not be defined.
#endif
 
#ifndef INT32_MAX
#ifdef _MSC_VER
#if _MSC_VER >= 1600
#define __STDC_LIMIT_MACROS
#include <stdint.h>
#else
typedef __int32 int_fast32_t;
typedef __int64 int64_t;
#endif
#else
#define __STDC_LIMIT_MACROS
#include <stdint.h>
#endif
#endif
 
#define FILL_N std::fill_n
#ifdef _MSC_VER
#if _MSC_VER < 1600
#undef FILL_N
#define FILL_N stdext::unchecked_fill_n
#endif
#endif
 
// Suppress warnings about (potentially) uninitialized variables.
#ifdef _MSC_VER
#pragma warning(disable : 4700)
#endif
 
#ifndef HAVE_DIAGNOSTIC
#if __GNUC__ > 4 || (__GNUC__ == 4 && (__GNUC_MINOR__ >= 6))
#define HAVE_DIAGNOSTIC 1
#endif
#endif
 
#ifndef HAVE_VISIBILITY
#if __GNUC__ >= 4
#define HAVE_VISIBILITY 1
#endif
#endif
 
/* Since the public interface is given by the Python respectively R interface,
 * we do not want other symbols than the interface initalization routines to be
 * visible in the shared object file. The "visibility" switch is a GCC concept.
 * Hiding symbols keeps the relocation table small and decreases startup time.
 * See http://gcc.gnu.org/wiki/Visibility
 */
#if HAVE_VISIBILITY
#pragma GCC visibility push(hidden)
#endif
 
typedef int_fast32_t t_index;
#ifndef INT32_MAX
#define MAX_INDEX 0x7fffffffL
#else
#define MAX_INDEX INT32_MAX
#endif
#if (LONG_MAX < MAX_INDEX)
#error The integer format "t_index" must not have a greater range than "long int".
#endif
#if (INT_MAX > MAX_INDEX)
#error The integer format "int" must not have a greater range than "t_index".
#endif
typedef double t_float;
 
/* Method codes.
 
   These codes must agree with the METHODS array in fastcluster.R and the
   dictionary mthidx in fastcluster.py.
*/
enum method_codes {
   // non-Euclidean methods
   METHOD_METR_SINGLE = 0,
   METHOD_METR_COMPLETE = 1,
   METHOD_METR_AVERAGE = 2,
   METHOD_METR_WEIGHTED = 3,
   METHOD_METR_WARD = 4,
   METHOD_METR_WARD_D = METHOD_METR_WARD,
   METHOD_METR_CENTROID = 5,
   METHOD_METR_MEDIAN = 6,
   METHOD_METR_WARD_D2 = 7,
 
   MIN_METHOD_CODE = 0,
   MAX_METHOD_CODE = 7
};
 
enum method_codes_vector {
   // Euclidean methods
   METHOD_VECTOR_SINGLE = 0,
   METHOD_VECTOR_WARD = 1,
   METHOD_VECTOR_CENTROID = 2,
   METHOD_VECTOR_MEDIAN = 3,
 
   MIN_METHOD_VECTOR_CODE = 0,
   MAX_METHOD_VECTOR_CODE = 3
};
 
// self-destructing array pointer
template <typename type>
class auto_array_ptr {
private:
   type *ptr;
   auto_array_ptr(auto_array_ptr const &);            // non construction-copyable
   auto_array_ptr &operator=(auto_array_ptr const &); // non copyable
public:
   auto_array_ptr() : ptr(NULL) {}
   template <typename index>
   auto_array_ptr(index const size) : ptr(new type[size])
   {
   }
   template <typename index, typename value>
   auto_array_ptr(index const size, value const val) : ptr(new type[size])
   {
      FILL_N(ptr, size, val);
   }
   ~auto_array_ptr() { delete[] ptr; }
   void free()
   {
      delete[] ptr;
      ptr = NULL;
   }
   template <typename index>
   void init(index const size)
   {
      ptr = new type[size];
   }
   template <typename index, typename value>
   void init(index const size, value const val)
   {
      init(size);
      FILL_N(ptr, size, val);
   }
   inline operator type *() const { return ptr; }
};
 
struct node {
   t_index node1, node2;
   t_float dist;
};
 
inline bool operator<(const node a, const node b)
{
   return (a.dist < b.dist);
}
 
class cluster_result {
private:
   auto_array_ptr<node> Z;
   t_index pos;
 
public:
   cluster_result(const t_index size) : Z(size), pos(0) {}
 
   void append(const t_index node1, const t_index node2, const t_float dist)
   {
      Z[pos].node1 = node1;
      Z[pos].node2 = node2;
      Z[pos].dist = dist;
      ++pos;
   }
 
   node *operator[](const t_index idx) const { return Z + idx; }
 
   /* Define several methods to postprocess the distances. All these functions
      are monotone, so they do not change the sorted order of distances. */
 
   void sqrt() const
   {
      for (node *ZZ = Z; ZZ != Z + pos; ++ZZ) {
         ZZ->dist = std::sqrt(ZZ->dist);
      }
   }
 
   void sqrt(const t_float) const
   { // ignore the argument
      sqrt();
   }
 
   void sqrtdouble(const t_float) const
   { // ignore the argument
      for (node *ZZ = Z; ZZ != Z + pos; ++ZZ) {
         ZZ->dist = std::sqrt(2 * ZZ->dist);
      }
   }
 
#ifdef R_pow
#define my_pow R_pow
#else
#define my_pow std::pow
#endif
 
   void power(const t_float p) const
   {
      t_float const q = 1 / p;
      for (node *ZZ = Z; ZZ != Z + pos; ++ZZ) {
         ZZ->dist = my_pow(ZZ->dist, q);
      }
   }
 
   void plusone(const t_float) const
   { // ignore the argument
      for (node *ZZ = Z; ZZ != Z + pos; ++ZZ) {
         ZZ->dist += 1;
      }
   }
 
   void divide(const t_float denom) const
   {
      for (node *ZZ = Z; ZZ != Z + pos; ++ZZ) {
         ZZ->dist /= denom;
      }
   }
};
 
class doubly_linked_list {
   /*
     Class for a doubly linked list. Initially, the list is the integer range
     [0, size]. We provide a forward iterator and a method to delete an index
     from the list.
 
     Typical use: for (i=L.start; L<size; i=L.succ[I])
     or
     for (i=somevalue; L<size; i=L.succ[I])
   */
public:
   t_index start;
   auto_array_ptr<t_index> succ;
 
private:
   auto_array_ptr<t_index> pred;
   // Not necessarily private, we just do not need it in this instance.
 
public:
   doubly_linked_list(const t_index size)
      // Initialize to the given size.
      : start(0), succ(size + 1), pred(size + 1)
   {
      for (t_index i = 0; i < size; ++i) {
         pred[i + 1] = i;
         succ[i] = i + 1;
      }
      // pred[0] is never accessed!
      // succ[size] is never accessed!
   }
 
   ~doubly_linked_list() {}
 
   void remove(const t_index idx)
   {
      // Remove an index from the list.
      if (idx == start) {
         start = succ[idx];
      } else {
         succ[pred[idx]] = succ[idx];
         pred[succ[idx]] = pred[idx];
      }
      succ[idx] = 0; // Mark as inactive
   }
 
   bool is_inactive(t_index idx) const { return (succ[idx] == 0); }
};
 
// Indexing functions
// D is the upper triangular part of a symmetric (NxN)-matrix
// We require r_ < c_ !
#define D_(r_, c_) (D[(static_cast<std::ptrdiff_t>(2 * N - 3 - (r_)) * (r_) >> 1) + (c_)-1])
// Z is an ((N-1)x4)-array
#define Z_(_r, _c) (Z[(_r)*4 + (_c)])
 
/*
  Lookup function for a union-find data structure.
 
  The function finds the root of idx by going iteratively through all
  parent elements until a root is found. An element i is a root if
  nodes[i] is zero. To make subsequent searches faster, the entry for
  idx and all its parents is updated with the root element.
 */
class union_find {
private:
   auto_array_ptr<t_index> parent;
   t_index nextparent;
 
public:
   union_find(const t_index size) : parent(size > 0 ? 2 * size - 1 : 0, 0), nextparent(size) {}
 
   t_index Find(t_index idx) const
   {
      // NOLINTNEXTLINE
      if (parent[idx] != 0) { // a → b
         t_index p = idx;
         idx = parent[idx];
         if (parent[idx] != 0) { // a → b → c
            do {
               idx = parent[idx];
            } while (parent[idx] != 0);
            do {
               t_index tmp = parent[p];
               parent[p] = idx;
               p = tmp;
            } while (parent[p] != idx);
         }
      }
      return idx;
   }
 
   void Union(const t_index node1, const t_index node2) { parent[node1] = parent[node2] = nextparent++; }
};
 
class nan_error {
};
#ifdef FE_INVALID
class fenv_error {
};
#endif
 
static void MST_linkage_core(const t_index N, const t_float *const D, cluster_result &Z2)
{
   /*
       N: integer, number of data points
       D: condensed distance matrix N*(N-1)/2
       Z2: output data structure
 
       The basis of this algorithm is an algorithm by Rohlf:
 
       F. James Rohlf, Hierarchical clustering using the minimum spanning tree,
       The Computer Journal, vol. 16, 1973, p. 93–95.
   */
   t_index i;
   t_index idx2;
   doubly_linked_list active_nodes(N);
   auto_array_ptr<t_float> d(N);
 
   t_index prev_node;
   t_float min;
 
   // first iteration
   idx2 = 1;
   min = std::numeric_limits<t_float>::infinity();
   for (i = 1; i < N; ++i) {
      d[i] = D[i - 1];
#if HAVE_DIAGNOSTIC
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wfloat-equal"
#endif
      if (d[i] < min) {
         min = d[i];
         idx2 = i;
      } else if (fc_isnan(d[i]))
         throw(nan_error());
#if HAVE_DIAGNOSTIC
#pragma GCC diagnostic pop
#endif
   }
   Z2.append(0, idx2, min);
 
   for (t_index j = 1; j < N - 1; ++j) {
      prev_node = idx2;
      active_nodes.remove(prev_node);
 
      idx2 = active_nodes.succ[0];
      min = d[idx2];
      for (i = idx2; i < prev_node; i = active_nodes.succ[i]) {
         t_float tmp = D_(i, prev_node);
#if HAVE_DIAGNOSTIC
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wfloat-equal"
#endif
         if (tmp < d[i])
            d[i] = tmp;
         else if (fc_isnan(tmp))
            throw(nan_error());
#if HAVE_DIAGNOSTIC
#pragma GCC diagnostic pop
#endif
         if (d[i] < min) {
            min = d[i];
            idx2 = i;
         }
      }
      for (; i < N; i = active_nodes.succ[i]) {
         t_float tmp = D_(prev_node, i);
#if HAVE_DIAGNOSTIC
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wfloat-equal"
#endif
         if (d[i] > tmp)
            d[i] = tmp;
         else if (fc_isnan(tmp))
            throw(nan_error());
#if HAVE_DIAGNOSTIC
#pragma GCC diagnostic pop
#endif
         if (d[i] < min) {
            min = d[i];
            idx2 = i;
         }
      }
      Z2.append(prev_node, idx2, min);
   }
}
 
/* Functions for the update of the dissimilarity array */
 
inline static void f_single(t_float *const b, const t_float a)
{
   if (*b > a)
      *b = a;
}
inline static void f_complete(t_float *const b, const t_float a)
{
   if (*b < a)
      *b = a;
}
inline static void f_average(t_float *const b, const t_float a, const t_float s, const t_float t)
{
   *b = s * a + t * (*b);
#ifndef FE_INVALID
#if HAVE_DIAGNOSTIC
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wfloat-equal"
#endif
   if (fc_isnan(*b)) {
      throw(nan_error());
   }
#if HAVE_DIAGNOSTIC
#pragma GCC diagnostic pop
#endif
#endif
}
inline static void f_weighted(t_float *const b, const t_float a)
{
   *b = (a + *b) * .5;
#ifndef FE_INVALID
#if HAVE_DIAGNOSTIC
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wfloat-equal"
#endif
   if (fc_isnan(*b)) {
      throw(nan_error());
   }
#if HAVE_DIAGNOSTIC
#pragma GCC diagnostic pop
#endif
#endif
}
inline static void
f_ward(t_float *const b, const t_float a, const t_float c, const t_float s, const t_float t, const t_float v)
{
   *b = ((v + s) * a - v * c + (v + t) * (*b)) / (s + t + v);
//*b = a+(*b)-(t*a+s*(*b)+v*c)/(s+t+v);
#ifndef FE_INVALID
#if HAVE_DIAGNOSTIC
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wfloat-equal"
#endif
   if (fc_isnan(*b)) {
      throw(nan_error());
   }
#if HAVE_DIAGNOSTIC
#pragma GCC diagnostic pop
#endif
#endif
}
inline static void f_centroid(t_float *const b, const t_float a, const t_float stc, const t_float s, const t_float t)
{
   *b = s * a - stc + t * (*b);
#ifndef FE_INVALID
   if (fc_isnan(*b)) {
      throw(nan_error());
   }
#if HAVE_DIAGNOSTIC
#pragma GCC diagnostic pop
#endif
#endif
}
inline static void f_median(t_float *const b, const t_float a, const t_float c_4)
{
   *b = (a + (*b)) * .5 - c_4;
#ifndef FE_INVALID
#if HAVE_DIAGNOSTIC
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wfloat-equal"
#endif
   if (fc_isnan(*b)) {
      throw(nan_error());
   }
#if HAVE_DIAGNOSTIC
#pragma GCC diagnostic pop
#endif
#endif
}
 
template <method_codes method, typename t_members>
static void NN_chain_core(const t_index N, t_float *const D, t_members *const members, cluster_result &Z2)
{
   /*
       N: integer
       D: condensed distance matrix N*(N-1)/2
       Z2: output data structure
 
       This is the NN-chain algorithm, described on page 86 in the following book:
 
       Fionn Murtagh, Multidimensional Clustering Algorithms,
       Vienna, Würzburg: Physica-Verlag, 1985.
   */
   t_index i;
 
   auto_array_ptr<t_index> NN_chain(N);
   t_index NN_chain_tip = 0;
 
   t_index idx1, idx2;
 
   t_float size1, size2;
   doubly_linked_list active_nodes(N);
 
   t_float min;
 
   for (t_float const *DD = D; DD != D + (static_cast<std::ptrdiff_t>(N) * (N - 1) >> 1); ++DD) {
#if HAVE_DIAGNOSTIC
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wfloat-equal"
#endif
      if (fc_isnan(*DD)) {
         throw(nan_error());
      }
#if HAVE_DIAGNOSTIC
#pragma GCC diagnostic pop
#endif
   }
 
#ifdef FE_INVALID
   if (feclearexcept(FE_INVALID))
      throw fenv_error();
#endif
 
   for (t_index j = 0; j < N - 1; ++j) {
      if (NN_chain_tip <= 3) {
         NN_chain[0] = idx1 = active_nodes.start;
         NN_chain_tip = 1;
 
         idx2 = active_nodes.succ[idx1];
         min = D_(idx1, idx2);
         for (i = active_nodes.succ[idx2]; i < N; i = active_nodes.succ[i]) {
            if (D_(idx1, i) < min) {
               min = D_(idx1, i);
               idx2 = i;
            }
         }
      } // a: idx1   b: idx2
      else {
         NN_chain_tip -= 3;
         idx1 = NN_chain[NN_chain_tip - 1];
         idx2 = NN_chain[NN_chain_tip];
         min = idx1 < idx2 ? D_(idx1, idx2) : D_(idx2, idx1);
      } // a: idx1   b: idx2
 
      do {
         NN_chain[NN_chain_tip] = idx2;
 
         for (i = active_nodes.start; i < idx2; i = active_nodes.succ[i]) {
            if (D_(i, idx2) < min) {
               min = D_(i, idx2);
               idx1 = i;
            }
         }
         for (i = active_nodes.succ[idx2]; i < N; i = active_nodes.succ[i]) {
            if (D_(idx2, i) < min) {
               min = D_(idx2, i);
               idx1 = i;
            }
         }
 
         idx2 = idx1;
         idx1 = NN_chain[NN_chain_tip++];
 
      } while (idx2 != NN_chain[NN_chain_tip - 2]);
 
      Z2.append(idx1, idx2, min);
 
      if (idx1 > idx2) {
         t_index tmp = idx1;
         idx1 = idx2;
         idx2 = tmp;
      }
 
      if (method == METHOD_METR_AVERAGE || method == METHOD_METR_WARD) {
         size1 = static_cast<t_float>(members[idx1]);
         size2 = static_cast<t_float>(members[idx2]);
         members[idx2] += members[idx1];
      }
 
      // Remove the smaller index from the valid indices (active_nodes).
      active_nodes.remove(idx1);
 
      switch (method) {
      case METHOD_METR_SINGLE:
         /*
         Single linkage.
 
         Characteristic: new distances are never longer than the old distances.
         */
         // Update the distance matrix in the range [start, idx1).
         for (i = active_nodes.start; i < idx1; i = active_nodes.succ[i])
            f_single(&D_(i, idx2), D_(i, idx1));
         // Update the distance matrix in the range (idx1, idx2).
         for (; i < idx2; i = active_nodes.succ[i])
            f_single(&D_(i, idx2), D_(idx1, i));
         // Update the distance matrix in the range (idx2, N).
         for (i = active_nodes.succ[idx2]; i < N; i = active_nodes.succ[i])
            f_single(&D_(idx2, i), D_(idx1, i));
         break;
 
      case METHOD_METR_COMPLETE:
         /*
         Complete linkage.
 
         Characteristic: new distances are never shorter than the old distances.
         */
         // Update the distance matrix in the range [start, idx1).
         for (i = active_nodes.start; i < idx1; i = active_nodes.succ[i])
            f_complete(&D_(i, idx2), D_(i, idx1));
         // Update the distance matrix in the range (idx1, idx2).
         for (; i < idx2; i = active_nodes.succ[i])
            f_complete(&D_(i, idx2), D_(idx1, i));
         // Update the distance matrix in the range (idx2, N).
         for (i = active_nodes.succ[idx2]; i < N; i = active_nodes.succ[i])
            f_complete(&D_(idx2, i), D_(idx1, i));
         break;
 
      case METHOD_METR_AVERAGE: {
         /*
         Average linkage.
 
         Shorter and longer distances can occur.
         */
         // Update the distance matrix in the range [start, idx1).
         t_float s = size1 / (size1 + size2);
         t_float t = size2 / (size1 + size2);
         for (i = active_nodes.start; i < idx1; i = active_nodes.succ[i])
            f_average(&D_(i, idx2), D_(i, idx1), s, t);
         // Update the distance matrix in the range (idx1, idx2).
         for (; i < idx2; i = active_nodes.succ[i])
            f_average(&D_(i, idx2), D_(idx1, i), s, t);
         // Update the distance matrix in the range (idx2, N).
         for (i = active_nodes.succ[idx2]; i < N; i = active_nodes.succ[i])
            f_average(&D_(idx2, i), D_(idx1, i), s, t);
         break;
      }
 
      case METHOD_METR_WEIGHTED:
         /*
         Weighted linkage.
 
         Shorter and longer distances can occur.
         */
         // Update the distance matrix in the range [start, idx1).
         for (i = active_nodes.start; i < idx1; i = active_nodes.succ[i])
            f_weighted(&D_(i, idx2), D_(i, idx1));
         // Update the distance matrix in the range (idx1, idx2).
         for (; i < idx2; i = active_nodes.succ[i])
            f_weighted(&D_(i, idx2), D_(idx1, i));
         // Update the distance matrix in the range (idx2, N).
         for (i = active_nodes.succ[idx2]; i < N; i = active_nodes.succ[i])
            f_weighted(&D_(idx2, i), D_(idx1, i));
         break;
 
      case METHOD_METR_WARD:
         /*
         Ward linkage.
 
         Shorter and longer distances can occur, not smaller than min(d1,d2)
         but maybe bigger than max(d1,d2).
         */
         // Update the distance matrix in the range [start, idx1).
         // t_float v = static_cast<t_float>(members[i]);
         for (i = active_nodes.start; i < idx1; i = active_nodes.succ[i])
            f_ward(&D_(i, idx2), D_(i, idx1), min, size1, size2, static_cast<t_float>(members[i]));
         // Update the distance matrix in the range (idx1, idx2).
         for (; i < idx2; i = active_nodes.succ[i])
            f_ward(&D_(i, idx2), D_(idx1, i), min, size1, size2, static_cast<t_float>(members[i]));
         // Update the distance matrix in the range (idx2, N).
         for (i = active_nodes.succ[idx2]; i < N; i = active_nodes.succ[i])
            f_ward(&D_(idx2, i), D_(idx1, i), min, size1, size2, static_cast<t_float>(members[i]));
         break;
 
      default: throw std::runtime_error(std::string("Invalid method."));
      }
   }
#ifdef FE_INVALID
   if (fetestexcept(FE_INVALID))
      throw fenv_error();
#endif
}
 
class binary_min_heap {
   /*
   Class for a binary min-heap. The data resides in an array A. The elements of
   A are not changed but two lists I and R of indices are generated which point
   to elements of A and backwards.
 
   The heap tree structure is
 
      H[2*i+1]     H[2*i+2]
          \            /
           \          /
            ≤        ≤
             \      /
              \    /
               H[i]
 
   where the children must be less or equal than their parent. Thus, H[0]
   contains the minimum. The lists I and R are made such that H[i] = A[I[i]]
   and R[I[i]] = i.
 
   This implementation is not designed to handle NaN values.
   */
private:
   t_float *const A;
   t_index size;
   auto_array_ptr<t_index> I;
   auto_array_ptr<t_index> R;
 
   // no default constructor
   binary_min_heap();
   // noncopyable
   binary_min_heap(binary_min_heap const &);
   binary_min_heap &operator=(binary_min_heap const &);
 
public:
   binary_min_heap(t_float *const A_, const t_index size_) : A(A_), size(size_), I(size), R(size)
   { // Allocate memory and initialize the lists I and R to the identity. This
      // does not make it a heap. Call heapify afterwards!
      for (t_index i = 0; i < size; ++i)
         R[i] = I[i] = i;
   }
 
   binary_min_heap(t_float *const A_, const t_index size1, const t_index size2, const t_index start)
      : A(A_), size(size1), I(size1), R(size2)
   { // Allocate memory and initialize the lists I and R to the identity. This
      // does not make it a heap. Call heapify afterwards!
      for (t_index i = 0; i < size; ++i) {
         R[i + start] = i;
         I[i] = i + start;
      }
   }
 
   ~binary_min_heap() {}
 
   void heapify()
   {
      // Arrange the indices I and R so that H[i] := A[I[i]] satisfies the heap
      // condition H[i] < H[2*i+1] and H[i] < H[2*i+2] for each i.
      //
      // Complexity: Θ(size)
      // Reference: Cormen, Leiserson, Rivest, Stein, Introduction to Algorithms,
      // 3rd ed., 2009, Section 6.3 “Building a heap”
      t_index idx;
      for (idx = (size >> 1); idx > 0;) {
         --idx;
         update_geq_(idx);
      }
   }
 
   inline t_index argmin() const
   {
      // Return the minimal element.
      return I[0];
   }
 
   void heap_pop()
   {
      // Remove the minimal element from the heap.
      --size;
      I[0] = I[size];
      R[I[0]] = 0;
      update_geq_(0);
   }
 
   void remove(t_index idx)
   {
      // Remove an element from the heap.
      --size;
      R[I[size]] = R[idx];
      I[R[idx]] = I[size];
      if (H(size) <= A[idx]) {
         update_leq_(R[idx]);
      } else {
         update_geq_(R[idx]);
      }
   }
 
   void replace(const t_index idxold, const t_index idxnew, const t_float val)
   {
      R[idxnew] = R[idxold];
      I[R[idxnew]] = idxnew;
      if (val <= A[idxold])
         update_leq(idxnew, val);
      else
         update_geq(idxnew, val);
   }
 
   void update(const t_index idx, const t_float val) const
   {
      // Update the element A[i] with val and re-arrange the indices to preserve
      // the heap condition.
      if (val <= A[idx])
         update_leq(idx, val);
      else
         update_geq(idx, val);
   }
 
   void update_leq(const t_index idx, const t_float val) const
   {
      // Use this when the new value is not more than the old value.
      A[idx] = val;
      update_leq_(R[idx]);
   }
 
   void update_geq(const t_index idx, const t_float val) const
   {
      // Use this when the new value is not less than the old value.
      A[idx] = val;
      update_geq_(R[idx]);
   }
 
private:
   void update_leq_(t_index i) const
   {
      t_index j;
      for (; (i > 0) && (H(i) < H(j = (i - 1) >> 1)); i = j)
         heap_swap(i, j);
   }
 
   void update_geq_(t_index i) const
   {
      t_index j;
      for (; (j = 2 * i + 1) < size; i = j) {
         if (H(j) >= H(i)) {
            ++j;
            if (j >= size || H(j) >= H(i))
               break;
         } else if (j + 1 < size && H(j + 1) < H(j))
            ++j;
         heap_swap(i, j);
      }
   }
 
   void heap_swap(const t_index i, const t_index j) const
   {
      // Swap two indices.
      t_index tmp = I[i];
      I[i] = I[j];
      I[j] = tmp;
      R[I[i]] = i;
      R[I[j]] = j;
   }
 
   inline t_float H(const t_index i) const { return A[I[i]]; }
};
 
template <method_codes method, typename t_members>
static void generic_linkage(const t_index N, t_float *const D, t_members *const members, cluster_result &Z2)
{
   /*
     N: integer, number of data points
     D: condensed distance matrix N*(N-1)/2
     Z2: output data structure
   */
 
   const t_index N_1 = N - 1;
   t_index i, j;       // loop variables
   t_index idx1, idx2; // row and column indices
 
   auto_array_ptr<t_index> n_nghbr(N_1); // array of nearest neighbors
   auto_array_ptr<t_float> mindist(N_1); // distances to the nearest neighbors
   auto_array_ptr<t_index> row_repr(N);  // row_repr[i]: node number that the
                                         // i-th row represents
   doubly_linked_list active_nodes(N);
   binary_min_heap nn_distances(&*mindist, N_1); // minimum heap structure for
                                                 // the distance to the nearest neighbor of each point
   t_index node1, node2;                         // node numbers in the output
   t_float size1, size2;                         // and their cardinalities
 
   t_float min; // minimum and row index for nearest-neighbor search
   t_index idx;
 
   for (i = 0; i < N; ++i)
      // Build a list of row ↔ node label assignments.
      // Initially i ↦ i
      row_repr[i] = i;
 
   // Initialize the minimal distances:
   // Find the nearest neighbor of each point.
   // n_nghbr[i] = argmin_{j>i} D(i,j) for i in range(N-1)
   t_float const *DD = D;
   for (i = 0; i < N_1; ++i) {
      min = std::numeric_limits<t_float>::infinity();
      for (idx = j = i + 1; j < N; ++j, ++DD) {
#if HAVE_DIAGNOSTIC
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wfloat-equal"
#endif
         if (*DD < min) {
            min = *DD;
            idx = j;
         } else if (fc_isnan(*DD))
            throw(nan_error());
      }
#if HAVE_DIAGNOSTIC
#pragma GCC diagnostic pop
#endif
      mindist[i] = min;
      n_nghbr[i] = idx;
   }
 
   // Put the minimal distances into a heap structure to make the repeated
   // global minimum searches fast.
   nn_distances.heapify();
 
#ifdef FE_INVALID
   if (feclearexcept(FE_INVALID))
      throw fenv_error();
#endif
 
   // Main loop: We have N-1 merging steps.
   for (i = 0; i < N_1; ++i) {
      /*
        Here is a special feature that allows fast bookkeeping and updates of the
        minimal distances.
 
        mindist[i] stores a lower bound on the minimum distance of the point i to
        all points of higher index:
 
            mindist[i] ≥ min_{j>i} D(i,j)
 
        Normally, we have equality. However, this minimum may become invalid due
        to the updates in the distance matrix. The rules are:
 
        1) If mindist[i] is equal to D(i, n_nghbr[i]), this is the correct
           minimum and n_nghbr[i] is a nearest neighbor.
 
        2) If mindist[i] is smaller than D(i, n_nghbr[i]), this might not be the
           correct minimum. The minimum needs to be recomputed.
 
        3) mindist[i] is never bigger than the true minimum. Hence, we never
           miss the true minimum if we take the smallest mindist entry,
           re-compute the value if necessary (thus maybe increasing it) and
           looking for the now smallest mindist entry until a valid minimal
           entry is found. This step is done in the lines below.
 
        The update process for D below takes care that these rules are
        fulfilled. This makes sure that the minima in the rows D(i,i+1:)of D are
        re-calculated when necessary but re-calculation is avoided whenever
        possible.
 
        The re-calculation of the minima makes the worst-case runtime of this
        algorithm cubic in N. We avoid this whenever possible, and in most cases
        the runtime appears to be quadratic.
      */
      idx1 = nn_distances.argmin();
      if (method != METHOD_METR_SINGLE) {
         while (mindist[idx1] < D_(idx1, n_nghbr[idx1])) {
            // Recompute the minimum mindist[idx1] and n_nghbr[idx1].
            n_nghbr[idx1] = j = active_nodes.succ[idx1]; // exists, maximally N-1
            min = D_(idx1, j);
            for (j = active_nodes.succ[j]; j < N; j = active_nodes.succ[j]) {
               if (D_(idx1, j) < min) {
                  min = D_(idx1, j);
                  n_nghbr[idx1] = j;
               }
            }
            /* Update the heap with the new true minimum and search for the
               (possibly different) minimal entry. */
            nn_distances.update_geq(idx1, min);
            idx1 = nn_distances.argmin();
         }
      }
 
      nn_distances.heap_pop(); // Remove the current minimum from the heap.
      idx2 = n_nghbr[idx1];
 
      // Write the newly found minimal pair of nodes to the output array.
      node1 = row_repr[idx1];
      node2 = row_repr[idx2];
 
      if (method == METHOD_METR_AVERAGE || method == METHOD_METR_WARD || method == METHOD_METR_CENTROID) {
         size1 = static_cast<t_float>(members[idx1]);
         size2 = static_cast<t_float>(members[idx2]);
         members[idx2] += members[idx1];
      }
      Z2.append(node1, node2, mindist[idx1]);
 
      // Remove idx1 from the list of active indices (active_nodes).
      active_nodes.remove(idx1);
      // Index idx2 now represents the new (merged) node with label N+i.
      row_repr[idx2] = N + i;
 
      // Update the distance matrix
      switch (method) {
      case METHOD_METR_SINGLE:
         /*
           Single linkage.
 
           Characteristic: new distances are never longer than the old distances.
         */
         // Update the distance matrix in the range [start, idx1).
         for (j = active_nodes.start; j < idx1; j = active_nodes.succ[j]) {
            f_single(&D_(j, idx2), D_(j, idx1));
            if (n_nghbr[j] == idx1)
               n_nghbr[j] = idx2;
         }
         // Update the distance matrix in the range (idx1, idx2).
         for (; j < idx2; j = active_nodes.succ[j]) {
            f_single(&D_(j, idx2), D_(idx1, j));
            // If the new value is below the old minimum in a row, update
            // the mindist and n_nghbr arrays.
            if (D_(j, idx2) < mindist[j]) {
               nn_distances.update_leq(j, D_(j, idx2));
               n_nghbr[j] = idx2;
            }
         }
         // Update the distance matrix in the range (idx2, N).
         // Recompute the minimum mindist[idx2] and n_nghbr[idx2].
         if (idx2 < N_1) {
            min = mindist[idx2];
            for (j = active_nodes.succ[idx2]; j < N; j = active_nodes.succ[j]) {
               f_single(&D_(idx2, j), D_(idx1, j));
               if (D_(idx2, j) < min) {
                  n_nghbr[idx2] = j;
                  min = D_(idx2, j);
               }
            }
            nn_distances.update_leq(idx2, min);
         }
         break;
 
      case METHOD_METR_COMPLETE:
         /*
           Complete linkage.
 
           Characteristic: new distances are never shorter than the old distances.
         */
         // Update the distance matrix in the range [start, idx1).
         for (j = active_nodes.start; j < idx1; j = active_nodes.succ[j]) {
            f_complete(&D_(j, idx2), D_(j, idx1));
            if (n_nghbr[j] == idx1)
               n_nghbr[j] = idx2;
         }
         // Update the distance matrix in the range (idx1, idx2).
         for (; j < idx2; j = active_nodes.succ[j])
            f_complete(&D_(j, idx2), D_(idx1, j));
         // Update the distance matrix in the range (idx2, N).
         for (j = active_nodes.succ[idx2]; j < N; j = active_nodes.succ[j])
            f_complete(&D_(idx2, j), D_(idx1, j));
         break;
 
      case METHOD_METR_AVERAGE: {
         /*
           Average linkage.
 
           Shorter and longer distances can occur.
         */
         // Update the distance matrix in the range [start, idx1).
         t_float s = size1 / (size1 + size2);
         t_float t = size2 / (size1 + size2);
         for (j = active_nodes.start; j < idx1; j = active_nodes.succ[j]) {
            f_average(&D_(j, idx2), D_(j, idx1), s, t);
            if (n_nghbr[j] == idx1)
               n_nghbr[j] = idx2;
         }
         // Update the distance matrix in the range (idx1, idx2).
         for (; j < idx2; j = active_nodes.succ[j]) {
            f_average(&D_(j, idx2), D_(idx1, j), s, t);
            if (D_(j, idx2) < mindist[j]) {
               nn_distances.update_leq(j, D_(j, idx2));
               n_nghbr[j] = idx2;
            }
         }
         // Update the distance matrix in the range (idx2, N).
         if (idx2 < N_1) {
            n_nghbr[idx2] = j = active_nodes.succ[idx2]; // exists, maximally N-1
            f_average(&D_(idx2, j), D_(idx1, j), s, t);
            min = D_(idx2, j);
            for (j = active_nodes.succ[j]; j < N; j = active_nodes.succ[j]) {
               f_average(&D_(idx2, j), D_(idx1, j), s, t);
               if (D_(idx2, j) < min) {
                  min = D_(idx2, j);
                  n_nghbr[idx2] = j;
               }
            }
            nn_distances.update(idx2, min);
         }
         break;
      }
 
      case METHOD_METR_WEIGHTED:
         /*
           Weighted linkage.
 
           Shorter and longer distances can occur.
         */
         // Update the distance matrix in the range [start, idx1).
         for (j = active_nodes.start; j < idx1; j = active_nodes.succ[j]) {
            f_weighted(&D_(j, idx2), D_(j, idx1));
            if (n_nghbr[j] == idx1)
               n_nghbr[j] = idx2;
         }
         // Update the distance matrix in the range (idx1, idx2).
         for (; j < idx2; j = active_nodes.succ[j]) {
            f_weighted(&D_(j, idx2), D_(idx1, j));
            if (D_(j, idx2) < mindist[j]) {
               nn_distances.update_leq(j, D_(j, idx2));
               n_nghbr[j] = idx2;
            }
         }
         // Update the distance matrix in the range (idx2, N).
         if (idx2 < N_1) {
            n_nghbr[idx2] = j = active_nodes.succ[idx2]; // exists, maximally N-1
            f_weighted(&D_(idx2, j), D_(idx1, j));
            min = D_(idx2, j);
            for (j = active_nodes.succ[j]; j < N; j = active_nodes.succ[j]) {
               f_weighted(&D_(idx2, j), D_(idx1, j));
               if (D_(idx2, j) < min) {
                  min = D_(idx2, j);
                  n_nghbr[idx2] = j;
               }
            }
            nn_distances.update(idx2, min);
         }
         break;
 
      case METHOD_METR_WARD:
         /*
           Ward linkage.
 
           Shorter and longer distances can occur, not smaller than min(d1,d2)
           but maybe bigger than max(d1,d2).
         */
         // Update the distance matrix in the range [start, idx1).
         for (j = active_nodes.start; j < idx1; j = active_nodes.succ[j]) {
            f_ward(&D_(j, idx2), D_(j, idx1), mindist[idx1], size1, size2, static_cast<t_float>(members[j]));
            if (n_nghbr[j] == idx1)
               n_nghbr[j] = idx2;
         }
         // Update the distance matrix in the range (idx1, idx2).
         for (; j < idx2; j = active_nodes.succ[j]) {
            f_ward(&D_(j, idx2), D_(idx1, j), mindist[idx1], size1, size2, static_cast<t_float>(members[j]));
            if (D_(j, idx2) < mindist[j]) {
               nn_distances.update_leq(j, D_(j, idx2));
               n_nghbr[j] = idx2;
            }
         }
         // Update the distance matrix in the range (idx2, N).
         if (idx2 < N_1) {
            n_nghbr[idx2] = j = active_nodes.succ[idx2]; // exists, maximally N-1
            f_ward(&D_(idx2, j), D_(idx1, j), mindist[idx1], size1, size2, static_cast<t_float>(members[j]));
            min = D_(idx2, j);
            for (j = active_nodes.succ[j]; j < N; j = active_nodes.succ[j]) {
               f_ward(&D_(idx2, j), D_(idx1, j), mindist[idx1], size1, size2, static_cast<t_float>(members[j]));
               if (D_(idx2, j) < min) {
                  min = D_(idx2, j);
                  n_nghbr[idx2] = j;
               }
            }
            nn_distances.update(idx2, min);
         }
         break;
 
      case METHOD_METR_CENTROID: {
         /*
           Centroid linkage.
 
           Shorter and longer distances can occur, not bigger than max(d1,d2)
           but maybe smaller than min(d1,d2).
         */
         // Update the distance matrix in the range [start, idx1).
         t_float s = size1 / (size1 + size2);
         t_float t = size2 / (size1 + size2);
         t_float stc = s * t * mindist[idx1];
         for (j = active_nodes.start; j < idx1; j = active_nodes.succ[j]) {
            f_centroid(&D_(j, idx2), D_(j, idx1), stc, s, t);
            if (D_(j, idx2) < mindist[j]) {
               nn_distances.update_leq(j, D_(j, idx2));
               n_nghbr[j] = idx2;
            } else if (n_nghbr[j] == idx1)
               n_nghbr[j] = idx2;
         }
         // Update the distance matrix in the range (idx1, idx2).
         for (; j < idx2; j = active_nodes.succ[j]) {
            f_centroid(&D_(j, idx2), D_(idx1, j), stc, s, t);
            if (D_(j, idx2) < mindist[j]) {
               nn_distances.update_leq(j, D_(j, idx2));
               n_nghbr[j] = idx2;
            }
         }
         // Update the distance matrix in the range (idx2, N).
         if (idx2 < N_1) {
            n_nghbr[idx2] = j = active_nodes.succ[idx2]; // exists, maximally N-1
            f_centroid(&D_(idx2, j), D_(idx1, j), stc, s, t);
            min = D_(idx2, j);
            for (j = active_nodes.succ[j]; j < N; j = active_nodes.succ[j]) {
               f_centroid(&D_(idx2, j), D_(idx1, j), stc, s, t);
               if (D_(idx2, j) < min) {
                  min = D_(idx2, j);
                  n_nghbr[idx2] = j;
               }
            }
            nn_distances.update(idx2, min);
         }
         break;
      }
 
      case METHOD_METR_MEDIAN: {
         /*
           Median linkage.
 
           Shorter and longer distances can occur, not bigger than max(d1,d2)
           but maybe smaller than min(d1,d2).
         */
         // Update the distance matrix in the range [start, idx1).
         t_float c_4 = mindist[idx1] * .25;
         for (j = active_nodes.start; j < idx1; j = active_nodes.succ[j]) {
            f_median(&D_(j, idx2), D_(j, idx1), c_4);
            if (D_(j, idx2) < mindist[j]) {
               nn_distances.update_leq(j, D_(j, idx2));
               n_nghbr[j] = idx2;
            } else if (n_nghbr[j] == idx1)
               n_nghbr[j] = idx2;
         }
         // Update the distance matrix in the range (idx1, idx2).
         for (; j < idx2; j = active_nodes.succ[j]) {
            f_median(&D_(j, idx2), D_(idx1, j), c_4);
            if (D_(j, idx2) < mindist[j]) {
               nn_distances.update_leq(j, D_(j, idx2));
               n_nghbr[j] = idx2;
            }
         }
         // Update the distance matrix in the range (idx2, N).
         if (idx2 < N_1) {
            n_nghbr[idx2] = j = active_nodes.succ[idx2]; // exists, maximally N-1
            f_median(&D_(idx2, j), D_(idx1, j), c_4);
            min = D_(idx2, j);
            for (j = active_nodes.succ[j]; j < N; j = active_nodes.succ[j]) {
               f_median(&D_(idx2, j), D_(idx1, j), c_4);
               if (D_(idx2, j) < min) {
                  min = D_(idx2, j);
                  n_nghbr[idx2] = j;
               }
            }
            nn_distances.update(idx2, min);
         }
         break;
      }
 
      default: throw std::runtime_error(std::string("Invalid method."));
      }
   }
#ifdef FE_INVALID
   if (fetestexcept(FE_INVALID))
      throw fenv_error();
#endif
}
 
/*
  Clustering methods for vector data
*/
 
template <typename t_dissimilarity>
static void MST_linkage_core_vector(const t_index N, t_dissimilarity &dist, cluster_result &Z2)
{
   /*
       N: integer, number of data points
       dist: function pointer to the metric
       Z2: output data structure
 
       The basis of this algorithm is an algorithm by Rohlf:
 
       F. James Rohlf, Hierarchical clustering using the minimum spanning tree,
       The Computer Journal, vol. 16, 1973, p. 93–95.
   */
   t_index i;
   t_index idx2;
   doubly_linked_list active_nodes(N);
   auto_array_ptr<t_float> d(N);
 
   t_index prev_node;
   t_float min;
 
   // first iteration
   idx2 = 1;
   min = std::numeric_limits<t_float>::infinity();
   for (i = 1; i < N; ++i) {
      d[i] = dist(0, i);
#if HAVE_DIAGNOSTIC
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wfloat-equal"
#endif
      if (d[i] < min) {
         min = d[i];
         idx2 = i;
      } else if (fc_isnan(d[i]))
         throw(nan_error());
#if HAVE_DIAGNOSTIC
#pragma GCC diagnostic pop
#endif
   }
 
   Z2.append(0, idx2, min);
 
   for (t_index j = 1; j < N - 1; ++j) {
      prev_node = idx2;
      active_nodes.remove(prev_node);
 
      idx2 = active_nodes.succ[0];
      min = d[idx2];
 
      for (i = idx2; i < N; i = active_nodes.succ[i]) {
         t_float tmp = dist(i, prev_node);
#if HAVE_DIAGNOSTIC
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wfloat-equal"
#endif
         if (d[i] > tmp)
            d[i] = tmp;
         else if (fc_isnan(tmp))
            throw(nan_error());
#if HAVE_DIAGNOSTIC
#pragma GCC diagnostic pop
#endif
         if (d[i] < min) {
            min = d[i];
            idx2 = i;
         }
      }
      Z2.append(prev_node, idx2, min);
   }
}
 
template <method_codes_vector method, typename t_dissimilarity>
static void generic_linkage_vector(const t_index N, t_dissimilarity &dist, cluster_result &Z2)
{
   /*
     N: integer, number of data points
     dist: function pointer to the metric
     Z2: output data structure
 
     This algorithm is valid for the distance update methods
     "Ward", "centroid" and "median" only!
   */
   const t_index N_1 = N - 1;
   t_index i, j;       // loop variables
   t_index idx1, idx2; // row and column indices
 
   auto_array_ptr<t_index> n_nghbr(N_1); // array of nearest neighbors
   auto_array_ptr<t_float> mindist(N_1); // distances to the nearest neighbors
   auto_array_ptr<t_index> row_repr(N);  // row_repr[i]: node number that the
                                         // i-th row represents
   doubly_linked_list active_nodes(N);
   binary_min_heap nn_distances(&*mindist, N_1); // minimum heap structure for
                                                 // the distance to the nearest neighbor of each point
   t_index node1, node2;                         // node numbers in the output
   t_float min;                                  // minimum and row index for nearest-neighbor search
 
   for (i = 0; i < N; ++i)
      // Build a list of row ↔ node label assignments.
      // Initially i ↦ i
      row_repr[i] = i;
 
   // Initialize the minimal distances:
   // Find the nearest neighbor of each point.
   // n_nghbr[i] = argmin_{j>i} D(i,j) for i in range(N-1)
   for (i = 0; i < N_1; ++i) {
      min = std::numeric_limits<t_float>::infinity();
      t_index idx;
      for (idx = j = i + 1; j < N; ++j) {
         t_float tmp;
         switch (method) {
         case METHOD_VECTOR_WARD: tmp = dist.ward_initial(i, j); break;
         default: tmp = dist.template sqeuclidean<true>(i, j);
         }
         if (tmp < min) {
            min = tmp;
            idx = j;
         }
      }
      switch (method) {
      case METHOD_VECTOR_WARD: mindist[i] = t_dissimilarity::ward_initial_conversion(min); break;
      default: mindist[i] = min;
      }
      n_nghbr[i] = idx;
   }
 
   // Put the minimal distances into a heap structure to make the repeated
   // global minimum searches fast.
   nn_distances.heapify();
 
   // Main loop: We have N-1 merging steps.
   for (i = 0; i < N_1; ++i) {
      idx1 = nn_distances.argmin();
 
      while (active_nodes.is_inactive(n_nghbr[idx1])) {
         // Recompute the minimum mindist[idx1] and n_nghbr[idx1].
         n_nghbr[idx1] = j = active_nodes.succ[idx1]; // exists, maximally N-1
         switch (method) {
         case METHOD_VECTOR_WARD:
            min = dist.ward(idx1, j);
            for (j = active_nodes.succ[j]; j < N; j = active_nodes.succ[j]) {
               t_float const tmp = dist.ward(idx1, j);
               if (tmp < min) {
                  min = tmp;
                  n_nghbr[idx1] = j;
               }
            }
            break;
         default:
            min = dist.template sqeuclidean<true>(idx1, j);
            for (j = active_nodes.succ[j]; j < N; j = active_nodes.succ[j]) {
               t_float const tmp = dist.template sqeuclidean<true>(idx1, j);
               if (tmp < min) {
                  min = tmp;
                  n_nghbr[idx1] = j;
               }
            }
         }
         /* Update the heap with the new true minimum and search for the (possibly
            different) minimal entry. */
         nn_distances.update_geq(idx1, min);
         idx1 = nn_distances.argmin();
      }
 
      nn_distances.heap_pop(); // Remove the current minimum from the heap.
      idx2 = n_nghbr[idx1];
 
      // Write the newly found minimal pair of nodes to the output array.
      node1 = row_repr[idx1];
      node2 = row_repr[idx2];
 
      Z2.append(node1, node2, mindist[idx1]);
 
      switch (method) {
      case METHOD_VECTOR_WARD:
      case METHOD_VECTOR_CENTROID: dist.merge_inplace(idx1, idx2); break;
      case METHOD_VECTOR_MEDIAN: dist.merge_inplace_weighted(idx1, idx2); break;
      default: throw std::runtime_error(std::string("Invalid method."));
      }
 
      // Index idx2 now represents the new (merged) node with label N+i.
      row_repr[idx2] = N + i;
      // Remove idx1 from the list of active indices (active_nodes).
      active_nodes.remove(idx1); // TBD later!!!
 
      // Update the distance matrix
      switch (method) {
      case METHOD_VECTOR_WARD:
         /*
           Ward linkage.
 
           Shorter and longer distances can occur, not smaller than min(d1,d2)
           but maybe bigger than max(d1,d2).
         */
         // Update the distance matrix in the range [start, idx1).
         for (j = active_nodes.start; j < idx1; j = active_nodes.succ[j]) {
            if (n_nghbr[j] == idx2) {
               n_nghbr[j] = idx1; // invalidate
            }
         }
         // Update the distance matrix in the range (idx1, idx2).
         for (; j < idx2; j = active_nodes.succ[j]) {
            t_float const tmp = dist.ward(j, idx2);
            if (tmp < mindist[j]) {
               nn_distances.update_leq(j, tmp);
               n_nghbr[j] = idx2;
            } else if (n_nghbr[j] == idx2) {
               n_nghbr[j] = idx1; // invalidate
            }
         }
         // Find the nearest neighbor for idx2.
         if (idx2 < N_1) {
            n_nghbr[idx2] = j = active_nodes.succ[idx2]; // exists, maximally N-1
            min = dist.ward(idx2, j);
            for (j = active_nodes.succ[j]; j < N; j = active_nodes.succ[j]) {
               t_float const tmp = dist.ward(idx2, j);
               if (tmp < min) {
                  min = tmp;
                  n_nghbr[idx2] = j;
               }
            }
            nn_distances.update(idx2, min);
         }
         break;
 
      default:
         /*
           Centroid and median linkage.
 
           Shorter and longer distances can occur, not bigger than max(d1,d2)
           but maybe smaller than min(d1,d2).
         */
         for (j = active_nodes.start; j < idx2; j = active_nodes.succ[j]) {
            t_float const tmp = dist.template sqeuclidean<true>(j, idx2);
            if (tmp < mindist[j]) {
               nn_distances.update_leq(j, tmp);
               n_nghbr[j] = idx2;
            } else if (n_nghbr[j] == idx2)
               n_nghbr[j] = idx1; // invalidate
         }
         // Find the nearest neighbor for idx2.
         if (idx2 < N_1) {
            n_nghbr[idx2] = j = active_nodes.succ[idx2]; // exists, maximally N-1
            min = dist.template sqeuclidean<true>(idx2, j);
            for (j = active_nodes.succ[j]; j < N; j = active_nodes.succ[j]) {
               t_float const tmp = dist.template sqeuclidean<true>(idx2, j);
               if (tmp < min) {
                  min = tmp;
                  n_nghbr[idx2] = j;
               }
            }
            nn_distances.update(idx2, min);
         }
      }
   }
}
 
template <method_codes_vector method, typename t_dissimilarity>
static void generic_linkage_vector_alternative(const t_index N, t_dissimilarity &dist, cluster_result &Z2)
{
   /*
     N: integer, number of data points
     dist: function pointer to the metric
     Z2: output data structure
 
     This algorithm is valid for the distance update methods
     "Ward", "centroid" and "median" only!
   */
   const t_index N_1 = N - 1;
   t_index i, j = 0;   // loop variables
   t_index idx1, idx2; // row and column indices
 
   auto_array_ptr<t_index> n_nghbr(2 * N - 2); // array of nearest neighbors
   auto_array_ptr<t_float> mindist(2 * N - 2); // distances to the nearest neighbors
 
   doubly_linked_list active_nodes(N + N_1);
   binary_min_heap nn_distances(&*mindist, N_1, 2 * N - 2,
                                1); // minimum heap
                                    // structure for the distance to the nearest neighbor of each point
 
   t_float min; // minimum for nearest-neighbor searches
 
   // Initialize the minimal distances:
   // Find the nearest neighbor of each point.
   // n_nghbr[i] = argmin_{j>i} D(i,j) for i in range(N-1)
   for (i = 1; i < N; ++i) {
      min = std::numeric_limits<t_float>::infinity();
      t_index idx;
      for (idx = j = 0; j < i; ++j) {
         t_float tmp;
         switch (method) {
         case METHOD_VECTOR_WARD: tmp = dist.ward_initial(i, j); break;
         default: tmp = dist.template sqeuclidean<true>(i, j);
         }
         if (tmp < min) {
            min = tmp;
            idx = j;
         }
      }
      switch (method) {
      case METHOD_VECTOR_WARD: mindist[i] = t_dissimilarity::ward_initial_conversion(min); break;
      default: mindist[i] = min;
      }
      n_nghbr[i] = idx;
   }
 
   // Put the minimal distances into a heap structure to make the repeated
   // global minimum searches fast.
   nn_distances.heapify();
 
   // Main loop: We have N-1 merging steps.
   for (i = N; i < N + N_1; ++i) {
      /*
        The bookkeeping is different from the "stored matrix approach" algorithm
        generic_linkage.
 
        mindist[i] stores a lower bound on the minimum distance of the point i to
        all points of *lower* index:
 
            mindist[i] ≥ min_{j<i} D(i,j)
 
        Moreover, new nodes do not re-use one of the old indices, but they are
        given a new, unique index (SciPy convention: initial nodes are 0,…,N−1,
        new nodes are N,…,2N−2).
 
        Invalid nearest neighbors are not recognized by the fact that the stored
        distance is smaller than the actual distance, but the list active_nodes
        maintains a flag whether a node is inactive. If n_nghbr[i] points to an
        active node, the entries nn_distances[i] and n_nghbr[i] are valid,
        otherwise they must be recomputed.
      */
      idx1 = nn_distances.argmin();
      while (active_nodes.is_inactive(n_nghbr[idx1])) {
         // Recompute the minimum mindist[idx1] and n_nghbr[idx1].
         n_nghbr[idx1] = j = active_nodes.start;
         switch (method) {
         case METHOD_VECTOR_WARD:
            min = dist.ward_extended(idx1, j);
            for (j = active_nodes.succ[j]; j < idx1; j = active_nodes.succ[j]) {
               t_float tmp = dist.ward_extended(idx1, j);
               if (tmp < min) {
                  min = tmp;
                  n_nghbr[idx1] = j;
               }
            }
            break;
         default:
            min = dist.sqeuclidean_extended(idx1, j);
            for (j = active_nodes.succ[j]; j < idx1; j = active_nodes.succ[j]) {
               t_float const tmp = dist.sqeuclidean_extended(idx1, j);
               if (tmp < min) {
                  min = tmp;
                  n_nghbr[idx1] = j;
               }
            }
         }
         /* Update the heap with the new true minimum and search for the (possibly
            different) minimal entry. */
         nn_distances.update_geq(idx1, min);
         idx1 = nn_distances.argmin();
      }
 
      idx2 = n_nghbr[idx1];
      active_nodes.remove(idx1);
      active_nodes.remove(idx2);
 
      Z2.append(idx1, idx2, mindist[idx1]);
 
      if (i < 2 * N_1) {
         switch (method) {
         case METHOD_VECTOR_WARD:
         case METHOD_VECTOR_CENTROID: dist.merge(idx1, idx2, i); break;
 
         case METHOD_VECTOR_MEDIAN: dist.merge_weighted(idx1, idx2, i); break;
 
         default: throw std::runtime_error(std::string("Invalid method."));
         }
 
         n_nghbr[i] = active_nodes.start;
         if (method == METHOD_VECTOR_WARD) {
            /*
              Ward linkage.
 
              Shorter and longer distances can occur, not smaller than min(d1,d2)
              but maybe bigger than max(d1,d2).
            */
            min = dist.ward_extended(active_nodes.start, i);
            for (j = active_nodes.succ[active_nodes.start]; j < i; j = active_nodes.succ[j]) {
               t_float tmp = dist.ward_extended(j, i);
               if (tmp < min) {
                  min = tmp;
                  n_nghbr[i] = j;
               }
            }
         } else {
            /*
              Centroid and median linkage.
 
              Shorter and longer distances can occur, not bigger than max(d1,d2)
              but maybe smaller than min(d1,d2).
            */
            min = dist.sqeuclidean_extended(active_nodes.start, i);
            for (j = active_nodes.succ[active_nodes.start]; j < i; j = active_nodes.succ[j]) {
               t_float tmp = dist.sqeuclidean_extended(j, i);
               if (tmp < min) {
                  min = tmp;
                  n_nghbr[i] = j;
               }
            }
         }
         if (idx2 < active_nodes.start) {
            nn_distances.remove(active_nodes.start);
         } else {
            nn_distances.remove(idx2);
         }
         nn_distances.replace(idx1, i, min);
      }
   }
}
 
#if HAVE_VISIBILITY
#pragma GCC visibility pop
#endif