scripts/builtin/outlierByIsolationForestApply.dml - systemds - Git at Google

 #-------------------------------------------------------------
 #
 # 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
 #
 # To compare implementations element by element, apply the same serialized iForestModel in each
 # implementation. Independently trained forests are generally different even when given the same
 # numeric seed because RNG and sampling semantics are implementation-specific.
 # ---------------------------------------------------------------------------------------------
 #
 # 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))
 }
	#-------------------------------------------------------------
	#
	# 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
	#
	# To compare implementations element by element, apply the same serialized iForestModel in each
	# implementation. Independently trained forests are generally different even when given the same
	# numeric seed because RNG and sampling semantics are implementation-specific.
	# ---------------------------------------------------------------------------------------------
	#
	# 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))
	}