diff --git a/scripts/builtin/outlierByIsolationForest.dml b/scripts/builtin/outlierByIsolationForest.dml new file mode 100644 index 00000000000..fa9d72b5352 --- /dev/null +++ b/scripts/builtin/outlierByIsolationForest.dml @@ -0,0 +1,521 @@ +#------------------------------------------------------------- +# +# 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 implements anomaly detection via isolation forest as described in +# [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. +# +# This function creates an iForest model for outlier detection. +# +# .. 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: +# --------------------------------------------------------------------------------------------- +# X Numerical feature matrix +# n_trees Number of iTrees to build +# subsampling_size Size of the subsample used to build each iTree. Values larger than nrow(X) +# are clamped to nrow(X) +# seed Seed for calls to `sample` and `rand`. -1 corresponds to a random seed +# --------------------------------------------------------------------------------------------- +# +# OUTPUT: +# --------------------------------------------------------------------------------------------- +# iForestModel The trained iForest model to be used in outlierByIsolationForestApply. +# The model is represented as a list with two entries: +# Entry 'model' (Matrix[Double]) - The iForest Model in linearized form (see m_iForest) +# Entry 'subsampling_size' (Integer) - The effective subsampling size used to build +# and score the model. +# --------------------------------------------------------------------------------------------- + +s_outlierByIsolationForest = function(Matrix[Double] X, Integer n_trees, Integer subsampling_size, Integer seed = -1) + return(List[Unknown] iForestModel) +{ + iForestModel = m_outlierByIsolationForest(X, n_trees, subsampling_size, seed) +} + +m_outlierByIsolationForest = function(Matrix[Double] X, Integer n_trees, Integer subsampling_size, Integer seed = -1) + return(List[Unknown] iForestModel) +{ + if (nrow(X) <= 1) + stop("outlierByIsolationForest: X must contain at least two rows.") + if (ncol(X) <= 0) + stop("outlierByIsolationForest: X must contain at least one feature column.") + if (n_trees <= 0) + stop("outlierByIsolationForest: n_trees must be positive.") + if (subsampling_size <= 1) + stop("outlierByIsolationForest: subsampling_size must be greater than 1.") + + # The effective size is part of the model because the same value normalizes + # path lengths when the forest is applied to new samples. + actual_subsampling_size = as.integer(min(subsampling_size, nrow(X))) + + M = m_iForest(X, n_trees, actual_subsampling_size, seed) + + iForestModel = list( + model = M, + subsampling_size = actual_subsampling_size + ) +} + +# This function implements isolation forest for numerical input features as +# described in [Liu2008]. +# +# The returned 'linearized' model is of type Matrix[Double] where each row +# corresponds to a linearized iTree (see m_iTree). Note that each tree in the +# model is padded with placeholder nodes such that each iTree has the same maximum depth. +# +# .. code-block:: +# +# For example, give a feature matrix with features [a,b,c,d] +# and the following iForest, M would look as follows: +# +# Level Tree 1 Tree 2 Node Depth +# ------------------------------------------------------------------- +# (L1) |d<5| |b<6| 0 +# / \ / \ +# (L2) 2 |a<7| 20 0 1 +# / \ +# (L3) 10 8 2 +# +# --> M := +# [[ 4, 5, 0, 2, 1, 7, -1, -1, -1, -1, 0, 10, 0, 8], (Tree 1) +# [ 2, 6, 0, 20, 0, 0, -1, -1, -1, -1, -1, -1, -1, -1]] (Tree 2) +# | (L1) | | (L2) | | (L3) | +# +# +# INPUT PARAMETERS: +# --------------------------------------------------------------------------------------------- +# NAME TYPE DEFAULT MEANING +# --------------------------------------------------------------------------------------------- +# X Matrix[Double] Numerical feature matrix +# n_trees Int Number of iTrees to build +# subsampling_size Int Size of the subsample to build iTrees with +# seed Int -1 Seed for calls to `sample` and `rand`. +# -1 corresponds to a random seed +# --------------------------------------------------------------------------------------------- +# OUTPUT PARAMETERS: +# --------------------------------------------------------------------------------------------- +# M Matrix containing the learned iForest in linearized form +# --------------------------------------------------------------------------------------------- + +m_iForest = function(Matrix[Double] X, Integer n_trees, Integer subsampling_size, Integer seed = -1) + return(Matrix[Double] M) +{ + if (n_trees <= 0) + stop("m_iForest: n_trees must be positive.") + if (subsampling_size <= 1 | subsampling_size > nrow(X)) + stop("m_iForest: subsampling_size must be in the interval [2, nrow(X)].") + + height_limit = ceil(log(subsampling_size, 2)) + tree_size = 2 * (2^(height_limit + 1) - 1) + + M = matrix(-1, rows = n_trees, cols = tree_size) + + # Each tree uses separate deterministic seeds for row sampling and splitting. + parfor (i_iTree in 1:n_trees, taskpartitioner="STATIC") { + tree_seed = ifelse(seed == -1, -1, seed + i_iTree) + X_subsample = s_sampleRows(X, subsampling_size, tree_seed) + + split_seed = ifelse(seed == -1, -1, seed + 1000003 * i_iTree) + M_tree = m_iTree(X_subsample, height_limit, split_seed) + + M[i_iTree, 1:ncol(M_tree)] = M_tree + } +} + +# This function implements isolation trees for numerical input features as +# described in [Liu2008]. +# +# The returned 'linearized' model is of type Matrix[Double] with exactly one row. +# Here, each node is represented by two consecutive entries in this row vector. +# Traversing the row vector from left to right corresponds to traversing the tree +# level-wise from top to bottom and left to right. If a node does not exist +# (e.g. because the parent node is already a leaf node), the node is still stored +# using placeholder values. +# Recall that for a binary tree with maximum depth `d`, the maximum number of nodes +# can be calculated by `2^(maximum depth + 1) - 1`. Hence, for a given maximum depth +# of an iTree, the row vector has exactly `2 * (2^(maximum depth + 1) - 1)` entries. +# +# There are three types of nodes that are represented in this model: +# - Internal Node +# A node that uses a "split feature" and corresponding "split value" to divide +# the data into two non-empty parts. The node is linearized in the following way: +# - Entry 1: Represents the index of the splitting feature in the feature matrix `X` +# - Entry 2: Represents splitting value +# +# - External Node +# A leaf node of the tree, It contains the "size" of the node. That is the +# number of remaining samples after splitting the feature matrix X by traversing +# the tree to this node. +# The node is linearized in the following way: +# - Entry 1: Always 0 - indicating an external node +# - Entry 2: The "size" of the node +# +# - Placeholder Node +# A node that is not present in the actual iTree and is used for "padding". +# Both entries are set to -1 +# +# .. code-block:: +# +# For example, give a feature matrix with features [a,b,c,d] +# and the following tree, M would look as follows: +# Level Tree Node Depth +# ------------------------------------------------- +# (L1) |d<5| 0 +# / \ +# (L2) 1 |a<7| 1 +# / \ +# (L3) 10 0 2 +# +# --> M := +# [[4, 5, 0, 1, 1, 7, -1, -1, -1, -1, 0, 10, 0, 0]] +# |(L1)| | (L2) | | (L3) | +# +# +# +# INPUT PARAMETERS: +# --------------------------------------------------------------------------------------------- +# NAME TYPE DEFAULT MEANING +# --------------------------------------------------------------------------------------------- +# X Matrix[Double] Numerical feature matrix +# max_depth Int Maximum depth of the learned tree where depth is the +# maximum number of edges from root to a leaf node +# seed Int -1 Seed for calls to `sample` and `rand`. +# -1 corresponds to a random seed +# --------------------------------------------------------------------------------------------- +# OUTPUT PARAMETERS: +# --------------------------------------------------------------------------------------------- +# M Matrix M containing the learned tree in linearized form +# --------------------------------------------------------------------------------------------- + +# A node can be split only if at least one feature has positive range. +s_hasSplittableFeature = function(Matrix[Double] X_node) + return(Boolean hasSplittableFeature) +{ + ranges = colMaxs(X_node) - colMins(X_node) + hasSplittableFeature = max(ranges) > 0 +} + +m_iTree = function(Matrix[Double] X, Integer max_depth, Integer seed = -1) + return(Matrix[Double] M) +{ + # These internal checks use warnings because m_iTree executes inside a parfor. + s_warning_assert_outlierByIsolationForest(max_depth > 0 & max_depth <= 32, + "m_iTree: max_depth must be in the interval [1, 32].") + s_warning_assert_outlierByIsolationForest(nrow(X) > 0, + "m_iTree: X must contain at least one row.") + + # Each node occupies two entries; the root has node ID 1 and depth 0. + M = matrix(-1, rows = 1, cols = 2 * (2^(max_depth + 1) - 1)) + + # Breadth-first construction replaces the recursive procedure from the paper. + # Queue entries contain the node ID and the data reaching that node. + node_queue = list(list(1, X)) + max_id = 1 + while (length(node_queue) > 0) { + [node_queue, queue_entry] = remove(node_queue, 1) + node = as.list(queue_entry) + node_id = as.scalar(node[1]) + X_node = as.matrix(node[2]) + + max_id = max(max_id, node_id) + + is_external_leaf = s_isExternalINode(X_node, node_id, max_depth) + if (is_external_leaf) { + M = s_addExternalINode(X_node, node_id, M) + } + else { + node_seed = ifelse(seed == -1, -1, seed + node_id) + + [split_feature, split_value] = s_drawSplitPoint(X_node, node_seed) + [left_id, X_left, right_id, X_right] = s_splitINode( + X_node, node_id, split_feature, split_value) + + # A defensive fallback keeps the serialized model traversable even if a + # random split happens to coincide with the minimum feature value. + if (nrow(X_left) == 0 | nrow(X_right) == 0) { + M = s_addExternalINode(X_node, node_id, M) + } + else { + M = s_addInternalINode(node_id, split_feature, split_value, M) + + node_queue = append(node_queue, list(left_id, X_left)) + node_queue = append(node_queue, list(right_id, X_right)) + } + } + } + + # Prune the individual tree; m_iForest pads all trees to a common width. + tree_depth = floor(log(max_id, 2)) + M = M[1, 1:(2 * (2^(tree_depth + 1) - 1))] +} + + +# Randomly draws a split point i.e. a feature and corresponding value to split a node by. +# +# INPUT PARAMETERS: +# --------------------------------------------------------------------------------------------- +# NAME TYPE DEFAULT MEANING +# --------------------------------------------------------------------------------------------- +# X Matrix[Double] Numerical feature matrix +# seed Int -1 Seed for calls to `sample` and `rand` +# -1 corresponds to a random seed +# +# --------------------------------------------------------------------------------------------- +# OUTPUT PARAMETERS: +# --------------------------------------------------------------------------------------------- +# split_feature Index of the feature used for splitting the node +# split_value Feature value used for splitting the node +# --------------------------------------------------------------------------------------------- + +# Draw only from features that can separate the current node. +s_drawSplitPoint = function(Matrix[Double] X, Integer seed = -1) + return(Integer split_feature, Double split_value) +{ + ranges = colMaxs(X) - colMins(X) + feature_ids = matrix(seq(1, ncol(X)), rows = 1, cols = ncol(X)) + + valid_features = removeEmpty( + target = feature_ids, + margin = "cols", + select = ranges > 0, + empty.return = FALSE + ) + + feature_idx = as.integer(as.scalar(sample(ncol(valid_features), 1, FALSE, seed))) + split_feature = as.integer(as.scalar(valid_features[1, feature_idx])) + + feature_min = min(X[, split_feature]) + feature_max = max(X[, split_feature]) + + split_value = as.scalar(rand( + rows = 1, + cols = 1, + min = feature_min, + max = feature_max, + seed = seed + )) +} + +# Adds an external (leaf) node to the linearized iTree model `M`. In the linearized form, +# each node is assigned two neighboring indices. For external nodes the value at the first +# index in M is always set to 0 while the value at the second index is set to the number of +# rows in the feature matrix corresponding to the node. +# +# INPUT PARAMETERS: +# --------------------------------------------------------------------------------------------- +# NAME TYPE DEFAULT MEANING +# --------------------------------------------------------------------------------------------- +# X_node Matrix[Double] Numerical feature matrix corresponding to the node +# node_id Int ID of the node +# M Matrix[Double] Linearized model to add the node to +# --------------------------------------------------------------------------------------------- +# OUTPUT PARAMETERS: +# --------------------------------------------------------------------------------------------- +# M The updated model +# --------------------------------------------------------------------------------------------- +s_addExternalINode = function(Matrix[Double] X_node, Integer node_id, Matrix[Double] M) + return(Matrix[Double] M) +{ + s_warning_assert_outlierByIsolationForest(node_id > 0, "s_addExternalINode: Requirement `node_id > 0` not satisfied!") + + node_start_index = 2 * node_id - 1 + M[, node_start_index] = 0 + M[, node_start_index + 1] = nrow(X_node) +} + +# Adds an internal node to the linearized iTree model `M`. In the linearized form, +# each node is assigned two neighboring indices. For internal nodes the value at the first +# index in M is set to index of the feature to split by while the value at the second index +# is set to the value to split the node by. +# +# INPUT PARAMETERS: +# --------------------------------------------------------------------------------------------- +# NAME TYPE DEFAULT MEANING +# --------------------------------------------------------------------------------------------- +# node_id Int ID of the node +# split_feature Int Index of the feature to split the node by +# split_value Double Value to split the node by +# M Matrix[Double] Linearized model to add the node to +# --------------------------------------------------------------------------------------------- +# OUTPUT PARAMETERS: +# --------------------------------------------------------------------------------------------- +# M The updated model +# --------------------------------------------------------------------------------------------- +s_addInternalINode = function(Integer node_id, Integer split_feature, Double split_value, Matrix[Double] M) + return(Matrix[Double] M) +{ + s_warning_assert_outlierByIsolationForest(node_id > 0, "s_addInternalINode: Requirement `node_id > 0` not satisfied!") + s_warning_assert_outlierByIsolationForest(split_feature > 0, "s_addInternalINode: Requirement `split_feature > 0` not satisfied!") + + node_start_index = 2 * node_id - 1 + M[, node_start_index] = split_feature + M[, node_start_index + 1] = split_value +} + +# Determines whether an iTree node is external based on its depth and data. +# +# INPUT PARAMETERS: +# --------------------------------------------------------------------------------------------- +# NAME TYPE DEFAULT MEANING +# --------------------------------------------------------------------------------------------- +# X_node Matrix[Double] Numerical feature matrix corresponding to the node +# node_id Int ID belonging to the node +# max_depth Int Maximum depth of the learned tree where depth is the +# maximum number of edges from root to a leaf node +# --------------------------------------------------------------------------------------------- +# OUTPUT PARAMETERS: +# --------------------------------------------------------------------------------------------- +# isExternalNode true if the node is an external (leaf) node, false otherwise. +# This is the case when a max depth is reached or the number of rows +# in the feature matrix corresponding to the node <= 1 +# --------------------------------------------------------------------------------------------- + +s_isExternalINode = function(Matrix[Double] X_node, Integer node_id, Integer max_depth) + return(Boolean isExternalNode) +{ + s_warning_assert_outlierByIsolationForest(max_depth > 0, + "s_isExternalINode: max_depth must be positive.") + s_warning_assert_outlierByIsolationForest(node_id > 0, + "s_isExternalINode: node_id must be positive.") + + node_depth = floor(log(node_id, 2)) + + isExternalNode = node_depth >= max_depth | nrow(X_node) <= 1 | !s_hasSplittableFeature(X_node) +} + + +# This function splits a node based on a given feature and value and returns the sub-matrices +# and IDs corresponding to the nodes resulting from the split. +# +# INPUT PARAMETERS: +# --------------------------------------------------------------------------------------------- +# NAME TYPE DEFAULT MEANING +# --------------------------------------------------------------------------------------------- +# X_node Matrix[Double] Numerical feature matrix corresponding +# node_id Int ID of the node to split +# split_feature Int Index of the feature to split the input matrix by +# split_value Double Value of the feature to split the input matrix by +# +# --------------------------------------------------------------------------------------------- +# OUTPUT PARAMETERS: +# --------------------------------------------------------------------------------------------- +# left_id ID of the resulting left node +# X_left Matrix corresponding to the left node resulting from the split with rows where +# value for feature `split_feature` < value `split_value` +# right_id ID of the resulting right node +# X_right Matrix corresponding to the right node resulting from the split with rows where +# value for feature `split_feature` >= value `split_value` +# --------------------------------------------------------------------------------------------- +s_splitINode = function(Matrix[Double] X_node, Integer node_id, Integer split_feature, Double split_value) + return(Integer left_id, Matrix[Double] X_left, Integer right_id, Matrix[Double] X_right) +{ + s_warning_assert_outlierByIsolationForest(nrow(X_node) > 0, "s_splitINode: Requirement `nrow(X_node) > 0` not satisfied!") + s_warning_assert_outlierByIsolationForest(node_id > 0, "s_splitINode: Requirement `nrow(X_node) > 0` not satisfied!") + s_warning_assert_outlierByIsolationForest(split_feature > 0, "s_splitINode: Requirement `split_feature > 0` not satisfied!") + + # Follow the original iTree convention: left is <, right is >=. + left_rows_mask = X_node[, split_feature] < split_value + + # In the linearized form of the iTree model, nodes need to be ordered by depth. + # Since iTrees are binary trees we can use 2*node_id/2*node_id+1 for left/right child ids + # to ensure that IDs are chosen accordingly. + left_id = 2 * node_id + X_left = removeEmpty(target=X_node, margin="rows", select=left_rows_mask, empty.return=FALSE) + + right_id = 2 * node_id + 1 + X_right = removeEmpty(target=X_node, margin="rows", select=!left_rows_mask, empty.return=FALSE) +} + +# Randomly samples `size` rows from a matrix X +# +# INPUT PARAMETERS: +# --------------------------------------------------------------------------------------------- +# NAME TYPE DEFAULT MEANING +# --------------------------------------------------------------------------------------------- +# X Matrix[Double] Matrix to sample rows from +# sample_size Int Number of rows to sample +# seed Int -1 Seed for calls to `sample` +# -1 corresponds to a random seed +# +# --------------------------------------------------------------------------------------------- +# OUTPUT PARAMETERS: +# --------------------------------------------------------------------------------------------- +# X_sampled Sampled rows from X +# --------------------------------------------------------------------------------------------- + +s_sampleRows = function(Matrix[Double] X, Integer size, Integer seed = -1) + return(Matrix[Double] X_extracted) +{ + s_warning_assert_outlierByIsolationForest(size > 0, + "s_sampleRows: size must be positive.") + s_warning_assert_outlierByIsolationForest(nrow(X) >= size, + "s_sampleRows: size must be <= nrow(X).") + + random_vector = rand(rows = nrow(X), cols = 1, seed = seed) + X_rand = cbind(X, random_vector) + + X_rand = order(target = X_rand, by = ncol(X_rand)) + + X_extracted = X_rand[1:size, 1:ncol(X)] +} + +# Prints a warning for internal invariant violations. Helpers called from the +# tree-building parfor cannot use `stop`; public inputs are validated before it. +# +# INPUT PARAMETERS: +# --------------------------------------------------------------------------------------------- +# NAME TYPE DEFAULT MEANING +# --------------------------------------------------------------------------------------------- +# assertion Boolean Assertion to check +# warning String Warning message to print if assertion is violated +# --------------------------------------------------------------------------------------------- +s_warning_assert_outlierByIsolationForest = function(Boolean assertion, String warning) +{ + if (!assertion) + print("outlierByIsolationForest: " + warning) +} diff --git a/scripts/builtin/outlierByIsolationForestApply.dml b/scripts/builtin/outlierByIsolationForestApply.dml new file mode 100644 index 00000000000..ba6c55f5c79 --- /dev/null +++ b/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)) +} diff --git a/scripts/staging/isolationForest/test/isolationForestTest.dml b/scripts/staging/isolationForest/test/isolationForestTest.dml index 9decfe6087d..1bae38aced1 100644 --- a/scripts/staging/isolationForest/test/isolationForestTest.dml +++ b/scripts/staging/isolationForest/test/isolationForestTest.dml @@ -702,6 +702,234 @@ parfor (i_run in 1:nr_runs) { print("\n") } +#------------------------------------------------------------- +# Isolation Forest - Edge Case Tests +# Following the structure of isolationForestTest.dml +#------------------------------------------------------------- + +print("===============================================================") +print("Isolation Forest - Edge Case Tests") +print("===============================================================") +print("") + +# Test data +X_100x5 = rand(rows=100, cols=5, seed=42) +X_zeros = matrix(0.0, rows=100, cols=5) +X_identical = matrix(5.0, rows=100, cols=5) + +print("===============================================================") +print("CATEGORY 1: Extreme Data Sizes") +print("===============================================================") +print("") + +# Test 1: Minimum dataset (2 rows) +print("Test 1: Minimum dataset (2 rows, 2 features)") +testname = "Minimum dataset (2x2)" +X_tiny = rand(rows=2, cols=2, seed=42) +model_tiny = iForest::outlierByIsolationForest(X=X_tiny, n_trees=10, subsampling_size=2) +scores_tiny = iForest::outlierByIsolationForestApply(iForestModel=model_tiny, X=X_tiny) +test_res = (nrow(scores_tiny) == 2) & (ncol(scores_tiny) == 1) +[test_cnt, fails] = record_test_result(testname, test_res, test_cnt, fails) +print("") + +# Test 2: Very small dataset (3 rows) +print("Test 2: Very small dataset (3 rows, 3 features)") +testname = "Very small dataset (3x3)" +X_mini = rand(rows=3, cols=3, seed=42) +model_mini = iForest::outlierByIsolationForest(X=X_mini, n_trees=5, subsampling_size=3) +scores_mini = iForest::outlierByIsolationForestApply(iForestModel=model_mini, X=X_mini) +test_res = (nrow(scores_mini) == 3) & (ncol(scores_mini) == 1) +[test_cnt, fails] = record_test_result(testname, test_res, test_cnt, fails) +print("") + +# Test 3: Large dataset +print("Test 3: Large dataset (5,000 rows, 10 features)") +testname = "Large dataset (5000x10)" +X_large = rand(rows=5000, cols=10, seed=42) +model_large = iForest::outlierByIsolationForest(X=X_large, n_trees=30, subsampling_size=256) +scores_large = iForest::outlierByIsolationForestApply(iForestModel=model_large, X=X_large[1:10,]) +test_res = (nrow(scores_large) == 10) +[test_cnt, fails] = record_test_result(testname, test_res, test_cnt, fails) +print("") + +print("===============================================================") +print("CATEGORY 2: Extreme Feature Counts") +print("===============================================================") +print("") + +# Test 4: Single feature +print("Test 4: Single feature dataset (100 rows, 1 feature)") +testname = "Single feature dataset" +X_one_feature = rand(rows=100, cols=1, seed=42) +model_one = iForest::outlierByIsolationForest(X=X_one_feature, n_trees=20, subsampling_size=50) +scores_one = iForest::outlierByIsolationForestApply(iForestModel=model_one, X=X_one_feature) +test_res = (nrow(scores_one) == 100) & (ncol(scores_one) == 1) +[test_cnt, fails] = record_test_result(testname, test_res, test_cnt, fails) +print("") + +# Test 5: High-dimensional dataset +print("Test 5: High-dimensional dataset (200 rows, 30 features)") +testname = "High-dimensional dataset (30 features)" +X_high_dim = rand(rows=200, cols=30, seed=42) +model_high = iForest::outlierByIsolationForest(X=X_high_dim, n_trees=20, subsampling_size=100) +scores_high = iForest::outlierByIsolationForestApply(iForestModel=model_high, X=X_high_dim[1:10,]) +test_res = (nrow(scores_high) == 10) +[test_cnt, fails] = record_test_result(testname, test_res, test_cnt, fails) +print("") + +print("===============================================================") +print("CATEGORY 3: Extreme Values") +print("===============================================================") +print("") + +# Test 6: All zeros +print("Test 6: All values are zero") +testname = "All zeros" +model_zeros = iForest::outlierByIsolationForest(X=X_zeros, n_trees=20, subsampling_size=50) +scores_zeros = iForest::outlierByIsolationForestApply(iForestModel=model_zeros, X=X_zeros) +mean_score = mean(scores_zeros) +# When all data is identical, algorithm can't split -> maximum path length +# -> LOW anomaly scores (all points are equally "normal") +test_res = (mean_score > 0.2) & (mean_score < 0.35) # Expect low scores for identical data +[test_cnt, fails] = record_test_result(testname, test_res, test_cnt, fails) +print(" Mean score: " + toString(mean_score) + " (low score expected - no variation to detect)") +print("") + +# Test 7: All identical values +print("Test 7: All identical values (all 5s)") +testname = "All identical values" +model_identical = iForest::outlierByIsolationForest(X=X_identical, n_trees=20, subsampling_size=50) +scores_identical = iForest::outlierByIsolationForestApply(iForestModel=model_identical, X=X_identical) +mean_score = mean(scores_identical) +# Same logic: identical data -> can't split -> maximum path -> low scores +test_res = (mean_score > 0.2) & (mean_score < 0.35) # Expect low scores for identical data +[test_cnt, fails] = record_test_result(testname, test_res, test_cnt, fails) +print(" Mean score: " + toString(mean_score) + " (low score expected - no variation to detect)") +print("") + +# Test 8: Negative values +print("Test 8: All negative values") +testname = "All negative values" +X_negative = rand(rows=100, cols=5, min=-100, max=-1, seed=42) +model_negative = iForest::outlierByIsolationForest(X=X_negative, n_trees=20, subsampling_size=50) +scores_negative = iForest::outlierByIsolationForestApply(iForestModel=model_negative, X=X_negative) +test_res = (nrow(scores_negative) == 100) +[test_cnt, fails] = record_test_result(testname, test_res, test_cnt, fails) +print("") + +# Test 9: Very large values +print("Test 9: Very large values (thousands)") +testname = "Very large values" +X_huge = rand(rows=100, cols=5, min=1000, max=10000, seed=42) +model_huge = iForest::outlierByIsolationForest(X=X_huge, n_trees=10, subsampling_size=50) +scores_huge = iForest::outlierByIsolationForestApply(iForestModel=model_huge, X=X_huge) +test_res = (nrow(scores_huge) == 100) +[test_cnt, fails] = record_test_result(testname, test_res, test_cnt, fails) +print("") + +# Test 10: Very small values +print("Test 10: Very small values (near zero)") +testname = "Very small values" +X_tiny_vals = rand(rows=100, cols=5, min=0.0001, max=0.001, seed=42) +model_tiny_vals = iForest::outlierByIsolationForest(X=X_tiny_vals, n_trees=10, subsampling_size=50) +scores_tiny_vals = iForest::outlierByIsolationForestApply(iForestModel=model_tiny_vals, X=X_tiny_vals) +test_res = (nrow(scores_tiny_vals) == 100) +[test_cnt, fails] = record_test_result(testname, test_res, test_cnt, fails) +print("") + +# Test 11: Mixed extreme values +print("Test 11: Mixed extreme values") +testname = "Mixed extreme values" +X_mixed = rand(rows=100, cols=5, min=-1000, max=1000, seed=42) +model_mixed = iForest::outlierByIsolationForest(X=X_mixed, n_trees=10, subsampling_size=50) +scores_mixed = iForest::outlierByIsolationForestApply(iForestModel=model_mixed, X=X_mixed) +test_res = (nrow(scores_mixed) == 100) +[test_cnt, fails] = record_test_result(testname, test_res, test_cnt, fails) +print("") + +print("===============================================================") +print("CATEGORY 4: Model Validation") +print("===============================================================") +print("") + +# NOTE: Seed reproducibility tests skipped due to integer overflow bug +# in seed generation (isolationForest.dml line 144) +# Bug: seeds = matrix(seq(1, n_trees), cols=n_trees, rows=1)*seed +# causes overflow when n_trees >= 10 and seed > 0 + +# Test 12: Model produces valid output +print("Test 12: Model produces valid anomaly scores") +testname = "Valid model output" +X_test = rand(rows=50, cols=5, seed=789) +# Use small n_trees and no seed to avoid bugs +model_test = iForest::outlierByIsolationForest(X=X_test, n_trees=5, subsampling_size=25) +scores_test = iForest::outlierByIsolationForestApply(iForestModel=model_test, X=X_test) +# Check we got scores for all rows +test_res = (nrow(scores_test) == 50) & (ncol(scores_test) == 1) +[test_cnt, fails] = record_test_result(testname, test_res, test_cnt, fails) +print(" Scores produced: " + toString(nrow(scores_test)) + " rows") +print("") + +# Test 13: Scores are in valid range +print("Test 13: Anomaly scores in valid range [0,1]") +testname = "Valid score range" +X_range = rand(rows=100, cols=5, seed=456) +model_range = iForest::outlierByIsolationForest(X=X_range, n_trees=5, subsampling_size=50) +scores_range = iForest::outlierByIsolationForestApply(iForestModel=model_range, X=X_range) +min_score = min(scores_range) +max_score = max(scores_range) +test_res = (min_score >= 0) & (max_score <= 1) +[test_cnt, fails] = record_test_result(testname, test_res, test_cnt, fails) +print(" Score range: [" + toString(min_score) + ", " + toString(max_score) + "]") +print("") + +print("===============================================================") +print("CATEGORY 5: Outlier Detection Scenarios") +print("===============================================================") +print("") + +# Test 14: Subtle outliers +print("Test 14: Subtle outliers (2 standard deviations)") +testname = "Subtle outliers" +X_normal = rand(rows=99, cols=5, pdf="normal", seed=42) +X_subtle_outlier = matrix(2.0, rows=1, cols=5) +X_with_subtle = rbind(X_normal, X_subtle_outlier) +model_subtle = iForest::outlierByIsolationForest(X=X_with_subtle, n_trees=30, subsampling_size=50) +scores_subtle = iForest::outlierByIsolationForestApply(iForestModel=model_subtle, X=X_with_subtle) +outlier_score = as.scalar(scores_subtle[100,1]) +normal_mean_score = mean(scores_subtle[1:99,]) +test_res = outlier_score > normal_mean_score +[test_cnt, fails] = record_test_result(testname, test_res, test_cnt, fails) +print(" Outlier score: " + toString(outlier_score) + ", Normal mean: " + toString(normal_mean_score)) +print("") + +# Test 15: Outlier in only one feature +print("Test 15: Outlier in only one feature") +testname = "One-dimensional outlier" +X_one_dim_outlier = rand(rows=100, cols=5, min=0, max=10, seed=42) +X_one_dim_outlier[1,1] = 100 +model_one_dim = iForest::outlierByIsolationForest(X=X_one_dim_outlier, n_trees=30, subsampling_size=50) +scores_one_dim = iForest::outlierByIsolationForestApply(iForestModel=model_one_dim, X=X_one_dim_outlier) +outlier_score = as.scalar(scores_one_dim[1,1]) +normal_mean = mean(scores_one_dim[2:100,]) +test_res = outlier_score > normal_mean +[test_cnt, fails] = record_test_result(testname, test_res, test_cnt, fails) +print(" Outlier score: " + toString(outlier_score) + ", Normal mean: " + toString(normal_mean)) +print("") + +print("===============================================================") +print("SUMMARY") +print("===============================================================") +succ_test_cnt = test_cnt - length(fails) +print(toString(succ_test_cnt) + "/" + toString(test_cnt) + " tests succeeded!") +if (length(fails) > 0) { + print("\nTests that failed:") + print(toString(fails)) +} + +print("===============================================================") + + print("===============================================================") print("TESTING FINISHED!") \ No newline at end of file diff --git a/src/main/java/org/apache/sysds/common/Builtins.java b/src/main/java/org/apache/sysds/common/Builtins.java index f5719641df7..2a9de5965df 100644 --- a/src/main/java/org/apache/sysds/common/Builtins.java +++ b/src/main/java/org/apache/sysds/common/Builtins.java @@ -267,6 +267,8 @@ public enum Builtins { OUTLIER_ARIMA("outlierByArima",true), OUTLIER_IQR("outlierByIQR", true), OUTLIER_IQR_APPLY("outlierByIQRApply", true), + OUTLIER_ISOLATION_FOREST("outlierByIsolationForest", true), + OUTLIER_ISOLATION_FOREST_APPLY("outlierByIsolationForestApply", true), OUTLIER_SD("outlierBySd", true), OUTLIER_SD_APPLY("outlierBySdApply", true), PAGERANK("pageRank", true), diff --git a/src/main/python/docs/source/api/operator/algorithms/outlierByIsolationForest.rst b/src/main/python/docs/source/api/operator/algorithms/outlierByIsolationForest.rst new file mode 100644 index 00000000000..4cc27950c73 --- /dev/null +++ b/src/main/python/docs/source/api/operator/algorithms/outlierByIsolationForest.rst @@ -0,0 +1,25 @@ +.. ------------------------------------------------------------- +.. +.. 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. +.. +.. ------------------------------------------------------------ + +outlierByIsolationForest +======================== + +.. autofunction:: systemds.operator.algorithm.outlierByIsolationForest \ No newline at end of file diff --git a/src/main/python/docs/source/api/operator/algorithms/outlierByIsolationForestApply.rst b/src/main/python/docs/source/api/operator/algorithms/outlierByIsolationForestApply.rst new file mode 100644 index 00000000000..aff908f4785 --- /dev/null +++ b/src/main/python/docs/source/api/operator/algorithms/outlierByIsolationForestApply.rst @@ -0,0 +1,25 @@ +.. ------------------------------------------------------------- +.. +.. 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. +.. +.. ------------------------------------------------------------ + +outlierByIsolationForestApply +============================= + +.. autofunction:: systemds.operator.algorithm.outlierByIsolationForestApply \ No newline at end of file diff --git a/src/main/python/systemds/operator/algorithm/__init__.py b/src/main/python/systemds/operator/algorithm/__init__.py index 9996c017f6b..94c2864b65a 100644 --- a/src/main/python/systemds/operator/algorithm/__init__.py +++ b/src/main/python/systemds/operator/algorithm/__init__.py @@ -159,6 +159,8 @@ from .builtin.outlierByArima import outlierByArima from .builtin.outlierByIQR import outlierByIQR from .builtin.outlierByIQRApply import outlierByIQRApply +from .builtin.outlierByIsolationForest import outlierByIsolationForest +from .builtin.outlierByIsolationForestApply import outlierByIsolationForestApply from .builtin.outlierBySd import outlierBySd from .builtin.outlierBySdApply import outlierBySdApply from .builtin.pageRank import pageRank @@ -359,6 +361,8 @@ 'outlierByArima', 'outlierByIQR', 'outlierByIQRApply', + 'outlierByIsolationForest', + 'outlierByIsolationForestApply', 'outlierBySd', 'outlierBySdApply', 'pageRank', diff --git a/src/main/python/systemds/operator/algorithm/builtin/outlierByIsolationForest.py b/src/main/python/systemds/operator/algorithm/builtin/outlierByIsolationForest.py new file mode 100644 index 00000000000..d6adf720ccf --- /dev/null +++ b/src/main/python/systemds/operator/algorithm/builtin/outlierByIsolationForest.py @@ -0,0 +1,86 @@ +# ------------------------------------------------------------- +# +# 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. +# +# ------------------------------------------------------------- + +# Autogenerated By : src/main/python/generator/generator.py +# Autogenerated From : scripts/builtin/outlierByIsolationForest.dml + +from typing import Dict, Iterable + +from systemds.operator import OperationNode, Matrix, Frame, List, MultiReturn, Scalar +from systemds.utils.consts import VALID_INPUT_TYPES + + +def outlierByIsolationForest(X: Matrix, + n_trees: int, + subsampling_size: int, + **kwargs: Dict[str, VALID_INPUT_TYPES]): + """ + Builtin function that implements anomaly detection via isolation forest as described in + [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. + + This function creates an iForest model for outlier detection. + + .. 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 + + + + + :param X: Numerical feature matrix + :param n_trees: Number of iTrees to build + :param subsampling_size: Size of the subsample to build iTrees with + :param seed: Seed for calls to `sample` and `rand`. -1 corresponds to a random seed + :return: The trained iForest model to be used in outlierByIsolationForestApply. + The model is represented as a list with two entries: + Entry 'model' (Matrix[Double]) - The iForest Model in linearized form (see m_iForest) + Entry 'subsampling_size' (Double) - The subsampling size used to build the model. + """ + + params_dict = {'X': X, 'n_trees': n_trees, 'subsampling_size': subsampling_size} + params_dict.update(kwargs) + return Matrix(X.sds_context, + 'outlierByIsolationForest', + named_input_nodes=params_dict) diff --git a/src/main/python/systemds/operator/algorithm/builtin/outlierByIsolationForestApply.py b/src/main/python/systemds/operator/algorithm/builtin/outlierByIsolationForestApply.py new file mode 100644 index 00000000000..2ecb612a3e5 --- /dev/null +++ b/src/main/python/systemds/operator/algorithm/builtin/outlierByIsolationForestApply.py @@ -0,0 +1,79 @@ +# ------------------------------------------------------------- +# +# 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. +# +# ------------------------------------------------------------- + +# Autogenerated By : src/main/python/generator/generator.py +# Autogenerated From : scripts/builtin/outlierByIsolationForestApply.dml + +from typing import Dict, Iterable + +from systemds.operator import OperationNode, Matrix, Frame, List, MultiReturn, Scalar +from systemds.utils.consts import VALID_INPUT_TYPES + + +def outlierByIsolationForestApply(iForestModel: List, + X: Matrix): + """ + 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 + + + + + :param iForestModel: The trained iForest model as returned by outlierByIsolationForest + :param X: Samples to calculate the anomaly score for + :return: Column vector of anomaly scores corresponding to the samples in X. + Samples with an anomaly score > 0.5 are generally considered to be outliers + """ + + params_dict = {'iForestModel': iForestModel, 'X': X} + return Matrix(iForestModel.sds_context, + 'outlierByIsolationForestApply', + named_input_nodes=params_dict) diff --git a/src/main/python/tests/auto_tests/test_outlierByIsolationForest.py b/src/main/python/tests/auto_tests/test_outlierByIsolationForest.py new file mode 100644 index 00000000000..41e0f4b5f31 --- /dev/null +++ b/src/main/python/tests/auto_tests/test_outlierByIsolationForest.py @@ -0,0 +1,56 @@ +# ------------------------------------------------------------- +# +# 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. +# +# ------------------------------------------------------------- + +# Autogenerated By : src/main/python/generator/generator.py +import unittest, contextlib, io + + +class TestOUTLIERBYISOLATIONFOREST(unittest.TestCase): + def test_outlierByIsolationForest(self): + # Example test case provided in python the code block + buf = io.StringIO() + with contextlib.redirect_stdout(buf): + 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) + + expected = """(2, 1) +True +True""" + self.assertEqual(buf.getvalue().strip(), expected) + + +if __name__ == '__main__': + unittest.main() diff --git a/src/main/python/tests/auto_tests/test_outlierByIsolationForestApply.py b/src/main/python/tests/auto_tests/test_outlierByIsolationForestApply.py new file mode 100644 index 00000000000..e0a9abb1a90 --- /dev/null +++ b/src/main/python/tests/auto_tests/test_outlierByIsolationForestApply.py @@ -0,0 +1,56 @@ +# ------------------------------------------------------------- +# +# 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. +# +# ------------------------------------------------------------- + +# Autogenerated By : src/main/python/generator/generator.py +import unittest, contextlib, io + + +class TestOUTLIERBYISOLATIONFORESTAPPLY(unittest.TestCase): + def test_outlierByIsolationForestApply(self): + # Example test case provided in python the code block + buf = io.StringIO() + with contextlib.redirect_stdout(buf): + 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) + + expected = """(2, 1) +True +True""" + self.assertEqual(buf.getvalue().strip(), expected) + + +if __name__ == '__main__': + unittest.main() diff --git a/src/test/java/org/apache/sysds/test/functions/builtin/part2/BuiltinIsolationForestTest.java b/src/test/java/org/apache/sysds/test/functions/builtin/part2/BuiltinIsolationForestTest.java new file mode 100644 index 00000000000..1086927a688 --- /dev/null +++ b/src/test/java/org/apache/sysds/test/functions/builtin/part2/BuiltinIsolationForestTest.java @@ -0,0 +1,418 @@ +/* + * 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. + */ + +package org.apache.sysds.test.functions.builtin.part2; + +import java.util.Arrays; +import java.util.Collection; +import java.util.HashMap; +import java.util.Locale; + +import org.apache.sysds.common.Types.ExecMode; +import org.apache.sysds.runtime.matrix.data.MatrixValue.CellIndex; +import org.apache.sysds.runtime.meta.MatrixCharacteristics; +import org.apache.sysds.test.AutomatedTestBase; +import org.apache.sysds.test.TestConfiguration; +import org.apache.sysds.test.TestUtils; +import org.junit.Assert; +import org.junit.Test; +import org.junit.runner.RunWith; +import org.junit.runners.Parameterized; + +@RunWith(value = Parameterized.class) + +@net.jcip.annotations.NotThreadSafe +public class BuiltinIsolationForestTest extends AutomatedTestBase +{ + private static final String TEST_NAME = "outlierByIsolationForestTest"; + private static final String TEST_DIR = "functions/builtin/"; + private static final String TEST_CLASS_DIR = TEST_DIR + BuiltinIsolationForestTest.class.getSimpleName() + "/"; + + private static final double SCORE_EPS = 1e-10; + private static final double R_SCORE_EPS = 1e-9; + private static final int DATA_SEED = 314159265; + private static final int FOREST_SEED = 42; + private static final int REFERENCE_SUBSAMPLING_SIZE = 4; + + private enum TestCase { + BASIC, + ANOMALY_RANKING, + SUBSAMPLING_CLAMP, + SINGLE_ROW_APPLY, + SINGLE_TREE, + CONSTANT_DATA, + MINIMUM_DATASET, + SINGLE_FEATURE, + HIGH_DIMENSIONAL, + PARTIALLY_CONSTANT + } + + private final TestCase testCase; + private final int numRows; + private final int numCols; + private final int numTrees; + private final int subsamplingSize; + + public BuiltinIsolationForestTest(TestCase testCase, int numRows, int numCols, + int numTrees, int subsamplingSize) + { + this.testCase = testCase; + this.numRows = numRows; + this.numCols = numCols; + this.numTrees = numTrees; + this.subsamplingSize = subsamplingSize; + } + + @Override + public void setUp() { + TestUtils.clearAssertionInformation(); + addTestConfiguration(TEST_NAME, new TestConfiguration(TEST_CLASS_DIR, TEST_NAME, + new String[] {"scores", "model", "subsampling_size", "same_seed_model", + "different_seed_model", "reference_scores", "dml_apply_runtime", "model_for_r"})); + } + + @Test + public void testIsolationForestSingleNode() { + runIsolationForestTest(ExecMode.SINGLE_NODE); + } + + @Test + public void testIsolationForestHybrid() { + runIsolationForestTest(ExecMode.HYBRID); + } + + private void runIsolationForestTest(ExecMode mode) { + ExecMode platformOld = setExecMode(mode); + + try { + loadTestConfiguration(getTestConfiguration(TEST_NAME)); + + double[][] X = createTrainingData(); + double[][] X_apply = createApplyData(X); + boolean runExtendedChecks = testCase == TestCase.BASIC && mode == ExecMode.SINGLE_NODE; + + writeInputMatrixWithMTD("X", X, true); + writeInputMatrixWithMTD("X_apply", X_apply, true); + writeInputMatrixWithMTD("ReferenceModel", referenceModel(), true); + writeInputMatrixWithMTD("ReferenceX", referenceSamples(), true); + + String home = SCRIPT_DIR + TEST_DIR; + fullDMLScriptName = home + TEST_NAME + ".dml"; + programArgs = new String[] {"-nvargs", + "X=" + input("X"), + "X_apply=" + input("X_apply"), + "reference_model=" + input("ReferenceModel"), + "reference_X=" + input("ReferenceX"), + "n_trees=" + numTrees, + "subsampling_size=" + subsamplingSize, + "seed=" + FOREST_SEED, + "check_seed=" + (runExtendedChecks ? 1 : 0), + "scores=" + output("scores"), + "model=" + output("model"), + "subsampling_size_out=" + output("subsampling_size"), + "same_seed_model=" + output("same_seed_model"), + "different_seed_model=" + output("different_seed_model"), + "reference_scores=" + output("reference_scores"), + "dml_apply_runtime=" + output("dml_apply_runtime"), + "model_for_r=" + output("model_for_r")}; + + runTest(true, EXCEPTION_NOT_EXPECTED, null, -1); + + int actualSubsamplingSize = Math.min(subsamplingSize, numRows); + HashMap scores = readDMLMatrixFromOutputDir("scores"); + HashMap model = readDMLMatrixFromOutputDir("model"); + + assertOutputDimensions(X_apply.length, actualSubsamplingSize); + assertScores(scores, X_apply.length); + assertModelStructure(model, X, actualSubsamplingSize); + assertScenarioSpecificResults(scores, X_apply.length); + assertReferenceScores(); + + if (runExtendedChecks) { + assertSeedBehavior(model); + assertRReference(scores, actualSubsamplingSize, X_apply.length, home); + } + } + finally { + resetExecMode(platformOld); + } + } + + private void assertOutputDimensions(int numApplyRows, int actualSubsamplingSize) { + MatrixCharacteristics scoreMeta = readDMLMetaDataFile("scores"); + Assert.assertEquals("One anomaly score is expected per input row", numApplyRows, scoreMeta.getRows()); + Assert.assertEquals("Anomaly scores must form a column vector", 1, scoreMeta.getCols()); + + int heightLimit = (int) Math.ceil(Math.log(actualSubsamplingSize) / Math.log(2)); + long expectedModelCols = 2L * ((1L << (heightLimit + 1)) - 1); + MatrixCharacteristics modelMeta = readDMLMetaDataFile("model"); + Assert.assertEquals("The model must contain exactly numTrees trees", numTrees, modelMeta.getRows()); + Assert.assertEquals("Every tree must be padded to the configured height limit", + expectedModelCols, modelMeta.getCols()); + + HashMap subsamplingSizeOut = readDMLScalarFromOutputDir("subsampling_size"); + Double actualValue = subsamplingSizeOut.get(new CellIndex(1, 1)); + Assert.assertNotNull("The effective subsampling size must be stored in the model", actualValue); + Assert.assertEquals("The stored value must reflect clamping to the number of training rows", + actualSubsamplingSize, actualValue.intValue()); + } + + private void assertScores(HashMap scores, int expectedRows) { + for (int row = 1; row <= expectedRows; row++) { + double score = valueAt(scores, row, 1); + Assert.assertTrue("Anomaly scores must be finite", Double.isFinite(score)); + Assert.assertTrue("Anomaly scores must be in the interval (0, 1]", score > 0 && score <= 1); + } + } + + private void assertModelStructure(HashMap model, double[][] X, + int actualSubsamplingSize) + { + MatrixCharacteristics modelMeta = readDMLMetaDataFile("model"); + int numNodes = (int) modelMeta.getCols() / 2; + int heightLimit = (int) Math.ceil(Math.log(actualSubsamplingSize) / Math.log(2)); + double[][] featureRanges = featureRanges(X); + + for (int tree = 1; tree <= numTrees; tree++) { + int leafSizeSum = 0; + + for (int nodeId = 1; nodeId <= numNodes; nodeId++) { + int nodeColumn = 2 * nodeId - 1; + double rawNodeType = valueAt(model, tree, nodeColumn); + int nodeType = (int) Math.rint(rawNodeType); + double nodeValue = valueAt(model, tree, nodeColumn + 1); + + Assert.assertEquals("Node types must be integer feature indices or the 0/-1 sentinels", + nodeType, rawNodeType, SCORE_EPS); + + if (nodeId > 1) { + int parentId = nodeId / 2; + int parentType = (int) Math.rint(valueAt(model, tree, 2 * parentId - 1)); + if (parentType > 0) + Assert.assertNotEquals("Both children of an internal node must be materialized", -1, nodeType); + else + Assert.assertEquals("Children below leaves and placeholders must remain placeholders", -1, nodeType); + } + + if (nodeType == -1) { + Assert.assertEquals("Both entries of a placeholder node must be -1", -1, nodeValue, SCORE_EPS); + } + else if (nodeType == 0) { + Assert.assertTrue("Every external node must contain at least one sample", nodeValue >= 1); + Assert.assertEquals("External-node sizes must be integral", Math.rint(nodeValue), nodeValue, SCORE_EPS); + leafSizeSum += (int) nodeValue; + } + else if (nodeType > 0) { + Assert.assertTrue("Split feature index exceeds the training width", nodeType <= numCols); + Assert.assertTrue("Split values must be finite", Double.isFinite(nodeValue)); + Assert.assertTrue("Split value is below the observed feature range", + nodeValue >= featureRanges[nodeType - 1][0]); + Assert.assertTrue("Split value is above the observed feature range", + nodeValue <= featureRanges[nodeType - 1][1]); + + int nodeDepth = 31 - Integer.numberOfLeadingZeros(nodeId); + Assert.assertTrue("Nodes at the height limit must be external", nodeDepth < heightLimit); + + if (testCase == TestCase.PARTIALLY_CONSTANT) + Assert.assertEquals("Constant columns must never be selected for a split", numCols, nodeType); + } + else { + Assert.fail("Invalid node type " + nodeType + " in tree " + tree); + } + } + + Assert.assertEquals("External nodes must partition the complete training subsample", + actualSubsamplingSize, leafSizeSum); + } + } + + private void assertScenarioSpecificResults(HashMap scores, int numApplyRows) { + if (testCase == TestCase.CONSTANT_DATA || testCase == TestCase.MINIMUM_DATASET) { + for (int row = 1; row <= numApplyRows; row++) + Assert.assertEquals("This scenario has the exact theoretical score 0.5", + 0.5, valueAt(scores, row, 1), SCORE_EPS); + } + else if (testCase == TestCase.ANOMALY_RANKING) { + double normalMean = (valueAt(scores, 1, 1) + valueAt(scores, 2, 1)) / 2; + double outlierMean = (valueAt(scores, 3, 1) + valueAt(scores, 4, 1)) / 2; + Assert.assertTrue("Points far outside the training range must rank above in-distribution points", + outlierMean > normalMean); + } + } + + private void assertReferenceScores() { + HashMap scores = readDMLMatrixFromOutputDir("reference_scores"); + double normalization = averagePathLength(REFERENCE_SUBSAMPLING_SIZE); + + Assert.assertEquals(isolationScore(1, normalization), valueAt(scores, 1, 1), SCORE_EPS); + Assert.assertEquals(isolationScore(2, normalization), valueAt(scores, 2, 1), SCORE_EPS); + Assert.assertEquals(isolationScore(3, normalization), valueAt(scores, 3, 1), SCORE_EPS); + } + + private void assertSeedBehavior(HashMap model) { + HashMap sameSeedModel = readDMLMatrixFromOutputDir("same_seed_model"); + HashMap differentSeedModel = readDMLMatrixFromOutputDir("different_seed_model"); + + Assert.assertEquals("The same input and seed must reproduce the exact forest", model, sameSeedModel); + Assert.assertNotEquals("Changing the seed must change at least one tree", model, differentSeedModel); + } + + // R scores the serialized SystemDS forest instead of training a second random + // forest whose RNG stream would not be element-wise comparable. + private void assertRReference(HashMap dmlScores, int actualSubsamplingSize, + int numApplyRows, String home) + { + fullRScriptName = home + TEST_NAME + ".R"; + rCmd = getRCmd( + output("model_for_r"), + input("X_apply.mtx"), + Integer.toString(actualSubsamplingSize), + expected("scores_R"), + expected("runtime_R")); + + runRScript(true); + + HashMap rScores = readRMatrixFromExpectedDir("scores_R"); + for (int row = 1; row <= numApplyRows; row++) + Assert.assertEquals("R and SystemDS must agree on the score for row " + row, + valueAt(dmlScores, row, 1), valueAt(rScores, row, 1), R_SCORE_EPS); + + double dmlRuntime = valueAt(readDMLScalarFromOutputDir("dml_apply_runtime"), 1, 1); + double rRuntime = valueAt(readRMatrixFromExpectedDir("runtime_R"), 1, 1); + Assert.assertTrue("SystemDS Apply runtime must be finite and non-negative", + Double.isFinite(dmlRuntime) && dmlRuntime >= 0); + Assert.assertTrue("R Apply runtime must be finite and non-negative", + Double.isFinite(rRuntime) && rRuntime >= 0); + + System.out.printf(Locale.ROOT, "Isolation Forest Apply runtime: SystemDS %.6f s, R %.6f s%n", + dmlRuntime, rRuntime); + } + + private double[][] createTrainingData() { + switch (testCase) { + case CONSTANT_DATA: + return constantMatrix(numRows, numCols, 7); + case MINIMUM_DATASET: + return new double[][] {{0.01, 0.02}, {0.03, 0.04}}; + case PARTIALLY_CONSTANT: + return partiallyConstantMatrix(numRows, numCols); + default: + return getRandomMatrix(numRows, numCols, -1, 1, 1, + DATA_SEED + testCase.ordinal()); + } + } + + private double[][] createApplyData(double[][] X) { + if (testCase == TestCase.ANOMALY_RANKING) { + double[][] result = new double[4][numCols]; + System.arraycopy(X[numRows / 3], 0, result[0], 0, numCols); + System.arraycopy(X[2 * numRows / 3], 0, result[1], 0, numCols); + Arrays.fill(result[2], -10); + Arrays.fill(result[3], 10); + return result; + } + else if (testCase == TestCase.SINGLE_ROW_APPLY) + return copyFirstRows(X, 1); + else if (testCase == TestCase.HIGH_DIMENSIONAL) + return copyFirstRows(X, 10); + + return X; + } + + private static double[][] constantMatrix(int rows, int cols, double value) { + double[][] result = new double[rows][cols]; + for (double[] row : result) + Arrays.fill(row, value); + return result; + } + + private static double[][] partiallyConstantMatrix(int rows, int cols) { + double[][] result = constantMatrix(rows, cols, 7); + for (int row = 0; row < rows; row++) + result[row][cols - 1] = (double) row / (rows - 1); + return result; + } + + private static double[][] copyFirstRows(double[][] X, int rows) { + double[][] result = new double[rows][X[0].length]; + for (int row = 0; row < rows; row++) + System.arraycopy(X[row], 0, result[row], 0, X[row].length); + return result; + } + + private static double[][] featureRanges(double[][] X) { + double[][] result = new double[X[0].length][2]; + for (double[] range : result) { + range[0] = Double.POSITIVE_INFINITY; + range[1] = Double.NEGATIVE_INFINITY; + } + + for (double[] row : X) { + for (int col = 0; col < row.length; col++) { + result[col][0] = Math.min(result[col][0], row[col]); + result[col][1] = Math.max(result[col][1], row[col]); + } + } + return result; + } + + private static double[][] referenceModel() { + return new double[][] {{ + 1, 0.5, + 0, 1, 1, 1.5, + -1, -1, -1, -1, 0, 1, 0, 2 + }}; + } + + private static double[][] referenceSamples() { + return new double[][] {{0.25}, {0.5}, {2}}; + } + + private static double averagePathLength(int n) { + double harmonic = 0; + for (int i = 1; i < n; i++) + harmonic += 1.0 / i; + return 2 * harmonic - 2.0 * (n - 1) / n; + } + + private static double isolationScore(int pathLength, double normalization) { + return Math.pow(2, -pathLength / normalization); + } + + private static double valueAt(HashMap matrix, int row, int col) { + return matrix.getOrDefault(new CellIndex(row, col), 0.0); + } + + @Parameterized.Parameters(name = "{0}: rows={1}, cols={2}, trees={3}, subsample={4}") + public static Collection data() { + // SCHEMA: TEST_CASE, #ROWS, #COLS, #TREES, REQUESTED_SUBSAMPLE_SIZE + return Arrays.asList(new Object[][] { + {TestCase.BASIC, 100, 3, 10, 20}, + {TestCase.ANOMALY_RANKING, 100, 3, 50, 64}, + {TestCase.SUBSAMPLING_CLAMP, 5, 3, 3, 256}, + {TestCase.SINGLE_ROW_APPLY, 20, 3, 5, 10}, + {TestCase.SINGLE_TREE, 20, 3, 1, 10}, + {TestCase.CONSTANT_DATA, 5, 3, 1, 5}, + {TestCase.MINIMUM_DATASET, 2, 2, 5, 2}, + {TestCase.SINGLE_FEATURE, 100, 1, 20, 50}, + {TestCase.HIGH_DIMENSIONAL, 200, 30, 20, 100}, + {TestCase.PARTIALLY_CONSTANT, 80, 5, 20, 40} + }); + } +} diff --git a/src/test/scripts/functions/builtin/outlierByIsolationForestTest.R b/src/test/scripts/functions/builtin/outlierByIsolationForestTest.R new file mode 100644 index 00000000000..168d2c71df6 --- /dev/null +++ b/src/test/scripts/functions/builtin/outlierByIsolationForestTest.R @@ -0,0 +1,129 @@ +#------------------------------------------------------------- +# +# 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. +# +#------------------------------------------------------------- + +# Independently applies a SystemDS Isolation Forest model in R. The model is +# trained by DML and serialized as MatrixMarket, so this script validates tree +# traversal and anomaly-score calculation without relying on matching random +# number generators between SystemDS and R. +# +# Arguments: +# 1. Linearized Isolation Forest model in MatrixMarket format +# 2. Samples to score in MatrixMarket format +# 3. Effective training subsampling size +# 4. Output anomaly scores in MatrixMarket format +# 5. Output Apply runtime in seconds as a 1x1 MatrixMarket matrix + +args <- commandArgs(trailingOnly = TRUE) +if (length(args) != 5) + stop("Expected model, samples, subsampling size, score output, and runtime output arguments.") + +suppressPackageStartupMessages(library("Matrix")) + +model <- as.matrix(readMM(args[1])) +X <- as.matrix(readMM(args[2])) +subsampling_size <- as.integer(args[3]) + +average_path_length <- function(n) { + if (n <= 1) + stop("average_path_length requires n > 1") + + if (n < 1000) + harmonic <- sum(1 / seq_len(n - 1)) + else + harmonic <- log(n - 1) + 0.57721566490153 + + 2 * harmonic - 2 * (n - 1) / n +} + +tree_path_length <- function(tree, x) { + node_id <- 1L + edges <- 0L + + repeat { + node_start <- 2L * node_id - 1L + if (node_start + 1L > length(tree)) + stop("Invalid iTree model: node index is out of bounds.") + + split_feature <- as.integer(round(tree[node_start])) + node_value <- tree[node_start + 1L] + + if (split_feature > 0L) { + if (split_feature > length(x)) + stop("Invalid iTree model: split feature exceeds the input width.") + + edges <- edges + 1L + if (x[split_feature] < node_value) + node_id <- 2L * node_id + else + node_id <- 2L * node_id + 1L + } + else if (split_feature == 0L) { + leaf_size <- as.integer(round(node_value)) + if (leaf_size < 1L) + stop("Invalid iTree model: external-node size must be positive.") + + if (leaf_size <= 1L) + return(as.numeric(edges)) + return(edges + average_path_length(leaf_size)) + } + else { + stop("Invalid iTree model: reached a placeholder node.") + } + } +} + +score_forest <- function(model, X, subsampling_size) { + if (nrow(model) < 1L) + stop("The model must contain at least one tree.") + if (nrow(X) < 1L) + stop("X must contain at least one row.") + if (subsampling_size <= 1L) + stop("subsampling_size must be greater than one.") + + height_limit <- ceiling(log(subsampling_size, base = 2)) + expected_columns <- 2 * (2^(height_limit + 1) - 1) + if (ncol(model) != expected_columns) + stop("The model has an invalid number of columns.") + + normalization <- average_path_length(subsampling_size) + num_samples <- nrow(X) + scores <- matrix(0, nrow = num_samples, ncol = 1L) + + for (sample_id in seq_len(num_samples)) { + path_sum <- 0 + for (tree_id in seq_len(nrow(model))) + path_sum <- path_sum + tree_path_length(model[tree_id, ], X[sample_id, ]) + + scores[sample_id, 1L] <- 2^(-(path_sum / nrow(model)) / normalization) + } + + scores +} + +timing <- system.time({ + scores <- score_forest(model, X, subsampling_size) +}) +apply_runtime <- unname(timing[["elapsed"]]) + +invisible(writeMM(Matrix(scores, sparse = TRUE), args[4])) +invisible(writeMM(Matrix(matrix(apply_runtime, nrow = 1L), sparse = TRUE), args[5])) + +cat(sprintf("R Isolation Forest Apply runtime: %.6f s\n", apply_runtime)) diff --git a/src/test/scripts/functions/builtin/outlierByIsolationForestTest.dml b/src/test/scripts/functions/builtin/outlierByIsolationForestTest.dml new file mode 100644 index 00000000000..b073074da93 --- /dev/null +++ b/src/test/scripts/functions/builtin/outlierByIsolationForestTest.dml @@ -0,0 +1,87 @@ +#------------------------------------------------------------- +# +# 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. +# +#------------------------------------------------------------- + +X = read($X) +X_apply = read($X_apply) + +iForestModel = outlierByIsolationForest( + X = X, + n_trees = $n_trees, + subsampling_size = $subsampling_size, + seed = $seed +) + +M = as.matrix(iForestModel["model"]) +actual_subsampling_size = as.scalar(iForestModel["subsampling_size"]) + +apply_start = time() +anomaly_scores = outlierByIsolationForestApply( + iForestModel = iForestModel, + X = X_apply +) +dml_apply_runtime = (time() - apply_start) / 1e9 + +# A repeated training run verifies deterministic seed handling. A second seed +# must produce a different model, proving that the seed is not simply ignored. +M_same_seed = M +M_different_seed = M +if ($check_seed == 1) { + sameSeedModel = outlierByIsolationForest( + X = X, + n_trees = $n_trees, + subsampling_size = $subsampling_size, + seed = $seed + ) + differentSeedModel = outlierByIsolationForest( + X = X, + n_trees = $n_trees, + subsampling_size = $subsampling_size, + seed = $seed + 1 + ) + M_same_seed = as.matrix(sameSeedModel["model"]) + M_different_seed = as.matrix(differentSeedModel["model"]) +} + +# This hand-built tree provides an implementation-independent oracle for the +# apply function, including the exact-threshold branch at x = 0.5. +reference_M = read($reference_model) +reference_X = read($reference_X) +referenceModel = list( + model = reference_M, + subsampling_size = 4 +) +reference_anomaly_scores = outlierByIsolationForestApply( + iForestModel = referenceModel, + X = reference_X +) + +write(anomaly_scores, $scores) +write(M, $model) +write(actual_subsampling_size, $subsampling_size_out) +write(M_same_seed, $same_seed_model) +write(M_different_seed, $different_seed_model) +write(reference_anomaly_scores, $reference_scores) +write(dml_apply_runtime, $dml_apply_runtime) + +# The R reference implementation scores this exact serialized forest. Writing +# only for the extended BASIC/SINGLE_NODE check avoids unnecessary test I/O. +if ($check_seed == 1) + write(M, $model_for_r, format = "mm")