ratio-markov 0.5.1

//! # Markov chain steady-state flow calculations module
//!
//! ## License
//!
//! This Source Code Form is subject to the terms of the Mozilla Public License, v. 2.0.
//! If a copy of the MPL was not distributed with this file,
//! You can obtain one at <https://mozilla.org/MPL/2.0/>.
//!
//! **Code examples both in the docstrings and rendered documentation are free to use.**

pub use crate::{DMatrix, DVector};

/// MarkovFlow struct contains all results of the Markov flow calculations based on an
/// adjacency matrix, evaporation constant, and scaling toggle.
///
/// # Example
/// ```
/// extern crate nalgebra as na;
/// use ratio_markov::{DMatrix, DVector, MarkovFlow};
/// let adj = DMatrix::from_row_slice(3, 3, &[2.0, 0.5, 0.7, 1.0, 2.0, 3.0, 1.0, 2.6, 0.1]);
/// let markov = MarkovFlow::new(&adj, 4.0, false);
/// ```
#[derive(Debug)]
pub struct MarkovFlow {
    /// Adjacency matrix with only positive values. Negatives are set to 0.0.
    pub adjacency_matrix: DMatrix<f64>,
    /// Dimension of this problem (size of adjacency matrix axes).
    pub dim: usize,
    /// Sink matrix. A probability matrix that has been provided with a evaporation sink.
    pub sink_matrix: DMatrix<f64>,
    /// Sensitivity matrix that shows how much flow from a source node ends up in another node if it were injected with 1.0.
    pub sensitivity_matrix: DMatrix<f64>,
    /// Influence matrix, percentage of outflow through j ends up in i.
    pub influence_matrix: DMatrix<f64>,
    /// Origin matrix that shows how much flow originates from a source node in a target node if its flow is 1.0.
    pub origin_matrix: DMatrix<f64>,
    /// Dependency matrix, percentage of inflow through j coming from i.
    pub dependency_matrix: DMatrix<f64>,
    /// Input vector (ones or scaled).
    pub in_vector: DVector<f64>,
    /// Flow vector through each node w.r.t. input vector.
    pub flow_vector: DVector<f64>,
    /// Flow matrix denoting the resulting flow over every edge when the system is injected with in_vector.
    pub flow_matrix: DMatrix<f64>,
}

impl MarkovFlow {
    /// New MarkovFlow struct containing all Markov chain calculations.
    ///
    /// # Arguments
    ///
    /// * `adj` - Square adjacency matrix with nodes in identical order on both axis.
    ///   Only positive adjacency values are used in the calculations.
    /// * `mu` - Evaporation constant. Factor by which flow is reduced after passing
    ///   through a node.
    /// * `scale` - Whether injected flow in each node should be scaled according to the
    ///   output adjacency of a node.
    pub fn new(adj: &DMatrix<f64>, mu: f64, scale: bool) -> MarkovFlow {
        let (rows, cols) = adj.shape();
        if rows != cols {
            panic!("Adjacency matrix has to be square.");
        }
        // No negative values beyond this point.
        let mut adjacency_matrix = adj.map(|x: f64| if x < 0.0 { 0.0 } else { x });
        adjacency_matrix.fill_diagonal(0.0);

        let sink_matrix = get_sink_matrix(&adjacency_matrix, mu);

        let sensitivity_matrix = get_sensitivity_matrix(&sink_matrix);
        let influence_matrix = sensitivity_matrix.view((0, 0), (rows, cols)).clone_owned();

        let origin_matrix = get_sensitivity_matrix(&sink_matrix.transpose());
        let dependency_matrix = origin_matrix.view((0, 0), (rows, cols)).clone_owned();

        let in_vector = get_in_vector(&adjacency_matrix, scale);
        let flow_vector = get_flow_vector(&sensitivity_matrix, &in_vector);
        let flow_matrix = get_flow_matrix(&sink_matrix, &in_vector);

        MarkovFlow {
            adjacency_matrix,
            dim: rows,
            sink_matrix,
            sensitivity_matrix,
            influence_matrix,
            origin_matrix,
            dependency_matrix,
            in_vector,
            flow_vector,
            flow_matrix,
        }
    }
}

/// get_sink_matrix calculates a column normalized matrix with an added evaporation sink
/// as the last entry. Input adjacency matrix 'adj' should be square and 'mu' should be a
/// number > 1.0.
pub fn get_sink_matrix(adj: &DMatrix<f64>, mu: f64) -> DMatrix<f64> {
    let cont = 1.0 / mu;
    let evap = 1.0 - cont;
    let (rows, cols) = adj.shape();
    if mu <= 1.0 {
        panic!("Mu should be >1!");
    }
    let mut sink_matrix = DMatrix::zeros(rows + 1, cols + 1);

    // Normalize according to the sums of each adjacency matrix column.
    // Also, multiply by the fraction that should continue (not evaporate).
    for j in 0..cols {
        let sum = adj.column(j).sum();
        // No outputs, only flow continuation through self.
        if sum.abs() < f64::EPSILON {
            sink_matrix[(j, j)] = cont;
        } else {
            // Output cell, normalize!
            let scale = cont / sum;
            for i in 0..rows {
                sink_matrix[(i, j)] = adj[(i, j)] * scale;
            }
        }
        // Set the evaporation row element.
        sink_matrix[(rows, j)] = evap;
    }

    sink_matrix
}

