ultrametric_matrix_tools/
utils.rs

1//! Usful funtions to construct and check propertyies of ultrametric matrices.
2
3use nalgebra::DMatrix;
4use ndarray::prelude::*;
5use numpy::{IntoPyArray, PyArray2, PyReadonlyArrayDyn};
6use pyo3::prelude::*;
7use rand::prelude::*;
8
9#[pymodule]
10fn utils(_py: Python, m: &PyModule) -> PyResult<()> {
11    #[pyfn(m)]
12    #[pyo3(name = "is_ultrametric")]
13    pub fn is_ultrametric_py(py_matrix: PyReadonlyArrayDyn<f64>) -> bool {
14        let size = py_matrix.shape()[0];
15        let py_array = py_matrix.as_array();
16        let mut matrix = DMatrix::<f64>::zeros(size, size);
17        for i in 0..size {
18            for j in 0..size {
19                matrix[(i, j)] = py_array[[i, j]];
20            }
21        }
22        return self::is_ultrametric(&matrix);
23    }
24    #[pyfn(m)]
25    #[pyo3(name = "random_ultrametric_matrix")]
26    pub fn random_ultrametric_matrix_py<'py>(py: Python<'py>, size: usize) -> &'py PyArray2<f64> {
27        let matrix = random_ultrametric_matrix(size);
28        let mut py_matrix = Array2::zeros((size, size));
29        for i in 0..size {
30            for j in 0..size {
31                py_matrix[[i, j]] = matrix[(i, j)];
32            }
33        }
34        return py_matrix.into_pyarray(py);
35    }
36
37    #[pyfn(m)]
38    #[pyo3(name = "random_special_ultrametric_matrix")]
39    pub fn random_special_ultrametric_matrix_py<'py>(
40        py: Python<'py>,
41        size: usize,
42    ) -> &'py PyArray2<f64> {
43        let matrix = random_special_ultrametric_matrix(size);
44        let mut py_matrix = Array2::zeros((size, size));
45        for i in 0..size {
46            for j in 0..size {
47                py_matrix[[i, j]] = matrix[(i, j)];
48            }
49        }
50        return py_matrix.into_pyarray(py);
51    }
52    Ok(())
53}
54
55/// Checks if a matrix is ultrametric.
56///
57/// To check the ultrametric property for a matrix, a similar algorithm as `from_matrix` for `UltrametricTree` is used. First, the symmetry is checked. If the matrix is symmetric, then it is recursively checked if a permutation of the matrix can be partitioned such that the off-diagonal blocks have one value and the diagonal blocks are ultrametric.
58///
59/// # Example:
60/// ```
61/// let ultrametric = ultrametric_matrix_tools::na::DMatrix::from_vec(4, 4,
62///     vec![0.0, 1.0, 3.0, 1.0, 1.0, 3.0, 1.0, 1.0, 3.0, 1.0, 5.0, 1.0, 1.0, 1.0, 1.0, 1.0]);
63/// let not_ultrametric = ultrametric_matrix_tools::na::DMatrix::from_vec(4, 4,
64///     vec![0.0, 1.0, 3.0, 2.0, 1.0, 3.0, 1.0, 1.0, 3.0, 1.0, 5.0, 1.0, 2.0, 1.0, 1.0, 1.0]);
65///
66/// assert_eq!(ultrametric_matrix_tools::utils::is_ultrametric(&ultrametric), true);
67/// assert_eq!(ultrametric_matrix_tools::utils::is_ultrametric(&not_ultrametric), false);
68/// ```
69pub fn is_ultrametric(matrix: &DMatrix<f64>) -> bool {
70    if !is_symmetric(matrix) {
71        return false;
72    }
73    let idx: Vec<usize> = (0..matrix.nrows()).collect();
74    if !is_submatrix_ultrametric(matrix, idx, 0.0) {
75        return false;
76    }
77    return true;
78}
79
80fn is_symmetric(matrix: &DMatrix<f64>) -> bool {
81    if !matrix.is_square() {
82        return false;
83    }
84    let size = matrix.shape().0;
85    for i in 1..size {
86        for j in (i + 1)..size {
87            if matrix[(i, j)] != matrix[(j, i)] {
88                return false;
89            }
90        }
91    }
92    return true;
93}
94
95fn is_submatrix_ultrametric(matrix: &DMatrix<f64>, idx: Vec<usize>, prev_value: f64) -> bool {
96    let first_i = idx[0];
97    if idx.len() == 1 {
98        return true;
99    }
100    let mut left_partition: Vec<usize> = vec![first_i];
101    let mut right_partition: Vec<usize> = Vec::new();
102    let mut min = f64::MAX;
103    for &i in &idx[1..] {
104        if min > matrix[(first_i, i)] {
105            min = matrix[(first_i, i)];
106            left_partition.extend(right_partition.iter());
107            right_partition.clear();
108            right_partition.push(i);
109        } else {
110            if min == matrix[(first_i, i)] {
111                right_partition.push(i);
112            } else {
113                left_partition.push(i);
114            }
115        }
116    }
117    if min < prev_value || !is_block_equal(matrix, &left_partition, &right_partition, min) {
118        return false;
119    }
120
121    let left = is_submatrix_ultrametric(matrix, left_partition, min);
122    let right = is_submatrix_ultrametric(matrix, right_partition, min);
123    return left && right;
124}
125
126fn is_block_equal(
127    matrix: &DMatrix<f64>,
128    left: &Vec<usize>,
129    right: &Vec<usize>,
130    value: f64,
131) -> bool {
132    for &i in left.iter() {
133        for &j in right.iter() {
134            if matrix[(i, j)] != value {
135                return false;
136            }
137        }
138    }
139    return true;
140}
141
142/// Constructs a random ultrametric matrix
143///
144/// The elements of the ultrametric matrix have an integer value between `1` and `size - 1`.
145///
146/// The construction is based on Theorem 1.1 in [(Fiedler, 2002)](fn.random_ultrametric_matrix.html#fn1)[^Fiedler, 2002]. The off-diagonal elements are constructed using the rules in Theorem 1.1. The diagonal elements are randomly chosen from `1` to `size - 1`.
147///
148/// # Example:
149/// ```
150/// let ultrametric_matrix = ultrametric_matrix_tools::utils::random_ultrametric_matrix(10);
151///
152/// assert_eq!(ultrametric_matrix_tools::utils::is_ultrametric(&ultrametric_matrix), true);
153/// ```
154///
155/// [^Fiedler, 2002]: [Fiedler, M., 2002. Remarks on Monge matrices. Mathematica Bohemica, 127(1), pp.27-32.](https://dml.cz/bitstream/handle/10338.dmlcz/133983/MathBohem_127-2002-1_3.pdf)
156pub fn random_ultrametric_matrix(size: usize) -> DMatrix<f64> {
157    let mut matrix = DMatrix::<f64>::zeros(size, size);
158    let mut rng = thread_rng();
159    for i in 1..size {
160        let elem = rng.gen_range(1..size) as f64;
161        matrix[(i - 1, i)] = elem;
162        matrix[(i, i - 1)] = elem;
163    }
164    for i in (0..(size - 2)).rev() {
165        for k in (i + 2)..size {
166            let elem = f64::min(matrix[(i, k - 1)], matrix[(i + 1, k)]);
167            matrix[(i, k)] = elem;
168            matrix[(k, i)] = elem;
169        }
170    }
171    let diag = random_vector(size);
172    for i in 0..size {
173        matrix[(i, i)] = diag[i];
174    }
175    permutate_matrix(&mut matrix);
176    return matrix;
177}
178
179/// Constructs a random special ultrametric matrix
180///
181/// The elements of the ultrametric matrix have an integer value between `1` and `size - 1`.
182///
183/// The construction is based on Theorem 1.1 in [(Fiedler, 2002)](fn.random_ultrametric_matrix.html#fn1)[^Fiedler, 2002]. The elements are constructed using the rules in Theorem 1.1.
184///
185/// # Example:
186/// ```
187/// let ultrametric_matrix = ultrametric_matrix_tools::utils::random_special_ultrametric_matrix(10);
188///
189/// assert_eq!(ultrametric_matrix_tools::utils::is_ultrametric(&ultrametric_matrix), true);
190/// ```
191///
192/// [^Fiedler, 2002]: [Fiedler, M., 2002. Remarks on Monge matrices. Mathematica Bohemica, 127(1), pp.27-32.](https://dml.cz/bitstream/handle/10338.dmlcz/133983/MathBohem_127-2002-1_3.pdf)
193pub fn random_special_ultrametric_matrix(size: usize) -> DMatrix<f64> {
194    let mut matrix = DMatrix::<f64>::zeros(size, size);
195    let mut rng = thread_rng();
196    for i in 1..size {
197        let elem = rng.gen_range(1..size) as f64;
198        matrix[(i - 1, i)] = elem;
199        matrix[(i, i - 1)] = elem;
200    }
201    for i in (0..(size - 2)).rev() {
202        for k in (i + 2)..size {
203            let elem = f64::min(matrix[(i, k - 1)], matrix[(i + 1, k)]);
204            matrix[(i, k)] = elem;
205            matrix[(k, i)] = elem;
206        }
207    }
208    for i in 1..(size - 1) {
209        matrix[(i, i)] = f64::max(matrix[(i - 1, i)], matrix[(i, i + 1)]);
210    }
211    matrix[(0, 0)] = matrix[(0, 1)];
212    matrix[(size - 1, size - 1)] = matrix[(size - 2, size - 1)];
213    permutate_matrix(&mut matrix);
214    return matrix;
215}
216
217fn random_vector(size: usize) -> Vec<f64> {
218    let mut rng = thread_rng();
219    let mut vector: Vec<f64> = Vec::new();
220    for _ in 0..size {
221        vector.push(rng.gen_range(1..size) as f64);
222    }
223    return vector;
224}
225
226fn permutate_matrix(matrix: &mut DMatrix<f64>) {
227    let size = matrix.shape().0;
228    let mut new_matrix = DMatrix::<f64>::zeros(size, size);
229    let mut rng: StdRng = SeedableRng::seed_from_u64(42);
230    let mut element_vector: Vec<usize> = (0..size).collect();
231    element_vector.shuffle(&mut rng);
232    for (old_i, &new_i) in element_vector.iter().enumerate() {
233        for (old_j, &new_j) in element_vector.iter().enumerate() {
234            new_matrix[(new_i, new_j)] = matrix[(old_i, old_j)];
235        }
236    }
237    *matrix = new_matrix;
238}
ultrametric_matrix_tools/utils.rs

ultrametric_matrix_tools/
utils.rs
