Aluode/PerceptionLabPortable
0
1"""Partition samples in the construction of a tree.
2
3This module contains the algorithms for moving sample indices to
4the left and right child node given a split determined by the
5splitting algorithm in `_splitter.pyx`.
6
7Partitioning is done in a way that is efficient for both dense data,
8and sparse data stored in a Compressed Sparse Column (CSC) format.
9"""
10# Authors: The scikit-learn developers
11# SPDX-License-Identifier: BSD-3-Clause
12
13from cython cimport final
14from libc.math cimport isnan, log2
15from libc.stdlib cimport qsort
16from libc.string cimport memcpy
17
18import numpy as np
19from scipy.sparse import issparse
20
21
22# Constant to switch between algorithm non zero value extract algorithm
23# in SparsePartitioner
24cdef float32_t EXTRACT_NNZ_SWITCH = 0.1
25
26# Allow for 32 bit float comparisons
27cdef float32_t INFINITY_32t = np.inf
28
29
30@final
31cdef class DensePartitioner:
32 """Partitioner specialized for dense data.
33
34 Note that this partitioner is agnostic to the splitting strategy (best vs. random).
35 """
36 def __init__(
37 self,
38 const float32_t[:, :] X,
39 intp_t[::1] samples,
40 float32_t[::1] feature_values,
41 const uint8_t[::1] missing_values_in_feature_mask,
42 ):
43 self.X = X
44 self.samples = samples
45 self.feature_values = feature_values
46 self.missing_values_in_feature_mask = missing_values_in_feature_mask
47
48 cdef inline void init_node_split(self, intp_t start, intp_t end) noexcept nogil:
49 """Initialize splitter at the beginning of node_split."""
50 self.start = start
51 self.end = end
52 self.n_missing = 0
53
54 cdef inline void sort_samples_and_feature_values(
55 self, intp_t current_feature
56 ) noexcept nogil:
57 """Simultaneously sort based on the feature_values.
58
59 Missing values are stored at the end of feature_values.
60 The number of missing values observed in feature_values is stored
61 in self.n_missing.
62 """
63 cdef:
64 intp_t i, current_end
65 float32_t[::1] feature_values = self.feature_values
66 const float32_t[:, :] X = self.X
67 intp_t[::1] samples = self.samples
68 intp_t n_missing = 0
69 const uint8_t[::1] missing_values_in_feature_mask = self.missing_values_in_feature_mask
70
71 # Sort samples along that feature; by copying the values into an array and
72 # sorting the array in a manner which utilizes the cache more effectively.
73 if missing_values_in_feature_mask is not None and missing_values_in_feature_mask[current_feature]:
74 i, current_end = self.start, self.end - 1
75 # Missing values are placed at the end and do not participate in the sorting.
76 while i <= current_end:
77 # Finds the right-most value that is not missing so that
78 # it can be swapped with missing values at its left.
79 if isnan(X[samples[current_end], current_feature]):
80 n_missing += 1
81 current_end -= 1
82 continue
83
84 # X[samples[current_end], current_feature] is a non-missing value
85 if isnan(X[samples[i], current_feature]):
86 samples[i], samples[current_end] = samples[current_end], samples[i]
87 n_missing += 1
88 current_end -= 1
89
90 feature_values[i] = X[samples[i], current_feature]
91 i += 1
92 else:
93 # When there are no missing values, we only need to copy the data into
94 # feature_values
95 for i in range(self.start, self.end):
96 feature_values[i] = X[samples[i], current_feature]
97
98 sort(&feature_values[self.start], &samples[self.start], self.end - self.start - n_missing)
99 self.n_missing = n_missing
100
101 cdef inline void find_min_max(
102 self,
103 intp_t current_feature,
104 float32_t* min_feature_value_out,
105 float32_t* max_feature_value_out,
106 ) noexcept nogil:
107 """Find the minimum and maximum value for current_feature.
108
109 Missing values are stored at the end of feature_values. The number of missing
110 values observed in feature_values is stored in self.n_missing.
111 """
112 cdef:
113 intp_t p, current_end
114 float32_t current_feature_value
115 const float32_t[:, :] X = self.X
116 intp_t[::1] samples = self.samples
117 float32_t min_feature_value = INFINITY_32t
118 float32_t max_feature_value = -INFINITY_32t
119 float32_t[::1] feature_values = self.feature_values
120 intp_t n_missing = 0
121 const uint8_t[::1] missing_values_in_feature_mask = self.missing_values_in_feature_mask
122
123 # We are copying the values into an array and finding min/max of the array in
124 # a manner which utilizes the cache more effectively. We need to also count
125 # the number of missing-values there are.
126 if missing_values_in_feature_mask is not None and missing_values_in_feature_mask[current_feature]:
127 p, current_end = self.start, self.end - 1
128 # Missing values are placed at the end and do not participate in the
129 # min/max calculation.
130 while p <= current_end:
131 # Finds the right-most value that is not missing so that
132 # it can be swapped with missing values towards its left.
133 if isnan(X[samples[current_end], current_feature]):
134 n_missing += 1
135 current_end -= 1
136 continue
137
138 # X[samples[current_end], current_feature] is a non-missing value
139 if isnan(X[samples[p], current_feature]):
140 samples[p], samples[current_end] = samples[current_end], samples[p]
141 n_missing += 1
142 current_end -= 1
143
144 current_feature_value = X[samples[p], current_feature]
145 feature_values[p] = current_feature_value
146 if current_feature_value < min_feature_value:
147 min_feature_value = current_feature_value
148 elif current_feature_value > max_feature_value:
149 max_feature_value = current_feature_value
150 p += 1
151 else:
152 min_feature_value = X[samples[self.start], current_feature]
153 max_feature_value = min_feature_value
154
155 feature_values[self.start] = min_feature_value
156 for p in range(self.start + 1, self.end):
157 current_feature_value = X[samples[p], current_feature]
158 feature_values[p] = current_feature_value
159
160 if current_feature_value < min_feature_value:
161 min_feature_value = current_feature_value
162 elif current_feature_value > max_feature_value:
163 max_feature_value = current_feature_value
164
165 min_feature_value_out[0] = min_feature_value
166 max_feature_value_out[0] = max_feature_value
167 self.n_missing = n_missing
168
169 cdef inline void next_p(self, intp_t* p_prev, intp_t* p) noexcept nogil:
170 """Compute the next p_prev and p for iterating over feature values.
171
172 The missing values are not included when iterating through the feature values.
173 """
174 cdef:
175 float32_t[::1] feature_values = self.feature_values
176 intp_t end_non_missing = self.end - self.n_missing
177
178 while (
179 p[0] + 1 < end_non_missing and
180 feature_values[p[0] + 1] <= feature_values[p[0]] + FEATURE_THRESHOLD
181 ):
182 p[0] += 1
183
184 p_prev[0] = p[0]
185
186 # By adding 1, we have
187 # (feature_values[p] >= end) or (feature_values[p] > feature_values[p - 1])
188 p[0] += 1
189
190 cdef inline intp_t partition_samples(
191 self,
192 float64_t current_threshold
193 ) noexcept nogil:
194 """Partition samples for feature_values at the current_threshold."""
195 cdef:
196 intp_t p = self.start
197 intp_t partition_end = self.end - self.n_missing
198 intp_t[::1] samples = self.samples
199 float32_t[::1] feature_values = self.feature_values
200
201 while p < partition_end:
202 if feature_values[p] <= current_threshold:
203 p += 1
204 else:
205 partition_end -= 1
206
207 feature_values[p], feature_values[partition_end] = (
208 feature_values[partition_end], feature_values[p]
209 )
210 samples[p], samples[partition_end] = samples[partition_end], samples[p]
211
212 return partition_end
213
214 cdef inline void partition_samples_final(
215 self,
216 intp_t best_pos,
217 float64_t best_threshold,
218 intp_t best_feature,
219 intp_t best_n_missing,
220 ) noexcept nogil:
221 """Partition samples for X at the best_threshold and best_feature.
222
223 If missing values are present, this method partitions `samples`
224 so that the `best_n_missing` missing values' indices are in the
225 right-most end of `samples`, that is `samples[end_non_missing:end]`.
226 """
227 cdef:
228 # Local invariance: start <= p <= partition_end <= end
229 intp_t start = self.start
230 intp_t p = start
231 intp_t end = self.end - 1
232 intp_t partition_end = end - best_n_missing
233 intp_t[::1] samples = self.samples
234 const float32_t[:, :] X = self.X
235 float32_t current_value
236
237 if best_n_missing != 0:
238 # Move samples with missing values to the end while partitioning the
239 # non-missing samples
240 while p < partition_end:
241 # Keep samples with missing values at the end
242 if isnan(X[samples[end], best_feature]):
243 end -= 1
244 continue
245
246 # Swap sample with missing values with the sample at the end
247 current_value = X[samples[p], best_feature]
248 if isnan(current_value):
249 samples[p], samples[end] = samples[end], samples[p]
250 end -= 1
251
252 # The swapped sample at the end is always a non-missing value, so
253 # we can continue the algorithm without checking for missingness.
254 current_value = X[samples[p], best_feature]
255
256 # Partition the non-missing samples
257 if current_value <= best_threshold:
258 p += 1
259 else:
260 samples[p], samples[partition_end] = samples[partition_end], samples[p]
261 partition_end -= 1
262 else:
263 # Partitioning routine when there are no missing values
264 while p < partition_end:
265 if X[samples[p], best_feature] <= best_threshold:
266 p += 1
267 else:
268 samples[p], samples[partition_end] = samples[partition_end], samples[p]
269 partition_end -= 1
270
271
272@final
273cdef class SparsePartitioner:
274 """Partitioner specialized for sparse CSC data.
275
276 Note that this partitioner is agnostic to the splitting strategy (best vs. random).
277 """
278 def __init__(
279 self,
280 object X,
281 intp_t[::1] samples,
282 intp_t n_samples,
283 float32_t[::1] feature_values,
284 const uint8_t[::1] missing_values_in_feature_mask,
285 ):
286 if not (issparse(X) and X.format == "csc"):
287 raise ValueError("X should be in csc format")
288
289 self.samples = samples
290 self.feature_values = feature_values
291
292 # Initialize X
293 cdef intp_t n_total_samples = X.shape[0]
294
295 self.X_data = X.data
296 self.X_indices = X.indices
297 self.X_indptr = X.indptr
298 self.n_total_samples = n_total_samples
299
300 # Initialize auxiliary array used to perform split
301 self.index_to_samples = np.full(n_total_samples, fill_value=-1, dtype=np.intp)
302 self.sorted_samples = np.empty(n_samples, dtype=np.intp)
303
304 cdef intp_t p
305 for p in range(n_samples):
306 self.index_to_samples[samples[p]] = p
307
308 self.missing_values_in_feature_mask = missing_values_in_feature_mask
309
310 cdef inline void init_node_split(self, intp_t start, intp_t end) noexcept nogil:
311 """Initialize splitter at the beginning of node_split."""
312 self.start = start
313 self.end = end
314 self.is_samples_sorted = 0
315 self.n_missing = 0
316
317 cdef inline void sort_samples_and_feature_values(
318 self,
319 intp_t current_feature
320 ) noexcept nogil:
321 """Simultaneously sort based on the feature_values."""
322 cdef:
323 float32_t[::1] feature_values = self.feature_values
324 intp_t[::1] index_to_samples = self.index_to_samples
325 intp_t[::1] samples = self.samples
326
327 self.extract_nnz(current_feature)
328 # Sort the positive and negative parts of `feature_values`
329 sort(&feature_values[self.start], &samples[self.start], self.end_negative - self.start)
330 if self.start_positive < self.end:
331 sort(
332 &feature_values[self.start_positive],
333 &samples[self.start_positive],
334 self.end - self.start_positive
335 )
336
337 # Update index_to_samples to take into account the sort
338 for p in range(self.start, self.end_negative):
339 index_to_samples[samples[p]] = p
340 for p in range(self.start_positive, self.end):
341 index_to_samples[samples[p]] = p
342
343 # Add one or two zeros in feature_values, if there is any
344 if self.end_negative < self.start_positive:
345 self.start_positive -= 1
346 feature_values[self.start_positive] = 0.
347
348 if self.end_negative != self.start_positive:
349 feature_values[self.end_negative] = 0.
350 self.end_negative += 1
351
352 # XXX: When sparse supports missing values, this should be set to the
353 # number of missing values for current_feature
354 self.n_missing = 0
355
356 cdef inline void find_min_max(
357 self,
358 intp_t current_feature,
359 float32_t* min_feature_value_out,
360 float32_t* max_feature_value_out,
361 ) noexcept nogil:
362 """Find the minimum and maximum value for current_feature."""
363 cdef:
364 intp_t p
365 float32_t current_feature_value, min_feature_value, max_feature_value
366 float32_t[::1] feature_values = self.feature_values
367
368 self.extract_nnz(current_feature)
369
370 if self.end_negative != self.start_positive:
371 # There is a zero
372 min_feature_value = 0
373 max_feature_value = 0
374 else:
375 min_feature_value = feature_values[self.start]
376 max_feature_value = min_feature_value
377
378 # Find min, max in feature_values[start:end_negative]
379 for p in range(self.start, self.end_negative):
380 current_feature_value = feature_values[p]
381
382 if current_feature_value < min_feature_value:
383 min_feature_value = current_feature_value
384 elif current_feature_value > max_feature_value:
385 max_feature_value = current_feature_value
386
387 # Update min, max given feature_values[start_positive:end]
388 for p in range(self.start_positive, self.end):
389 current_feature_value = feature_values[p]
390
391 if current_feature_value < min_feature_value:
392 min_feature_value = current_feature_value
393 elif current_feature_value > max_feature_value:
394 max_feature_value = current_feature_value
395
396 min_feature_value_out[0] = min_feature_value
397 max_feature_value_out[0] = max_feature_value
398
399 cdef inline void next_p(self, intp_t* p_prev, intp_t* p) noexcept nogil:
400 """Compute the next p_prev and p for iterating over feature values."""
401 cdef:
402 intp_t p_next
403 float32_t[::1] feature_values = self.feature_values
404
405 if p[0] + 1 != self.end_negative:
406 p_next = p[0] + 1
407 else:
408 p_next = self.start_positive
409
410 while (p_next < self.end and
411 feature_values[p_next] <= feature_values[p[0]] + FEATURE_THRESHOLD):
412 p[0] = p_next
413 if p[0] + 1 != self.end_negative:
414 p_next = p[0] + 1
415 else:
416 p_next = self.start_positive
417
418 p_prev[0] = p[0]
419 p[0] = p_next
420
421 cdef inline intp_t partition_samples(
422 self,
423 float64_t current_threshold
424 ) noexcept nogil:
425 """Partition samples for feature_values at the current_threshold."""
426 return self._partition(current_threshold, self.start_positive)
427
428 cdef inline void partition_samples_final(
429 self,
430 intp_t best_pos,
431 float64_t best_threshold,
432 intp_t best_feature,
433 intp_t n_missing,
434 ) noexcept nogil:
435 """Partition samples for X at the best_threshold and best_feature."""
436 self.extract_nnz(best_feature)
437 self._partition(best_threshold, best_pos)
438
439 cdef inline intp_t _partition(self, float64_t threshold, intp_t zero_pos) noexcept nogil:
440 """Partition samples[start:end] based on threshold."""
441 cdef:
442 intp_t p, partition_end
443 intp_t[::1] index_to_samples = self.index_to_samples
444 float32_t[::1] feature_values = self.feature_values
445 intp_t[::1] samples = self.samples
446
447 if threshold < 0.:
448 p = self.start
449 partition_end = self.end_negative
450 elif threshold > 0.:
451 p = self.start_positive
452 partition_end = self.end
453 else:
454 # Data are already split
455 return zero_pos
456
457 while p < partition_end:
458 if feature_values[p] <= threshold:
459 p += 1
460
461 else:
462 partition_end -= 1
463
464 feature_values[p], feature_values[partition_end] = (
465 feature_values[partition_end], feature_values[p]
466 )
467 sparse_swap(index_to_samples, samples, p, partition_end)
468
469 return partition_end
470
471 cdef inline void extract_nnz(self, intp_t feature) noexcept nogil:
472 """Extract and partition values for a given feature.
473
474 The extracted values are partitioned between negative values
475 feature_values[start:end_negative[0]] and positive values
476 feature_values[start_positive[0]:end].
477 The samples and index_to_samples are modified according to this
478 partition.
479
480 The extraction corresponds to the intersection between the arrays
481 X_indices[indptr_start:indptr_end] and samples[start:end].
482 This is done efficiently using either an index_to_samples based approach
483 or binary search based approach.
484
485 Parameters
486 ----------
487 feature : intp_t,
488 Index of the feature we want to extract non zero value.
489 """
490 cdef intp_t[::1] samples = self.samples
491 cdef float32_t[::1] feature_values = self.feature_values
492 cdef intp_t indptr_start = self.X_indptr[feature],
493 cdef intp_t indptr_end = self.X_indptr[feature + 1]
494 cdef intp_t n_indices = <intp_t>(indptr_end - indptr_start)
495 cdef intp_t n_samples = self.end - self.start
496 cdef intp_t[::1] index_to_samples = self.index_to_samples
497 cdef intp_t[::1] sorted_samples = self.sorted_samples
498 cdef const int32_t[::1] X_indices = self.X_indices
499 cdef const float32_t[::1] X_data = self.X_data
500
501 # Use binary search if n_samples * log(n_indices) <
502 # n_indices and index_to_samples approach otherwise.
503 # O(n_samples * log(n_indices)) is the running time of binary
504 # search and O(n_indices) is the running time of index_to_samples
505 # approach.
506 if ((1 - self.is_samples_sorted) * n_samples * log2(n_samples) +
507 n_samples * log2(n_indices) < EXTRACT_NNZ_SWITCH * n_indices):
508 extract_nnz_binary_search(X_indices, X_data,
509 indptr_start, indptr_end,
510 samples, self.start, self.end,
511 index_to_samples,
512 feature_values,
513 &self.end_negative, &self.start_positive,
514 sorted_samples, &self.is_samples_sorted)
515
516 # Using an index to samples technique to extract non zero values
517 # index_to_samples is a mapping from X_indices to samples
518 else:
519 extract_nnz_index_to_samples(X_indices, X_data,
520 indptr_start, indptr_end,
521 samples, self.start, self.end,
522 index_to_samples,
523 feature_values,
524 &self.end_negative, &self.start_positive)
525
526
527cdef int compare_SIZE_t(const void* a, const void* b) noexcept nogil:
528 """Comparison function for sort.
529
530 This must return an `int` as it is used by stdlib's qsort, which expects
531 an `int` return value.
532 """
533 return <int>((<intp_t*>a)[0] - (<intp_t*>b)[0])
534
535
536cdef inline void binary_search(const int32_t[::1] sorted_array,
537 int32_t start, int32_t end,
538 intp_t value, intp_t* index,
539 int32_t* new_start) noexcept nogil:
540 """Return the index of value in the sorted array.
541
542 If not found, return -1. new_start is the last pivot + 1
543 """
544 cdef int32_t pivot
545 index[0] = -1
546 while start < end:
547 pivot = start + (end - start) / 2
548
549 if sorted_array[pivot] == value:
550 index[0] = pivot
551 start = pivot + 1
552 break
553
554 if sorted_array[pivot] < value:
555 start = pivot + 1
556 else:
557 end = pivot
558 new_start[0] = start
559
560
561cdef inline void extract_nnz_index_to_samples(const int32_t[::1] X_indices,
562 const float32_t[::1] X_data,
563 int32_t indptr_start,
564 int32_t indptr_end,
565 intp_t[::1] samples,
566 intp_t start,
567 intp_t end,
568 intp_t[::1] index_to_samples,
569 float32_t[::1] feature_values,
570 intp_t* end_negative,
571 intp_t* start_positive) noexcept nogil:
572 """Extract and partition values for a feature using index_to_samples.
573
574 Complexity is O(indptr_end - indptr_start).
575 """
576 cdef int32_t k
577 cdef intp_t index
578 cdef intp_t end_negative_ = start
579 cdef intp_t start_positive_ = end
580
581 for k in range(indptr_start, indptr_end):
582 if start <= index_to_samples[X_indices[k]] < end:
583 if X_data[k] > 0:
584 start_positive_ -= 1
585 feature_values[start_positive_] = X_data[k]
586 index = index_to_samples[X_indices[k]]
587 sparse_swap(index_to_samples, samples, index, start_positive_)
588
589 elif X_data[k] < 0:
590 feature_values[end_negative_] = X_data[k]
591 index = index_to_samples[X_indices[k]]
592 sparse_swap(index_to_samples, samples, index, end_negative_)
593 end_negative_ += 1
594
595 # Returned values
596 end_negative[0] = end_negative_
597 start_positive[0] = start_positive_
598
599
600cdef inline void extract_nnz_binary_search(const int32_t[::1] X_indices,
601 const float32_t[::1] X_data,
602 int32_t indptr_start,
603 int32_t indptr_end,
604 intp_t[::1] samples,
605 intp_t start,
606 intp_t end,
607 intp_t[::1] index_to_samples,
608 float32_t[::1] feature_values,
609 intp_t* end_negative,
610 intp_t* start_positive,
611 intp_t[::1] sorted_samples,
612 bint* is_samples_sorted) noexcept nogil:
613 """Extract and partition values for a given feature using binary search.
614
615 If n_samples = end - start and n_indices = indptr_end - indptr_start,
616 the complexity is
617
618 O((1 - is_samples_sorted[0]) * n_samples * log(n_samples) +
619 n_samples * log(n_indices)).
620 """
621 cdef intp_t n_samples
622
623 if not is_samples_sorted[0]:
624 n_samples = end - start
625 memcpy(&sorted_samples[start], &samples[start],
626 n_samples * sizeof(intp_t))
627 qsort(&sorted_samples[start], n_samples, sizeof(intp_t),
628 compare_SIZE_t)
629 is_samples_sorted[0] = 1
630
631 while (indptr_start < indptr_end and
632 sorted_samples[start] > X_indices[indptr_start]):
633 indptr_start += 1
634
635 while (indptr_start < indptr_end and
636 sorted_samples[end - 1] < X_indices[indptr_end - 1]):
637 indptr_end -= 1
638
639 cdef intp_t p = start
640 cdef intp_t index
641 cdef intp_t k
642 cdef intp_t end_negative_ = start
643 cdef intp_t start_positive_ = end
644
645 while (p < end and indptr_start < indptr_end):
646 # Find index of sorted_samples[p] in X_indices
647 binary_search(X_indices, indptr_start, indptr_end,
648 sorted_samples[p], &k, &indptr_start)
649
650 if k != -1:
651 # If k != -1, we have found a non zero value
652
653 if X_data[k] > 0:
654 start_positive_ -= 1
655 feature_values[start_positive_] = X_data[k]
656 index = index_to_samples[X_indices[k]]
657 sparse_swap(index_to_samples, samples, index, start_positive_)
658
659 elif X_data[k] < 0:
660 feature_values[end_negative_] = X_data[k]
661 index = index_to_samples[X_indices[k]]
662 sparse_swap(index_to_samples, samples, index, end_negative_)
663 end_negative_ += 1
664 p += 1
665
666 # Returned values
667 end_negative[0] = end_negative_
668 start_positive[0] = start_positive_
669
670
671cdef inline void sparse_swap(intp_t[::1] index_to_samples, intp_t[::1] samples,
672 intp_t pos_1, intp_t pos_2) noexcept nogil:
673 """Swap sample pos_1 and pos_2 preserving sparse invariant."""
674 samples[pos_1], samples[pos_2] = samples[pos_2], samples[pos_1]
675 index_to_samples[samples[pos_1]] = pos_1
676 index_to_samples[samples[pos_2]] = pos_2
677
678
679cdef inline void shift_missing_values_to_left_if_required(
680 SplitRecord* best,
681 intp_t[::1] samples,
682 intp_t end,
683) noexcept nogil:
684 """Shift missing value sample indices to the left of the split if required.
685
686 Note: this should always be called at the very end because it will
687 move samples around, thereby affecting the criterion.
688 This affects the computation of the children impurity, which affects
689 the computation of the next node.
690 """
691 cdef intp_t i, p, current_end
692 # The partitioner partitions the data such that the missing values are in
693 # samples[-n_missing:] for the criterion to consume. If the missing values
694 # are going to the right node, then the missing values are already in the
695 # correct position. If the missing values go left, then we move the missing
696 # values to samples[best.pos:best.pos+n_missing] and update `best.pos`.
697 if best.n_missing > 0 and best.missing_go_to_left:
698 for p in range(best.n_missing):
699 i = best.pos + p
700 current_end = end - 1 - p
701 samples[i], samples[current_end] = samples[current_end], samples[i]
702 best.pos += best.n_missing
703
704
705def _py_sort(float32_t[::1] feature_values, intp_t[::1] samples, intp_t n):
706 """Used for testing sort."""
707 sort(&feature_values[0], &samples[0], n)
708
709
710# Sort n-element arrays pointed to by feature_values and samples, simultaneously,
711# by the values in feature_values. Algorithm: Introsort (Musser, SP&E, 1997).
712cdef inline void sort(float32_t* feature_values, intp_t* samples, intp_t n) noexcept nogil:
713 if n == 0:
714 return
715 cdef intp_t maxd = 2 * <intp_t>log2(n)
716 introsort(feature_values, samples, n, maxd)
717
718
719cdef inline void swap(float32_t* feature_values, intp_t* samples,
720 intp_t i, intp_t j) noexcept nogil:
721 # Helper for sort
722 feature_values[i], feature_values[j] = feature_values[j], feature_values[i]
723 samples[i], samples[j] = samples[j], samples[i]
724
725
726cdef inline float32_t median3(float32_t* feature_values, intp_t n) noexcept nogil:
727 # Median of three pivot selection, after Bentley and McIlroy (1993).
728 # Engineering a sort function. SP&E. Requires 8/3 comparisons on average.
729 cdef float32_t a = feature_values[0], b = feature_values[n / 2], c = feature_values[n - 1]
730 if a < b:
731 if b < c:
732 return b
733 elif a < c:
734 return c
735 else:
736 return a
737 elif b < c:
738 if a < c:
739 return a
740 else:
741 return c
742 else:
743 return b
744
745
746# Introsort with median of 3 pivot selection and 3-way partition function
747# (robust to repeated elements, e.g. lots of zero features).
748cdef void introsort(float32_t* feature_values, intp_t *samples,
749 intp_t n, intp_t maxd) noexcept nogil:
750 cdef float32_t pivot
751 cdef intp_t i, l, r
752
753 while n > 1:
754 if maxd <= 0: # max depth limit exceeded ("gone quadratic")
755 heapsort(feature_values, samples, n)
756 return
757 maxd -= 1
758
759 pivot = median3(feature_values, n)
760
761 # Three-way partition.
762 i = l = 0
763 r = n
764 while i < r:
765 if feature_values[i] < pivot:
766 swap(feature_values, samples, i, l)
767 i += 1
768 l += 1
769 elif feature_values[i] > pivot:
770 r -= 1
771 swap(feature_values, samples, i, r)
772 else:
773 i += 1
774
775 introsort(feature_values, samples, l, maxd)
776 feature_values += r
777 samples += r
778 n -= r
779
780
781cdef inline void sift_down(float32_t* feature_values, intp_t* samples,
782 intp_t start, intp_t end) noexcept nogil:
783 # Restore heap order in feature_values[start:end] by moving the max element to start.
784 cdef intp_t child, maxind, root
785
786 root = start
787 while True:
788 child = root * 2 + 1
789
790 # find max of root, left child, right child
791 maxind = root
792 if child < end and feature_values[maxind] < feature_values[child]:
793 maxind = child
794 if child + 1 < end and feature_values[maxind] < feature_values[child + 1]:
795 maxind = child + 1
796
797 if maxind == root:
798 break
799 else:
800 swap(feature_values, samples, root, maxind)
801 root = maxind
802
803
804cdef void heapsort(float32_t* feature_values, intp_t* samples, intp_t n) noexcept nogil:
805 cdef intp_t start, end
806
807 # heapify
808 start = (n - 2) / 2
809 end = n
810 while True:
811 sift_down(feature_values, samples, start, end)
812 if start == 0:
813 break
814 start -= 1
815
816 # sort by shrinking the heap, putting the max element immediately after it
817 end = n - 1
818 while end > 0:
819 swap(feature_values, samples, 0, end)
820 sift_down(feature_values, samples, 0, end)
821 end = end - 1
822