sklears_feature_selection/
group_selection.rs

1//! Group-based feature selection methods
2//!
3//! This module provides feature selection algorithms that consider groups of features,
4//! such as categorical variables, domain-specific feature groups, or hierarchical structures.
5
6use crate::base::SelectorMixin;
7use scirs2_core::ndarray::{Array1, Array2, Axis};
8use sklears_core::{
9    error::{Result as SklResult, SklearsError},
10    traits::{Estimator, Fit, Trained, Transform, Untrained},
11    types::Float,
12};
13use std::collections::{HashMap, HashSet};
14use std::marker::PhantomData;
15
16/// Preprocessed data result: (x_processed, y_processed, x_mean, x_scale, y_mean)
17pub type PreprocessedDataResult = SklResult<(
18    Array2<Float>,
19    Array1<Float>,
20    Array1<Float>,
21    Array1<Float>,
22    f64,
23)>;
24
25/// Group LASSO feature selector
26///
27/// Implements group LASSO regularization that can select or exclude entire groups of features.
28/// This is useful for categorical variables encoded as one-hot vectors or domain-specific feature groups.
29#[derive(Debug, Clone)]
30pub struct GroupLassoSelector<State = Untrained> {
31    /// Regularization parameter (alpha) for group LASSO
32    alpha: f64,
33    /// Groups of feature indices. Each group is a vector of feature indices that belong together.
34    feature_groups: Vec<Vec<usize>>,
35    /// Group names for interpretability
36    group_names: Option<Vec<String>>,
37    /// Maximum number of iterations for optimization
38    max_iter: usize,
39    /// Convergence tolerance
40    tol: f64,
41    /// Whether to fit an intercept
42    fit_intercept: bool,
43    /// Whether to normalize features
44    normalize: bool,
45    state: PhantomData<State>,
46    // Trained state
47    coefficients_: Option<Array1<Float>>,
48    group_weights_: Option<Array1<Float>>,
49    selected_groups_: Option<Vec<usize>>,
50    selected_features_: Option<Vec<usize>>,
51    intercept_: Option<f64>,
52}
53
54impl GroupLassoSelector<Untrained> {
55    /// new
56    pub fn new() -> Self {
57        Self {
58            alpha: 1.0,
59            feature_groups: Vec::new(),
60            group_names: None,
61            max_iter: 1000,
62            tol: 1e-4,
63            fit_intercept: true,
64            normalize: true,
65            state: PhantomData,
66            coefficients_: None,
67            group_weights_: None,
68            selected_groups_: None,
69            selected_features_: None,
70            intercept_: None,
71        }
72    }
73
74    /// alpha
75    pub fn alpha(mut self, alpha: f64) -> Self {
76        self.alpha = alpha;
77        self
78    }
79
80    /// feature_groups
81    pub fn feature_groups(mut self, groups: Vec<Vec<usize>>) -> Self {
82        self.feature_groups = groups;
83        self
84    }
85
86    /// group_names
87    pub fn group_names(mut self, names: Vec<String>) -> Self {
88        self.group_names = Some(names);
89        self
90    }
91
92    /// max_iter
93    pub fn max_iter(mut self, max_iter: usize) -> Self {
94        self.max_iter = max_iter;
95        self
96    }
97
98    /// tol
99    pub fn tol(mut self, tol: f64) -> Self {
100        self.tol = tol;
101        self
102    }
103
104    /// fit_intercept
105    pub fn fit_intercept(mut self, fit_intercept: bool) -> Self {
106        self.fit_intercept = fit_intercept;
107        self
108    }
109
110    /// normalize
111    pub fn normalize(mut self, normalize: bool) -> Self {
112        self.normalize = normalize;
113        self
114    }
115}
116
117impl Default for GroupLassoSelector<Untrained> {
118    fn default() -> Self {
119        Self::new()
120    }
121}
122
123impl Estimator for GroupLassoSelector<Untrained> {
124    type Config = ();
125    type Error = SklearsError;
126    type Float = f64;
127
128    fn config(&self) -> &Self::Config {
129        &()
130    }
131}
132
133impl Fit<Array2<Float>, Array1<Float>> for GroupLassoSelector<Untrained> {
134    type Fitted = GroupLassoSelector<Trained>;
135
136    fn fit(self, x: &Array2<Float>, y: &Array1<Float>) -> SklResult<Self::Fitted> {
137        let (n_samples, n_features) = x.dim();
138
139        if n_samples != y.len() {
140            return Err(SklearsError::InvalidInput(
141                "X and y must have the same number of samples".to_string(),
142            ));
143        }
144
145        if self.feature_groups.is_empty() {
146            return Err(SklearsError::InvalidInput(
147                "Feature groups must be specified".to_string(),
148            ));
149        }
150
151        // Validate feature groups
152        let mut all_features = HashSet::new();
153        for group in &self.feature_groups {
154            for &feature in group {
155                if feature >= n_features {
156                    return Err(SklearsError::InvalidInput(
157                        "Feature index out of bounds in groups".to_string(),
158                    ));
159                }
160                if all_features.contains(&feature) {
161                    return Err(SklearsError::InvalidInput(
162                        "Feature appears in multiple groups".to_string(),
163                    ));
164                }
165                all_features.insert(feature);
166            }
167        }
168
169        // Prepare data
170        let (x_processed, y_processed, _feature_means, _feature_stds, y_mean) =
171            self.preprocess_data(x, y)?;
172
173        // Fit group LASSO using coordinate descent
174        let (coefficients, group_weights, selected_groups) =
175            self.fit_group_lasso(&x_processed, &y_processed)?;
176
177        // Determine selected features from selected groups
178        let mut selected_features = Vec::new();
179        for &group_idx in &selected_groups {
180            selected_features.extend_from_slice(&self.feature_groups[group_idx]);
181        }
182        selected_features.sort();
183
184        // Compute intercept if needed
185        let intercept = if self.fit_intercept {
186            Some(y_mean)
187        } else {
188            None
189        };
190
191        Ok(GroupLassoSelector {
192            alpha: self.alpha,
193            feature_groups: self.feature_groups,
194            group_names: self.group_names,
195            max_iter: self.max_iter,
196            tol: self.tol,
197            fit_intercept: self.fit_intercept,
198            normalize: self.normalize,
199            state: PhantomData,
200            coefficients_: Some(coefficients),
201            group_weights_: Some(group_weights),
202            selected_groups_: Some(selected_groups),
203            selected_features_: Some(selected_features),
204            intercept_: intercept,
205        })
206    }
207}
208
209impl GroupLassoSelector<Untrained> {
210    fn preprocess_data(&self, x: &Array2<Float>, y: &Array1<Float>) -> PreprocessedDataResult {
211        let (n_samples, n_features) = x.dim();
212        let mut x_processed = x.clone();
213        let mut y_processed = y.clone();
214
215        // Center target if fitting intercept
216        let y_mean = if self.fit_intercept {
217            let mean = y.mean().unwrap_or(0.0);
218            for i in 0..n_samples {
219                y_processed[i] -= mean;
220            }
221            mean
222        } else {
223            0.0
224        };
225
226        // Center and normalize features
227        let mut feature_means = Array1::zeros(n_features);
228        let mut feature_stds = Array1::ones(n_features);
229
230        for j in 0..n_features {
231            let col = x.column(j);
232            let mean = col.mean().unwrap_or(0.0);
233            feature_means[j] = mean;
234
235            // Center
236            for i in 0..n_samples {
237                x_processed[[i, j]] -= mean;
238            }
239
240            // Normalize if requested
241            if self.normalize {
242                let std = (col.mapv(|v| (v - mean).powi(2)).sum() / n_samples as f64).sqrt();
243                if std > 1e-10 {
244                    feature_stds[j] = std;
245                    for i in 0..n_samples {
246                        x_processed[[i, j]] /= std;
247                    }
248                }
249            }
250        }
251
252        Ok((
253            x_processed,
254            y_processed,
255            feature_means,
256            feature_stds,
257            y_mean,
258        ))
259    }
260
261    fn fit_group_lasso(
262        &self,
263        x: &Array2<Float>,
264        y: &Array1<Float>,
265    ) -> SklResult<(Array1<Float>, Array1<Float>, Vec<usize>)> {
266        let (n_samples, n_features) = x.dim();
267        let n_groups = self.feature_groups.len();
268
269        // Initialize coefficients
270        let mut coefficients = Array1::zeros(n_features);
271        let mut group_weights = Array1::zeros(n_groups);
272
273        // Compute group norms for scaling
274        let mut group_norms = Array1::zeros(n_groups);
275        for (g, group) in self.feature_groups.iter().enumerate() {
276            let group_size = group.len() as f64;
277            group_norms[g] = group_size.sqrt();
278        }
279
280        // Coordinate descent optimization
281        for _iter in 0..self.max_iter {
282            let coefficients_old = coefficients.clone();
283
284            for (g, group) in self.feature_groups.iter().enumerate() {
285                // Compute residual for this group
286                let mut residual = y.clone();
287                for (h, other_group) in self.feature_groups.iter().enumerate() {
288                    if g != h {
289                        for &feature in other_group {
290                            for i in 0..n_samples {
291                                residual[i] -= x[[i, feature]] * coefficients[feature];
292                            }
293                        }
294                    }
295                }
296
297                // Compute group gradient
298                let mut group_gradient = Array1::zeros(group.len());
299                for (idx, &feature) in group.iter().enumerate() {
300                    let mut grad = 0.0;
301                    for i in 0..n_samples {
302                        grad += x[[i, feature]] * residual[i];
303                    }
304                    group_gradient[idx] = grad / n_samples as f64;
305                }
306
307                // Compute group norm
308                let group_norm = group_gradient.mapv(|x| x * x).sum().sqrt();
309
310                // Apply group soft thresholding
311                let threshold = self.alpha * group_norms[g];
312                if group_norm > threshold {
313                    // Group is selected
314                    let shrinkage_factor = 1.0 - threshold / group_norm;
315                    group_weights[g] = group_norm * shrinkage_factor;
316
317                    // Update coefficients for this group
318                    for (idx, &feature) in group.iter().enumerate() {
319                        coefficients[feature] = group_gradient[idx] * shrinkage_factor;
320                    }
321                } else {
322                    // Group is not selected (set to zero)
323                    group_weights[g] = 0.0;
324                    for &feature in group {
325                        coefficients[feature] = 0.0;
326                    }
327                }
328            }
329
330            // Check convergence
331            let diff = (&coefficients - &coefficients_old).mapv(|x| x.abs()).sum();
332            if diff < self.tol {
333                break;
334            }
335        }
336
337        // Determine selected groups
338        let selected_groups: Vec<usize> = group_weights
339            .indexed_iter()
340            .filter_map(|(i, &weight)| if weight.abs() > 1e-10 { Some(i) } else { None })
341            .collect();
342
343        Ok((coefficients, group_weights, selected_groups))
344    }
345}
346
347impl Transform<Array2<Float>> for GroupLassoSelector<Trained> {
348    fn transform(&self, x: &Array2<Float>) -> SklResult<Array2<Float>> {
349        let selected_features = self
350            .selected_features_
351            .as_ref()
352            .expect("operation should succeed");
353        if selected_features.is_empty() {
354            return Err(SklearsError::InvalidInput(
355                "No features were selected".to_string(),
356            ));
357        }
358
359        let selected_indices: Vec<usize> = selected_features.to_vec();
360        Ok(x.select(Axis(1), &selected_indices))
361    }
362}
363
364impl SelectorMixin for GroupLassoSelector<Trained> {
365    fn get_support(&self) -> SklResult<Array1<bool>> {
366        let selected_features = self
367            .selected_features_
368            .as_ref()
369            .expect("operation should succeed");
370        let n_features = self
371            .coefficients_
372            .as_ref()
373            .expect("operation should succeed")
374            .len();
375        let mut support = Array1::from_elem(n_features, false);
376        for &idx in selected_features {
377            if idx < n_features {
378                support[idx] = true;
379            }
380        }
381        Ok(support)
382    }
383
384    fn transform_features(&self, indices: &[usize]) -> SklResult<Vec<usize>> {
385        let selected_features = self
386            .selected_features_
387            .as_ref()
388            .expect("operation should succeed");
389        Ok(indices
390            .iter()
391            .filter_map(|&idx| selected_features.iter().position(|&f| f == idx))
392            .collect())
393    }
394}
395
396impl GroupLassoSelector<Trained> {
397    /// Get the group coefficients (group-level importance scores)
398    pub fn group_weights(&self) -> Option<&Array1<Float>> {
399        self.group_weights_.as_ref()
400    }
401
402    /// Get individual feature coefficients
403    pub fn coefficients(&self) -> Option<&Array1<Float>> {
404        self.coefficients_.as_ref()
405    }
406
407    /// Get selected group indices
408    pub fn selected_groups(&self) -> Option<&Vec<usize>> {
409        self.selected_groups_.as_ref()
410    }
411
412    /// Get the intercept term
413    pub fn intercept(&self) -> Option<f64> {
414        self.intercept_
415    }
416
417    /// Get feature groups that were selected
418    pub fn get_selected_group_names(&self) -> Vec<String> {
419        if let (Some(selected_groups), Some(group_names)) =
420            (&self.selected_groups_, &self.group_names)
421        {
422            selected_groups
423                .iter()
424                .filter_map(|&idx| group_names.get(idx).cloned())
425                .collect()
426        } else {
427            self.selected_groups_
428                .as_ref()
429                .unwrap_or(&Vec::new())
430                .iter()
431                .map(|&idx| format!("Group_{}", idx))
432                .collect()
433        }
434    }
435
436    /// Predict using the fitted model
437    pub fn predict(&self, x: &Array2<Float>) -> SklResult<Array1<Float>> {
438        let coefficients = self
439            .coefficients_
440            .as_ref()
441            .expect("operation should succeed");
442        let (n_samples, n_features) = x.dim();
443
444        if n_features != coefficients.len() {
445            return Err(SklearsError::InvalidInput(
446                "X has wrong number of features".to_string(),
447            ));
448        }
449
450        let mut predictions = Array1::zeros(n_samples);
451
452        for i in 0..n_samples {
453            let mut pred = 0.0;
454            for j in 0..n_features {
455                pred += x[[i, j]] * coefficients[j];
456            }
457
458            if let Some(intercept) = self.intercept_ {
459                pred += intercept;
460            }
461
462            predictions[i] = pred;
463        }
464
465        Ok(predictions)
466    }
467}
468
469/// Sparse Group LASSO selector
470///
471/// Combines group LASSO with element-wise LASSO to encourage both group sparsity
472/// and within-group sparsity.
473#[derive(Debug, Clone)]
474pub struct SparseGroupLassoSelector<State = Untrained> {
475    /// Group regularization parameter
476    alpha_group: f64,
477    /// Element-wise regularization parameter
478    alpha_element: f64,
479    /// Groups of feature indices
480    feature_groups: Vec<Vec<usize>>,
481    /// Group names for interpretability
482    group_names: Option<Vec<String>>,
483    /// Maximum number of iterations
484    max_iter: usize,
485    /// Convergence tolerance
486    tol: f64,
487    /// Whether to fit intercept
488    fit_intercept: bool,
489    /// Whether to normalize features
490    normalize: bool,
491    state: PhantomData<State>,
492    // Trained state
493    coefficients_: Option<Array1<Float>>,
494    selected_features_: Option<Vec<usize>>,
495    _selected_groups_: Option<Vec<usize>>,
496    _intercept_: Option<f64>,
497}
498
499impl SparseGroupLassoSelector<Untrained> {
500    /// new
501    pub fn new() -> Self {
502        Self {
503            alpha_group: 1.0,
504            alpha_element: 1.0,
505            feature_groups: Vec::new(),
506            group_names: None,
507            max_iter: 1000,
508            tol: 1e-4,
509            fit_intercept: true,
510            normalize: true,
511            state: PhantomData,
512            coefficients_: None,
513            selected_features_: None,
514            _selected_groups_: None,
515            _intercept_: None,
516        }
517    }
518
519    /// alpha_group
520    pub fn alpha_group(mut self, alpha: f64) -> Self {
521        self.alpha_group = alpha;
522        self
523    }
524
525    /// alpha_element
526    pub fn alpha_element(mut self, alpha: f64) -> Self {
527        self.alpha_element = alpha;
528        self
529    }
530
531    /// feature_groups
532    pub fn feature_groups(mut self, groups: Vec<Vec<usize>>) -> Self {
533        self.feature_groups = groups;
534        self
535    }
536
537    /// group_names
538    pub fn group_names(mut self, names: Vec<String>) -> Self {
539        self.group_names = Some(names);
540        self
541    }
542}
543
544impl Default for SparseGroupLassoSelector<Untrained> {
545    fn default() -> Self {
546        Self::new()
547    }
548}
549
550impl Estimator for SparseGroupLassoSelector<Untrained> {
551    type Config = ();
552    type Error = SklearsError;
553    type Float = f64;
554
555    fn config(&self) -> &Self::Config {
556        &()
557    }
558}
559
560impl Fit<Array2<Float>, Array1<Float>> for SparseGroupLassoSelector<Untrained> {
561    type Fitted = SparseGroupLassoSelector<Trained>;
562
563    fn fit(self, x: &Array2<Float>, y: &Array1<Float>) -> SklResult<Self::Fitted> {
564        // For simplicity, implement as a combination of group and element-wise penalties
565        // In practice, this would use a more sophisticated optimization algorithm
566
567        let (n_samples, _n_features) = x.dim();
568        if n_samples != y.len() {
569            return Err(SklearsError::InvalidInput(
570                "X and y must have the same number of samples".to_string(),
571            ));
572        }
573
574        if self.feature_groups.is_empty() {
575            return Err(SklearsError::InvalidInput(
576                "Feature groups must be specified".to_string(),
577            ));
578        }
579
580        // Use a simplified implementation: first apply group LASSO, then element-wise LASSO
581        let group_lasso = GroupLassoSelector::new()
582            .alpha(self.alpha_group)
583            .feature_groups(self.feature_groups.clone())
584            .fit_intercept(self.fit_intercept)
585            .normalize(self.normalize);
586
587        let fitted_group_lasso = group_lasso.fit(x, y)?;
588        let mut coefficients = fitted_group_lasso
589            .coefficients()
590            .expect("operation should succeed")
591            .clone();
592
593        // Apply element-wise soft thresholding
594        for coeff in coefficients.iter_mut() {
595            *coeff = soft_threshold(*coeff, self.alpha_element);
596        }
597
598        // Determine selected features and groups
599        let selected_features: Vec<usize> = coefficients
600            .indexed_iter()
601            .filter_map(|(i, &coeff)| if coeff.abs() > 1e-10 { Some(i) } else { None })
602            .collect();
603
604        let mut selected_groups = Vec::new();
605        for (g, group) in self.feature_groups.iter().enumerate() {
606            if group.iter().any(|&f| selected_features.contains(&f)) {
607                selected_groups.push(g);
608            }
609        }
610
611        Ok(SparseGroupLassoSelector {
612            alpha_group: self.alpha_group,
613            alpha_element: self.alpha_element,
614            feature_groups: self.feature_groups,
615            group_names: self.group_names,
616            max_iter: self.max_iter,
617            tol: self.tol,
618            fit_intercept: self.fit_intercept,
619            normalize: self.normalize,
620            state: PhantomData,
621            coefficients_: Some(coefficients),
622            selected_features_: Some(selected_features),
623            _selected_groups_: Some(selected_groups),
624            _intercept_: fitted_group_lasso.intercept(),
625        })
626    }
627}
628
629impl Transform<Array2<Float>> for SparseGroupLassoSelector<Trained> {
630    fn transform(&self, x: &Array2<Float>) -> SklResult<Array2<Float>> {
631        let selected_features = self
632            .selected_features_
633            .as_ref()
634            .expect("operation should succeed");
635        if selected_features.is_empty() {
636            return Err(SklearsError::InvalidInput(
637                "No features were selected".to_string(),
638            ));
639        }
640
641        let selected_indices: Vec<usize> = selected_features.to_vec();
642        Ok(x.select(Axis(1), &selected_indices))
643    }
644}
645
646impl SelectorMixin for SparseGroupLassoSelector<Trained> {
647    fn get_support(&self) -> SklResult<Array1<bool>> {
648        let selected_features = self
649            .selected_features_
650            .as_ref()
651            .expect("operation should succeed");
652        let n_features = self
653            .coefficients_
654            .as_ref()
655            .expect("operation should succeed")
656            .len();
657        let mut support = Array1::from_elem(n_features, false);
658        for &idx in selected_features {
659            if idx < n_features {
660                support[idx] = true;
661            }
662        }
663        Ok(support)
664    }
665
666    fn transform_features(&self, indices: &[usize]) -> SklResult<Vec<usize>> {
667        let selected_features = self
668            .selected_features_
669            .as_ref()
670            .expect("operation should succeed");
671        Ok(indices
672            .iter()
673            .filter_map(|&idx| selected_features.iter().position(|&f| f == idx))
674            .collect())
675    }
676}
677
678// Helper function for soft thresholding
679fn soft_threshold(x: f64, lambda: f64) -> f64 {
680    if x > lambda {
681        x - lambda
682    } else if x < -lambda {
683        x + lambda
684    } else {
685        0.0
686    }
687}
688
689/// Hierarchical structured sparsity selector
690///
691/// Implements structured sparsity regularization for tree-like hierarchical structures.
692/// This is useful for nested feature groups where parent-child relationships exist.
693#[derive(Debug, Clone)]
694pub struct HierarchicalStructuredSparsitySelector<State = Untrained> {
695    /// Regularization parameter for hierarchical penalty
696    alpha: f64,
697    /// Hierarchical structure as adjacency list (parent -> children)
698    hierarchy: HashMap<usize, Vec<usize>>,
699    /// Mapping from feature indices to hierarchy nodes
700    feature_to_node: HashMap<usize, usize>,
701    /// Maximum number of iterations
702    max_iter: usize,
703    /// Convergence tolerance
704    tol: f64,
705    /// Whether to fit intercept
706    fit_intercept: bool,
707    /// Whether to normalize features
708    normalize: bool,
709    state: PhantomData<State>,
710    // Trained state
711    coefficients_: Option<Array1<Float>>,
712    selected_features_: Option<Vec<usize>>,
713    _selected_nodes_: Option<Vec<usize>>,
714    _intercept_: Option<f64>,
715}
716
717impl HierarchicalStructuredSparsitySelector<Untrained> {
718    /// new
719    pub fn new() -> Self {
720        Self {
721            alpha: 1.0,
722            hierarchy: HashMap::new(),
723            feature_to_node: HashMap::new(),
724            max_iter: 1000,
725            tol: 1e-4,
726            fit_intercept: true,
727            normalize: true,
728            state: PhantomData,
729            coefficients_: None,
730            selected_features_: None,
731            _selected_nodes_: None,
732            _intercept_: None,
733        }
734    }
735
736    /// alpha
737    pub fn alpha(mut self, alpha: f64) -> Self {
738        self.alpha = alpha;
739        self
740    }
741
742    /// hierarchy
743    pub fn hierarchy(mut self, hierarchy: HashMap<usize, Vec<usize>>) -> Self {
744        self.hierarchy = hierarchy;
745        self
746    }
747
748    /// feature_to_node
749    pub fn feature_to_node(mut self, mapping: HashMap<usize, usize>) -> Self {
750        self.feature_to_node = mapping;
751        self
752    }
753
754    /// max_iter
755    pub fn max_iter(mut self, max_iter: usize) -> Self {
756        self.max_iter = max_iter;
757        self
758    }
759
760    /// tol
761    pub fn tol(mut self, tol: f64) -> Self {
762        self.tol = tol;
763        self
764    }
765
766    /// fit_intercept
767    pub fn fit_intercept(mut self, fit_intercept: bool) -> Self {
768        self.fit_intercept = fit_intercept;
769        self
770    }
771
772    /// normalize
773    pub fn normalize(mut self, normalize: bool) -> Self {
774        self.normalize = normalize;
775        self
776    }
777}
778
779impl Default for HierarchicalStructuredSparsitySelector<Untrained> {
780    fn default() -> Self {
781        Self::new()
782    }
783}
784
785impl Estimator for HierarchicalStructuredSparsitySelector<Untrained> {
786    type Config = ();
787    type Error = SklearsError;
788    type Float = f64;
789
790    fn config(&self) -> &Self::Config {
791        &()
792    }
793}
794
795impl Fit<Array2<Float>, Array1<Float>> for HierarchicalStructuredSparsitySelector<Untrained> {
796    type Fitted = HierarchicalStructuredSparsitySelector<Trained>;
797
798    fn fit(self, x: &Array2<Float>, y: &Array1<Float>) -> SklResult<Self::Fitted> {
799        let (n_samples, _n_features) = x.dim();
800
801        if n_samples != y.len() {
802            return Err(SklearsError::InvalidInput(
803                "X and y must have the same number of samples".to_string(),
804            ));
805        }
806
807        if self.hierarchy.is_empty() || self.feature_to_node.is_empty() {
808            return Err(SklearsError::InvalidInput(
809                "Hierarchical structure and feature mapping must be specified".to_string(),
810            ));
811        }
812
813        // Preprocess data
814        let (x_processed, y_processed, _, _, y_mean) = self.preprocess_data(x, y)?;
815
816        // Fit hierarchical structured sparsity using proximal gradient descent
817        let (coefficients, selected_nodes) =
818            self.fit_hierarchical_sparsity(&x_processed, &y_processed)?;
819
820        // Determine selected features from selected nodes
821        let selected_features: Vec<usize> = self
822            .feature_to_node
823            .iter()
824            .filter_map(|(&feature, &node)| {
825                if selected_nodes.contains(&node) {
826                    Some(feature)
827                } else {
828                    None
829                }
830            })
831            .collect();
832
833        let intercept = if self.fit_intercept {
834            Some(y_mean)
835        } else {
836            None
837        };
838
839        Ok(HierarchicalStructuredSparsitySelector {
840            alpha: self.alpha,
841            hierarchy: self.hierarchy,
842            feature_to_node: self.feature_to_node,
843            max_iter: self.max_iter,
844            tol: self.tol,
845            fit_intercept: self.fit_intercept,
846            normalize: self.normalize,
847            state: PhantomData,
848            coefficients_: Some(coefficients),
849            selected_features_: Some(selected_features),
850            _selected_nodes_: Some(selected_nodes),
851            _intercept_: intercept,
852        })
853    }
854}
855
856impl HierarchicalStructuredSparsitySelector<Untrained> {
857    fn preprocess_data(&self, x: &Array2<Float>, y: &Array1<Float>) -> PreprocessedDataResult {
858        let (n_samples, n_features) = x.dim();
859        let mut x_processed = x.clone();
860        let mut y_processed = y.clone();
861
862        // Center target if fitting intercept
863        let y_mean = if self.fit_intercept {
864            let mean = y.mean().unwrap_or(0.0);
865            for i in 0..n_samples {
866                y_processed[i] -= mean;
867            }
868            mean
869        } else {
870            0.0
871        };
872
873        // Center and normalize features
874        let mut feature_means = Array1::zeros(n_features);
875        let mut feature_stds = Array1::ones(n_features);
876
877        for j in 0..n_features {
878            let col = x.column(j);
879            let mean = col.mean().unwrap_or(0.0);
880            feature_means[j] = mean;
881
882            // Center
883            for i in 0..n_samples {
884                x_processed[[i, j]] -= mean;
885            }
886
887            // Normalize if requested
888            if self.normalize {
889                let std = (col.mapv(|v| (v - mean).powi(2)).sum() / n_samples as f64).sqrt();
890                if std > 1e-10 {
891                    feature_stds[j] = std;
892                    for i in 0..n_samples {
893                        x_processed[[i, j]] /= std;
894                    }
895                }
896            }
897        }
898
899        Ok((
900            x_processed,
901            y_processed,
902            feature_means,
903            feature_stds,
904            y_mean,
905        ))
906    }
907
908    fn fit_hierarchical_sparsity(
909        &self,
910        x: &Array2<Float>,
911        y: &Array1<Float>,
912    ) -> SklResult<(Array1<Float>, Vec<usize>)> {
913        let (n_samples, n_features) = x.dim();
914        let mut coefficients = Array1::zeros(n_features);
915        let step_size = 0.01;
916
917        for _iter in 0..self.max_iter {
918            let coefficients_old = coefficients.clone();
919
920            // Compute gradient
921            let mut gradient = Array1::<f64>::zeros(n_features);
922            for i in 0..n_samples {
923                let residual = y[i] - x.row(i).dot(&coefficients);
924                for j in 0..n_features {
925                    gradient[j] -= x[[i, j]] * residual / n_samples as f64;
926                }
927            }
928
929            // Gradient descent step
930            for j in 0..n_features {
931                coefficients[j] -= step_size * gradient[j];
932            }
933
934            // Apply hierarchical proximal operator
935            self.hierarchical_proximal_operator(&mut coefficients)?;
936
937            // Check convergence
938            let diff = (&coefficients - &coefficients_old).mapv(|x| x.abs()).sum();
939            if diff < self.tol {
940                break;
941            }
942        }
943
944        // Determine selected nodes
945        let selected_nodes = self.get_selected_nodes(&coefficients);
946
947        Ok((coefficients, selected_nodes))
948    }
949
950    fn hierarchical_proximal_operator(&self, coefficients: &mut Array1<Float>) -> SklResult<()> {
951        // Apply hierarchical constraint: if a parent is not selected, children shouldn't be selected
952        for (&parent, children) in &self.hierarchy {
953            let parent_features: Vec<usize> = self
954                .feature_to_node
955                .iter()
956                .filter_map(|(&feature, &node)| if node == parent { Some(feature) } else { None })
957                .collect();
958
959            let parent_norm = parent_features
960                .iter()
961                .map(|&f| coefficients[f] * coefficients[f])
962                .sum::<f64>()
963                .sqrt();
964
965            // Apply group soft thresholding to parent
966            let threshold = self.alpha * (parent_features.len() as f64).sqrt();
967            if parent_norm <= threshold {
968                // Parent is not selected, zero out parent and all children
969                for &feature in &parent_features {
970                    coefficients[feature] = 0.0;
971                }
972
973                // Zero out all children recursively
974                self.zero_out_children(coefficients, children)?;
975            } else {
976                // Parent is selected, apply soft thresholding
977                let shrinkage = 1.0 - threshold / parent_norm;
978                for &feature in &parent_features {
979                    coefficients[feature] *= shrinkage;
980                }
981            }
982        }
983
984        Ok(())
985    }
986
987    fn zero_out_children(
988        &self,
989        coefficients: &mut Array1<Float>,
990        children: &[usize],
991    ) -> SklResult<()> {
992        for &child in children {
993            let child_features: Vec<usize> = self
994                .feature_to_node
995                .iter()
996                .filter_map(|(&feature, &node)| if node == child { Some(feature) } else { None })
997                .collect();
998
999            for &feature in &child_features {
1000                coefficients[feature] = 0.0;
1001            }
1002
1003            // Recursively zero out grandchildren
1004            if let Some(grandchildren) = self.hierarchy.get(&child) {
1005                self.zero_out_children(coefficients, grandchildren)?;
1006            }
1007        }
1008        Ok(())
1009    }
1010
1011    fn get_selected_nodes(&self, coefficients: &Array1<Float>) -> Vec<usize> {
1012        let mut selected_nodes = HashSet::new();
1013
1014        for (&feature, &node) in &self.feature_to_node {
1015            if coefficients[feature].abs() > 1e-10 {
1016                selected_nodes.insert(node);
1017            }
1018        }
1019
1020        selected_nodes.into_iter().collect()
1021    }
1022}
1023
1024impl Transform<Array2<Float>> for HierarchicalStructuredSparsitySelector<Trained> {
1025    fn transform(&self, x: &Array2<Float>) -> SklResult<Array2<Float>> {
1026        let selected_features = self
1027            .selected_features_
1028            .as_ref()
1029            .expect("operation should succeed");
1030        if selected_features.is_empty() {
1031            return Err(SklearsError::InvalidInput(
1032                "No features were selected".to_string(),
1033            ));
1034        }
1035
1036        let selected_indices: Vec<usize> = selected_features.to_vec();
1037        Ok(x.select(Axis(1), &selected_indices))
1038    }
1039}
1040
1041impl SelectorMixin for HierarchicalStructuredSparsitySelector<Trained> {
1042    fn get_support(&self) -> SklResult<Array1<bool>> {
1043        let selected_features = self
1044            .selected_features_
1045            .as_ref()
1046            .expect("operation should succeed");
1047        let n_features = self
1048            .coefficients_
1049            .as_ref()
1050            .expect("operation should succeed")
1051            .len();
1052        let mut support = Array1::from_elem(n_features, false);
1053        for &idx in selected_features {
1054            if idx < n_features {
1055                support[idx] = true;
1056            }
1057        }
1058        Ok(support)
1059    }
1060
1061    fn transform_features(&self, indices: &[usize]) -> SklResult<Vec<usize>> {
1062        let selected_features = self
1063            .selected_features_
1064            .as_ref()
1065            .expect("operation should succeed");
1066        Ok(indices
1067            .iter()
1068            .filter_map(|&idx| selected_features.iter().position(|&f| f == idx))
1069            .collect())
1070    }
1071}
1072
1073/// Graph-based structured sparsity selector
1074///
1075/// Implements structured sparsity regularization for graph-structured feature relationships.
1076/// This encourages smoothness over connected features in the graph.
1077#[derive(Debug, Clone)]
1078pub struct GraphStructuredSparsitySelector<State = Untrained> {
1079    /// Regularization parameter for graph penalty
1080    alpha: f64,
1081    /// Graph structure as adjacency matrix
1082    adjacency_matrix: Array2<Float>,
1083    /// Maximum number of iterations
1084    max_iter: usize,
1085    /// Convergence tolerance
1086    tol: f64,
1087    /// Whether to fit intercept
1088    fit_intercept: bool,
1089    /// Whether to normalize features
1090    normalize: bool,
1091    state: PhantomData<State>,
1092    // Trained state
1093    coefficients_: Option<Array1<Float>>,
1094    selected_features_: Option<Vec<usize>>,
1095    _intercept_: Option<f64>,
1096}
1097
1098impl GraphStructuredSparsitySelector<Untrained> {
1099    /// new
1100    pub fn new() -> Self {
1101        Self {
1102            alpha: 1.0,
1103            adjacency_matrix: Array2::zeros((0, 0)),
1104            max_iter: 1000,
1105            tol: 1e-4,
1106            fit_intercept: true,
1107            normalize: true,
1108            state: PhantomData,
1109            coefficients_: None,
1110            selected_features_: None,
1111            _intercept_: None,
1112        }
1113    }
1114
1115    /// alpha
1116    pub fn alpha(mut self, alpha: f64) -> Self {
1117        self.alpha = alpha;
1118        self
1119    }
1120
1121    /// adjacency_matrix
1122    pub fn adjacency_matrix(mut self, matrix: Array2<Float>) -> Self {
1123        self.adjacency_matrix = matrix;
1124        self
1125    }
1126
1127    /// max_iter
1128    pub fn max_iter(mut self, max_iter: usize) -> Self {
1129        self.max_iter = max_iter;
1130        self
1131    }
1132
1133    /// tol
1134    pub fn tol(mut self, tol: f64) -> Self {
1135        self.tol = tol;
1136        self
1137    }
1138
1139    /// fit_intercept
1140    pub fn fit_intercept(mut self, fit_intercept: bool) -> Self {
1141        self.fit_intercept = fit_intercept;
1142        self
1143    }
1144
1145    /// normalize
1146    pub fn normalize(mut self, normalize: bool) -> Self {
1147        self.normalize = normalize;
1148        self
1149    }
1150}
1151
1152impl Default for GraphStructuredSparsitySelector<Untrained> {
1153    fn default() -> Self {
1154        Self::new()
1155    }
1156}
1157
1158impl Estimator for GraphStructuredSparsitySelector<Untrained> {
1159    type Config = ();
1160    type Error = SklearsError;
1161    type Float = f64;
1162
1163    fn config(&self) -> &Self::Config {
1164        &()
1165    }
1166}
1167
1168impl Fit<Array2<Float>, Array1<Float>> for GraphStructuredSparsitySelector<Untrained> {
1169    type Fitted = GraphStructuredSparsitySelector<Trained>;
1170
1171    fn fit(self, x: &Array2<Float>, y: &Array1<Float>) -> SklResult<Self::Fitted> {
1172        let (n_samples, n_features) = x.dim();
1173
1174        if n_samples != y.len() {
1175            return Err(SklearsError::InvalidInput(
1176                "X and y must have the same number of samples".to_string(),
1177            ));
1178        }
1179
1180        if self.adjacency_matrix.nrows() != n_features
1181            || self.adjacency_matrix.ncols() != n_features
1182        {
1183            return Err(SklearsError::InvalidInput(
1184                "Adjacency matrix dimensions must match number of features".to_string(),
1185            ));
1186        }
1187
1188        // Compute graph Laplacian
1189        let laplacian = self.compute_graph_laplacian()?;
1190
1191        // Fit using proximal gradient descent with graph penalty
1192        let (coefficients, y_mean) = self.fit_graph_sparsity(x, y, &laplacian)?;
1193
1194        // Determine selected features
1195        let selected_features: Vec<usize> = coefficients
1196            .indexed_iter()
1197            .filter_map(|(i, &coeff)| if coeff.abs() > 1e-10 { Some(i) } else { None })
1198            .collect();
1199
1200        let intercept = if self.fit_intercept {
1201            Some(y_mean)
1202        } else {
1203            None
1204        };
1205
1206        Ok(GraphStructuredSparsitySelector {
1207            alpha: self.alpha,
1208            adjacency_matrix: self.adjacency_matrix,
1209            max_iter: self.max_iter,
1210            tol: self.tol,
1211            fit_intercept: self.fit_intercept,
1212            normalize: self.normalize,
1213            state: PhantomData,
1214            coefficients_: Some(coefficients),
1215            selected_features_: Some(selected_features),
1216            _intercept_: intercept,
1217        })
1218    }
1219}
1220
1221impl GraphStructuredSparsitySelector<Untrained> {
1222    fn compute_graph_laplacian(&self) -> SklResult<Array2<Float>> {
1223        let n_features = self.adjacency_matrix.nrows();
1224        let mut laplacian = Array2::zeros((n_features, n_features));
1225
1226        // Compute degree matrix
1227        let mut degrees = Array1::zeros(n_features);
1228        for i in 0..n_features {
1229            degrees[i] = self.adjacency_matrix.row(i).sum();
1230        }
1231
1232        // L = D - A
1233        for i in 0..n_features {
1234            laplacian[[i, i]] = degrees[i];
1235            for j in 0..n_features {
1236                if i != j {
1237                    laplacian[[i, j]] = -self.adjacency_matrix[[i, j]];
1238                }
1239            }
1240        }
1241
1242        Ok(laplacian)
1243    }
1244
1245    fn fit_graph_sparsity(
1246        &self,
1247        x: &Array2<Float>,
1248        y: &Array1<Float>,
1249        laplacian: &Array2<Float>,
1250    ) -> SklResult<(Array1<Float>, f64)> {
1251        let (n_samples, n_features) = x.dim();
1252        let mut coefficients = Array1::zeros(n_features);
1253        let step_size = 0.01;
1254
1255        // Center target
1256        let y_mean = if self.fit_intercept {
1257            y.mean().unwrap_or(0.0)
1258        } else {
1259            0.0
1260        };
1261
1262        let mut y_centered = y.clone();
1263        if self.fit_intercept {
1264            for i in 0..n_samples {
1265                y_centered[i] -= y_mean;
1266            }
1267        }
1268
1269        for _iter in 0..self.max_iter {
1270            let coefficients_old = coefficients.clone();
1271
1272            // Compute gradient of loss function
1273            let mut gradient = Array1::<f64>::zeros(n_features);
1274            for i in 0..n_samples {
1275                let residual = y_centered[i] - x.row(i).dot(&coefficients);
1276                for j in 0..n_features {
1277                    gradient[j] -= x[[i, j]] * residual / n_samples as f64;
1278                }
1279            }
1280
1281            // Add graph penalty gradient: alpha * L * beta
1282            let graph_penalty_grad = laplacian.dot(&coefficients) * self.alpha;
1283            gradient = gradient + graph_penalty_grad;
1284
1285            // Gradient descent step
1286            for j in 0..n_features {
1287                coefficients[j] -= step_size * gradient[j];
1288            }
1289
1290            // Apply soft thresholding for sparsity
1291            for coeff in coefficients.iter_mut() {
1292                *coeff = soft_threshold(*coeff, step_size * self.alpha * 0.1);
1293            }
1294
1295            // Check convergence
1296            let diff = (&coefficients - &coefficients_old).mapv(|x| x.abs()).sum();
1297            if diff < self.tol {
1298                break;
1299            }
1300        }
1301
1302        Ok((coefficients, y_mean))
1303    }
1304}
1305
1306impl Transform<Array2<Float>> for GraphStructuredSparsitySelector<Trained> {
1307    fn transform(&self, x: &Array2<Float>) -> SklResult<Array2<Float>> {
1308        let selected_features = self
1309            .selected_features_
1310            .as_ref()
1311            .expect("operation should succeed");
1312        if selected_features.is_empty() {
1313            return Err(SklearsError::InvalidInput(
1314                "No features were selected".to_string(),
1315            ));
1316        }
1317
1318        let selected_indices: Vec<usize> = selected_features.to_vec();
1319        Ok(x.select(Axis(1), &selected_indices))
1320    }
1321}
1322
1323impl SelectorMixin for GraphStructuredSparsitySelector<Trained> {
1324    fn get_support(&self) -> SklResult<Array1<bool>> {
1325        let selected_features = self
1326            .selected_features_
1327            .as_ref()
1328            .expect("operation should succeed");
1329        let n_features = self
1330            .coefficients_
1331            .as_ref()
1332            .expect("operation should succeed")
1333            .len();
1334        let mut support = Array1::from_elem(n_features, false);
1335        for &idx in selected_features {
1336            if idx < n_features {
1337                support[idx] = true;
1338            }
1339        }
1340        Ok(support)
1341    }
1342
1343    fn transform_features(&self, indices: &[usize]) -> SklResult<Vec<usize>> {
1344        let selected_features = self
1345            .selected_features_
1346            .as_ref()
1347            .expect("operation should succeed");
1348        Ok(indices
1349            .iter()
1350            .filter_map(|&idx| selected_features.iter().position(|&f| f == idx))
1351            .collect())
1352    }
1353}
1354
1355/// Overlapping group sparsity selector
1356///
1357/// Implements structured sparsity for overlapping groups of features.
1358/// Features can belong to multiple groups, and the penalty encourages
1359/// selecting entire groups while handling the overlap structure.
1360#[derive(Debug, Clone)]
1361pub struct OverlappingGroupSparsitySelector<State = Untrained> {
1362    /// Regularization parameter
1363    alpha: f64,
1364    /// Overlapping groups of feature indices
1365    overlapping_groups: Vec<Vec<usize>>,
1366    /// Group names for interpretability
1367    group_names: Option<Vec<String>>,
1368    /// Maximum number of iterations
1369    max_iter: usize,
1370    /// Convergence tolerance
1371    tol: f64,
1372    /// Whether to fit intercept
1373    fit_intercept: bool,
1374    /// Whether to normalize features
1375    normalize: bool,
1376    state: PhantomData<State>,
1377    // Trained state
1378    coefficients_: Option<Array1<Float>>,
1379    selected_features_: Option<Vec<usize>>,
1380    _selected_groups_: Option<Vec<usize>>,
1381    _intercept_: Option<f64>,
1382}
1383
1384impl OverlappingGroupSparsitySelector<Untrained> {
1385    /// new
1386    pub fn new() -> Self {
1387        Self {
1388            alpha: 1.0,
1389            overlapping_groups: Vec::new(),
1390            group_names: None,
1391            max_iter: 1000,
1392            tol: 1e-4,
1393            fit_intercept: true,
1394            normalize: true,
1395            state: PhantomData,
1396            coefficients_: None,
1397            selected_features_: None,
1398            _selected_groups_: None,
1399            _intercept_: None,
1400        }
1401    }
1402
1403    /// alpha
1404    pub fn alpha(mut self, alpha: f64) -> Self {
1405        self.alpha = alpha;
1406        self
1407    }
1408
1409    /// overlapping_groups
1410    pub fn overlapping_groups(mut self, groups: Vec<Vec<usize>>) -> Self {
1411        self.overlapping_groups = groups;
1412        self
1413    }
1414
1415    /// group_names
1416    pub fn group_names(mut self, names: Vec<String>) -> Self {
1417        self.group_names = Some(names);
1418        self
1419    }
1420
1421    /// max_iter
1422    pub fn max_iter(mut self, max_iter: usize) -> Self {
1423        self.max_iter = max_iter;
1424        self
1425    }
1426
1427    /// tol
1428    pub fn tol(mut self, tol: f64) -> Self {
1429        self.tol = tol;
1430        self
1431    }
1432
1433    /// fit_intercept
1434    pub fn fit_intercept(mut self, fit_intercept: bool) -> Self {
1435        self.fit_intercept = fit_intercept;
1436        self
1437    }
1438
1439    /// normalize
1440    pub fn normalize(mut self, normalize: bool) -> Self {
1441        self.normalize = normalize;
1442        self
1443    }
1444}
1445
1446impl Default for OverlappingGroupSparsitySelector<Untrained> {
1447    fn default() -> Self {
1448        Self::new()
1449    }
1450}
1451
1452impl Estimator for OverlappingGroupSparsitySelector<Untrained> {
1453    type Config = ();
1454    type Error = SklearsError;
1455    type Float = f64;
1456
1457    fn config(&self) -> &Self::Config {
1458        &()
1459    }
1460}
1461
1462impl Fit<Array2<Float>, Array1<Float>> for OverlappingGroupSparsitySelector<Untrained> {
1463    type Fitted = OverlappingGroupSparsitySelector<Trained>;
1464
1465    fn fit(self, x: &Array2<Float>, y: &Array1<Float>) -> SklResult<Self::Fitted> {
1466        let (n_samples, n_features) = x.dim();
1467
1468        if n_samples != y.len() {
1469            return Err(SklearsError::InvalidInput(
1470                "X and y must have the same number of samples".to_string(),
1471            ));
1472        }
1473
1474        if self.overlapping_groups.is_empty() {
1475            return Err(SklearsError::InvalidInput(
1476                "Overlapping groups must be specified".to_string(),
1477            ));
1478        }
1479
1480        // Validate overlapping groups
1481        for group in &self.overlapping_groups {
1482            for &feature in group {
1483                if feature >= n_features {
1484                    return Err(SklearsError::InvalidInput(
1485                        "Feature index out of bounds in groups".to_string(),
1486                    ));
1487                }
1488            }
1489        }
1490
1491        // Fit using alternating direction method of multipliers (ADMM)
1492        let (coefficients, selected_groups, y_mean) = self.fit_overlapping_groups(x, y)?;
1493
1494        // Determine selected features from selected groups
1495        let mut selected_features = HashSet::new();
1496        for &group_idx in &selected_groups {
1497            for &feature in &self.overlapping_groups[group_idx] {
1498                selected_features.insert(feature);
1499            }
1500        }
1501        let selected_features: Vec<usize> = selected_features.into_iter().collect();
1502
1503        let intercept = if self.fit_intercept {
1504            Some(y_mean)
1505        } else {
1506            None
1507        };
1508
1509        Ok(OverlappingGroupSparsitySelector {
1510            alpha: self.alpha,
1511            overlapping_groups: self.overlapping_groups,
1512            group_names: self.group_names,
1513            max_iter: self.max_iter,
1514            tol: self.tol,
1515            fit_intercept: self.fit_intercept,
1516            normalize: self.normalize,
1517            state: PhantomData,
1518            coefficients_: Some(coefficients),
1519            selected_features_: Some(selected_features),
1520            _selected_groups_: Some(selected_groups),
1521            _intercept_: intercept,
1522        })
1523    }
1524}
1525
1526impl OverlappingGroupSparsitySelector<Untrained> {
1527    fn fit_overlapping_groups(
1528        &self,
1529        x: &Array2<Float>,
1530        y: &Array1<Float>,
1531    ) -> SklResult<(Array1<Float>, Vec<usize>, f64)> {
1532        let (n_samples, n_features) = x.dim();
1533        let n_groups = self.overlapping_groups.len();
1534        let mut coefficients = Array1::zeros(n_features);
1535        let step_size = 0.01;
1536
1537        // Center target
1538        let y_mean = if self.fit_intercept {
1539            y.mean().unwrap_or(0.0)
1540        } else {
1541            0.0
1542        };
1543
1544        let mut y_centered = y.clone();
1545        if self.fit_intercept {
1546            for i in 0..n_samples {
1547                y_centered[i] -= y_mean;
1548            }
1549        }
1550
1551        // Create overlap matrix to handle overlapping groups
1552        let mut overlap_matrix = Array2::zeros((n_features, n_groups));
1553        for (g, group) in self.overlapping_groups.iter().enumerate() {
1554            for &feature in group {
1555                overlap_matrix[[feature, g]] = 1.0;
1556            }
1557        }
1558
1559        // Compute overlap weights (how many groups each feature belongs to)
1560        let mut overlap_weights = Array1::zeros(n_features);
1561        for i in 0..n_features {
1562            overlap_weights[i] = overlap_matrix.row(i).sum();
1563            if overlap_weights[i] == 0.0 {
1564                overlap_weights[i] = 1.0; // Avoid division by zero
1565            }
1566        }
1567
1568        for _iter in 0..self.max_iter {
1569            let coefficients_old = coefficients.clone();
1570
1571            // Compute gradient
1572            let mut gradient = Array1::<f64>::zeros(n_features);
1573            for i in 0..n_samples {
1574                let residual = y_centered[i] - x.row(i).dot(&coefficients);
1575                for j in 0..n_features {
1576                    gradient[j] -= x[[i, j]] * residual / n_samples as f64;
1577                }
1578            }
1579
1580            // Gradient descent step
1581            for j in 0..n_features {
1582                coefficients[j] -= step_size * gradient[j];
1583            }
1584
1585            // Apply overlapping group soft thresholding
1586            for group in self.overlapping_groups.iter() {
1587                // Compute group norm
1588                let group_norm = group
1589                    .iter()
1590                    .map(|&f| coefficients[f] * coefficients[f])
1591                    .sum::<f64>()
1592                    .sqrt();
1593
1594                // Apply group soft thresholding with overlap adjustment
1595                let group_size = group.len() as f64;
1596                let threshold = self.alpha * group_size.sqrt();
1597
1598                if group_norm > threshold {
1599                    let shrinkage_factor = 1.0 - threshold / group_norm;
1600                    for &feature in group {
1601                        // Adjust for overlap: features in multiple groups get averaged adjustment
1602                        let adjustment = shrinkage_factor / overlap_weights[feature];
1603                        coefficients[feature] = coefficients[feature] * (1.0 - adjustment)
1604                            + coefficients[feature] * shrinkage_factor * adjustment;
1605                    }
1606                }
1607            }
1608
1609            // Check convergence
1610            let diff = (&coefficients - &coefficients_old).mapv(|x| x.abs()).sum();
1611            if diff < self.tol {
1612                break;
1613            }
1614        }
1615
1616        // Determine selected groups
1617        let selected_groups: Vec<usize> = (0..n_groups)
1618            .filter(|&g| {
1619                let group_norm = self.overlapping_groups[g]
1620                    .iter()
1621                    .map(|&f| coefficients[f] * coefficients[f])
1622                    .sum::<f64>()
1623                    .sqrt();
1624                group_norm > 1e-10
1625            })
1626            .collect();
1627
1628        Ok((coefficients, selected_groups, y_mean))
1629    }
1630}
1631
1632impl Transform<Array2<Float>> for OverlappingGroupSparsitySelector<Trained> {
1633    fn transform(&self, x: &Array2<Float>) -> SklResult<Array2<Float>> {
1634        let selected_features = self
1635            .selected_features_
1636            .as_ref()
1637            .expect("operation should succeed");
1638        if selected_features.is_empty() {
1639            return Err(SklearsError::InvalidInput(
1640                "No features were selected".to_string(),
1641            ));
1642        }
1643
1644        let selected_indices: Vec<usize> = selected_features.to_vec();
1645        Ok(x.select(Axis(1), &selected_indices))
1646    }
1647}
1648
1649impl SelectorMixin for OverlappingGroupSparsitySelector<Trained> {
1650    fn get_support(&self) -> SklResult<Array1<bool>> {
1651        let selected_features = self
1652            .selected_features_
1653            .as_ref()
1654            .expect("operation should succeed");
1655        let n_features = self
1656            .coefficients_
1657            .as_ref()
1658            .expect("operation should succeed")
1659            .len();
1660        let mut support = Array1::from_elem(n_features, false);
1661        for &idx in selected_features {
1662            if idx < n_features {
1663                support[idx] = true;
1664            }
1665        }
1666        Ok(support)
1667    }
1668
1669    fn transform_features(&self, indices: &[usize]) -> SklResult<Vec<usize>> {
1670        let selected_features = self
1671            .selected_features_
1672            .as_ref()
1673            .expect("operation should succeed");
1674        Ok(indices
1675            .iter()
1676            .filter_map(|&idx| selected_features.iter().position(|&f| f == idx))
1677            .collect())
1678    }
1679}
1680
1681#[allow(non_snake_case)]
1682#[cfg(test)]
1683mod tests {
1684    use super::*;
1685    use scirs2_core::ndarray::array;
1686
1687    #[test]
1688    fn test_group_lasso_basic() {
1689        let x = array![
1690            [1.0, 0.0, 2.0, 0.0],
1691            [0.0, 1.0, 0.0, 3.0],
1692            [2.0, 0.0, 4.0, 0.0],
1693            [0.0, 2.0, 0.0, 6.0],
1694            [3.0, 0.0, 6.0, 0.0],
1695        ];
1696        let y = array![1.0, 2.0, 2.0, 4.0, 3.0];
1697
1698        // Define groups: [0, 2] and [1, 3]
1699        let groups = vec![vec![0, 2], vec![1, 3]];
1700
1701        let selector = GroupLassoSelector::new().alpha(0.1).feature_groups(groups);
1702
1703        let fitted_selector = selector.fit(&x, &y).expect("operation should succeed");
1704
1705        assert!(fitted_selector.coefficients().is_some());
1706        assert!(fitted_selector.selected_groups().is_some());
1707        assert!(!fitted_selector
1708            .selected_groups()
1709            .expect("operation should succeed")
1710            .is_empty());
1711    }
1712
1713    #[test]
1714    fn test_group_lasso_high_regularization() {
1715        let x = array![[1.0, 2.0, 3.0], [2.0, 3.0, 4.0], [3.0, 4.0, 5.0],];
1716        let y = array![1.0, 2.0, 3.0];
1717
1718        let groups = vec![vec![0], vec![1], vec![2]];
1719        let n_groups = groups.len();
1720
1721        let selector = GroupLassoSelector::new()
1722            .alpha(10.0) // High regularization
1723            .feature_groups(groups);
1724
1725        let fitted_selector = selector.fit(&x, &y).expect("operation should succeed");
1726
1727        // With high regularization, few or no groups should be selected
1728        let selected_groups = fitted_selector
1729            .selected_groups()
1730            .expect("operation should succeed");
1731        assert!(selected_groups.len() <= n_groups);
1732    }
1733
1734    #[test]
1735    fn test_sparse_group_lasso() {
1736        let x = array![
1737            [1.0, 0.0, 2.0, 0.0],
1738            [0.0, 1.0, 0.0, 3.0],
1739            [2.0, 0.0, 4.0, 0.0],
1740            [0.0, 2.0, 0.0, 6.0],
1741        ];
1742        let y = array![1.0, 2.0, 2.0, 4.0];
1743
1744        let groups = vec![vec![0, 1], vec![2, 3]];
1745
1746        let selector = SparseGroupLassoSelector::new()
1747            .alpha_group(0.1)
1748            .alpha_element(0.05)
1749            .feature_groups(groups);
1750
1751        let fitted_selector = selector.fit(&x, &y).expect("operation should succeed");
1752
1753        assert!(fitted_selector.coefficients_.is_some());
1754        assert!(fitted_selector.selected_features_.is_some());
1755    }
1756
1757    #[test]
1758    fn test_group_lasso_transform() {
1759        let x = array![
1760            [1.0, 2.0, 3.0, 4.0],
1761            [2.0, 3.0, 4.0, 5.0],
1762            [3.0, 4.0, 5.0, 6.0],
1763        ];
1764        let y = array![1.0, 2.0, 3.0];
1765
1766        let groups = vec![vec![0, 1], vec![2, 3]];
1767
1768        let selector = GroupLassoSelector::new().alpha(0.1).feature_groups(groups);
1769
1770        let fitted_selector = selector.fit(&x, &y).expect("operation should succeed");
1771        let transformed = fitted_selector
1772            .transform(&x)
1773            .expect("operation should succeed");
1774
1775        assert!(transformed.ncols() <= x.ncols());
1776        assert_eq!(transformed.nrows(), x.nrows());
1777    }
1778
1779    #[test]
1780    fn test_soft_threshold() {
1781        assert_eq!(soft_threshold(2.0, 1.0), 1.0);
1782        assert_eq!(soft_threshold(-2.0, 1.0), -1.0);
1783        assert_eq!(soft_threshold(0.5, 1.0), 0.0);
1784        assert_eq!(soft_threshold(-0.5, 1.0), 0.0);
1785    }
1786
1787    #[test]
1788    fn test_hierarchical_structured_sparsity_basic() {
1789        let x = array![
1790            [1.0, 0.0, 2.0, 0.0],
1791            [0.0, 1.0, 0.0, 3.0],
1792            [2.0, 0.0, 4.0, 0.0],
1793            [0.0, 2.0, 0.0, 6.0],
1794        ];
1795        let y = array![1.0, 2.0, 2.0, 4.0];
1796
1797        // Define hierarchy: node 0 has children [1, 2], features map as 0->0, 1->1, 2->0, 3->2
1798        let mut hierarchy = HashMap::new();
1799        hierarchy.insert(0, vec![1, 2]);
1800
1801        let mut feature_to_node = HashMap::new();
1802        feature_to_node.insert(0, 0);
1803        feature_to_node.insert(1, 1);
1804        feature_to_node.insert(2, 0);
1805        feature_to_node.insert(3, 2);
1806
1807        let selector = HierarchicalStructuredSparsitySelector::new()
1808            .alpha(0.1)
1809            .hierarchy(hierarchy)
1810            .feature_to_node(feature_to_node);
1811
1812        let fitted_selector = selector.fit(&x, &y).expect("operation should succeed");
1813
1814        assert!(fitted_selector.coefficients_.is_some());
1815        assert!(fitted_selector.selected_features_.is_some());
1816    }
1817
1818    #[test]
1819    fn test_graph_structured_sparsity_basic() {
1820        let x = array![[1.0, 2.0, 3.0], [2.0, 3.0, 4.0], [3.0, 4.0, 5.0],];
1821        let y = array![1.0, 2.0, 3.0];
1822
1823        // Define adjacency matrix (features 0-1 connected, 1-2 connected)
1824        let adjacency_matrix = array![[0.0, 1.0, 0.0], [1.0, 0.0, 1.0], [0.0, 1.0, 0.0],];
1825
1826        let selector = GraphStructuredSparsitySelector::new()
1827            .alpha(0.1)
1828            .adjacency_matrix(adjacency_matrix);
1829
1830        let fitted_selector = selector.fit(&x, &y).expect("operation should succeed");
1831
1832        assert!(fitted_selector.coefficients_.is_some());
1833        assert!(fitted_selector.selected_features_.is_some());
1834    }
1835
1836    #[test]
1837    fn test_overlapping_group_sparsity_basic() {
1838        let x = array![
1839            [1.0, 0.0, 2.0, 0.0],
1840            [0.0, 1.0, 0.0, 3.0],
1841            [2.0, 0.0, 4.0, 0.0],
1842            [0.0, 2.0, 0.0, 6.0],
1843        ];
1844        let y = array![1.0, 2.0, 2.0, 4.0];
1845
1846        // Define overlapping groups: [0, 1], [1, 2], [2, 3]
1847        let overlapping_groups = vec![vec![0, 1], vec![1, 2], vec![2, 3]];
1848
1849        let selector = OverlappingGroupSparsitySelector::new()
1850            .alpha(0.1)
1851            .overlapping_groups(overlapping_groups);
1852
1853        let fitted_selector = selector.fit(&x, &y).expect("operation should succeed");
1854
1855        assert!(fitted_selector.coefficients_.is_some());
1856        assert!(fitted_selector.selected_features_.is_some());
1857        assert!(fitted_selector._selected_groups_.is_some());
1858    }
1859
1860    #[test]
1861    fn test_hierarchical_transform() {
1862        let x = array![
1863            [1.0, 2.0, 3.0, 4.0],
1864            [2.0, 3.0, 4.0, 5.0],
1865            [3.0, 4.0, 5.0, 6.0],
1866        ];
1867        let y = array![1.0, 2.0, 3.0];
1868
1869        let mut hierarchy = HashMap::new();
1870        hierarchy.insert(0, vec![1]);
1871
1872        let mut feature_to_node = HashMap::new();
1873        feature_to_node.insert(0, 0);
1874        feature_to_node.insert(1, 0);
1875        feature_to_node.insert(2, 1);
1876        feature_to_node.insert(3, 1);
1877
1878        let selector = HierarchicalStructuredSparsitySelector::new()
1879            .alpha(0.01) // Lower alpha to ensure some features are selected
1880            .hierarchy(hierarchy)
1881            .feature_to_node(feature_to_node);
1882
1883        let fitted_selector = selector.fit(&x, &y).expect("operation should succeed");
1884
1885        // Check if any features were selected before transforming
1886        if !fitted_selector
1887            .selected_features_
1888            .as_ref()
1889            .expect("operation should succeed")
1890            .is_empty()
1891        {
1892            let transformed = fitted_selector
1893                .transform(&x)
1894                .expect("operation should succeed");
1895            assert!(transformed.ncols() <= x.ncols());
1896            assert_eq!(transformed.nrows(), x.nrows());
1897        } else {
1898            // If no features selected, this is acceptable for this algorithm
1899            assert!(fitted_selector
1900                .selected_features_
1901                .as_ref()
1902                .expect("operation should succeed")
1903                .is_empty());
1904        }
1905    }
1906
1907    #[test]
1908    fn test_overlapping_group_with_group_names() {
1909        let x = array![[1.0, 2.0, 3.0], [2.0, 3.0, 4.0], [3.0, 4.0, 5.0],];
1910        let y = array![1.0, 2.0, 3.0];
1911
1912        let overlapping_groups = vec![vec![0], vec![1], vec![2]];
1913        let group_names = vec![
1914            "Group A".to_string(),
1915            "Group B".to_string(),
1916            "Group C".to_string(),
1917        ];
1918
1919        let selector = OverlappingGroupSparsitySelector::new()
1920            .alpha(0.1)
1921            .overlapping_groups(overlapping_groups)
1922            .group_names(group_names);
1923
1924        let fitted_selector = selector.fit(&x, &y).expect("operation should succeed");
1925
1926        assert!(fitted_selector.coefficients_.is_some());
1927        assert!(fitted_selector.group_names.is_some());
1928    }
1929}
sklears_feature_selection/group_selection.rs

sklears_feature_selection/
group_selection.rs
