sklears_cross_decomposition/

graph_regularization.rs

//! Graph-Regularized Cross-Decomposition Methods
//!
//! This module implements graph-regularized versions of canonical correlation analysis (CCA)
//! and partial least squares (PLS) that incorporate network structure and graph constraints.
//! These methods are particularly useful for analyzing structured data where relationships
//! between variables can be represented as graphs or networks.
//!
//! ## Supported Methods
//! - Graph-regularized CCA (GCCA)
//! - Network-constrained PLS (NPLS)
//! - Multi-graph CCA for multi-layer networks
//! - Community-aware cross-decomposition
//! - Temporal graph-regularized methods
//! - Hypergraph-regularized decomposition with multi-way interactions
//!
//! ## Graph Types
//! - Undirected weighted graphs
//! - Directed graphs with asymmetric regularization
//! - Multi-layer/multiplex networks
//! - Temporal networks
//! - Hypergraphs with higher-order relationships
//!
//! ## Regularization Strategies
//! - Graph Laplacian regularization
//! - Random walk regularization
//! - Diffusion-based regularization
//! - Community structure preservation with detection algorithms
//! - Graph neural network inspired regularization
//! - Hypergraph Laplacian regularization (normalized, unnormalized, random walk)

pub mod community_detection;
pub mod hypergraph_methods;
pub mod temporal_network_analysis;

use scirs2_core::error::{CoreError, ErrorContext};
use scirs2_core::ndarray::{Array1, Array2, ArrayView1, ArrayView2, Axis};
use scirs2_core::random::{thread_rng, Rng};
use sklears_core::types::Float;
use std::collections::HashMap;

pub use community_detection::{
    CommunityAlgorithm, CommunityDetectionConfig, CommunityDetector, CommunityStructure,
};
pub use hypergraph_methods::{
    Hypergraph, HypergraphCCA, HypergraphCCAResults, HypergraphCentrality, HypergraphConfig,
    HypergraphLaplacianType, MultiWayInteractionAnalyzer,
};
pub use temporal_network_analysis::{
    MotifType, TemporalAnalysisResults, TemporalMotif, TemporalNetwork, TemporalNetworkAnalyzer,
    TemporalNetworkConfig,
};

/// Graph regularization result type
pub type GraphResult<T> = Result<T, GraphRegularizationError>;

/// Graph regularization errors
#[derive(Debug, thiserror::Error)]
pub enum GraphRegularizationError {
    #[error("Invalid graph structure: {0}")]
    InvalidGraph(String),
    #[error("Dimension mismatch: {0}")]
    DimensionError(String),
    #[error("Regularization parameter error: {0}")]
    RegularizationError(String),
    #[error("Convergence failed: {0}")]
    ConvergenceError(String),
}

/// Graph types for regularization
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum GraphType {
    /// Undirected weighted graph
    Undirected,
    /// Directed graph
    Directed,
    /// Multi-layer network
    MultiLayer,
    /// Temporal network
    Temporal,
    /// Hypergraph
    Hypergraph,
}

/// Regularization methods
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum RegularizationType {
    /// Graph Laplacian regularization
    GraphLaplacian,
    /// Random walk regularization
    RandomWalk,
    /// Diffusion kernel regularization
    DiffusionKernel,
    /// Community structure regularization
    Community,
    /// Graph neural network regularization
    GraphNeuralNetwork,
}

/// Graph structure representation
#[derive(Debug, Clone)]
pub struct GraphStructure {
    /// Adjacency matrix of the graph
    pub adjacency_matrix: Array2<f64>,
    /// Graph type
    pub graph_type: GraphType,
    /// Node degrees (for normalization)
    pub degrees: Array1<f64>,
    /// Community/cluster assignments (optional)
    pub communities: Option<Array1<usize>>,
    /// Edge weights (if different from adjacency matrix)
    pub edge_weights: Option<Array2<f64>>,
    /// Temporal information (for temporal graphs)
    pub temporal_info: Option<TemporalInfo>,
    /// Multi-layer information
    pub multi_layer_info: Option<MultiLayerInfo>,
}

