microsoft/qdk
Publicmirrored from https://github.com/microsoft/qdkAvailable
source/paulimer/tests/bitmatrix_test.rs
422lines · modecode
| 1 | // Copyright (c) Microsoft Corporation. |
| 2 | // Licensed under the MIT License. |
| 3 | |
| 4 | use paulimer::bits::bitmatrix::{directly_summed, kernel_basis_matrix, rref_with_transforms}; |
| 5 | use paulimer::bits::tiny_matrix::{tiny_matrix_from_bitmatrix, tiny_matrix_rref}; |
| 6 | use paulimer::bits::{BitMatrix, BitVec, Bitwise, BitwiseBinaryOps, WORD_COUNT_DEFAULT}; |
| 7 | use proptest::prelude::*; |
| 8 | use rand::prelude::*; |
| 9 | use rustc_hash::FxHashSet; |
| 10 | use sorted_iter::assume::AssumeSortedByItemExt; |
| 11 | use sorted_iter::SortedIterator; |
| 12 | use std::collections::BTreeMap; |
| 13 | use std::str::FromStr; |
| 14 | |
| 15 | proptest! { |
| 16 | #[test] |
| 17 | fn shape(rowcount in 0..100usize, columncount in 0..100usize) { |
| 18 | let matrix = BitMatrix::<WORD_COUNT_DEFAULT>::with_shape(rowcount, columncount); |
| 19 | assert_eq!(matrix.rowcount(), rowcount); |
| 20 | assert_eq!(matrix.columncount(), columncount); |
| 21 | assert_eq!(matrix.shape(), (rowcount, columncount)); |
| 22 | } |
| 23 | |
| 24 | #[test] |
| 25 | fn zeros(rowcount in 0..100usize, columncount in 0..100usize) { |
| 26 | let matrix = BitMatrix::<WORD_COUNT_DEFAULT>::zeros(rowcount, columncount); |
| 27 | for irow in 0..matrix.rowcount() { |
| 28 | for icol in 0..matrix.columncount() { |
| 29 | assert!(!matrix[(irow, icol)]); |
| 30 | } |
| 31 | } |
| 32 | } |
| 33 | |
| 34 | #[test] |
| 35 | fn indexing(matrix in arbitrary_bitmatrix(100)) { |
| 36 | for irow in 0..matrix.rowcount() { |
| 37 | for icol in 0..matrix.columncount() { |
| 38 | assert_eq!(matrix[(irow, icol)], matrix[[irow, icol]]); |
| 39 | } |
| 40 | } |
| 41 | } |
| 42 | |
| 43 | #[test] |
| 44 | fn clone(matrix in arbitrary_bitmatrix(100)) { |
| 45 | assert_eq!(matrix, matrix.clone()); |
| 46 | } |
| 47 | |
| 48 | #[test] |
| 49 | fn swap_rows(matrix in nonempty_bitmatrix(100), raw_row_indexes in (0..100usize, 0..100usize)) { |
| 50 | let row_indexes = [raw_row_indexes.0 % matrix.rowcount(), raw_row_indexes.1 % matrix.rowcount()]; |
| 51 | let mut swapped = matrix.clone(); |
| 52 | swapped.swap_rows(row_indexes[0], row_indexes[1]); |
| 53 | for column_index in 0..matrix.columncount() { |
| 54 | assert_eq!(matrix[[row_indexes[0], column_index]], swapped[[row_indexes[1], column_index]]); |
| 55 | } |
| 56 | let row_indexes = row_indexes.into_iter().collect::<rustc_hash::FxHashSet<usize>>(); |
| 57 | for row_index in (0..matrix.rowcount()).collect::<FxHashSet<usize>>().difference(&row_indexes) { |
| 58 | for column_index in 0..matrix.columncount() { |
| 59 | assert_eq!(matrix[[*row_index, column_index]], swapped[[*row_index, column_index]]); |
| 60 | } |
| 61 | } |
| 62 | } |
| 63 | |
| 64 | #[test] |
| 65 | fn swap_columns(matrix in nonempty_bitmatrix(100), raw_column_indexes in (0..100usize, 0..100usize)) { |
| 66 | let column_indexes = [raw_column_indexes.0 % matrix.columncount(), raw_column_indexes.1 % matrix.columncount()]; |
| 67 | let mut swapped = matrix.clone(); |
| 68 | swapped.swap_columns(column_indexes[0], column_indexes[1]); |
| 69 | for row_index in 0..matrix.rowcount() { |
| 70 | assert_eq!(matrix[[row_index, column_indexes[0]]], swapped[[row_index, column_indexes[1]]]); |
| 71 | } |
| 72 | let column_indexes = column_indexes.into_iter().collect::<rustc_hash::FxHashSet<usize>>(); |
| 73 | for column_index in (0..matrix.columncount()).collect::<FxHashSet<usize>>().difference(&column_indexes) { |
| 74 | for row_index in 0..matrix.rowcount() { |
| 75 | assert_eq!(matrix[[row_index, *column_index]], swapped[[row_index, *column_index]]); |
| 76 | } |
| 77 | } |
| 78 | } |
| 79 | |
| 80 | #[test] |
| 81 | fn addition((left, right) in equal_shape_bitmatrices(100)) { |
| 82 | let sum = &left + &right; |
| 83 | for irow in 0..left.rowcount() { |
| 84 | for icol in 0..right.columncount() { |
| 85 | let index = (irow, icol); |
| 86 | assert_eq!(sum[index], left[index] ^ right[index]); |
| 87 | } |
| 88 | } |
| 89 | assert_eq!(sum, &right + &left); |
| 90 | } |
| 91 | |
| 92 | #[test] |
| 93 | fn addition_inplace((mut left, right) in equal_shape_bitmatrices(100)) { |
| 94 | let sum = &left + &right; |
| 95 | left += &right; |
| 96 | assert_eq!(sum, left); |
| 97 | } |
| 98 | |
| 99 | #[test] |
| 100 | fn xor((left, right) in equal_shape_bitmatrices(100)) { |
| 101 | assert_eq!(&left ^ &right, &left + &right); |
| 102 | } |
| 103 | |
| 104 | #[test] |
| 105 | fn xor_inplace((mut left, right) in equal_shape_bitmatrices(100)) { |
| 106 | let xor = &left ^ &right; |
| 107 | left ^= &right; |
| 108 | assert_eq!(xor, left); |
| 109 | } |
| 110 | |
| 111 | #[test] |
| 112 | fn and((left, right) in equal_shape_bitmatrices(100)) { |
| 113 | let and = &left & &right; |
| 114 | for irow in 0..left.rowcount() { |
| 115 | for icol in 0..left.columncount() { |
| 116 | let index = (irow, icol); |
| 117 | assert_eq!(and[index], left[index] & right[index]); |
| 118 | } |
| 119 | } |
| 120 | assert_eq!(and, &right & &left); |
| 121 | } |
| 122 | |
| 123 | |
| 124 | #[test] |
| 125 | fn and_inplace((mut left, right) in equal_shape_bitmatrices(100)) { |
| 126 | let and = &left & &right; |
| 127 | left &= &right; |
| 128 | assert_eq!(and, left); |
| 129 | } |
| 130 | |
| 131 | #[test] |
| 132 | fn equality(left in arbitrary_bitmatrix(100), right in arbitrary_bitmatrix(100)) { |
| 133 | let mut are_equal = left.shape() == right.shape(); |
| 134 | if are_equal { |
| 135 | for irow in 0..left.rowcount() { |
| 136 | for icol in 0..right.columncount() { |
| 137 | let index = (irow, icol); |
| 138 | are_equal &= left[index] == right[index]; |
| 139 | } |
| 140 | } |
| 141 | } |
| 142 | assert_eq!(left == right, are_equal); |
| 143 | } |
| 144 | |
| 145 | #[test] |
| 146 | fn transpose(matrix in arbitrary_bitmatrix(100)) { |
| 147 | let transposed = matrix.transposed(); |
| 148 | for row in 0..matrix.rowcount() { |
| 149 | for column in 0..matrix.columncount() { |
| 150 | assert_eq!(matrix[(row, column)], transposed[(column, row)]); |
| 151 | } |
| 152 | } |
| 153 | } |
| 154 | |
| 155 | #[test] |
| 156 | fn inverse(matrix in invertible_bitmatrix(100)) { |
| 157 | let inverted = matrix.inverted(); |
| 158 | let identity = BitMatrix::identity(matrix.rowcount()); |
| 159 | assert_eq!(&matrix * &inverted, identity); |
| 160 | } |
| 161 | |
| 162 | #[test] |
| 163 | fn echelon_form(matrix in arbitrary_bitmatrix(100)) { |
| 164 | let mut echeloned = matrix.clone(); |
| 165 | let profile = echeloned.echelonize(); |
| 166 | assert!(is_rref(&echeloned, &profile)); |
| 167 | assert!(preserves_rowspan_of(&matrix, &echeloned)); |
| 168 | } |
| 169 | |
| 170 | #[test] |
| 171 | fn tiny_matrix_echelon_form(matrix in fixed_size_bitmatrix(32,60)) { |
| 172 | let mut echeloned = matrix.clone(); |
| 173 | let _ = echeloned.echelonize(); |
| 174 | let mut tiny1 = tiny_matrix_from_bitmatrix::<32>(&matrix); |
| 175 | tiny_matrix_rref::<32,60>(&mut tiny1); |
| 176 | let tiny2 = tiny_matrix_from_bitmatrix::<32>(&echeloned); |
| 177 | assert_eq!(tiny1,tiny2); |
| 178 | } |
| 179 | |
| 180 | #[test] |
| 181 | fn direct_sum(left in arbitrary_bitmatrix(100), right in arbitrary_bitmatrix(100)) { |
| 182 | let summed = directly_summed([&left, &right]); |
| 183 | let expected_shape = (left.rowcount() + right.rowcount(), left.columncount() + right.columncount()); |
| 184 | assert_eq!(expected_shape, summed.shape()); |
| 185 | for row_index in 0..left.rowcount() { |
| 186 | for column_index in 0..left.columncount() { |
| 187 | assert_eq!(left[(row_index, column_index)], summed[(row_index, column_index)]); |
| 188 | } |
| 189 | for column_index in left.columncount()..summed.columncount() { |
| 190 | assert!(!summed[(row_index, column_index)]); |
| 191 | } |
| 192 | } |
| 193 | for row_index in 0..right.rowcount() { |
| 194 | for column_index in 0..right.columncount() { |
| 195 | assert_eq!(right[(row_index, column_index)], summed[(left.rowcount() + row_index, left.columncount() + column_index)]); |
| 196 | } |
| 197 | for column_index in 0..left.columncount() { |
| 198 | assert!(!summed[(left.rowcount() + row_index, column_index)]); |
| 199 | } |
| 200 | } |
| 201 | } |
| 202 | |
| 203 | } |
| 204 | |
| 205 | macro_rules! bitmatrix{ |
| 206 | ($($t:tt)+) => { |
| 207 | $crate::BitMatrix::<{paulimer::bits::WORD_COUNT_DEFAULT}>::from_str(stringify!($($t)+)).unwrap() |
| 208 | }; |
| 209 | } |
| 210 | |
| 211 | prop_compose! { |
| 212 | fn arbitrary_bitmatrix(max_dimension: usize)(shape in (0..=max_dimension, 0..=max_dimension)) -> BitMatrix { |
| 213 | random_bitmatrix(shape.0, shape.1) |
| 214 | } |
| 215 | } |
| 216 | |
| 217 | prop_compose! { |
| 218 | fn fixed_size_bitmatrix(row_count: usize, column_count: usize)(_ in 0..column_count) -> BitMatrix { |
| 219 | random_bitmatrix(row_count, column_count) |
| 220 | } |
| 221 | } |
| 222 | |
| 223 | prop_compose! { |
| 224 | fn invertible_bitmatrix(max_dimension: usize)(dimension in 1..=max_dimension) -> BitMatrix { |
| 225 | let mut matrix = BitMatrix::identity(dimension); |
| 226 | for _ in 0..dimension^2 { |
| 227 | let from_index = rand::rng().random_range(0..dimension); |
| 228 | let to_index = rand::rng().random_range(0..dimension); |
| 229 | if from_index != to_index { |
| 230 | matrix.add_into_row(to_index, from_index); |
| 231 | } |
| 232 | } |
| 233 | for _ in 0..dimension.pow(2) { |
| 234 | let from_index = rand::rng().random_range(0..dimension); |
| 235 | let to_index = rand::rng().random_range(0..dimension); |
| 236 | matrix.swap_rows(from_index, to_index); |
| 237 | } |
| 238 | matrix |
| 239 | } |
| 240 | } |
| 241 | |
| 242 | prop_compose! { |
| 243 | fn nonempty_bitmatrix(max_dimension: usize)(shape in (1..=max_dimension, 1..=max_dimension)) -> BitMatrix { |
| 244 | random_bitmatrix(shape.0, shape.1) |
| 245 | } |
| 246 | } |
| 247 | |
| 248 | prop_compose! { |
| 249 | fn equal_shape_bitmatrices(max_dimension: usize)(shape in (1..=max_dimension, 1..=max_dimension)) -> (BitMatrix, BitMatrix) { |
| 250 | (random_bitmatrix(shape.0, shape.1), random_bitmatrix(shape.0, shape.1)) |
| 251 | } |
| 252 | } |
| 253 | |
| 254 | // #[test] |
| 255 | // fn reduce() { |
| 256 | // for _ in 0..100 { |
| 257 | // let array = random_bitmatrix(100, 100); |
| 258 | // let reduced = rref(array); |
| 259 | // assert!(is_rref(&reduced)); |
| 260 | // } |
| 261 | |
| 262 | // for _ in 0..100 { |
| 263 | // let array = random_bitmatrix(50, 100); |
| 264 | // let (reduced, profile) = rref_with_rank_profile(array); |
| 265 | // assert_eq!(profile.len(), reduced.rowcount()); |
| 266 | // assert!(is_rref(&reduced)); |
| 267 | // } |
| 268 | |
| 269 | // { |
| 270 | // let matrix = bitmatrix!( |
| 271 | // |10 011 01| |
| 272 | // |.. 111 01| |
| 273 | // |.. ... 10|); |
| 274 | // assert!(is_rref(&matrix)); |
| 275 | // let (reduced, profile) = rref_with_rank_profile(matrix); |
| 276 | // assert!(is_rref(&reduced)); |
| 277 | // assert_eq!(profile, vec![0, 2, 5]); |
| 278 | // } |
| 279 | // } |
| 280 | |
| 281 | #[test] |
| 282 | fn reduce_with_transforms() { |
| 283 | for _ in 0..100 { |
| 284 | check_rref_with_transforms_on_random_matrix(100, 100); |
| 285 | } |
| 286 | for _ in 0..100 { |
| 287 | check_rref_with_transforms_on_random_matrix(50, 100); |
| 288 | } |
| 289 | } |
| 290 | |
| 291 | fn check_rref_with_transforms_on_random_matrix(nrows: usize, ncols: usize) { |
| 292 | let array = random_bitmatrix(nrows, ncols); |
| 293 | let (reduced, t, t_inv_t, profile) = rref_with_transforms(array.clone()); |
| 294 | assert!(is_rref(&reduced, &profile)); |
| 295 | assert_eq!(t.dot(&array), reduced); |
| 296 | assert_eq!( |
| 297 | t.dot(&t_inv_t.transposed()), |
| 298 | BitMatrix::identity(array.rowcount()) |
| 299 | ); |
| 300 | } |
| 301 | |
| 302 | #[test] |
| 303 | fn test_dot() { |
| 304 | println!("0"); |
| 305 | let x = bitmatrix!( |
| 306 | |01| |
| 307 | |10|); |
| 308 | let id = bitmatrix!( |
| 309 | |10| |
| 310 | |01|); |
| 311 | println!("1"); |
| 312 | assert_eq!(x.dot(&x), id); |
| 313 | assert_eq!(x.dot(&id), x); |
| 314 | assert_eq!(id.dot(&x), x); |
| 315 | |
| 316 | // multiplication is associative |
| 317 | println!("2"); |
| 318 | for _ in 0..100 { |
| 319 | let a = random_bitmatrix(10, 10); |
| 320 | let b = random_bitmatrix(10, 10); |
| 321 | let c = random_bitmatrix(10, 10); |
| 322 | assert_eq!((a.dot(&b)).dot(&c), a.dot(&b.dot(&c))); |
| 323 | } |
| 324 | |
| 325 | println!("3"); |
| 326 | // multiplication by zero is zero |
| 327 | for _ in 0..100 { |
| 328 | let a = random_bitmatrix(10, 10); |
| 329 | let z = BitMatrix::zeros(10, 10); |
| 330 | assert_eq!(a.dot(&z), z); |
| 331 | } |
| 332 | |
| 333 | // multiplication by id |
| 334 | for _ in 0..100 { |
| 335 | let a = random_bitmatrix(3, 3); |
| 336 | let id = BitMatrix::identity(3); |
| 337 | assert_eq!(a.dot(&id), a); |
| 338 | } |
| 339 | } |
| 340 | |
| 341 | #[test] |
| 342 | fn test_kernel_basis() { |
| 343 | let num_cols = 100; |
| 344 | for _ in 0..100 { |
| 345 | let mut matrix = random_bitmatrix(50, 100); |
| 346 | let rrp = matrix.echelonize(); |
| 347 | let mut kernel_basis_matrix = kernel_basis_matrix(&matrix); |
| 348 | let prod = matrix.dot(&kernel_basis_matrix.transposed()); |
| 349 | assert!(prod.is_zero()); |
| 350 | let rrpc = kernel_basis_matrix.echelonize(); |
| 351 | assert_eq!(rrp.len() + rrpc.len(), num_cols); |
| 352 | } |
| 353 | } |
| 354 | |
| 355 | fn preserves_rowspan_of(matrix: &BitMatrix, rref_matrix: &BitMatrix) -> bool { |
| 356 | let profile = fast_profile_of(rref_matrix); |
| 357 | let mut profile_rows = BTreeMap::new(); |
| 358 | for (row_index, column_index) in profile.iter().enumerate() { |
| 359 | profile_rows.insert(column_index, row_index); |
| 360 | } |
| 361 | for row in matrix.rows() { |
| 362 | let mut reduced = BitVec::<WORD_COUNT_DEFAULT>::from_view(&row); |
| 363 | let support = row |
| 364 | .support() |
| 365 | .assume_sorted_by_item() |
| 366 | .intersection(profile.iter().copied().assume_sorted_by_item()); |
| 367 | |
| 368 | for column_index in support { |
| 369 | let row_index = profile_rows[&column_index]; |
| 370 | let rref_row = BitVec::<WORD_COUNT_DEFAULT>::from_view(&rref_matrix.row(row_index)); |
| 371 | reduced.bitxor_assign(&rref_row); |
| 372 | } |
| 373 | if reduced.weight() > 0 { |
| 374 | return false; |
| 375 | } |
| 376 | } |
| 377 | true |
| 378 | } |
| 379 | |
| 380 | fn is_rref(matrix: &BitMatrix, with_profile: &[usize]) -> bool { |
| 381 | let expected_profile = fast_profile_of(matrix); |
| 382 | (expected_profile == with_profile) && columns_are_pivots_of(matrix, with_profile) |
| 383 | } |
| 384 | |
| 385 | fn columns_are_pivots_of(matrix: &BitMatrix, column_indexes: &[usize]) -> bool { |
| 386 | for &column_index in column_indexes { |
| 387 | let column = matrix.column(column_index); |
| 388 | if column.weight() != 1 { |
| 389 | return false; |
| 390 | } |
| 391 | } |
| 392 | true |
| 393 | } |
| 394 | |
| 395 | fn fast_profile_of(matrix: &BitMatrix) -> Vec<usize> { |
| 396 | let mut profile = vec![]; |
| 397 | for row_index in 0..matrix.rowcount() { |
| 398 | let row = matrix.row(row_index); |
| 399 | let pivot = row.into_iter().position(|bit| bit); |
| 400 | if pivot.is_none() { |
| 401 | break; |
| 402 | } |
| 403 | profile.push(pivot.unwrap()); |
| 404 | } |
| 405 | profile |
| 406 | } |
| 407 | |
| 408 | fn random_bitmatrix(rowcount: usize, columncount: usize) -> BitMatrix { |
| 409 | let mut matrix = BitMatrix::with_shape(rowcount, columncount); |
| 410 | let mut bits = std::iter::from_fn(move || Some(rand::rng().random::<bool>())); |
| 411 | for row_index in 0..rowcount { |
| 412 | for column_index in 0..columncount { |
| 413 | matrix.set((row_index, column_index), bits.next().expect("boom")); |
| 414 | } |
| 415 | } |
| 416 | for _ in 0..rowcount { |
| 417 | let from_index = rand::rng().random_range(0..rowcount); |
| 418 | let to_index = rand::rng().random_range(0..rowcount); |
| 419 | matrix.swap_rows(from_index, to_index); |
| 420 | } |
| 421 | matrix |
| 422 | } |
| 423 | |