diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index a329c4362..741b9b3b2 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -66,6 +66,11 @@ jobs: run: | sudo apt-get update sudo apt-get install -y libopenblas-dev libatlas-base-dev liblapack-dev gfortran libgmp-dev libmpfr-dev libsuitesparse-dev ccache libmpc-dev + sudo apt-get install -y gcc g++ + + - name: show-gcc + run: | + gcc --version - name: Install Python packages run: | diff --git a/.spin/cmds.py b/.spin/cmds.py index 1a5fc67e4..eb1fd8406 100644 --- a/.spin/cmds.py +++ b/.spin/cmds.py @@ -7,6 +7,8 @@ from spin import util from spin.cmds import meson +PROJECT_MODULE = "sktree" + def get_git_revision_hash(submodule) -> str: return subprocess.check_output(["git", "rev-parse", f"@:{submodule}"]).decode("ascii").strip() diff --git a/DEVELOPING.md b/DEVELOPING.md index a9d64457c..8956e1251 100644 --- a/DEVELOPING.md +++ b/DEVELOPING.md @@ -76,4 +76,16 @@ Verify that installations work as expected on your machine. twine upload dist/* -4. Update version number on ``meson.build`` and ``_version.py`` to the relevant version. \ No newline at end of file +4. Update version number on ``meson.build`` and ``_version.py`` to the relevant version. + +# Comparing branches using ASV benchmarks + +We use asv to compare performance between the current branch and main branch. For +example, one can run: + + asv continuous --verbose --split --bench ObliqueRandomForest origin/main constantsv2 + +# Extending Cython API + +Due to the current state of scikit-learn's internal Cython code for trees, we have to instead leverage a fork of scikit-learn at https://github.com/neurodata/scikit-learn when +extending the decision tree model API of scikit-learn. Specifically, we extend the Python and Cython API of the tree submodule in scikit-learn in our submodule, so we can introduce the tree models housed in this package. Thus these extend the functionality of decision-tree based models in a way that is not possible yet in scikit-learn itself. As one example, we introduce an abstract API to allow users to implement their own oblique splits. Our plan in the future is to benchmark these functionalities and introduce them upstream to scikit-learn where applicable and inclusion criterion are met. diff --git a/README.md b/README.md index 470641ef1..f1483986f 100644 --- a/README.md +++ b/README.md @@ -5,6 +5,7 @@ [![codecov](https://codecov.io/gh/neurodata/scikit-tree/branch/main/graph/badge.svg?token=H1reh7Qwf4)](https://codecov.io/gh/neurodata/scikit-tree) [![PyPI Download count](https://img.shields.io/pypi/dm/scikit-tree.svg)](https://pypistats.org/packages/scikit-tree) [![Latest PyPI release](https://img.shields.io/pypi/v/scikit-tree.svg)](https://pypi.org/project/scikit-tree/) +[![Benchmark](https://img.shields.io/badge/Benchmarked%20by-asv-blue)](https://img.shields.io/badge/Benchmarked%20by-asv-blue) scikit-tree =========== @@ -15,9 +16,6 @@ Tree-models have withstood the test of time, and are consistently used for moder We welcome contributions for modern tree-based algorithms. We use Cython to achieve fast C/C++ speeds, while abiding by a scikit-learn compatible (tested) API. Moreover, our Cython internals are easily extensible because they follow the internal Cython API of scikit-learn as well. -Due to the current state of scikit-learn's internal Cython code for trees, we have to instead leverage a fork of scikit-learn at https://github.com/neurodata/scikit-learn when -extending the decision tree model API of scikit-learn. Specifically, we extend the Python and Cython API of the tree submodule in scikit-learn in our submodule, so we can introduce the tree models housed in this package. Thus these extend the functionality of decision-tree based models in a way that is not possible yet in scikit-learn itself. As one example, we introduce an abstract API to allow users to implement their own oblique splits. Our plan in the future is to benchmark these functionalities and introduce them upstream to scikit-learn where applicable and inclusion criterion are met. - Documentation ============= @@ -25,7 +23,7 @@ See here for the documentation for our dev version: https://docs.neurodata.io/sc Why oblique trees and why trees beyond those in scikit-learn? ============================================================= -In 2001, Leo Breiman proposed two types of Random Forests. One was known as ``Forest-RI``, which is the axis-aligned traditional random forest. One was known as ``Forest-RC``, which is the random oblique linear combinations random forest. This leveraged random combinations of features to perform splits. [MORF](1) builds upon ``Forest-RC`` by proposing additional functions to combine features. Other modern tree variants such as Canonical Correlation Forests (CCF), or unsupervised random forests are also important at solving real-world problems using robust decision tree models. +In 2001, Leo Breiman proposed two types of Random Forests. One was known as ``Forest-RI``, which is the axis-aligned traditional random forest. One was known as ``Forest-RC``, which is the random oblique linear combinations random forest. This leveraged random combinations of features to perform splits. A generalization of ``Forest-RC`` known as [MORF](1) builds upon ``Forest-RC`` by proposing additional functions to combine features. Other modern tree variants such as extended isolation forests, or unsupervised random forests are also important at solving real-world problems using robust decision tree models. Installation ============ diff --git a/build_requirements.txt b/build_requirements.txt index 360c66711..167b9f41e 100644 --- a/build_requirements.txt +++ b/build_requirements.txt @@ -1,6 +1,6 @@ meson meson-python -cython>=3.0 +cython<3.0 ninja numpy scikit-learn>=1.3 diff --git a/doc/whats_new/v0.2.rst b/doc/whats_new/v0.2.rst index ae0e75305..80a333bd2 100644 --- a/doc/whats_new/v0.2.rst +++ b/doc/whats_new/v0.2.rst @@ -31,6 +31,7 @@ Changelog - |Efficiency| Around 1.5-2x speed improvement for unsupervised forests, by `Adam Li`_ (:pr:`114`) - |API| Allow ``sqrt`` and ``log2`` keywords to be used for ``min_samples_split`` parameter in unsupervised forests, by `Adam Li`_ (:pr:`114`) - |Feature| Implement extended isolation forest, by `Adam Li`_ (:pr:`101`) +- |Feature| Track constant columns within the tree to prevent splitting or evaluating columns in oblique splits that are constant, by `Adam Li`_ (:pr:`121`) Code and Documentation Contributors ----------------------------------- diff --git a/pyproject.toml b/pyproject.toml index 7bb840d16..73990e254 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -77,7 +77,7 @@ requires = [ "wheel", "setuptools<=65.5", "packaging", - "Cython>=0.29.24", + "Cython==0.29.36", "scikit-learn>=1.3", "lazy_loader>=0.1", diff --git a/sktree/_lib/sklearn_fork b/sktree/_lib/sklearn_fork index 68015082c..63e72410d 160000 --- a/sktree/_lib/sklearn_fork +++ b/sktree/_lib/sklearn_fork @@ -1 +1 @@ -Subproject commit 68015082cc740d7859fad964a77f9e684544d868 +Subproject commit 63e72410d7803d90aacb49645a1acd28d05ab0c6 diff --git a/sktree/tests/test_supervised_forest.py b/sktree/tests/test_supervised_forest.py index 045653490..75cb7d926 100644 --- a/sktree/tests/test_supervised_forest.py +++ b/sktree/tests/test_supervised_forest.py @@ -193,6 +193,7 @@ def _trunk(n, p=10, random_state=None): ) def test_sklearn_compatible_estimator(estimator, check): # TODO: remove when we can replicate the CI error... + # this seems to be due to a compiler issue since it is not replicable on MacOSx if isinstance( estimator, ( @@ -207,9 +208,9 @@ def test_sklearn_compatible_estimator(estimator, check): def test_oblique_forest_sparse_parity(): # Sparse parity dataset - n = 1000 + n = 500 X, y = _sparse_parity(n, random_state=0) - n_test = 0.1 + n_test = 0.2 X_train, X_test, y_train, y_test = train_test_split( X, y, @@ -217,7 +218,7 @@ def test_oblique_forest_sparse_parity(): random_state=0, ) - rc_clf = ObliqueRandomForestClassifier(max_features=None, random_state=0) + rc_clf = ObliqueRandomForestClassifier(random_state=0) rc_clf.fit(X_train, y_train) y_hat = rc_clf.predict(X_test) rc_accuracy = accuracy_score(y_test, y_hat) @@ -227,9 +228,11 @@ def test_oblique_forest_sparse_parity(): y_hat = ri_clf.predict(X_test) ri_accuracy = accuracy_score(y_test, y_hat) - assert ri_accuracy < rc_accuracy - assert ri_accuracy > 0.45 - assert rc_accuracy > 0.5 + assert ( + ri_accuracy < rc_accuracy + ), f"Oblique forest: {rc_accuracy} < Axis-aligned forest: {ri_accuracy}" + assert ri_accuracy > 0.45, f"Axis-aligned forest: {ri_accuracy} < 0.45" + assert rc_accuracy > 0.5, f"Oblique forest: {rc_accuracy} < 0.5" def test_oblique_forest_orthant(): @@ -258,7 +261,9 @@ def test_oblique_forest_orthant(): y_hat = ri_clf.predict(X_test) ri_accuracy = accuracy_score(y_test, y_hat) - assert rc_accuracy >= ri_accuracy + assert ( + rc_accuracy >= ri_accuracy + ), f"Oblique forest: {rc_accuracy} < Axis-aligned forest: {ri_accuracy}" assert ri_accuracy > 0.84 assert rc_accuracy > 0.85 diff --git a/sktree/tree/_oblique_splitter.pxd b/sktree/tree/_oblique_splitter.pxd index c9c833d90..98ceb1b7e 100644 --- a/sktree/tree/_oblique_splitter.pxd +++ b/sktree/tree/_oblique_splitter.pxd @@ -11,6 +11,7 @@ import numpy as np cimport numpy as cnp +from libcpp.unordered_map cimport unordered_map from libcpp.vector cimport vector from .._lib.sklearn.tree._criterion cimport Criterion @@ -46,6 +47,13 @@ cdef class BaseObliqueSplitter(Splitter): cdef vector[vector[DTYPE_t]] proj_mat_weights # nonzero weights of sparse proj_mat matrix cdef vector[vector[SIZE_t]] proj_mat_indices # nonzero indices of sparse proj_mat matrix + # keep a hashmap of every projection vector indices sampled + cdef unordered_map[size_t, bint] proj_vec_hash + + cdef unordered_map[SIZE_t, DTYPE_t] min_val_map + cdef unordered_map[SIZE_t, DTYPE_t] max_val_map + cdef unordered_map[SIZE_t, bint] constant_col_map + # TODO: assumes all oblique splitters only work with dense data cdef const DTYPE_t[:, :] X @@ -58,7 +66,8 @@ cdef class BaseObliqueSplitter(Splitter): cdef void sample_proj_mat( self, vector[vector[DTYPE_t]]& proj_mat_weights, - vector[vector[SIZE_t]]& proj_mat_indices + vector[vector[SIZE_t]]& proj_mat_indices, + SIZE_t n_known_constants, ) noexcept nogil # Redefined here since the new logic requires calling sample_proj_mat @@ -76,7 +85,8 @@ cdef class BaseObliqueSplitter(Splitter): const SIZE_t[:] samples, DTYPE_t[:] feature_values, vector[DTYPE_t]* proj_vec_weights, # weights of the vector (max_features,) - vector[SIZE_t]* proj_vec_indices # indices of the features (max_features,) + vector[SIZE_t]* proj_vec_indices, # indices of the features (max_features,) + SIZE_t* n_known_constants ) noexcept nogil cdef int node_split( @@ -97,6 +107,7 @@ cdef class ObliqueSplitter(BaseObliqueSplitter): cdef public double feature_combinations # Number of features to combine cdef SIZE_t n_non_zeros # Number of non-zero features cdef SIZE_t[::1] indices_to_sample # an array of indices to sample of size mtry X n_features + cdef SIZE_t floor_feature_combinations # All oblique splitters (i.e. non-axis aligned splitters) require a # function to sample a projection matrix that is applied to the feature matrix @@ -104,5 +115,6 @@ cdef class ObliqueSplitter(BaseObliqueSplitter): cdef void sample_proj_mat( self, vector[vector[DTYPE_t]]& proj_mat_weights, - vector[vector[SIZE_t]]& proj_mat_indices + vector[vector[SIZE_t]]& proj_mat_indices, + SIZE_t n_known_constants, ) noexcept nogil diff --git a/sktree/tree/_oblique_splitter.pyx b/sktree/tree/_oblique_splitter.pyx index bc29cafa7..66dc12bcb 100644 --- a/sktree/tree/_oblique_splitter.pyx +++ b/sktree/tree/_oblique_splitter.pyx @@ -81,7 +81,8 @@ cdef class BaseObliqueSplitter(Splitter): cdef void sample_proj_mat( self, vector[vector[DTYPE_t]]& proj_mat_weights, - vector[vector[SIZE_t]]& proj_mat_indices + vector[vector[SIZE_t]]& proj_mat_indices, + SIZE_t n_known_constants, ) noexcept nogil: """ Sample the projection vector. @@ -94,14 +95,15 @@ cdef class BaseObliqueSplitter(Splitter): return sizeof(ObliqueSplitRecord) - cdef inline void compute_features_over_samples( + cdef void compute_features_over_samples( self, SIZE_t start, SIZE_t end, const SIZE_t[:] samples, DTYPE_t[:] feature_values, vector[DTYPE_t]* proj_vec_weights, # weights of the vector (max_features,) - vector[SIZE_t]* proj_vec_indices # indices of the features (max_features,) + vector[SIZE_t]* proj_vec_indices, # indices of the features (max_features,) + SIZE_t* n_known_constants ) noexcept nogil: """Compute the feature values for the samples[start:end] range. @@ -164,7 +166,7 @@ cdef class BaseObliqueSplitter(Splitter): _init_split(&best_split, end) # Sample the projection matrix - self.sample_proj_mat(self.proj_mat_weights, self.proj_mat_indices) + self.sample_proj_mat(self.proj_mat_weights, self.proj_mat_indices, 0) # For every vector in the projection matrix for feat_i in range(max_features): @@ -186,7 +188,8 @@ cdef class BaseObliqueSplitter(Splitter): samples, feature_values, &self.proj_mat_weights[feat_i], - &self.proj_mat_indices[feat_i] + &self.proj_mat_indices[feat_i], + n_constant_features ) # Sort the samples @@ -326,6 +329,7 @@ cdef class ObliqueSplitter(BaseObliqueSplitter): # Oblique tree parameters self.feature_combinations = feature_combinations + self.floor_feature_combinations = (self.feature_combinations) # or max w/ 1... self.n_non_zeros = max(int(self.max_features * self.feature_combinations), 1) @@ -351,6 +355,16 @@ cdef class ObliqueSplitter(BaseObliqueSplitter): self.indices_to_sample = np.arange(self.max_features * self.n_features, dtype=np.intp) + # re-initialize the hashmap for looking at constant features + cdef unordered_map[SIZE_t, DTYPE_t] min_val_map + cdef unordered_map[SIZE_t, DTYPE_t] max_val_map + self.min_val_map = min_val_map + self.max_val_map = max_val_map + + # re-initialize the hashmap for projection vector hash + cdef unordered_map[size_t, bint] proj_vec_hash + self.proj_vec_hash = proj_vec_hash + # XXX: Just to initialize stuff # self.feature_weights = np.ones((self.n_features,), dtype=DTYPE_t) / self.n_features return 0 @@ -358,7 +372,8 @@ cdef class ObliqueSplitter(BaseObliqueSplitter): cdef void sample_proj_mat( self, vector[vector[DTYPE_t]]& proj_mat_weights, - vector[vector[SIZE_t]]& proj_mat_indices + vector[vector[SIZE_t]]& proj_mat_indices, + SIZE_t n_known_constants, ) noexcept nogil: """Sample oblique projection matrix. @@ -385,6 +400,7 @@ cdef class ObliqueSplitter(BaseObliqueSplitter): cdef SIZE_t n_non_zeros = self.n_non_zeros cdef UINT32_t* random_state = &self.rand_r_state + cdef SIZE_t col_idx cdef int i, feat_i, proj_i, rand_vec_index cdef DTYPE_t weight @@ -406,12 +422,13 @@ cdef class ObliqueSplitter(BaseObliqueSplitter): # get the projection index and feature index proj_i = rand_vec_index // n_features - feat_i = rand_vec_index % n_features + feat_i = rand_vec_index % (n_features - n_known_constants) + col_idx = self.features[feat_i + n_known_constants] # sample a random weight weight = 1 if (rand_int(0, 2, random_state) == 1) else -1 - proj_mat_indices[proj_i].push_back(feat_i) # Store index of nonzero + proj_mat_indices[proj_i].push_back(col_idx) # Store index of nonzero proj_mat_weights[proj_i].push_back(weight) # Store weight of nonzero cdef class BestObliqueSplitter(ObliqueSplitter): @@ -428,6 +445,168 @@ cdef class BestObliqueSplitter(ObliqueSplitter): self.feature_combinations, ), self.__getstate__()) + cdef int node_split( + self, + double impurity, + SplitRecord* split, + SIZE_t* n_constant_features, + double lower_bound, + double upper_bound, + ) except -1 nogil: + """Find the best_split split on node samples[start:end] + Returns -1 in case of failure to allocate memory (and raise MemoryError) + or 0 otherwise. + """ + # typecast the pointer to an ObliqueSplitRecord + cdef ObliqueSplitRecord* oblique_split = (split) + + # use index sample array to allow quick re-ordering of samples + cdef SIZE_t[::1] samples = self.samples + cdef SIZE_t start = self.start + cdef SIZE_t end = self.end + cdef SIZE_t n_known_constants = n_constant_features[0] + + # The number of sampled feature values (i.e. mtry) + cdef SIZE_t n_visited_features = 0 + + # pointer array to store feature values to split on + cdef DTYPE_t[::1] feature_values = self.feature_values + cdef SIZE_t max_features = self.max_features + cdef SIZE_t min_samples_leaf = self.min_samples_leaf + + # keep track of split record for current_split node and the best_split split + # found among the sampled projection vectors + cdef ObliqueSplitRecord best_split, current_split + cdef double current_proxy_improvement = -INFINITY + cdef double best_proxy_improvement = -INFINITY + + cdef SIZE_t p # index over computed features and start/end + cdef SIZE_t partition_end + cdef DTYPE_t temp_d # to compute a projection feature value + + # instantiate the split records + _init_split(&best_split, end) + + # Sample the projection matrix + self.sample_proj_mat(self.proj_mat_weights, self.proj_mat_indices, n_known_constants) + + # Note: Compared to axis-aligned sampling, we do not keep track of "constant" features + # values that are drawn and enforce that we continue drawing features even if + # we go over the `max_features` mtry. This implicitly assumes that if we sample a new + # projection vector, the probability of sampling a constant feature projection is relatively low. + while (self.n_features > n_known_constants and # Stop early if remaining features + # are constant, or + # if we have reached max_features mtry + n_visited_features < max_features): + + if self.proj_mat_weights[n_visited_features].empty() or self.proj_mat_indices[n_visited_features].empty(): + # increment the mtry + n_visited_features += 1 + continue + + # XXX: 'feature' is not actually used in oblique split records + # Just indicates which split was sampled + current_split.feature = n_visited_features + current_split.proj_vec_weights = &self.proj_mat_weights[n_visited_features] + current_split.proj_vec_indices = &self.proj_mat_indices[n_visited_features] + + # Compute linear combination of features + self.compute_features_over_samples( + start, + end, + samples, + feature_values, + &self.proj_mat_weights[n_visited_features], + &self.proj_mat_indices[n_visited_features], + &n_known_constants + ) + + # Sort the samples + sort(&feature_values[start], &samples[start], end - start) + + # increment the mtry + n_visited_features += 1 + + # Evaluate all splits + self.criterion.reset() + p = start + while p < end: + while (p + 1 < end and feature_values[p + 1] <= feature_values[p] + FEATURE_THRESHOLD): + p += 1 + + p += 1 + + if p < end: + current_split.pos = p + + # Reject if min_samples_leaf is not guaranteed + if (((current_split.pos - start) < min_samples_leaf) or + ((end - current_split.pos) < min_samples_leaf)): + continue + + self.criterion.update(current_split.pos) + # Reject if min_weight_leaf is not satisfied + if self.check_postsplit_conditions() == 1: + continue + + current_proxy_improvement = \ + self.criterion.proxy_impurity_improvement() + + if current_proxy_improvement > best_proxy_improvement: + best_proxy_improvement = current_proxy_improvement + # sum of halves is used to avoid infinite value + current_split.threshold = feature_values[p - 1] / 2.0 + feature_values[p] / 2.0 + + if ( + (current_split.threshold == feature_values[p]) or + (current_split.threshold == INFINITY) or + (current_split.threshold == -INFINITY) + ): + current_split.threshold = feature_values[p - 1] + + best_split = current_split # copy + + # Reorganize into samples[start:best_split.pos] + samples[best_split.pos:end] + if best_split.pos < end: + partition_end = end + p = start + + while p < partition_end: + # Account for projection vector + temp_d = 0.0 + for j in range(best_split.proj_vec_indices.size()): + temp_d += self.X[samples[p], deref(best_split.proj_vec_indices)[j]] *\ + deref(best_split.proj_vec_weights)[j] + + if temp_d <= best_split.threshold: + p += 1 + + else: + partition_end -= 1 + samples[p], samples[partition_end] = \ + samples[partition_end], samples[p] + + self.criterion.reset() + self.criterion.update(best_split.pos) + self.criterion.children_impurity(&best_split.impurity_left, + &best_split.impurity_right) + best_split.improvement = self.criterion.impurity_improvement( + impurity, best_split.impurity_left, best_split.impurity_right) + + # keep track of the known constants at each depth of the tree + n_constant_features[0] = n_known_constants + + # Return values + deref(oblique_split).proj_vec_indices = best_split.proj_vec_indices + deref(oblique_split).proj_vec_weights = best_split.proj_vec_weights + deref(oblique_split).feature = best_split.feature + deref(oblique_split).pos = best_split.pos + deref(oblique_split).threshold = best_split.threshold + deref(oblique_split).improvement = best_split.improvement + deref(oblique_split).impurity_left = best_split.impurity_left + deref(oblique_split).impurity_right = best_split.impurity_right + return 0 + cdef class RandomObliqueSplitter(ObliqueSplitter): def __reduce__(self): """Enable pickling the splitter.""" @@ -505,10 +684,14 @@ cdef class RandomObliqueSplitter(ObliqueSplitter): # typecast the pointer to an ObliqueSplitRecord cdef ObliqueSplitRecord* oblique_split = (split) - # Draw random splits and pick the best_split + # use index sample array to allow quick re-ordering of samples cdef SIZE_t[::1] samples = self.samples cdef SIZE_t start = self.start cdef SIZE_t end = self.end + cdef SIZE_t n_known_constants = n_constant_features[0] + + # The number of sampled feature values (i.e. mtry) + cdef SIZE_t n_visited_features = 0 cdef UINT32_t* random_state = &self.rand_r_state # pointer array to store feature values to split on @@ -524,42 +707,38 @@ cdef class RandomObliqueSplitter(ObliqueSplitter): cdef double best_proxy_improvement = -INFINITY cdef SIZE_t p - cdef SIZE_t feat_i cdef SIZE_t partition_end cdef DTYPE_t temp_d # to compute a projection feature value cdef DTYPE_t min_feature_value cdef DTYPE_t max_feature_value - # Number of features discovered to be constant during the split search - # cdef SIZE_t n_found_constants = 0 - # cdef SIZE_t n_known_constants = n_constant_features[0] - # n_total_constants = n_known_constants + n_found_constants - # cdef SIZE_t n_total_constants = n_known_constants - cdef SIZE_t n_visited_features = 0 - # instantiate the split records _init_split(&best_split, end) # Sample the projection matrix - self.sample_proj_mat(self.proj_mat_weights, self.proj_mat_indices) - - # For every vector in the projection matrix - for feat_i in range(max_features): - # Break if already reached max_features - if n_visited_features >= max_features: - break - # Skip features known to be constant - # if feat_i < n_total_constants: - # continue - # Projection vector has no nonzeros - if self.proj_mat_weights[feat_i].empty(): + self.sample_proj_mat(self.proj_mat_weights, self.proj_mat_indices, n_known_constants) + + # Note: Compared to axis-aligned sampling, we do not keep track of "constant" features + # values that are drawn and enforce that we continue drawing features even if + # we go over the `max_features` mtry. This implicitly assumes that if we sample a new + # projection vector, the probability of sampling a constant feature projection is relatively low. + while (self.n_features > n_known_constants and # Stop early if remaining features + # are constant, or + # if we have reached max_features mtry + n_visited_features < max_features): + + if self.proj_mat_weights[n_visited_features].empty() or self.proj_mat_indices[n_visited_features].empty(): + # increment the mtry + n_visited_features += 1 continue + # sample a projection vector for the current mtry + # self.sample_proj_vector(self.proj_mat_weights, self.proj_mat_indices, n_known_constants, n_visited_features) # XXX: 'feature' is not actually used in oblique split records # Just indicates which split was sampled - current_split.feature = feat_i - current_split.proj_vec_weights = &self.proj_mat_weights[feat_i] - current_split.proj_vec_indices = &self.proj_mat_indices[feat_i] + current_split.feature = n_visited_features + current_split.proj_vec_weights = &self.proj_mat_weights[n_visited_features] + current_split.proj_vec_indices = &self.proj_mat_indices[n_visited_features] # Compute linear combination of features self.compute_features_over_samples( @@ -567,19 +746,18 @@ cdef class RandomObliqueSplitter(ObliqueSplitter): end, samples, feature_values, - &self.proj_mat_weights[feat_i], - &self.proj_mat_indices[feat_i] + &self.proj_mat_weights[n_visited_features], + &self.proj_mat_indices[n_visited_features], + &n_known_constants ) + # increment the count after all data structures that require this index + # are accessed and used + n_visited_features += 1 + # find min, max of the feature_values self.find_min_max(feature_values, &min_feature_value, &max_feature_value) - # XXX: Add logic to keep track of constant features if they exist - # if max_feature_value <= min_feature_value + FEATURE_THRESHOLD: - # n_found_constants += 1 - # n_total_constants += 1 - # continue - # Draw a random threshold current_split.threshold = rand_uniform( min_feature_value, @@ -613,8 +791,6 @@ cdef class RandomObliqueSplitter(ObliqueSplitter): best_proxy_improvement = current_proxy_improvement best_split = current_split # copy - n_visited_features += 1 - # Reorganize into samples[start:best_split.pos] + samples[best_split.pos:end] if best_split.pos < end: partition_end = end @@ -642,6 +818,9 @@ cdef class RandomObliqueSplitter(ObliqueSplitter): best_split.improvement = self.criterion.impurity_improvement( impurity, best_split.impurity_left, best_split.impurity_right) + # keep track of the known constants at each depth of the tree + n_constant_features[0] = n_known_constants + # Return values deref(oblique_split).proj_vec_indices = best_split.proj_vec_indices deref(oblique_split).proj_vec_weights = best_split.proj_vec_weights diff --git a/sktree/tree/_oblique_tree.pxd b/sktree/tree/_oblique_tree.pxd index 9f4906a4d..9b7d7aa54 100644 --- a/sktree/tree/_oblique_tree.pxd +++ b/sktree/tree/_oblique_tree.pxd @@ -37,7 +37,7 @@ cdef class ObliqueTree(Tree): SplitRecord* split_node, Node *node, SIZE_t node_id - ) nogil except -1 + ) except -1 nogil cdef DTYPE_t _compute_feature( self, const DTYPE_t[:, :] X_ndarray, diff --git a/sktree/tree/_utils.pxd b/sktree/tree/_utils.pxd index fc93163b2..a1efb4a13 100644 --- a/sktree/tree/_utils.pxd +++ b/sktree/tree/_utils.pxd @@ -1,6 +1,7 @@ import numpy as np cimport numpy as cnp +from libcpp.vector cimport vector cnp.import_array() @@ -23,3 +24,7 @@ cpdef ravel_multi_index(SIZE_t[:] coords, const SIZE_t[:] shape) cdef void unravel_index_cython(SIZE_t index, const SIZE_t[:] shape, SIZE_t[:] coords) noexcept nogil cdef SIZE_t ravel_multi_index_cython(SIZE_t[:] coords, const SIZE_t[:] shape) noexcept nogil + +cdef SIZE_t vector_hash( + const vector[SIZE_t]& v +) noexcept nogil diff --git a/sktree/tree/_utils.pyx b/sktree/tree/_utils.pyx index 1b7b6c293..e931f230e 100644 --- a/sktree/tree/_utils.pyx +++ b/sktree/tree/_utils.pyx @@ -145,3 +145,14 @@ cdef SIZE_t ravel_multi_index_cython(SIZE_t[:] coords, const SIZE_t[:] shape) no flat_index *= shape[i + 1] return flat_index + +cdef SIZE_t vector_hash( + const vector[SIZE_t]& v +) noexcept nogil: + """Hash a vector of size_t.""" + cdef SIZE_t seed = v.size() + cdef SIZE_t i, hash_val + for i in range(v.size()): + seed = seed ^ (v[i] + 0x9e3779b9 + (seed << 6) + (seed >> 2)) + hash_val = seed + return hash_val diff --git a/sktree/tree/manifold/_morf_splitter.pxd b/sktree/tree/manifold/_morf_splitter.pxd index 9387d6797..2c1642190 100644 --- a/sktree/tree/manifold/_morf_splitter.pxd +++ b/sktree/tree/manifold/_morf_splitter.pxd @@ -80,7 +80,8 @@ cdef class PatchSplitter(BaseObliqueSplitter): cdef void sample_proj_mat( self, vector[vector[DTYPE_t]]& proj_mat_weights, - vector[vector[SIZE_t]]& proj_mat_indices + vector[vector[SIZE_t]]& proj_mat_indices, + SIZE_t n_known_constants, ) noexcept nogil @@ -100,5 +101,6 @@ cdef class GaussianKernelSplitter(PatchSplitter): cdef void sample_proj_mat( self, vector[vector[DTYPE_t]]& proj_mat_weights, - vector[vector[SIZE_t]]& proj_mat_indices + vector[vector[SIZE_t]]& proj_mat_indices, + SIZE_t n_known_constants, ) noexcept nogil diff --git a/sktree/tree/manifold/_morf_splitter.pyx b/sktree/tree/manifold/_morf_splitter.pyx index 039d58008..b7e1197ea 100644 --- a/sktree/tree/manifold/_morf_splitter.pyx +++ b/sktree/tree/manifold/_morf_splitter.pyx @@ -82,7 +82,8 @@ cdef class PatchSplitter(BaseObliqueSplitter): cdef void sample_proj_mat( self, vector[vector[DTYPE_t]]& proj_mat_weights, - vector[vector[SIZE_t]]& proj_mat_indices + vector[vector[SIZE_t]]& proj_mat_indices, + SIZE_t n_known_constants, ) noexcept nogil: """ Sample the projection vector. @@ -271,7 +272,8 @@ cdef class BestPatchSplitter(BaseDensePatchSplitter): cdef void sample_proj_mat( self, vector[vector[DTYPE_t]]& proj_mat_weights, - vector[vector[SIZE_t]]& proj_mat_indices + vector[vector[SIZE_t]]& proj_mat_indices, + SIZE_t n_known_constants, ) noexcept nogil: """Sample projection matrix using a contiguous patch. @@ -410,7 +412,8 @@ cdef class BestPatchSplitter(BaseDensePatchSplitter): const SIZE_t[:] samples, DTYPE_t[:] feature_values, vector[DTYPE_t]* proj_vec_weights, # weights of the vector (max_features,) - vector[SIZE_t]* proj_vec_indices # indices of the features (max_features,) + vector[SIZE_t]* proj_vec_indices, # indices of the features (max_features,) + SIZE_t* n_known_constants ) noexcept nogil: """Compute the feature values for the samples[start:end] range. @@ -450,7 +453,7 @@ cdef class BestPatchSplitterTester(BestPatchSplitter): patch_dims = np.array(self.patch_dims_buff, dtype=np.intp) return top_left_patch_seed, patch_size, patch_dims - cpdef sample_projection_vector( + cpdef sample_projection_vector_py( self, SIZE_t proj_i, SIZE_t patch_size, @@ -492,7 +495,7 @@ cdef class BestPatchSplitterTester(BestPatchSplitter): cdef SIZE_t i, j # sample projection matrix in C/C++ - self.sample_proj_mat(proj_mat_weights, proj_mat_indices) + self.sample_proj_mat(proj_mat_weights, proj_mat_indices, 0) # convert the projection matrix to something that can be used in Python proj_vecs = np.zeros((self.max_features, self.n_features), dtype=np.float64) diff --git a/sktree/tree/meson.build b/sktree/tree/meson.build index 09ac76013..56090f51b 100644 --- a/sktree/tree/meson.build +++ b/sktree/tree/meson.build @@ -29,7 +29,6 @@ py3.install_sources( subdir: 'sktree/tree' # Folder relative to site-packages to install to ) -# TODO: comment in if we include tests subdir('tests') subdir('unsupervised') subdir('manifold') diff --git a/sktree/tree/tests/test_tree.py b/sktree/tree/tests/test_tree.py index 5bcd23c23..480800bcd 100644 --- a/sktree/tree/tests/test_tree.py +++ b/sktree/tree/tests/test_tree.py @@ -223,16 +223,11 @@ def test_pickle_splitters(): ] ) def test_sklearn_compatible_estimator(estimator, check): - # TODO: remove when we can replicate the CI error... - if isinstance( - estimator, (PatchObliqueDecisionTreeClassifier, ExtraObliqueDecisionTreeClassifier) - ) and check.func.__name__ in ["check_fit_score_takes_y"]: - pytest.skip() check(estimator) @pytest.mark.parametrize("Tree", CLF_TREES.values()) -def test_oblique_tree_sampling(Tree, random_state=0): +def test_oblique_tree_sampling_iris(Tree, random_state=0): """Test Oblique Decision Trees. Oblique trees can sample more candidate splits than @@ -253,9 +248,11 @@ def test_oblique_tree_sampling(Tree, random_state=0): tree_rc = Tree(random_state=random_state, max_features=n_features * 2) ri_cv_scores = cross_val_score(tree_ri, X, y, scoring="accuracy", cv=10, error_score="raise") rc_cv_scores = cross_val_score(tree_rc, X, y, scoring="accuracy", cv=10, error_score="raise") - assert rc_cv_scores.mean() > ri_cv_scores.mean() - assert rc_cv_scores.std() < ri_cv_scores.std() - assert rc_cv_scores.mean() > 0.91 + assert ( + rc_cv_scores.mean() > ri_cv_scores.mean() + ), f"{rc_cv_scores.mean()} <= {ri_cv_scores.mean()}" + assert rc_cv_scores.std() < ri_cv_scores.std(), f"{rc_cv_scores.std()} >= {ri_cv_scores.std()}" + assert rc_cv_scores.mean() > 0.91, f"{rc_cv_scores.mean()} <= 0.91" @pytest.mark.parametrize("Tree", OBLIQUE_TREES.values()) @@ -483,7 +480,9 @@ def test_diabetes_underfit(name, Tree, criterion, max_depth, metric, max_loss): reg = Tree(criterion=criterion, max_depth=max_depth, max_features=10, random_state=1234) reg.fit(diabetes.data, diabetes.target) loss = metric(diabetes.target, reg.predict(diabetes.data)) - assert 0.0 <= loss < max_loss + assert ( + 0 <= loss < max_loss + ), f"Failed with {name}, criterion = {criterion} and loss = {loss} vs {max_loss}" def test_numerical_stability(): @@ -521,4 +520,7 @@ def test_balance_property(criterion, Tree): X, y = diabetes.data, diabetes.target reg = Tree(criterion=criterion) reg.fit(X, y) - assert np.sum(reg.predict(X)) == pytest.approx(np.sum(y)) + y_pred_sum = np.sum(reg.predict(X)) + assert y_pred_sum == pytest.approx( + np.sum(y) + ), f"Failed with {criterion} and {Tree}: {y_pred_sum} != {np.sum(y)}" diff --git a/sktree/tree/tests/test_utils.py b/sktree/tree/tests/test_utils.py index 3ba69adbc..a0842d6ec 100644 --- a/sktree/tree/tests/test_utils.py +++ b/sktree/tree/tests/test_utils.py @@ -122,7 +122,7 @@ def test_best_patch_splitter_contiguous(): top_left_patch_seed, patch_size, patch_dims = splitter.sample_top_left_seed_cpdef() # sample a bunch of projection vectors - proj_vec = splitter.sample_projection_vector( + proj_vec = splitter.sample_projection_vector_py( proj_i, patch_size, top_left_patch_seed, patch_dims ) # print(proj_vec.reshape(data_dims)) diff --git a/sktree/tree/unsupervised/_unsup_criterion.pyx b/sktree/tree/unsupervised/_unsup_criterion.pyx index 87874b800..4ca6cea9e 100644 --- a/sktree/tree/unsupervised/_unsup_criterion.pyx +++ b/sktree/tree/unsupervised/_unsup_criterion.pyx @@ -288,7 +288,7 @@ cdef class TwoMeans(UnsupervisedCriterion): to compute a weighted sample mean and weighted sample variance. The weighted sample variance has two approaches to compute an unbiased estimate. The first is using "frequencies" and the second is using "reliability" weights. Currently, - we have impmlemented the frequencies approach. + we have implemented the frequencies approach. # TODO: implement reliability weighting diff --git a/sktree/tree/unsupervised/_unsup_oblique_splitter.pxd b/sktree/tree/unsupervised/_unsup_oblique_splitter.pxd index 07a7049bb..17a7956fd 100644 --- a/sktree/tree/unsupervised/_unsup_oblique_splitter.pxd +++ b/sktree/tree/unsupervised/_unsup_oblique_splitter.pxd @@ -1,6 +1,7 @@ import numpy as np cimport numpy as cnp +from libcpp.unordered_map cimport unordered_map from libcpp.vector cimport vector from ..._lib.sklearn.tree._splitter cimport SplitRecord @@ -44,16 +45,30 @@ cdef class UnsupervisedObliqueSplitter(UnsupervisedSplitter): cdef vector[vector[SIZE_t]] proj_mat_indices # nonzero indices of sparse proj_mat matrix cdef SIZE_t[::1] indices_to_sample # an array of indices to sample of size mtry X n_features + # keep a hashmap of every projection vector indices sampled + cdef unordered_map[size_t, bint] proj_vec_hash + + cdef unordered_map[SIZE_t, DTYPE_t] min_val_map + cdef unordered_map[SIZE_t, DTYPE_t] max_val_map + cdef unordered_map[SIZE_t, bint] constant_col_map + # All oblique splitters (i.e. non-axis aligned splitters) require a # function to sample a projection matrix that is applied to the feature matrix # to quickly obtain the sampled projections for candidate splits. - cdef void sample_proj_mat(self, - vector[vector[DTYPE_t]]& proj_mat_weights, - vector[vector[SIZE_t]]& proj_mat_indices) noexcept nogil + cdef void sample_proj_mat( + self, + vector[vector[DTYPE_t]]& proj_mat_weights, + vector[vector[SIZE_t]]& proj_mat_indices, + SIZE_t n_known_constants + ) noexcept nogil # Redefined here since the new logic requires calling sample_proj_mat - cdef int node_reset(self, SIZE_t start, SIZE_t end, - double* weighted_n_node_samples) except -1 nogil + cdef int node_reset( + self, + SIZE_t start, + SIZE_t end, + double* weighted_n_node_samples + ) except -1 nogil cdef int node_split( self, @@ -77,5 +92,6 @@ cdef class UnsupervisedObliqueSplitter(UnsupervisedSplitter): const SIZE_t[:] samples, DTYPE_t[:] feature_values, vector[DTYPE_t]* proj_vec_weights, # weights of the vector (max_features,) - vector[SIZE_t]* proj_vec_indices # indices of the features (max_features,) + vector[SIZE_t]* proj_vec_indices, # indices of the features (max_features,) + SIZE_t* n_known_constants ) noexcept nogil diff --git a/sktree/tree/unsupervised/_unsup_oblique_splitter.pyx b/sktree/tree/unsupervised/_unsup_oblique_splitter.pyx index e65612095..3bcb571ca 100644 --- a/sktree/tree/unsupervised/_unsup_oblique_splitter.pyx +++ b/sktree/tree/unsupervised/_unsup_oblique_splitter.pyx @@ -155,7 +155,8 @@ cdef class UnsupervisedObliqueSplitter(UnsupervisedSplitter): cdef void sample_proj_mat(self, vector[vector[DTYPE_t]]& proj_mat_weights, - vector[vector[SIZE_t]]& proj_mat_indices) noexcept nogil: + vector[vector[SIZE_t]]& proj_mat_indices, + SIZE_t n_known_constants) noexcept nogil: """ Sample the projection vector. This is a placeholder method. @@ -166,14 +167,15 @@ cdef class UnsupervisedObliqueSplitter(UnsupervisedSplitter): """Get size of a pointer to record for ObliqueSplitter.""" return sizeof(ObliqueSplitRecord) - cdef inline void compute_features_over_samples( + cdef void compute_features_over_samples( self, SIZE_t start, SIZE_t end, const SIZE_t[:] samples, DTYPE_t[:] feature_values, vector[DTYPE_t]* proj_vec_weights, # weights of the vector (max_features,) - vector[SIZE_t]* proj_vec_indices # indices of the features (max_features,) + vector[SIZE_t]* proj_vec_indices, # indices of the features (max_features,) + SIZE_t* n_known_constants ) noexcept nogil: """Compute the feature values for the samples[start:end] range. @@ -196,23 +198,58 @@ cdef class UnsupervisedObliqueSplitter(UnsupervisedSplitter): feature_values[idx] = 0.0 feature_values[idx] += self.X[samples[idx], col_idx] * col_weight + # keep track of the min/max of X[samples[:], col_idx] + if (self.X[samples[idx], col_idx] < self.min_val_map[col_idx] or + self.min_val_map.find(col_idx) == self.min_val_map.end()): + self.min_val_map[col_idx] = self.X[samples[idx], col_idx] + if (self.X[samples[idx], col_idx] > self.max_val_map[col_idx] or + self.max_val_map.find(col_idx) == self.max_val_map.end()): + self.max_val_map[col_idx] = self.X[samples[idx], col_idx] + + if self.max_val_map[col_idx] <= self.min_val_map[col_idx] + FEATURE_THRESHOLD: + self.constant_features[col_idx] = 1 + + # move features pointer around to make sure we keep track of the constant features + self.features[col_idx], self.features[n_known_constants[0]] = \ + self.features[n_known_constants[0]], self.features[col_idx] + + # increment the number of known constants + n_known_constants[0] += 1 cdef class BestObliqueUnsupervisedSplitter(UnsupervisedObliqueSplitter): # NOTE: vectors are passed by value, so & is needed to pass by reference cdef void sample_proj_mat( self, vector[vector[DTYPE_t]]& proj_mat_weights, - vector[vector[SIZE_t]]& proj_mat_indices + vector[vector[SIZE_t]]& proj_mat_indices, + SIZE_t n_known_constants, ) noexcept nogil: - """ - Sparse Oblique Projection matrix. + """Sample oblique projection matrix. + Randomly sample features to put in randomly sampled projection vectors - weight = 1 or -1 with probability 0.5 + weight = 1 or -1 with probability 0.5. + + Note: vectors are passed by value, so & is needed to pass by reference. + + Parameters + ---------- + proj_mat_weights : vector of vectors reference + The memory address of projection matrix non-zero weights. + proj_mat_indices : vector of vectors reference + The memory address of projection matrix non-zero indices. + + Notes + ----- + Note that grid_size must be larger than or equal to n_non_zeros because + it is assumed ``feature_combinations`` is forced to be smaller than + ``n_features`` before instantiating an oblique splitter. """ + cdef SIZE_t n_features = self.n_features cdef SIZE_t n_non_zeros = self.n_non_zeros cdef UINT32_t* random_state = &self.rand_r_state + cdef SIZE_t col_idx cdef int i, feat_i, proj_i, rand_vec_index cdef DTYPE_t weight @@ -234,12 +271,13 @@ cdef class BestObliqueUnsupervisedSplitter(UnsupervisedObliqueSplitter): # get the projection index and feature index proj_i = rand_vec_index // n_features - feat_i = rand_vec_index % n_features + feat_i = rand_vec_index % (n_features - n_known_constants) + col_idx = self.features[feat_i + n_known_constants] # sample a random weight weight = 1 if (rand_int(0, 2, random_state) == 1) else -1 - proj_mat_indices[proj_i].push_back(feat_i) # Store index of nonzero + proj_mat_indices[proj_i].push_back(col_idx) # Store index of nonzero proj_mat_weights[proj_i].push_back(weight) # Store weight of nonzero cdef int node_split( @@ -263,6 +301,8 @@ cdef class BestObliqueUnsupervisedSplitter(UnsupervisedObliqueSplitter): cdef SIZE_t start = self.start cdef SIZE_t end = self.end + cdef SIZE_t n_known_constants = n_constant_features[0] + # pointer array to store feature values to split on cdef DTYPE_t[::1] feature_values = self.feature_values cdef SIZE_t max_features = self.max_features @@ -283,7 +323,7 @@ cdef class BestObliqueUnsupervisedSplitter(UnsupervisedObliqueSplitter): _init_split(&best_split, end) # Sample the projection matrix - self.sample_proj_mat(self.proj_mat_weights, self.proj_mat_indices) + self.sample_proj_mat(self.proj_mat_weights, self.proj_mat_indices, n_known_constants) # For every vector in the projection matrix for feat_i in range(max_features): @@ -305,7 +345,8 @@ cdef class BestObliqueUnsupervisedSplitter(UnsupervisedObliqueSplitter): samples, feature_values, &self.proj_mat_weights[feat_i], - &self.proj_mat_indices[feat_i] + &self.proj_mat_indices[feat_i], + &n_known_constants ) # Sort the samples