[DRAFT] [SYSTEMDS-3858] New HDBSCAN Builtin Function

src/test/scripts/functions/builtin/hdbscan.dml

@@ -0,0 +1,20 @@
+#-------------------------------------------------------------
+#
+# Licensed to the Apache Software Foundation (ASF) under one
+# or more contributor license agreements.  See the NOTICE file
+# distributed with this work for additional information
+# regarding copyright ownership.  The ASF licenses this file
+# to you under the Apache License, Version 2.0 (the
+# "License"); you may not use this file except in compliance
+# with the License.  You may obtain a copy of the License at
+#
+#   http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing,
+# software distributed under the License is distributed on an
+# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+# KIND, either express or implied.  See the License for the
+# specific language governing permissions and limitations
+# under the License.
+#
+#-------------------------------------------------------------

test_build_hierarchy.dml

@@ -0,0 +1,121 @@
+source("scripts/builtin/hdbscan.dml") as hdb
+
+#        4
+#      / | \
+#     /  |  \
+#    /  (2)  (5)
+#    |   |    \
+#    |   |     \
+#   (2)  1-(3)--3
+#    |   |    /
+#    \   |   (4)
+#     \  (1) /
+#      \ | /
+#       \|/
+#        2
+
+distances = matrix(0, rows=4, cols=4)
+distances[1,2] = 1 
+distances[2,1] = 1
+
+distances[1,3] = 3 
+distances[3,1] = 3
+
+distances[1,4] = 2  
+distances[4,1] = 2
+
+distances[2,3] = 4 
+distances[3,2] = 4
+
+distances[2,4] = 2
+distances[4,2] = 2
+
+distances[3,4] = 5 
+distances[4,3] = 5
+
+[edges, weights] = hdb::buildMST(distances, 4)
+
+[hierarchy, sizes] = hdb::buildHierarchy(edges, weights, 4)
+
+print("Hierarchy (format: [cluster1, cluster2, merge_distance]):")
+print(toString(hierarchy))
+
+# Should have n-1 merge operations for n nodes
+num_merges = nrow(hierarchy)
+print("Number of merges: " + num_merges + " (should be 3)")
+test1 = (num_merges == 3)
+
+# Merge distances should be in ascending order (or equal)
+# Because we process edges from low weight to high weight
+dist1 = as.scalar(hierarchy[1,3])
+dist2 = as.scalar(hierarchy[2,3])
+dist3 = as.scalar(hierarchy[3,3])
+print("\nMerge distances: [" + dist1 + ", " + dist2 + ", " + dist3 + "]" + " (Should be in ascending order)")
+test2 = (dist1 <= dist2) & (dist2 <= dist3)
+
+# Cluster sizes should increase
+size1 = as.scalar(sizes[1])
+size2 = as.scalar(sizes[2])
+size3 = as.scalar(sizes[3])
+print("\nCluster sizes: [" + size1 + ", " + size2 + ", " + size3 + "]" + " (Should be increasing)")
+test3 = (size1 <= size2) & (size2 <= size3)
+
+# Final size should equal total number of nodes
+print("Final cluster size: " + size3 + " (should be 4)")
+test4 = (size3 == 4)
+
+# First merge should be size 2
+print("First merge size: " + size1 + " (should be 2)")
+test5 = (size1 == 2)
+
+# New classic-linkage checks
+n = 4
+hasInternal = (sum(hierarchy[,1] > n) + sum(hierarchy[,2] > n)) > 0
+print("Has internal node ids (>n): " + hasInternal + " (should be true)")
+test6 = hasInternal
+
+# Check that child ids are within valid range
+maxChild1 = max(hierarchy[,1])
+maxChild2 = max(hierarchy[,2])
+maxChild = maxChild1
+if(maxChild2 > maxChild) { maxChild = maxChild2 }
+
+print("Max child id: " + maxChild + " (should be <= " + (2*n-2) + ")")
+test7 = (maxChild <= (2*n-2))
+
+# Check that internal ids are “created in order”
+test8 = TRUE
+for(r in 1:(n-1)) {
+    child1 = as.scalar(hierarchy[r,1])
+    child2 = as.scalar(hierarchy[r,2])
+    newId = n + r
+    ok = (child1 < newId) & (child2 < newId)
+    test8 = test8 & ok
+}
+print("Children reference only existing nodes: " + test8 + " (should be true)")
+
+# Recompute node sizes from hierarchy and verify
+nodeSize = matrix(0, rows=2*n-1, cols=1)
+for(j in 1:n) { nodeSize[j,1] = 1 }
+
+test9 = TRUE
+for(r in 1:(n-1)) {
+    left = as.integer(as.scalar(hierarchy[r,1]))
+    right = as.integer(as.scalar(hierarchy[r,2]))
+    newId = n + r
+
+    expected = as.scalar(nodeSize[left,1]) + as.scalar(nodeSize[right,1])
+    nodeSize[newId,1] = expected
+
+    ok = (as.scalar(sizes[r,1]) == expected)
+    test9 = test9 & ok
+}
+print("sizes[r] equals sum of child sizes: " + test9 + " (should be true)")
+
+test_pass = test1 & test2 & test3 & test4 & test5 & test6 & test7 & test8 & test9
+
+if(test_pass) {
+    print("Passed")
+} else {
+    print("Failed")
+}

test_build_mst.dml