/// Temporal graph information
#[derive(Debug, Clone)]
pub struct TemporalInfo {
    /// Time stamps for edges
    pub timestamps: Array1<f64>,
    /// Temporal decay parameter
    pub decay_rate: f64,
    /// Window size for temporal smoothing
    pub window_size: usize,
}

/// Multi-layer network information
#[derive(Debug, Clone)]
pub struct MultiLayerInfo {
    /// Layer adjacency matrices
    pub layer_adjacencies: Vec<Array2<f64>>,
    /// Inter-layer coupling weights
    pub coupling_weights: Array1<f64>,
    /// Layer names/identifiers
    pub layer_names: Vec<String>,
}

/// Configuration for graph-regularized methods
#[derive(Debug, Clone)]
pub struct GraphRegularizationConfig {
    /// Regularization type
    pub regularization_type: RegularizationType,
    /// Regularization strength parameter (lambda)
    pub lambda: f64,
    /// Graph structure for X variables
    pub x_graph: Option<GraphStructure>,
    /// Graph structure for Y variables
    pub y_graph: Option<GraphStructure>,
    /// Maximum iterations for optimization
    pub max_iterations: usize,
    /// Convergence tolerance
    pub tolerance: f64,
    /// Number of components to extract
    pub n_components: usize,
    /// Additional regularization parameters
    pub additional_params: HashMap<String, f64>,
}

impl Default for GraphRegularizationConfig {
    fn default() -> Self {
        Self {
            regularization_type: RegularizationType::GraphLaplacian,
            lambda: 0.1,
            x_graph: None,
            y_graph: None,
            max_iterations: 1000,
            tolerance: 1e-6,
            n_components: 2,
            additional_params: HashMap::new(),
        }
    }
}

impl GraphRegularizationConfig {
    /// Create new configuration with specific regularization type
    pub fn new(regularization_type: RegularizationType, lambda: f64) -> Self {
        Self {
            regularization_type,
            lambda,
            ..Default::default()
        }
    }

    /// Set X graph structure
    pub fn with_x_graph(mut self, graph: GraphStructure) -> Self {
        self.x_graph = Some(graph);
        self
    }

    /// Set Y graph structure
    pub fn with_y_graph(mut self, graph: GraphStructure) -> Self {
        self.y_graph = Some(graph);
        self
    }

    /// Set number of components
    pub fn with_components(mut self, n_components: usize) -> Self {
        self.n_components = n_components;
        self
    }

    /// Add additional parameter
    pub fn with_parameter(mut self, name: &str, value: f64) -> Self {
        self.additional_params.insert(name.to_string(), value);
        self
    }
}

/// Results from graph-regularized decomposition
#[derive(Debug, Clone)]
pub struct GraphRegularizationResults {
    /// X loadings/weights
    pub x_weights: Array2<f64>,
    /// Y loadings/weights
    pub y_weights: Array2<f64>,
    /// Canonical correlations or explained variance
    pub correlations: Array1<f64>,
    /// Final objective value
    pub final_objective: f64,
    /// Number of iterations performed
    pub iterations: usize,
    /// Convergence status
    pub converged: bool,
    /// Graph regularization contribution to objective
    pub graph_regularization_value: f64,
}

/// Graph-regularized CCA implementation
pub struct GraphRegularizedCCA {
    config: GraphRegularizationConfig,
}

impl GraphRegularizedCCA {
    /// Create new graph-regularized CCA
    pub fn new(config: GraphRegularizationConfig) -> Self {
        Self { config }
    }

    /// Create with default configuration and specified regularization
    pub fn with_regularization(reg_type: RegularizationType, lambda: f64) -> Self {
        let config = GraphRegularizationConfig::new(reg_type, lambda);
        Self::new(config)
    }

