hypergraphx 0.0.5

use crate::prelude::*;
use nalgebra::{Const, CsMatrix, DMatrix, Dyn, Matrix};

/// The root trait for all algebraic operations on graphs.
/// Provides methods to convert graphs into various matrix representations.
pub trait MatrixRepresentation<'a, D: GraphType>: GraphBasics<'a> {
    /// A |V| x |E| binary incidence matrix, where |V| is the number of nodes and |E| is the number of edges.
    /// If M_{ij} = 1, then node i is contained in edge j.
    fn binary_incidence_matrix(&'a self) -> CsMatrix<i8>;
    /// A standard |V| x |V| adjacency matrix, where |V| is the number of nodes.
    /// This does not preserve hypergraph structure, and cannot be transformed back into the original hypergraph.
    ///
    /// Note: Adjacency tensors are defined for uniform hypergraphs, and are implemented separately.
    fn adjacency_matrix(&'a self) -> DMatrix<usize>;
    /// Diag(degree_sequence) - adjacency_matrix
    fn laplacian_matrix(&'a self) -> DMatrix<usize>;
    /// The default implementation of this method calls `into_adjacency_matrix` and transposes the result.
    /// All graph types provided by this library implement this method separately, without
    /// the overhead of `transpose`.
    fn dual_adjacency_matrix(&'a self) -> DMatrix<usize> {
        self.adjacency_matrix().transpose()
    }
    fn dual_laplacian_matrix(&'a self) -> DMatrix<usize> {
        self.laplacian_matrix().transpose()
    }

    /// Eigenvalues of the adjacency matrix.
    fn alphas(&'a self) -> Vec<f64> {
        self.adjacency_matrix()
            .map(|x| x as f64)
            .eigenvalues()
            .unwrap()
            .into_iter()
            .map(|x| *x)
            .collect()
    }

    /// Eigenvalues of the laplacian matrix.
    fn lambdas(&'a self) -> Vec<f64> {
        self.laplacian_matrix()
            .map(|x| x as f64)
            .eigenvalues()
            .unwrap()
            .into_iter()
            .map(|x| *x)
            .collect()
    }
}

impl<'a, T> MatrixRepresentation<'a, Undirected> for T
where
    T: GraphBasics<'a> + GraphProperties<'a>,
{
    fn binary_incidence_matrix(&'a self) -> CsMatrix<i8> {
        let mut irows = Vec::new();
        let mut icols = Vec::new();

        for edge_index in 0..self.edge_count() {
            for &node_index in &self.contained_nodes(edge_index.into()).unwrap() {
                irows.push(node_index.into());
                icols.push(edge_index);
            }
        }

        let matrix = CsMatrix::from_triplet(
            self.node_count(),
            self.edge_count(),
            &irows,
            &icols,
            &vec![1; irows.len()],
        );
        matrix
    }

    fn adjacency_matrix(&'a self) -> DMatrix<usize> {
        let l = self.node_count();
        let mut mat = DMatrix::<usize>::zeros(l, l);

        for (node_index, _) in self.nodes().enumerate() {
            let n = self.neighbours(node_index.into()).unwrap().len();
            for neighbour in 0..n {
                mat[(node_index, neighbour)] += 1;
            }
        }

        mat
    }

    fn dual_adjacency_matrix(&'a self) -> DMatrix<usize> {
        let l = self.node_count();
        let mut mat = DMatrix::<usize>::zeros(l, l);

        for (node_index, _) in self.nodes().enumerate() {
            let n = self.neighbours(node_index.into()).unwrap().len();
            for neighbour in 0..n {
                mat[(neighbour, node_index)] += 1;
            }
        }

        mat
    }

    fn laplacian_matrix(&'a self) -> DMatrix<usize> {
        let dseq = Matrix::<usize, Dyn, Const<1>, _>::from(self.degree_sequence());
        let out = DMatrix::<usize>::from_diagonal(&dseq);

        return out - self.adjacency_matrix();
    }

    fn dual_laplacian_matrix(&'a self) -> DMatrix<usize> {
        let dseq = Matrix::<usize, Dyn, Const<1>, _>::from(self.degree_sequence());
        let out = DMatrix::<usize>::from_diagonal(&dseq);

        return out - self.dual_adjacency_matrix();
    }
}

impl<'a, T> MatrixRepresentation<'a, Directed> for T
where
    T: GraphBasics<'a> + DiGraphProperties<'a>,
{
    fn binary_incidence_matrix(&'a self) -> CsMatrix<i8> {
        let mut irows = Vec::new();
        let mut icols = Vec::new();
        let mut vals = Vec::new();

        for edge_index in 0..self.edge_count() {
            for &node_index in self.source(edge_index.into()).unwrap().iter() {
                irows.push(node_index.into());
                icols.push(edge_index);
                vals.push(1);
            }
            for &node_index in self.target(edge_index.into()).unwrap().iter() {
                irows.push(node_index.into());
                icols.push(edge_index);
                vals.push(-1);
            }
        }

        let matrix =
            CsMatrix::from_triplet(self.node_count(), self.edge_count(), &irows, &icols, &vals);
        matrix
    }
    fn adjacency_matrix(&'a self) -> DMatrix<usize> {
        let l = self.node_count();
        let mut mat = DMatrix::<usize>::zeros(l, l);

        for node_index in 0..self.node_count() {
            let n = self.out_neighbours(node_index.into()).unwrap();
            for neighbour in n {
                mat[(node_index, neighbour.into())] += 1;
            }
        }

        mat
    }

    fn dual_adjacency_matrix(&'a self) -> DMatrix<usize> {
        let l = self.node_count();
        let mut mat = DMatrix::<usize>::zeros(l, l);

        for (node_index, _) in self.nodes().enumerate() {
            let n = self.out_neighbours(node_index.into()).unwrap();
            for neighbour in n {
                mat[(neighbour.into(), node_index)] += 1;
            }
        }

        mat
    }
    fn laplacian_matrix(&'a self) -> DMatrix<usize> {
        let dseq = Matrix::<usize, Dyn, Const<1>, _>::from(self.out_degree_sequence());
        let out = DMatrix::<usize>::from_diagonal(&dseq);

        return out - self.adjacency_matrix();
    }

    fn dual_laplacian_matrix(&'a self) -> DMatrix<usize> {
        let dseq = Matrix::<usize, Dyn, Const<1>, _>::from(self.out_degree_sequence());
        let out = DMatrix::<usize>::from_diagonal(&dseq);

        return out - self.dual_adjacency_matrix();
    }
}