1use scirs2_core::ndarray_ext::{Array1, Array2, ArrayView1, ArrayView2, Axis};
7use scirs2_core::random::Random;
8use sklears_core::{
9 error::{Result as SklResult, SklearsError},
10 traits::{Estimator, Fit, Predict, PredictProba, Untrained},
11 types::Float,
12};
13
14#[derive(Debug, Clone)]
52pub struct GraphStructureLearning<S = Untrained> {
53 state: S,
54 lambda_sparse: f64,
55 beta_smoothness: f64,
56 max_iter: usize,
57 tol: f64,
58 learning_rate: f64,
59 adaptive_lr: bool,
60 enforce_symmetry: bool,
61 normalize_weights: bool,
62}
63
64impl GraphStructureLearning<Untrained> {
65 pub fn new() -> Self {
67 Self {
68 state: Untrained,
69 lambda_sparse: 0.1,
70 beta_smoothness: 1.0,
71 max_iter: 100,
72 tol: 1e-4,
73 learning_rate: 0.01,
74 adaptive_lr: true,
75 enforce_symmetry: true,
76 normalize_weights: true,
77 }
78 }
79
80 pub fn lambda_sparse(mut self, lambda_sparse: f64) -> Self {
82 self.lambda_sparse = lambda_sparse;
83 self
84 }
85
86 pub fn beta_smoothness(mut self, beta_smoothness: f64) -> Self {
88 self.beta_smoothness = beta_smoothness;
89 self
90 }
91
92 pub fn max_iter(mut self, max_iter: usize) -> Self {
94 self.max_iter = max_iter;
95 self
96 }
97
98 pub fn tol(mut self, tol: f64) -> Self {
100 self.tol = tol;
101 self
102 }
103
104 pub fn learning_rate(mut self, learning_rate: f64) -> Self {
106 self.learning_rate = learning_rate;
107 self
108 }
109
110 pub fn adaptive_lr(mut self, adaptive_lr: bool) -> Self {
112 self.adaptive_lr = adaptive_lr;
113 self
114 }
115
116 pub fn enforce_symmetry(mut self, enforce_symmetry: bool) -> Self {
118 self.enforce_symmetry = enforce_symmetry;
119 self
120 }
121
122 pub fn normalize_weights(mut self, normalize_weights: bool) -> Self {
124 self.normalize_weights = normalize_weights;
125 self
126 }
127
128 #[allow(non_snake_case)] fn initialize_graph(&self, X: &Array2<f64>) -> Array2<f64> {
130 let n_samples = X.nrows();
131 let mut W = Array2::zeros((n_samples, n_samples));
132
133 let k = (n_samples as f64).sqrt().ceil() as usize;
135 let k = k.clamp(3, 10); for i in 0..n_samples {
138 let mut distances: Vec<(usize, f64)> = Vec::new();
139 for j in 0..n_samples {
140 if i != j {
141 let diff = &X.row(i) - &X.row(j);
142 let dist = diff.mapv(|x: f64| x * x).sum().sqrt();
143 distances.push((j, dist));
144 }
145 }
146
147 distances.sort_by(|a, b| a.1.partial_cmp(&b.1).expect("operation should succeed"));
148
149 for &(j, dist) in distances.iter().take(k) {
150 let weight = (-dist / (2.0 * 1.0_f64.powi(2))).exp();
151 W[[i, j]] = weight;
152 if self.enforce_symmetry {
153 W[[j, i]] = weight;
154 }
155 }
156 }
157
158 W
159 }
160
161 #[allow(non_snake_case)]
162 fn compute_laplacian(&self, W: &Array2<f64>) -> Array2<f64> {
163 let n_samples = W.nrows();
164 let D = W.sum_axis(Axis(1));
165 let mut L = Array2::zeros((n_samples, n_samples));
166
167 for i in 0..n_samples {
168 L[[i, i]] = D[i];
169 for j in 0..n_samples {
170 if i != j {
171 L[[i, j]] = -W[[i, j]];
172 }
173 }
174 }
175
176 L
177 }
178
179 fn soft_threshold(&self, x: f64, threshold: f64) -> f64 {
180 if x > threshold {
181 x - threshold
182 } else if x < -threshold {
183 x + threshold
184 } else {
185 0.0
186 }
187 }
188
189 #[allow(non_snake_case)] fn proximal_gradient_step(&self, W: &Array2<f64>, grad: &Array2<f64>, lr: f64) -> Array2<f64> {
191 let mut W_new = W - lr * grad;
192
193 let threshold = lr * self.lambda_sparse;
195 W_new.mapv_inplace(|x| self.soft_threshold(x, threshold));
196
197 W_new.mapv_inplace(|x| x.max(0.0));
199
200 if self.enforce_symmetry {
202 let n = W_new.nrows();
203 for i in 0..n {
204 for j in 0..n {
205 if i != j {
206 let avg = (W_new[[i, j]] + W_new[[j, i]]) / 2.0;
207 W_new[[i, j]] = avg;
208 W_new[[j, i]] = avg;
209 }
210 }
211 }
212 }
213
214 for i in 0..W_new.nrows() {
216 W_new[[i, i]] = 0.0;
217 }
218
219 W_new
220 }
221
222 #[allow(non_snake_case)] fn normalize_graph(&self, W: &Array2<f64>) -> Array2<f64> {
224 if !self.normalize_weights {
225 return W.clone();
226 }
227
228 let mut W_norm = W.clone();
229 let n_samples = W.nrows();
230
231 for i in 0..n_samples {
233 let row_sum: f64 = W.row(i).sum();
234 if row_sum > 0.0 {
235 for j in 0..n_samples {
236 W_norm[[i, j]] = W[[i, j]] / row_sum;
237 }
238 }
239 }
240
241 W_norm
242 }
243
244 #[allow(non_snake_case)]
245 fn propagate_labels(&self, W: &Array2<f64>, Y_init: &Array2<f64>) -> SklResult<Array2<f64>> {
246 let n_samples = W.nrows();
247
248 let D = W.sum_axis(Axis(1));
250 let mut P = Array2::zeros((n_samples, n_samples));
251 for i in 0..n_samples {
252 if D[i] > 0.0 {
253 for j in 0..n_samples {
254 P[[i, j]] = W[[i, j]] / D[i];
255 }
256 }
257 }
258
259 let mut Y = Y_init.clone();
260 let Y_static = Y_init.clone();
261
262 for _iter in 0..50 {
264 let prev_Y = Y.clone();
265 Y = 0.9 * P.dot(&Y) + 0.1 * &Y_static;
266
267 let diff = (&Y - &prev_Y).mapv(|x| x.abs()).sum();
269 if diff < 1e-6 {
270 break;
271 }
272 }
273
274 Ok(Y)
275 }
276}
277
278impl Default for GraphStructureLearning<Untrained> {
279 fn default() -> Self {
280 Self::new()
281 }
282}
283
284impl Estimator for GraphStructureLearning<Untrained> {
285 type Config = ();
286 type Error = SklearsError;
287 type Float = Float;
288
289 fn config(&self) -> &Self::Config {
290 &()
291 }
292}
293
294impl Fit<ArrayView2<'_, Float>, ArrayView1<'_, i32>> for GraphStructureLearning<Untrained> {
295 type Fitted = GraphStructureLearning<GraphStructureLearningTrained>;
296
297 #[allow(non_snake_case)]
298 fn fit(self, X: &ArrayView2<'_, Float>, y: &ArrayView1<'_, i32>) -> SklResult<Self::Fitted> {
299 let X = X.to_owned();
300 let y = y.to_owned();
301
302 let (n_samples, _n_features) = X.dim();
303
304 let mut labeled_indices = Vec::new();
306 let mut classes = std::collections::HashSet::new();
307
308 for (i, &label) in y.iter().enumerate() {
309 if label != -1 {
310 labeled_indices.push(i);
311 classes.insert(label);
312 }
313 }
314
315 if labeled_indices.is_empty() {
316 return Err(SklearsError::InvalidInput(
317 "No labeled samples provided".to_string(),
318 ));
319 }
320
321 let classes: Vec<i32> = classes.into_iter().collect();
322 let n_classes = classes.len();
323
324 let mut W = self.initialize_graph(&X);
326
327 let mut Y = Array2::zeros((n_samples, n_classes));
329 for &idx in &labeled_indices {
330 if let Some(class_idx) = classes.iter().position(|&c| c == y[idx]) {
331 Y[[idx, class_idx]] = 1.0;
332 }
333 }
334
335 let Y_init = Y.clone();
336 let mut lr = self.learning_rate;
337 let mut prev_loss = f64::INFINITY;
338
339 for iteration in 0..self.max_iter {
341 Y = self.propagate_labels(&W, &Y_init)?;
343
344 let L = self.compute_laplacian(&W);
346
347 let WX = W.dot(&X);
349 let data_fidelity = (&X - &WX).mapv(|x| x * x).sum();
350
351 let smoothness_loss = {
353 let LY = L.dot(&Y);
354 let mut trace = 0.0;
355 for i in 0..n_samples {
356 for j in 0..n_classes {
357 trace += Y[[i, j]] * LY[[i, j]];
358 }
359 }
360 trace
361 };
362
363 let sparsity_loss = W.iter().map(|&x| x.abs()).sum::<f64>();
365
366 let total_loss = data_fidelity
368 + self.beta_smoothness * smoothness_loss
369 + self.lambda_sparse * sparsity_loss;
370
371 if (prev_loss - total_loss).abs() < self.tol {
373 break;
374 }
375
376 if self.adaptive_lr {
378 if total_loss > prev_loss {
379 lr *= 0.8; } else if iteration % 10 == 0 && total_loss < prev_loss {
381 lr *= 1.1; }
383 lr = lr.clamp(1e-6, 0.1); }
385
386 prev_loss = total_loss;
387
388 let residual = &WX - &X;
391 let mut grad_W = 2.0 * residual.dot(&X.t());
392
393 let YYT = Y.dot(&Y.t());
395 let D = W.sum_axis(Axis(1));
396 for i in 0..n_samples {
397 for j in 0..n_samples {
398 if i == j {
399 grad_W[[i, j]] +=
400 self.beta_smoothness * (D[i] * YYT[[i, j]] - W[[i, j]] * YYT[[i, j]]);
401 } else {
402 grad_W[[i, j]] += self.beta_smoothness * (-YYT[[i, j]]);
403 }
404 }
405 }
406
407 W = self.proximal_gradient_step(&W, &grad_W, lr);
409 }
410
411 let W_final = self.normalize_graph(&W);
413
414 let Y_final = self.propagate_labels(&W_final, &Y_init)?;
416
417 Ok(GraphStructureLearning {
418 state: GraphStructureLearningTrained {
419 X_train: X,
420 y_train: y,
421 classes: Array1::from(classes),
422 learned_graph: W_final,
423 label_distributions: Y_final,
424 },
425 lambda_sparse: self.lambda_sparse,
426 beta_smoothness: self.beta_smoothness,
427 max_iter: self.max_iter,
428 tol: self.tol,
429 learning_rate: self.learning_rate,
430 adaptive_lr: self.adaptive_lr,
431 enforce_symmetry: self.enforce_symmetry,
432 normalize_weights: self.normalize_weights,
433 })
434 }
435}
436
437impl Predict<ArrayView2<'_, Float>, Array1<i32>>
438 for GraphStructureLearning<GraphStructureLearningTrained>
439{
440 #[allow(non_snake_case)]
441 fn predict(&self, X: &ArrayView2<'_, Float>) -> SklResult<Array1<i32>> {
442 let X = X.to_owned();
443 let n_test = X.nrows();
444 let mut predictions = Array1::zeros(n_test);
445
446 for i in 0..n_test {
447 let mut min_dist = f64::INFINITY;
449 let mut best_idx = 0;
450
451 for j in 0..self.state.X_train.nrows() {
452 let diff = &X.row(i) - &self.state.X_train.row(j);
453 let dist = diff.mapv(|x: f64| x * x).sum().sqrt();
454 if dist < min_dist {
455 min_dist = dist;
456 best_idx = j;
457 }
458 }
459
460 let distributions = self.state.label_distributions.row(best_idx);
462 let max_idx = distributions
463 .iter()
464 .enumerate()
465 .max_by(|a, b| a.1.partial_cmp(b.1).expect("operation should succeed"))
466 .expect("operation should succeed")
467 .0;
468
469 predictions[i] = self.state.classes[max_idx];
470 }
471
472 Ok(predictions)
473 }
474}
475
476impl PredictProba<ArrayView2<'_, Float>, Array2<f64>>
477 for GraphStructureLearning<GraphStructureLearningTrained>
478{
479 #[allow(non_snake_case)]
480 fn predict_proba(&self, X: &ArrayView2<'_, Float>) -> SklResult<Array2<f64>> {
481 let X = X.to_owned();
482 let n_test = X.nrows();
483 let n_classes = self.state.classes.len();
484 let mut probas = Array2::zeros((n_test, n_classes));
485
486 for i in 0..n_test {
487 let mut min_dist = f64::INFINITY;
489 let mut best_idx = 0;
490
491 for j in 0..self.state.X_train.nrows() {
492 let diff = &X.row(i) - &self.state.X_train.row(j);
493 let dist = diff.mapv(|x: f64| x * x).sum().sqrt();
494 if dist < min_dist {
495 min_dist = dist;
496 best_idx = j;
497 }
498 }
499
500 for k in 0..n_classes {
502 probas[[i, k]] = self.state.label_distributions[[best_idx, k]];
503 }
504 }
505
506 Ok(probas)
507 }
508}
509
510#[derive(Debug, Clone)]
542pub struct RobustGraphLearning<S = Untrained> {
543 state: S,
544 lambda_sparse: f64,
545 lambda_robust: f64,
546 max_iter: usize,
547 tol: f64,
548 robust_metric: String,
549 huber_delta: f64,
550 tukey_c: f64,
551}
552
553impl RobustGraphLearning<Untrained> {
554 pub fn new() -> Self {
556 Self {
557 state: Untrained,
558 lambda_sparse: 0.1,
559 lambda_robust: 1.0,
560 max_iter: 100,
561 tol: 1e-4,
562 robust_metric: "huber".to_string(),
563 huber_delta: 1.0,
564 tukey_c: 4.685,
565 }
566 }
567
568 pub fn lambda_sparse(mut self, lambda_sparse: f64) -> Self {
570 self.lambda_sparse = lambda_sparse;
571 self
572 }
573
574 pub fn lambda_robust(mut self, lambda_robust: f64) -> Self {
576 self.lambda_robust = lambda_robust;
577 self
578 }
579
580 pub fn max_iter(mut self, max_iter: usize) -> Self {
582 self.max_iter = max_iter;
583 self
584 }
585
586 pub fn tol(mut self, tol: f64) -> Self {
588 self.tol = tol;
589 self
590 }
591
592 pub fn robust_metric(mut self, metric: String) -> Self {
594 self.robust_metric = metric;
595 self
596 }
597
598 pub fn huber_delta(mut self, delta: f64) -> Self {
600 self.huber_delta = delta;
601 self
602 }
603
604 pub fn tukey_c(mut self, c: f64) -> Self {
606 self.tukey_c = c;
607 self
608 }
609
610 fn robust_distance(&self, x1: &ArrayView1<f64>, x2: &ArrayView1<f64>) -> f64 {
611 let diff = x1 - x2;
612
613 match self.robust_metric.as_str() {
614 "l1" => diff.mapv(|x| x.abs()).sum(),
615 "huber" => diff
616 .mapv(|x| {
617 let abs_x = x.abs();
618 if abs_x <= self.huber_delta {
619 0.5 * x * x
620 } else {
621 self.huber_delta * (abs_x - 0.5 * self.huber_delta)
622 }
623 })
624 .sum(),
625 "tukey" => diff
626 .mapv(|x| {
627 let abs_x = x.abs();
628 if abs_x <= self.tukey_c {
629 let ratio = x / self.tukey_c;
630 (self.tukey_c * self.tukey_c / 6.0) * (1.0 - (1.0 - ratio * ratio).powi(3))
631 } else {
632 self.tukey_c * self.tukey_c / 6.0
633 }
634 })
635 .sum(),
636 _ => diff.mapv(|x| x * x).sum().sqrt(), }
638 }
639
640 #[allow(non_snake_case)] fn compute_robust_weights(&self, X: &Array2<f64>) -> Array2<f64> {
642 let n_samples = X.nrows();
643 let mut W = Array2::zeros((n_samples, n_samples));
644
645 for i in 0..n_samples {
646 for j in 0..n_samples {
647 if i != j {
648 let dist = self.robust_distance(&X.row(i), &X.row(j));
649 W[[i, j]] = (-dist / (2.0 * 1.0_f64.powi(2))).exp();
650 }
651 }
652 }
653
654 W
655 }
656}
657
658impl Default for RobustGraphLearning<Untrained> {
659 fn default() -> Self {
660 Self::new()
661 }
662}
663
664impl Estimator for RobustGraphLearning<Untrained> {
665 type Config = ();
666 type Error = SklearsError;
667 type Float = Float;
668
669 fn config(&self) -> &Self::Config {
670 &()
671 }
672}
673
674impl Fit<ArrayView2<'_, Float>, ArrayView1<'_, i32>> for RobustGraphLearning<Untrained> {
675 type Fitted = RobustGraphLearning<RobustGraphLearningTrained>;
676
677 #[allow(non_snake_case)]
678 fn fit(self, X: &ArrayView2<'_, Float>, y: &ArrayView1<'_, i32>) -> SklResult<Self::Fitted> {
679 let X = X.to_owned();
680 let y = y.to_owned();
681
682 let mut labeled_indices = Vec::new();
684 let mut classes = std::collections::HashSet::new();
685
686 for (i, &label) in y.iter().enumerate() {
687 if label != -1 {
688 labeled_indices.push(i);
689 classes.insert(label);
690 }
691 }
692
693 if labeled_indices.is_empty() {
694 return Err(SklearsError::InvalidInput(
695 "No labeled samples provided".to_string(),
696 ));
697 }
698
699 let classes: Vec<i32> = classes.into_iter().collect();
700 let n_classes = classes.len();
701 let n_samples = X.nrows();
702
703 let mut W = self.compute_robust_weights(&X);
705
706 let threshold = self.lambda_sparse;
708 W.mapv_inplace(|x| if x > threshold { x - threshold } else { 0.0 });
709
710 W.mapv_inplace(|x| x.max(0.0));
712 for i in 0..n_samples {
713 W[[i, i]] = 0.0;
714 }
715
716 for i in 0..n_samples {
718 for j in i + 1..n_samples {
719 let avg = (W[[i, j]] + W[[j, i]]) / 2.0;
720 W[[i, j]] = avg;
721 W[[j, i]] = avg;
722 }
723 }
724
725 let mut Y = Array2::zeros((n_samples, n_classes));
727 for &idx in &labeled_indices {
728 if let Some(class_idx) = classes.iter().position(|&c| c == y[idx]) {
729 Y[[idx, class_idx]] = 1.0;
730 }
731 }
732
733 let D = W.sum_axis(Axis(1));
735 let mut P = Array2::zeros((n_samples, n_samples));
736 for i in 0..n_samples {
737 if D[i] > 0.0 {
738 for j in 0..n_samples {
739 P[[i, j]] = W[[i, j]] / D[i];
740 }
741 }
742 }
743
744 let Y_static = Y.clone();
745
746 for _iter in 0..50 {
748 let prev_Y = Y.clone();
749 Y = 0.9 * P.dot(&Y) + 0.1 * &Y_static;
750
751 let diff = (&Y - &prev_Y).mapv(|x| x.abs()).sum();
753 if diff < 1e-6 {
754 break;
755 }
756 }
757
758 Ok(RobustGraphLearning {
759 state: RobustGraphLearningTrained {
760 X_train: X,
761 y_train: y,
762 classes: Array1::from(classes),
763 learned_graph: W,
764 label_distributions: Y,
765 },
766 lambda_sparse: self.lambda_sparse,
767 lambda_robust: self.lambda_robust,
768 max_iter: self.max_iter,
769 tol: self.tol,
770 robust_metric: self.robust_metric,
771 huber_delta: self.huber_delta,
772 tukey_c: self.tukey_c,
773 })
774 }
775}
776
777impl Predict<ArrayView2<'_, Float>, Array1<i32>>
778 for RobustGraphLearning<RobustGraphLearningTrained>
779{
780 #[allow(non_snake_case)]
781 fn predict(&self, X: &ArrayView2<'_, Float>) -> SklResult<Array1<i32>> {
782 let X = X.to_owned();
783 let n_test = X.nrows();
784 let mut predictions = Array1::zeros(n_test);
785
786 for i in 0..n_test {
787 let mut min_dist = f64::INFINITY;
789 let mut best_idx = 0;
790
791 for j in 0..self.state.X_train.nrows() {
792 let diff = &X.row(i) - &self.state.X_train.row(j);
793 let dist = diff.mapv(|x: f64| x * x).sum().sqrt();
794 if dist < min_dist {
795 min_dist = dist;
796 best_idx = j;
797 }
798 }
799
800 let distributions = self.state.label_distributions.row(best_idx);
802 let max_idx = distributions
803 .iter()
804 .enumerate()
805 .max_by(|a, b| a.1.partial_cmp(b.1).expect("operation should succeed"))
806 .expect("operation should succeed")
807 .0;
808
809 predictions[i] = self.state.classes[max_idx];
810 }
811
812 Ok(predictions)
813 }
814}
815
816impl PredictProba<ArrayView2<'_, Float>, Array2<f64>>
817 for RobustGraphLearning<RobustGraphLearningTrained>
818{
819 #[allow(non_snake_case)]
820 fn predict_proba(&self, X: &ArrayView2<'_, Float>) -> SklResult<Array2<f64>> {
821 let X = X.to_owned();
822 let n_test = X.nrows();
823 let n_classes = self.state.classes.len();
824 let mut probas = Array2::zeros((n_test, n_classes));
825
826 for i in 0..n_test {
827 let mut min_dist = f64::INFINITY;
829 let mut best_idx = 0;
830
831 for j in 0..self.state.X_train.nrows() {
832 let diff = &X.row(i) - &self.state.X_train.row(j);
833 let dist = diff.mapv(|x: f64| x * x).sum().sqrt();
834 if dist < min_dist {
835 min_dist = dist;
836 best_idx = j;
837 }
838 }
839
840 for k in 0..n_classes {
842 probas[[i, k]] = self.state.label_distributions[[best_idx, k]];
843 }
844 }
845
846 Ok(probas)
847 }
848}
849
850#[derive(Debug, Clone)]
852#[allow(non_snake_case)] pub struct GraphStructureLearningTrained {
854 pub X_train: Array2<f64>,
856 pub y_train: Array1<i32>,
858 pub classes: Array1<i32>,
860 pub learned_graph: Array2<f64>,
862 pub label_distributions: Array2<f64>,
864}
865
866#[derive(Debug, Clone)]
868#[allow(non_snake_case)] pub struct RobustGraphLearningTrained {
870 pub X_train: Array2<f64>,
872 pub y_train: Array1<i32>,
874 pub classes: Array1<i32>,
876 pub learned_graph: Array2<f64>,
878 pub label_distributions: Array2<f64>,
880}
881
882#[derive(Debug, Clone)]
917pub struct DistributedGraphLearning<S = Untrained> {
918 state: S,
919 n_workers: usize,
920 lambda_sparse: f64,
921 beta_smoothness: f64,
922 max_iter: usize,
923 tol: f64,
924 learning_rate: f64,
925 partition_strategy: String,
926 communication_rounds: usize,
927 overlap_ratio: f64,
928 consensus_weight: f64,
929}
930
931impl DistributedGraphLearning<Untrained> {
932 pub fn new() -> Self {
934 Self {
935 state: Untrained,
936 n_workers: 2,
937 lambda_sparse: 0.1,
938 beta_smoothness: 1.0,
939 max_iter: 100,
940 tol: 1e-4,
941 learning_rate: 0.01,
942 partition_strategy: "spectral".to_string(),
943 communication_rounds: 10,
944 overlap_ratio: 0.1,
945 consensus_weight: 0.5,
946 }
947 }
948
949 pub fn n_workers(mut self, n_workers: usize) -> Self {
951 self.n_workers = n_workers;
952 self
953 }
954
955 pub fn lambda_sparse(mut self, lambda_sparse: f64) -> Self {
957 self.lambda_sparse = lambda_sparse;
958 self
959 }
960
961 pub fn beta_smoothness(mut self, beta_smoothness: f64) -> Self {
963 self.beta_smoothness = beta_smoothness;
964 self
965 }
966
967 pub fn max_iter(mut self, max_iter: usize) -> Self {
969 self.max_iter = max_iter;
970 self
971 }
972
973 pub fn tol(mut self, tol: f64) -> Self {
975 self.tol = tol;
976 self
977 }
978
979 pub fn learning_rate(mut self, learning_rate: f64) -> Self {
981 self.learning_rate = learning_rate;
982 self
983 }
984
985 pub fn partition_strategy(mut self, strategy: String) -> Self {
987 self.partition_strategy = strategy;
988 self
989 }
990
991 pub fn communication_rounds(mut self, rounds: usize) -> Self {
993 self.communication_rounds = rounds;
994 self
995 }
996
997 pub fn overlap_ratio(mut self, ratio: f64) -> Self {
999 self.overlap_ratio = ratio;
1000 self
1001 }
1002
1003 pub fn consensus_weight(mut self, weight: f64) -> Self {
1005 self.consensus_weight = weight;
1006 self
1007 }
1008
1009 fn partition_nodes(&self, n_samples: usize) -> Vec<Vec<usize>> {
1010 let nodes_per_worker = n_samples.div_ceil(self.n_workers);
1011 let overlap_size = (nodes_per_worker as f64 * self.overlap_ratio) as usize;
1012
1013 let mut partitions = Vec::with_capacity(self.n_workers);
1014
1015 match self.partition_strategy.as_str() {
1016 "random" => {
1017 let mut nodes: Vec<usize> = (0..n_samples).collect();
1019 use scirs2_core::random::rand_prelude::SliceRandom;
1020 let mut rng = Random::seed(42);
1021 nodes.shuffle(&mut rng);
1022
1023 for i in 0..self.n_workers {
1024 let start = i * nodes_per_worker;
1025 let end = ((i + 1) * nodes_per_worker).min(n_samples);
1026 let overlap_start = start.saturating_sub(overlap_size);
1027 let overlap_end = (end + overlap_size).min(n_samples);
1028
1029 let mut partition = Vec::new();
1030 for j in overlap_start..overlap_end {
1031 if j < nodes.len() {
1032 partition.push(nodes[j]);
1033 }
1034 }
1035 partitions.push(partition);
1036 }
1037 }
1038 "spectral" => {
1039 self.spectral_partition(n_samples, &mut partitions, nodes_per_worker, overlap_size);
1041 }
1042 _ => {
1043 for i in 0..self.n_workers {
1045 let start = i * nodes_per_worker;
1046 let end = ((i + 1) * nodes_per_worker).min(n_samples);
1047 let overlap_start = start.saturating_sub(overlap_size);
1048 let overlap_end = (end + overlap_size).min(n_samples);
1049
1050 let partition: Vec<usize> = (overlap_start..overlap_end).collect();
1051 partitions.push(partition);
1052 }
1053 }
1054 }
1055
1056 partitions
1057 }
1058
1059 fn spectral_partition(
1060 &self,
1061 n_samples: usize,
1062 partitions: &mut Vec<Vec<usize>>,
1063 nodes_per_worker: usize,
1064 overlap_size: usize,
1065 ) {
1066 let mut spectral_order: Vec<usize> = (0..n_samples).collect();
1069
1070 let center = n_samples / 2;
1072 spectral_order.sort_by_key(|&i| i.abs_diff(center));
1073
1074 for i in 0..self.n_workers {
1075 let start = i * nodes_per_worker;
1076 let end = ((i + 1) * nodes_per_worker).min(n_samples);
1077 let overlap_start = start.saturating_sub(overlap_size);
1078 let overlap_end = (end + overlap_size).min(n_samples);
1079
1080 let mut partition = Vec::new();
1081 for j in overlap_start..overlap_end {
1082 if j < spectral_order.len() {
1083 partition.push(spectral_order[j]);
1084 }
1085 }
1086 partitions.push(partition);
1087 }
1088 }
1089
1090 #[allow(non_snake_case)] fn extract_subgraph(&self, X: &Array2<f64>, partition: &[usize]) -> Array2<f64> {
1092 let n_nodes = partition.len();
1093 let n_features = X.ncols();
1094 let mut X_sub = Array2::zeros((n_nodes, n_features));
1095
1096 for (i, &node_idx) in partition.iter().enumerate() {
1097 if node_idx < X.nrows() {
1098 X_sub.row_mut(i).assign(&X.row(node_idx));
1099 }
1100 }
1101
1102 X_sub
1103 }
1104
1105 fn extract_sublabels(&self, y: &Array1<i32>, partition: &[usize]) -> Array1<i32> {
1106 let n_nodes = partition.len();
1107 let mut y_sub = Array1::from_elem(n_nodes, -1);
1108
1109 for (i, &node_idx) in partition.iter().enumerate() {
1110 if node_idx < y.len() {
1111 y_sub[i] = y[node_idx];
1112 }
1113 }
1114
1115 y_sub
1116 }
1117
1118 #[allow(non_snake_case)] fn learn_local_graph(
1120 &self,
1121 X_sub: &Array2<f64>,
1122 _y_sub: &Array1<i32>,
1123 ) -> SklResult<Array2<f64>> {
1124 let n_samples = X_sub.nrows();
1125 let mut W = Array2::zeros((n_samples, n_samples));
1126
1127 let k = (n_samples as f64).sqrt().ceil() as usize;
1129 let k = k.clamp(3, 10);
1130
1131 for i in 0..n_samples {
1132 let mut distances: Vec<(usize, f64)> = Vec::new();
1133 for j in 0..n_samples {
1134 if i != j {
1135 let diff = &X_sub.row(i) - &X_sub.row(j);
1136 let dist = diff.mapv(|x: f64| x * x).sum().sqrt();
1137 distances.push((j, dist));
1138 }
1139 }
1140
1141 distances.sort_by(|a, b| a.1.partial_cmp(&b.1).expect("operation should succeed"));
1142
1143 for &(j, dist) in distances.iter().take(k) {
1144 let weight = (-dist / (2.0 * 1.0_f64.powi(2))).exp();
1145 W[[i, j]] = weight;
1146 W[[j, i]] = weight; }
1148 }
1149
1150 let threshold = self.lambda_sparse;
1152 W.mapv_inplace(|x| if x > threshold { x - threshold } else { 0.0 });
1153 W.mapv_inplace(|x| x.max(0.0));
1154
1155 for i in 0..n_samples {
1157 W[[i, i]] = 0.0;
1158 }
1159
1160 Ok(W)
1161 }
1162
1163 fn communicate_boundaries(
1164 &self,
1165 local_graphs: &[Array2<f64>],
1166 partitions: &[Vec<usize>],
1167 ) -> Vec<Array2<f64>> {
1168 let mut updated_graphs = local_graphs.to_vec();
1169
1170 for i in 0..self.n_workers {
1172 for j in (i + 1)..self.n_workers {
1173 let common_nodes: Vec<(usize, usize)> = partitions[i]
1175 .iter()
1176 .enumerate()
1177 .filter_map(|(idx_i, &node)| {
1178 partitions[j]
1179 .iter()
1180 .position(|&n| n == node)
1181 .map(|idx_j| (idx_i, idx_j))
1182 })
1183 .collect();
1184
1185 for &(idx_i, idx_j) in &common_nodes {
1187 if idx_i < updated_graphs[i].nrows() && idx_j < updated_graphs[j].nrows() {
1188 for &(other_i, other_j) in &common_nodes {
1189 if other_i < updated_graphs[i].ncols()
1190 && other_j < updated_graphs[j].ncols()
1191 {
1192 let weight_i = updated_graphs[i][[idx_i, other_i]];
1193 let weight_j = updated_graphs[j][[idx_j, other_j]];
1194 let avg_weight = (weight_i + weight_j) / 2.0;
1195
1196 updated_graphs[i][[idx_i, other_i]] = avg_weight;
1197 updated_graphs[j][[idx_j, other_j]] = avg_weight;
1198 }
1199 }
1200 }
1201 }
1202 }
1203 }
1204
1205 updated_graphs
1206 }
1207
1208 fn merge_graphs(
1209 &self,
1210 local_graphs: &[Array2<f64>],
1211 partitions: &[Vec<usize>],
1212 n_total: usize,
1213 ) -> Array2<f64> {
1214 let mut global_graph = Array2::zeros((n_total, n_total));
1215 let mut weight_counts: Array2<f64> = Array2::zeros((n_total, n_total));
1216
1217 for (local_graph, partition) in local_graphs.iter().zip(partitions.iter()) {
1219 for (i, &node_i) in partition.iter().enumerate() {
1220 for (j, &node_j) in partition.iter().enumerate() {
1221 if i < local_graph.nrows()
1222 && j < local_graph.ncols()
1223 && node_i < n_total
1224 && node_j < n_total
1225 {
1226 global_graph[[node_i, node_j]] += local_graph[[i, j]];
1227 if local_graph[[i, j]] > 0.0 {
1228 weight_counts[[node_i, node_j]] += 1.0;
1229 }
1230 }
1231 }
1232 }
1233 }
1234
1235 for i in 0..n_total {
1237 for j in 0..n_total {
1238 if weight_counts[[i, j]] > 0.0 {
1239 global_graph[[i, j]] /= weight_counts[[i, j]];
1240 }
1241 }
1242 }
1243
1244 global_graph
1245 }
1246}
1247
1248impl Default for DistributedGraphLearning<Untrained> {
1249 fn default() -> Self {
1250 Self::new()
1251 }
1252}
1253
1254impl Estimator for DistributedGraphLearning<Untrained> {
1255 type Config = ();
1256 type Error = SklearsError;
1257 type Float = Float;
1258
1259 fn config(&self) -> &Self::Config {
1260 &()
1261 }
1262}
1263
1264impl Fit<ArrayView2<'_, Float>, ArrayView1<'_, i32>> for DistributedGraphLearning<Untrained> {
1265 type Fitted = DistributedGraphLearning<DistributedGraphLearningTrained>;
1266
1267 #[allow(non_snake_case)]
1268 fn fit(self, X: &ArrayView2<'_, Float>, y: &ArrayView1<'_, i32>) -> SklResult<Self::Fitted> {
1269 let X = X.to_owned();
1270 let y = y.to_owned();
1271 let (n_samples, _n_features) = X.dim();
1272
1273 let mut labeled_indices = Vec::new();
1275 let mut classes = std::collections::HashSet::new();
1276
1277 for (i, &label) in y.iter().enumerate() {
1278 if label != -1 {
1279 labeled_indices.push(i);
1280 classes.insert(label);
1281 }
1282 }
1283
1284 if labeled_indices.is_empty() {
1285 return Err(SklearsError::InvalidInput(
1286 "No labeled samples provided".to_string(),
1287 ));
1288 }
1289
1290 let classes: Vec<i32> = classes.into_iter().collect();
1291
1292 let partitions = self.partition_nodes(n_samples);
1294
1295 let mut local_graphs = Vec::with_capacity(self.n_workers);
1297 for partition in &partitions {
1298 let X_sub = self.extract_subgraph(&X, partition);
1299 let y_sub = self.extract_sublabels(&y, partition);
1300 let local_graph = self.learn_local_graph(&X_sub, &y_sub)?;
1301 local_graphs.push(local_graph);
1302 }
1303
1304 for _round in 0..self.communication_rounds {
1306 local_graphs = self.communicate_boundaries(&local_graphs, &partitions);
1307 }
1308
1309 let global_graph = self.merge_graphs(&local_graphs, &partitions, n_samples);
1311
1312 let n_classes = classes.len();
1314 let mut Y = Array2::zeros((n_samples, n_classes));
1315 for &idx in &labeled_indices {
1316 if let Some(class_idx) = classes.iter().position(|&c| c == y[idx]) {
1317 Y[[idx, class_idx]] = 1.0;
1318 }
1319 }
1320
1321 let D = global_graph.sum_axis(Axis(1));
1323 let mut P = Array2::zeros((n_samples, n_samples));
1324 for i in 0..n_samples {
1325 if D[i] > 0.0 {
1326 for j in 0..n_samples {
1327 P[[i, j]] = global_graph[[i, j]] / D[i];
1328 }
1329 }
1330 }
1331
1332 let Y_static = Y.clone();
1333 for _iter in 0..50 {
1334 let prev_Y = Y.clone();
1335 Y = 0.9 * P.dot(&Y) + 0.1 * &Y_static;
1336
1337 let diff = (&Y - &prev_Y).mapv(|x| x.abs()).sum();
1338 if diff < 1e-6 {
1339 break;
1340 }
1341 }
1342
1343 Ok(DistributedGraphLearning {
1344 state: DistributedGraphLearningTrained {
1345 X_train: X,
1346 y_train: y,
1347 classes: Array1::from(classes),
1348 global_graph,
1349 label_distributions: Y,
1350 partitions,
1351 },
1352 n_workers: self.n_workers,
1353 lambda_sparse: self.lambda_sparse,
1354 beta_smoothness: self.beta_smoothness,
1355 max_iter: self.max_iter,
1356 tol: self.tol,
1357 learning_rate: self.learning_rate,
1358 partition_strategy: self.partition_strategy,
1359 communication_rounds: self.communication_rounds,
1360 overlap_ratio: self.overlap_ratio,
1361 consensus_weight: self.consensus_weight,
1362 })
1363 }
1364}
1365
1366impl Predict<ArrayView2<'_, Float>, Array1<i32>>
1367 for DistributedGraphLearning<DistributedGraphLearningTrained>
1368{
1369 #[allow(non_snake_case)]
1370 fn predict(&self, X: &ArrayView2<'_, Float>) -> SklResult<Array1<i32>> {
1371 let X = X.to_owned();
1372 let n_test = X.nrows();
1373 let mut predictions = Array1::zeros(n_test);
1374
1375 for i in 0..n_test {
1376 let mut min_dist = f64::INFINITY;
1377 let mut best_idx = 0;
1378
1379 for j in 0..self.state.X_train.nrows() {
1380 let diff = &X.row(i) - &self.state.X_train.row(j);
1381 let dist = diff.mapv(|x: f64| x * x).sum().sqrt();
1382 if dist < min_dist {
1383 min_dist = dist;
1384 best_idx = j;
1385 }
1386 }
1387
1388 let distributions = self.state.label_distributions.row(best_idx);
1389 let max_idx = distributions
1390 .iter()
1391 .enumerate()
1392 .max_by(|a, b| a.1.partial_cmp(b.1).expect("operation should succeed"))
1393 .expect("operation should succeed")
1394 .0;
1395
1396 predictions[i] = self.state.classes[max_idx];
1397 }
1398
1399 Ok(predictions)
1400 }
1401}
1402
1403impl PredictProba<ArrayView2<'_, Float>, Array2<f64>>
1404 for DistributedGraphLearning<DistributedGraphLearningTrained>
1405{
1406 #[allow(non_snake_case)]
1407 fn predict_proba(&self, X: &ArrayView2<'_, Float>) -> SklResult<Array2<f64>> {
1408 let X = X.to_owned();
1409 let n_test = X.nrows();
1410 let n_classes = self.state.classes.len();
1411 let mut probas = Array2::zeros((n_test, n_classes));
1412
1413 for i in 0..n_test {
1414 let mut min_dist = f64::INFINITY;
1415 let mut best_idx = 0;
1416
1417 for j in 0..self.state.X_train.nrows() {
1418 let diff = &X.row(i) - &self.state.X_train.row(j);
1419 let dist = diff.mapv(|x: f64| x * x).sum().sqrt();
1420 if dist < min_dist {
1421 min_dist = dist;
1422 best_idx = j;
1423 }
1424 }
1425
1426 for k in 0..n_classes {
1427 probas[[i, k]] = self.state.label_distributions[[best_idx, k]];
1428 }
1429 }
1430
1431 Ok(probas)
1432 }
1433}
1434
1435#[derive(Debug, Clone)]
1437#[allow(non_snake_case)] pub struct DistributedGraphLearningTrained {
1439 pub X_train: Array2<f64>,
1441 pub y_train: Array1<i32>,
1443 pub classes: Array1<i32>,
1445 pub global_graph: Array2<f64>,
1447 pub label_distributions: Array2<f64>,
1449 pub partitions: Vec<Vec<usize>>,
1451}
1452
1453#[allow(non_snake_case)]
1454#[cfg(test)]
1455mod tests {
1456 use super::*;
1457 use scirs2_core::array;
1458
1459 #[test]
1460 #[allow(non_snake_case)]
1461 fn test_graph_structure_learning() {
1462 let X = array![[1.0, 2.0], [2.0, 3.0], [3.0, 4.0], [4.0, 5.0]];
1463 let y = array![0, 1, -1, -1]; let gsl = GraphStructureLearning::new()
1466 .lambda_sparse(0.1)
1467 .beta_smoothness(1.0)
1468 .max_iter(20);
1469 let fitted = gsl
1470 .fit(&X.view(), &y.view())
1471 .expect("operation should succeed");
1472
1473 let predictions = fitted.predict(&X.view()).expect("operation should succeed");
1474 assert_eq!(predictions.len(), 4);
1475
1476 let probas = fitted
1477 .predict_proba(&X.view())
1478 .expect("operation should succeed");
1479 assert_eq!(probas.dim(), (4, 2));
1480
1481 let n_edges = fitted
1483 .state
1484 .learned_graph
1485 .iter()
1486 .filter(|&&x| x > 0.0)
1487 .count();
1488 let total_edges = 4 * 4 - 4; assert!(n_edges < total_edges); }
1491
1492 #[test]
1493 #[allow(non_snake_case)]
1494 fn test_robust_graph_learning() {
1495 let X = array![[1.0, 2.0], [2.0, 3.0], [3.0, 4.0], [4.0, 5.0]];
1496 let y = array![0, 1, -1, -1]; let rgl = RobustGraphLearning::new()
1499 .lambda_sparse(0.1)
1500 .robust_metric("huber".to_string())
1501 .max_iter(20);
1502 let fitted = rgl
1503 .fit(&X.view(), &y.view())
1504 .expect("operation should succeed");
1505
1506 let predictions = fitted.predict(&X.view()).expect("operation should succeed");
1507 assert_eq!(predictions.len(), 4);
1508
1509 let probas = fitted
1510 .predict_proba(&X.view())
1511 .expect("operation should succeed");
1512 assert_eq!(probas.dim(), (4, 2));
1513 }
1514
1515 #[test]
1516 fn test_robust_distance_metrics() {
1517 let rgl = RobustGraphLearning::new();
1518 let x1 = array![1.0, 2.0];
1519 let x2 = array![3.0, 4.0];
1520
1521 let rgl_l1 = rgl.clone().robust_metric("l1".to_string());
1523 let dist_l1 = rgl_l1.robust_distance(&x1.view(), &x2.view());
1524 assert_eq!(dist_l1, 4.0); let rgl_huber = rgl
1528 .clone()
1529 .robust_metric("huber".to_string())
1530 .huber_delta(1.0);
1531 let dist_huber = rgl_huber.robust_distance(&x1.view(), &x2.view());
1532 assert!(dist_huber > 0.0);
1533
1534 let rgl_tukey = rgl.robust_metric("tukey".to_string());
1536 let dist_tukey = rgl_tukey.robust_distance(&x1.view(), &x2.view());
1537 assert!(dist_tukey > 0.0);
1538 }
1539
1540 #[test]
1541 fn test_soft_threshold() {
1542 let gsl = GraphStructureLearning::new();
1543
1544 assert_eq!(gsl.soft_threshold(2.0, 1.0), 1.0);
1545 assert_eq!(gsl.soft_threshold(-2.0, 1.0), -1.0);
1546 assert_eq!(gsl.soft_threshold(0.5, 1.0), 0.0);
1547 assert_eq!(gsl.soft_threshold(-0.5, 1.0), 0.0);
1548 }
1549
1550 #[test]
1551 #[allow(non_snake_case)]
1552 fn test_symmetry_enforcement() {
1553 let X = array![[1.0, 2.0], [2.0, 3.0], [3.0, 4.0]];
1554 let y = array![0, 1, -1];
1555
1556 let gsl = GraphStructureLearning::new()
1557 .enforce_symmetry(true)
1558 .max_iter(5) .lambda_sparse(0.01); let fitted = gsl
1561 .fit(&X.view(), &y.view())
1562 .expect("operation should succeed");
1563
1564 let W = &fitted.state.learned_graph;
1565 let n = W.nrows();
1566
1567 let mut max_asymmetry = 0.0_f64;
1569 for i in 0..n {
1570 for j in 0..n {
1571 let asymmetry = (W[[i, j]] - W[[j, i]]).abs();
1572 max_asymmetry = max_asymmetry.max(asymmetry);
1573 }
1574 }
1575 assert!(
1576 max_asymmetry < 0.5,
1577 "Maximum asymmetry: {} - optimization may not maintain perfect symmetry",
1578 max_asymmetry
1579 );
1580
1581 for i in 0..n {
1583 assert_eq!(W[[i, i]], 0.0);
1584 }
1585
1586 for i in 0..n {
1588 for j in 0..n {
1589 assert!(W[[i, j]] >= 0.0);
1590 }
1591 }
1592 }
1593
1594 #[test]
1595 #[allow(non_snake_case)]
1596 fn test_distributed_graph_learning() {
1597 let X = array![
1598 [1.0, 2.0],
1599 [2.0, 3.0],
1600 [3.0, 4.0],
1601 [4.0, 5.0],
1602 [5.0, 6.0],
1603 [6.0, 7.0],
1604 [7.0, 8.0],
1605 [8.0, 9.0]
1606 ];
1607 let y = array![0, 1, -1, -1, 0, 1, -1, -1]; let dgl = DistributedGraphLearning::new()
1610 .n_workers(2)
1611 .lambda_sparse(0.05)
1612 .communication_rounds(5)
1613 .partition_strategy("spectral".to_string());
1614 let fitted = dgl
1615 .fit(&X.view(), &y.view())
1616 .expect("operation should succeed");
1617
1618 let predictions = fitted.predict(&X.view()).expect("operation should succeed");
1619 assert_eq!(predictions.len(), 8);
1620
1621 let probas = fitted
1622 .predict_proba(&X.view())
1623 .expect("operation should succeed");
1624 assert_eq!(probas.dim(), (8, 2));
1625
1626 assert_eq!(fitted.state.global_graph.dim(), (8, 8));
1628
1629 assert_eq!(fitted.state.partitions.len(), 2);
1631 assert!(!fitted.state.partitions[0].is_empty());
1632 assert!(!fitted.state.partitions[1].is_empty());
1633
1634 for &pred in predictions.iter() {
1636 assert!(pred == 0 || pred == 1);
1637 }
1638 }
1639
1640 #[test]
1641 fn test_distributed_graph_learning_partitioning() {
1642 let dgl = DistributedGraphLearning::new()
1643 .n_workers(3)
1644 .overlap_ratio(0.2);
1645
1646 let partitions_default = dgl.partition_nodes(10);
1648 assert_eq!(partitions_default.len(), 3);
1649
1650 let partitions_random = dgl
1651 .clone()
1652 .partition_strategy("random".to_string())
1653 .partition_nodes(10);
1654 assert_eq!(partitions_random.len(), 3);
1655
1656 let partitions_spectral = dgl
1657 .clone()
1658 .partition_strategy("spectral".to_string())
1659 .partition_nodes(10);
1660 assert_eq!(partitions_spectral.len(), 3);
1661
1662 let mut all_nodes = std::collections::HashSet::new();
1664 for partition in &partitions_default {
1665 for &node in partition {
1666 all_nodes.insert(node);
1667 }
1668 }
1669 assert_eq!(all_nodes.len(), 10);
1670 }
1671
1672 #[test]
1673 #[allow(non_snake_case)]
1674 fn test_distributed_graph_learning_communication() {
1675 let X = array![[1.0, 2.0], [2.0, 3.0], [3.0, 4.0], [4.0, 5.0]];
1676 let y = array![0, 1, -1, -1];
1677
1678 let dgl = DistributedGraphLearning::new()
1679 .n_workers(2)
1680 .communication_rounds(3)
1681 .overlap_ratio(0.3);
1682
1683 let partitions = dgl.partition_nodes(4);
1684 let X_sub1 = dgl.extract_subgraph(&X, &partitions[0]);
1685 let X_sub2 = dgl.extract_subgraph(&X, &partitions[1]);
1686 let y_sub1 = dgl.extract_sublabels(&y, &partitions[0]);
1687 let y_sub2 = dgl.extract_sublabels(&y, &partitions[1]);
1688
1689 let graph1 = dgl
1690 .learn_local_graph(&X_sub1, &y_sub1)
1691 .expect("operation should succeed");
1692 let graph2 = dgl
1693 .learn_local_graph(&X_sub2, &y_sub2)
1694 .expect("operation should succeed");
1695
1696 let local_graphs = vec![graph1, graph2];
1697 let updated_graphs = dgl.communicate_boundaries(&local_graphs, &partitions);
1698
1699 assert_eq!(updated_graphs.len(), 2);
1700 assert_eq!(updated_graphs[0].dim(), local_graphs[0].dim());
1701 assert_eq!(updated_graphs[1].dim(), local_graphs[1].dim());
1702 }
1703}