1 //! This crate implements various functions that help speed up dynamic
2 //! programming, most importantly the SMAWK algorithm for finding row
3 //! or column minima in a totally monotone matrix with *m* rows and
4 //! *n* columns in time O(*m* + *n*). This is much better than the
5 //! brute force solution which would take O(*mn*). When *m* and *n*
6 //! are of the same order, this turns a quadratic function into a
7 //! linear function.
8 //!
9 //! # Examples
10 //!
11 //! Computing the column minima of an *m* ✕ *n* Monge matrix can be
12 //! done efficiently with `smawk_column_minima`:
13 //!
14 //! ```
15 //! use smawk::{Matrix, smawk_column_minima};
16 //!
17 //! let matrix = vec![
18 //! vec![3, 2, 4, 5, 6],
19 //! vec![2, 1, 3, 3, 4],
20 //! vec![2, 1, 3, 3, 4],
21 //! vec![3, 2, 4, 3, 4],
22 //! vec![4, 3, 2, 1, 1],
23 //! ];
24 //! let minima = vec![1, 1, 4, 4, 4];
25 //! assert_eq!(smawk_column_minima(&matrix), minima);
26 //! ```
27 //!
28 //! The `minima` vector gives the index of the minimum value per
29 //! column, so `minima[0] == 1` since the minimum value in the first
30 //! column is 2 (row 1). Note that the smallest row index is returned.
31 //!
32 //! # Definitions
33 //!
34 //! Some of the functions in this crate only work on matrices that are
35 //! *totally monotone*, which we will define below.
36 //!
37 //! ## Monotone Matrices
38 //!
39 //! We start with a helper definition. Given an *m* ✕ *n* matrix `M`,
40 //! we say that `M` is *monotone* when the minimum value of row `i` is
41 //! found to the left of the minimum value in row `i'` where `i < i'`.
42 //!
43 //! More formally, if we let `rm(i)` denote the column index of the
44 //! left-most minimum value in row `i`, then we have
45 //!
46 //! ```text
47 //! rm(0) ≤ rm(1) ≤ ... ≤ rm(m - 1)
48 //! ```
49 //!
50 //! This means that as you go down the rows from top to bottom, the
51 //! row-minima proceed from left to right.
52 //!
53 //! The algorithms in this crate deal with finding such row- and
54 //! column-minima.
55 //!
56 //! ## Totally Monotone Matrices
57 //!
58 //! We say that a matrix `M` is *totally monotone* when every
59 //! sub-matrix is monotone. A sub-matrix is formed by the intersection
60 //! of any two rows `i < i'` and any two columns `j < j'`.
61 //!
62 //! This is often expressed as via this equivalent condition:
63 //!
64 //! ```text
65 //! M[i, j] > M[i, j'] => M[i', j] > M[i', j']
66 //! ```
67 //!
68 //! for all `i < i'` and `j < j'`.
69 //!
70 //! ## Monge Property for Matrices
71 //!
72 //! A matrix `M` is said to fulfill the *Monge property* if
73 //!
74 //! ```text
75 //! M[i, j] + M[i', j'] ≤ M[i, j'] + M[i', j]
76 //! ```
77 //!
78 //! for all `i < i'` and `j < j'`. This says that given any rectangle
79 //! in the matrix, the sum of the top-left and bottom-right corners is
80 //! less than or equal to the sum of the bottom-left and upper-right
81 //! corners.
82 //!
83 //! All Monge matrices are totally monotone, so it is enough to
84 //! establish that the Monge property holds in order to use a matrix
85 //! with the functions in this crate. If your program is dealing with
86 //! unknown inputs, it can use [`monge::is_monge`] to verify that a
87 //! matrix is a Monge matrix.
88
89 #![doc(html_root_url = "https://docs.rs/smawk/0.3.1")]
90
91 #[cfg(feature = "ndarray")]
92 pub mod brute_force;
93 pub mod monge;
94 #[cfg(feature = "ndarray")]
95 pub mod recursive;
96
97 /// Minimal matrix trait for two-dimensional arrays.
98 ///
99 /// This provides the functionality needed to represent a read-only
100 /// numeric matrix. You can query the size of the matrix and access
101 /// elements. Modeled after [`ndarray::Array2`] from the [ndarray
102 /// crate](https://crates.io/crates/ndarray).
103 ///
104 /// Enable the `ndarray` Cargo feature if you want to use it with
105 /// `ndarray::Array2`.
106 pub trait Matrix<T: Copy> {
107 /// Return the number of rows.
nrows(&self) -> usize108 fn nrows(&self) -> usize;
109 /// Return the number of columns.
ncols(&self) -> usize110 fn ncols(&self) -> usize;
111 /// Return a matrix element.
index(&self, row: usize, column: usize) -> T112 fn index(&self, row: usize, column: usize) -> T;
113 }
114
115 /// Simple and inefficient matrix representation used for doctest
116 /// examples and simple unit tests.
117 ///
118 /// You should prefer implementing it yourself, or you can enable the
119 /// `ndarray` Cargo feature and use the provided implementation for
120 /// [`ndarray::Array2`].
121 impl<T: Copy> Matrix<T> for Vec<Vec<T>> {
nrows(&self) -> usize122 fn nrows(&self) -> usize {
123 self.len()
124 }
ncols(&self) -> usize125 fn ncols(&self) -> usize {
126 self[0].len()
127 }
index(&self, row: usize, column: usize) -> T128 fn index(&self, row: usize, column: usize) -> T {
129 self[row][column]
130 }
131 }
132
133 /// Adapting [`ndarray::Array2`] to the `Matrix` trait.
134 ///
135 /// **Note: this implementation is only available if you enable the
136 /// `ndarray` Cargo feature.**
137 #[cfg(feature = "ndarray")]
138 impl<T: Copy> Matrix<T> for ndarray::Array2<T> {
139 #[inline]
nrows(&self) -> usize140 fn nrows(&self) -> usize {
141 self.nrows()
142 }
143 #[inline]
ncols(&self) -> usize144 fn ncols(&self) -> usize {
145 self.ncols()
146 }
147 #[inline]
index(&self, row: usize, column: usize) -> T148 fn index(&self, row: usize, column: usize) -> T {
149 self[[row, column]]
150 }
151 }
152
153 /// Compute row minima in O(*m* + *n*) time.
154 ///
155 /// This implements the SMAWK algorithm for finding row minima in a
156 /// totally monotone matrix.
157 ///
158 /// The SMAWK algorithm is from Agarwal, Klawe, Moran, Shor, and
159 /// Wilbur, *Geometric applications of a matrix searching algorithm*,
160 /// Algorithmica 2, pp. 195-208 (1987) and the code here is a
161 /// translation [David Eppstein's Python code][pads].
162 ///
163 /// [pads]: https://github.com/jfinkels/PADS/blob/master/pads/smawk.py
164 ///
165 /// Running time on an *m* ✕ *n* matrix: O(*m* + *n*).
166 ///
167 /// # Panics
168 ///
169 /// It is an error to call this on a matrix with zero columns.
smawk_row_minima<T: PartialOrd + Copy, M: Matrix<T>>(matrix: &M) -> Vec<usize>170 pub fn smawk_row_minima<T: PartialOrd + Copy, M: Matrix<T>>(matrix: &M) -> Vec<usize> {
171 // Benchmarking shows that SMAWK performs roughly the same on row-
172 // and column-major matrices.
173 let mut minima = vec![0; matrix.nrows()];
174 smawk_inner(
175 &|j, i| matrix.index(i, j),
176 &(0..matrix.ncols()).collect::<Vec<_>>(),
177 &(0..matrix.nrows()).collect::<Vec<_>>(),
178 &mut minima,
179 );
180 minima
181 }
182
183 /// Compute column minima in O(*m* + *n*) time.
184 ///
185 /// This implements the SMAWK algorithm for finding column minima in a
186 /// totally monotone matrix.
187 ///
188 /// The SMAWK algorithm is from Agarwal, Klawe, Moran, Shor, and
189 /// Wilbur, *Geometric applications of a matrix searching algorithm*,
190 /// Algorithmica 2, pp. 195-208 (1987) and the code here is a
191 /// translation [David Eppstein's Python code][pads].
192 ///
193 /// [pads]: https://github.com/jfinkels/PADS/blob/master/pads/smawk.py
194 ///
195 /// Running time on an *m* ✕ *n* matrix: O(*m* + *n*).
196 ///
197 /// # Panics
198 ///
199 /// It is an error to call this on a matrix with zero rows.
smawk_column_minima<T: PartialOrd + Copy, M: Matrix<T>>(matrix: &M) -> Vec<usize>200 pub fn smawk_column_minima<T: PartialOrd + Copy, M: Matrix<T>>(matrix: &M) -> Vec<usize> {
201 let mut minima = vec![0; matrix.ncols()];
202 smawk_inner(
203 &|i, j| matrix.index(i, j),
204 &(0..matrix.nrows()).collect::<Vec<_>>(),
205 &(0..matrix.ncols()).collect::<Vec<_>>(),
206 &mut minima,
207 );
208 minima
209 }
210
211 /// Compute column minima in the given area of the matrix. The
212 /// `minima` slice is updated inplace.
smawk_inner<T: PartialOrd + Copy, M: Fn(usize, usize) -> T>( matrix: &M, rows: &[usize], cols: &[usize], mut minima: &mut [usize], )213 fn smawk_inner<T: PartialOrd + Copy, M: Fn(usize, usize) -> T>(
214 matrix: &M,
215 rows: &[usize],
216 cols: &[usize],
217 mut minima: &mut [usize],
218 ) {
219 if cols.is_empty() {
220 return;
221 }
222
223 let mut stack = Vec::with_capacity(cols.len());
224 for r in rows {
225 // TODO: use stack.last() instead of stack.is_empty() etc
226 while !stack.is_empty()
227 && matrix(stack[stack.len() - 1], cols[stack.len() - 1])
228 > matrix(*r, cols[stack.len() - 1])
229 {
230 stack.pop();
231 }
232 if stack.len() != cols.len() {
233 stack.push(*r);
234 }
235 }
236 let rows = &stack;
237
238 let mut odd_cols = Vec::with_capacity(1 + cols.len() / 2);
239 for (idx, c) in cols.iter().enumerate() {
240 if idx % 2 == 1 {
241 odd_cols.push(*c);
242 }
243 }
244
245 smawk_inner(matrix, rows, &odd_cols, &mut minima);
246
247 let mut r = 0;
248 for (c, &col) in cols.iter().enumerate().filter(|(c, _)| c % 2 == 0) {
249 let mut row = rows[r];
250 let last_row = if c == cols.len() - 1 {
251 rows[rows.len() - 1]
252 } else {
253 minima[cols[c + 1]]
254 };
255 let mut pair = (matrix(row, col), row);
256 while row != last_row {
257 r += 1;
258 row = rows[r];
259 if (matrix(row, col), row) < pair {
260 pair = (matrix(row, col), row);
261 }
262 }
263 minima[col] = pair.1;
264 }
265 }
266
267 /// Compute upper-right column minima in O(*m* + *n*) time.
268 ///
269 /// The input matrix must be totally monotone.
270 ///
271 /// The function returns a vector of `(usize, T)`. The `usize` in the
272 /// tuple at index `j` tells you the row of the minimum value in
273 /// column `j` and the `T` value is minimum value itself.
274 ///
275 /// The algorithm only considers values above the main diagonal, which
276 /// means that it computes values `v(j)` where:
277 ///
278 /// ```text
279 /// v(0) = initial
280 /// v(j) = min { M[i, j] | i < j } for j > 0
281 /// ```
282 ///
283 /// If we let `r(j)` denote the row index of the minimum value in
284 /// column `j`, the tuples in the result vector become `(r(j), M[r(j),
285 /// j])`.
286 ///
287 /// The algorithm is an *online* algorithm, in the sense that `matrix`
288 /// function can refer back to previously computed column minima when
289 /// determining an entry in the matrix. The guarantee is that we only
290 /// call `matrix(i, j)` after having computed `v(i)`. This is
291 /// reflected in the `&[(usize, T)]` argument to `matrix`, which grows
292 /// as more and more values are computed.
online_column_minima<T: Copy + PartialOrd, M: Fn(&[(usize, T)], usize, usize) -> T>( initial: T, size: usize, matrix: M, ) -> Vec<(usize, T)>293 pub fn online_column_minima<T: Copy + PartialOrd, M: Fn(&[(usize, T)], usize, usize) -> T>(
294 initial: T,
295 size: usize,
296 matrix: M,
297 ) -> Vec<(usize, T)> {
298 let mut result = vec![(0, initial)];
299
300 // State used by the algorithm.
301 let mut finished = 0;
302 let mut base = 0;
303 let mut tentative = 0;
304
305 // Shorthand for evaluating the matrix. We need a macro here since
306 // we don't want to borrow the result vector.
307 macro_rules! m {
308 ($i:expr, $j:expr) => {{
309 assert!($i < $j, "(i, j) not above diagonal: ({}, {})", $i, $j);
310 assert!(
311 $i < size && $j < size,
312 "(i, j) out of bounds: ({}, {}), size: {}",
313 $i,
314 $j,
315 size
316 );
317 matrix(&result[..finished + 1], $i, $j)
318 }};
319 }
320
321 // Keep going until we have finished all size columns. Since the
322 // columns are zero-indexed, we're done when finished == size - 1.
323 while finished < size - 1 {
324 // First case: we have already advanced past the previous
325 // tentative value. We make a new tentative value by applying
326 // smawk_inner to the largest square submatrix that fits under
327 // the base.
328 let i = finished + 1;
329 if i > tentative {
330 let rows = (base..finished + 1).collect::<Vec<_>>();
331 tentative = std::cmp::min(finished + rows.len(), size - 1);
332 let cols = (finished + 1..tentative + 1).collect::<Vec<_>>();
333 let mut minima = vec![0; tentative + 1];
334 smawk_inner(&|i, j| m![i, j], &rows, &cols, &mut minima);
335 for col in cols {
336 let row = minima[col];
337 let v = m![row, col];
338 if col >= result.len() {
339 result.push((row, v));
340 } else if v < result[col].1 {
341 result[col] = (row, v);
342 }
343 }
344 finished = i;
345 continue;
346 }
347
348 // Second case: the new column minimum is on the diagonal. All
349 // subsequent ones will be at least as low, so we can clear
350 // out all our work from higher rows. As in the fourth case,
351 // the loss of tentative is amortized against the increase in
352 // base.
353 let diag = m![i - 1, i];
354 if diag < result[i].1 {
355 result[i] = (i - 1, diag);
356 base = i - 1;
357 tentative = i;
358 finished = i;
359 continue;
360 }
361
362 // Third case: row i-1 does not supply a column minimum in any
363 // column up to tentative. We simply advance finished while
364 // maintaining the invariant.
365 if m![i - 1, tentative] >= result[tentative].1 {
366 finished = i;
367 continue;
368 }
369
370 // Fourth and final case: a new column minimum at tentative.
371 // This allows us to make progress by incorporating rows prior
372 // to finished into the base. The base invariant holds because
373 // these rows cannot supply any later column minima. The work
374 // done when we last advanced tentative (and undone by this
375 // step) can be amortized against the increase in base.
376 base = i - 1;
377 tentative = i;
378 finished = i;
379 }
380
381 result
382 }
383
384 #[cfg(test)]
385 mod tests {
386 use super::*;
387
388 #[test]
smawk_1x1()389 fn smawk_1x1() {
390 let matrix = vec![vec![2]];
391 assert_eq!(smawk_row_minima(&matrix), vec![0]);
392 assert_eq!(smawk_column_minima(&matrix), vec![0]);
393 }
394
395 #[test]
smawk_2x1()396 fn smawk_2x1() {
397 let matrix = vec![
398 vec![3], //
399 vec![2],
400 ];
401 assert_eq!(smawk_row_minima(&matrix), vec![0, 0]);
402 assert_eq!(smawk_column_minima(&matrix), vec![1]);
403 }
404
405 #[test]
smawk_1x2()406 fn smawk_1x2() {
407 let matrix = vec![vec![2, 1]];
408 assert_eq!(smawk_row_minima(&matrix), vec![1]);
409 assert_eq!(smawk_column_minima(&matrix), vec![0, 0]);
410 }
411
412 #[test]
smawk_2x2()413 fn smawk_2x2() {
414 let matrix = vec![
415 vec![3, 2], //
416 vec![2, 1],
417 ];
418 assert_eq!(smawk_row_minima(&matrix), vec![1, 1]);
419 assert_eq!(smawk_column_minima(&matrix), vec![1, 1]);
420 }
421
422 #[test]
smawk_3x3()423 fn smawk_3x3() {
424 let matrix = vec![
425 vec![3, 4, 4], //
426 vec![3, 4, 4],
427 vec![2, 3, 3],
428 ];
429 assert_eq!(smawk_row_minima(&matrix), vec![0, 0, 0]);
430 assert_eq!(smawk_column_minima(&matrix), vec![2, 2, 2]);
431 }
432
433 #[test]
smawk_4x4()434 fn smawk_4x4() {
435 let matrix = vec![
436 vec![4, 5, 5, 5], //
437 vec![2, 3, 3, 3],
438 vec![2, 3, 3, 3],
439 vec![2, 2, 2, 2],
440 ];
441 assert_eq!(smawk_row_minima(&matrix), vec![0, 0, 0, 0]);
442 assert_eq!(smawk_column_minima(&matrix), vec![1, 3, 3, 3]);
443 }
444
445 #[test]
smawk_5x5()446 fn smawk_5x5() {
447 let matrix = vec![
448 vec![3, 2, 4, 5, 6],
449 vec![2, 1, 3, 3, 4],
450 vec![2, 1, 3, 3, 4],
451 vec![3, 2, 4, 3, 4],
452 vec![4, 3, 2, 1, 1],
453 ];
454 assert_eq!(smawk_row_minima(&matrix), vec![1, 1, 1, 1, 3]);
455 assert_eq!(smawk_column_minima(&matrix), vec![1, 1, 4, 4, 4]);
456 }
457
458 #[test]
online_1x1()459 fn online_1x1() {
460 let matrix = vec![vec![0]];
461 let minima = vec![(0, 0)];
462 assert_eq!(online_column_minima(0, 1, |_, i, j| matrix[i][j]), minima);
463 }
464
465 #[test]
online_2x2()466 fn online_2x2() {
467 let matrix = vec![
468 vec![0, 2], //
469 vec![0, 0],
470 ];
471 let minima = vec![(0, 0), (0, 2)];
472 assert_eq!(online_column_minima(0, 2, |_, i, j| matrix[i][j]), minima);
473 }
474
475 #[test]
online_3x3()476 fn online_3x3() {
477 let matrix = vec![
478 vec![0, 4, 4], //
479 vec![0, 0, 4],
480 vec![0, 0, 0],
481 ];
482 let minima = vec![(0, 0), (0, 4), (0, 4)];
483 assert_eq!(online_column_minima(0, 3, |_, i, j| matrix[i][j]), minima);
484 }
485
486 #[test]
online_4x4()487 fn online_4x4() {
488 let matrix = vec![
489 vec![0, 5, 5, 5], //
490 vec![0, 0, 3, 3],
491 vec![0, 0, 0, 3],
492 vec![0, 0, 0, 0],
493 ];
494 let minima = vec![(0, 0), (0, 5), (1, 3), (1, 3)];
495 assert_eq!(online_column_minima(0, 4, |_, i, j| matrix[i][j]), minima);
496 }
497
498 #[test]
online_5x5()499 fn online_5x5() {
500 let matrix = vec![
501 vec![0, 2, 4, 6, 7],
502 vec![0, 0, 3, 4, 5],
503 vec![0, 0, 0, 3, 4],
504 vec![0, 0, 0, 0, 4],
505 vec![0, 0, 0, 0, 0],
506 ];
507 let minima = vec![(0, 0), (0, 2), (1, 3), (2, 3), (2, 4)];
508 assert_eq!(online_column_minima(0, 5, |_, i, j| matrix[i][j]), minima);
509 }
510
511 #[test]
smawk_works_with_partial_ord()512 fn smawk_works_with_partial_ord() {
513 let matrix = vec![
514 vec![3.0, 2.0], //
515 vec![2.0, 1.0],
516 ];
517 assert_eq!(smawk_row_minima(&matrix), vec![1, 1]);
518 assert_eq!(smawk_column_minima(&matrix), vec![1, 1]);
519 }
520
521 #[test]
online_works_with_partial_ord()522 fn online_works_with_partial_ord() {
523 let matrix = vec![
524 vec![0.0, 2.0], //
525 vec![0.0, 0.0],
526 ];
527 let minima = vec![(0, 0.0), (0, 2.0)];
528 assert_eq!(online_column_minima(0.0, 2, |_, i:usize, j:usize| matrix[i][j]), minima);
529 }
530
531 }
532