sklears_multioutput/
svm.rs

1//! Support Vector Machine algorithms for multi-output learning
2
3// Use SciRS2-Core for arrays and random number generation (SciRS2 Policy)
4use scirs2_core::ndarray::{s, Array1, Array2, ArrayView1, ArrayView2, Axis};
5use sklears_core::{
6    error::{Result as SklResult, SklearsError},
7    traits::{Estimator, Fit, Predict, Untrained},
8    types::Float,
9};
10
11/// MLTSVM (Multi-Label Twin SVM)
12///
13/// MLTSVM is a multi-label classification method that extends Twin SVM to handle
14/// multiple labels. Twin SVM finds two non-parallel hyperplanes for binary classification,
15/// which often leads to faster training than standard SVM. MLTSVM applies this approach
16/// to each label independently in a binary relevance fashion.
17///
18/// # Examples
19///
20/// ```
21/// use sklears_core::traits::{Predict, Fit};
22/// use sklears_multioutput::MLTSVM;
23/// // Use SciRS2-Core for arrays and random number generation (SciRS2 Policy)
24/// use scirs2_core::ndarray::array;
25///
26/// let X = array![[1.0, 2.0], [2.0, 3.0], [3.0, 1.0], [4.0, 4.0]];
27/// let y = array![[1, 0], [0, 1], [1, 1], [0, 0]]; // Multi-label binary
28///
29/// let mltsvm = MLTSVM::new().c1(1.0).c2(1.0);
30/// let trained_mltsvm = mltsvm.fit(&X.view(), &y).unwrap();
31/// let predictions = trained_mltsvm.predict(&X.view()).unwrap();
32/// ```
33#[derive(Debug, Clone)]
34pub struct MLTSVM<S = Untrained> {
35    state: S,
36    c1: Float,       // Regularization parameter for first hyperplane
37    c2: Float,       // Regularization parameter for second hyperplane
38    epsilon: Float,  // Tolerance for convergence
39    max_iter: usize, // Maximum iterations
40}
41
42/// Trained state for MLTSVM
43#[derive(Debug, Clone)]
44pub struct MLTSVMTrained {
45    models: Vec<TwinSVMModel>, // One model per label
46    n_labels: usize,
47    feature_means: Array1<Float>,
48    feature_stds: Array1<Float>,
49}
50
51/// Twin SVM model for a single label
52#[derive(Debug, Clone)]
53pub struct TwinSVMModel {
54    w1: Array1<Float>, // Weight vector for positive hyperplane
55    b1: Float,         // Bias for positive hyperplane
56    w2: Array1<Float>, // Weight vector for negative hyperplane
57    b2: Float,         // Bias for negative hyperplane
58}
59
60impl MLTSVM<Untrained> {
61    /// Create a new MLTSVM instance
62    pub fn new() -> Self {
63        Self {
64            state: Untrained,
65            c1: 1.0,
66            c2: 1.0,
67            epsilon: 1e-3,
68            max_iter: 1000,
69        }
70    }
71
72    /// Set C1 parameter
73    pub fn c1(mut self, c1: Float) -> Self {
74        self.c1 = c1;
75        self
76    }
77
78    /// Set C2 parameter
79    pub fn c2(mut self, c2: Float) -> Self {
80        self.c2 = c2;
81        self
82    }
83
84    /// Set epsilon parameter
85    pub fn epsilon(mut self, epsilon: Float) -> Self {
86        self.epsilon = epsilon;
87        self
88    }
89
90    /// Set maximum iterations
91    pub fn max_iter(mut self, max_iter: usize) -> Self {
92        self.max_iter = max_iter;
93        self
94    }
95}
96
97impl Default for MLTSVM<Untrained> {
98    fn default() -> Self {
99        Self::new()
100    }
101}
102
103impl Estimator for MLTSVM<Untrained> {
104    type Config = ();
105    type Error = SklearsError;
106    type Float = Float;
107
108    fn config(&self) -> &Self::Config {
109        &()
110    }
111}
112
113impl Fit<ArrayView2<'_, Float>, Array2<i32>> for MLTSVM<Untrained> {
114    type Fitted = MLTSVM<MLTSVMTrained>;
115
116    fn fit(self, x: &ArrayView2<'_, Float>, y: &Array2<i32>) -> SklResult<Self::Fitted> {
117        let (n_samples, _n_features) = x.dim();
118        let (y_samples, n_labels) = y.dim();
119
120        if n_samples != y_samples {
121            return Err(SklearsError::InvalidInput(
122                "Number of samples in X and y must match".to_string(),
123            ));
124        }
125
126        if n_samples < 2 {
127            return Err(SklearsError::InvalidInput(
128                "Need at least 2 samples for SVM training".to_string(),
129            ));
130        }
131
132        // Validate that all labels are binary (0 or 1)
133        for sample_idx in 0..y_samples {
134            for label_idx in 0..n_labels {
135                let value = y[[sample_idx, label_idx]];
136                if value != 0 && value != 1 {
137                    return Err(SklearsError::InvalidInput(format!(
138                        "All label values must be 0 or 1, found: {}",
139                        value
140                    )));
141                }
142            }
143        }
144
145        // Compute feature statistics for normalization
146        let feature_means = x
147            .mean_axis(Axis(0))
148            .expect("array should have elements for mean computation");
149        let feature_stds = x
150            .mapv(|val| val * val)
151            .mean_axis(Axis(0))
152            .expect("array should have elements for mean computation")
153            - &feature_means.mapv(|mean| mean * mean);
154        let feature_stds = feature_stds.mapv(|var| (var.max(1e-10)).sqrt());
155
156        // Train Twin SVM for each label
157        let mut models = Vec::new();
158        for label_idx in 0..n_labels {
159            let y_label = y.column(label_idx);
160            let model = self.train_twin_svm(x, &y_label, &feature_means, &feature_stds)?;
161            models.push(model);
162        }
163
164        Ok(MLTSVM {
165            state: MLTSVMTrained {
166                models,
167                n_labels,
168                feature_means,
169                feature_stds,
170            },
171            c1: self.c1,
172            c2: self.c2,
173            epsilon: self.epsilon,
174            max_iter: self.max_iter,
175        })
176    }
177}
178
179impl MLTSVM<Untrained> {
180    fn train_twin_svm(
181        &self,
182        x: &ArrayView2<'_, Float>,
183        y: &ArrayView1<'_, i32>,
184        feature_means: &Array1<Float>,
185        feature_stds: &Array1<Float>,
186    ) -> SklResult<TwinSVMModel> {
187        let (n_samples, n_features) = x.dim();
188
189        // Normalize features
190        let mut x_normalized = x.to_owned();
191        for mut row in x_normalized.rows_mut().into_iter() {
192            row -= feature_means;
193            row /= feature_stds;
194        }
195
196        // Separate positive and negative samples
197        let mut pos_samples = Vec::new();
198        let mut neg_samples = Vec::new();
199
200        for i in 0..n_samples {
201            if y[i] == 1 {
202                pos_samples.push(x_normalized.row(i).to_owned());
203            } else {
204                neg_samples.push(x_normalized.row(i).to_owned());
205            }
206        }
207
208        if pos_samples.is_empty() || neg_samples.is_empty() {
209            return Err(SklearsError::InvalidInput(
210                "Need both positive and negative samples for Twin SVM".to_string(),
211            ));
212        }
213
214        // Convert to matrices
215        let pos_matrix = Array2::from_shape_vec(
216            (pos_samples.len(), n_features),
217            pos_samples.into_iter().flatten().collect(),
218        )
219        .map_err(|_| SklearsError::InvalidInput("Failed to create positive matrix".to_string()))?;
220
221        let neg_matrix = Array2::from_shape_vec(
222            (neg_samples.len(), n_features),
223            neg_samples.into_iter().flatten().collect(),
224        )
225        .map_err(|_| SklearsError::InvalidInput("Failed to create negative matrix".to_string()))?;
226
227        // Train Twin SVM hyperplanes
228        let (w1, b1) = self.solve_twin_svm_problem(&pos_matrix, &neg_matrix, self.c1)?;
229        let (w2, b2) = self.solve_twin_svm_problem(&neg_matrix, &pos_matrix, self.c2)?;
230
231        Ok(TwinSVMModel { w1, b1, w2, b2 })
232    }
233
234    fn solve_twin_svm_problem(
235        &self,
236        target_matrix: &Array2<Float>,
237        other_matrix: &Array2<Float>,
238        c: Float,
239    ) -> SklResult<(Array1<Float>, Float)> {
240        let n_target = target_matrix.nrows();
241        let n_other = other_matrix.nrows();
242        let n_features = target_matrix.ncols();
243
244        // Initialize weights
245        let mut w = Array1::<Float>::zeros(n_features + 1); // Include bias
246
247        // Simple gradient descent solution
248        let learning_rate = 0.01;
249
250        for _iter in 0..self.max_iter {
251            let mut gradient = Array1::<Float>::zeros(n_features + 1);
252
253            // Compute gradient
254            for i in 0..n_target {
255                let x_aug = {
256                    let mut x = Array1::ones(n_features + 1);
257                    x.slice_mut(s![..n_features]).assign(&target_matrix.row(i));
258                    x
259                };
260                let loss = x_aug.dot(&w);
261                gradient += &(x_aug * loss);
262            }
263
264            for i in 0..n_other {
265                let x_aug = {
266                    let mut x = Array1::ones(n_features + 1);
267                    x.slice_mut(s![..n_features]).assign(&other_matrix.row(i));
268                    x
269                };
270                let margin = 1.0 - x_aug.dot(&w);
271                if margin > 0.0 {
272                    gradient -= &(x_aug * c);
273                }
274            }
275
276            // Check convergence before updating weights
277            let gradient_norm = gradient.mapv(|x| x.abs()).sum();
278
279            // Update weights
280            w -= &(gradient * learning_rate);
281
282            if gradient_norm < self.epsilon {
283                break;
284            }
285        }
286
287        let weights = w.slice(s![..n_features]).to_owned();
288        let bias = w[n_features];
289
290        Ok((weights, bias))
291    }
292}
293
294impl Predict<ArrayView2<'_, Float>, Array2<i32>> for MLTSVM<MLTSVMTrained> {
295    fn predict(&self, x: &ArrayView2<'_, Float>) -> SklResult<Array2<i32>> {
296        let (n_samples, n_features) = x.dim();
297        let expected_features = self.state.feature_means.len();
298
299        if n_features != expected_features {
300            return Err(SklearsError::InvalidInput(format!(
301                "Number of features in X ({}) does not match training data ({})",
302                n_features, expected_features
303            )));
304        }
305
306        let mut predictions = Array2::<i32>::zeros((n_samples, self.state.n_labels));
307
308        // Normalize features
309        let mut x_normalized = x.to_owned();
310        for mut row in x_normalized.rows_mut().into_iter() {
311            row -= &self.state.feature_means;
312            row /= &self.state.feature_stds;
313        }
314
315        for label_idx in 0..self.state.n_labels {
316            let model = &self.state.models[label_idx];
317
318            for sample_idx in 0..n_samples {
319                let x_sample = x_normalized.row(sample_idx);
320
321                // Compute distances to both hyperplanes
322                let dist1 = (x_sample.dot(&model.w1) + model.b1).abs();
323                let dist2 = (x_sample.dot(&model.w2) + model.b2).abs();
324
325                // Predict based on closer hyperplane
326                predictions[[sample_idx, label_idx]] = if dist1 < dist2 { 1 } else { 0 };
327            }
328        }
329
330        Ok(predictions)
331    }
332}
333
334impl MLTSVM<MLTSVMTrained> {
335    /// Get the number of labels
336    pub fn n_labels(&self) -> usize {
337        self.state.n_labels
338    }
339
340    /// Get decision function values (distances to hyperplanes)
341    pub fn decision_function(&self, x: &ArrayView2<'_, Float>) -> SklResult<Array2<Float>> {
342        let (n_samples, _n_features) = x.dim();
343        let mut decision_values = Array2::<Float>::zeros((n_samples, self.state.n_labels));
344
345        // Normalize features
346        let mut x_normalized = x.to_owned();
347        for mut row in x_normalized.rows_mut().into_iter() {
348            row -= &self.state.feature_means;
349            row /= &self.state.feature_stds;
350        }
351
352        for label_idx in 0..self.state.n_labels {
353            let model = &self.state.models[label_idx];
354
355            for sample_idx in 0..n_samples {
356                let x_sample = x_normalized.row(sample_idx);
357
358                // Compute distances to both hyperplanes and use the difference
359                let dist1 = x_sample.dot(&model.w1) + model.b1;
360                let dist2 = x_sample.dot(&model.w2) + model.b2;
361
362                // Decision function is the difference (positive means class 1)
363                decision_values[[sample_idx, label_idx]] = dist1 - dist2;
364            }
365        }
366
367        Ok(decision_values)
368    }
369}
370
371/// RankSVM for Multi-Label Classification
372///
373/// RankSVM is a ranking-based approach for multi-label classification that optimizes
374/// ranking loss functions. It learns to rank labels by their relevance scores and
375/// can handle both label ranking and threshold selection for multi-label prediction.
376#[derive(Debug, Clone)]
377pub struct RankSVM<S = Untrained> {
378    state: S,
379    c: Float,                              // Regularization parameter
380    epsilon: Float,                        // Tolerance for convergence
381    max_iter: usize,                       // Maximum iterations
382    threshold_strategy: ThresholdStrategy, // How to determine prediction thresholds
383}
384
385/// Threshold strategy for RankSVM
386#[derive(Debug, Clone)]
387pub enum ThresholdStrategy {
388    /// Use fixed threshold for all labels
389    Fixed(Float),
390    /// Optimize threshold to maximize F1 score for each label
391    OptimizeF1,
392    /// Use top-k labels (fixed number of labels per sample)
393    TopK(usize),
394}
395
396/// Trained state for RankSVM
397#[derive(Debug, Clone)]
398pub struct RankSVMTrained {
399    models: Vec<RankingSVMModel>, // One model per label
400    thresholds: Vec<Float>,       // Prediction thresholds for each label
401    n_labels: usize,
402    feature_means: Array1<Float>,
403    feature_stds: Array1<Float>,
404}
405
406/// Single ranking SVM model for one label
407#[derive(Debug, Clone)]
408pub struct RankingSVMModel {
409    weights: Array1<Float>,
410    bias: Float,
411}
412
413impl RankSVM<Untrained> {
414    /// Create a new RankSVM instance
415    pub fn new() -> Self {
416        Self {
417            state: Untrained,
418            c: 1.0,
419            epsilon: 1e-3,
420            max_iter: 1000,
421            threshold_strategy: ThresholdStrategy::Fixed(0.0),
422        }
423    }
424
425    /// Set regularization parameter
426    pub fn c(mut self, c: Float) -> Self {
427        self.c = c;
428        self
429    }
430
431    /// Set convergence tolerance
432    pub fn epsilon(mut self, epsilon: Float) -> Self {
433        self.epsilon = epsilon;
434        self
435    }
436
437    /// Set maximum iterations
438    pub fn max_iter(mut self, max_iter: usize) -> Self {
439        self.max_iter = max_iter;
440        self
441    }
442
443    /// Set threshold strategy
444    pub fn threshold_strategy(mut self, strategy: ThresholdStrategy) -> Self {
445        self.threshold_strategy = strategy;
446        self
447    }
448}
449
450impl Default for RankSVM<Untrained> {
451    fn default() -> Self {
452        Self::new()
453    }
454}
455
456impl Estimator for RankSVM<Untrained> {
457    type Config = ();
458    type Error = SklearsError;
459    type Float = Float;
460
461    fn config(&self) -> &Self::Config {
462        &()
463    }
464}
465
466impl Fit<ArrayView2<'_, Float>, Array2<i32>> for RankSVM<Untrained> {
467    type Fitted = RankSVM<RankSVMTrained>;
468
469    fn fit(self, x: &ArrayView2<'_, Float>, y: &Array2<i32>) -> SklResult<Self::Fitted> {
470        let (n_samples, _n_features) = x.dim();
471        let (y_samples, n_labels) = y.dim();
472
473        if n_samples != y_samples {
474            return Err(SklearsError::InvalidInput(
475                "Number of samples in X and y must match".to_string(),
476            ));
477        }
478
479        // Validate that all labels are binary (0 or 1)
480        for sample_idx in 0..y_samples {
481            for label_idx in 0..n_labels {
482                let value = y[[sample_idx, label_idx]];
483                if value != 0 && value != 1 {
484                    return Err(SklearsError::InvalidInput(format!(
485                        "All label values must be 0 or 1, found: {}",
486                        value
487                    )));
488                }
489            }
490        }
491
492        // Compute feature statistics
493        let feature_means = x.mean_axis(Axis(0)).ok_or_else(|| {
494            SklearsError::InvalidInput("Cannot compute feature means from input data".to_string())
495        })?;
496
497        let squared_means = x.mapv(|val| val * val).mean_axis(Axis(0)).ok_or_else(|| {
498            SklearsError::InvalidInput("Cannot compute squared means from input data".to_string())
499        })?;
500
501        let feature_stds = squared_means - &feature_means.mapv(|mean| mean * mean);
502        let feature_stds = feature_stds.mapv(|var| (var.max(1e-10)).sqrt());
503
504        // Train ranking SVM for each label
505        let mut models = Vec::new();
506        for label_idx in 0..n_labels {
507            let y_label = y.column(label_idx);
508            let model = self.train_ranking_svm(x, &y_label, &feature_means, &feature_stds)?;
509            models.push(model);
510        }
511
512        // Determine thresholds
513        let thresholds = match &self.threshold_strategy {
514            ThresholdStrategy::Fixed(threshold) => vec![*threshold; n_labels],
515            ThresholdStrategy::OptimizeF1 => {
516                self.optimize_f1_thresholds(x, y, &models, &feature_means, &feature_stds)?
517            }
518            ThresholdStrategy::TopK(_) => vec![0.0; n_labels], // No threshold needed for TopK
519        };
520
521        Ok(RankSVM {
522            state: RankSVMTrained {
523                models,
524                thresholds,
525                n_labels,
526                feature_means,
527                feature_stds,
528            },
529            c: self.c,
530            epsilon: self.epsilon,
531            max_iter: self.max_iter,
532            threshold_strategy: self.threshold_strategy,
533        })
534    }
535}
536
537impl RankSVM<Untrained> {
538    fn train_ranking_svm(
539        &self,
540        x: &ArrayView2<'_, Float>,
541        y: &ArrayView1<'_, i32>,
542        feature_means: &Array1<Float>,
543        feature_stds: &Array1<Float>,
544    ) -> SklResult<RankingSVMModel> {
545        let (n_samples, n_features) = x.dim();
546
547        // Normalize features
548        let mut x_normalized = x.to_owned();
549        for mut row in x_normalized.rows_mut().into_iter() {
550            row -= feature_means;
551            row /= feature_stds;
552        }
553
554        // Initialize weights and bias
555        let mut weights = Array1::<Float>::zeros(n_features);
556        let mut bias = 0.0;
557
558        let learning_rate = 0.01;
559
560        // Gradient descent optimization
561        for _iter in 0..self.max_iter {
562            let mut weight_gradient = Array1::<Float>::zeros(n_features);
563            let mut bias_gradient = 0.0;
564
565            // Create ranking pairs
566            for i in 0..n_samples {
567                for j in 0..n_samples {
568                    if y[i] > y[j] {
569                        // i should be ranked higher than j
570                        let x_i = x_normalized.row(i);
571                        let x_j = x_normalized.row(j);
572                        let x_diff = &x_i.to_owned() - &x_j.to_owned();
573
574                        let score_diff = x_diff.dot(&weights) + bias;
575                        let margin = 1.0 - score_diff;
576
577                        if margin > 0.0 {
578                            // Hinge loss gradient
579                            weight_gradient -= &(x_diff * self.c);
580                            bias_gradient -= self.c;
581                        }
582                    }
583                }
584            }
585
586            // L2 regularization
587            weight_gradient += &(&weights * 2.0);
588
589            // Check convergence before updating parameters
590            let gradient_norm = weight_gradient.mapv(|x| x.abs()).sum();
591
592            // Update parameters
593            weights -= &(weight_gradient * learning_rate);
594            bias -= bias_gradient * learning_rate;
595
596            if gradient_norm < self.epsilon {
597                break;
598            }
599        }
600
601        Ok(RankingSVMModel { weights, bias })
602    }
603
604    fn optimize_f1_thresholds(
605        &self,
606        x: &ArrayView2<'_, Float>,
607        y: &Array2<i32>,
608        models: &[RankingSVMModel],
609        feature_means: &Array1<Float>,
610        feature_stds: &Array1<Float>,
611    ) -> SklResult<Vec<Float>> {
612        let mut thresholds = Vec::new();
613
614        for (label_idx, model) in models.iter().enumerate().take(y.ncols()) {
615            let y_true = y.column(label_idx);
616            let scores = self.predict_scores_single_label(x, model, feature_means, feature_stds)?;
617
618            let threshold = self.find_optimal_f1_threshold(&y_true, &scores)?;
619            thresholds.push(threshold);
620        }
621
622        Ok(thresholds)
623    }
624
625    fn predict_scores_single_label(
626        &self,
627        x: &ArrayView2<'_, Float>,
628        model: &RankingSVMModel,
629        feature_means: &Array1<Float>,
630        feature_stds: &Array1<Float>,
631    ) -> SklResult<Array1<Float>> {
632        let (n_samples, _) = x.dim();
633        let mut scores = Array1::<Float>::zeros(n_samples);
634
635        for i in 0..n_samples {
636            let x_sample = x.row(i);
637            let x_normalized = (&x_sample.to_owned() - feature_means) / feature_stds;
638            scores[i] = x_normalized.dot(&model.weights) + model.bias;
639        }
640
641        Ok(scores)
642    }
643
644    fn find_optimal_f1_threshold(
645        &self,
646        y_true: &ArrayView1<'_, i32>,
647        scores: &Array1<Float>,
648    ) -> SklResult<Float> {
649        let mut score_threshold_pairs: Vec<(Float, i32)> = scores
650            .iter()
651            .zip(y_true.iter())
652            .map(|(&score, &label)| (score, label))
653            .collect();
654
655        score_threshold_pairs
656            .sort_by(|a, b| a.0.partial_cmp(&b.0).expect("operation should succeed"));
657
658        let mut best_f1 = 0.0;
659        let mut best_threshold = 0.0;
660
661        // Try each unique score as a threshold
662        for &(threshold, _) in &score_threshold_pairs {
663            let mut tp = 0;
664            let mut fp = 0;
665            let mut fn_count = 0;
666
667            for (&score, &true_label) in scores.iter().zip(y_true.iter()) {
668                let predicted = if score >= threshold { 1 } else { 0 };
669
670                match (true_label, predicted) {
671                    (1, 1) => tp += 1,
672                    (0, 1) => fp += 1,
673                    (1, 0) => fn_count += 1,
674                    _ => {}
675                }
676            }
677
678            let precision = if tp + fp > 0 {
679                tp as Float / (tp + fp) as Float
680            } else {
681                0.0
682            };
683            let recall = if tp + fn_count > 0 {
684                tp as Float / (tp + fn_count) as Float
685            } else {
686                0.0
687            };
688            let f1 = if precision + recall > 0.0 {
689                2.0 * precision * recall / (precision + recall)
690            } else {
691                0.0
692            };
693
694            if f1 > best_f1 {
695                best_f1 = f1;
696                best_threshold = threshold;
697            }
698        }
699
700        Ok(best_threshold)
701    }
702}
703
704impl Predict<ArrayView2<'_, Float>, Array2<i32>> for RankSVM<RankSVMTrained> {
705    fn predict(&self, x: &ArrayView2<'_, Float>) -> SklResult<Array2<i32>> {
706        let (n_samples, n_features) = x.dim();
707        let expected_features = self.state.feature_means.len();
708
709        if n_features != expected_features {
710            return Err(SklearsError::InvalidInput(format!(
711                "Number of features in X ({}) does not match training data ({})",
712                n_features, expected_features
713            )));
714        }
715
716        let mut predictions = Array2::<i32>::zeros((n_samples, self.state.n_labels));
717
718        match &self.threshold_strategy {
719            ThresholdStrategy::TopK(k) => {
720                // For TopK, rank all labels and select top k
721                for sample_idx in 0..n_samples {
722                    let mut scores = Vec::new();
723                    for label_idx in 0..self.state.n_labels {
724                        let x_sample = x.row(sample_idx);
725                        let x_normalized = (&x_sample.to_owned() - &self.state.feature_means)
726                            / &self.state.feature_stds;
727                        let score = x_normalized.dot(&self.state.models[label_idx].weights)
728                            + self.state.models[label_idx].bias;
729                        scores.push((score, label_idx));
730                    }
731
732                    scores.sort_by(|a, b| b.0.partial_cmp(&a.0).expect("operation should succeed"));
733
734                    for &(_, label_idx) in scores.iter().take(*k) {
735                        predictions[[sample_idx, label_idx]] = 1;
736                    }
737                }
738            }
739            _ => {
740                // For fixed or optimized thresholds
741                for label_idx in 0..self.state.n_labels {
742                    let threshold = self.state.thresholds[label_idx];
743
744                    for sample_idx in 0..n_samples {
745                        let x_sample = x.row(sample_idx);
746                        let x_normalized = (&x_sample.to_owned() - &self.state.feature_means)
747                            / &self.state.feature_stds;
748                        let score = x_normalized.dot(&self.state.models[label_idx].weights)
749                            + self.state.models[label_idx].bias;
750
751                        predictions[[sample_idx, label_idx]] =
752                            if score >= threshold { 1 } else { 0 };
753                    }
754                }
755            }
756        }
757
758        Ok(predictions)
759    }
760}
761
762impl RankSVM<RankSVMTrained> {
763    /// Get the number of labels
764    pub fn n_labels(&self) -> usize {
765        self.state.n_labels
766    }
767
768    /// Get decision function values (ranking scores)
769    pub fn decision_function(&self, x: &ArrayView2<'_, Float>) -> SklResult<Array2<Float>> {
770        let (n_samples, n_features) = x.dim();
771        let expected_features = self.state.feature_means.len();
772
773        if n_features != expected_features {
774            return Err(SklearsError::InvalidInput(format!(
775                "Number of features in X ({}) does not match training data ({})",
776                n_features, expected_features
777            )));
778        }
779
780        let mut decision_values = Array2::<Float>::zeros((n_samples, self.state.n_labels));
781
782        for sample_idx in 0..n_samples {
783            for label_idx in 0..self.state.n_labels {
784                let x_sample = x.row(sample_idx);
785                let x_normalized =
786                    (&x_sample.to_owned() - &self.state.feature_means) / &self.state.feature_stds;
787                let score = x_normalized.dot(&self.state.models[label_idx].weights)
788                    + self.state.models[label_idx].bias;
789                decision_values[[sample_idx, label_idx]] = score;
790            }
791        }
792
793        Ok(decision_values)
794    }
795
796    /// Get ranking predictions (label indices ordered by relevance)
797    pub fn predict_ranking(&self, x: &ArrayView2<'_, Float>) -> SklResult<Array2<usize>> {
798        let (n_samples, n_features) = x.dim();
799        let expected_features = self.state.feature_means.len();
800
801        if n_features != expected_features {
802            return Err(SklearsError::InvalidInput(format!(
803                "Number of features in X ({}) does not match training data ({})",
804                n_features, expected_features
805            )));
806        }
807
808        let mut rankings = Array2::<usize>::zeros((n_samples, self.state.n_labels));
809
810        for sample_idx in 0..n_samples {
811            let mut scores = Vec::new();
812            for label_idx in 0..self.state.n_labels {
813                let x_sample = x.row(sample_idx);
814                let x_normalized =
815                    (&x_sample.to_owned() - &self.state.feature_means) / &self.state.feature_stds;
816                let score = x_normalized.dot(&self.state.models[label_idx].weights)
817                    + self.state.models[label_idx].bias;
818                scores.push((score, label_idx));
819            }
820
821            // Sort by score descending
822            scores.sort_by(|a, b| b.0.partial_cmp(&a.0).expect("operation should succeed"));
823
824            // Assign rankings
825            for (rank, &(_score, label_idx)) in scores.iter().enumerate() {
826                rankings[[sample_idx, rank]] = label_idx;
827            }
828        }
829
830        Ok(rankings)
831    }
832
833    /// Get the thresholds used for prediction
834    pub fn thresholds(&self) -> &Vec<Float> {
835        &self.state.thresholds
836    }
837}
838
839/// Multi-Output Support Vector Machine
840///
841/// A multi-output support vector machine that handles multiple regression or
842/// classification targets simultaneously by training separate SVM models for each output.
843#[derive(Debug, Clone)]
844pub struct MultiOutputSVM<S = Untrained> {
845    state: S,
846    kernel: SVMKernel,
847    c: Float,
848    epsilon: Float,
849    gamma: Option<Float>,
850}
851
852/// SVM Kernel types
853#[derive(Debug, Clone, Copy, PartialEq)]
854pub enum SVMKernel {
855    /// Linear kernel: K(x, y) = x^T y
856    Linear,
857    /// Polynomial kernel: K(x, y) = (gamma * x^T y + coef0)^degree
858    Polynomial {
859        degree: i32,
860        gamma: Float,
861        coef0: Float,
862    },
863    /// Radial Basis Function kernel: K(x, y) = exp(-gamma * ||x - y||^2)
864    Rbf { gamma: Float },
865    /// Sigmoid kernel: K(x, y) = tanh(gamma * x^T y + coef0)
866    Sigmoid { gamma: Float, coef0: Float },
867}
868
869/// Trained state for MultiOutputSVM
870#[derive(Debug, Clone)]
871pub struct MultiOutputSVMTrained {
872    models: Vec<SVMModel>,
873    n_outputs: usize,
874    feature_means: Array1<Float>,
875    feature_stds: Array1<Float>,
876}
877
878/// Single SVM model for one output
879#[derive(Debug, Clone)]
880pub struct SVMModel {
881    support_vectors: Array2<Float>,
882    support_coefficients: Array1<Float>,
883    bias: Float,
884    kernel: SVMKernel,
885}
886
887impl MultiOutputSVM<Untrained> {
888    /// Create a new MultiOutputSVM instance
889    pub fn new() -> Self {
890        Self {
891            state: Untrained,
892            kernel: SVMKernel::Rbf { gamma: 1.0 },
893            c: 1.0,
894            epsilon: 1e-3,
895            gamma: None,
896        }
897    }
898
899    /// Set the kernel
900    pub fn kernel(mut self, kernel: SVMKernel) -> Self {
901        self.kernel = kernel;
902        self
903    }
904
905    /// Set regularization parameter
906    pub fn c(mut self, c: Float) -> Self {
907        self.c = c;
908        self
909    }
910
911    /// Set tolerance for stopping criterion
912    pub fn epsilon(mut self, epsilon: Float) -> Self {
913        self.epsilon = epsilon;
914        self
915    }
916
917    /// Set gamma parameter (will override kernel-specific gamma)
918    pub fn gamma(mut self, gamma: Float) -> Self {
919        self.gamma = Some(gamma);
920        self
921    }
922}
923
924impl Default for MultiOutputSVM<Untrained> {
925    fn default() -> Self {
926        Self::new()
927    }
928}
929
930impl Estimator for MultiOutputSVM<Untrained> {
931    type Config = ();
932    type Error = SklearsError;
933    type Float = Float;
934
935    fn config(&self) -> &Self::Config {
936        &()
937    }
938}
939
940impl Fit<ArrayView2<'_, Float>, ArrayView2<'_, Float>> for MultiOutputSVM<Untrained> {
941    type Fitted = MultiOutputSVM<MultiOutputSVMTrained>;
942
943    fn fit(self, x: &ArrayView2<'_, Float>, y: &ArrayView2<'_, Float>) -> SklResult<Self::Fitted> {
944        let (n_samples, _n_features) = x.dim();
945        let (y_samples, n_outputs) = y.dim();
946
947        if n_samples != y_samples {
948            return Err(SklearsError::InvalidInput(
949                "Number of samples in X and y must match".to_string(),
950            ));
951        }
952
953        // Compute feature statistics
954        let feature_means = x
955            .mean_axis(Axis(0))
956            .expect("array should have elements for mean computation");
957        let feature_stds = x
958            .mapv(|val| val * val)
959            .mean_axis(Axis(0))
960            .expect("array should have elements for mean computation")
961            - &feature_means.mapv(|mean| mean * mean);
962        let feature_stds = feature_stds.mapv(|var| (var.max(1e-10)).sqrt());
963
964        // Update kernel gamma if specified
965        let kernel = if let Some(gamma) = self.gamma {
966            match self.kernel {
967                SVMKernel::Rbf { .. } => SVMKernel::Rbf { gamma },
968                SVMKernel::Polynomial { degree, coef0, .. } => SVMKernel::Polynomial {
969                    degree,
970                    gamma,
971                    coef0,
972                },
973                SVMKernel::Sigmoid { coef0, .. } => SVMKernel::Sigmoid { gamma, coef0 },
974                other => other,
975            }
976        } else {
977            self.kernel
978        };
979
980        // Train one SVM for each output
981        let mut models = Vec::new();
982        for output_idx in 0..n_outputs {
983            let y_output = y.column(output_idx);
984            let model =
985                self.train_single_svm(x, &y_output, &feature_means, &feature_stds, kernel)?;
986            models.push(model);
987        }
988
989        Ok(MultiOutputSVM {
990            state: MultiOutputSVMTrained {
991                models,
992                n_outputs,
993                feature_means,
994                feature_stds,
995            },
996            kernel,
997            c: self.c,
998            epsilon: self.epsilon,
999            gamma: self.gamma,
1000        })
1001    }
1002}
1003
1004impl MultiOutputSVM<Untrained> {
1005    fn train_single_svm(
1006        &self,
1007        x: &ArrayView2<'_, Float>,
1008        y: &ArrayView1<'_, Float>,
1009        feature_means: &Array1<Float>,
1010        feature_stds: &Array1<Float>,
1011        kernel: SVMKernel,
1012    ) -> SklResult<SVMModel> {
1013        let (n_samples, _n_features) = x.dim();
1014
1015        // Normalize features
1016        let mut x_normalized = x.to_owned();
1017        for mut row in x_normalized.rows_mut().into_iter() {
1018            row -= feature_means;
1019            row /= feature_stds;
1020        }
1021
1022        // For simplicity, we'll implement a basic SVM using all samples as support vectors
1023        // In a real implementation, you'd use SMO or other optimization algorithms
1024        let support_vectors = x_normalized.clone();
1025        let mut support_coefficients = Array1::<Float>::zeros(n_samples);
1026
1027        // Simple heuristic: coefficients proportional to target values
1028        let y_mean = y
1029            .mean()
1030            .expect("array should have elements for mean computation");
1031        for i in 0..n_samples {
1032            support_coefficients[i] = (y[i] - y_mean) / self.c;
1033        }
1034
1035        let bias = y_mean;
1036
1037        Ok(SVMModel {
1038            support_vectors,
1039            support_coefficients,
1040            bias,
1041            kernel,
1042        })
1043    }
1044}
1045
1046impl Predict<ArrayView2<'_, Float>, Array2<Float>> for MultiOutputSVM<MultiOutputSVMTrained> {
1047    fn predict(&self, x: &ArrayView2<'_, Float>) -> SklResult<Array2<Float>> {
1048        let (n_samples, _) = x.dim();
1049        let mut predictions = Array2::<Float>::zeros((n_samples, self.state.n_outputs));
1050
1051        // Normalize input features
1052        let mut x_normalized = x.to_owned();
1053        for mut row in x_normalized.rows_mut().into_iter() {
1054            row -= &self.state.feature_means;
1055            row /= &self.state.feature_stds;
1056        }
1057
1058        for output_idx in 0..self.state.n_outputs {
1059            let model = &self.state.models[output_idx];
1060
1061            for sample_idx in 0..n_samples {
1062                let x_sample = x_normalized.row(sample_idx);
1063                let mut prediction = model.bias;
1064
1065                // Compute kernel sum
1066                for (sv_idx, support_vector) in model.support_vectors.rows().into_iter().enumerate()
1067                {
1068                    let kernel_value =
1069                        compute_kernel_value(&x_sample, &support_vector, model.kernel);
1070                    prediction += model.support_coefficients[sv_idx] * kernel_value;
1071                }
1072
1073                predictions[[sample_idx, output_idx]] = prediction;
1074            }
1075        }
1076
1077        Ok(predictions)
1078    }
1079}
1080
1081/// Compute kernel value between two vectors
1082fn compute_kernel_value(
1083    x1: &ArrayView1<Float>,
1084    x2: &ArrayView1<Float>,
1085    kernel: SVMKernel,
1086) -> Float {
1087    match kernel {
1088        SVMKernel::Linear => x1.dot(x2),
1089        SVMKernel::Polynomial {
1090            degree,
1091            gamma,
1092            coef0,
1093        } => (gamma * x1.dot(x2) + coef0).powi(degree),
1094        SVMKernel::Rbf { gamma } => {
1095            let dist_sq = x1
1096                .iter()
1097                .zip(x2.iter())
1098                .map(|(a, b)| (a - b).powi(2))
1099                .sum::<Float>();
1100            (-gamma * dist_sq).exp()
1101        }
1102        SVMKernel::Sigmoid { gamma, coef0 } => (gamma * x1.dot(x2) + coef0).tanh(),
1103    }
1104}
sklears_multioutput/svm.rs

sklears_multioutput/
svm.rs
