dragonpilot beta3

date: 2023-07-26T22:20:36
commit: c6d842c412052be1985b63d683c63be9dcb2b0eb

third_party/cluster/fastcluster.cpp

@@ -0,0 +1,218 @@
+//
+// 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;
+  }
+}