text: //
  // C++ standalone verion of fastcluster by Daniel Müllner
  //
  // Copyright: Christoph Dalitz, 2018
  //            Daniel Müllner, 2011
  // License:   BSD style license
  //            (see the file LICENSE for details)
  //
  
  
  #include <vector>
  #include <algorithm>
  #include <cmath>
  
  
  extern "C" {
  #include "fastcluster.h"
  }
  
  // Code by Daniel Müllner
  // workaround to make it usable as a standalone version (without R)
  bool fc_isnan(double x) { return false; }
  #include "fastcluster_dm.cpp"
  #include "fastcluster_R_dm.cpp"
  
  extern "C" {
  //
  // Assigns cluster labels (0, ..., nclust-1) to the n points such
  // that the cluster result is split into nclust clusters.
  //
  // Input arguments:
  //   n      = number of observables
  //   merge  = clustering result in R format
  //   nclust = number of clusters
  // Output arguments:
  //   labels = allocated integer array of size n for result
  //
    void cutree_k(int n, const int* merge, int nclust, int* labels) {
  
      int k,m1,m2,j,l;
  
      if (nclust > n || nclust < 2) {
        for (j=0; j<n; j++) labels[j] = 0;
        return;
      }
  
      // assign to each observable the number of its last merge step
      // beware: indices of observables in merge start at 1 (R convention)
      std::vector<int> last_merge(n, 0);
      for (k=1; k<=(n-nclust); k++) {
        // (m1,m2) = merge[k,]
        m1 = merge[k-1];
        m2 = merge[n-1+k-1];
        if (m1 < 0 && m2 < 0) { // both single observables
          last_merge[-m1-1] = last_merge[-m2-1] = k;
        }
        else if (m1 < 0 || m2 < 0) { // one is a cluster
          if(m1 < 0) { j = -m1; m1 = m2; } else j = -m2;
          // merging single observable and cluster
          for(l = 0; l < n; l++)
            if (last_merge[l] == m1)
              last_merge[l] = k;
          last_merge[j-1] = k;
        }
        else { // both cluster
          for(l=0; l < n; l++) {
            if( last_merge[l] == m1 || last_merge[l] == m2 )
              last_merge[l] = k;
          }
        }
      }
  
      // assign cluster labels
      int label = 0;
      std::vector<int> z(n,-1);
      for (j=0; j<n; j++) {
        if (last_merge[j] == 0) { // still singleton
          labels[j] = label++;
        } else {
          if (z[last_merge[j]] < 0) {
            z[last_merge[j]] = label++;
          }
          labels[j] = z[last_merge[j]];
        }
      }
    }
  
    //
    // Assigns cluster labels (0, ..., nclust-1) to the n points such
    // that the hierarchical clustering is stopped when cluster distance >= cdist
    //
    // Input arguments:
    //   n      = number of observables
    //   merge  = clustering result in R format
    //   height = cluster distance at each merge step
    //   cdist  = cutoff cluster distance
    // Output arguments:
    //   labels = allocated integer array of size n for result
    //
    void cutree_cdist(int n, const int* merge, double* height, double cdist, int* labels) {
  
      int k;
  
      for (k=0; k<(n-1); k++) {
        if (height[k] >= cdist) {
          break;
        }
      }
      cutree_k(n, merge, n-k, labels);
    }
  
  
    //
    // Hierarchical clustering with one of Daniel Muellner's fast algorithms
    //
    // Input arguments:
    //   n       = number of observables
    //   distmat = condensed distance matrix, i.e. an n*(n-1)/2 array representing
    //             the upper triangle (without diagonal elements) of the distance
    //             matrix, e.g. for n=4:
    //               d00 d01 d02 d03
    //               d10 d11 d12 d13   ->  d01 d02 d03 d12 d13 d23
    //               d20 d21 d22 d23
    //               d30 d31 d32 d33
    //   method  = cluster metric (see enum method_code)
    // Output arguments:
    //   merge   = allocated (n-1)x2 matrix (2*(n-1) array) for storing result.
    //             Result follows R hclust convention:
    //              - observabe indices start with one
    //              - merge[i][] contains the merged nodes in step i
    //              - merge[i][j] is negative when the node is an atom
    //   height  = allocated (n-1) array with distances at each merge step
    // Return code:
    //   0 = ok
    //   1 = invalid method
    //
    int hclust_fast(int n, double* distmat, int method, int* merge, double* height) {
  
      // call appropriate culstering function
      cluster_result Z2(n-1);
      if (method == HCLUST_METHOD_SINGLE) {
        // single link
        MST_linkage_core(n, distmat, Z2);
      }
      else if (method == HCLUST_METHOD_COMPLETE) {
        // complete link
        NN_chain_core<METHOD_METR_COMPLETE, t_float>(n, distmat, NULL, Z2);
      }
      else if (method == HCLUST_METHOD_AVERAGE) {
        // best average distance
        double* members = new double[n];
        for (int i=0; i<n; i++) members[i] = 1;
        NN_chain_core<METHOD_METR_AVERAGE, t_float>(n, distmat, members, Z2);
        delete[] members;
      }
      else if (method == HCLUST_METHOD_MEDIAN) {
        // best median distance (beware: O(n^3))
        generic_linkage<METHOD_METR_MEDIAN, t_float>(n, distmat, NULL, Z2);
      }
      else if (method == HCLUST_METHOD_CENTROID) {
        // best centroid distance (beware: O(n^3))
        double* members = new double[n];
        for (int i=0; i<n; i++) members[i] = 1;
        generic_linkage<METHOD_METR_CENTROID, t_float>(n, distmat, members, Z2);
        delete[] members;
      }
      else {
        return 1;
      }
  
      int* order = new int[n];
      if (method == HCLUST_METHOD_MEDIAN || method == HCLUST_METHOD_CENTROID) {
        generate_R_dendrogram<true>(merge, height, order, Z2, n);
      } else {
        generate_R_dendrogram<false>(merge, height, order, Z2, n);
      }
      delete[] order; // only needed for visualization
  
      return 0;
    }
  
  
    // Build condensed distance matrix
    // Input arguments:
    //   n  = number of observables
    //   m  = dimension of observable
    // Output arguments:
    //   out = allocated integer array of size n * (n - 1) / 2 for result
    void hclust_pdist(int n, int m, double* pts, double* out) {
      int ii = 0;
      for (int i = 0; i < n; i++) {
        for (int j = i + 1; j < n; j++) {
          // Compute euclidian distance
          double d = 0;
          for (int k = 0; k < m; k ++) {
            double error = pts[i * m + k] - pts[j * m + k];
            d += (error * error);
          }
          out[ii] = d;//sqrt(d);
          ii++;
        }
      }
    }
  
    void cluster_points_centroid(int n, int m, double* pts, double dist, int* idx) {
      double* pdist = new double[n * (n - 1) / 2];
      int* merge = new int[2 * (n - 1)];
      double* height = new double[n - 1];
  
      hclust_pdist(n, m, pts, pdist);
      hclust_fast(n, pdist, HCLUST_METHOD_CENTROID, merge, height);
      cutree_cdist(n, merge, height, dist, idx);
  
      delete[] pdist;
      delete[] merge;
      delete[] height;
    }
  }
  