[SYSTEMDS-3166] Refactor Isolation Forest builtin and add edge-case tests

scripts/builtin/outlierByIsolationForestApply.dml

@@ -0,0 +1,260 @@
+#-------------------------------------------------------------
+#
+# 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.
+#
+#-------------------------------------------------------------
+
+# Builtin function that calculates the anomaly score as described in [Liu2008]
+# for a set of samples `X` based on an iForest model.
+#
+# [Liu2008]:
+#   Liu, F. T., Ting, K. M., & Zhou, Z. H.
+#   (2008, December).
+#   Isolation forest.
+#   In 2008 eighth ieee international conference on data mining (pp. 413-422).
+#   IEEE.
+#
+# .. code-block:: python
+#
+#   >>> import numpy as np
+#   >>> from systemds.context import SystemDSContext
+#   >>> from systemds.operator.algorithm import outlierByIsolationForest, outlierByIsolationForestApply
+#   >>> with SystemDSContext() as sds:
+#   ...     # Create training data: 20 points clustered near origin
+#   ...     X_train = sds.from_numpy(np.array([
+#   ...         [0.0, 0.0], [0.1, 0.1], [0.2, 0.2], [0.3, 0.3], [0.4, 0.4],
+#   ...         [0.5, 0.5], [0.6, 0.6], [0.7, 0.7], [0.8, 0.8], [0.9, 0.9],
+#   ...         [1.0, 1.0], [1.1, 1.1], [1.2, 1.2], [1.3, 1.3], [1.4, 1.4],
+#   ...         [1.5, 1.5], [1.6, 1.6], [1.7, 1.7], [1.8, 1.8], [1.9, 1.9]
+#   ...     ]))
+#   ...     model = outlierByIsolationForest(X_train, n_trees=100, subsampling_size=10, seed=42)
+#   ...     X_test = sds.from_numpy(np.array([[1.0, 1.0], [100.0, 100.0]]))
+#   ...     scores = outlierByIsolationForestApply(model, X_test).compute()
+#   ...     print(scores.shape)
+#   ...     print(scores[1, 0] > scores[0, 0])
+#   ...     print(scores[1, 0] > 0.5)
+#   (2, 1)
+#   True
+#   True
+#
+#
+# INPUT:
+# ---------------------------------------------------------------------------------------------
+# iForestModel  The trained iForest model as returned by outlierByIsolationForest
+# X             Samples to calculate the anomaly score for. X must contain every feature referenced
+#               by an internal model node
+# ---------------------------------------------------------------------------------------------
+#
+# OUTPUT:
+# ---------------------------------------------------------------------------------------------
+# anomaly_scores   Column vector of anomaly scores corresponding to the samples in X.
+#                  Samples with an anomaly score > 0.5 are generally considered to be outliers
+# ---------------------------------------------------------------------------------------------
+
+s_outlierByIsolationForestApply = function(List[Unknown] iForestModel, Matrix[Double] X)
+  return(Matrix[Double] anomaly_scores)
+{
+  anomaly_scores = m_outlierByIsolationForestApply(iForestModel, X)
+}
+
+m_outlierByIsolationForestApply = function(List[Unknown] iForestModel, Matrix[Double] X)
+  return(Matrix[Double] anomaly_scores)
+{
+  if (nrow(X) < 1)
+    stop("outlierByIsolationForestApply: X must contain at least one row.")
+
+  M = as.matrix(iForestModel["model"])
+  subsampling_size = as.integer(as.scalar(iForestModel["subsampling_size"]))
+
+  if (subsampling_size <= 1)
+    stop("outlierByIsolationForestApply: model subsampling_size must be greater than 1.")
+  if (nrow(M) < 1)
+    stop("outlierByIsolationForestApply: model must contain at least one tree.")
+
+  height_limit = ceil(log(subsampling_size, 2))
+  tree_size = 2 * (2^(height_limit + 1) - 1)
+
+  if (ncol(M) != tree_size)
+    stop("outlierByIsolationForestApply: model has an invalid number of columns.")
+
+  anomaly_scores = matrix(0, rows = nrow(X), cols = 1)
+
+  for (i_x in 1:nrow(X))
+    anomaly_scores[i_x, 1] = m_score(M, X[i_x, ], subsampling_size)
+}
+
+# Calculates the PathLength as defined in [Liu2008] based on a sample x
+#
+# INPUT PARAMETERS:
+# ---------------------------------------------------------------------------------------------
+# NAME          TYPE             DEFAULT      MEANING
+# ---------------------------------------------------------------------------------------------
+# M             Matrix[Double]                The linearized iTree model
+# x             Matrix[Double]                The sample to calculate the PathLength
+#
+# ---------------------------------------------------------------------------------------------
+# OUTPUT PARAMETERS:
+# ---------------------------------------------------------------------------------------------
+# PathLength  The PathLength for the sample
+# ---------------------------------------------------------------------------------------------
+m_PathLength = function(Matrix[Double] M, Matrix[Double] x)
+  return(Double PathLength)
+{
+  [nrEdgesTraversed, externalNodeSize] = s_traverseITree(M, x)
+
+  if (externalNodeSize <= 1) {
+    PathLength = nrEdgesTraversed
+  }
+  else {
+    PathLength = nrEdgesTraversed + s_cn(externalNodeSize)
+  }
+}
+
+
+# Traverses an iTree based on a sample x
+#
+# INPUT PARAMETERS:
+# ---------------------------------------------------------------------------------------------
+# NAME          TYPE             DEFAULT      MEANING
+# ---------------------------------------------------------------------------------------------
+# M             Matrix[Double]                The linearized iTree model to traverse
+# x             Matrix[Double]                The sample to traverse the iTree with
+#
+# ---------------------------------------------------------------------------------------------
+# OUTPUT PARAMETERS:
+# ---------------------------------------------------------------------------------------------
+# nrEdgesTraversed         The number of edges traversed until an external node was reached
+# externalNodeSize   The size of the external node assigned during training
+# ---------------------------------------------------------------------------------------------
+
+s_traverseITree = function(Matrix[Double] M, Matrix[Double] x)
+  return(Integer nrEdgesTraversed, Integer externalNodeSize)
+{
+  if (nrow(x) != 1)
+    stop("s_traverseITree: x must be a single row.")
+
+  nrEdgesTraversed = 0
+  externalNodeSize = 0
+  is_external_node = FALSE
+  node_id = 1
+
+  while (!is_external_node) {
+    node_start_idx = node_id * 2 - 1
+
+    if (node_start_idx + 1 > ncol(M))
+      stop("s_traverseITree: invalid iTree model; node index is out of bounds.")
+
+    split_feature = as.integer(as.scalar(M[1, node_start_idx]))
+    node_value = as.scalar(M[1, node_start_idx + 1])
+
+    if (split_feature > 0) {
+      if (split_feature > ncol(x))
+        stop("s_traverseITree: model split feature exceeds the input width.")
+
+      nrEdgesTraversed = nrEdgesTraversed + 1
+      x_val = as.scalar(x[1, split_feature])
+
+      # Training uses the same < / >= partition at every internal node.
+      if (x_val < node_value)
+        node_id = node_id * 2
+      else
+        node_id = node_id * 2 + 1
+    }
+    else if (split_feature == 0) {
+      if (node_value < 1)
+        stop("s_traverseITree: invalid iTree model; external-node size must be positive.")
+
+      externalNodeSize = as.integer(node_value)
+      is_external_node = TRUE
+    }
+    else {
+      stop("s_traverseITree: invalid iTree model; reached a placeholder node.")
+    }
+  }
+}
+
+
+# This function gives the average path length of unsuccessful search in BST `c(n)`
+# for `n` nodes as given in [Liu2008]. This function is used to normalize the path length
+#
+# INPUT PARAMETERS:
+# ---------------------------------------------------------------------------------------------
+# NAME          TYPE             DEFAULT      MEANING
+# ---------------------------------------------------------------------------------------------
+# n             Int                           Number of samples in the external node for which c(n)
+#                                             should be calculated
+# ---------------------------------------------------------------------------------------------
+# OUTPUT PARAMETERS:
+# ---------------------------------------------------------------------------------------------
+# cn   Value for c(n)
+# ---------------------------------------------------------------------------------------------
+s_cn = function(Integer n)
+  return(Double cn)
+{
+  if (n <= 1)
+    stop("s_cn: n must be greater than 1.")
+
+  # The logarithmic approximation has noticeable error for small n, so H(n-1)
+  # is evaluated directly below 1000.
+  if (n < 1000) {
+    indices = seq(1, n - 1)
+    H_nminus1 = sum(1 / indices)
+  }
+  else {
+    # Euler–Mascheroni's constant
+    eulergamma = 0.57721566490153
+    H_nminus1 = log(n - 1) + eulergamma
+  }
+
+  cn = 2 * H_nminus1 - 2 * (n - 1) / n
+}
+
+# Scores a sample `x` according to the score function `s(x, n)` described in [Liu2008].
+#
+# INPUT PARAMETERS:
+# ---------------------------------------------------------------------------------------------
+# NAME          TYPE             DEFAULT      MEANING
+# ---------------------------------------------------------------------------------------------
+# M             Matrix[Double]                 iForest model used to score
+# x             Matrix[Double]                 Sample to be scored
+# n             Int                            Subsample size the iTrees were built from
+# ---------------------------------------------------------------------------------------------
+# OUTPUT PARAMETERS:
+# ---------------------------------------------------------------------------------------------
+# score   The anomaly score for x
+# ---------------------------------------------------------------------------------------------
+m_score = function(Matrix[Double] M, Matrix[Double] x, Integer n)
+  return(Double score)
+{
+  if (n <= 1)
+    stop("m_score: n must be greater than 1.")
+  if (nrow(x) != 1)
+    stop("m_score: x must be a single row.")
+  if (nrow(M) < 1)
+    stop("m_score: model must contain at least one tree.")
+
+  # Only the mean path length is required, so no intermediate vector is built.
+  path_length_sum = 0.0
+
+  for (i_iTree in 1:nrow(M))
+    path_length_sum = path_length_sum + m_PathLength(M[i_iTree, ], x)
+
+  avg_path_length = path_length_sum / nrow(M)
+
+  score = 2^-(avg_path_length / s_cn(n))
+}