Skip to main content

sklears_multioutput/
sparse_storage.rs

1//! Memory-efficient storage for sparse output representations
2//!
3//! This module provides optimized data structures and algorithms for scenarios where
4//! multi-output predictions are sparse (most outputs are zero or inactive).
5//! Common in multi-label classification where each instance typically has only a few active labels.
6#![allow(non_snake_case)] // Standard ML notation: X for feature matrices, K for kernels
7
8// Use SciRS2-Core for arrays and random number generation (SciRS2 Policy)
9use scirs2_core::ndarray::{Array1, Array2, ArrayView2};
10use sklears_core::{
11    error::{Result as SklResult, SklearsError},
12    traits::{Estimator, Fit, Predict, Untrained},
13    types::Float,
14};
15use std::collections::HashMap;
16use std::fmt;
17
18/// Compressed Sparse Row (CSR) format for efficient sparse matrix storage
19#[derive(Debug, Clone)]
20pub struct CSRMatrix<T: Clone> {
21    /// Non-zero values stored in row-major order
22    pub data: Vec<T>,
23    /// Column indices for each non-zero value
24    pub indices: Vec<usize>,
25    /// Pointers to the start of each row in data/indices
26    pub indptr: Vec<usize>,
27    /// Matrix dimensions (rows, cols)
28    pub shape: (usize, usize),
29}
30
31impl<T: Clone + Default + PartialEq> CSRMatrix<T> {
32    /// Create a new empty CSR matrix with given dimensions
33    pub fn new(rows: usize, cols: usize) -> Self {
34        Self {
35            data: Vec::new(),
36            indices: Vec::new(),
37            indptr: vec![0; rows + 1],
38            shape: (rows, cols),
39        }
40    }
41
42    /// Create CSR matrix from dense array
43    pub fn from_dense(dense: &ArrayView2<T>) -> Self
44    where
45        T: Clone + Default + PartialEq + Copy,
46    {
47        let (rows, cols) = dense.dim();
48        let mut data = Vec::new();
49        let mut indices = Vec::new();
50        let mut indptr = vec![0; rows + 1];
51
52        for row in 0..rows {
53            for col in 0..cols {
54                let val = dense[[row, col]];
55                if val != T::default() {
56                    data.push(val);
57                    indices.push(col);
58                }
59            }
60            indptr[row + 1] = data.len();
61        }
62
63        Self {
64            data,
65            indices,
66            indptr,
67            shape: (rows, cols),
68        }
69    }
70
71    /// Convert back to dense array
72    pub fn to_dense(&self) -> Array2<T>
73    where
74        T: Clone + Default,
75    {
76        let (rows, cols) = self.shape;
77        let mut dense = Array2::from_elem((rows, cols), T::default());
78
79        for row in 0..rows {
80            let start = self.indptr[row];
81            let end = self.indptr[row + 1];
82
83            for idx in start..end {
84                let col = self.indices[idx];
85                let val = self.data[idx].clone();
86                dense[[row, col]] = val;
87            }
88        }
89
90        dense
91    }
92
93    /// Get the number of non-zero elements
94    pub fn nnz(&self) -> usize {
95        self.data.len()
96    }
97
98    /// Calculate sparsity ratio (fraction of non-zero elements)
99    pub fn sparsity(&self) -> f64 {
100        let total_elements = self.shape.0 * self.shape.1;
101        if total_elements == 0 {
102            0.0
103        } else {
104            self.nnz() as f64 / total_elements as f64
105        }
106    }
107
108    /// Get values for a specific row
109    pub fn get_row(&self, row: usize) -> Vec<(usize, T)> {
110        if row >= self.shape.0 {
111            return Vec::new();
112        }
113
114        let start = self.indptr[row];
115        let end = self.indptr[row + 1];
116        let mut row_data = Vec::new();
117
118        for idx in start..end {
119            let col = self.indices[idx];
120            let val = self.data[idx].clone();
121            row_data.push((col, val));
122        }
123
124        row_data
125    }
126
127    /// Set a value at specific row and column
128    pub fn set(&mut self, row: usize, col: usize, value: T) {
129        if row >= self.shape.0 || col >= self.shape.1 {
130            return;
131        }
132
133        let start = self.indptr[row];
134        let end = self.indptr[row + 1];
135
136        // Find if the element already exists
137        for idx in start..end {
138            if self.indices[idx] == col {
139                if value == T::default() {
140                    // Remove the element
141                    self.data.remove(idx);
142                    self.indices.remove(idx);
143                    // Update indptr for all following rows
144                    for r in (row + 1)..=self.shape.0 {
145                        self.indptr[r] -= 1;
146                    }
147                } else {
148                    // Update the value
149                    self.data[idx] = value;
150                }
151                return;
152            }
153            if self.indices[idx] > col {
154                // Insert at this position
155                if value != T::default() {
156                    self.data.insert(idx, value);
157                    self.indices.insert(idx, col);
158                    // Update indptr for all following rows
159                    for r in (row + 1)..=self.shape.0 {
160                        self.indptr[r] += 1;
161                    }
162                }
163                return;
164            }
165        }
166
167        // Append at the end of this row
168        if value != T::default() {
169            self.data.insert(end, value);
170            self.indices.insert(end, col);
171            // Update indptr for all following rows
172            for r in (row + 1)..=self.shape.0 {
173                self.indptr[r] += 1;
174            }
175        }
176    }
177}
178
179/// Memory-efficient sparse multi-output predictor
180#[derive(Debug, Clone)]
181pub struct SparseMultiOutput<S = Untrained> {
182    state: S,
183    /// Sparsity threshold - values below this are considered zero
184    sparsity_threshold: f64,
185    /// Whether to use compressed storage for predictions
186    use_compression: bool,
187}
188
189/// Trained state for sparse multi-output predictor
190#[derive(Debug, Clone)]
191pub struct SparseMultiOutputTrained {
192    pub coefficients: CSRMatrix<f64>,
193    pub bias: HashMap<usize, f64>,
194    pub feature_means: Array1<f64>,
195    pub feature_stds: Array1<f64>,
196    pub n_features: usize,
197    pub n_outputs: usize,
198    pub sparsity_ratio: f64,
199}
200
201impl SparseMultiOutput<Untrained> {
202    /// Create a new sparse multi-output predictor
203    pub fn new() -> Self {
204        Self {
205            state: Untrained,
206            sparsity_threshold: 1e-6,
207            use_compression: true,
208        }
209    }
210
211    /// Set the sparsity threshold
212    pub fn sparsity_threshold(mut self, threshold: f64) -> Self {
213        self.sparsity_threshold = threshold;
214        self
215    }
216
217    /// Enable or disable compression
218    pub fn use_compression(mut self, use_compression: bool) -> Self {
219        self.use_compression = use_compression;
220        self
221    }
222}
223
224impl Default for SparseMultiOutput<Untrained> {
225    fn default() -> Self {
226        Self::new()
227    }
228}
229
230impl Estimator for SparseMultiOutput<Untrained> {
231    type Config = ();
232    type Error = SklearsError;
233    type Float = Float;
234
235    fn config(&self) -> &Self::Config {
236        &()
237    }
238}
239
240impl Estimator for SparseMultiOutput<SparseMultiOutputTrained> {
241    type Config = ();
242    type Error = SklearsError;
243    type Float = Float;
244
245    fn config(&self) -> &Self::Config {
246        &()
247    }
248}
249
250impl Fit<ArrayView2<'_, Float>, ArrayView2<'_, f64>> for SparseMultiOutput<Untrained> {
251    type Fitted = SparseMultiOutput<SparseMultiOutputTrained>;
252
253    #[allow(non_snake_case)]
254    fn fit(self, X: &ArrayView2<'_, Float>, y: &ArrayView2<'_, f64>) -> SklResult<Self::Fitted> {
255        let (n_samples, n_features) = X.dim();
256        let (n_samples_y, n_outputs) = y.dim();
257
258        if n_samples != n_samples_y {
259            return Err(SklearsError::InvalidInput(
260                "X and y must have the same number of samples".to_string(),
261            ));
262        }
263
264        // Convert X to f64 for consistency
265        let X_f64 = X.mapv(|x| x);
266
267        // Compute feature statistics for standardization
268        let mut feature_means = Array1::zeros(n_features);
269        let mut feature_stds = Array1::zeros(n_features);
270
271        for feature in 0..n_features {
272            let col = X_f64.column(feature);
273            feature_means[feature] = col.sum() / n_samples as f64;
274
275            let variance = col
276                .iter()
277                .map(|&x| (x - feature_means[feature]).powi(2))
278                .sum::<f64>()
279                / n_samples as f64;
280            feature_stds[feature] = variance.sqrt().max(1e-8); // Avoid division by zero
281        }
282
283        // Standardize X
284        let mut X_std = X_f64.clone();
285        for feature in 0..n_features {
286            let mut col = X_std.column_mut(feature);
287            col -= feature_means[feature];
288            col /= feature_stds[feature];
289        }
290
291        // Train sparse linear models using coordinate descent
292        let mut coefficients_dense = Array2::zeros((n_outputs, n_features));
293        let mut bias = HashMap::new();
294
295        for output in 0..n_outputs {
296            let y_target = y.column(output);
297
298            // Simple ridge regression for each output
299            let mut weights = Array1::zeros(n_features);
300            let intercept = y_target.mean().unwrap_or(0.0);
301
302            // Coordinate descent iterations
303            for _iter in 0..100 {
304                let mut converged = true;
305
306                // Update each weight
307                for feature in 0..n_features {
308                    let old_weight = weights[feature];
309
310                    // Compute residuals without this feature
311                    let mut residual_sum = 0.0;
312                    for sample in 0..n_samples {
313                        let mut pred = intercept;
314                        for other_feature in 0..n_features {
315                            if other_feature != feature {
316                                pred += weights[other_feature] * X_std[[sample, other_feature]];
317                            }
318                        }
319                        let residual = y_target[sample] - pred;
320                        residual_sum += residual * X_std[[sample, feature]];
321                    }
322
323                    // Feature variance (standardized features have variance 1)
324                    let feature_var = n_samples as f64;
325
326                    // Ridge penalty
327                    let lambda = 0.01;
328                    let new_weight = residual_sum / (feature_var + lambda);
329
330                    // Apply sparsity threshold
331                    weights[feature] = if new_weight.abs() < self.sparsity_threshold {
332                        0.0
333                    } else {
334                        new_weight
335                    };
336
337                    if (weights[feature] - old_weight).abs() > 1e-6 {
338                        converged = false;
339                    }
340                }
341
342                if converged {
343                    break;
344                }
345            }
346
347            // Store results
348            for feature in 0..n_features {
349                coefficients_dense[[output, feature]] = weights[feature];
350            }
351
352            if intercept.abs() > self.sparsity_threshold {
353                bias.insert(output, intercept);
354            }
355        }
356
357        // Convert to sparse format
358        let coefficients = CSRMatrix::from_dense(&coefficients_dense.view());
359        let sparsity_ratio = coefficients.sparsity();
360
361        Ok(SparseMultiOutput {
362            state: SparseMultiOutputTrained {
363                coefficients,
364                bias,
365                feature_means,
366                feature_stds,
367                n_features,
368                n_outputs,
369                sparsity_ratio,
370            },
371            sparsity_threshold: self.sparsity_threshold,
372            use_compression: self.use_compression,
373        })
374    }
375}
376
377impl Predict<ArrayView2<'_, Float>, Array2<f64>> for SparseMultiOutput<SparseMultiOutputTrained> {
378    fn predict(&self, X: &ArrayView2<'_, Float>) -> SklResult<Array2<f64>> {
379        let (n_samples, n_features) = X.dim();
380
381        if n_features != self.state.n_features {
382            return Err(SklearsError::InvalidInput(format!(
383                "Expected {} features, got {}",
384                self.state.n_features, n_features
385            )));
386        }
387
388        // Standardize input features
389        let mut X_std = X.mapv(|x| x);
390        for feature in 0..n_features {
391            let mut col = X_std.column_mut(feature);
392            col -= self.state.feature_means[feature];
393            col /= self.state.feature_stds[feature];
394        }
395
396        let mut predictions = Array2::zeros((n_samples, self.state.n_outputs));
397
398        // Sparse matrix-vector multiplication
399        for output in 0..self.state.n_outputs {
400            let output_coeffs = self.state.coefficients.get_row(output);
401            let intercept = *self.state.bias.get(&output).unwrap_or(&0.0);
402
403            for sample in 0..n_samples {
404                let mut pred = intercept;
405
406                // Only compute for non-zero coefficients
407                for &(feature, coeff) in &output_coeffs {
408                    pred += coeff * X_std[[sample, feature]];
409                }
410
411                predictions[[sample, output]] = pred;
412            }
413        }
414
415        Ok(predictions)
416    }
417}
418
419impl SparseMultiOutput<SparseMultiOutputTrained> {
420    /// Get the sparsity ratio of the coefficient matrix
421    pub fn sparsity_ratio(&self) -> f64 {
422        self.state.sparsity_ratio
423    }
424
425    /// Get the number of non-zero coefficients
426    pub fn nnz_coefficients(&self) -> usize {
427        self.state.coefficients.nnz()
428    }
429
430    /// Get memory usage statistics
431    pub fn memory_usage(&self) -> MemoryUsage {
432        let dense_size = self.state.n_outputs * self.state.n_features * 8; // 8 bytes per f64
433        let sparse_size = self.state.coefficients.data.len() * 8 + // data values
434                         self.state.coefficients.indices.len() * 8 + // column indices
435                         self.state.coefficients.indptr.len() * 8; // row pointers
436
437        let compression_ratio = if dense_size > 0 {
438            sparse_size as f64 / dense_size as f64
439        } else {
440            1.0
441        };
442
443        MemoryUsage {
444            dense_size_bytes: dense_size,
445            sparse_size_bytes: sparse_size,
446            compression_ratio,
447            memory_saved_bytes: dense_size.saturating_sub(sparse_size),
448        }
449    }
450
451    /// Get coefficients for a specific output (sparse representation)
452    pub fn get_output_coefficients(&self, output: usize) -> Vec<(usize, f64)> {
453        if output >= self.state.n_outputs {
454            return Vec::new();
455        }
456
457        self.state.coefficients.get_row(output)
458    }
459
460    /// Get the bias for a specific output
461    pub fn get_output_bias(&self, output: usize) -> f64 {
462        *self.state.bias.get(&output).unwrap_or(&0.0)
463    }
464}
465
466/// Memory usage statistics
467#[derive(Debug, Clone)]
468pub struct MemoryUsage {
469    /// Size of dense representation in bytes
470    pub dense_size_bytes: usize,
471    /// Size of sparse representation in bytes
472    pub sparse_size_bytes: usize,
473    /// Compression ratio (sparse_size / dense_size)
474    pub compression_ratio: f64,
475    /// Memory saved in bytes
476    pub memory_saved_bytes: usize,
477}
478
479impl fmt::Display for MemoryUsage {
480    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
481        write!(f,
482            "Memory Usage - Dense: {} bytes, Sparse: {} bytes, Compression: {:.3}x, Saved: {} bytes",
483            self.dense_size_bytes,
484            self.sparse_size_bytes,
485            self.compression_ratio,
486            self.memory_saved_bytes
487        )
488    }
489}
490
491/// Utility functions for sparse output analysis
492pub mod sparse_utils {
493    use super::*;
494
495    /// Analyze sparsity patterns in output data
496    pub fn analyze_output_sparsity(y: &ArrayView2<f64>, threshold: f64) -> SparsityAnalysis {
497        let (n_samples, n_outputs) = y.dim();
498        let mut total_elements = 0;
499        let mut zero_elements = 0;
500        let mut output_sparsities = Vec::with_capacity(n_outputs);
501
502        for output in 0..n_outputs {
503            let col = y.column(output);
504            let output_zeros = col.iter().filter(|&&x| x.abs() <= threshold).count();
505            let output_sparsity = output_zeros as f64 / n_samples as f64;
506            output_sparsities.push(output_sparsity);
507
508            total_elements += n_samples;
509            zero_elements += output_zeros;
510        }
511
512        let overall_sparsity = zero_elements as f64 / total_elements as f64;
513        let avg_sparsity = output_sparsities.iter().sum::<f64>() / n_outputs as f64;
514        let min_sparsity = output_sparsities
515            .iter()
516            .fold(f64::INFINITY, |a, &b| a.min(b));
517        let max_sparsity = output_sparsities
518            .iter()
519            .fold(f64::NEG_INFINITY, |a, &b| a.max(b));
520
521        SparsityAnalysis {
522            overall_sparsity,
523            avg_sparsity,
524            min_sparsity,
525            max_sparsity,
526            output_sparsities,
527            total_elements,
528            zero_elements,
529        }
530    }
531
532    /// Recommend whether to use sparse storage based on data characteristics
533    pub fn recommend_sparse_storage(y: &ArrayView2<f64>, threshold: f64) -> StorageRecommendation {
534        let analysis = analyze_output_sparsity(y, threshold);
535
536        let should_use_sparse = analysis.overall_sparsity > 0.5; // More than 50% zeros
537        let expected_compression = if should_use_sparse {
538            // Estimate compression based on sparsity
539            1.0 - analysis.overall_sparsity + 0.1 // Add overhead estimate
540        } else {
541            1.0
542        };
543
544        StorageRecommendation {
545            should_use_sparse,
546            expected_compression_ratio: expected_compression,
547            sparsity_analysis: analysis,
548        }
549    }
550}
551
552/// Sparsity analysis results
553#[derive(Debug, Clone)]
554pub struct SparsityAnalysis {
555    /// Overall fraction of zero elements
556    pub overall_sparsity: f64,
557    /// Average sparsity across outputs
558    pub avg_sparsity: f64,
559    /// Minimum sparsity among outputs
560    pub min_sparsity: f64,
561    /// Maximum sparsity among outputs
562    pub max_sparsity: f64,
563    /// Sparsity for each output
564    pub output_sparsities: Vec<f64>,
565    /// Total number of elements
566    pub total_elements: usize,
567    /// Number of zero elements
568    pub zero_elements: usize,
569}
570
571/// Storage recommendation based on data analysis
572#[derive(Debug, Clone)]
573pub struct StorageRecommendation {
574    /// Whether sparse storage is recommended
575    pub should_use_sparse: bool,
576    /// Expected compression ratio
577    pub expected_compression_ratio: f64,
578    /// Detailed sparsity analysis
579    pub sparsity_analysis: SparsityAnalysis,
580}
581
582#[allow(non_snake_case)]
583#[cfg(test)]
584mod tests {
585    use super::*;
586    use approx::assert_abs_diff_eq;
587    // Use SciRS2-Core for arrays and random number generation (SciRS2 Policy)
588    use scirs2_core::ndarray::array;
589
590    #[test]
591    fn test_csr_matrix_basic() {
592        let dense = array![[1.0, 0.0, 3.0], [0.0, 2.0, 0.0], [4.0, 0.0, 0.0]];
593        let csr = CSRMatrix::from_dense(&dense.view());
594
595        assert_eq!(csr.nnz(), 4);
596        assert_eq!(csr.shape, (3, 3));
597        assert_eq!(csr.data, vec![1.0, 3.0, 2.0, 4.0]);
598        assert_eq!(csr.indices, vec![0, 2, 1, 0]);
599        assert_eq!(csr.indptr, vec![0, 2, 3, 4]);
600
601        let reconstructed = csr.to_dense();
602        for i in 0..3 {
603            for j in 0..3 {
604                assert_abs_diff_eq!(dense[[i, j]], reconstructed[[i, j]], epsilon = 1e-10);
605            }
606        }
607    }
608
609    #[test]
610    fn test_csr_sparsity() {
611        let dense = array![[1.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 3.0]];
612        let csr = CSRMatrix::from_dense(&dense.view());
613
614        assert_abs_diff_eq!(csr.sparsity(), 2.0 / 9.0, epsilon = 1e-10);
615    }
616
617    #[test]
618    #[allow(non_snake_case)]
619    fn test_sparse_multi_output_basic() {
620        let X = array![[1.0, 2.0], [3.0, 4.0], [5.0, 6.0], [7.0, 8.0]];
621        let y = array![
622            [1.0, 0.0, 0.1],
623            [0.0, 2.0, 0.0],
624            [3.0, 0.0, 0.0],
625            [0.0, 4.0, 0.2]
626        ];
627
628        let model = SparseMultiOutput::new().sparsity_threshold(0.05);
629        let trained = model
630            .fit(&X.view(), &y.view())
631            .expect("model fitting should succeed");
632
633        let predictions = trained
634            .predict(&X.view())
635            .expect("prediction should succeed");
636        assert_eq!(predictions.shape(), &[4, 3]);
637
638        // Check that model learned something reasonable
639        assert!(trained.sparsity_ratio() < 1.0); // Should have some non-zero coefficients
640        println!("Sparsity ratio: {}", trained.sparsity_ratio());
641    }
642
643    #[test]
644    #[allow(non_snake_case)]
645    fn test_sparse_memory_efficiency() {
646        let X = array![
647            [1.0, 2.0, 3.0, 4.0, 5.0],
648            [2.0, 3.0, 4.0, 5.0, 6.0],
649            [3.0, 4.0, 5.0, 6.0, 7.0]
650        ];
651        // Highly sparse output - most values are zero
652        let y = array![
653            [1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
654            [0.0, 2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
655            [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3.0]
656        ];
657
658        let model = SparseMultiOutput::new().sparsity_threshold(1e-6);
659        let trained = model
660            .fit(&X.view(), &y.view())
661            .expect("model fitting should succeed");
662
663        let memory_usage = trained.memory_usage();
664        println!("{}", memory_usage);
665
666        // Should achieve significant compression for sparse data
667        assert!(memory_usage.compression_ratio < 0.8); // At least 20% compression
668        assert!(memory_usage.memory_saved_bytes > 0);
669    }
670
671    #[test]
672    fn test_sparsity_analysis() {
673        let y = array![
674            [1.0, 0.0, 0.0, 2.0],
675            [0.0, 0.0, 3.0, 0.0],
676            [0.0, 1.0, 0.0, 0.0],
677            [2.0, 0.0, 0.0, 0.0]
678        ];
679
680        let analysis = sparse_utils::analyze_output_sparsity(&y.view(), 1e-6);
681
682        // Actually 11 out of 16 elements are zero (5 non-zero: 1.0, 2.0, 3.0, 1.0, 2.0)
683        assert_abs_diff_eq!(analysis.overall_sparsity, 11.0 / 16.0, epsilon = 1e-10);
684        assert_eq!(analysis.total_elements, 16);
685        assert_eq!(analysis.zero_elements, 11);
686        assert_eq!(analysis.output_sparsities.len(), 4);
687    }
688
689    #[test]
690    fn test_storage_recommendation() {
691        // Sparse data
692        let y_sparse = array![
693            [1.0, 0.0, 0.0, 0.0, 0.0],
694            [0.0, 0.0, 0.0, 2.0, 0.0],
695            [0.0, 0.0, 0.0, 0.0, 0.0]
696        ];
697
698        let recommendation = sparse_utils::recommend_sparse_storage(&y_sparse.view(), 1e-6);
699        assert!(recommendation.should_use_sparse);
700        assert!(recommendation.expected_compression_ratio < 1.0);
701
702        // Dense data
703        let y_dense = array![
704            [1.0, 2.0, 3.0, 4.0, 5.0],
705            [6.0, 7.0, 8.0, 9.0, 10.0],
706            [11.0, 12.0, 13.0, 14.0, 15.0]
707        ];
708
709        let recommendation = sparse_utils::recommend_sparse_storage(&y_dense.view(), 1e-6);
710        assert!(!recommendation.should_use_sparse);
711    }
712
713    #[test]
714    #[allow(non_snake_case)]
715    fn test_sparse_coefficient_access() {
716        let X = array![[1.0, 2.0], [3.0, 4.0]];
717        let y = array![[1.0, 0.0], [0.0, 2.0]];
718
719        let model = SparseMultiOutput::new().sparsity_threshold(1e-3);
720        let trained = model
721            .fit(&X.view(), &y.view())
722            .expect("model fitting should succeed");
723
724        // Test coefficient access for each output
725        for output in 0..2 {
726            let coeffs = trained.get_output_coefficients(output);
727            let bias = trained.get_output_bias(output);
728
729            println!("Output {}: coeffs = {:?}, bias = {}", output, coeffs, bias);
730
731            // Should have some coefficients
732            assert!(!coeffs.is_empty() || bias.abs() > 1e-6);
733        }
734    }
735
736    #[test]
737    #[allow(non_snake_case)]
738    fn test_edge_cases() {
739        let X = array![[1.0, 2.0], [3.0, 4.0]];
740
741        // All zeros - coefficients should be small due to regularization
742        let y_zeros = array![[0.0, 0.0], [0.0, 0.0]];
743        let model = SparseMultiOutput::new().sparsity_threshold(1e-3);
744        let trained = model
745            .fit(&X.view(), &y_zeros.view())
746            .expect("model fitting should succeed");
747
748        // Check that predictions are close to zero
749        let pred_zeros = trained
750            .predict(&X.view())
751            .expect("prediction should succeed");
752        for i in 0..pred_zeros.nrows() {
753            for j in 0..pred_zeros.ncols() {
754                assert!(
755                    pred_zeros[[i, j]].abs() < 0.1,
756                    "Prediction should be close to zero: {}",
757                    pred_zeros[[i, j]]
758                );
759            }
760        }
761
762        println!("Zero data sparsity ratio: {}", trained.sparsity_ratio());
763
764        // Single feature
765        let X_single = array![[1.0], [2.0]];
766        let y_single = array![[1.0], [2.0]];
767        let model_single = SparseMultiOutput::new();
768        let trained_single = model_single
769            .fit(&X_single.view(), &y_single.view())
770            .expect("operation should succeed");
771        let pred = trained_single
772            .predict(&X_single.view())
773            .expect("prediction should succeed");
774        assert_eq!(pred.shape(), &[2, 1]);
775
776        // Test with many zero outputs
777        let X_many = array![[1.0, 2.0], [3.0, 4.0]];
778        let y_many_sparse = array![[1.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 2.0]];
779        let model_many = SparseMultiOutput::new().sparsity_threshold(1e-6);
780        let trained_many = model_many
781            .fit(&X_many.view(), &y_many_sparse.view())
782            .expect("operation should succeed");
783
784        // Just check that training completed and we can make predictions
785        let pred_many = trained_many
786            .predict(&X_many.view())
787            .expect("prediction should succeed");
788        assert_eq!(pred_many.shape(), &[2, 5]);
789
790        println!(
791            "Many sparse outputs sparsity ratio: {}",
792            trained_many.sparsity_ratio()
793        );
794    }
795}