1#![allow(non_snake_case)] use 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#[derive(Debug, Clone)]
20pub struct CSRMatrix<T: Clone> {
21 pub data: Vec<T>,
23 pub indices: Vec<usize>,
25 pub indptr: Vec<usize>,
27 pub shape: (usize, usize),
29}
30
31impl<T: Clone + Default + PartialEq> CSRMatrix<T> {
32 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 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 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 pub fn nnz(&self) -> usize {
95 self.data.len()
96 }
97
98 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 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 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 for idx in start..end {
138 if self.indices[idx] == col {
139 if value == T::default() {
140 self.data.remove(idx);
142 self.indices.remove(idx);
143 for r in (row + 1)..=self.shape.0 {
145 self.indptr[r] -= 1;
146 }
147 } else {
148 self.data[idx] = value;
150 }
151 return;
152 }
153 if self.indices[idx] > col {
154 if value != T::default() {
156 self.data.insert(idx, value);
157 self.indices.insert(idx, col);
158 for r in (row + 1)..=self.shape.0 {
160 self.indptr[r] += 1;
161 }
162 }
163 return;
164 }
165 }
166
167 if value != T::default() {
169 self.data.insert(end, value);
170 self.indices.insert(end, col);
171 for r in (row + 1)..=self.shape.0 {
173 self.indptr[r] += 1;
174 }
175 }
176 }
177}
178
179#[derive(Debug, Clone)]
181pub struct SparseMultiOutput<S = Untrained> {
182 state: S,
183 sparsity_threshold: f64,
185 use_compression: bool,
187}
188
189#[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 pub fn new() -> Self {
204 Self {
205 state: Untrained,
206 sparsity_threshold: 1e-6,
207 use_compression: true,
208 }
209 }
210
211 pub fn sparsity_threshold(mut self, threshold: f64) -> Self {
213 self.sparsity_threshold = threshold;
214 self
215 }
216
217 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 let X_f64 = X.mapv(|x| x);
266
267 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); }
282
283 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 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 let mut weights = Array1::zeros(n_features);
300 let intercept = y_target.mean().unwrap_or(0.0);
301
302 for _iter in 0..100 {
304 let mut converged = true;
305
306 for feature in 0..n_features {
308 let old_weight = weights[feature];
309
310 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 let feature_var = n_samples as f64;
325
326 let lambda = 0.01;
328 let new_weight = residual_sum / (feature_var + lambda);
329
330 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 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 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 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 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 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 pub fn sparsity_ratio(&self) -> f64 {
422 self.state.sparsity_ratio
423 }
424
425 pub fn nnz_coefficients(&self) -> usize {
427 self.state.coefficients.nnz()
428 }
429
430 pub fn memory_usage(&self) -> MemoryUsage {
432 let dense_size = self.state.n_outputs * self.state.n_features * 8; let sparse_size = self.state.coefficients.data.len() * 8 + self.state.coefficients.indices.len() * 8 + self.state.coefficients.indptr.len() * 8; 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 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 pub fn get_output_bias(&self, output: usize) -> f64 {
462 *self.state.bias.get(&output).unwrap_or(&0.0)
463 }
464}
465
466#[derive(Debug, Clone)]
468pub struct MemoryUsage {
469 pub dense_size_bytes: usize,
471 pub sparse_size_bytes: usize,
473 pub compression_ratio: f64,
475 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
491pub mod sparse_utils {
493 use super::*;
494
495 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 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; let expected_compression = if should_use_sparse {
538 1.0 - analysis.overall_sparsity + 0.1 } 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#[derive(Debug, Clone)]
554pub struct SparsityAnalysis {
555 pub overall_sparsity: f64,
557 pub avg_sparsity: f64,
559 pub min_sparsity: f64,
561 pub max_sparsity: f64,
563 pub output_sparsities: Vec<f64>,
565 pub total_elements: usize,
567 pub zero_elements: usize,
569}
570
571#[derive(Debug, Clone)]
573pub struct StorageRecommendation {
574 pub should_use_sparse: bool,
576 pub expected_compression_ratio: f64,
578 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::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 assert!(trained.sparsity_ratio() < 1.0); 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 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 assert!(memory_usage.compression_ratio < 0.8); 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 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 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 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 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 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 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 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 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 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 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}