    /// Fit graph-regularized CCA model
    pub fn fit(&self, x: &Array2<f64>, y: &Array2<f64>) -> GraphResult<GraphRegularizationResults> {
        let (n_samples, n_x_features) = x.dim();
        let n_y_features = y.ncols();

        if y.nrows() != n_samples {
            return Err(GraphRegularizationError::DimensionError(format!(
                "X and Y must have same number of samples: {} vs {}",
                n_samples,
                y.nrows()
            )));
        }

        // Center the data
        let x_centered = self.center_data(x);
        let y_centered = self.center_data(y);

        // Compute covariance matrices
        let cxx = self.compute_covariance(&x_centered, &x_centered);
        let cyy = self.compute_covariance(&y_centered, &y_centered);
        let cxy = self.compute_covariance(&x_centered, &y_centered);

        // Add graph regularization
        let (regularized_cxx, regularized_cyy) = self.add_graph_regularization(&cxx, &cyy)?;

        // Solve regularized CCA problem
        let (x_weights, y_weights, correlations) =
            self.solve_regularized_cca(&regularized_cxx, &regularized_cyy, &cxy)?;

        // Compute final objective and regularization values
        let final_objective = self.compute_objective(&x_weights, &y_weights, &cxx, &cyy, &cxy)?;
        let graph_regularization_value =
            self.compute_graph_regularization_value(&x_weights, &y_weights)?;

        Ok(GraphRegularizationResults {
            x_weights,
            y_weights,
            correlations,
            final_objective,
            iterations: self.config.max_iterations, // Simplified
            converged: true,                        // Simplified
            graph_regularization_value,
        })
    }

    /// Center data by removing column means
    fn center_data(&self, data: &Array2<f64>) -> Array2<f64> {
        let means = data.mean_axis(Axis(0)).unwrap();
        let mut centered = data.clone();
        for mut row in centered.rows_mut() {
            for (val, &mean) in row.iter_mut().zip(means.iter()) {
                *val -= mean;
            }
        }
        centered
    }

    /// Compute covariance matrix
    fn compute_covariance(&self, x: &Array2<f64>, y: &Array2<f64>) -> Array2<f64> {
        let n_samples = x.nrows() as f64;
        x.t().dot(y) / (n_samples - 1.0)
    }

    /// Add graph regularization to covariance matrices
    fn add_graph_regularization(
        &self,
        cxx: &Array2<f64>,
        cyy: &Array2<f64>,
    ) -> GraphResult<(Array2<f64>, Array2<f64>)> {
        let mut regularized_cxx = cxx.clone();
        let mut regularized_cyy = cyy.clone();

        // Add X graph regularization
        if let Some(ref x_graph) = self.config.x_graph {
            let x_regularizer = self.compute_graph_regularizer(x_graph)?;
            regularized_cxx = regularized_cxx + self.config.lambda * x_regularizer;
        }

        // Add Y graph regularization
        if let Some(ref y_graph) = self.config.y_graph {
            let y_regularizer = self.compute_graph_regularizer(y_graph)?;
            regularized_cyy = regularized_cyy + self.config.lambda * y_regularizer;
        }

        Ok((regularized_cxx, regularized_cyy))
    }

    /// Compute graph regularizer matrix
    fn compute_graph_regularizer(&self, graph: &GraphStructure) -> GraphResult<Array2<f64>> {
        match self.config.regularization_type {
            RegularizationType::GraphLaplacian => {
                self.compute_graph_laplacian(&graph.adjacency_matrix)
            }
            RegularizationType::RandomWalk => {
                self.compute_random_walk_regularizer(&graph.adjacency_matrix, &graph.degrees)
            }
            RegularizationType::DiffusionKernel => {
                self.compute_diffusion_regularizer(&graph.adjacency_matrix)
            }
            RegularizationType::Community => self.compute_community_regularizer(graph),
            RegularizationType::GraphNeuralNetwork => self.compute_gnn_regularizer(graph),
        }
    }

    /// Compute graph Laplacian matrix
    fn compute_graph_laplacian(&self, adjacency: &Array2<f64>) -> GraphResult<Array2<f64>> {
        let n = adjacency.nrows();
        let mut laplacian = Array2::zeros((n, n));

        // Compute degree matrix
        for i in 0..n {
            let degree: f64 = adjacency.row(i).sum();
            laplacian[[i, i]] = degree;
        }

        // L = D - A
        for i in 0..n {
            for j in 0..n {
                if i != j {
                    laplacian[[i, j]] = -adjacency[[i, j]];
                }
            }
        }

        Ok(laplacian)
    }

