diff --git a/.github/workflows/emscripten.yml b/.github/workflows/emscripten.yml index 906c093a35232..8a2863fd062e2 100644 --- a/.github/workflows/emscripten.yml +++ b/.github/workflows/emscripten.yml @@ -73,8 +73,9 @@ jobs: SKLEARN_SKIP_OPENMP_TEST: "true" SKLEARN_SKIP_NETWORK_TESTS: 1 CIBW_TEST_REQUIRES: "pytest pandas" - # -s pytest argument is needed to avoid an issue in pytest output capturing with Pyodide - CIBW_TEST_COMMAND: "python -m pytest -sra --pyargs sklearn --durations 20 --showlocals" + # pytest -s argument is needed to avoid an issue in pytest output capturing with Pyodide + # -p no:cacheprovider is needed to avoid cleaning up pytest cache failing with PermissionError + CIBW_TEST_COMMAND: "python -m pytest -sra -p no:cacheprovider --pyargs sklearn --durations 20 --showlocals" - name: Upload wheel artifact uses: actions/upload-artifact@v7 diff --git a/sklearn/linear_model/_cd_fast.pyx b/sklearn/linear_model/_cd_fast.pyx index b3e9cb712ce21..d730049909c4a 100644 --- a/sklearn/linear_model/_cd_fast.pyx +++ b/sklearn/linear_model/_cd_fast.pyx @@ -607,7 +607,7 @@ cdef (floating, floating) gap_enet_sparse( bint positive, bint gap_smaller_eps, ) noexcept nogil: - """Compute dual gap for use in sparse_enet_coordinate_descent. + """Compute dual gap for use in enet_coordinate_descent_sparse. alpha > 0: formulation A of the duality gap alpha = 0 & beta > 0: formulation B of the duality gap @@ -705,7 +705,7 @@ cdef (floating, floating) gap_enet_sparse( return gap, dual_norm_XtA -def sparse_enet_coordinate_descent( +def enet_coordinate_descent_sparse( floating[::1] w, floating alpha, floating beta, diff --git a/sklearn/linear_model/_coordinate_descent.py b/sklearn/linear_model/_coordinate_descent.py index 98df190307ed6..1eaf274a871fa 100644 --- a/sklearn/linear_model/_coordinate_descent.py +++ b/sklearn/linear_model/_coordinate_descent.py @@ -5,6 +5,7 @@ import sys import warnings from abc import ABC, abstractmethod +from enum import Enum from functools import partial from numbers import Integral, Real @@ -52,6 +53,13 @@ ) +class CD_Algo(Enum): + ENET_CD = 0 + ENET_CD_GRAM = 1 + ENET_CD_SPARSE = 2 + ENET_CD_MULTITASK = 3 + + def _set_order(X, y, order="C"): """Change the order of X and y if necessary. @@ -781,12 +789,54 @@ def enet_path( X_indices = None X_indptr = None + if multi_output: + algo = CD_Algo.ENET_CD_MULTITASK + elif X_is_sparse: + algo = CD_Algo.ENET_CD_SPARSE + elif isinstance(precompute, np.ndarray): + algo = CD_Algo.ENET_CD_GRAM + # We expect precompute to be already Fortran ordered when bypassing checks + if check_input: + precompute = check_array(precompute, dtype=X.dtype.type, order="C") + else: # precompute is False + algo = CD_Algo.ENET_CD + + params = dict( + max_iter=max_iter, + tol=tol, + rng=rng, + random=random, + do_screening=do_screening, + early_stopping=early_stopping, + ) + for i, alpha in enumerate(alphas): # account for n_samples scaling in objectives between here and cd_fast l1_reg = alpha * l1_ratio * n_samples l2_reg = alpha * (1.0 - l1_ratio) * n_samples - if not multi_output and X_is_sparse: - model = cd_fast.sparse_enet_coordinate_descent( + if algo == CD_Algo.ENET_CD: + model = cd_fast.enet_coordinate_descent( + w=coef_, + alpha=l1_reg, + beta=l2_reg, + X=X, + y=y, + positive=positive, + **params, + ) + elif algo == CD_Algo.ENET_CD_GRAM: + model = cd_fast.enet_coordinate_descent_gram( + w=coef_, + alpha=l1_reg, + beta=l2_reg, + Q=precompute, # the gram matrix + q=Xy, + y=y, + positive=positive, + **params, + ) + elif algo == CD_Algo.ENET_CD_SPARSE: + model = cd_fast.enet_coordinate_descent_sparse( w=coef_, alpha=l1_reg, beta=l2_reg, @@ -796,15 +846,10 @@ def enet_path( y=y, sample_weight=sample_weight, X_mean=X_sparse_scaling, - max_iter=max_iter, - tol=tol, - rng=rng, - random=random, positive=positive, - do_screening=do_screening, - early_stopping=early_stopping, + **params, ) - elif multi_output: + elif algo == CD_Algo.ENET_CD_MULTITASK: model = cd_fast.enet_coordinate_descent_multi_task( W=coef_, alpha=l1_reg, @@ -817,53 +862,9 @@ def enet_path( Y=y, sample_weight=sample_weight, X_mean=X_sparse_scaling, - max_iter=max_iter, - tol=tol, - rng=rng, - random=random, - do_screening=do_screening, - early_stopping=early_stopping, - ) - elif isinstance(precompute, np.ndarray): - # We expect precompute to be already Fortran ordered when bypassing - # checks - if check_input: - precompute = check_array(precompute, dtype=X.dtype.type, order="C") - model = cd_fast.enet_coordinate_descent_gram( - coef_, - l1_reg, - l2_reg, - precompute, - Xy, - y, - max_iter, - tol, - rng, - random, - positive, - do_screening, - early_stopping, - ) - elif precompute is False: - model = cd_fast.enet_coordinate_descent( - coef_, - l1_reg, - l2_reg, - X, - y, - max_iter, - tol, - rng, - random, - positive, - do_screening, - early_stopping, - ) - else: - raise ValueError( - "Precompute should be one of True, False, 'auto' or array-like. Got %r" - % precompute + **params, ) + coef_, dual_gap_, eps_, n_iter_ = model coefs[..., i] = coef_ # we correct the scale of the returned dual gap, as the objective diff --git a/sklearn/linear_model/tests/test_coordinate_descent.py b/sklearn/linear_model/tests/test_coordinate_descent.py index cee527a8d1502..b978a70e8d2d1 100644 --- a/sklearn/linear_model/tests/test_coordinate_descent.py +++ b/sklearn/linear_model/tests/test_coordinate_descent.py @@ -141,7 +141,7 @@ def zc(): Xs = sparse_csc_type(X) for do_screening in [True, False]: coef_3 = zc() - cd_fast.sparse_enet_coordinate_descent( + cd_fast.enet_coordinate_descent_sparse( w=coef_3, alpha=alpha, X_data=Xs.data, @@ -196,7 +196,7 @@ def test_cython_solver_early_stopping(cd): elif cd == "sparse_enet": Xs = sparse.csc_matrix(X) cd_solve = partial( - cd_fast.sparse_enet_coordinate_descent, + cd_fast.enet_coordinate_descent_sparse, X_data=Xs.data, X_indices=Xs.indices, X_indptr=Xs.indptr, diff --git a/sklearn/linear_model/tests/test_sparse_coordinate_descent.py b/sklearn/linear_model/tests/test_sparse_coordinate_descent.py index 0e34e8b2db4c3..a060b5b4be1b4 100644 --- a/sklearn/linear_model/tests/test_sparse_coordinate_descent.py +++ b/sklearn/linear_model/tests/test_sparse_coordinate_descent.py @@ -371,7 +371,7 @@ def test_same_multiple_output_sparse_dense(coo_container): @pytest.mark.parametrize("csc_container", CSC_CONTAINERS) -def test_sparse_enet_coordinate_descent(csc_container): +def test_enet_coordinate_descent_sparse(csc_container): """Test that a warning is issued if model does not converge""" clf = Lasso( alpha=1e-10, fit_intercept=False, warm_start=True, max_iter=2, tol=1e-10 diff --git a/sklearn/tree/_partitioner.pxd b/sklearn/tree/_partitioner.pxd index d09f10323ccd9..d528d651dfedb 100644 --- a/sklearn/tree/_partitioner.pxd +++ b/sklearn/tree/_partitioner.pxd @@ -141,12 +141,64 @@ cdef class DensePartitioner: float32_t* min_feature_value_out, float32_t* max_feature_value_out, ) noexcept nogil - cdef void next_p( + cdef inline void next_p( self, intp_t* p_prev, intp_t* p, - bint missing_go_to_left - ) noexcept nogil + bint missing_go_to_left, + ) noexcept nogil: + """ + Compute the next p_prev and p for iterating over feature values. + + This method is used inside the best-split search function pass which starts + by setting p = start at the beginning of each search pass and calls + this method repeatedly with the same missing_go_to_left as for that pass. + The expected layout of self.feature_values[start:end] is: + - first pass (missing_go_to_left=False): after + sort_samples_and_feature_values(), non-missing values are sorted and + missing values are grouped at the right; + - second pass (missing_go_to_left=True): after + shift_missing_to_the_left(), missing values are grouped at the left. + + Given that layout, this method advances p to the next valid split + position while skipping ties up to FEATURE_THRESHOLD: + - if missing_go_to_left: iterate p in [start + n_missing + 1, end) + - otherwise: iterate p in [start, end - n_missing]. + The special case p == end - n_missing corresponds to "all non-missing + values on the left and all missing values on the right". The next + call then sets p to end to terminate the search loop. + """ + cdef intp_t end_non_missing = ( + self.end if missing_go_to_left + else self.end - self.n_missing) + + # First pass special end marker: include "all non-missing left, all missing right" + if p[0] == end_non_missing and not missing_go_to_left: + p[0] = self.end + p_prev[0] = self.end + return + + # Second pass starts with missing on the left; jump to first non-missing + if missing_go_to_left and p[0] == self.start: + p[0] = self.start + self.n_missing + + # Move to next candidate split position + p[0] += 1 + + if self.n_categories_current > 0: + while ( + p[0] < end_non_missing and + self.feature_values[p[0]] == self.feature_values[p[0] - 1] + ): + p[0] += 1 + else: + while ( + p[0] < end_non_missing and + self.feature_values[p[0]] <= self.feature_values[p[0] - 1] + FEATURE_THRESHOLD + ): + p[0] += 1 + + p_prev[0] = p[0] - 1 cdef intp_t partition_samples( self, float64_t current_threshold, @@ -157,12 +209,33 @@ cdef class DensePartitioner: const SplitRecord* best_split, ) noexcept nogil - cdef void cat_position_to_split_bitset( + cdef inline void cat_position_to_split_bitset( self, intp_t position, bint missing_go_to_left, - BITSET_DTYPE_C left_cat_bitset - ) noexcept nogil + BITSET_DTYPE_C left_cat_bitset, + ) noexcept nogil: + """Convert a categorical split position into a bitset.""" + cdef: + intp_t n_left_non_missing = position - self.start + intp_t offset = 0 + intp_t r + intp_t c + + if missing_go_to_left: + n_left_non_missing -= self.n_missing + + init_bitset(left_cat_bitset) + + if n_left_non_missing <= 0: + return + + for r in range(self.n_categories_current): + c = self.sorted_cat[r] + set_bitset(left_cat_bitset, c) + offset += self.counts[c] + if offset >= n_left_non_missing: + break cdef void sort_categories( self, intp_t nc @@ -217,12 +290,34 @@ cdef class SparsePartitioner: float32_t* min_feature_value_out, float32_t* max_feature_value_out, ) noexcept nogil - cdef void next_p( + cdef inline void next_p( self, intp_t* p_prev, intp_t* p, - bint missing_go_to_left - ) noexcept nogil + bint missing_go_to_left, + ) noexcept nogil: + """Compute the next p_prev and p for iterating over feature values. + + The missing_go_to_left argument is ignored for sparse data because + sparse partitioning does not currently support missing values. + """ + cdef intp_t p_next + + if p[0] + 1 != self.end_negative: + p_next = p[0] + 1 + else: + p_next = self.start_positive + + while (p_next < self.end and + self.feature_values[p_next] <= self.feature_values[p[0]] + FEATURE_THRESHOLD): + p[0] = p_next + if p[0] + 1 != self.end_negative: + p_next = p[0] + 1 + else: + p_next = self.start_positive + + p_prev[0] = p[0] + p[0] = p_next cdef intp_t partition_samples( self, float64_t current_threshold, @@ -233,12 +328,14 @@ cdef class SparsePartitioner: const SplitRecord* best_split, ) noexcept nogil - cdef void cat_position_to_split_bitset( + cdef inline void cat_position_to_split_bitset( self, intp_t position, bint missing_go_to_left, - BITSET_DTYPE_C left_cat_bitset - ) noexcept nogil + BITSET_DTYPE_C left_cat_bitset, + ) noexcept nogil: + # Sparse categorical features are rejected before split search. + init_bitset(left_cat_bitset) cdef void extract_nnz( self, intp_t feature diff --git a/sklearn/tree/_partitioner.pyx b/sklearn/tree/_partitioner.pyx index 13e573431aed0..86d788d197466 100644 --- a/sklearn/tree/_partitioner.pyx +++ b/sklearn/tree/_partitioner.pyx @@ -19,7 +19,6 @@ import numpy as np cimport numpy as cnp cnp.import_array() from scipy.sparse import issparse -from sklearn.utils._bitset cimport BITSET_DTYPE_C, init_bitset from sklearn.tree._utils cimport goes_left, MAX_NUM_CATEGORIES from sklearn.tree._splitter cimport SplitRecord from sklearn.utils._sorting cimport simultaneous_sort @@ -299,65 +298,6 @@ cdef class DensePartitioner: max_feature_value_out[0] = max_feature_value self.n_missing = n_missing - cdef inline void next_p( - self, - intp_t* p_prev, - intp_t* p, - bint missing_go_to_left, - ) noexcept nogil: - """ - Compute the next p_prev and p for iterating over feature values. - - This method is used inside the best-split search function pass which starts - by setting p = start at the beginning of each search pass and calls - this method repeatedly with the same missing_go_to_left as for that pass. - The expected layout of self.feature_values[start:end] is: - - first pass (missing_go_to_left=False): after - sort_samples_and_feature_values(), non-missing values are sorted and - missing values are grouped at the right; - - second pass (missing_go_to_left=True): after - shift_missing_to_the_left(), missing values are grouped at the left. - - Given that layout, this method advances p to the next valid split - position while skipping ties up to FEATURE_THRESHOLD: - - if missing_go_to_left: iterate p in [start + n_missing + 1, end) - - otherwise: iterate p in [start, end - n_missing]. - The special case p == end - n_missing corresponds to "all non-missing - values on the left and all missing values on the right". The next - call then sets p to end to terminate the search loop. - """ - cdef intp_t end_non_missing = ( - self.end if missing_go_to_left - else self.end - self.n_missing) - - # First pass special end marker: include "all non-missing left, all missing right" - if p[0] == end_non_missing and not missing_go_to_left: - p[0] = self.end - p_prev[0] = self.end - return - - # Second pass starts with missing on the left; jump to first non-missing - if missing_go_to_left and p[0] == self.start: - p[0] = self.start + self.n_missing - - # Move to next candidate split position - p[0] += 1 - - if self.n_categories_current > 0: - while ( - p[0] < end_non_missing and - self.feature_values[p[0]] == self.feature_values[p[0] - 1] - ): - p[0] += 1 - else: - while ( - p[0] < end_non_missing and - self.feature_values[p[0]] <= self.feature_values[p[0] - 1] + FEATURE_THRESHOLD - ): - p[0] += 1 - - p_prev[0] = p[0] - 1 - cdef inline intp_t partition_samples( self, float64_t threshold, @@ -426,34 +366,6 @@ cdef class DensePartitioner: samples[partition_start], samples[partition_end] = ( samples[partition_end], samples[partition_start]) - cdef inline void cat_position_to_split_bitset( - self, - intp_t position, - bint missing_go_to_left, - BITSET_DTYPE_C left_cat_bitset, - ) noexcept nogil: - """Convert a categorical split position into a bitset.""" - cdef: - intp_t n_left_non_missing = position - self.start - intp_t offset = 0 - intp_t r - intp_t c - - if missing_go_to_left: - n_left_non_missing -= self.n_missing - - init_bitset(left_cat_bitset) - - if n_left_non_missing <= 0: - return - - for r in range(self.n_categories_current): - c = self.sorted_cat[r] - set_bitset(left_cat_bitset, c) - offset += self.counts[c] - if offset >= n_left_non_missing: - break - @final cdef class SparsePartitioner: @@ -598,35 +510,6 @@ cdef class SparsePartitioner: min_feature_value_out[0] = min_feature_value max_feature_value_out[0] = max_feature_value - cdef inline void next_p( - self, - intp_t* p_prev, - intp_t* p, - bint missing_go_to_left, - ) noexcept nogil: - """Compute the next p_prev and p for iterating over feature values. - - The missing_go_to_left argument is ignored for sparse data because - sparse partitioning does not currently support missing values. - """ - cdef intp_t p_next - - if p[0] + 1 != self.end_negative: - p_next = p[0] + 1 - else: - p_next = self.start_positive - - while (p_next < self.end and - self.feature_values[p_next] <= self.feature_values[p[0]] + FEATURE_THRESHOLD): - p[0] = p_next - if p[0] + 1 != self.end_negative: - p_next = p[0] + 1 - else: - p_next = self.start_positive - - p_prev[0] = p[0] - p[0] = p_next - cdef inline intp_t partition_samples( self, float64_t current_threshold, @@ -682,15 +565,6 @@ cdef class SparsePartitioner: return partition_end - cdef inline void cat_position_to_split_bitset( - self, - intp_t position, - bint missing_go_to_left, - BITSET_DTYPE_C left_cat_bitset, - ) noexcept nogil: - # Sparse categorical features are rejected before split search. - init_bitset(left_cat_bitset) - cdef inline void extract_nnz(self, intp_t feature) noexcept nogil: """Extract and partition values for a given feature. diff --git a/sklearn/tree/_utils.pxd b/sklearn/tree/_utils.pxd index f37593169004c..81d7bad628fcb 100644 --- a/sklearn/tree/_utils.pxd +++ b/sklearn/tree/_utils.pxd @@ -2,11 +2,15 @@ # SPDX-License-Identifier: BSD-3-Clause # See _utils.pyx for details. +from libc.math cimport isnan +from libc.math cimport log as ln +from libc.stdlib cimport realloc + cimport numpy as cnp from sklearn.neighbors._quad_tree cimport Cell from sklearn.utils._typedefs cimport float32_t, float64_t, intp_t, uint8_t, int32_t, uint32_t, uint64_t -from sklearn.utils._bitset cimport BITSET_DTYPE_C, BITSET_INNER_DTYPE_C, N_BITSETS - +from sklearn.utils._bitset cimport BITSET_DTYPE_C, BITSET_INNER_DTYPE_C, N_BITSETS, in_bitset +from sklearn.utils._random cimport our_rand_r cdef enum: MAX_NUM_CATEGORIES = N_BITSETS @@ -31,13 +35,19 @@ cdef struct Node: uint8_t missing_go_to_left # Whether features have missing values -cdef bint goes_left( +cdef inline bint goes_left( float64_t threshold, const BITSET_INNER_DTYPE_C* left_cat_bitset, bint missing_go_to_left, bint is_categorical, float32_t value, -) noexcept nogil +) noexcept nogil: + if isnan(value): + return missing_go_to_left + elif is_categorical: + return in_bitset(left_cat_bitset, value) + else: + return value <= threshold cdef enum: @@ -67,22 +77,32 @@ ctypedef fused realloc_ptr: (Cell*) (Node**) -cdef int safe_realloc(realloc_ptr* p, size_t nelems) except -1 nogil +cdef int safe_realloc(realloc_ptr* p, size_t nelems) except -1 nogil -cdef cnp.ndarray sizet_ptr_to_ndarray(intp_t* data, intp_t size) +cdef inline cnp.ndarray sizet_ptr_to_ndarray(intp_t* data, intp_t size): + """Return copied data as 1D numpy array of intp's.""" + cdef cnp.npy_intp shape[1] + shape[0] = size + return cnp.PyArray_SimpleNewFromData(1, shape, cnp.NPY_INTP, data).copy() -cdef intp_t rand_int(intp_t low, intp_t high, - uint32_t* random_state) noexcept nogil +cdef inline intp_t rand_int(intp_t low, intp_t high, + uint32_t* random_state) noexcept nogil: + """Generate a random integer in [low; end).""" + return low + our_rand_r(random_state) % (high - low) -cdef float64_t rand_uniform(float64_t low, float64_t high, - uint32_t* random_state) noexcept nogil +cdef inline float64_t rand_uniform(float64_t low, float64_t high, + uint32_t* random_state) noexcept nogil: + """Generate a random float64_t in [low; high).""" + return ((high - low) * our_rand_r(random_state) / + RAND_R_MAX) + low -cdef float64_t log(float64_t x) noexcept nogil +cdef inline float64_t log(float64_t x) noexcept nogil: + return ln(x) / ln(2.0) cdef class WeightedFenwickTree: cdef intp_t size # number of leaves (ranks) diff --git a/sklearn/tree/_utils.pyx b/sklearn/tree/_utils.pyx index c3aa0b39e190f..9c9b0fa9ec72a 100644 --- a/sklearn/tree/_utils.pyx +++ b/sklearn/tree/_utils.pyx @@ -2,30 +2,11 @@ # SPDX-License-Identifier: BSD-3-Clause from libc.stdlib cimport free -from libc.stdlib cimport realloc -from libc.math cimport log as ln -from libc.math cimport isnan from libc.string cimport memset cimport numpy as cnp cnp.import_array() -from sklearn.utils._bitset cimport BITSET_INNER_DTYPE_C, in_bitset -from sklearn.utils._random cimport our_rand_r - -cdef inline bint goes_left( - float64_t threshold, - const BITSET_INNER_DTYPE_C* left_cat_bitset, - bint missing_go_to_left, - bint is_categorical, - float32_t value, -) noexcept nogil: - if isnan(value): - return missing_go_to_left - elif is_categorical: - return in_bitset(left_cat_bitset, value) - else: - return value <= threshold # ============================================================================= # Helper functions @@ -57,30 +38,6 @@ def _realloc_test(): assert False -cdef inline cnp.ndarray sizet_ptr_to_ndarray(intp_t* data, intp_t size): - """Return copied data as 1D numpy array of intp's.""" - cdef cnp.npy_intp shape[1] - shape[0] = size - return cnp.PyArray_SimpleNewFromData(1, shape, cnp.NPY_INTP, data).copy() - - -cdef inline intp_t rand_int(intp_t low, intp_t high, - uint32_t* random_state) noexcept nogil: - """Generate a random integer in [low; end).""" - return low + our_rand_r(random_state) % (high - low) - - -cdef inline float64_t rand_uniform(float64_t low, float64_t high, - uint32_t* random_state) noexcept nogil: - """Generate a random float64_t in [low; high).""" - return ((high - low) * our_rand_r(random_state) / - RAND_R_MAX) + low - - -cdef inline float64_t log(float64_t x) noexcept nogil: - return ln(x) / ln(2.0) - - cdef class WeightedFenwickTree: """ Fenwick tree (Binary Indexed Tree) specialized for maintaining: diff --git a/sklearn/utils/_heap.pxd b/sklearn/utils/_heap.pxd index 44293d5c2ef62..7583da88acf66 100644 --- a/sklearn/utils/_heap.pxd +++ b/sklearn/utils/_heap.pxd @@ -5,10 +5,83 @@ from cython cimport floating from sklearn.utils._typedefs cimport intp_t -cdef int heap_push( +cdef inline int heap_push( floating* values, intp_t* indices, intp_t size, floating val, intp_t val_idx, -) noexcept nogil +) noexcept nogil: + """Push a tuple (val, val_idx) onto a fixed-size max-heap. + + The max-heap is represented as a Structure of Arrays where: + - values is the array containing the data to construct the heap with + - indices is the array containing the indices (meta-data) of each value + + Notes + ----- + Arrays are manipulated via a pointer to there first element and their size + as to ease the processing of dynamically allocated buffers. + + For instance, in pseudo-code: + + values = [1.2, 0.4, 0.1], + indices = [42, 1, 5], + heap_push( + values=values, + indices=indices, + size=3, + val=0.2, + val_idx=4, + ) + + will modify values and indices inplace, giving at the end of the call: + + values == [0.4, 0.2, 0.1] + indices == [1, 4, 5] + + """ + cdef: + intp_t current_idx, left_child_idx, right_child_idx, swap_idx + + # Check if val should be in heap + if val >= values[0]: + return 0 + + # Insert val at position zero + values[0] = val + indices[0] = val_idx + + # Descend the heap, swapping values until the max heap criterion is met + current_idx = 0 + while True: + left_child_idx = 2 * current_idx + 1 + right_child_idx = left_child_idx + 1 + + if left_child_idx >= size: + break + elif right_child_idx >= size: + if values[left_child_idx] > val: + swap_idx = left_child_idx + else: + break + elif values[left_child_idx] >= values[right_child_idx]: + if val < values[left_child_idx]: + swap_idx = left_child_idx + else: + break + else: + if val < values[right_child_idx]: + swap_idx = right_child_idx + else: + break + + values[current_idx] = values[swap_idx] + indices[current_idx] = indices[swap_idx] + + current_idx = swap_idx + + values[current_idx] = val + indices[current_idx] = val_idx + + return 0 diff --git a/sklearn/utils/_heap.pyx b/sklearn/utils/_heap.pyx deleted file mode 100644 index 2e39118d10a7c..0000000000000 --- a/sklearn/utils/_heap.pyx +++ /dev/null @@ -1,85 +0,0 @@ -from cython cimport floating - -from sklearn.utils._typedefs cimport intp_t - - -cdef inline int heap_push( - floating* values, - intp_t* indices, - intp_t size, - floating val, - intp_t val_idx, -) noexcept nogil: - """Push a tuple (val, val_idx) onto a fixed-size max-heap. - - The max-heap is represented as a Structure of Arrays where: - - values is the array containing the data to construct the heap with - - indices is the array containing the indices (meta-data) of each value - - Notes - ----- - Arrays are manipulated via a pointer to there first element and their size - as to ease the processing of dynamically allocated buffers. - - For instance, in pseudo-code: - - values = [1.2, 0.4, 0.1], - indices = [42, 1, 5], - heap_push( - values=values, - indices=indices, - size=3, - val=0.2, - val_idx=4, - ) - - will modify values and indices inplace, giving at the end of the call: - - values == [0.4, 0.2, 0.1] - indices == [1, 4, 5] - - """ - cdef: - intp_t current_idx, left_child_idx, right_child_idx, swap_idx - - # Check if val should be in heap - if val >= values[0]: - return 0 - - # Insert val at position zero - values[0] = val - indices[0] = val_idx - - # Descend the heap, swapping values until the max heap criterion is met - current_idx = 0 - while True: - left_child_idx = 2 * current_idx + 1 - right_child_idx = left_child_idx + 1 - - if left_child_idx >= size: - break - elif right_child_idx >= size: - if values[left_child_idx] > val: - swap_idx = left_child_idx - else: - break - elif values[left_child_idx] >= values[right_child_idx]: - if val < values[left_child_idx]: - swap_idx = left_child_idx - else: - break - else: - if val < values[right_child_idx]: - swap_idx = right_child_idx - else: - break - - values[current_idx] = values[swap_idx] - indices[current_idx] = indices[swap_idx] - - current_idx = swap_idx - - values[current_idx] = val - indices[current_idx] = val_idx - - return 0 diff --git a/sklearn/utils/meson.build b/sklearn/utils/meson.build index 71b98c088d4c1..dd12e74baeb33 100644 --- a/sklearn/utils/meson.build +++ b/sklearn/utils/meson.build @@ -48,7 +48,6 @@ utils_extension_metadata = { '_openmp_helpers': {'sources': [cython_gen.process('_openmp_helpers.pyx')], 'dependencies': [openmp_dep]}, '_random': {'sources': [cython_gen.process('_random.pyx')]}, '_typedefs': {'sources': [cython_gen.process('_typedefs.pyx')]}, - '_heap': {'sources': [cython_gen.process('_heap.pyx')]}, '_sorting': {'sources': [cython_gen.process('_sorting.pyx')]}, '_vector_sentinel': {'sources': [cython_gen_cpp.process('_vector_sentinel.pyx')],