Aluode/PerceptionLabPortable
0
1"""
2This file contains preprocessing tools based on polynomials.
3"""
4
5# Authors: The scikit-learn developers
6# SPDX-License-Identifier: BSD-3-Clause
7
8import collections
9from itertools import chain, combinations
10from itertools import combinations_with_replacement as combinations_w_r
11from numbers import Integral
12
13import numpy as np
14from scipy import sparse
15from scipy.interpolate import BSpline
16from scipy.special import comb
17
18from ..base import BaseEstimator, TransformerMixin, _fit_context
19from ..utils import check_array
20from ..utils._param_validation import Interval, StrOptions
21from ..utils.fixes import parse_version, sp_version
22from ..utils.stats import _weighted_percentile
23from ..utils.validation import (
24 FLOAT_DTYPES,
25 _check_feature_names_in,
26 _check_sample_weight,
27 check_is_fitted,
28 validate_data,
29)
30from ._csr_polynomial_expansion import (
31 _calc_expanded_nnz,
32 _calc_total_nnz,
33 _csr_polynomial_expansion,
34)
35
36__all__ = [
37 "PolynomialFeatures",
38 "SplineTransformer",
39]
40
41
42def _create_expansion(X, interaction_only, deg, n_features, cumulative_size=0):
43 """Helper function for creating and appending sparse expansion matrices"""
44
45 total_nnz = _calc_total_nnz(X.indptr, interaction_only, deg)
46 expanded_col = _calc_expanded_nnz(n_features, interaction_only, deg)
47
48 if expanded_col == 0:
49 return None
50 # This only checks whether each block needs 64bit integers upon
51 # expansion. We prefer to keep int32 indexing where we can,
52 # since currently SciPy's CSR construction downcasts when possible,
53 # so we prefer to avoid an unnecessary cast. The dtype may still
54 # change in the concatenation process if needed.
55 # See: https://github.com/scipy/scipy/issues/16569
56 max_indices = expanded_col - 1
57 max_indptr = total_nnz
58 max_int32 = np.iinfo(np.int32).max
59 needs_int64 = max(max_indices, max_indptr) > max_int32
60 index_dtype = np.int64 if needs_int64 else np.int32
61
62 # Result of the expansion, modified in place by the
63 # `_csr_polynomial_expansion` routine.
64 expanded_data = np.empty(shape=total_nnz, dtype=X.data.dtype)
65 expanded_indices = np.empty(shape=total_nnz, dtype=index_dtype)
66 expanded_indptr = np.empty(shape=X.indptr.shape[0], dtype=index_dtype)
67 _csr_polynomial_expansion(
68 X.data,
69 X.indices,
70 X.indptr,
71 X.shape[1],
72 expanded_data,
73 expanded_indices,
74 expanded_indptr,
75 interaction_only,
76 deg,
77 )
78 return sparse.csr_matrix(
79 (expanded_data, expanded_indices, expanded_indptr),
80 shape=(X.indptr.shape[0] - 1, expanded_col),
81 dtype=X.dtype,
82 )
83
84
85class PolynomialFeatures(TransformerMixin, BaseEstimator):
86 """Generate polynomial and interaction features.
87
88 Generate a new feature matrix consisting of all polynomial combinations
89 of the features with degree less than or equal to the specified degree.
90 For example, if an input sample is two dimensional and of the form
91 [a, b], the degree-2 polynomial features are [1, a, b, a^2, ab, b^2].
92
93 Read more in the :ref:`User Guide <polynomial_features>`.
94
95 Parameters
96 ----------
97 degree : int or tuple (min_degree, max_degree), default=2
98 If a single int is given, it specifies the maximal degree of the
99 polynomial features. If a tuple `(min_degree, max_degree)` is passed,
100 then `min_degree` is the minimum and `max_degree` is the maximum
101 polynomial degree of the generated features. Note that `min_degree=0`
102 and `min_degree=1` are equivalent as outputting the degree zero term is
103 determined by `include_bias`.
104
105 interaction_only : bool, default=False
106 If `True`, only interaction features are produced: features that are
107 products of at most `degree` *distinct* input features, i.e. terms with
108 power of 2 or higher of the same input feature are excluded:
109
110 - included: `x[0]`, `x[1]`, `x[0] * x[1]`, etc.
111 - excluded: `x[0] ** 2`, `x[0] ** 2 * x[1]`, etc.
112
113 include_bias : bool, default=True
114 If `True` (default), then include a bias column, the feature in which
115 all polynomial powers are zero (i.e. a column of ones - acts as an
116 intercept term in a linear model).
117
118 order : {'C', 'F'}, default='C'
119 Order of output array in the dense case. `'F'` order is faster to
120 compute, but may slow down subsequent estimators.
121
122 .. versionadded:: 0.21
123
124 Attributes
125 ----------
126 powers_ : ndarray of shape (`n_output_features_`, `n_features_in_`)
127 `powers_[i, j]` is the exponent of the jth input in the ith output.
128
129 n_features_in_ : int
130 Number of features seen during :term:`fit`.
131
132 .. versionadded:: 0.24
133
134 feature_names_in_ : ndarray of shape (`n_features_in_`,)
135 Names of features seen during :term:`fit`. Defined only when `X`
136 has feature names that are all strings.
137
138 .. versionadded:: 1.0
139
140 n_output_features_ : int
141 The total number of polynomial output features. The number of output
142 features is computed by iterating over all suitably sized combinations
143 of input features.
144
145 See Also
146 --------
147 SplineTransformer : Transformer that generates univariate B-spline bases
148 for features.
149
150 Notes
151 -----
152 Be aware that the number of features in the output array scales
153 polynomially in the number of features of the input array, and
154 exponentially in the degree. High degrees can cause overfitting.
155
156 See :ref:`examples/linear_model/plot_polynomial_interpolation.py
157 <sphx_glr_auto_examples_linear_model_plot_polynomial_interpolation.py>`
158
159 Examples
160 --------
161 >>> import numpy as np
162 >>> from sklearn.preprocessing import PolynomialFeatures
163 >>> X = np.arange(6).reshape(3, 2)
164 >>> X
165 array([[0, 1],
166 [2, 3],
167 [4, 5]])
168 >>> poly = PolynomialFeatures(2)
169 >>> poly.fit_transform(X)
170 array([[ 1., 0., 1., 0., 0., 1.],
171 [ 1., 2., 3., 4., 6., 9.],
172 [ 1., 4., 5., 16., 20., 25.]])
173 >>> poly = PolynomialFeatures(interaction_only=True)
174 >>> poly.fit_transform(X)
175 array([[ 1., 0., 1., 0.],
176 [ 1., 2., 3., 6.],
177 [ 1., 4., 5., 20.]])
178 """
179
180 _parameter_constraints: dict = {
181 "degree": [Interval(Integral, 0, None, closed="left"), "array-like"],
182 "interaction_only": ["boolean"],
183 "include_bias": ["boolean"],
184 "order": [StrOptions({"C", "F"})],
185 }
186
187 def __init__(
188 self, degree=2, *, interaction_only=False, include_bias=True, order="C"
189 ):
190 self.degree = degree
191 self.interaction_only = interaction_only
192 self.include_bias = include_bias
193 self.order = order
194
195 @staticmethod
196 def _combinations(
197 n_features, min_degree, max_degree, interaction_only, include_bias
198 ):
199 comb = combinations if interaction_only else combinations_w_r
200 start = max(1, min_degree)
201 iter = chain.from_iterable(
202 comb(range(n_features), i) for i in range(start, max_degree + 1)
203 )
204 if include_bias:
205 iter = chain(comb(range(n_features), 0), iter)
206 return iter
207
208 @staticmethod
209 def _num_combinations(
210 n_features, min_degree, max_degree, interaction_only, include_bias
211 ):
212 """Calculate number of terms in polynomial expansion
213
214 This should be equivalent to counting the number of terms returned by
215 _combinations(...) but much faster.
216 """
217
218 if interaction_only:
219 combinations = sum(
220 [
221 comb(n_features, i, exact=True)
222 for i in range(max(1, min_degree), min(max_degree, n_features) + 1)
223 ]
224 )
225 else:
226 combinations = comb(n_features + max_degree, max_degree, exact=True) - 1
227 if min_degree > 0:
228 d = min_degree - 1
229 combinations -= comb(n_features + d, d, exact=True) - 1
230
231 if include_bias:
232 combinations += 1
233
234 return combinations
235
236 @property
237 def powers_(self):
238 """Exponent for each of the inputs in the output."""
239 check_is_fitted(self)
240
241 combinations = self._combinations(
242 n_features=self.n_features_in_,
243 min_degree=self._min_degree,
244 max_degree=self._max_degree,
245 interaction_only=self.interaction_only,
246 include_bias=self.include_bias,
247 )
248 return np.vstack(
249 [np.bincount(c, minlength=self.n_features_in_) for c in combinations]
250 )
251
252 def get_feature_names_out(self, input_features=None):
253 """Get output feature names for transformation.
254
255 Parameters
256 ----------
257 input_features : array-like of str or None, default=None
258 Input features.
259
260 - If `input_features is None`, then `feature_names_in_` is
261 used as feature names in. If `feature_names_in_` is not defined,
262 then the following input feature names are generated:
263 `["x0", "x1", ..., "x(n_features_in_ - 1)"]`.
264 - If `input_features` is an array-like, then `input_features` must
265 match `feature_names_in_` if `feature_names_in_` is defined.
266
267 Returns
268 -------
269 feature_names_out : ndarray of str objects
270 Transformed feature names.
271 """
272 powers = self.powers_
273 input_features = _check_feature_names_in(self, input_features)
274 feature_names = []
275 for row in powers:
276 inds = np.where(row)[0]
277 if len(inds):
278 name = " ".join(
279 (
280 "%s^%d" % (input_features[ind], exp)
281 if exp != 1
282 else input_features[ind]
283 )
284 for ind, exp in zip(inds, row[inds])
285 )
286 else:
287 name = "1"
288 feature_names.append(name)
289 return np.asarray(feature_names, dtype=object)
290
291 @_fit_context(prefer_skip_nested_validation=True)
292 def fit(self, X, y=None):
293 """
294 Compute number of output features.
295
296 Parameters
297 ----------
298 X : {array-like, sparse matrix} of shape (n_samples, n_features)
299 The data.
300
301 y : Ignored
302 Not used, present here for API consistency by convention.
303
304 Returns
305 -------
306 self : object
307 Fitted transformer.
308 """
309 _, n_features = validate_data(self, X, accept_sparse=True).shape
310
311 if isinstance(self.degree, Integral):
312 if self.degree == 0 and not self.include_bias:
313 raise ValueError(
314 "Setting degree to zero and include_bias to False would result in"
315 " an empty output array."
316 )
317
318 self._min_degree = 0
319 self._max_degree = self.degree
320 elif (
321 isinstance(self.degree, collections.abc.Iterable) and len(self.degree) == 2
322 ):
323 self._min_degree, self._max_degree = self.degree
324 if not (
325 isinstance(self._min_degree, Integral)
326 and isinstance(self._max_degree, Integral)
327 and self._min_degree >= 0
328 and self._min_degree <= self._max_degree
329 ):
330 raise ValueError(
331 "degree=(min_degree, max_degree) must "
332 "be non-negative integers that fulfil "
333 "min_degree <= max_degree, got "
334 f"{self.degree}."
335 )
336 elif self._max_degree == 0 and not self.include_bias:
337 raise ValueError(
338 "Setting both min_degree and max_degree to zero and include_bias to"
339 " False would result in an empty output array."
340 )
341 else:
342 raise ValueError(
343 "degree must be a non-negative int or tuple "
344 "(min_degree, max_degree), got "
345 f"{self.degree}."
346 )
347
348 self.n_output_features_ = self._num_combinations(
349 n_features=n_features,
350 min_degree=self._min_degree,
351 max_degree=self._max_degree,
352 interaction_only=self.interaction_only,
353 include_bias=self.include_bias,
354 )
355 if self.n_output_features_ > np.iinfo(np.intp).max:
356 msg = (
357 "The output that would result from the current configuration would"
358 f" have {self.n_output_features_} features which is too large to be"
359 f" indexed by {np.intp().dtype.name}. Please change some or all of the"
360 " following:\n- The number of features in the input, currently"
361 f" {n_features=}\n- The range of degrees to calculate, currently"
362 f" [{self._min_degree}, {self._max_degree}]\n- Whether to include only"
363 f" interaction terms, currently {self.interaction_only}\n- Whether to"
364 f" include a bias term, currently {self.include_bias}."
365 )
366 if (
367 np.intp == np.int32
368 and self.n_output_features_ <= np.iinfo(np.int64).max
369 ): # pragma: nocover
370 msg += (
371 "\nNote that the current Python runtime has a limited 32 bit "
372 "address space and that this configuration would have been "
373 "admissible if run on a 64 bit Python runtime."
374 )
375 raise ValueError(msg)
376 # We also record the number of output features for
377 # _min_degree = 0
378 self._n_out_full = self._num_combinations(
379 n_features=n_features,
380 min_degree=0,
381 max_degree=self._max_degree,
382 interaction_only=self.interaction_only,
383 include_bias=self.include_bias,
384 )
385
386 return self
387
388 def transform(self, X):
389 """Transform data to polynomial features.
390
391 Parameters
392 ----------
393 X : {array-like, sparse matrix} of shape (n_samples, n_features)
394 The data to transform, row by row.
395
396 Prefer CSR over CSC for sparse input (for speed), but CSC is
397 required if the degree is 4 or higher. If the degree is less than
398 4 and the input format is CSC, it will be converted to CSR, have
399 its polynomial features generated, then converted back to CSC.
400
401 If the degree is 2 or 3, the method described in "Leveraging
402 Sparsity to Speed Up Polynomial Feature Expansions of CSR Matrices
403 Using K-Simplex Numbers" by Andrew Nystrom and John Hughes is
404 used, which is much faster than the method used on CSC input. For
405 this reason, a CSC input will be converted to CSR, and the output
406 will be converted back to CSC prior to being returned, hence the
407 preference of CSR.
408
409 Returns
410 -------
411 XP : {ndarray, sparse matrix} of shape (n_samples, NP)
412 The matrix of features, where `NP` is the number of polynomial
413 features generated from the combination of inputs. If a sparse
414 matrix is provided, it will be converted into a sparse
415 `csr_matrix`.
416 """
417 check_is_fitted(self)
418
419 X = validate_data(
420 self,
421 X,
422 order="F",
423 dtype=FLOAT_DTYPES,
424 reset=False,
425 accept_sparse=("csr", "csc"),
426 )
427
428 n_samples, n_features = X.shape
429 max_int32 = np.iinfo(np.int32).max
430 if sparse.issparse(X) and X.format == "csr":
431 if self._max_degree > 3:
432 return self.transform(X.tocsc()).tocsr()
433 to_stack = []
434 if self.include_bias:
435 to_stack.append(
436 sparse.csr_matrix(np.ones(shape=(n_samples, 1), dtype=X.dtype))
437 )
438 if self._min_degree <= 1 and self._max_degree > 0:
439 to_stack.append(X)
440
441 cumulative_size = sum(mat.shape[1] for mat in to_stack)
442 for deg in range(max(2, self._min_degree), self._max_degree + 1):
443 expanded = _create_expansion(
444 X=X,
445 interaction_only=self.interaction_only,
446 deg=deg,
447 n_features=n_features,
448 cumulative_size=cumulative_size,
449 )
450 if expanded is not None:
451 to_stack.append(expanded)
452 cumulative_size += expanded.shape[1]
453 if len(to_stack) == 0:
454 # edge case: deal with empty matrix
455 XP = sparse.csr_matrix((n_samples, 0), dtype=X.dtype)
456 else:
457 # `scipy.sparse.hstack` breaks in scipy<1.9.2
458 # when `n_output_features_ > max_int32`
459 all_int32 = all(mat.indices.dtype == np.int32 for mat in to_stack)
460 if (
461 sp_version < parse_version("1.9.2")
462 and self.n_output_features_ > max_int32
463 and all_int32
464 ):
465 raise ValueError( # pragma: no cover
466 "In scipy versions `<1.9.2`, the function `scipy.sparse.hstack`"
467 " produces negative columns when:\n1. The output shape contains"
468 " `n_cols` too large to be represented by a 32bit signed"
469 " integer.\n2. All sub-matrices to be stacked have indices of"
470 " dtype `np.int32`.\nTo avoid this error, either use a version"
471 " of scipy `>=1.9.2` or alter the `PolynomialFeatures`"
472 " transformer to produce fewer than 2^31 output features"
473 )
474 XP = sparse.hstack(to_stack, dtype=X.dtype, format="csr")
475 elif sparse.issparse(X) and X.format == "csc" and self._max_degree < 4:
476 return self.transform(X.tocsr()).tocsc()
477 elif sparse.issparse(X):
478 combinations = self._combinations(
479 n_features=n_features,
480 min_degree=self._min_degree,
481 max_degree=self._max_degree,
482 interaction_only=self.interaction_only,
483 include_bias=self.include_bias,
484 )
485 columns = []
486 for combi in combinations:
487 if combi:
488 out_col = 1
489 for col_idx in combi:
490 out_col = X[:, [col_idx]].multiply(out_col)
491 columns.append(out_col)
492 else:
493 bias = sparse.csc_matrix(np.ones((X.shape[0], 1)))
494 columns.append(bias)
495 XP = sparse.hstack(columns, dtype=X.dtype).tocsc()
496 else:
497 # Do as if _min_degree = 0 and cut down array after the
498 # computation, i.e. use _n_out_full instead of n_output_features_.
499 XP = np.empty(
500 shape=(n_samples, self._n_out_full), dtype=X.dtype, order=self.order
501 )
502
503 # What follows is a faster implementation of:
504 # for i, comb in enumerate(combinations):
505 # XP[:, i] = X[:, comb].prod(1)
506 # This implementation uses two optimisations.
507 # First one is broadcasting,
508 # multiply ([X1, ..., Xn], X1) -> [X1 X1, ..., Xn X1]
509 # multiply ([X2, ..., Xn], X2) -> [X2 X2, ..., Xn X2]
510 # ...
511 # multiply ([X[:, start:end], X[:, start]) -> ...
512 # Second optimisation happens for degrees >= 3.
513 # Xi^3 is computed reusing previous computation:
514 # Xi^3 = Xi^2 * Xi.
515
516 # degree 0 term
517 if self.include_bias:
518 XP[:, 0] = 1
519 current_col = 1
520 else:
521 current_col = 0
522
523 if self._max_degree == 0:
524 return XP
525
526 # degree 1 term
527 XP[:, current_col : current_col + n_features] = X
528 index = list(range(current_col, current_col + n_features))
529 current_col += n_features
530 index.append(current_col)
531
532 # loop over degree >= 2 terms
533 for _ in range(2, self._max_degree + 1):
534 new_index = []
535 end = index[-1]
536 for feature_idx in range(n_features):
537 start = index[feature_idx]
538 new_index.append(current_col)
539 if self.interaction_only:
540 start += index[feature_idx + 1] - index[feature_idx]
541 next_col = current_col + end - start
542 if next_col <= current_col:
543 break
544 # XP[:, start:end] are terms of degree d - 1
545 # that exclude feature #feature_idx.
546 np.multiply(
547 XP[:, start:end],
548 X[:, feature_idx : feature_idx + 1],
549 out=XP[:, current_col:next_col],
550 casting="no",
551 )
552 current_col = next_col
553
554 new_index.append(current_col)
555 index = new_index
556
557 if self._min_degree > 1:
558 n_XP, n_Xout = self._n_out_full, self.n_output_features_
559 if self.include_bias:
560 Xout = np.empty(
561 shape=(n_samples, n_Xout), dtype=XP.dtype, order=self.order
562 )
563 Xout[:, 0] = 1
564 Xout[:, 1:] = XP[:, n_XP - n_Xout + 1 :]
565 else:
566 Xout = XP[:, n_XP - n_Xout :].copy()
567 XP = Xout
568 return XP
569
570 def __sklearn_tags__(self):
571 tags = super().__sklearn_tags__()
572 tags.input_tags.sparse = True
573 return tags
574
575
576class SplineTransformer(TransformerMixin, BaseEstimator):
577 """Generate univariate B-spline bases for features.
578
579 Generate a new feature matrix consisting of
580 `n_splines=n_knots + degree - 1` (`n_knots - 1` for
581 `extrapolation="periodic"`) spline basis functions
582 (B-splines) of polynomial order=`degree` for each feature.
583
584 In order to learn more about the SplineTransformer class go to:
585 :ref:`sphx_glr_auto_examples_applications_plot_cyclical_feature_engineering.py`
586
587 Read more in the :ref:`User Guide <spline_transformer>`.
588
589 .. versionadded:: 1.0
590
591 Parameters
592 ----------
593 n_knots : int, default=5
594 Number of knots of the splines if `knots` equals one of
595 {'uniform', 'quantile'}. Must be larger or equal 2. Ignored if `knots`
596 is array-like.
597
598 degree : int, default=3
599 The polynomial degree of the spline basis. Must be a non-negative
600 integer.
601
602 knots : {'uniform', 'quantile'} or array-like of shape \
603 (n_knots, n_features), default='uniform'
604 Set knot positions such that first knot <= features <= last knot.
605
606 - If 'uniform', `n_knots` number of knots are distributed uniformly
607 from min to max values of the features.
608 - If 'quantile', they are distributed uniformly along the quantiles of
609 the features.
610 - If an array-like is given, it directly specifies the sorted knot
611 positions including the boundary knots. Note that, internally,
612 `degree` number of knots are added before the first knot, the same
613 after the last knot.
614
615 extrapolation : {'error', 'constant', 'linear', 'continue', 'periodic'}, \
616 default='constant'
617 If 'error', values outside the min and max values of the training
618 features raises a `ValueError`. If 'constant', the value of the
619 splines at minimum and maximum value of the features is used as
620 constant extrapolation. If 'linear', a linear extrapolation is used.
621 If 'continue', the splines are extrapolated as is, i.e. option
622 `extrapolate=True` in :class:`scipy.interpolate.BSpline`. If
623 'periodic', periodic splines with a periodicity equal to the distance
624 between the first and last knot are used. Periodic splines enforce
625 equal function values and derivatives at the first and last knot.
626 For example, this makes it possible to avoid introducing an arbitrary
627 jump between Dec 31st and Jan 1st in spline features derived from a
628 naturally periodic "day-of-year" input feature. In this case it is
629 recommended to manually set the knot values to control the period.
630
631 include_bias : bool, default=True
632 If False, then the last spline element inside the data range
633 of a feature is dropped. As B-splines sum to one over the spline basis
634 functions for each data point, they implicitly include a bias term,
635 i.e. a column of ones. It acts as an intercept term in a linear models.
636
637 order : {'C', 'F'}, default='C'
638 Order of output array in the dense case. `'F'` order is faster to compute, but
639 may slow down subsequent estimators.
640
641 sparse_output : bool, default=False
642 Will return sparse CSR matrix if set True else will return an array.
643
644 .. versionadded:: 1.2
645
646 Attributes
647 ----------
648 bsplines_ : list of shape (n_features,)
649 List of BSplines objects, one for each feature.
650
651 n_features_in_ : int
652 The total number of input features.
653
654 feature_names_in_ : ndarray of shape (`n_features_in_`,)
655 Names of features seen during :term:`fit`. Defined only when `X`
656 has feature names that are all strings.
657
658 .. versionadded:: 1.0
659
660 n_features_out_ : int
661 The total number of output features, which is computed as
662 `n_features * n_splines`, where `n_splines` is
663 the number of bases elements of the B-splines,
664 `n_knots + degree - 1` for non-periodic splines and
665 `n_knots - 1` for periodic ones.
666 If `include_bias=False`, then it is only
667 `n_features * (n_splines - 1)`.
668
669 See Also
670 --------
671 KBinsDiscretizer : Transformer that bins continuous data into intervals.
672
673 PolynomialFeatures : Transformer that generates polynomial and interaction
674 features.
675
676 Notes
677 -----
678 High degrees and a high number of knots can cause overfitting.
679
680 See :ref:`examples/linear_model/plot_polynomial_interpolation.py
681 <sphx_glr_auto_examples_linear_model_plot_polynomial_interpolation.py>`.
682
683 Examples
684 --------
685 >>> import numpy as np
686 >>> from sklearn.preprocessing import SplineTransformer
687 >>> X = np.arange(6).reshape(6, 1)
688 >>> spline = SplineTransformer(degree=2, n_knots=3)
689 >>> spline.fit_transform(X)
690 array([[0.5 , 0.5 , 0. , 0. ],
691 [0.18, 0.74, 0.08, 0. ],
692 [0.02, 0.66, 0.32, 0. ],
693 [0. , 0.32, 0.66, 0.02],
694 [0. , 0.08, 0.74, 0.18],
695 [0. , 0. , 0.5 , 0.5 ]])
696 """
697
698 _parameter_constraints: dict = {
699 "n_knots": [Interval(Integral, 2, None, closed="left")],
700 "degree": [Interval(Integral, 0, None, closed="left")],
701 "knots": [StrOptions({"uniform", "quantile"}), "array-like"],
702 "extrapolation": [
703 StrOptions({"error", "constant", "linear", "continue", "periodic"})
704 ],
705 "include_bias": ["boolean"],
706 "order": [StrOptions({"C", "F"})],
707 "sparse_output": ["boolean"],
708 }
709
710 def __init__(
711 self,
712 n_knots=5,
713 degree=3,
714 *,
715 knots="uniform",
716 extrapolation="constant",
717 include_bias=True,
718 order="C",
719 sparse_output=False,
720 ):
721 self.n_knots = n_knots
722 self.degree = degree
723 self.knots = knots
724 self.extrapolation = extrapolation
725 self.include_bias = include_bias
726 self.order = order
727 self.sparse_output = sparse_output
728
729 @staticmethod
730 def _get_base_knot_positions(X, n_knots=10, knots="uniform", sample_weight=None):
731 """Calculate base knot positions.
732
733 Base knots such that first knot <= feature <= last knot. For the
734 B-spline construction with scipy.interpolate.BSpline, 2*degree knots
735 beyond the base interval are added.
736
737 Returns
738 -------
739 knots : ndarray of shape (n_knots, n_features), dtype=np.float64
740 Knot positions (points) of base interval.
741 """
742 if knots == "quantile":
743 percentile_ranks = 100 * np.linspace(
744 start=0, stop=1, num=n_knots, dtype=np.float64
745 )
746
747 if sample_weight is None:
748 knots = np.percentile(X, percentile_ranks, axis=0)
749 else:
750 knots = np.array(
751 [
752 _weighted_percentile(X, sample_weight, percentile_rank)
753 for percentile_rank in percentile_ranks
754 ]
755 )
756
757 else:
758 # knots == 'uniform':
759 # Note that the variable `knots` has already been validated and
760 # `else` is therefore safe.
761 # Disregard observations with zero weight.
762 mask = slice(None, None, 1) if sample_weight is None else sample_weight > 0
763 x_min = np.amin(X[mask], axis=0)
764 x_max = np.amax(X[mask], axis=0)
765
766 knots = np.linspace(
767 start=x_min,
768 stop=x_max,
769 num=n_knots,
770 endpoint=True,
771 dtype=np.float64,
772 )
773
774 return knots
775
776 def get_feature_names_out(self, input_features=None):
777 """Get output feature names for transformation.
778
779 Parameters
780 ----------
781 input_features : array-like of str or None, default=None
782 Input features.
783
784 - If `input_features` is `None`, then `feature_names_in_` is
785 used as feature names in. If `feature_names_in_` is not defined,
786 then the following input feature names are generated:
787 `["x0", "x1", ..., "x(n_features_in_ - 1)"]`.
788 - If `input_features` is an array-like, then `input_features` must
789 match `feature_names_in_` if `feature_names_in_` is defined.
790
791 Returns
792 -------
793 feature_names_out : ndarray of str objects
794 Transformed feature names.
795 """
796 check_is_fitted(self, "n_features_in_")
797 n_splines = self.bsplines_[0].c.shape[1]
798
799 input_features = _check_feature_names_in(self, input_features)
800 feature_names = []
801 for i in range(self.n_features_in_):
802 for j in range(n_splines - 1 + self.include_bias):
803 feature_names.append(f"{input_features[i]}_sp_{j}")
804 return np.asarray(feature_names, dtype=object)
805
806 @_fit_context(prefer_skip_nested_validation=True)
807 def fit(self, X, y=None, sample_weight=None):
808 """Compute knot positions of splines.
809
810 Parameters
811 ----------
812 X : array-like of shape (n_samples, n_features)
813 The data.
814
815 y : None
816 Ignored.
817
818 sample_weight : array-like of shape (n_samples,), default = None
819 Individual weights for each sample. Used to calculate quantiles if
820 `knots="quantile"`. For `knots="uniform"`, zero weighted
821 observations are ignored for finding the min and max of `X`.
822
823 Returns
824 -------
825 self : object
826 Fitted transformer.
827 """
828 X = validate_data(
829 self,
830 X,
831 reset=True,
832 accept_sparse=False,
833 ensure_min_samples=2,
834 ensure_2d=True,
835 )
836 if sample_weight is not None:
837 sample_weight = _check_sample_weight(sample_weight, X, dtype=X.dtype)
838
839 _, n_features = X.shape
840
841 if isinstance(self.knots, str):
842 base_knots = self._get_base_knot_positions(
843 X, n_knots=self.n_knots, knots=self.knots, sample_weight=sample_weight
844 )
845 else:
846 base_knots = check_array(self.knots, dtype=np.float64)
847 if base_knots.shape[0] < 2:
848 raise ValueError("Number of knots, knots.shape[0], must be >= 2.")
849 elif base_knots.shape[1] != n_features:
850 raise ValueError("knots.shape[1] == n_features is violated.")
851 elif not np.all(np.diff(base_knots, axis=0) > 0):
852 raise ValueError("knots must be sorted without duplicates.")
853
854 # number of knots for base interval
855 n_knots = base_knots.shape[0]
856
857 if self.extrapolation == "periodic" and n_knots <= self.degree:
858 raise ValueError(
859 "Periodic splines require degree < n_knots. Got n_knots="
860 f"{n_knots} and degree={self.degree}."
861 )
862
863 # number of splines basis functions
864 if self.extrapolation != "periodic":
865 n_splines = n_knots + self.degree - 1
866 else:
867 # periodic splines have self.degree less degrees of freedom
868 n_splines = n_knots - 1
869
870 degree = self.degree
871 n_out = n_features * n_splines
872 # We have to add degree number of knots below, and degree number knots
873 # above the base knots in order to make the spline basis complete.
874 if self.extrapolation == "periodic":
875 # For periodic splines the spacing of the first / last degree knots
876 # needs to be a continuation of the spacing of the last / first
877 # base knots.
878 period = base_knots[-1] - base_knots[0]
879 knots = np.r_[
880 base_knots[-(degree + 1) : -1] - period,
881 base_knots,
882 base_knots[1 : (degree + 1)] + period,
883 ]
884
885 else:
886 # Eilers & Marx in "Flexible smoothing with B-splines and
887 # penalties" https://doi.org/10.1214/ss/1038425655 advice
888 # against repeating first and last knot several times, which
889 # would have inferior behaviour at boundaries if combined with
890 # a penalty (hence P-Spline). We follow this advice even if our
891 # splines are unpenalized. Meaning we do not:
892 # knots = np.r_[
893 # np.tile(base_knots.min(axis=0), reps=[degree, 1]),
894 # base_knots,
895 # np.tile(base_knots.max(axis=0), reps=[degree, 1])
896 # ]
897 # Instead, we reuse the distance of the 2 fist/last knots.
898 dist_min = base_knots[1] - base_knots[0]
899 dist_max = base_knots[-1] - base_knots[-2]
900
901 knots = np.r_[
902 np.linspace(
903 base_knots[0] - degree * dist_min,
904 base_knots[0] - dist_min,
905 num=degree,
906 ),
907 base_knots,
908 np.linspace(
909 base_knots[-1] + dist_max,
910 base_knots[-1] + degree * dist_max,
911 num=degree,
912 ),
913 ]
914
915 # With a diagonal coefficient matrix, we get back the spline basis
916 # elements, i.e. the design matrix of the spline.
917 # Note, BSpline appreciates C-contiguous float64 arrays as c=coef.
918 coef = np.eye(n_splines, dtype=np.float64)
919 if self.extrapolation == "periodic":
920 coef = np.concatenate((coef, coef[:degree, :]))
921
922 extrapolate = self.extrapolation in ["periodic", "continue"]
923
924 bsplines = [
925 BSpline.construct_fast(
926 knots[:, i], coef, self.degree, extrapolate=extrapolate
927 )
928 for i in range(n_features)
929 ]
930 self.bsplines_ = bsplines
931
932 self.n_features_out_ = n_out - n_features * (1 - self.include_bias)
933 return self
934
935 def transform(self, X):
936 """Transform each feature data to B-splines.
937
938 Parameters
939 ----------
940 X : array-like of shape (n_samples, n_features)
941 The data to transform.
942
943 Returns
944 -------
945 XBS : {ndarray, sparse matrix} of shape (n_samples, n_features * n_splines)
946 The matrix of features, where n_splines is the number of bases
947 elements of the B-splines, n_knots + degree - 1.
948 """
949 check_is_fitted(self)
950
951 X = validate_data(self, X, reset=False, accept_sparse=False, ensure_2d=True)
952
953 n_samples, n_features = X.shape
954 n_splines = self.bsplines_[0].c.shape[1]
955 degree = self.degree
956
957 # TODO: Remove this condition, once scipy 1.10 is the minimum version.
958 # Only scipy => 1.10 supports design_matrix(.., extrapolate=..).
959 # The default (implicit in scipy < 1.10) is extrapolate=False.
960 scipy_1_10 = sp_version >= parse_version("1.10.0")
961 # Note: self.bsplines_[0].extrapolate is True for extrapolation in
962 # ["periodic", "continue"]
963 if scipy_1_10:
964 use_sparse = self.sparse_output
965 kwargs_extrapolate = {"extrapolate": self.bsplines_[0].extrapolate}
966 else:
967 use_sparse = self.sparse_output and not self.bsplines_[0].extrapolate
968 kwargs_extrapolate = dict()
969
970 # Note that scipy BSpline returns float64 arrays and converts input
971 # x=X[:, i] to c-contiguous float64.
972 n_out = self.n_features_out_ + n_features * (1 - self.include_bias)
973 if X.dtype in FLOAT_DTYPES:
974 dtype = X.dtype
975 else:
976 dtype = np.float64
977 if use_sparse:
978 output_list = []
979 else:
980 XBS = np.zeros((n_samples, n_out), dtype=dtype, order=self.order)
981
982 for i in range(n_features):
983 spl = self.bsplines_[i]
984
985 if self.extrapolation in ("continue", "error", "periodic"):
986 if self.extrapolation == "periodic":
987 # With periodic extrapolation we map x to the segment
988 # [spl.t[k], spl.t[n]].
989 # This is equivalent to BSpline(.., extrapolate="periodic")
990 # for scipy>=1.0.0.
991 n = spl.t.size - spl.k - 1
992 # Assign to new array to avoid inplace operation
993 x = spl.t[spl.k] + (X[:, i] - spl.t[spl.k]) % (
994 spl.t[n] - spl.t[spl.k]
995 )
996 else:
997 x = X[:, i]
998
999 if use_sparse:
1000 XBS_sparse = BSpline.design_matrix(
1001 x, spl.t, spl.k, **kwargs_extrapolate
1002 )
1003 if self.extrapolation == "periodic":
1004 # See the construction of coef in fit. We need to add the last
1005 # degree spline basis function to the first degree ones and
1006 # then drop the last ones.
1007 # Note: See comment about SparseEfficiencyWarning below.
1008 XBS_sparse = XBS_sparse.tolil()
1009 XBS_sparse[:, :degree] += XBS_sparse[:, -degree:]
1010 XBS_sparse = XBS_sparse[:, :-degree]
1011 else:
1012 XBS[:, (i * n_splines) : ((i + 1) * n_splines)] = spl(x)
1013 else: # extrapolation in ("constant", "linear")
1014 xmin, xmax = spl.t[degree], spl.t[-degree - 1]
1015 # spline values at boundaries
1016 f_min, f_max = spl(xmin), spl(xmax)
1017 mask = (xmin <= X[:, i]) & (X[:, i] <= xmax)
1018 if use_sparse:
1019 mask_inv = ~mask
1020 x = X[:, i].copy()
1021 # Set some arbitrary values outside boundary that will be reassigned
1022 # later.
1023 x[mask_inv] = spl.t[self.degree]
1024 XBS_sparse = BSpline.design_matrix(x, spl.t, spl.k)
1025 # Note: Without converting to lil_matrix we would get:
1026 # scipy.sparse._base.SparseEfficiencyWarning: Changing the sparsity
1027 # structure of a csr_matrix is expensive. lil_matrix is more
1028 # efficient.
1029 if np.any(mask_inv):
1030 XBS_sparse = XBS_sparse.tolil()
1031 XBS_sparse[mask_inv, :] = 0
1032 else:
1033 XBS[mask, (i * n_splines) : ((i + 1) * n_splines)] = spl(X[mask, i])
1034
1035 # Note for extrapolation:
1036 # 'continue' is already returned as is by scipy BSplines
1037 if self.extrapolation == "error":
1038 # BSpline with extrapolate=False does not raise an error, but
1039 # outputs np.nan.
1040 if (use_sparse and np.any(np.isnan(XBS_sparse.data))) or (
1041 not use_sparse
1042 and np.any(
1043 np.isnan(XBS[:, (i * n_splines) : ((i + 1) * n_splines)])
1044 )
1045 ):
1046 raise ValueError(
1047 "X contains values beyond the limits of the knots."
1048 )
1049 elif self.extrapolation == "constant":
1050 # Set all values beyond xmin and xmax to the value of the
1051 # spline basis functions at those two positions.
1052 # Only the first degree and last degree number of splines
1053 # have non-zero values at the boundaries.
1054
1055 mask = X[:, i] < xmin
1056 if np.any(mask):
1057 if use_sparse:
1058 # Note: See comment about SparseEfficiencyWarning above.
1059 XBS_sparse = XBS_sparse.tolil()
1060 XBS_sparse[mask, :degree] = f_min[:degree]
1061
1062 else:
1063 XBS[mask, (i * n_splines) : (i * n_splines + degree)] = f_min[
1064 :degree
1065 ]
1066
1067 mask = X[:, i] > xmax
1068 if np.any(mask):
1069 if use_sparse:
1070 # Note: See comment about SparseEfficiencyWarning above.
1071 XBS_sparse = XBS_sparse.tolil()
1072 XBS_sparse[mask, -degree:] = f_max[-degree:]
1073 else:
1074 XBS[
1075 mask,
1076 ((i + 1) * n_splines - degree) : ((i + 1) * n_splines),
1077 ] = f_max[-degree:]
1078
1079 elif self.extrapolation == "linear":
1080 # Continue the degree first and degree last spline bases
1081 # linearly beyond the boundaries, with slope = derivative at
1082 # the boundary.
1083 # Note that all others have derivative = value = 0 at the
1084 # boundaries.
1085
1086 # spline derivatives = slopes at boundaries
1087 fp_min, fp_max = spl(xmin, nu=1), spl(xmax, nu=1)
1088 # Compute the linear continuation.
1089 if degree <= 1:
1090 # For degree=1, the derivative of 2nd spline is not zero at
1091 # boundary. For degree=0 it is the same as 'constant'.
1092 degree += 1
1093 for j in range(degree):
1094 mask = X[:, i] < xmin
1095 if np.any(mask):
1096 linear_extr = f_min[j] + (X[mask, i] - xmin) * fp_min[j]
1097 if use_sparse:
1098 # Note: See comment about SparseEfficiencyWarning above.
1099 XBS_sparse = XBS_sparse.tolil()
1100 XBS_sparse[mask, j] = linear_extr
1101 else:
1102 XBS[mask, i * n_splines + j] = linear_extr
1103
1104 mask = X[:, i] > xmax
1105 if np.any(mask):
1106 k = n_splines - 1 - j
1107 linear_extr = f_max[k] + (X[mask, i] - xmax) * fp_max[k]
1108 if use_sparse:
1109 # Note: See comment about SparseEfficiencyWarning above.
1110 XBS_sparse = XBS_sparse.tolil()
1111 XBS_sparse[mask, k : k + 1] = linear_extr[:, None]
1112 else:
1113 XBS[mask, i * n_splines + k] = linear_extr
1114
1115 if use_sparse:
1116 XBS_sparse = XBS_sparse.tocsr()
1117 output_list.append(XBS_sparse)
1118
1119 if use_sparse:
1120 # TODO: Remove this conditional error when the minimum supported version of
1121 # SciPy is 1.9.2
1122 # `scipy.sparse.hstack` breaks in scipy<1.9.2
1123 # when `n_features_out_ > max_int32`
1124 max_int32 = np.iinfo(np.int32).max
1125 all_int32 = True
1126 for mat in output_list:
1127 all_int32 &= mat.indices.dtype == np.int32
1128 if (
1129 sp_version < parse_version("1.9.2")
1130 and self.n_features_out_ > max_int32
1131 and all_int32
1132 ):
1133 raise ValueError(
1134 "In scipy versions `<1.9.2`, the function `scipy.sparse.hstack`"
1135 " produces negative columns when:\n1. The output shape contains"
1136 " `n_cols` too large to be represented by a 32bit signed"
1137 " integer.\n. All sub-matrices to be stacked have indices of"
1138 " dtype `np.int32`.\nTo avoid this error, either use a version"
1139 " of scipy `>=1.9.2` or alter the `SplineTransformer`"
1140 " transformer to produce fewer than 2^31 output features"
1141 )
1142 XBS = sparse.hstack(output_list, format="csr")
1143 elif self.sparse_output:
1144 # TODO: Remove ones scipy 1.10 is the minimum version. See comments above.
1145 XBS = sparse.csr_matrix(XBS)
1146
1147 if self.include_bias:
1148 return XBS
1149 else:
1150 # We throw away one spline basis per feature.
1151 # We chose the last one.
1152 indices = [j for j in range(XBS.shape[1]) if (j + 1) % n_splines != 0]
1153 return XBS[:, indices]
1154