@@ -0,0 +1,46 @@
+source("scripts/builtin/hdbscan.dml") as hdb
+
+#        4
+#      / | \
+#     /  |  \
+#    /  (2)  (5)
+#    |   |    \
+#    |   |     \
+#   (2)  1-(3)--3
+#    |   |    /
+#    \   |   (4)
+#     \  (1) /
+#      \ | /
+#       \|/
+#        2
+
+distances = matrix(0, rows=4, cols=4)
+distances[1,2] = 1 
+distances[2,1] = 1
+
+distances[1,3] = 3 
+distances[3,1] = 3
+
+distances[1,4] = 2  
+distances[4,1] = 2
+
+distances[2,3] = 4 
+distances[3,2] = 4
+
+distances[2,4] = 2
+distances[4,2] = 2
+
+distances[3,4] = 5 
+distances[4,3] = 5
+
+[edges, weights] = hdb::buildMST(distances, 4)
+
+totalWeight = sum(weights)
+
+test_pass = (nrow(edges) == 3) & (totalWeight == 6)
+
+if(test_pass) {
+    print("Passed")
+} else {
+    print("Failed")
+}

test_cluster_model.dml

@@ -0,0 +1,129 @@
+source("scripts/builtin/hdbscan.dml") as hdb
+
+# 3 clear clusters
+# A: points around (0, 0)
+# B: points around (10, 10)
+# C: points around (20, 0)
+
+n = 12
+d = 2
+X = matrix(0, rows=n, cols=d)
+
+# A
+X[1,] = matrix("0.0 0.0", rows=1, cols=2)
+X[2,] = matrix("0.5 0.5", rows=1, cols=2)
+X[3,] = matrix("-0.5 0.5", rows=1, cols=2)
+X[4,] = matrix("0.0 -0.5", rows=1, cols=2)
+
+# B
+X[5,] = matrix("10.0 10.0", rows=1, cols=2)
+X[6,] = matrix("10.5 10.5", rows=1, cols=2)
+X[7,] = matrix("9.5 10.5", rows=1, cols=2)
+X[8,] = matrix("10.0 9.5", rows=1, cols=2)
+
+# C
+X[9,] = matrix("20.0 0.0", rows=1, cols=2)
+X[10,] = matrix("20.5 0.5", rows=1, cols=2)
+X[11,] = matrix("19.5 0.5", rows=1, cols=2)
+X[12,] = matrix("20.0 -0.5", rows=1, cols=2)
+
+print(toString(X))
+
+
+# get distances
+distances = hdb::dist(X)
+
+
+# get core distances
+minPts = 3
+coreDistances = matrix(0, rows=n, cols=1)
+for(i in 1:n) {
+    kthDist = hdb::computeKthSmallest(t(distances[i,]), minPts)
+    coreDistances[i] = kthDist
+}
+
+
+# get mutual reachability
+mutualReach = hdb::computeMutualReachability(distances, coreDistances)
+
+
+# get MST
+[edges, weights] = hdb::buildMST(mutualReach, n)
+
+
+# get hierarchy
+[hierarchy, sizes] = hdb::buildHierarchy(edges, weights, n)
+
+
+# get stable clusters
+minClSize = 3
+[labels, stabilities, clusterToNode] = hdb::extractStableClusters(hierarchy, weights, n, minClSize)
+expected_labels = matrix("1 1 1 1 2 2 2 2 3 3 3 3", rows=12, cols=1)
+labels_match = (min(labels == expected_labels) == 1)
+if (labels_match) {
+    print("Pass: labels match.")
+} else {
+    print("Fail: labels don't match.")
+}
+print("Cluster labels:")
+print(toString(labels))
+
+
+# get cluster model
+[centroids, clusterInfo] = hdb::buildClusterModel(X, labels, stabilities, clusterToNode)
+
+print("\nCentroids:")
+print(toString(centroids))
+
+print("\nInfo [size, stability]:")
+print(toString(clusterInfo))
+
+
+# check results
+numClusters = nrow(centroids)
+print("\nNumber of clusters found: " + numClusters)
+
+
+# should find 3 clusters
+test1 = (numClusters == 3)
+print("Found 3 clusters: " + test1)
+
+
+# each cluster should have 4 points
+allSizesFour = TRUE
+for(c in 1:numClusters) {
+    size = as.scalar(clusterInfo[c,1])
+    allSizesFour = allSizesFour & (size == 4)
+}
+print("All clusters have size 4: " + allSizesFour)
+
+
+test_A = min(sqrt(rowSums((centroids - matrix("0 0.125", 1, 2))^2))) < 0.001
+test_B = min(sqrt(rowSums((centroids - matrix("10 10.125", 1, 2))^2))) < 0.001
+test_C = min(sqrt(rowSums((centroids - matrix("20 0.125", 1, 2))^2))) < 0.001
+all_found = test_A & test_B & test_C
+print("Expected centroids near expected positions: " + all_found)
+
+
+# no noise (all assigned to clusters)
+numNoise = sum(labels == -1)
+test4 = (numNoise == 0)
+print("No noise points: " + test4)
+
+
+# stabilities are populated (not 0)
+test_stability = TRUE
+for(c in 1:numClusters) {
+    stab = as.scalar(clusterInfo[c,2])
+    test_stability = test_stability & (stab > 0)
+}
+print("Stabilities populated: " + test_stability)
+
+
+test_pass = test1 & allSizesFour & all_found & test4 & test_stability
+
+if(test_pass) {
+    print("\nAll tests passed")
+} else {
+    print("\nTests failed")
+}