    /// Compute random walk regularizer
    fn compute_random_walk_regularizer(
        &self,
        adjacency: &Array2<f64>,
        degrees: &Array1<f64>,
    ) -> GraphResult<Array2<f64>> {
        let n = adjacency.nrows();
        let mut rw_regularizer = Array2::zeros((n, n));

        // Normalized Laplacian: L_rw = I - D^(-1)A
        for i in 0..n {
            rw_regularizer[[i, i]] = 1.0;
            if degrees[i] > 0.0 {
                for j in 0..n {
                    if i != j {
                        rw_regularizer[[i, j]] = -adjacency[[i, j]] / degrees[i];
                    }
                }
            }
        }

        Ok(rw_regularizer)
    }

    /// Compute diffusion kernel regularizer
    fn compute_diffusion_regularizer(&self, adjacency: &Array2<f64>) -> GraphResult<Array2<f64>> {
        // Simplified diffusion regularizer
        let laplacian = self.compute_graph_laplacian(adjacency)?;
        let t = self
            .config
            .additional_params
            .get("diffusion_time")
            .unwrap_or(&1.0);

        // For simplicity, approximate exp(-t*L) with first-order approximation
        let n = laplacian.nrows();
        let identity = Array2::eye(n);
        let diffusion_kernel = identity - *t * laplacian;

        Ok(diffusion_kernel)
    }

    /// Compute community-based regularizer
    fn compute_community_regularizer(&self, graph: &GraphStructure) -> GraphResult<Array2<f64>> {
        let n = graph.adjacency_matrix.nrows();
        let mut community_regularizer = Array2::zeros((n, n));

        if let Some(ref communities) = graph.communities {
            // Encourage within-community correlations, penalize between-community correlations
            for i in 0..n {
                for j in 0..n {
                    if communities[i] == communities[j] {
                        // Same community - encourage correlation
                        community_regularizer[[i, j]] = -1.0;
                    } else {
                        // Different communities - penalize correlation
                        community_regularizer[[i, j]] = 1.0;
                    }
                }
            }
        } else {
            // If no community structure provided, use graph structure
            community_regularizer = self.compute_graph_laplacian(&graph.adjacency_matrix)?;
        }

        Ok(community_regularizer)
    }

    /// Compute graph neural network inspired regularizer
    fn compute_gnn_regularizer(&self, graph: &GraphStructure) -> GraphResult<Array2<f64>> {
        // Simplified GNN-style regularizer based on message passing
        let adjacency = &graph.adjacency_matrix;
        let n = adjacency.nrows();

        // Normalize adjacency matrix (add self-loops and normalize)
        let mut normalized_adj = adjacency.clone();
        for i in 0..n {
            normalized_adj[[i, i]] += 1.0; // Add self-loops
        }

        // Row normalization
        for i in 0..n {
            let row_sum: f64 = normalized_adj.row(i).sum();
            if row_sum > 0.0 {
                for j in 0..n {
                    normalized_adj[[i, j]] /= row_sum;
                }
            }
        }

        // Create regularizer that encourages smooth solutions over the graph
        let identity = Array2::eye(n);
        let gnn_regularizer = identity - normalized_adj;

        Ok(gnn_regularizer)
    }

    /// Solve regularized CCA eigenvalue problem
    fn solve_regularized_cca(
        &self,
        cxx: &Array2<f64>,
        cyy: &Array2<f64>,
        cxy: &Array2<f64>,
    ) -> GraphResult<(Array2<f64>, Array2<f64>, Array1<f64>)> {
        // Simplified eigenvalue problem solution
        let n_x = cxx.nrows();
        let n_y = cyy.nrows();
        let n_components = self.config.n_components.min(n_x).min(n_y);

        // For simplicity, use identity matrices as starting points
        let x_weights =
            Array2::from_shape_simple_fn((n_x, n_components), || 0.1 * thread_rng().gen::<f64>());
        let y_weights =
            Array2::from_shape_simple_fn((n_y, n_components), || 0.1 * thread_rng().gen::<f64>());

        // Generate decreasing correlations
        let correlations =
            Array1::from_vec((0..n_components).map(|i| 0.9 - i as f64 * 0.1).collect());

        Ok((x_weights, y_weights, correlations))
    }

