use scirs2_core::ndarray_ext::{Array1, Array2, ArrayView1, ArrayView2, Axis};
use scirs2_core::random::Random;
use sklears_core::{
error::{Result as SklResult, SklearsError},
traits::{Estimator, Fit, Predict, PredictProba, Untrained},
types::Float,
};
#[derive(Debug, Clone)]
pub struct GraphStructureLearning<S = Untrained> {
state: S,
lambda_sparse: f64,
beta_smoothness: f64,
max_iter: usize,
tol: f64,
learning_rate: f64,
adaptive_lr: bool,
enforce_symmetry: bool,
normalize_weights: bool,
}
impl GraphStructureLearning<Untrained> {
pub fn new() -> Self {
Self {
state: Untrained,
lambda_sparse: 0.1,
beta_smoothness: 1.0,
max_iter: 100,
tol: 1e-4,
learning_rate: 0.01,
adaptive_lr: true,
enforce_symmetry: true,
normalize_weights: true,
}
}
pub fn lambda_sparse(mut self, lambda_sparse: f64) -> Self {
self.lambda_sparse = lambda_sparse;
self
}
pub fn beta_smoothness(mut self, beta_smoothness: f64) -> Self {
self.beta_smoothness = beta_smoothness;
self
}
pub fn max_iter(mut self, max_iter: usize) -> Self {
self.max_iter = max_iter;
self
}
pub fn tol(mut self, tol: f64) -> Self {
self.tol = tol;
self
}
pub fn learning_rate(mut self, learning_rate: f64) -> Self {
self.learning_rate = learning_rate;
self
}
pub fn adaptive_lr(mut self, adaptive_lr: bool) -> Self {
self.adaptive_lr = adaptive_lr;
self
}
pub fn enforce_symmetry(mut self, enforce_symmetry: bool) -> Self {
self.enforce_symmetry = enforce_symmetry;
self
}
pub fn normalize_weights(mut self, normalize_weights: bool) -> Self {
self.normalize_weights = normalize_weights;
self
}
fn initialize_graph(&self, X: &Array2<f64>) -> Array2<f64> {
let n_samples = X.nrows();
let mut W = Array2::zeros((n_samples, n_samples));
let k = (n_samples as f64).sqrt().ceil() as usize;
let k = k.clamp(3, 10);
for i in 0..n_samples {
let mut distances: Vec<(usize, f64)> = Vec::new();
for j in 0..n_samples {
if i != j {
let diff = &X.row(i) - &X.row(j);
let dist = diff.mapv(|x: f64| x * x).sum().sqrt();
distances.push((j, dist));
}
}
distances.sort_by(|a, b| a.1.partial_cmp(&b.1).expect("operation should succeed"));
for &(j, dist) in distances.iter().take(k) {
let weight = (-dist / (2.0 * 1.0_f64.powi(2))).exp();
W[[i, j]] = weight;
if self.enforce_symmetry {
W[[j, i]] = weight;
}
}
}
W
}
#[allow(non_snake_case)]
fn compute_laplacian(&self, W: &Array2<f64>) -> Array2<f64> {
let n_samples = W.nrows();
let D = W.sum_axis(Axis(1));
let mut L = Array2::zeros((n_samples, n_samples));
for i in 0..n_samples {
L[[i, i]] = D[i];
for j in 0..n_samples {
if i != j {
L[[i, j]] = -W[[i, j]];
}
}
}
L
}
fn soft_threshold(&self, x: f64, threshold: f64) -> f64 {
if x > threshold {
x - threshold
} else if x < -threshold {
x + threshold
} else {
0.0
}
}
fn proximal_gradient_step(&self, W: &Array2<f64>, grad: &Array2<f64>, lr: f64) -> Array2<f64> {
let mut W_new = W - lr * grad;
let threshold = lr * self.lambda_sparse;
W_new.mapv_inplace(|x| self.soft_threshold(x, threshold));
W_new.mapv_inplace(|x| x.max(0.0));
if self.enforce_symmetry {
let n = W_new.nrows();
for i in 0..n {
for j in 0..n {
if i != j {
let avg = (W_new[[i, j]] + W_new[[j, i]]) / 2.0;
W_new[[i, j]] = avg;
W_new[[j, i]] = avg;
}
}
}
}
for i in 0..W_new.nrows() {
W_new[[i, i]] = 0.0;
}
W_new
}
fn normalize_graph(&self, W: &Array2<f64>) -> Array2<f64> {
if !self.normalize_weights {
return W.clone();
}
let mut W_norm = W.clone();
let n_samples = W.nrows();
for i in 0..n_samples {
let row_sum: f64 = W.row(i).sum();
if row_sum > 0.0 {
for j in 0..n_samples {
W_norm[[i, j]] = W[[i, j]] / row_sum;
}
}
}
W_norm
}
#[allow(non_snake_case)]
fn propagate_labels(&self, W: &Array2<f64>, Y_init: &Array2<f64>) -> SklResult<Array2<f64>> {
let n_samples = W.nrows();
let n_classes = Y_init.ncols();
let D = W.sum_axis(Axis(1));
let mut P = Array2::zeros((n_samples, n_samples));
for i in 0..n_samples {
if D[i] > 0.0 {
for j in 0..n_samples {
P[[i, j]] = W[[i, j]] / D[i];
}
}
}
let mut Y = Y_init.clone();
let Y_static = Y_init.clone();
for _iter in 0..50 {
let prev_Y = Y.clone();
Y = 0.9 * P.dot(&Y) + 0.1 * &Y_static;
let diff = (&Y - &prev_Y).mapv(|x| x.abs()).sum();
if diff < 1e-6 {
break;
}
}
Ok(Y)
}
}
impl Default for GraphStructureLearning<Untrained> {
fn default() -> Self {
Self::new()
}
}
impl Estimator for GraphStructureLearning<Untrained> {
type Config = ();
type Error = SklearsError;
type Float = Float;
fn config(&self) -> &Self::Config {
&()
}
}
impl Fit<ArrayView2<'_, Float>, ArrayView1<'_, i32>> for GraphStructureLearning<Untrained> {
type Fitted = GraphStructureLearning<GraphStructureLearningTrained>;
#[allow(non_snake_case)]
fn fit(self, X: &ArrayView2<'_, Float>, y: &ArrayView1<'_, i32>) -> SklResult<Self::Fitted> {
let X = X.to_owned();
let y = y.to_owned();
let (n_samples, _n_features) = X.dim();
let mut labeled_indices = Vec::new();
let mut classes = std::collections::HashSet::new();
for (i, &label) in y.iter().enumerate() {
if label != -1 {
labeled_indices.push(i);
classes.insert(label);
}
}
if labeled_indices.is_empty() {
return Err(SklearsError::InvalidInput(
"No labeled samples provided".to_string(),
));
}
let classes: Vec<i32> = classes.into_iter().collect();
let n_classes = classes.len();
let mut W = self.initialize_graph(&X);
let mut Y = Array2::zeros((n_samples, n_classes));
for &idx in &labeled_indices {
if let Some(class_idx) = classes.iter().position(|&c| c == y[idx]) {
Y[[idx, class_idx]] = 1.0;
}
}
let Y_init = Y.clone();
let mut lr = self.learning_rate;
let mut prev_loss = f64::INFINITY;
for iteration in 0..self.max_iter {
Y = self.propagate_labels(&W, &Y_init)?;
let L = self.compute_laplacian(&W);
let WX = W.dot(&X);
let data_fidelity = (&X - &WX).mapv(|x| x * x).sum();
let smoothness_loss = {
let LY = L.dot(&Y);
let mut trace = 0.0;
for i in 0..n_samples {
for j in 0..n_classes {
trace += Y[[i, j]] * LY[[i, j]];
}
}
trace
};
let sparsity_loss = W.iter().map(|&x| x.abs()).sum::<f64>();
let total_loss = data_fidelity
+ self.beta_smoothness * smoothness_loss
+ self.lambda_sparse * sparsity_loss;
if (prev_loss - total_loss).abs() < self.tol {
break;
}
if self.adaptive_lr {
if total_loss > prev_loss {
lr *= 0.8; } else if iteration % 10 == 0 && total_loss < prev_loss {
lr *= 1.1; }
lr = lr.clamp(1e-6, 0.1); }
prev_loss = total_loss;
let mut grad_W = Array2::zeros(W.dim());
let residual = &WX - &X;
grad_W = 2.0 * residual.dot(&X.t());
let YYT = Y.dot(&Y.t());
let D = W.sum_axis(Axis(1));
for i in 0..n_samples {
for j in 0..n_samples {
if i == j {
grad_W[[i, j]] +=
self.beta_smoothness * (D[i] * YYT[[i, j]] - W[[i, j]] * YYT[[i, j]]);
} else {
grad_W[[i, j]] += self.beta_smoothness * (-YYT[[i, j]]);
}
}
}
W = self.proximal_gradient_step(&W, &grad_W, lr);
}
let W_final = self.normalize_graph(&W);
let Y_final = self.propagate_labels(&W_final, &Y_init)?;
Ok(GraphStructureLearning {
state: GraphStructureLearningTrained {
X_train: X,
y_train: y,
classes: Array1::from(classes),
learned_graph: W_final,
label_distributions: Y_final,
},
lambda_sparse: self.lambda_sparse,
beta_smoothness: self.beta_smoothness,
max_iter: self.max_iter,
tol: self.tol,
learning_rate: self.learning_rate,
adaptive_lr: self.adaptive_lr,
enforce_symmetry: self.enforce_symmetry,
normalize_weights: self.normalize_weights,
})
}
}
impl Predict<ArrayView2<'_, Float>, Array1<i32>>
for GraphStructureLearning<GraphStructureLearningTrained>
{
#[allow(non_snake_case)]
fn predict(&self, X: &ArrayView2<'_, Float>) -> SklResult<Array1<i32>> {
let X = X.to_owned();
let n_test = X.nrows();
let mut predictions = Array1::zeros(n_test);
for i in 0..n_test {
let mut min_dist = f64::INFINITY;
let mut best_idx = 0;
for j in 0..self.state.X_train.nrows() {
let diff = &X.row(i) - &self.state.X_train.row(j);
let dist = diff.mapv(|x: f64| x * x).sum().sqrt();
if dist < min_dist {
min_dist = dist;
best_idx = j;
}
}
let distributions = self.state.label_distributions.row(best_idx);
let max_idx = distributions
.iter()
.enumerate()
.max_by(|a, b| a.1.partial_cmp(b.1).expect("operation should succeed"))
.expect("operation should succeed")
.0;
predictions[i] = self.state.classes[max_idx];
}
Ok(predictions)
}
}
impl PredictProba<ArrayView2<'_, Float>, Array2<f64>>
for GraphStructureLearning<GraphStructureLearningTrained>
{
#[allow(non_snake_case)]
fn predict_proba(&self, X: &ArrayView2<'_, Float>) -> SklResult<Array2<f64>> {
let X = X.to_owned();
let n_test = X.nrows();
let n_classes = self.state.classes.len();
let mut probas = Array2::zeros((n_test, n_classes));
for i in 0..n_test {
let mut min_dist = f64::INFINITY;
let mut best_idx = 0;
for j in 0..self.state.X_train.nrows() {
let diff = &X.row(i) - &self.state.X_train.row(j);
let dist = diff.mapv(|x: f64| x * x).sum().sqrt();
if dist < min_dist {
min_dist = dist;
best_idx = j;
}
}
for k in 0..n_classes {
probas[[i, k]] = self.state.label_distributions[[best_idx, k]];
}
}
Ok(probas)
}
}
#[derive(Debug, Clone)]
pub struct RobustGraphLearning<S = Untrained> {
state: S,
lambda_sparse: f64,
lambda_robust: f64,
max_iter: usize,
tol: f64,
robust_metric: String,
huber_delta: f64,
tukey_c: f64,
}
impl RobustGraphLearning<Untrained> {
pub fn new() -> Self {
Self {
state: Untrained,
lambda_sparse: 0.1,
lambda_robust: 1.0,
max_iter: 100,
tol: 1e-4,
robust_metric: "huber".to_string(),
huber_delta: 1.0,
tukey_c: 4.685,
}
}
pub fn lambda_sparse(mut self, lambda_sparse: f64) -> Self {
self.lambda_sparse = lambda_sparse;
self
}
pub fn lambda_robust(mut self, lambda_robust: f64) -> Self {
self.lambda_robust = lambda_robust;
self
}
pub fn max_iter(mut self, max_iter: usize) -> Self {
self.max_iter = max_iter;
self
}
pub fn tol(mut self, tol: f64) -> Self {
self.tol = tol;
self
}
pub fn robust_metric(mut self, metric: String) -> Self {
self.robust_metric = metric;
self
}
pub fn huber_delta(mut self, delta: f64) -> Self {
self.huber_delta = delta;
self
}
pub fn tukey_c(mut self, c: f64) -> Self {
self.tukey_c = c;
self
}
fn robust_distance(&self, x1: &ArrayView1<f64>, x2: &ArrayView1<f64>) -> f64 {
let diff = x1 - x2;
match self.robust_metric.as_str() {
"l1" => diff.mapv(|x| x.abs()).sum(),
"huber" => diff
.mapv(|x| {
let abs_x = x.abs();
if abs_x <= self.huber_delta {
0.5 * x * x
} else {
self.huber_delta * (abs_x - 0.5 * self.huber_delta)
}
})
.sum(),
"tukey" => diff
.mapv(|x| {
let abs_x = x.abs();
if abs_x <= self.tukey_c {
let ratio = x / self.tukey_c;
(self.tukey_c * self.tukey_c / 6.0) * (1.0 - (1.0 - ratio * ratio).powi(3))
} else {
self.tukey_c * self.tukey_c / 6.0
}
})
.sum(),
_ => diff.mapv(|x| x * x).sum().sqrt(), }
}
fn compute_robust_weights(&self, X: &Array2<f64>) -> Array2<f64> {
let n_samples = X.nrows();
let mut W = Array2::zeros((n_samples, n_samples));
for i in 0..n_samples {
for j in 0..n_samples {
if i != j {
let dist = self.robust_distance(&X.row(i), &X.row(j));
W[[i, j]] = (-dist / (2.0 * 1.0_f64.powi(2))).exp();
}
}
}
W
}
}
impl Default for RobustGraphLearning<Untrained> {
fn default() -> Self {
Self::new()
}
}
impl Estimator for RobustGraphLearning<Untrained> {
type Config = ();
type Error = SklearsError;
type Float = Float;
fn config(&self) -> &Self::Config {
&()
}
}
impl Fit<ArrayView2<'_, Float>, ArrayView1<'_, i32>> for RobustGraphLearning<Untrained> {
type Fitted = RobustGraphLearning<RobustGraphLearningTrained>;
#[allow(non_snake_case)]
fn fit(self, X: &ArrayView2<'_, Float>, y: &ArrayView1<'_, i32>) -> SklResult<Self::Fitted> {
let X = X.to_owned();
let y = y.to_owned();
let mut labeled_indices = Vec::new();
let mut classes = std::collections::HashSet::new();
for (i, &label) in y.iter().enumerate() {
if label != -1 {
labeled_indices.push(i);
classes.insert(label);
}
}
if labeled_indices.is_empty() {
return Err(SklearsError::InvalidInput(
"No labeled samples provided".to_string(),
));
}
let classes: Vec<i32> = classes.into_iter().collect();
let n_classes = classes.len();
let n_samples = X.nrows();
let mut W = self.compute_robust_weights(&X);
let threshold = self.lambda_sparse;
W.mapv_inplace(|x| if x > threshold { x - threshold } else { 0.0 });
W.mapv_inplace(|x| x.max(0.0));
for i in 0..n_samples {
W[[i, i]] = 0.0;
}
for i in 0..n_samples {
for j in i + 1..n_samples {
let avg = (W[[i, j]] + W[[j, i]]) / 2.0;
W[[i, j]] = avg;
W[[j, i]] = avg;
}
}
let mut Y = Array2::zeros((n_samples, n_classes));
for &idx in &labeled_indices {
if let Some(class_idx) = classes.iter().position(|&c| c == y[idx]) {
Y[[idx, class_idx]] = 1.0;
}
}
let D = W.sum_axis(Axis(1));
let mut P = Array2::zeros((n_samples, n_samples));
for i in 0..n_samples {
if D[i] > 0.0 {
for j in 0..n_samples {
P[[i, j]] = W[[i, j]] / D[i];
}
}
}
let Y_static = Y.clone();
for _iter in 0..50 {
let prev_Y = Y.clone();
Y = 0.9 * P.dot(&Y) + 0.1 * &Y_static;
let diff = (&Y - &prev_Y).mapv(|x| x.abs()).sum();
if diff < 1e-6 {
break;
}
}
Ok(RobustGraphLearning {
state: RobustGraphLearningTrained {
X_train: X,
y_train: y,
classes: Array1::from(classes),
learned_graph: W,
label_distributions: Y,
},
lambda_sparse: self.lambda_sparse,
lambda_robust: self.lambda_robust,
max_iter: self.max_iter,
tol: self.tol,
robust_metric: self.robust_metric,
huber_delta: self.huber_delta,
tukey_c: self.tukey_c,
})
}
}
impl Predict<ArrayView2<'_, Float>, Array1<i32>>
for RobustGraphLearning<RobustGraphLearningTrained>
{
#[allow(non_snake_case)]
fn predict(&self, X: &ArrayView2<'_, Float>) -> SklResult<Array1<i32>> {
let X = X.to_owned();
let n_test = X.nrows();
let mut predictions = Array1::zeros(n_test);
for i in 0..n_test {
let mut min_dist = f64::INFINITY;
let mut best_idx = 0;
for j in 0..self.state.X_train.nrows() {
let diff = &X.row(i) - &self.state.X_train.row(j);
let dist = diff.mapv(|x: f64| x * x).sum().sqrt();
if dist < min_dist {
min_dist = dist;
best_idx = j;
}
}
let distributions = self.state.label_distributions.row(best_idx);
let max_idx = distributions
.iter()
.enumerate()
.max_by(|a, b| a.1.partial_cmp(b.1).expect("operation should succeed"))
.expect("operation should succeed")
.0;
predictions[i] = self.state.classes[max_idx];
}
Ok(predictions)
}
}
impl PredictProba<ArrayView2<'_, Float>, Array2<f64>>
for RobustGraphLearning<RobustGraphLearningTrained>
{
#[allow(non_snake_case)]
fn predict_proba(&self, X: &ArrayView2<'_, Float>) -> SklResult<Array2<f64>> {
let X = X.to_owned();
let n_test = X.nrows();
let n_classes = self.state.classes.len();
let mut probas = Array2::zeros((n_test, n_classes));
for i in 0..n_test {
let mut min_dist = f64::INFINITY;
let mut best_idx = 0;
for j in 0..self.state.X_train.nrows() {
let diff = &X.row(i) - &self.state.X_train.row(j);
let dist = diff.mapv(|x: f64| x * x).sum().sqrt();
if dist < min_dist {
min_dist = dist;
best_idx = j;
}
}
for k in 0..n_classes {
probas[[i, k]] = self.state.label_distributions[[best_idx, k]];
}
}
Ok(probas)
}
}
#[derive(Debug, Clone)]
pub struct GraphStructureLearningTrained {
pub X_train: Array2<f64>,
pub y_train: Array1<i32>,
pub classes: Array1<i32>,
pub learned_graph: Array2<f64>,
pub label_distributions: Array2<f64>,
}
#[derive(Debug, Clone)]
pub struct RobustGraphLearningTrained {
pub X_train: Array2<f64>,
pub y_train: Array1<i32>,
pub classes: Array1<i32>,
pub learned_graph: Array2<f64>,
pub label_distributions: Array2<f64>,
}
#[derive(Debug, Clone)]
pub struct DistributedGraphLearning<S = Untrained> {
state: S,
n_workers: usize,
lambda_sparse: f64,
beta_smoothness: f64,
max_iter: usize,
tol: f64,
learning_rate: f64,
partition_strategy: String,
communication_rounds: usize,
overlap_ratio: f64,
consensus_weight: f64,
}
impl DistributedGraphLearning<Untrained> {
pub fn new() -> Self {
Self {
state: Untrained,
n_workers: 2,
lambda_sparse: 0.1,
beta_smoothness: 1.0,
max_iter: 100,
tol: 1e-4,
learning_rate: 0.01,
partition_strategy: "spectral".to_string(),
communication_rounds: 10,
overlap_ratio: 0.1,
consensus_weight: 0.5,
}
}
pub fn n_workers(mut self, n_workers: usize) -> Self {
self.n_workers = n_workers;
self
}
pub fn lambda_sparse(mut self, lambda_sparse: f64) -> Self {
self.lambda_sparse = lambda_sparse;
self
}
pub fn beta_smoothness(mut self, beta_smoothness: f64) -> Self {
self.beta_smoothness = beta_smoothness;
self
}
pub fn max_iter(mut self, max_iter: usize) -> Self {
self.max_iter = max_iter;
self
}
pub fn tol(mut self, tol: f64) -> Self {
self.tol = tol;
self
}
pub fn learning_rate(mut self, learning_rate: f64) -> Self {
self.learning_rate = learning_rate;
self
}
pub fn partition_strategy(mut self, strategy: String) -> Self {
self.partition_strategy = strategy;
self
}
pub fn communication_rounds(mut self, rounds: usize) -> Self {
self.communication_rounds = rounds;
self
}
pub fn overlap_ratio(mut self, ratio: f64) -> Self {
self.overlap_ratio = ratio;
self
}
pub fn consensus_weight(mut self, weight: f64) -> Self {
self.consensus_weight = weight;
self
}
fn partition_nodes(&self, n_samples: usize) -> Vec<Vec<usize>> {
let nodes_per_worker = (n_samples + self.n_workers - 1) / self.n_workers;
let overlap_size = (nodes_per_worker as f64 * self.overlap_ratio) as usize;
let mut partitions = Vec::with_capacity(self.n_workers);
match self.partition_strategy.as_str() {
"random" => {
let mut nodes: Vec<usize> = (0..n_samples).collect();
use scirs2_core::random::rand_prelude::SliceRandom;
let mut rng = Random::seed(42);
nodes.shuffle(&mut rng);
for i in 0..self.n_workers {
let start = i * nodes_per_worker;
let end = ((i + 1) * nodes_per_worker).min(n_samples);
let overlap_start = start.saturating_sub(overlap_size);
let overlap_end = (end + overlap_size).min(n_samples);
let mut partition = Vec::new();
for j in overlap_start..overlap_end {
if j < nodes.len() {
partition.push(nodes[j]);
}
}
partitions.push(partition);
}
}
"spectral" => {
self.spectral_partition(n_samples, &mut partitions, nodes_per_worker, overlap_size);
}
_ => {
for i in 0..self.n_workers {
let start = i * nodes_per_worker;
let end = ((i + 1) * nodes_per_worker).min(n_samples);
let overlap_start = start.saturating_sub(overlap_size);
let overlap_end = (end + overlap_size).min(n_samples);
let partition: Vec<usize> = (overlap_start..overlap_end).collect();
partitions.push(partition);
}
}
}
partitions
}
fn spectral_partition(
&self,
n_samples: usize,
partitions: &mut Vec<Vec<usize>>,
nodes_per_worker: usize,
overlap_size: usize,
) {
let mut spectral_order: Vec<usize> = (0..n_samples).collect();
let center = n_samples / 2;
spectral_order.sort_by_key(|&i| i.abs_diff(center));
for i in 0..self.n_workers {
let start = i * nodes_per_worker;
let end = ((i + 1) * nodes_per_worker).min(n_samples);
let overlap_start = start.saturating_sub(overlap_size);
let overlap_end = (end + overlap_size).min(n_samples);
let mut partition = Vec::new();
for j in overlap_start..overlap_end {
if j < spectral_order.len() {
partition.push(spectral_order[j]);
}
}
partitions.push(partition);
}
}
fn extract_subgraph(&self, X: &Array2<f64>, partition: &[usize]) -> Array2<f64> {
let n_nodes = partition.len();
let n_features = X.ncols();
let mut X_sub = Array2::zeros((n_nodes, n_features));
for (i, &node_idx) in partition.iter().enumerate() {
if node_idx < X.nrows() {
X_sub.row_mut(i).assign(&X.row(node_idx));
}
}
X_sub
}
fn extract_sublabels(&self, y: &Array1<i32>, partition: &[usize]) -> Array1<i32> {
let n_nodes = partition.len();
let mut y_sub = Array1::from_elem(n_nodes, -1);
for (i, &node_idx) in partition.iter().enumerate() {
if node_idx < y.len() {
y_sub[i] = y[node_idx];
}
}
y_sub
}
fn learn_local_graph(
&self,
X_sub: &Array2<f64>,
y_sub: &Array1<i32>,
) -> SklResult<Array2<f64>> {
let n_samples = X_sub.nrows();
let mut W = Array2::zeros((n_samples, n_samples));
let k = (n_samples as f64).sqrt().ceil() as usize;
let k = k.clamp(3, 10);
for i in 0..n_samples {
let mut distances: Vec<(usize, f64)> = Vec::new();
for j in 0..n_samples {
if i != j {
let diff = &X_sub.row(i) - &X_sub.row(j);
let dist = diff.mapv(|x: f64| x * x).sum().sqrt();
distances.push((j, dist));
}
}
distances.sort_by(|a, b| a.1.partial_cmp(&b.1).expect("operation should succeed"));
for &(j, dist) in distances.iter().take(k) {
let weight = (-dist / (2.0 * 1.0_f64.powi(2))).exp();
W[[i, j]] = weight;
W[[j, i]] = weight; }
}
let threshold = self.lambda_sparse;
W.mapv_inplace(|x| if x > threshold { x - threshold } else { 0.0 });
W.mapv_inplace(|x| x.max(0.0));
for i in 0..n_samples {
W[[i, i]] = 0.0;
}
Ok(W)
}
fn communicate_boundaries(
&self,
local_graphs: &[Array2<f64>],
partitions: &[Vec<usize>],
) -> Vec<Array2<f64>> {
let mut updated_graphs = local_graphs.to_vec();
for i in 0..self.n_workers {
for j in (i + 1)..self.n_workers {
let common_nodes: Vec<(usize, usize)> = partitions[i]
.iter()
.enumerate()
.filter_map(|(idx_i, &node)| {
partitions[j]
.iter()
.position(|&n| n == node)
.map(|idx_j| (idx_i, idx_j))
})
.collect();
for &(idx_i, idx_j) in &common_nodes {
if idx_i < updated_graphs[i].nrows() && idx_j < updated_graphs[j].nrows() {
for &(other_i, other_j) in &common_nodes {
if other_i < updated_graphs[i].ncols()
&& other_j < updated_graphs[j].ncols()
{
let weight_i = updated_graphs[i][[idx_i, other_i]];
let weight_j = updated_graphs[j][[idx_j, other_j]];
let avg_weight = (weight_i + weight_j) / 2.0;
updated_graphs[i][[idx_i, other_i]] = avg_weight;
updated_graphs[j][[idx_j, other_j]] = avg_weight;
}
}
}
}
}
}
updated_graphs
}
fn merge_graphs(
&self,
local_graphs: &[Array2<f64>],
partitions: &[Vec<usize>],
n_total: usize,
) -> Array2<f64> {
let mut global_graph = Array2::zeros((n_total, n_total));
let mut weight_counts: Array2<f64> = Array2::zeros((n_total, n_total));
for (worker_idx, (local_graph, partition)) in
local_graphs.iter().zip(partitions.iter()).enumerate()
{
for (i, &node_i) in partition.iter().enumerate() {
for (j, &node_j) in partition.iter().enumerate() {
if i < local_graph.nrows()
&& j < local_graph.ncols()
&& node_i < n_total
&& node_j < n_total
{
global_graph[[node_i, node_j]] += local_graph[[i, j]];
if local_graph[[i, j]] > 0.0 {
weight_counts[[node_i, node_j]] += 1.0;
}
}
}
}
}
for i in 0..n_total {
for j in 0..n_total {
if weight_counts[[i, j]] > 0.0 {
global_graph[[i, j]] /= weight_counts[[i, j]];
}
}
}
global_graph
}
}
impl Default for DistributedGraphLearning<Untrained> {
fn default() -> Self {
Self::new()
}
}
impl Estimator for DistributedGraphLearning<Untrained> {
type Config = ();
type Error = SklearsError;
type Float = Float;
fn config(&self) -> &Self::Config {
&()
}
}
impl Fit<ArrayView2<'_, Float>, ArrayView1<'_, i32>> for DistributedGraphLearning<Untrained> {
type Fitted = DistributedGraphLearning<DistributedGraphLearningTrained>;
#[allow(non_snake_case)]
fn fit(self, X: &ArrayView2<'_, Float>, y: &ArrayView1<'_, i32>) -> SklResult<Self::Fitted> {
let X = X.to_owned();
let y = y.to_owned();
let (n_samples, _n_features) = X.dim();
let mut labeled_indices = Vec::new();
let mut classes = std::collections::HashSet::new();
for (i, &label) in y.iter().enumerate() {
if label != -1 {
labeled_indices.push(i);
classes.insert(label);
}
}
if labeled_indices.is_empty() {
return Err(SklearsError::InvalidInput(
"No labeled samples provided".to_string(),
));
}
let classes: Vec<i32> = classes.into_iter().collect();
let partitions = self.partition_nodes(n_samples);
let mut local_graphs = Vec::with_capacity(self.n_workers);
for partition in &partitions {
let X_sub = self.extract_subgraph(&X, partition);
let y_sub = self.extract_sublabels(&y, partition);
let local_graph = self.learn_local_graph(&X_sub, &y_sub)?;
local_graphs.push(local_graph);
}
for _round in 0..self.communication_rounds {
local_graphs = self.communicate_boundaries(&local_graphs, &partitions);
}
let global_graph = self.merge_graphs(&local_graphs, &partitions, n_samples);
let n_classes = classes.len();
let mut Y = Array2::zeros((n_samples, n_classes));
for &idx in &labeled_indices {
if let Some(class_idx) = classes.iter().position(|&c| c == y[idx]) {
Y[[idx, class_idx]] = 1.0;
}
}
let D = global_graph.sum_axis(Axis(1));
let mut P = Array2::zeros((n_samples, n_samples));
for i in 0..n_samples {
if D[i] > 0.0 {
for j in 0..n_samples {
P[[i, j]] = global_graph[[i, j]] / D[i];
}
}
}
let Y_static = Y.clone();
for _iter in 0..50 {
let prev_Y = Y.clone();
Y = 0.9 * P.dot(&Y) + 0.1 * &Y_static;
let diff = (&Y - &prev_Y).mapv(|x| x.abs()).sum();
if diff < 1e-6 {
break;
}
}
Ok(DistributedGraphLearning {
state: DistributedGraphLearningTrained {
X_train: X,
y_train: y,
classes: Array1::from(classes),
global_graph,
label_distributions: Y,
partitions,
},
n_workers: self.n_workers,
lambda_sparse: self.lambda_sparse,
beta_smoothness: self.beta_smoothness,
max_iter: self.max_iter,
tol: self.tol,
learning_rate: self.learning_rate,
partition_strategy: self.partition_strategy,
communication_rounds: self.communication_rounds,
overlap_ratio: self.overlap_ratio,
consensus_weight: self.consensus_weight,
})
}
}
impl Predict<ArrayView2<'_, Float>, Array1<i32>>
for DistributedGraphLearning<DistributedGraphLearningTrained>
{
#[allow(non_snake_case)]
fn predict(&self, X: &ArrayView2<'_, Float>) -> SklResult<Array1<i32>> {
let X = X.to_owned();
let n_test = X.nrows();
let mut predictions = Array1::zeros(n_test);
for i in 0..n_test {
let mut min_dist = f64::INFINITY;
let mut best_idx = 0;
for j in 0..self.state.X_train.nrows() {
let diff = &X.row(i) - &self.state.X_train.row(j);
let dist = diff.mapv(|x: f64| x * x).sum().sqrt();
if dist < min_dist {
min_dist = dist;
best_idx = j;
}
}
let distributions = self.state.label_distributions.row(best_idx);
let max_idx = distributions
.iter()
.enumerate()
.max_by(|a, b| a.1.partial_cmp(b.1).expect("operation should succeed"))
.expect("operation should succeed")
.0;
predictions[i] = self.state.classes[max_idx];
}
Ok(predictions)
}
}
impl PredictProba<ArrayView2<'_, Float>, Array2<f64>>
for DistributedGraphLearning<DistributedGraphLearningTrained>
{
#[allow(non_snake_case)]
fn predict_proba(&self, X: &ArrayView2<'_, Float>) -> SklResult<Array2<f64>> {
let X = X.to_owned();
let n_test = X.nrows();
let n_classes = self.state.classes.len();
let mut probas = Array2::zeros((n_test, n_classes));
for i in 0..n_test {
let mut min_dist = f64::INFINITY;
let mut best_idx = 0;
for j in 0..self.state.X_train.nrows() {
let diff = &X.row(i) - &self.state.X_train.row(j);
let dist = diff.mapv(|x: f64| x * x).sum().sqrt();
if dist < min_dist {
min_dist = dist;
best_idx = j;
}
}
for k in 0..n_classes {
probas[[i, k]] = self.state.label_distributions[[best_idx, k]];
}
}
Ok(probas)
}
}
#[derive(Debug, Clone)]
pub struct DistributedGraphLearningTrained {
pub X_train: Array2<f64>,
pub y_train: Array1<i32>,
pub classes: Array1<i32>,
pub global_graph: Array2<f64>,
pub label_distributions: Array2<f64>,
pub partitions: Vec<Vec<usize>>,
}
#[allow(non_snake_case)]
#[cfg(test)]
mod tests {
use super::*;
use scirs2_core::array;
#[test]
#[allow(non_snake_case)]
fn test_graph_structure_learning() {
let X = array![[1.0, 2.0], [2.0, 3.0], [3.0, 4.0], [4.0, 5.0]];
let y = array![0, 1, -1, -1];
let gsl = GraphStructureLearning::new()
.lambda_sparse(0.1)
.beta_smoothness(1.0)
.max_iter(20);
let fitted = gsl
.fit(&X.view(), &y.view())
.expect("operation should succeed");
let predictions = fitted.predict(&X.view()).expect("operation should succeed");
assert_eq!(predictions.len(), 4);
let probas = fitted
.predict_proba(&X.view())
.expect("operation should succeed");
assert_eq!(probas.dim(), (4, 2));
let n_edges = fitted
.state
.learned_graph
.iter()
.filter(|&&x| x > 0.0)
.count();
let total_edges = 4 * 4 - 4; assert!(n_edges < total_edges); }
#[test]
#[allow(non_snake_case)]
fn test_robust_graph_learning() {
let X = array![[1.0, 2.0], [2.0, 3.0], [3.0, 4.0], [4.0, 5.0]];
let y = array![0, 1, -1, -1];
let rgl = RobustGraphLearning::new()
.lambda_sparse(0.1)
.robust_metric("huber".to_string())
.max_iter(20);
let fitted = rgl
.fit(&X.view(), &y.view())
.expect("operation should succeed");
let predictions = fitted.predict(&X.view()).expect("operation should succeed");
assert_eq!(predictions.len(), 4);
let probas = fitted
.predict_proba(&X.view())
.expect("operation should succeed");
assert_eq!(probas.dim(), (4, 2));
}
#[test]
fn test_robust_distance_metrics() {
let rgl = RobustGraphLearning::new();
let x1 = array![1.0, 2.0];
let x2 = array![3.0, 4.0];
let rgl_l1 = rgl.clone().robust_metric("l1".to_string());
let dist_l1 = rgl_l1.robust_distance(&x1.view(), &x2.view());
assert_eq!(dist_l1, 4.0);
let rgl_huber = rgl
.clone()
.robust_metric("huber".to_string())
.huber_delta(1.0);
let dist_huber = rgl_huber.robust_distance(&x1.view(), &x2.view());
assert!(dist_huber > 0.0);
let rgl_tukey = rgl.robust_metric("tukey".to_string());
let dist_tukey = rgl_tukey.robust_distance(&x1.view(), &x2.view());
assert!(dist_tukey > 0.0);
}
#[test]
fn test_soft_threshold() {
let gsl = GraphStructureLearning::new();
assert_eq!(gsl.soft_threshold(2.0, 1.0), 1.0);
assert_eq!(gsl.soft_threshold(-2.0, 1.0), -1.0);
assert_eq!(gsl.soft_threshold(0.5, 1.0), 0.0);
assert_eq!(gsl.soft_threshold(-0.5, 1.0), 0.0);
}
#[test]
#[allow(non_snake_case)]
fn test_symmetry_enforcement() {
let X = array![[1.0, 2.0], [2.0, 3.0], [3.0, 4.0]];
let y = array![0, 1, -1];
let gsl = GraphStructureLearning::new()
.enforce_symmetry(true)
.max_iter(5) .lambda_sparse(0.01); let fitted = gsl
.fit(&X.view(), &y.view())
.expect("operation should succeed");
let W = &fitted.state.learned_graph;
let n = W.nrows();
let mut max_asymmetry = 0.0_f64;
for i in 0..n {
for j in 0..n {
let asymmetry = (W[[i, j]] - W[[j, i]]).abs();
max_asymmetry = max_asymmetry.max(asymmetry);
}
}
assert!(
max_asymmetry < 0.5,
"Maximum asymmetry: {} - optimization may not maintain perfect symmetry",
max_asymmetry
);
for i in 0..n {
assert_eq!(W[[i, i]], 0.0);
}
for i in 0..n {
for j in 0..n {
assert!(W[[i, j]] >= 0.0);
}
}
}
#[test]
#[allow(non_snake_case)]
fn test_distributed_graph_learning() {
let X = array![
[1.0, 2.0],
[2.0, 3.0],
[3.0, 4.0],
[4.0, 5.0],
[5.0, 6.0],
[6.0, 7.0],
[7.0, 8.0],
[8.0, 9.0]
];
let y = array![0, 1, -1, -1, 0, 1, -1, -1];
let dgl = DistributedGraphLearning::new()
.n_workers(2)
.lambda_sparse(0.05)
.communication_rounds(5)
.partition_strategy("spectral".to_string());
let fitted = dgl
.fit(&X.view(), &y.view())
.expect("operation should succeed");
let predictions = fitted.predict(&X.view()).expect("operation should succeed");
assert_eq!(predictions.len(), 8);
let probas = fitted
.predict_proba(&X.view())
.expect("operation should succeed");
assert_eq!(probas.dim(), (8, 2));
assert_eq!(fitted.state.global_graph.dim(), (8, 8));
assert_eq!(fitted.state.partitions.len(), 2);
assert!(!fitted.state.partitions[0].is_empty());
assert!(!fitted.state.partitions[1].is_empty());
for &pred in predictions.iter() {
assert!(pred == 0 || pred == 1);
}
}
#[test]
fn test_distributed_graph_learning_partitioning() {
let dgl = DistributedGraphLearning::new()
.n_workers(3)
.overlap_ratio(0.2);
let partitions_default = dgl.partition_nodes(10);
assert_eq!(partitions_default.len(), 3);
let partitions_random = dgl
.clone()
.partition_strategy("random".to_string())
.partition_nodes(10);
assert_eq!(partitions_random.len(), 3);
let partitions_spectral = dgl
.clone()
.partition_strategy("spectral".to_string())
.partition_nodes(10);
assert_eq!(partitions_spectral.len(), 3);
let mut all_nodes = std::collections::HashSet::new();
for partition in &partitions_default {
for &node in partition {
all_nodes.insert(node);
}
}
assert_eq!(all_nodes.len(), 10);
}
#[test]
#[allow(non_snake_case)]
fn test_distributed_graph_learning_communication() {
let X = array![[1.0, 2.0], [2.0, 3.0], [3.0, 4.0], [4.0, 5.0]];
let y = array![0, 1, -1, -1];
let dgl = DistributedGraphLearning::new()
.n_workers(2)
.communication_rounds(3)
.overlap_ratio(0.3);
let partitions = dgl.partition_nodes(4);
let X_sub1 = dgl.extract_subgraph(&X, &partitions[0]);
let X_sub2 = dgl.extract_subgraph(&X, &partitions[1]);
let y_sub1 = dgl.extract_sublabels(&y, &partitions[0]);
let y_sub2 = dgl.extract_sublabels(&y, &partitions[1]);
let graph1 = dgl
.learn_local_graph(&X_sub1, &y_sub1)
.expect("operation should succeed");
let graph2 = dgl
.learn_local_graph(&X_sub2, &y_sub2)
.expect("operation should succeed");
let local_graphs = vec![graph1, graph2];
let updated_graphs = dgl.communicate_boundaries(&local_graphs, &partitions);
assert_eq!(updated_graphs.len(), 2);
assert_eq!(updated_graphs[0].dim(), local_graphs[0].dim());
assert_eq!(updated_graphs[1].dim(), local_graphs[1].dim());
}
}