1use scirs2_core::ndarray::{Array1, Array2};
4use sklears_core::{
5 error::{Result, SklearsError},
6 prelude::{Fit, Transform},
7 traits::{Estimator, Trained, Untrained},
8 types::Float,
9};
10use std::marker::PhantomData;
11
12#[derive(Debug, Clone)]
14pub enum NormalizationMethod {
16 None,
18 L2,
20 L1,
22 Max,
24 Standard,
26}
27
28#[derive(Debug, Clone)]
30pub enum CoefficientMethod {
32 Multinomial,
34 Unit,
36 SqrtMultinomial,
38}
39
40#[derive(Debug, Clone)]
70pub struct HomogeneousPolynomialFeatures<State = Untrained> {
72 pub degree: u32,
74 pub interaction_only: bool,
76 pub normalization: NormalizationMethod,
78 pub coefficient_method: CoefficientMethod,
80
81 n_input_features_: Option<usize>,
83 n_output_features_: Option<usize>,
84 power_combinations_: Option<Vec<Vec<u32>>>,
85 coefficients_: Option<Vec<Float>>,
86 normalization_params_: Option<(Array1<Float>, Array1<Float>)>, _state: PhantomData<State>,
89}
90
91impl HomogeneousPolynomialFeatures<Untrained> {
92 pub fn new(degree: u32) -> Self {
94 Self {
95 degree,
96 interaction_only: false,
97 normalization: NormalizationMethod::None,
98 coefficient_method: CoefficientMethod::Unit,
99 n_input_features_: None,
100 n_output_features_: None,
101 power_combinations_: None,
102 coefficients_: None,
103 normalization_params_: None,
104 _state: PhantomData,
105 }
106 }
107
108 pub fn interaction_only(mut self, interaction_only: bool) -> Self {
110 self.interaction_only = interaction_only;
111 self
112 }
113
114 pub fn normalization(mut self, method: NormalizationMethod) -> Self {
116 self.normalization = method;
117 self
118 }
119
120 pub fn coefficient_method(mut self, method: CoefficientMethod) -> Self {
122 self.coefficient_method = method;
123 self
124 }
125}
126
127impl Estimator for HomogeneousPolynomialFeatures<Untrained> {
128 type Config = ();
129 type Error = SklearsError;
130 type Float = Float;
131
132 fn config(&self) -> &Self::Config {
133 &()
134 }
135}
136
137impl Fit<Array2<Float>, ()> for HomogeneousPolynomialFeatures<Untrained> {
138 type Fitted = HomogeneousPolynomialFeatures<Trained>;
139
140 fn fit(self, x: &Array2<Float>, _y: &()) -> Result<Self::Fitted> {
141 let (_, n_features) = x.dim();
142
143 if self.degree == 0 {
144 return Err(SklearsError::InvalidInput(
145 "degree must be positive".to_string(),
146 ));
147 }
148
149 let power_combinations = self.generate_homogeneous_combinations(n_features)?;
151
152 let coefficients = self.compute_coefficients(&power_combinations)?;
154
155 let n_output_features = power_combinations.len();
156
157 let normalization_params = match self.normalization {
159 NormalizationMethod::Standard => {
160 Some(self.compute_normalization_params(x, &power_combinations, &coefficients)?)
161 }
162 _ => None,
163 };
164
165 Ok(HomogeneousPolynomialFeatures {
166 degree: self.degree,
167 interaction_only: self.interaction_only,
168 normalization: self.normalization,
169 coefficient_method: self.coefficient_method,
170 n_input_features_: Some(n_features),
171 n_output_features_: Some(n_output_features),
172 power_combinations_: Some(power_combinations),
173 coefficients_: Some(coefficients),
174 normalization_params_: normalization_params,
175 _state: PhantomData,
176 })
177 }
178}
179
180impl HomogeneousPolynomialFeatures<Untrained> {
181 fn generate_homogeneous_combinations(&self, n_features: usize) -> Result<Vec<Vec<u32>>> {
183 let mut combinations = Vec::new();
184 let mut current_combination = vec![0; n_features];
185
186 self.generate_combinations_recursive(
187 n_features,
188 self.degree,
189 0,
190 &mut current_combination,
191 &mut combinations,
192 );
193
194 if self.interaction_only {
196 combinations.retain(|combination| self.is_valid_for_interaction_only(combination));
197 }
198
199 Ok(combinations)
200 }
201
202 fn generate_combinations_recursive(
204 &self,
205 n_features: usize,
206 remaining_degree: u32,
207 feature_idx: usize,
208 current: &mut Vec<u32>,
209 combinations: &mut Vec<Vec<u32>>,
210 ) {
211 if feature_idx == n_features {
212 if remaining_degree == 0 {
213 combinations.push(current.clone());
214 }
215 return;
216 }
217
218 for power in 0..=remaining_degree {
220 current[feature_idx] = power;
221 self.generate_combinations_recursive(
222 n_features,
223 remaining_degree - power,
224 feature_idx + 1,
225 current,
226 combinations,
227 );
228 }
229 current[feature_idx] = 0;
230 }
231
232 fn is_valid_for_interaction_only(&self, combination: &[u32]) -> bool {
234 let non_zero_count = combination.iter().filter(|&&p| p > 0).count();
235 let max_power = combination.iter().max().unwrap_or(&0);
236
237 *max_power == 1 && non_zero_count >= 2
241 }
242
243 fn compute_coefficients(&self, combinations: &[Vec<u32>]) -> Result<Vec<Float>> {
245 let mut coefficients = Vec::new();
246
247 for combination in combinations {
248 let coeff = match self.coefficient_method {
249 CoefficientMethod::Unit => 1.0,
250 CoefficientMethod::Multinomial => self.compute_multinomial_coefficient(combination),
251 CoefficientMethod::SqrtMultinomial => {
252 self.compute_multinomial_coefficient(combination).sqrt()
253 }
254 };
255 coefficients.push(coeff);
256 }
257
258 Ok(coefficients)
259 }
260
261 fn compute_multinomial_coefficient(&self, powers: &[u32]) -> Float {
263 let total_degree = powers.iter().sum::<u32>();
264
265 if total_degree == 0 {
266 return 1.0;
267 }
268
269 let numerator = self.factorial(total_degree);
271 let mut denominator = 1.0;
272
273 for &power in powers {
274 if power > 0 {
275 denominator *= self.factorial(power);
276 }
277 }
278
279 numerator / denominator
280 }
281
282 fn factorial(&self, n: u32) -> Float {
284 if n <= 1 {
285 1.0
286 } else {
287 (1..=n).map(|i| i as Float).product()
288 }
289 }
290
291 fn compute_normalization_params(
293 &self,
294 x: &Array2<Float>,
295 combinations: &[Vec<u32>],
296 coefficients: &[Float],
297 ) -> Result<(Array1<Float>, Array1<Float>)> {
298 let (n_samples, _) = x.dim();
299 let n_features = combinations.len();
300
301 let mut means = Array1::zeros(n_features);
302 let mut stds = Array1::zeros(n_features);
303
304 for i in 0..n_samples {
306 for (j, (combination, &coeff)) in
307 combinations.iter().zip(coefficients.iter()).enumerate()
308 {
309 let feature_value = self.compute_polynomial_value(&x.row(i), combination) * coeff;
310 means[j] += feature_value;
311 }
312 }
313 means /= n_samples as Float;
314
315 for i in 0..n_samples {
317 for (j, (combination, &coeff)) in
318 combinations.iter().zip(coefficients.iter()).enumerate()
319 {
320 let feature_value = self.compute_polynomial_value(&x.row(i), combination) * coeff;
321 let diff = feature_value - means[j];
322 stds[j] += diff * diff;
323 }
324 }
325 stds = stds.mapv(|var: Float| (var / ((n_samples - 1) as Float)).sqrt().max(1e-12));
326
327 Ok((means, stds))
328 }
329
330 fn compute_polynomial_value(
332 &self,
333 sample: &scirs2_core::ndarray::ArrayView1<Float>,
334 powers: &[u32],
335 ) -> Float {
336 let mut value = 1.0;
337 for (i, &power) in powers.iter().enumerate() {
338 if power > 0 && i < sample.len() {
339 value *= sample[i].powi(power as i32);
340 }
341 }
342 value
343 }
344}
345
346impl Transform<Array2<Float>, Array2<Float>> for HomogeneousPolynomialFeatures<Trained> {
347 fn transform(&self, x: &Array2<Float>) -> Result<Array2<Float>> {
348 let (n_samples, n_features) = x.dim();
349 let n_input_features = self.n_input_features_.expect("operation should succeed");
350 let n_output_features = self.n_output_features_.expect("operation should succeed");
351 let combinations = self
352 .power_combinations_
353 .as_ref()
354 .expect("operation should succeed");
355 let coefficients = self
356 .coefficients_
357 .as_ref()
358 .expect("operation should succeed");
359
360 if n_features != n_input_features {
361 return Err(SklearsError::InvalidInput(format!(
362 "X has {} features, but HomogeneousPolynomialFeatures was fitted with {} features",
363 n_features, n_input_features
364 )));
365 }
366
367 let mut result = Array2::zeros((n_samples, n_output_features));
368
369 for i in 0..n_samples {
371 for (j, (combination, &coeff)) in
372 combinations.iter().zip(coefficients.iter()).enumerate()
373 {
374 let mut feature_value = coeff;
375 for (k, &power) in combination.iter().enumerate() {
376 if power > 0 {
377 feature_value *= x[[i, k]].powi(power as i32);
378 }
379 }
380 result[[i, j]] = feature_value;
381 }
382 }
383
384 match &self.normalization {
386 NormalizationMethod::None => {}
387 NormalizationMethod::L2 => {
388 for mut row in result.rows_mut() {
389 let norm = (row.dot(&row)).sqrt();
390 if norm > 1e-12 {
391 row /= norm;
392 }
393 }
394 }
395 NormalizationMethod::L1 => {
396 for mut row in result.rows_mut() {
397 let norm = row.mapv(|v| v.abs()).sum();
398 if norm > 1e-12 {
399 row /= norm;
400 }
401 }
402 }
403 NormalizationMethod::Max => {
404 for mut row in result.rows_mut() {
405 let max_val = row.mapv(|v| v.abs()).fold(0.0_f64, |a: Float, &b| a.max(b));
406 if max_val > 1e-12 {
407 row /= max_val;
408 }
409 }
410 }
411 NormalizationMethod::Standard => {
412 if let Some((ref means, ref stds)) = self.normalization_params_ {
413 for i in 0..n_samples {
414 for j in 0..n_output_features {
415 result[[i, j]] = (result[[i, j]] - means[j]) / stds[j];
416 }
417 }
418 }
419 }
420 }
421
422 Ok(result)
423 }
424}
425
426impl HomogeneousPolynomialFeatures<Trained> {
427 pub fn n_input_features(&self) -> usize {
429 self.n_input_features_.expect("operation should succeed")
430 }
431
432 pub fn n_output_features(&self) -> usize {
434 self.n_output_features_.expect("operation should succeed")
435 }
436
437 pub fn power_combinations(&self) -> &[Vec<u32>] {
439 self.power_combinations_
440 .as_ref()
441 .expect("operation should succeed")
442 }
443
444 pub fn coefficients(&self) -> &[Float] {
446 self.coefficients_
447 .as_ref()
448 .expect("operation should succeed")
449 }
450
451 pub fn normalization_params(&self) -> Option<&(Array1<Float>, Array1<Float>)> {
453 self.normalization_params_.as_ref()
454 }
455
456 pub fn count_homogeneous_terms(
458 degree: u32,
459 n_features: usize,
460 interaction_only: bool,
461 ) -> usize {
462 if degree == 0 {
463 return if interaction_only { 0 } else { 1 };
464 }
465
466 if interaction_only {
467 if degree > n_features as u32 {
469 return 0;
470 }
471 Self::binomial_coefficient(n_features, degree as usize)
473 } else {
474 Self::binomial_coefficient(n_features + degree as usize - 1, degree as usize)
476 }
477 }
478
479 fn binomial_coefficient(n: usize, k: usize) -> usize {
481 if k > n {
482 return 0;
483 }
484 if k == 0 || k == n {
485 return 1;
486 }
487
488 let k = k.min(n - k); let mut result = 1;
490
491 for i in 0..k {
492 result = result * (n - i) / (i + 1);
493 }
494
495 result
496 }
497}
498
499#[allow(non_snake_case)]
500#[cfg(test)]
501mod tests {
502 use super::*;
503 use approx::assert_abs_diff_eq;
504 use scirs2_core::ndarray::array;
505
506 #[test]
507 fn test_homogeneous_polynomial_basic() {
508 let x = array![[1.0, 2.0], [3.0, 4.0]];
509
510 let homo_poly = HomogeneousPolynomialFeatures::new(2);
511 let fitted = homo_poly.fit(&x, &()).expect("operation should succeed");
512 let x_transformed = fitted.transform(&x).expect("operation should succeed");
513
514 assert_eq!(x_transformed.nrows(), 2);
515
516 assert_eq!(x_transformed.ncols(), 3);
518
519 assert_abs_diff_eq!(x_transformed[[0, 0]], 4.0, epsilon = 1e-10); assert_abs_diff_eq!(x_transformed[[0, 1]], 2.0, epsilon = 1e-10); assert_abs_diff_eq!(x_transformed[[0, 2]], 1.0, epsilon = 1e-10); }
525
526 #[test]
527 fn test_homogeneous_polynomial_interaction_only() {
528 let x = array![[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]];
529
530 let homo_poly = HomogeneousPolynomialFeatures::new(2).interaction_only(true);
531 let fitted = homo_poly.fit(&x, &()).expect("operation should succeed");
532 let x_transformed = fitted.transform(&x).expect("operation should succeed");
533
534 assert_eq!(x_transformed.nrows(), 2);
535
536 assert_eq!(x_transformed.ncols(), 3);
538
539 assert_abs_diff_eq!(x_transformed[[0, 0]], 6.0, epsilon = 1e-10); assert_abs_diff_eq!(x_transformed[[0, 1]], 3.0, epsilon = 1e-10); assert_abs_diff_eq!(x_transformed[[0, 2]], 2.0, epsilon = 1e-10); }
545
546 #[test]
547 fn test_homogeneous_polynomial_degree_3() {
548 let x = array![[1.0, 2.0]];
549
550 let homo_poly = HomogeneousPolynomialFeatures::new(3);
551 let fitted = homo_poly.fit(&x, &()).expect("operation should succeed");
552 let x_transformed = fitted.transform(&x).expect("operation should succeed");
553
554 assert_eq!(x_transformed.ncols(), 4);
556
557 assert_abs_diff_eq!(x_transformed[[0, 0]], 8.0, epsilon = 1e-10); assert_abs_diff_eq!(x_transformed[[0, 1]], 4.0, epsilon = 1e-10); assert_abs_diff_eq!(x_transformed[[0, 2]], 2.0, epsilon = 1e-10); assert_abs_diff_eq!(x_transformed[[0, 3]], 1.0, epsilon = 1e-10); }
564
565 #[test]
566 fn test_homogeneous_polynomial_multinomial_coefficients() {
567 let x = array![[1.0, 1.0]];
568
569 let homo_poly = HomogeneousPolynomialFeatures::new(2)
570 .coefficient_method(CoefficientMethod::Multinomial);
571 let fitted = homo_poly.fit(&x, &()).expect("operation should succeed");
572 let x_transformed = fitted.transform(&x).expect("operation should succeed");
573
574 assert_abs_diff_eq!(x_transformed[[0, 0]], 1.0, epsilon = 1e-10);
579 assert_abs_diff_eq!(x_transformed[[0, 1]], 2.0, epsilon = 1e-10);
580 assert_abs_diff_eq!(x_transformed[[0, 2]], 1.0, epsilon = 1e-10);
581 }
582
583 #[test]
584 fn test_homogeneous_polynomial_l2_normalization() {
585 let x = array![[3.0, 4.0]]; let homo_poly =
588 HomogeneousPolynomialFeatures::new(2).normalization(NormalizationMethod::L2);
589 let fitted = homo_poly.fit(&x, &()).expect("operation should succeed");
590 let x_transformed = fitted.transform(&x).expect("operation should succeed");
591
592 let row_norm = (x_transformed.row(0).dot(&x_transformed.row(0))).sqrt();
594 assert_abs_diff_eq!(row_norm, 1.0, epsilon = 1e-10);
595 }
596
597 #[test]
598 fn test_homogeneous_polynomial_l1_normalization() {
599 let x = array![[2.0, 2.0]]; let homo_poly =
602 HomogeneousPolynomialFeatures::new(2).normalization(NormalizationMethod::L1);
603 let fitted = homo_poly.fit(&x, &()).expect("operation should succeed");
604 let x_transformed = fitted.transform(&x).expect("operation should succeed");
605
606 let row_norm = x_transformed.row(0).mapv(|v| v.abs()).sum();
608 assert_abs_diff_eq!(row_norm, 1.0, epsilon = 1e-10);
609 }
610
611 #[test]
612 fn test_homogeneous_polynomial_standard_normalization() {
613 let x = array![[1.0, 2.0], [3.0, 4.0], [5.0, 6.0]];
614
615 let homo_poly =
616 HomogeneousPolynomialFeatures::new(2).normalization(NormalizationMethod::Standard);
617 let fitted = homo_poly.fit(&x, &()).expect("operation should succeed");
618 let x_transformed = fitted.transform(&x).expect("operation should succeed");
619
620 for j in 0..x_transformed.ncols() {
622 let column = x_transformed.column(j);
623 let mean = column.sum() / column.len() as Float;
624 let variance = column.mapv(|v| (v - mean).powi(2)).sum() / (column.len() - 1) as Float;
625
626 assert_abs_diff_eq!(mean, 0.0, epsilon = 1e-10);
627 assert_abs_diff_eq!(variance, 1.0, epsilon = 1e-10);
628 }
629 }
630
631 #[test]
632 fn test_homogeneous_polynomial_count_terms() {
633 assert_eq!(
635 HomogeneousPolynomialFeatures::<Trained>::count_homogeneous_terms(2, 2, false),
636 3
637 );
638 assert_eq!(
639 HomogeneousPolynomialFeatures::<Trained>::count_homogeneous_terms(2, 3, false),
640 6
641 );
642 assert_eq!(
643 HomogeneousPolynomialFeatures::<Trained>::count_homogeneous_terms(3, 2, false),
644 4
645 );
646
647 assert_eq!(
649 HomogeneousPolynomialFeatures::<Trained>::count_homogeneous_terms(2, 3, true),
650 3
651 );
652 assert_eq!(
653 HomogeneousPolynomialFeatures::<Trained>::count_homogeneous_terms(3, 4, true),
654 4
655 );
656 }
657
658 #[test]
659 fn test_homogeneous_polynomial_zero_degree() {
660 let x = array![[1.0, 2.0]];
661 let homo_poly = HomogeneousPolynomialFeatures::new(0);
662 let result = homo_poly.fit(&x, &());
663 assert!(result.is_err());
664 }
665
666 #[test]
667 fn test_homogeneous_polynomial_feature_mismatch() {
668 let x_train = array![[1.0, 2.0], [3.0, 4.0]];
669 let x_test = array![[1.0, 2.0, 3.0]]; let homo_poly = HomogeneousPolynomialFeatures::new(2);
672 let fitted = homo_poly
673 .fit(&x_train, &())
674 .expect("operation should succeed");
675 let result = fitted.transform(&x_test);
676 assert!(result.is_err());
677 }
678
679 #[test]
680 fn test_homogeneous_polynomial_single_feature() {
681 let x = array![[2.0], [3.0]];
682
683 let homo_poly = HomogeneousPolynomialFeatures::new(3);
684 let fitted = homo_poly.fit(&x, &()).expect("operation should succeed");
685 let x_transformed = fitted.transform(&x).expect("operation should succeed");
686
687 assert_eq!(x_transformed.shape(), &[2, 1]);
689 assert_abs_diff_eq!(x_transformed[[0, 0]], 8.0, epsilon = 1e-10); assert_abs_diff_eq!(x_transformed[[1, 0]], 27.0, epsilon = 1e-10); }
692
693 #[test]
694 fn test_homogeneous_polynomial_degree_1() {
695 let x = array![[1.0, 2.0, 3.0]];
696
697 let homo_poly = HomogeneousPolynomialFeatures::new(1);
698 let fitted = homo_poly.fit(&x, &()).expect("operation should succeed");
699 let x_transformed = fitted.transform(&x).expect("operation should succeed");
700
701 assert_eq!(x_transformed.shape(), &[1, 3]);
703 assert_abs_diff_eq!(x_transformed[[0, 0]], 3.0, epsilon = 1e-10); assert_abs_diff_eq!(x_transformed[[0, 1]], 2.0, epsilon = 1e-10); assert_abs_diff_eq!(x_transformed[[0, 2]], 1.0, epsilon = 1e-10); }
707
708 #[test]
709 fn test_homogeneous_polynomial_interaction_high_degree() {
710 let x = array![[1.0, 2.0]];
711
712 let homo_poly = HomogeneousPolynomialFeatures::new(3).interaction_only(true);
714 let fitted = homo_poly.fit(&x, &()).expect("operation should succeed");
715 let x_transformed = fitted.transform(&x).expect("operation should succeed");
716
717 assert_eq!(x_transformed.ncols(), 0);
719 }
720}