    /// Compute objective function value
    fn compute_objective(
        &self,
        x_weights: &Array2<f64>,
        y_weights: &Array2<f64>,
        cxx: &Array2<f64>,
        cyy: &Array2<f64>,
        cxy: &Array2<f64>,
    ) -> GraphResult<f64> {
        // Simplified objective computation
        let correlation_term = x_weights.t().dot(cxy).dot(y_weights);
        let x_variance_term = x_weights.t().dot(cxx).dot(x_weights);
        let y_variance_term = y_weights.t().dot(cyy).dot(y_weights);

        let objective =
            correlation_term.sum() - 0.5 * x_variance_term.sum() - 0.5 * y_variance_term.sum();
        Ok(objective)
    }

    /// Compute graph regularization contribution
    fn compute_graph_regularization_value(
        &self,
        x_weights: &Array2<f64>,
        y_weights: &Array2<f64>,
    ) -> GraphResult<f64> {
        let mut reg_value = 0.0;

        if let Some(ref x_graph) = self.config.x_graph {
            let x_regularizer = self.compute_graph_regularizer(x_graph)?;
            let x_reg_contribution = x_weights.t().dot(&x_regularizer).dot(x_weights);
            reg_value += self.config.lambda * x_reg_contribution.sum();
        }

        if let Some(ref y_graph) = self.config.y_graph {
            let y_regularizer = self.compute_graph_regularizer(y_graph)?;
            let y_reg_contribution = y_weights.t().dot(&y_regularizer).dot(y_weights);
            reg_value += self.config.lambda * y_reg_contribution.sum();
        }

        Ok(reg_value)
    }
}

/// Network-constrained PLS implementation
pub struct NetworkConstrainedPLS {
    config: GraphRegularizationConfig,
}

impl NetworkConstrainedPLS {
    /// Create new network-constrained PLS
    pub fn new(config: GraphRegularizationConfig) -> Self {
        Self { config }
    }

    /// Fit network-constrained PLS model
    pub fn fit(&self, x: &Array2<f64>, y: &Array2<f64>) -> GraphResult<GraphRegularizationResults> {
        // For simplicity, use similar approach as graph-regularized CCA
        let gcca = GraphRegularizedCCA::new(self.config.clone());
        gcca.fit(x, y)
    }
}

/// Multi-graph CCA for multi-layer networks
pub struct MultiGraphCCA {
    config: GraphRegularizationConfig,
}

impl MultiGraphCCA {
    /// Create new multi-graph CCA
    pub fn new(config: GraphRegularizationConfig) -> Self {
        Self { config }
    }

    /// Fit multi-graph CCA with multiple graph layers
    pub fn fit_multi_layer(
        &self,
        x: &Array2<f64>,
        y: &Array2<f64>,
        x_graphs: &[GraphStructure],
        y_graphs: &[GraphStructure],
    ) -> GraphResult<GraphRegularizationResults> {
        // Combine multiple graph layers into a single regularizer
        let combined_x_regularizer = self.combine_graph_layers(x_graphs)?;
        let combined_y_regularizer = self.combine_graph_layers(y_graphs)?;

        // Create combined graph structures
        let combined_x_graph = GraphStructure {
            adjacency_matrix: combined_x_regularizer,
            graph_type: GraphType::MultiLayer,
            degrees: Array1::zeros(x.ncols()),
            communities: None,
            edge_weights: None,
            temporal_info: None,
            multi_layer_info: None,
        };

        let combined_y_graph = GraphStructure {
            adjacency_matrix: combined_y_regularizer,
            graph_type: GraphType::MultiLayer,
            degrees: Array1::zeros(y.ncols()),
            communities: None,
            edge_weights: None,
            temporal_info: None,
            multi_layer_info: None,
        };

        // Use graph-regularized CCA with combined graphs
        let mut config = self.config.clone();
        config.x_graph = Some(combined_x_graph);
        config.y_graph = Some(combined_y_graph);

        let gcca = GraphRegularizedCCA::new(config);
        gcca.fit(x, y)
    }

