aboutsummaryrefslogtreecommitdiff
path: root/vendor/smawk/src
diff options
context:
space:
mode:
Diffstat (limited to 'vendor/smawk/src')
-rw-r--r--vendor/smawk/src/brute_force.rs150
-rw-r--r--vendor/smawk/src/lib.rs570
-rw-r--r--vendor/smawk/src/monge.rs121
-rw-r--r--vendor/smawk/src/recursive.rs191
4 files changed, 1032 insertions, 0 deletions
diff --git a/vendor/smawk/src/brute_force.rs b/vendor/smawk/src/brute_force.rs
new file mode 100644
index 0000000..1ec0ca3
--- /dev/null
+++ b/vendor/smawk/src/brute_force.rs
@@ -0,0 +1,150 @@
+//! Brute-force algorithm for finding column minima.
+//!
+//! The functions here are mostly meant to be used for testing
+//! correctness of the SMAWK implementation.
+//!
+//! **Note: this module is only available if you enable the `ndarray`
+//! Cargo feature.**
+
+use ndarray::{Array2, ArrayView1};
+
+/// Compute lane minimum by brute force.
+///
+/// This does a simple scan through the lane (row or column).
+#[inline]
+pub fn lane_minimum<T: Ord>(lane: ArrayView1<'_, T>) -> usize {
+ lane.iter()
+ .enumerate()
+ .min_by_key(|&(idx, elem)| (elem, idx))
+ .map(|(idx, _)| idx)
+ .expect("empty lane in matrix")
+}
+
+/// Compute row minima by brute force in O(*mn*) time.
+///
+/// This function implements a simple brute-force approach where each
+/// matrix row is scanned completely. This means that the function
+/// works on all matrices, not just Monge matrices.
+///
+/// # Examples
+///
+/// ```
+/// let matrix = ndarray::arr2(&[[4, 2, 4, 3],
+/// [5, 3, 5, 3],
+/// [5, 3, 3, 1]]);
+/// assert_eq!(smawk::brute_force::row_minima(&matrix),
+/// vec![1, 1, 3]);
+/// ```
+///
+/// # Panics
+///
+/// It is an error to call this on a matrix with zero columns.
+pub fn row_minima<T: Ord>(matrix: &Array2<T>) -> Vec<usize> {
+ matrix.rows().into_iter().map(lane_minimum).collect()
+}
+
+/// Compute column minima by brute force in O(*mn*) time.
+///
+/// This function implements a simple brute-force approach where each
+/// matrix column is scanned completely. This means that the function
+/// works on all matrices, not just Monge matrices.
+///
+/// # Examples
+///
+/// ```
+/// let matrix = ndarray::arr2(&[[4, 2, 4, 3],
+/// [5, 3, 5, 3],
+/// [5, 3, 3, 1]]);
+/// assert_eq!(smawk::brute_force::column_minima(&matrix),
+/// vec![0, 0, 2, 2]);
+/// ```
+///
+/// # Panics
+///
+/// It is an error to call this on a matrix with zero rows.
+pub fn column_minima<T: Ord>(matrix: &Array2<T>) -> Vec<usize> {
+ matrix.columns().into_iter().map(lane_minimum).collect()
+}
+
+#[cfg(test)]
+mod tests {
+ use super::*;
+ use ndarray::arr2;
+
+ #[test]
+ fn brute_force_1x1() {
+ let matrix = arr2(&[[2]]);
+ let minima = vec![0];
+ assert_eq!(row_minima(&matrix), minima);
+ assert_eq!(column_minima(&matrix.reversed_axes()), minima);
+ }
+
+ #[test]
+ fn brute_force_2x1() {
+ let matrix = arr2(&[
+ [3], //
+ [2],
+ ]);
+ let minima = vec![0, 0];
+ assert_eq!(row_minima(&matrix), minima);
+ assert_eq!(column_minima(&matrix.reversed_axes()), minima);
+ }
+
+ #[test]
+ fn brute_force_1x2() {
+ let matrix = arr2(&[[2, 1]]);
+ let minima = vec![1];
+ assert_eq!(row_minima(&matrix), minima);
+ assert_eq!(column_minima(&matrix.reversed_axes()), minima);
+ }
+
+ #[test]
+ fn brute_force_2x2() {
+ let matrix = arr2(&[
+ [3, 2], //
+ [2, 1],
+ ]);
+ let minima = vec![1, 1];
+ assert_eq!(row_minima(&matrix), minima);
+ assert_eq!(column_minima(&matrix.reversed_axes()), minima);
+ }
+
+ #[test]
+ fn brute_force_3x3() {
+ let matrix = arr2(&[
+ [3, 4, 4], //
+ [3, 4, 4],
+ [2, 3, 3],
+ ]);
+ let minima = vec![0, 0, 0];
+ assert_eq!(row_minima(&matrix), minima);
+ assert_eq!(column_minima(&matrix.reversed_axes()), minima);
+ }
+
+ #[test]
+ fn brute_force_4x4() {
+ let matrix = arr2(&[
+ [4, 5, 5, 5], //
+ [2, 3, 3, 3],
+ [2, 3, 3, 3],
+ [2, 2, 2, 2],
+ ]);
+ let minima = vec![0, 0, 0, 0];
+ assert_eq!(row_minima(&matrix), minima);
+ assert_eq!(column_minima(&matrix.reversed_axes()), minima);
+ }
+
+ #[test]
+ fn brute_force_5x5() {
+ let matrix = arr2(&[
+ [3, 2, 4, 5, 6],
+ [2, 1, 3, 3, 4],
+ [2, 1, 3, 3, 4],
+ [3, 2, 4, 3, 4],
+ [4, 3, 2, 1, 1],
+ ]);
+ let minima = vec![1, 1, 1, 1, 3];
+ assert_eq!(row_minima(&matrix), minima);
+ assert_eq!(column_minima(&matrix.reversed_axes()), minima);
+ }
+}
diff --git a/vendor/smawk/src/lib.rs b/vendor/smawk/src/lib.rs
new file mode 100644
index 0000000..367d033
--- /dev/null
+++ b/vendor/smawk/src/lib.rs
@@ -0,0 +1,570 @@
+//! This crate implements various functions that help speed up dynamic
+//! programming, most importantly the SMAWK algorithm for finding row
+//! or column minima in a totally monotone matrix with *m* rows and
+//! *n* columns in time O(*m* + *n*). This is much better than the
+//! brute force solution which would take O(*mn*). When *m* and *n*
+//! are of the same order, this turns a quadratic function into a
+//! linear function.
+//!
+//! # Examples
+//!
+//! Computing the column minima of an *m* ✕ *n* Monge matrix can be
+//! done efficiently with `smawk::column_minima`:
+//!
+//! ```
+//! use smawk::Matrix;
+//!
+//! let matrix = vec![
+//! vec![3, 2, 4, 5, 6],
+//! vec![2, 1, 3, 3, 4],
+//! vec![2, 1, 3, 3, 4],
+//! vec![3, 2, 4, 3, 4],
+//! vec![4, 3, 2, 1, 1],
+//! ];
+//! let minima = vec![1, 1, 4, 4, 4];
+//! assert_eq!(smawk::column_minima(&matrix), minima);
+//! ```
+//!
+//! The `minima` vector gives the index of the minimum value per
+//! column, so `minima[0] == 1` since the minimum value in the first
+//! column is 2 (row 1). Note that the smallest row index is returned.
+//!
+//! # Definitions
+//!
+//! Some of the functions in this crate only work on matrices that are
+//! *totally monotone*, which we will define below.
+//!
+//! ## Monotone Matrices
+//!
+//! We start with a helper definition. Given an *m* ✕ *n* matrix `M`,
+//! we say that `M` is *monotone* when the minimum value of row `i` is
+//! found to the left of the minimum value in row `i'` where `i < i'`.
+//!
+//! More formally, if we let `rm(i)` denote the column index of the
+//! left-most minimum value in row `i`, then we have
+//!
+//! ```text
+//! rm(0) ≤ rm(1) ≤ ... ≤ rm(m - 1)
+//! ```
+//!
+//! This means that as you go down the rows from top to bottom, the
+//! row-minima proceed from left to right.
+//!
+//! The algorithms in this crate deal with finding such row- and
+//! column-minima.
+//!
+//! ## Totally Monotone Matrices
+//!
+//! We say that a matrix `M` is *totally monotone* when every
+//! sub-matrix is monotone. A sub-matrix is formed by the intersection
+//! of any two rows `i < i'` and any two columns `j < j'`.
+//!
+//! This is often expressed as via this equivalent condition:
+//!
+//! ```text
+//! M[i, j] > M[i, j'] => M[i', j] > M[i', j']
+//! ```
+//!
+//! for all `i < i'` and `j < j'`.
+//!
+//! ## Monge Property for Matrices
+//!
+//! A matrix `M` is said to fulfill the *Monge property* if
+//!
+//! ```text
+//! M[i, j] + M[i', j'] ≤ M[i, j'] + M[i', j]
+//! ```
+//!
+//! for all `i < i'` and `j < j'`. This says that given any rectangle
+//! in the matrix, the sum of the top-left and bottom-right corners is
+//! less than or equal to the sum of the bottom-left and upper-right
+//! corners.
+//!
+//! All Monge matrices are totally monotone, so it is enough to
+//! establish that the Monge property holds in order to use a matrix
+//! with the functions in this crate. If your program is dealing with
+//! unknown inputs, it can use [`monge::is_monge`] to verify that a
+//! matrix is a Monge matrix.
+
+#![doc(html_root_url = "https://docs.rs/smawk/0.3.2")]
+// The s! macro from ndarray uses unsafe internally, so we can only
+// forbid unsafe code when building with the default features.
+#![cfg_attr(not(feature = "ndarray"), forbid(unsafe_code))]
+
+#[cfg(feature = "ndarray")]
+pub mod brute_force;
+pub mod monge;
+#[cfg(feature = "ndarray")]
+pub mod recursive;
+
+/// Minimal matrix trait for two-dimensional arrays.
+///
+/// This provides the functionality needed to represent a read-only
+/// numeric matrix. You can query the size of the matrix and access
+/// elements. Modeled after [`ndarray::Array2`] from the [ndarray
+/// crate](https://crates.io/crates/ndarray).
+///
+/// Enable the `ndarray` Cargo feature if you want to use it with
+/// `ndarray::Array2`.
+pub trait Matrix<T: Copy> {
+ /// Return the number of rows.
+ fn nrows(&self) -> usize;
+ /// Return the number of columns.
+ fn ncols(&self) -> usize;
+ /// Return a matrix element.
+ fn index(&self, row: usize, column: usize) -> T;
+}
+
+/// Simple and inefficient matrix representation used for doctest
+/// examples and simple unit tests.
+///
+/// You should prefer implementing it yourself, or you can enable the
+/// `ndarray` Cargo feature and use the provided implementation for
+/// [`ndarray::Array2`].
+impl<T: Copy> Matrix<T> for Vec<Vec<T>> {
+ fn nrows(&self) -> usize {
+ self.len()
+ }
+ fn ncols(&self) -> usize {
+ self[0].len()
+ }
+ fn index(&self, row: usize, column: usize) -> T {
+ self[row][column]
+ }
+}
+
+/// Adapting [`ndarray::Array2`] to the `Matrix` trait.
+///
+/// **Note: this implementation is only available if you enable the
+/// `ndarray` Cargo feature.**
+#[cfg(feature = "ndarray")]
+impl<T: Copy> Matrix<T> for ndarray::Array2<T> {
+ #[inline]
+ fn nrows(&self) -> usize {
+ self.nrows()
+ }
+ #[inline]
+ fn ncols(&self) -> usize {
+ self.ncols()
+ }
+ #[inline]
+ fn index(&self, row: usize, column: usize) -> T {
+ self[[row, column]]
+ }
+}
+
+/// Compute row minima in O(*m* + *n*) time.
+///
+/// This implements the [SMAWK algorithm] for efficiently finding row
+/// minima in a totally monotone matrix.
+///
+/// The SMAWK algorithm is from Agarwal, Klawe, Moran, Shor, and
+/// Wilbur, *Geometric applications of a matrix searching algorithm*,
+/// Algorithmica 2, pp. 195-208 (1987) and the code here is a
+/// translation [David Eppstein's Python code][pads].
+///
+/// Running time on an *m* ✕ *n* matrix: O(*m* + *n*).
+///
+/// # Examples
+///
+/// ```
+/// use smawk::Matrix;
+/// let matrix = vec![vec![4, 2, 4, 3],
+/// vec![5, 3, 5, 3],
+/// vec![5, 3, 3, 1]];
+/// assert_eq!(smawk::row_minima(&matrix),
+/// vec![1, 1, 3]);
+/// ```
+///
+/// # Panics
+///
+/// It is an error to call this on a matrix with zero columns.
+///
+/// [pads]: https://github.com/jfinkels/PADS/blob/master/pads/smawk.py
+/// [SMAWK algorithm]: https://en.wikipedia.org/wiki/SMAWK_algorithm
+pub fn row_minima<T: PartialOrd + Copy, M: Matrix<T>>(matrix: &M) -> Vec<usize> {
+ // Benchmarking shows that SMAWK performs roughly the same on row-
+ // and column-major matrices.
+ let mut minima = vec![0; matrix.nrows()];
+ smawk_inner(
+ &|j, i| matrix.index(i, j),
+ &(0..matrix.ncols()).collect::<Vec<_>>(),
+ &(0..matrix.nrows()).collect::<Vec<_>>(),
+ &mut minima,
+ );
+ minima
+}
+
+#[deprecated(since = "0.3.2", note = "Please use `row_minima` instead.")]
+pub fn smawk_row_minima<T: PartialOrd + Copy, M: Matrix<T>>(matrix: &M) -> Vec<usize> {
+ row_minima(matrix)
+}
+
+/// Compute column minima in O(*m* + *n*) time.
+///
+/// This implements the [SMAWK algorithm] for efficiently finding
+/// column minima in a totally monotone matrix.
+///
+/// The SMAWK algorithm is from Agarwal, Klawe, Moran, Shor, and
+/// Wilbur, *Geometric applications of a matrix searching algorithm*,
+/// Algorithmica 2, pp. 195-208 (1987) and the code here is a
+/// translation [David Eppstein's Python code][pads].
+///
+/// Running time on an *m* ✕ *n* matrix: O(*m* + *n*).
+///
+/// # Examples
+///
+/// ```
+/// use smawk::Matrix;
+/// let matrix = vec![vec![4, 2, 4, 3],
+/// vec![5, 3, 5, 3],
+/// vec![5, 3, 3, 1]];
+/// assert_eq!(smawk::column_minima(&matrix),
+/// vec![0, 0, 2, 2]);
+/// ```
+///
+/// # Panics
+///
+/// It is an error to call this on a matrix with zero rows.
+///
+/// [SMAWK algorithm]: https://en.wikipedia.org/wiki/SMAWK_algorithm
+/// [pads]: https://github.com/jfinkels/PADS/blob/master/pads/smawk.py
+pub fn column_minima<T: PartialOrd + Copy, M: Matrix<T>>(matrix: &M) -> Vec<usize> {
+ let mut minima = vec![0; matrix.ncols()];
+ smawk_inner(
+ &|i, j| matrix.index(i, j),
+ &(0..matrix.nrows()).collect::<Vec<_>>(),
+ &(0..matrix.ncols()).collect::<Vec<_>>(),
+ &mut minima,
+ );
+ minima
+}
+
+#[deprecated(since = "0.3.2", note = "Please use `column_minima` instead.")]
+pub fn smawk_column_minima<T: PartialOrd + Copy, M: Matrix<T>>(matrix: &M) -> Vec<usize> {
+ column_minima(matrix)
+}
+
+/// Compute column minima in the given area of the matrix. The
+/// `minima` slice is updated inplace.
+fn smawk_inner<T: PartialOrd + Copy, M: Fn(usize, usize) -> T>(
+ matrix: &M,
+ rows: &[usize],
+ cols: &[usize],
+ minima: &mut [usize],
+) {
+ if cols.is_empty() {
+ return;
+ }
+
+ let mut stack = Vec::with_capacity(cols.len());
+ for r in rows {
+ // TODO: use stack.last() instead of stack.is_empty() etc
+ while !stack.is_empty()
+ && matrix(stack[stack.len() - 1], cols[stack.len() - 1])
+ > matrix(*r, cols[stack.len() - 1])
+ {
+ stack.pop();
+ }
+ if stack.len() != cols.len() {
+ stack.push(*r);
+ }
+ }
+ let rows = &stack;
+
+ let mut odd_cols = Vec::with_capacity(1 + cols.len() / 2);
+ for (idx, c) in cols.iter().enumerate() {
+ if idx % 2 == 1 {
+ odd_cols.push(*c);
+ }
+ }
+
+ smawk_inner(matrix, rows, &odd_cols, minima);
+
+ let mut r = 0;
+ for (c, &col) in cols.iter().enumerate().filter(|(c, _)| c % 2 == 0) {
+ let mut row = rows[r];
+ let last_row = if c == cols.len() - 1 {
+ rows[rows.len() - 1]
+ } else {
+ minima[cols[c + 1]]
+ };
+ let mut pair = (matrix(row, col), row);
+ while row != last_row {
+ r += 1;
+ row = rows[r];
+ if (matrix(row, col), row) < pair {
+ pair = (matrix(row, col), row);
+ }
+ }
+ minima[col] = pair.1;
+ }
+}
+
+/// Compute upper-right column minima in O(*m* + *n*) time.
+///
+/// The input matrix must be totally monotone.
+///
+/// The function returns a vector of `(usize, T)`. The `usize` in the
+/// tuple at index `j` tells you the row of the minimum value in
+/// column `j` and the `T` value is minimum value itself.
+///
+/// The algorithm only considers values above the main diagonal, which
+/// means that it computes values `v(j)` where:
+///
+/// ```text
+/// v(0) = initial
+/// v(j) = min { M[i, j] | i < j } for j > 0
+/// ```
+///
+/// If we let `r(j)` denote the row index of the minimum value in
+/// column `j`, the tuples in the result vector become `(r(j), M[r(j),
+/// j])`.
+///
+/// The algorithm is an *online* algorithm, in the sense that `matrix`
+/// function can refer back to previously computed column minima when
+/// determining an entry in the matrix. The guarantee is that we only
+/// call `matrix(i, j)` after having computed `v(i)`. This is
+/// reflected in the `&[(usize, T)]` argument to `matrix`, which grows
+/// as more and more values are computed.
+pub fn online_column_minima<T: Copy + PartialOrd, M: Fn(&[(usize, T)], usize, usize) -> T>(
+ initial: T,
+ size: usize,
+ matrix: M,
+) -> Vec<(usize, T)> {
+ let mut result = vec![(0, initial)];
+
+ // State used by the algorithm.
+ let mut finished = 0;
+ let mut base = 0;
+ let mut tentative = 0;
+
+ // Shorthand for evaluating the matrix. We need a macro here since
+ // we don't want to borrow the result vector.
+ macro_rules! m {
+ ($i:expr, $j:expr) => {{
+ assert!($i < $j, "(i, j) not above diagonal: ({}, {})", $i, $j);
+ assert!(
+ $i < size && $j < size,
+ "(i, j) out of bounds: ({}, {}), size: {}",
+ $i,
+ $j,
+ size
+ );
+ matrix(&result[..finished + 1], $i, $j)
+ }};
+ }
+
+ // Keep going until we have finished all size columns. Since the
+ // columns are zero-indexed, we're done when finished == size - 1.
+ while finished < size - 1 {
+ // First case: we have already advanced past the previous
+ // tentative value. We make a new tentative value by applying
+ // smawk_inner to the largest square submatrix that fits under
+ // the base.
+ let i = finished + 1;
+ if i > tentative {
+ let rows = (base..finished + 1).collect::<Vec<_>>();
+ tentative = std::cmp::min(finished + rows.len(), size - 1);
+ let cols = (finished + 1..tentative + 1).collect::<Vec<_>>();
+ let mut minima = vec![0; tentative + 1];
+ smawk_inner(&|i, j| m![i, j], &rows, &cols, &mut minima);
+ for col in cols {
+ let row = minima[col];
+ let v = m![row, col];
+ if col >= result.len() {
+ result.push((row, v));
+ } else if v < result[col].1 {
+ result[col] = (row, v);
+ }
+ }
+ finished = i;
+ continue;
+ }
+
+ // Second case: the new column minimum is on the diagonal. All
+ // subsequent ones will be at least as low, so we can clear
+ // out all our work from higher rows. As in the fourth case,
+ // the loss of tentative is amortized against the increase in
+ // base.
+ let diag = m![i - 1, i];
+ if diag < result[i].1 {
+ result[i] = (i - 1, diag);
+ base = i - 1;
+ tentative = i;
+ finished = i;
+ continue;
+ }
+
+ // Third case: row i-1 does not supply a column minimum in any
+ // column up to tentative. We simply advance finished while
+ // maintaining the invariant.
+ if m![i - 1, tentative] >= result[tentative].1 {
+ finished = i;
+ continue;
+ }
+
+ // Fourth and final case: a new column minimum at tentative.
+ // This allows us to make progress by incorporating rows prior
+ // to finished into the base. The base invariant holds because
+ // these rows cannot supply any later column minima. The work
+ // done when we last advanced tentative (and undone by this
+ // step) can be amortized against the increase in base.
+ base = i - 1;
+ tentative = i;
+ finished = i;
+ }
+
+ result
+}
+
+#[cfg(test)]
+mod tests {
+ use super::*;
+
+ #[test]
+ fn smawk_1x1() {
+ let matrix = vec![vec![2]];
+ assert_eq!(row_minima(&matrix), vec![0]);
+ assert_eq!(column_minima(&matrix), vec![0]);
+ }
+
+ #[test]
+ fn smawk_2x1() {
+ let matrix = vec![
+ vec![3], //
+ vec![2],
+ ];
+ assert_eq!(row_minima(&matrix), vec![0, 0]);
+ assert_eq!(column_minima(&matrix), vec![1]);
+ }
+
+ #[test]
+ fn smawk_1x2() {
+ let matrix = vec![vec![2, 1]];
+ assert_eq!(row_minima(&matrix), vec![1]);
+ assert_eq!(column_minima(&matrix), vec![0, 0]);
+ }
+
+ #[test]
+ fn smawk_2x2() {
+ let matrix = vec![
+ vec![3, 2], //
+ vec![2, 1],
+ ];
+ assert_eq!(row_minima(&matrix), vec![1, 1]);
+ assert_eq!(column_minima(&matrix), vec![1, 1]);
+ }
+
+ #[test]
+ fn smawk_3x3() {
+ let matrix = vec![
+ vec![3, 4, 4], //
+ vec![3, 4, 4],
+ vec![2, 3, 3],
+ ];
+ assert_eq!(row_minima(&matrix), vec![0, 0, 0]);
+ assert_eq!(column_minima(&matrix), vec![2, 2, 2]);
+ }
+
+ #[test]
+ fn smawk_4x4() {
+ let matrix = vec![
+ vec![4, 5, 5, 5], //
+ vec![2, 3, 3, 3],
+ vec![2, 3, 3, 3],
+ vec![2, 2, 2, 2],
+ ];
+ assert_eq!(row_minima(&matrix), vec![0, 0, 0, 0]);
+ assert_eq!(column_minima(&matrix), vec![1, 3, 3, 3]);
+ }
+
+ #[test]
+ fn smawk_5x5() {
+ let matrix = vec![
+ vec![3, 2, 4, 5, 6],
+ vec![2, 1, 3, 3, 4],
+ vec![2, 1, 3, 3, 4],
+ vec![3, 2, 4, 3, 4],
+ vec![4, 3, 2, 1, 1],
+ ];
+ assert_eq!(row_minima(&matrix), vec![1, 1, 1, 1, 3]);
+ assert_eq!(column_minima(&matrix), vec![1, 1, 4, 4, 4]);
+ }
+
+ #[test]
+ fn online_1x1() {
+ let matrix = vec![vec![0]];
+ let minima = vec![(0, 0)];
+ assert_eq!(online_column_minima(0, 1, |_, i, j| matrix[i][j]), minima);
+ }
+
+ #[test]
+ fn online_2x2() {
+ let matrix = vec![
+ vec![0, 2], //
+ vec![0, 0],
+ ];
+ let minima = vec![(0, 0), (0, 2)];
+ assert_eq!(online_column_minima(0, 2, |_, i, j| matrix[i][j]), minima);
+ }
+
+ #[test]
+ fn online_3x3() {
+ let matrix = vec![
+ vec![0, 4, 4], //
+ vec![0, 0, 4],
+ vec![0, 0, 0],
+ ];
+ let minima = vec![(0, 0), (0, 4), (0, 4)];
+ assert_eq!(online_column_minima(0, 3, |_, i, j| matrix[i][j]), minima);
+ }
+
+ #[test]
+ fn online_4x4() {
+ let matrix = vec![
+ vec![0, 5, 5, 5], //
+ vec![0, 0, 3, 3],
+ vec![0, 0, 0, 3],
+ vec![0, 0, 0, 0],
+ ];
+ let minima = vec![(0, 0), (0, 5), (1, 3), (1, 3)];
+ assert_eq!(online_column_minima(0, 4, |_, i, j| matrix[i][j]), minima);
+ }
+
+ #[test]
+ fn online_5x5() {
+ let matrix = vec![
+ vec![0, 2, 4, 6, 7],
+ vec![0, 0, 3, 4, 5],
+ vec![0, 0, 0, 3, 4],
+ vec![0, 0, 0, 0, 4],
+ vec![0, 0, 0, 0, 0],
+ ];
+ let minima = vec![(0, 0), (0, 2), (1, 3), (2, 3), (2, 4)];
+ assert_eq!(online_column_minima(0, 5, |_, i, j| matrix[i][j]), minima);
+ }
+
+ #[test]
+ fn smawk_works_with_partial_ord() {
+ let matrix = vec![
+ vec![3.0, 2.0], //
+ vec![2.0, 1.0],
+ ];
+ assert_eq!(row_minima(&matrix), vec![1, 1]);
+ assert_eq!(column_minima(&matrix), vec![1, 1]);
+ }
+
+ #[test]
+ fn online_works_with_partial_ord() {
+ let matrix = vec![
+ vec![0.0, 2.0], //
+ vec![0.0, 0.0],
+ ];
+ let minima = vec![(0, 0.0), (0, 2.0)];
+ assert_eq!(
+ online_column_minima(0.0, 2, |_, i: usize, j: usize| matrix[i][j]),
+ minima
+ );
+ }
+}
diff --git a/vendor/smawk/src/monge.rs b/vendor/smawk/src/monge.rs
new file mode 100644
index 0000000..dbc80e1
--- /dev/null
+++ b/vendor/smawk/src/monge.rs
@@ -0,0 +1,121 @@
+//! Functions for generating and checking Monge arrays.
+//!
+//! The functions here are mostly meant to be used for testing
+//! correctness of the SMAWK implementation.
+
+use crate::Matrix;
+use std::num::Wrapping;
+use std::ops::Add;
+
+/// Verify that a matrix is a Monge matrix.
+///
+/// A [Monge matrix] \(or array) is a matrix where the following
+/// inequality holds:
+///
+/// ```text
+/// M[i, j] + M[i', j'] <= M[i, j'] + M[i', j] for all i < i', j < j'
+/// ```
+///
+/// The inequality says that the sum of the main diagonal is less than
+/// the sum of the antidiagonal. Checking this condition is done by
+/// checking *n* ✕ *m* submatrices, so the running time is O(*mn*).
+///
+/// [Monge matrix]: https://en.wikipedia.org/wiki/Monge_array
+pub fn is_monge<T: Ord + Copy, M: Matrix<T>>(matrix: &M) -> bool
+where
+ Wrapping<T>: Add<Output = Wrapping<T>>,
+{
+ /// Returns `Ok(a + b)` if the computation can be done without
+ /// overflow, otherwise `Err(a + b - T::MAX - 1)` is returned.
+ fn checked_add<T: Ord + Copy>(a: Wrapping<T>, b: Wrapping<T>) -> Result<T, T>
+ where
+ Wrapping<T>: Add<Output = Wrapping<T>>,
+ {
+ let sum = a + b;
+ if sum < a {
+ Err(sum.0)
+ } else {
+ Ok(sum.0)
+ }
+ }
+
+ (0..matrix.nrows() - 1)
+ .flat_map(|row| (0..matrix.ncols() - 1).map(move |col| (row, col)))
+ .all(|(row, col)| {
+ let top_left = Wrapping(matrix.index(row, col));
+ let top_right = Wrapping(matrix.index(row, col + 1));
+ let bot_left = Wrapping(matrix.index(row + 1, col));
+ let bot_right = Wrapping(matrix.index(row + 1, col + 1));
+
+ match (
+ checked_add(top_left, bot_right),
+ checked_add(bot_left, top_right),
+ ) {
+ (Ok(a), Ok(b)) => a <= b, // No overflow.
+ (Err(a), Err(b)) => a <= b, // Double overflow.
+ (Ok(_), Err(_)) => true, // Anti-diagonal overflow.
+ (Err(_), Ok(_)) => false, // Main diagonal overflow.
+ }
+ })
+}
+
+#[cfg(test)]
+mod tests {
+ use super::*;
+
+ #[test]
+ fn is_monge_handles_overflow() {
+ // The x + y <= z + w computations will overflow for an u8
+ // matrix unless is_monge is careful.
+ let matrix: Vec<Vec<u8>> = vec![
+ vec![200, 200, 200, 200],
+ vec![200, 200, 200, 200],
+ vec![200, 200, 200, 200],
+ ];
+ assert!(is_monge(&matrix));
+ }
+
+ #[test]
+ fn monge_constant_rows() {
+ let matrix = vec![
+ vec![42, 42, 42, 42],
+ vec![0, 0, 0, 0],
+ vec![100, 100, 100, 100],
+ vec![1000, 1000, 1000, 1000],
+ ];
+ assert!(is_monge(&matrix));
+ }
+
+ #[test]
+ fn monge_constant_cols() {
+ let matrix = vec![
+ vec![42, 0, 100, 1000],
+ vec![42, 0, 100, 1000],
+ vec![42, 0, 100, 1000],
+ vec![42, 0, 100, 1000],
+ ];
+ assert!(is_monge(&matrix));
+ }
+
+ #[test]
+ fn monge_upper_right() {
+ let matrix = vec![
+ vec![10, 10, 42, 42, 42],
+ vec![10, 10, 42, 42, 42],
+ vec![10, 10, 10, 10, 10],
+ vec![10, 10, 10, 10, 10],
+ ];
+ assert!(is_monge(&matrix));
+ }
+
+ #[test]
+ fn monge_lower_left() {
+ let matrix = vec![
+ vec![10, 10, 10, 10, 10],
+ vec![10, 10, 10, 10, 10],
+ vec![42, 42, 42, 10, 10],
+ vec![42, 42, 42, 10, 10],
+ ];
+ assert!(is_monge(&matrix));
+ }
+}
diff --git a/vendor/smawk/src/recursive.rs b/vendor/smawk/src/recursive.rs
new file mode 100644
index 0000000..9df8b9c
--- /dev/null
+++ b/vendor/smawk/src/recursive.rs
@@ -0,0 +1,191 @@
+//! Recursive algorithm for finding column minima.
+//!
+//! The functions here are mostly meant to be used for testing
+//! correctness of the SMAWK implementation.
+//!
+//! **Note: this module is only available if you enable the `ndarray`
+//! Cargo feature.**
+
+use ndarray::{s, Array2, ArrayView2, Axis};
+
+/// Compute row minima in O(*m* + *n* log *m*) time.
+///
+/// This function computes row minima in a totally monotone matrix
+/// using a recursive algorithm.
+///
+/// Running time on an *m* ✕ *n* matrix: O(*m* + *n* log *m*).
+///
+/// # Examples
+///
+/// ```
+/// let matrix = ndarray::arr2(&[[4, 2, 4, 3],
+/// [5, 3, 5, 3],
+/// [5, 3, 3, 1]]);
+/// assert_eq!(smawk::recursive::row_minima(&matrix),
+/// vec![1, 1, 3]);
+/// ```
+///
+/// # Panics
+///
+/// It is an error to call this on a matrix with zero columns.
+pub fn row_minima<T: Ord>(matrix: &Array2<T>) -> Vec<usize> {
+ let mut minima = vec![0; matrix.nrows()];
+ recursive_inner(matrix.view(), &|| Direction::Row, 0, &mut minima);
+ minima
+}
+
+/// Compute column minima in O(*n* + *m* log *n*) time.
+///
+/// This function computes column minima in a totally monotone matrix
+/// using a recursive algorithm.
+///
+/// Running time on an *m* ✕ *n* matrix: O(*n* + *m* log *n*).
+///
+/// # Examples
+///
+/// ```
+/// let matrix = ndarray::arr2(&[[4, 2, 4, 3],
+/// [5, 3, 5, 3],
+/// [5, 3, 3, 1]]);
+/// assert_eq!(smawk::recursive::column_minima(&matrix),
+/// vec![0, 0, 2, 2]);
+/// ```
+///
+/// # Panics
+///
+/// It is an error to call this on a matrix with zero rows.
+pub fn column_minima<T: Ord>(matrix: &Array2<T>) -> Vec<usize> {
+ let mut minima = vec![0; matrix.ncols()];
+ recursive_inner(matrix.view(), &|| Direction::Column, 0, &mut minima);
+ minima
+}
+
+/// The type of minima (row or column) we compute.
+enum Direction {
+ Row,
+ Column,
+}
+
+/// Compute the minima along the given direction (`Direction::Row` for
+/// row minima and `Direction::Column` for column minima).
+///
+/// The direction is given as a generic function argument to allow
+/// monomorphization to kick in. The function calls will be inlined
+/// and optimized away and the result is that the compiler generates
+/// differnet code for finding row and column minima.
+fn recursive_inner<T: Ord, F: Fn() -> Direction>(
+ matrix: ArrayView2<'_, T>,
+ dir: &F,
+ offset: usize,
+ minima: &mut [usize],
+) {
+ if matrix.is_empty() {
+ return;
+ }
+
+ let axis = match dir() {
+ Direction::Row => Axis(0),
+ Direction::Column => Axis(1),
+ };
+ let mid = matrix.len_of(axis) / 2;
+ let min_idx = crate::brute_force::lane_minimum(matrix.index_axis(axis, mid));
+ minima[mid] = offset + min_idx;
+
+ if mid == 0 {
+ return; // Matrix has a single row or column, so we're done.
+ }
+
+ let top_left = match dir() {
+ Direction::Row => matrix.slice(s![..mid, ..(min_idx + 1)]),
+ Direction::Column => matrix.slice(s![..(min_idx + 1), ..mid]),
+ };
+ let bot_right = match dir() {
+ Direction::Row => matrix.slice(s![(mid + 1).., min_idx..]),
+ Direction::Column => matrix.slice(s![min_idx.., (mid + 1)..]),
+ };
+ recursive_inner(top_left, dir, offset, &mut minima[..mid]);
+ recursive_inner(bot_right, dir, offset + min_idx, &mut minima[mid + 1..]);
+}
+
+#[cfg(test)]
+mod tests {
+ use super::*;
+ use ndarray::arr2;
+
+ #[test]
+ fn recursive_1x1() {
+ let matrix = arr2(&[[2]]);
+ let minima = vec![0];
+ assert_eq!(row_minima(&matrix), minima);
+ assert_eq!(column_minima(&matrix.reversed_axes()), minima);
+ }
+
+ #[test]
+ fn recursive_2x1() {
+ let matrix = arr2(&[
+ [3], //
+ [2],
+ ]);
+ let minima = vec![0, 0];
+ assert_eq!(row_minima(&matrix), minima);
+ assert_eq!(column_minima(&matrix.reversed_axes()), minima);
+ }
+
+ #[test]
+ fn recursive_1x2() {
+ let matrix = arr2(&[[2, 1]]);
+ let minima = vec![1];
+ assert_eq!(row_minima(&matrix), minima);
+ assert_eq!(column_minima(&matrix.reversed_axes()), minima);
+ }
+
+ #[test]
+ fn recursive_2x2() {
+ let matrix = arr2(&[
+ [3, 2], //
+ [2, 1],
+ ]);
+ let minima = vec![1, 1];
+ assert_eq!(row_minima(&matrix), minima);
+ assert_eq!(column_minima(&matrix.reversed_axes()), minima);
+ }
+
+ #[test]
+ fn recursive_3x3() {
+ let matrix = arr2(&[
+ [3, 4, 4], //
+ [3, 4, 4],
+ [2, 3, 3],
+ ]);
+ let minima = vec![0, 0, 0];
+ assert_eq!(row_minima(&matrix), minima);
+ assert_eq!(column_minima(&matrix.reversed_axes()), minima);
+ }
+
+ #[test]
+ fn recursive_4x4() {
+ let matrix = arr2(&[
+ [4, 5, 5, 5], //
+ [2, 3, 3, 3],
+ [2, 3, 3, 3],
+ [2, 2, 2, 2],
+ ]);
+ let minima = vec![0, 0, 0, 0];
+ assert_eq!(row_minima(&matrix), minima);
+ assert_eq!(column_minima(&matrix.reversed_axes()), minima);
+ }
+
+ #[test]
+ fn recursive_5x5() {
+ let matrix = arr2(&[
+ [3, 2, 4, 5, 6],
+ [2, 1, 3, 3, 4],
+ [2, 1, 3, 3, 4],
+ [3, 2, 4, 3, 4],
+ [4, 3, 2, 1, 1],
+ ]);
+ let minima = vec![1, 1, 1, 1, 3];
+ assert_eq!(row_minima(&matrix), minima);
+ assert_eq!(column_minima(&matrix.reversed_axes()), minima);
+ }
+}
Copyright © 2024 Valentin A. Popov. All Rights Reverse Engineered.