Aluode/PerceptionLabPortable
0
1# Authors: The scikit-learn developers
2# SPDX-License-Identifier: BSD-3-Clause
3
4"""
5Random utility function
6=======================
7This module complements missing features of ``numpy.random``.
8
9The module contains:
10 * Several algorithms to sample integers without replacement.
11 * Fast rand_r alternative based on xor shifts
12"""
13import numpy as np
14from . import check_random_state
15
16from ._typedefs cimport intp_t
17
18
19cdef uint32_t DEFAULT_SEED = 1
20
21
22# Compatibility type to always accept the default int type used by NumPy, both
23# before and after NumPy 2. On Windows, `long` does not always match `inp_t`.
24# See the comments in the `sample_without_replacement` Python function for more
25# details.
26ctypedef fused default_int:
27 intp_t
28 long
29
30
31cpdef _sample_without_replacement_check_input(default_int n_population,
32 default_int n_samples):
33 """ Check that input are consistent for sample_without_replacement"""
34 if n_population < 0:
35 raise ValueError('n_population should be greater than 0, got %s.'
36 % n_population)
37
38 if n_samples > n_population:
39 raise ValueError('n_population should be greater or equal than '
40 'n_samples, got n_samples > n_population (%s > %s)'
41 % (n_samples, n_population))
42
43
44cpdef _sample_without_replacement_with_tracking_selection(
45 default_int n_population,
46 default_int n_samples,
47 random_state=None):
48 r"""Sample integers without replacement.
49
50 Select n_samples integers from the set [0, n_population) without
51 replacement.
52
53 Time complexity:
54 - Worst-case: unbounded
55 - Average-case:
56 O(O(np.random.randint) * \sum_{i=1}^n_samples 1 /
57 (1 - i / n_population)))
58 <= O(O(np.random.randint) *
59 n_population * ln((n_population - 2)
60 /(n_population - 1 - n_samples)))
61 <= O(O(np.random.randint) *
62 n_population * 1 / (1 - n_samples / n_population))
63
64 Space complexity of O(n_samples) in a python set.
65
66
67 Parameters
68 ----------
69 n_population : int
70 The size of the set to sample from.
71
72 n_samples : int
73 The number of integer to sample.
74
75 random_state : int, RandomState instance or None, default=None
76 If int, random_state is the seed used by the random number generator;
77 If RandomState instance, random_state is the random number generator;
78 If None, the random number generator is the RandomState instance used
79 by `np.random`.
80
81 Returns
82 -------
83 out : ndarray of shape (n_samples,)
84 The sampled subsets of integer.
85 """
86 _sample_without_replacement_check_input(n_population, n_samples)
87
88 cdef default_int i
89 cdef default_int j
90 cdef default_int[::1] out = np.empty((n_samples, ), dtype=int)
91
92 rng = check_random_state(random_state)
93 rng_randint = rng.randint
94
95 # The following line of code are heavily inspired from python core,
96 # more precisely of random.sample.
97 cdef set selected = set()
98
99 for i in range(n_samples):
100 j = rng_randint(n_population)
101 while j in selected:
102 j = rng_randint(n_population)
103 selected.add(j)
104 out[i] = j
105
106 return np.asarray(out)
107
108
109cpdef _sample_without_replacement_with_pool(default_int n_population,
110 default_int n_samples,
111 random_state=None):
112 """Sample integers without replacement.
113
114 Select n_samples integers from the set [0, n_population) without
115 replacement.
116
117 Time complexity: O(n_population + O(np.random.randint) * n_samples)
118
119 Space complexity of O(n_population + n_samples).
120
121
122 Parameters
123 ----------
124 n_population : int
125 The size of the set to sample from.
126
127 n_samples : int
128 The number of integer to sample.
129
130 random_state : int, RandomState instance or None, default=None
131 If int, random_state is the seed used by the random number generator;
132 If RandomState instance, random_state is the random number generator;
133 If None, the random number generator is the RandomState instance used
134 by `np.random`.
135
136 Returns
137 -------
138 out : ndarray of shape (n_samples,)
139 The sampled subsets of integer.
140 """
141 _sample_without_replacement_check_input(n_population, n_samples)
142
143 cdef default_int i
144 cdef default_int j
145 cdef default_int[::1] out = np.empty((n_samples,), dtype=int)
146 cdef default_int[::1] pool = np.empty((n_population,), dtype=int)
147
148 rng = check_random_state(random_state)
149 rng_randint = rng.randint
150
151 # Initialize the pool
152 for i in range(n_population):
153 pool[i] = i
154
155 # The following line of code are heavily inspired from python core,
156 # more precisely of random.sample.
157 for i in range(n_samples):
158 j = rng_randint(n_population - i) # invariant: non-selected at [0,n-i)
159 out[i] = pool[j]
160 pool[j] = pool[n_population - i - 1] # move non-selected item into vacancy
161
162 return np.asarray(out)
163
164
165cpdef _sample_without_replacement_with_reservoir_sampling(
166 default_int n_population,
167 default_int n_samples,
168 random_state=None
169):
170 """Sample integers without replacement.
171
172 Select n_samples integers from the set [0, n_population) without
173 replacement.
174
175 Time complexity of
176 O((n_population - n_samples) * O(np.random.randint) + n_samples)
177 Space complexity of O(n_samples)
178
179
180 Parameters
181 ----------
182 n_population : int
183 The size of the set to sample from.
184
185 n_samples : int
186 The number of integer to sample.
187
188 random_state : int, RandomState instance or None, default=None
189 If int, random_state is the seed used by the random number generator;
190 If RandomState instance, random_state is the random number generator;
191 If None, the random number generator is the RandomState instance used
192 by `np.random`.
193
194 Returns
195 -------
196 out : ndarray of shape (n_samples,)
197 The sampled subsets of integer. The order of the items is not
198 necessarily random. Use a random permutation of the array if the order
199 of the items has to be randomized.
200 """
201 _sample_without_replacement_check_input(n_population, n_samples)
202
203 cdef default_int i
204 cdef default_int j
205 cdef default_int[::1] out = np.empty((n_samples, ), dtype=int)
206
207 rng = check_random_state(random_state)
208 rng_randint = rng.randint
209
210 # This cython implementation is based on the one of Robert Kern:
211 # http://mail.scipy.org/pipermail/numpy-discussion/2010-December/
212 # 054289.html
213 #
214 for i in range(n_samples):
215 out[i] = i
216
217 for i from n_samples <= i < n_population:
218 j = rng_randint(0, i + 1)
219 if j < n_samples:
220 out[j] = i
221
222 return np.asarray(out)
223
224
225cdef _sample_without_replacement(default_int n_population,
226 default_int n_samples,
227 method="auto",
228 random_state=None):
229 """Sample integers without replacement.
230
231 Private function for the implementation, see sample_without_replacement
232 documentation for more details.
233 """
234 _sample_without_replacement_check_input(n_population, n_samples)
235
236 all_methods = ("auto", "tracking_selection", "reservoir_sampling", "pool")
237
238 ratio = <double> n_samples / n_population if n_population != 0.0 else 1.0
239
240 # Check ratio and use permutation unless ratio < 0.01 or ratio > 0.99
241 if method == "auto" and ratio > 0.01 and ratio < 0.99:
242 rng = check_random_state(random_state)
243 return rng.permutation(n_population)[:n_samples]
244
245 if method == "auto" or method == "tracking_selection":
246 # TODO the pool based method can also be used.
247 # however, it requires special benchmark to take into account
248 # the memory requirement of the array vs the set.
249
250 # The value 0.2 has been determined through benchmarking.
251 if ratio < 0.2:
252 return _sample_without_replacement_with_tracking_selection(
253 n_population, n_samples, random_state)
254 else:
255 return _sample_without_replacement_with_reservoir_sampling(
256 n_population, n_samples, random_state)
257
258 elif method == "reservoir_sampling":
259 return _sample_without_replacement_with_reservoir_sampling(
260 n_population, n_samples, random_state)
261
262 elif method == "pool":
263 return _sample_without_replacement_with_pool(n_population, n_samples,
264 random_state)
265 else:
266 raise ValueError('Expected a method name in %s, got %s. '
267 % (all_methods, method))
268
269
270def sample_without_replacement(
271 object n_population, object n_samples, method="auto", random_state=None):
272 """Sample integers without replacement.
273
274 Select n_samples integers from the set [0, n_population) without
275 replacement.
276
277
278 Parameters
279 ----------
280 n_population : int
281 The size of the set to sample from.
282
283 n_samples : int
284 The number of integer to sample.
285
286 random_state : int, RandomState instance or None, default=None
287 If int, random_state is the seed used by the random number generator;
288 If RandomState instance, random_state is the random number generator;
289 If None, the random number generator is the RandomState instance used
290 by `np.random`.
291
292 method : {"auto", "tracking_selection", "reservoir_sampling", "pool"}, \
293 default='auto'
294 If method == "auto", the ratio of n_samples / n_population is used
295 to determine which algorithm to use:
296 If ratio is between 0 and 0.01, tracking selection is used.
297 If ratio is between 0.01 and 0.99, numpy.random.permutation is used.
298 If ratio is greater than 0.99, reservoir sampling is used.
299 The order of the selected integers is undefined. If a random order is
300 desired, the selected subset should be shuffled.
301
302 If method =="tracking_selection", a set based implementation is used
303 which is suitable for `n_samples` <<< `n_population`.
304
305 If method == "reservoir_sampling", a reservoir sampling algorithm is
306 used which is suitable for high memory constraint or when
307 O(`n_samples`) ~ O(`n_population`).
308 The order of the selected integers is undefined. If a random order is
309 desired, the selected subset should be shuffled.
310
311 If method == "pool", a pool based algorithm is particularly fast, even
312 faster than the tracking selection method. However, a vector containing
313 the entire population has to be initialized.
314 If n_samples ~ n_population, the reservoir sampling method is faster.
315
316 Returns
317 -------
318 out : ndarray of shape (n_samples,)
319 The sampled subsets of integer. The subset of selected integer might
320 not be randomized, see the method argument.
321
322 Examples
323 --------
324 >>> from sklearn.utils.random import sample_without_replacement
325 >>> sample_without_replacement(10, 5, random_state=42)
326 array([8, 1, 5, 0, 7])
327 """
328 cdef:
329 intp_t n_pop_intp, n_samples_intp
330 long n_pop_long, n_samples_long
331
332 # On most platforms `np.int_ is np.intp`. However, before NumPy 2 the
333 # default integer `np.int_` was a long which is 32bit on 64bit windows
334 # while `intp` is 64bit on 64bit platforms and 32bit on 32bit ones.
335 if np.int_ is np.intp:
336 # Branch always taken on NumPy >=2 (or when not on 64bit windows).
337 # Cython has different rules for conversion of values to integers.
338 # For NumPy <1.26.2 AND Cython 3, this first branch requires `int()`
339 # called explicitly to allow e.g. floats.
340 n_pop_intp = int(n_population)
341 n_samples_intp = int(n_samples)
342 return _sample_without_replacement(
343 n_pop_intp, n_samples_intp, method, random_state)
344 else:
345 # Branch taken on 64bit windows with Numpy<2.0 where `long` is 32bit
346 n_pop_long = n_population
347 n_samples_long = n_samples
348 return _sample_without_replacement(
349 n_pop_long, n_samples_long, method, random_state)
350
351
352def _our_rand_r_py(seed):
353 """Python utils to test the our_rand_r function"""
354 cdef uint32_t my_seed = seed
355 return our_rand_r(&my_seed)
356