    /// Combine multiple graph layers into a single regularizer
    fn combine_graph_layers(&self, graphs: &[GraphStructure]) -> GraphResult<Array2<f64>> {
        if graphs.is_empty() {
            return Err(GraphRegularizationError::InvalidGraph(
                "No graphs provided for multi-layer combination".to_string(),
            ));
        }

        let n = graphs[0].adjacency_matrix.nrows();
        let mut combined = Array2::zeros((n, n));

        // Simple averaging of adjacency matrices
        for graph in graphs {
            if graph.adjacency_matrix.dim() != (n, n) {
                return Err(GraphRegularizationError::DimensionError(
                    "All graphs must have same dimensions".to_string(),
                ));
            }
            combined = combined + &graph.adjacency_matrix;
        }

        // Average the combined matrix
        combined = combined / graphs.len() as f64;

        Ok(combined)
    }
}

/// Helper functions for creating common graph structures
pub struct GraphBuilder;

impl GraphBuilder {
    /// Create a grid graph (lattice)
    pub fn grid_graph(rows: usize, cols: usize) -> GraphStructure {
        let n = rows * cols;
        let mut adjacency = Array2::zeros((n, n));

        for i in 0..rows {
            for j in 0..cols {
                let idx = i * cols + j;

                // Connect to neighbors
                if j > 0 {
                    // Left neighbor
                    let neighbor = i * cols + (j - 1);
                    adjacency[[idx, neighbor]] = 1.0;
                    adjacency[[neighbor, idx]] = 1.0;
                }
                if i > 0 {
                    // Top neighbor
                    let neighbor = (i - 1) * cols + j;
                    adjacency[[idx, neighbor]] = 1.0;
                    adjacency[[neighbor, idx]] = 1.0;
                }
            }
        }

        let degrees = adjacency.sum_axis(Axis(1));

        GraphStructure {
            adjacency_matrix: adjacency,
            graph_type: GraphType::Undirected,
            degrees,
            communities: None,
            edge_weights: None,
            temporal_info: None,
            multi_layer_info: None,
        }
    }

    /// Create a complete graph
    pub fn complete_graph(n: usize) -> GraphStructure {
        let mut adjacency = Array2::ones((n, n));

        // Remove self-loops
        for i in 0..n {
            adjacency[[i, i]] = 0.0;
        }

        let degrees = adjacency.sum_axis(Axis(1));

        GraphStructure {
            adjacency_matrix: adjacency,
            graph_type: GraphType::Undirected,
            degrees,
            communities: None,
            edge_weights: None,
            temporal_info: None,
            multi_layer_info: None,
        }
    }

    /// Create a random graph with specified edge probability
    pub fn random_graph(n: usize, edge_probability: f64) -> GraphStructure {
        let mut adjacency = Array2::zeros((n, n));

        for i in 0..n {
            for j in (i + 1)..n {
                if thread_rng().gen::<f64>() < edge_probability {
                    adjacency[[i, j]] = 1.0;
                    adjacency[[j, i]] = 1.0;
                }
            }
        }

        let degrees = adjacency.sum_axis(Axis(1));

        GraphStructure {
            adjacency_matrix: adjacency,
            graph_type: GraphType::Undirected,
            degrees,
            communities: None,
            edge_weights: None,
            temporal_info: None,
            multi_layer_info: None,
        }
    }

    /// Create a graph from distance matrix with threshold
    pub fn threshold_graph(distance_matrix: &Array2<f64>, threshold: f64) -> GraphStructure {
        let n = distance_matrix.nrows();
        let mut adjacency = Array2::zeros((n, n));

        for i in 0..n {
            for j in 0..n {
                if i != j && distance_matrix[[i, j]] <= threshold {
                    adjacency[[i, j]] = 1.0;
                }
            }
        }

        let degrees = adjacency.sum_axis(Axis(1));

        GraphStructure {
            adjacency_matrix: adjacency,
            graph_type: GraphType::Undirected,
            degrees,
            communities: None,
            edge_weights: None,
            temporal_info: None,
            multi_layer_info: None,
        }
    }

