Min Plus Convolution
(library/convolution/min_plus_convolution.hpp)

• View this file on GitHub
• Last update: 2023-09-06 20:32:52+09:00
• Include: #include "library/convolution/min_plus_convolution.hpp"

Min Plus Convolution

Depends on

• Monotone Minima (library/algorithm/monotone_minima.hpp)

Verified with

• test/src/convolution/min_plus_convolution/min_plus_convolution_convex_arbitrary.test.cpp
• test/src/convolution/min_plus_convolution/min_plus_convolution_convex_convex.test.cpp

Code

#ifndef SUISEN_MIN_PLUS_CONVOLUTION
#define SUISEN_MIN_PLUS_CONVOLUTION

#include <cassert>
#include <cstddef>
#include <limits>
#include <vector>

#include "library/algorithm/monotone_minima.hpp"

namespace suisen {
    template <typename T>
    std::vector<T> min_plus_convolution_convex_convex(const std::vector<T> &a, const std::vector<T> &b) {
        const int n = a.size(), m = b.size();
        assert(n and m);
        // check if convex
        for (int i = 2; i < n; ++i) assert(a[i - 1] - a[i - 2] <= a[i] - a[i - 1]);
        // check if convex
        for (int j = 2; j < m; ++j) assert(b[j - 1] - b[j - 2] <= b[j] - b[j - 1]);

        std::vector<T> c(n + m - 1);
        c[0] = a[0] + b[0];
        for (int k = 0, i = 0; k < n + m - 2; ++k) {
            int j = k - i;
            if (j == m - 1 or (i < n - 1 and a[i + 1] + b[j] < a[i] + b[j + 1])) {
                c[k + 1] = a[++i] + b[j];
            } else {
                c[k + 1] = a[i] + b[++j];
            }
        }
        return c;
    }

    template <typename T>
    std::vector<T> min_plus_convolution_convex_arbitrary(const std::vector<T> &a, const std::vector<T> &b) {
        const int n = a.size(), m = b.size();
        assert(n and m);
        // check if convex
        for (int i = 2; i < n; ++i) assert(a[i - 1] - a[i - 2] <= a[i] - a[i - 1]);

        std::vector<int> argmin = monotone_minima(
            n + m - 1, m,
            [&](int k, int j1, int j2) {
                const int i1 = k - j1, i2 = k - j2;
                // upper right triangle
                if (i2 < 0) return true;
                // lower left triangle
                if (i1 >= n) return false;
                return a[i1] + b[j1] < a[i2] + b[j2];
            }
        );
        std::vector<T> c(n + m - 1);
        for (int k = 0; k < n + m - 1; ++k) {
            const int j = argmin[k], i = k - j;
            c[k] = a[i] + b[j];
        }
        return c;
    }
} // namespace suisen

#endif // SUISEN_MIN_PLUS_CONVOLUTION

#line 1 "library/convolution/min_plus_convolution.hpp"



#include <cassert>
#include <cstddef>
#include <limits>
#include <vector>

#line 1 "library/algorithm/monotone_minima.hpp"



#line 7 "library/algorithm/monotone_minima.hpp"

namespace suisen {
    /**
     * @param n # rows
     * @param m # cols
     * @param compare (row, col1, col2 (> col1)) -> bool (= A(row, col1) <= A(row, col2))
     * @return std::vector<int> res s.t. res[i] = argmin_j f(i,j)
     */
    template <typename Compare, std::enable_if_t<std::is_invocable_r_v<bool, Compare, size_t, size_t, size_t>, std::nullptr_t> = nullptr>
    std::vector<int> monotone_minima(size_t n, size_t m, const Compare &compare) {
        std::vector<int> res(n);
        auto rec = [&](auto rec, size_t u, size_t d, size_t l, size_t r) -> void {
            if (u == d) return;
            assert(l < r);
            const size_t row = (u + d) >> 1;
            size_t argmin = l;
            for (size_t col = l + 1; col < r; ++col) if (not compare(row, argmin, col)) argmin = col;
            res[row] = argmin;
            rec(rec, u, row, l, argmin + 1);
            rec(rec, row + 1, d, argmin, r);
        };
        rec(rec, 0, n, 0, m);
        return res;
    }

    /**
     * @param n # rows
     * @param m # cols
     * @param matrix (row, col) -> value
     * @return std::vector<int> res s.t. res[i] = argmin_j f(i,j)
     */
    template <typename Matrix, std::enable_if_t<std::is_invocable_v<Matrix, size_t, size_t>, std::nullptr_t> = nullptr>
    std::vector<int> monotone_minima(size_t n, size_t m, const Matrix &matrix) {
        return monotone_minima(n, m, [&matrix](size_t i, size_t j1, size_t j2) { return matrix(i, j1) <= matrix(i, j2); });
    }
} // namespace suisen



#line 10 "library/convolution/min_plus_convolution.hpp"

namespace suisen {
    template <typename T>
    std::vector<T> min_plus_convolution_convex_convex(const std::vector<T> &a, const std::vector<T> &b) {
        const int n = a.size(), m = b.size();
        assert(n and m);
        // check if convex
        for (int i = 2; i < n; ++i) assert(a[i - 1] - a[i - 2] <= a[i] - a[i - 1]);
        // check if convex
        for (int j = 2; j < m; ++j) assert(b[j - 1] - b[j - 2] <= b[j] - b[j - 1]);

        std::vector<T> c(n + m - 1);
        c[0] = a[0] + b[0];
        for (int k = 0, i = 0; k < n + m - 2; ++k) {
            int j = k - i;
            if (j == m - 1 or (i < n - 1 and a[i + 1] + b[j] < a[i] + b[j + 1])) {
                c[k + 1] = a[++i] + b[j];
            } else {
                c[k + 1] = a[i] + b[++j];
            }
        }
        return c;
    }

    template <typename T>
    std::vector<T> min_plus_convolution_convex_arbitrary(const std::vector<T> &a, const std::vector<T> &b) {
        const int n = a.size(), m = b.size();
        assert(n and m);
        // check if convex
        for (int i = 2; i < n; ++i) assert(a[i - 1] - a[i - 2] <= a[i] - a[i - 1]);

        std::vector<int> argmin = monotone_minima(
            n + m - 1, m,
            [&](int k, int j1, int j2) {
                const int i1 = k - j1, i2 = k - j2;
                // upper right triangle
                if (i2 < 0) return true;
                // lower left triangle
                if (i1 >= n) return false;
                return a[i1] + b[j1] < a[i2] + b[j2];
            }
        );
        std::vector<T> c(n + m - 1);
        for (int k = 0; k < n + m - 1; ++k) {
            const int j = argmin[k], i = k - j;
            c[k] = a[i] + b[j];
        }
        return c;
    }
} // namespace suisen

Back to top page