/// get_sensitivity_matrix calculates the sensitivity matrix based on the sink matrix.
pub fn get_sensitivity_matrix(sink_matrix: &DMatrix<f64>) -> DMatrix<f64> {
    let (rows, cols) = sink_matrix.shape();
    let mut sensitivity_matrix = DMatrix::identity(rows, cols) - sink_matrix;
    let inverted = sensitivity_matrix.try_inverse_mut();
    if !inverted {
        panic!("Could not invert the given matrix.");
    }
    sensitivity_matrix
}

/// get_in_vector calculates the input vector based on the shape of the adjacency matrix
/// and optional scaling. Without scaling, all nodes except the final one (sink) of the
/// sink matrix should be supplied with 1.0 input flow. With scaling, input flow equals
/// the column sum of the adjacency matrix.
pub fn get_in_vector(adjacency_matrix: &DMatrix<f64>, scale: bool) -> DVector<f64> {
    let rows = adjacency_matrix.shape().0;
    let mut in_vector = DVector::from_element(rows + 1, 1.0);
    if scale {
        let sums = adjacency_matrix.row_sum();
        for idx in 0..(rows) {
            in_vector[idx] = sums[idx];
        }
    } else {
        in_vector[rows] = 0.0;
    }
    in_vector
}

/// get_flow_vector calculates the flow over every node according to the input vector and
/// sensitivity matrix. It equals the (outer) product of the sensitivity matrix with the
/// input vector.
pub fn get_flow_vector(
    sensitivity_matrix: &DMatrix<f64>,
    in_vector: &DVector<f64>,
) -> DVector<f64> {
    sensitivity_matrix * in_vector
}

/// get_flow_matrix calculates the flow matrix after the system is injected with an
/// input vector. Each entry denotes the steady state flow over the corresponding edge.
pub fn get_flow_matrix(sink_matrix: &DMatrix<f64>, in_vector: &DVector<f64>) -> DMatrix<f64> {
    let (rows, cols) = sink_matrix.shape();
    DMatrix::from_fn(rows, cols, |i, j| sink_matrix[(i, j)] * in_vector[j])
}

#[cfg(test)]
mod tests {
    use super::*;
    #[test]
    fn example() {
        let adj = DMatrix::from_row_slice(3, 3, &[2.0, 0.5, 0.7, 1.0, 2.0, 3.0, 1.0, 2.6, 0.1]);
        let flow = MarkovFlow::new(&adj, 4.0, false);
        // Check the sink matrix.
        assert_eq!(flow.sink_matrix[(0, 0)], 0.0);
        assert_eq!(flow.sink_matrix[(3, 0)], 0.75);
        let reference = DMatrix::from_row_slice(
            4,
            1,
            &[1.12594131394443, 1.43287457803168, 1.4411841080238899, 3.0],
        );
        dbg!(&flow.flow_vector);
        assert!(
            flow.flow_vector
                .relative_eq(&reference, 0.000000000000001, 0.000000000000001)
        );
    }

    #[test]
    fn climate_control() {
        let adj = DMatrix::from_row_slice(
            16,
            16,
            &[
                0., 0., 0., 0., 0., 0., 6., 0., 0., 0., 6., 0., -1., 0., 2., 0., 0., 0., 4., 0.,
                0., 2., 0., 0., 0., 4., 0., 0., 0., 0., 0., 0., 0., 4., 0., 4., 0., 2., 4., 0., 0.,
                0., 0., 0., 0., 0., 4., 4., 0., 0., 4., 0., 8., 2., 0., 0., 0., 4., 0., 0., 0., 0.,
                0., 0., 0., 0., 0., 8., 0., 2., 0., 0., 0., 8., 4., 4., 0., 0., 0., 0., 0., 2., 2.,
                2., 2., 0., 2., 0., 2., 0., 0., 0., 0., 0., 2., 2., 6., 0., 4., 0., 0., 2., 0., 8.,
                0., 0., 8., 0., 0., 0., 4., 0., 0., 0., 0., 0., 0., 0., 8., 0., 8., 0., 4., 0., 0.,
                0., 0., 0., 0., 0., 0., 0., 0., 2., 0., 8., 0., 0., 0., 0., 0., 8., 0., 0., 0., 4.,
                0., 4., 8., 0., 0., 0., 0., 0., 4., 4., 0., 0., 0., 0., 6., 0., 0., 0., 4., 0., 8.,
                4., 0., 4., 0., -1., 0., 0., 0., 0., 0., 0., 0., 0., 4., 0., 0., 0., 0., 4., -1.,
                0., 2., 0., 0., 0., -1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 2., 0., 0., 0.,
                0., 0., 0., 0., 0., 0., 0., 0., 0., 8., 0., 0., 0., 0., 0., 0., 0., 2., 0., 4., 0.,
                0., 2., 4., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 4., 0., 0., 2., 0., 0., 0.,
                0., 0., 0., 0., 0., 0., 0.,
            ],
        );
        let flow = MarkovFlow::new(&adj, 2.0, false);
        let flow_vec_ref = DMatrix::from_row_slice(
            17,
            1,
            &[
                1.66245393, 1.57043341, 2.5625501, 1.96454112, 2.56170014, 2.05911155, 2.75477062,
                2.03070853, 2.2891907, 2.50535804, 2.309491, 2.00614073, 1.20061407, 1.50870904,
                1.65257344, 1.36165357, 16.,
            ],
        );
        dbg!(&flow.flow_vector);
        assert!(
            flow.flow_vector
                .relative_eq(&flow_vec_ref, 0.00000001, 0.00000001)
        );
    }
}