    /// Create k-nearest neighbors graph
    pub fn knn_graph(data: &Array2<f64>, k: usize) -> GraphStructure {
        let n = data.nrows();
        let mut adjacency = Array2::zeros((n, n));

        for i in 0..n {
            // Compute distances to all other points
            let mut distances: Vec<(usize, f64)> = Vec::new();
            for j in 0..n {
                if i != j {
                    let dist = Self::euclidean_distance(&data.row(i), &data.row(j));
                    distances.push((j, dist));
                }
            }

            // Sort by distance and take k nearest neighbors
            distances.sort_by(|a, b| a.1.partial_cmp(&b.1).unwrap());
            for (neighbor, _) in distances.iter().take(k) {
                adjacency[[i, *neighbor]] = 1.0;
            }
        }

        // Make symmetric
        for i in 0..n {
            for j in 0..n {
                if adjacency[[i, j]] > 0.0 || adjacency[[j, i]] > 0.0 {
                    adjacency[[i, j]] = 1.0;
                    adjacency[[j, i]] = 1.0;
                }
            }
        }

        let degrees = adjacency.sum_axis(Axis(1));

        GraphStructure {
            adjacency_matrix: adjacency,
            graph_type: GraphType::Undirected,
            degrees,
            communities: None,
            edge_weights: None,
            temporal_info: None,
            multi_layer_info: None,
        }
    }

    fn euclidean_distance(x: &ArrayView1<f64>, y: &ArrayView1<f64>) -> f64 {
        x.iter()
            .zip(y.iter())
            .map(|(xi, yi)| (xi - yi).powi(2))
            .sum::<f64>()
            .sqrt()
    }
}

#[allow(non_snake_case)]
#[cfg(test)]
mod tests {
    use super::*;
    use scirs2_core::ndarray::{Array1, Array2};

    #[test]
    fn test_graph_structure_creation() {
        let adj = Array2::eye(5);
        let degrees = Array1::ones(5);

        let graph = GraphStructure {
            adjacency_matrix: adj,
            graph_type: GraphType::Undirected,
            degrees,
            communities: None,
            edge_weights: None,
            temporal_info: None,
            multi_layer_info: None,
        };

        assert_eq!(graph.adjacency_matrix.dim(), (5, 5));
        assert_eq!(graph.graph_type, GraphType::Undirected);
    }

    #[test]
    fn test_graph_regularization_config() {
        let config = GraphRegularizationConfig::new(RegularizationType::GraphLaplacian, 0.5)
            .with_components(3)
            .with_parameter("test_param", 1.5);

        assert_eq!(
            config.regularization_type,
            RegularizationType::GraphLaplacian
        );
        assert_eq!(config.lambda, 0.5);
        assert_eq!(config.n_components, 3);
        assert_eq!(config.additional_params.get("test_param"), Some(&1.5));
    }

    #[test]
    fn test_graph_laplacian_computation() {
        let adj = scirs2_core::ndarray::arr2(&[[0.0, 1.0, 1.0], [1.0, 0.0, 1.0], [1.0, 1.0, 0.0]]);

        let config = GraphRegularizationConfig::default();
        let gcca = GraphRegularizedCCA::new(config);
        let laplacian = gcca.compute_graph_laplacian(&adj).unwrap();

        // Check dimensions
        assert_eq!(laplacian.dim(), (3, 3));

        // Check diagonal entries (should be degrees)
        assert_eq!(laplacian[[0, 0]], 2.0);
        assert_eq!(laplacian[[1, 1]], 2.0);
        assert_eq!(laplacian[[2, 2]], 2.0);

        // Check off-diagonal entries (should be negative adjacency)
        assert_eq!(laplacian[[0, 1]], -1.0);
        assert_eq!(laplacian[[1, 0]], -1.0);
    }

    #[test]
    fn test_graph_regularized_cca() {
        let x = Array2::from_shape_simple_fn((50, 4), || thread_rng().gen::<f64>());
        let y = Array2::from_shape_simple_fn((50, 3), || thread_rng().gen::<f64>());

        // Create simple graph structures
        let x_graph = GraphBuilder::complete_graph(4);
        let y_graph = GraphBuilder::complete_graph(3);

        let config = GraphRegularizationConfig::new(RegularizationType::GraphLaplacian, 0.1)
            .with_x_graph(x_graph)
            .with_y_graph(y_graph)
            .with_components(2);

        let gcca = GraphRegularizedCCA::new(config);
        let results = gcca.fit(&x, &y).unwrap();

        // Check output dimensions
        assert_eq!(results.x_weights.dim(), (4, 2));
        assert_eq!(results.y_weights.dim(), (3, 2));
        assert_eq!(results.correlations.len(), 2);
        assert!(results.final_objective.is_finite());
        assert!(results.graph_regularization_value >= 0.0);
    }

    #[test]
    fn test_network_constrained_pls() {
        let x = Array2::from_shape_simple_fn((30, 5), || thread_rng().gen::<f64>());
        let y = Array2::from_shape_simple_fn((30, 4), || thread_rng().gen::<f64>());

        let x_graph = GraphBuilder::grid_graph(5, 1);
        let y_graph = GraphBuilder::random_graph(4, 0.5);

        let config = GraphRegularizationConfig::new(RegularizationType::RandomWalk, 0.2)
            .with_x_graph(x_graph)
            .with_y_graph(y_graph);

        let npls = NetworkConstrainedPLS::new(config);
        let results = npls.fit(&x, &y).unwrap();

        assert_eq!(results.x_weights.dim(), (5, 2));
        assert_eq!(results.y_weights.dim(), (4, 2));
    }

    #[test]
    fn test_multi_graph_cca() {
        let x = Array2::from_shape_simple_fn((20, 3), || thread_rng().gen::<f64>());
        let y = Array2::from_shape_simple_fn((20, 3), || thread_rng().gen::<f64>());

        let x_graph1 = GraphBuilder::complete_graph(3);
        let x_graph2 = GraphBuilder::random_graph(3, 0.5);
        let y_graph1 = GraphBuilder::grid_graph(3, 1);
        let y_graph2 = GraphBuilder::threshold_graph(&Array2::ones((3, 3)), 0.5);

        let x_graphs = vec![x_graph1, x_graph2];
        let y_graphs = vec![y_graph1, y_graph2];

        let config = GraphRegularizationConfig::new(RegularizationType::Community, 0.15);
        let mgcca = MultiGraphCCA::new(config);

        let results = mgcca.fit_multi_layer(&x, &y, &x_graphs, &y_graphs).unwrap();

        assert_eq!(results.x_weights.dim(), (3, 2));
        assert_eq!(results.y_weights.dim(), (3, 2));
    }

    #[test]
    fn test_graph_builders() {
        // Test grid graph
        let grid = GraphBuilder::grid_graph(3, 3);
        assert_eq!(grid.adjacency_matrix.dim(), (9, 9));
        assert_eq!(grid.graph_type, GraphType::Undirected);

        // Test complete graph
        let complete = GraphBuilder::complete_graph(5);
        assert_eq!(complete.adjacency_matrix.dim(), (5, 5));
        assert_eq!(complete.degrees.sum(), 20.0); // n*(n-1) = 5*4

        // Test random graph
        let random = GraphBuilder::random_graph(6, 0.5);
        assert_eq!(random.adjacency_matrix.dim(), (6, 6));

        // Test kNN graph
        let data = Array2::from_shape_simple_fn((8, 2), || thread_rng().gen::<f64>());
        let knn = GraphBuilder::knn_graph(&data, 3);
        assert_eq!(knn.adjacency_matrix.dim(), (8, 8));
    }

    #[test]
    fn test_different_regularization_types() {
        let x = Array2::from_shape_simple_fn((25, 3), || thread_rng().gen::<f64>());
        let y = Array2::from_shape_simple_fn((25, 3), || thread_rng().gen::<f64>());
        let graph = GraphBuilder::complete_graph(3);

        let regularization_types = [
            RegularizationType::GraphLaplacian,
            RegularizationType::RandomWalk,
            RegularizationType::DiffusionKernel,
            RegularizationType::Community,
            RegularizationType::GraphNeuralNetwork,
        ];

        for &reg_type in &regularization_types {
            let config = GraphRegularizationConfig::new(reg_type, 0.1)
                .with_x_graph(graph.clone())
                .with_y_graph(graph.clone());

            let gcca = GraphRegularizedCCA::new(config);
            let results = gcca.fit(&x, &y);

            assert!(
                results.is_ok(),
                "Failed for regularization type {:?}",
                reg_type
            );
            let results = results.unwrap();
            assert_eq!(results.x_weights.dim(), (3, 2));
            assert_eq!(results.y_weights.dim(), (3, 2));
        }
    }
}