From 8bf58e729cfa515ed4f925b4bf777889207f4b16 Mon Sep 17 00:00:00 2001 From: cassmasschelein Date: Sat, 23 Nov 2024 19:24:51 -0700 Subject: [PATCH] fixed both ryser algs, added tests, updated tuning --- Makefile | 8 +- include/permanent/opt.h | 49 ++++-- include/permanent/ryser.h | 185 +++++++++----------- include/permanent/tuning.default.h | 29 +++- pyproject.toml | 3 + src/tuning.cc | 127 +++++++++----- tests/test_permanent.py | 264 ++++------------------------- tools/generate_tuning_header.py | 147 ++++++++-------- 8 files changed, 350 insertions(+), 462 deletions(-) diff --git a/Makefile b/Makefile index 213f262..7c4f05f 100644 --- a/Makefile +++ b/Makefile @@ -26,9 +26,9 @@ clean: install: _build @cmake --install build +build: + @cmake -B build $(CMAKE_FLAGS) + _build: build @cmake --build build --target all - @test -e build/compile_commands.json && ln -sf build/compile_commands.json . - -build: - @cmake -B build $(CMAKE_FLAGS) . + @if test -e build/compile_commands.json; then ln -sf build/compile_commands.json .; fi diff --git a/include/permanent/opt.h b/include/permanent/opt.h index 00d42b0..24e3b00 100644 --- a/include/permanent/opt.h +++ b/include/permanent/opt.h @@ -20,22 +20,49 @@ namespace permanent { template result_t opt_square(const size_t m, const size_t n, const T *ptr) { - if (n <= PARAM_4) { - return ryser_square(m, n, ptr); - } else if (n * PARAM_1 + PARAM_2 + PARAM_3 > 0) { - return combinatoric_square(m, n, ptr); - } else { - return glynn_square(m, n, ptr); + if (n <= PARAM_8) { // N ≤ 13 + if (n <= PARAM_4) { + return combinatoric_square(m, n, ptr); + } else { + // Check first hyperplane: -6.67(m/n) + (-0.35)(n) + 6.49 = 0 + // Using stored parameters for flexibility + const double ratio = static_cast(m) / n; + if (PARAM_1 * ratio + PARAM_2 * n + PARAM_3 > 0) { + return combinatoric_square(m, n, ptr); + } else { + return glynn_square(m, n, ptr); + } + } + } else { // N > 13 + // Check second hyperplane: -10.76(m/n) + 0.09(n) + 1.96 = 0 + const double ratio = static_cast(m) / n; + if (PARAM_5 * ratio + PARAM_6 * n + PARAM_7 > 0) { + return glynn_square(m, n, ptr); + } else { + return ryser_square(m, n, ptr); + } } } template result_t opt_rectangular(const size_t m, const size_t n, const T *ptr) { - if (n * PARAM_1 + PARAM_2 * m / n + PARAM_3 > 0) { - return combinatoric_rectangular(m, n, ptr); - } else { - return glynn_rectangular(m, n, ptr); + if (n <= PARAM_8) { // N ≤ 13 + // First hyperplane for small matrices + const double ratio = static_cast(m) / n; + if (PARAM_1 * ratio + PARAM_2 * n + PARAM_3 > 0) { + return combinatoric_rectangular(m, n, ptr); + } else { + return glynn_rectangular(m, n, ptr); + } + } else { // N > 13 + // Second hyperplane for large matrices + const double ratio = static_cast(m) / n; + if (PARAM_5 * ratio + PARAM_6 * n + PARAM_7 > 0) { + return glynn_rectangular(m, n, ptr); + } else { + return ryser_rectangular(m, n, ptr); + } } } @@ -47,4 +74,4 @@ result_t opt(const size_t m, const size_t n, const T *ptr) } // namespace permanent -#endif // permanent_opt_h_ +#endif // permanent_opt_h_ \ No newline at end of file diff --git a/include/permanent/ryser.h b/include/permanent/ryser.h index aa03ab8..2deb6e2 100644 --- a/include/permanent/ryser.h +++ b/include/permanent/ryser.h @@ -8,9 +8,40 @@ #include #include +#include namespace { +/* Helper functions for rectangle ryser. */ +inline int gray_code (int n) { + return n ^ (n >> 1); +} + +inline int count_bits (int n) { + int res = 0; + for (; n; n >>= 1) + res += n & 1; + return res; +} + +inline void all_combinations(int n, int k, std::vector >& k_perms) { + k_perms.clear(); // Clear k_perms to make sure it starts empty + + for (int i = 0; i < (1 << n); i++) { + int cur = gray_code(i); + if (count_bits(cur) == k) { + std::vector vec; // Initialize vector for new combination + + for (int j = 0; j < n; j++) { + if (cur & (1 << j)) { + vec.push_back(j); // Add element to vector + } + } + k_perms.push_back(vec); // Store the combination + } + } +} + template auto parity(const T &x) { @@ -95,120 +126,72 @@ namespace permanent { template result_t ryser_square(const size_t m, const size_t n, const T *ptr) { - (void)n; - - size_t i, j; - size_t k; - result_t rowsum; - result_t out = 0; - - /* Iterate over c = pow(2, m) submatrices (equal to (1 << m)) submatrices. - */ - - size_t c = 1UL << m; - - /* Loop over columns of submatrix; compute product of row sums. */ - - for (k = 0; k < c; ++k) { - result_t rowsumprod = 1.0; - for (i = 0; i < m; ++i) { - /* Loop over rows of submatrix; compute row sum. */ - - rowsum = 0.0; - for (j = 0; j < m; ++j) { - /* Add element to row sum if the row index is in the - * characteristic * vector of the submatrix, which is the binary - * vector given by k. */ - - if (k & (1UL << j)) { - rowsum += ptr[m * i + j]; + (void)n; + + size_t i, j, k; + T rowsum; // Use T instead of result_t for intermediate calculations + T out = 0; + size_t c = 1UL << m; + + for (k = 0; k < c; ++k) { + T rowsumprod = 1.0; + for (i = 0; i < m; ++i) { + rowsum = 0.0; + for (j = 0; j < m; ++j) { + if (k & (1UL << j)) { + rowsum += ptr[m * i + j]; + } + } + rowsumprod *= rowsum; } - } - - /* Update product of row sums. */ - - rowsumprod *= rowsum; + out += rowsumprod * (1 - ((__builtin_popcountll(k) & 1) << 1)); } - /* Add term multiplied by the parity of the characteristic vector. */ - - out += rowsumprod * parity(std::popcount(k)); - } - - /* Return answer with the correct sign (times -1 for odd m). */ - - return out * parity(m); + return static_cast>(out * ((m % 2 == 1) ? -1 : 1)); } template result_t ryser_rectangular(const size_t m, const size_t n, const T *ptr) { - size_t falling_fact[128]; - size_t perm[128]; - size_t inv_perm[128]; - falling_fact[0] = 0; - - size_t i, j, k; - - int value_sign = 1; - - result_t sum_of_matrix_vals; - result_t sum_over_k_vals = 0.0; - result_t vec[64]; - - /* Dealing with a rectangle. Can't use bit hacking trick here. */ - for (k = 0; k < m; ++k) { - /* Store the binomial coefficient for this k value bin_c. */ - size_t counter = 0; // Count how many permutations you have generated - size_t bin_c = binomial((n - m + k), k); - result_t prod_of_cols = 1; - result_t result = 0; - - /* (Re)initialize the set to permute for this k value. */ - init_perm(n, falling_fact, perm, inv_perm); + size_t i, j, k, bin; + int sign = 1; + + result_t colprod, matsum, permsum; + result_t out = 0.0; + result_t vec[64]; + std::vector > k_perms; + + /* Iterate over subsets from size 0 to size m */ + + for (k = 0; k < m; ++k) { + all_combinations(n, m - k, k_perms); + + bin = binomial((n - m + k), k); + permsum = 0.0; + + /* Compute permanents of each submatrix */ + for (const auto& combination : k_perms) { + for (i = 0; i < m; ++i) { + matsum = 0.0; + for (j = 0; j < (m - k); ++j) + matsum += ptr[i * n + combination[j]]; + vec[i] = matsum; + } + + colprod = 1.0; + for (i = 0; i < m; ++i) + colprod *= vec[i]; + permsum += colprod * sign * bin; + } - /* sort up to position u + 1 where u = min(m - k, n - 1). */ - size_t sort_up_to = n - 1; + /* Add term to result */ - if ((m - k) < sort_up_to) { - sort_up_to = (m - k); - } + out += permsum; - for (i = 0; i < m; ++i) { - sum_of_matrix_vals = 0.0; - for (j = 0; j < sort_up_to; ++j) { - sum_of_matrix_vals += ptr[i * n + perm[j]]; - } - vec[i] = sum_of_matrix_vals; - } - for (i = 0; i < m; ++i) { - prod_of_cols *= vec[i]; + sign *= -1; } - result += value_sign * bin_c * prod_of_cols; - counter += 1; - - /* Iterate over second to last permutations of the set. */ - while (gen_next_perm(falling_fact, perm, inv_perm, n, sort_up_to)) { - for (i = 0; i < m; ++i) { - sum_of_matrix_vals = 0.0; - for (j = 0; j < m - k; ++j) { - sum_of_matrix_vals += ptr[i * n + perm[j]]; - } - vec[i] = sum_of_matrix_vals; - } - prod_of_cols = 1.0; - for (i = 0; i < m; ++i) { - prod_of_cols *= vec[i]; - } - - result += value_sign * bin_c * prod_of_cols; - counter += 1; - } - sum_over_k_vals += (result / counter) * binomial(n, (m - k)); - value_sign *= -1; - } - return sum_over_k_vals; + return out; } template diff --git a/include/permanent/tuning.default.h b/include/permanent/tuning.default.h index a73d11b..ac8e9a1 100644 --- a/include/permanent/tuning.default.h +++ b/include/permanent/tuning.default.h @@ -8,10 +8,17 @@ namespace permanent { template struct _tuning_params_t { - static constexpr double PARAM_1 = -5.72098e-01; - static constexpr double PARAM_2 = -2.20142e+01; - static constexpr double PARAM_3 = +1.52979e+01; - static constexpr double PARAM_4 = +3.00000e+00; + // First hyperplane parameters (N ≤ 13): -6.67(m/n) + (-0.35)(n) + 6.49 = 0 + static constexpr double PARAM_1 = -6.67000e+00; // coefficient of m/n + static constexpr double PARAM_2 = -3.50000e-01; // coefficient of n + static constexpr double PARAM_3 = +6.49000e+00; // constant term + static constexpr double PARAM_4 = +4.00000e+00; // Combinatorial crossover limit + + // Second hyperplane parameters (N > 13): -10.76(m/n) + 0.09(n) + 1.96 = 0 + static constexpr double PARAM_5 = -1.07600e+01; // coefficient of m/n + static constexpr double PARAM_6 = +9.00000e-02; // coefficient of n + static constexpr double PARAM_7 = +1.96000e+00; // constant term + static constexpr double PARAM_8 = +1.30000e+01; // Combinatorial limit (N=13) }; template @@ -26,6 +33,18 @@ static constexpr double PARAM_3 = _tuning_params_t::PARAM_3; template static constexpr double PARAM_4 = _tuning_params_t::PARAM_4; +template +static constexpr double PARAM_5 = _tuning_params_t::PARAM_5; + +template +static constexpr double PARAM_6 = _tuning_params_t::PARAM_6; + +template +static constexpr double PARAM_7 = _tuning_params_t::PARAM_7; + +template +static constexpr double PARAM_8 = _tuning_params_t::PARAM_8; + } // namespace permanent -#endif // permanent_tuning_h_ +#endif // permanent_tuning_h_ \ No newline at end of file diff --git a/pyproject.toml b/pyproject.toml index 09d5480..abb624b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -38,6 +38,9 @@ test = ["pytest"] [tool.scikit-build] wheel.expand-macos-universal-tags = true +[tool.scikit-build] +cmake.args = ["-DPERMANENT_TUNE=ON"] + [tool.pytest.ini_options] minversion = "6.0" addopts = ["-ra", "--showlocals", "--strict-markers", "--strict-config"] diff --git a/src/tuning.cc b/src/tuning.cc index b7d6f64..5617924 100644 --- a/src/tuning.cc +++ b/src/tuning.cc @@ -28,14 +28,16 @@ constexpr char CSV_HEADER[] = "M/N,N,Combn,Glynn,Ryser,Fastest"; constexpr size_t NUM_TRIALS = 5; -constexpr size_t MAX_ROWS = 16; +constexpr size_t MAX_ROWS = 26; -constexpr size_t MAX_COLS = 16; +constexpr size_t MAX_COLS = 26; constexpr size_t DATA_POINTS = MAX_ROWS * MAX_COLS; constexpr double TOLERANCE = 0.0001; +constexpr size_t RUN_COMBINATORIAL_UNTIL = 13; + namespace { template @@ -51,13 +53,18 @@ int generate_tuning_data(std::ofstream &csv_file, const size_t m, const size_t n // Solve the permanent using each algorithm NUM_TRIALS number of times result_t soln_combn, soln_glynn, soln_ryser; std::clock_t begin, end; + bool run_combinatorial = (m <= RUN_COMBINATORIAL_UNTIL); if (m == n) { for (size_t i = 0; i != NUM_TRIALS; ++i) { - begin = std::clock(); - soln_combn = permanent::combinatoric_square(m, n, array); - ensure(soln_combn); - end = std::clock(); - time_combn[i] = static_cast(end - begin); + if (run_combinatorial) { + begin = std::clock(); + soln_combn = permanent::combinatoric_square(m, n, array); + ensure(soln_combn); + end = std::clock(); + time_combn[i] = static_cast(end - begin); + } else { + time_combn[i] = std::numeric_limits::infinity(); + } begin = std::clock(); soln_glynn = permanent::glynn_square(m, n, array); @@ -71,29 +78,46 @@ int generate_tuning_data(std::ofstream &csv_file, const size_t m, const size_t n end = std::clock(); time_ryser[i] = static_cast(end - begin); - if ((std::fabs(soln_combn - soln_glynn) / - std::min(std::fabs(soln_combn), std::fabs(soln_glynn)) > - TOLERANCE) || - (std::fabs(soln_ryser - soln_combn) / - std::min(std::fabs(soln_ryser), std::fabs(soln_combn)) > - TOLERANCE) || - (std::fabs(soln_glynn - soln_ryser) / - std::min(std::fabs(soln_glynn), std::fabs(soln_ryser)) > - TOLERANCE)) { - std::cerr << "Bad permanent values:" - << "\nCombn: " << soln_combn << "\nGlynn: " << soln_glynn - << "\nRyser: " << soln_ryser << std::endl; - csv_file.close(); - return 3; + if (run_combinatorial) { + if ((std::fabs(soln_combn - soln_glynn) / + std::min(std::fabs(soln_combn), std::fabs(soln_glynn)) > + TOLERANCE) || + (std::fabs(soln_ryser - soln_combn) / + std::min(std::fabs(soln_ryser), std::fabs(soln_combn)) > + TOLERANCE) || + (std::fabs(soln_glynn - soln_ryser) / + std::min(std::fabs(soln_glynn), std::fabs(soln_ryser)) > + TOLERANCE)) { + std::cerr << "Bad permanent values:" + << "\nCombn: " << soln_combn << "\nGlynn: " << soln_glynn + << "\nRyser: " << soln_ryser << std::endl; + csv_file.close(); + return 3; + } + } else { + // Only check Glynn vs Ryser when not running combinatorial + if (std::fabs(soln_glynn - soln_ryser) / + std::min(std::fabs(soln_glynn), std::fabs(soln_ryser)) > + TOLERANCE) { + std::cerr << "Bad permanent values:" + << "\nGlynn: " << soln_glynn + << "\nRyser: " << soln_ryser << std::endl; + csv_file.close(); + return 3; + } } } } else { for (size_t i = 0; i != NUM_TRIALS; ++i) { - begin = std::clock(); - soln_combn = permanent::combinatoric_rectangular(m, n, array); - ensure(soln_combn); - end = std::clock(); - time_combn[i] = static_cast(end - begin); + if (run_combinatorial) { + begin = std::clock(); + soln_combn = permanent::combinatoric_rectangular(m, n, array); + ensure(soln_combn); + end = std::clock(); + time_combn[i] = static_cast(end - begin); + } else { + time_combn[i] = std::numeric_limits::infinity(); + } begin = std::clock(); soln_glynn = permanent::glynn_rectangular(m, n, array); @@ -107,20 +131,33 @@ int generate_tuning_data(std::ofstream &csv_file, const size_t m, const size_t n end = std::clock(); time_ryser[i] = static_cast(end - begin); - if ((std::fabs(soln_combn - soln_glynn) / - std::min(std::fabs(soln_combn), std::fabs(soln_glynn)) > - TOLERANCE) || - (std::fabs(soln_ryser - soln_combn) / - std::min(std::fabs(soln_ryser), std::fabs(soln_combn)) > - TOLERANCE) || - (std::fabs(soln_glynn - soln_ryser) / - std::min(std::fabs(soln_glynn), std::fabs(soln_ryser)) > - TOLERANCE)) { - std::cerr << std::scientific << "Bad permanent values:" - << "\nCombn: " << soln_combn << "\nGlynn: " << soln_glynn - << "\nRyser: " << soln_ryser << std::endl; - csv_file.close(); - return 3; + if (run_combinatorial) { + if ((std::fabs(soln_combn - soln_glynn) / + std::min(std::fabs(soln_combn), std::fabs(soln_glynn)) > + TOLERANCE) || + (std::fabs(soln_ryser - soln_combn) / + std::min(std::fabs(soln_ryser), std::fabs(soln_combn)) > + TOLERANCE) || + (std::fabs(soln_glynn - soln_ryser) / + std::min(std::fabs(soln_glynn), std::fabs(soln_ryser)) > + TOLERANCE)) { + std::cerr << std::scientific << "Bad permanent values:" + << "\nCombn: " << soln_combn << "\nGlynn: " << soln_glynn + << "\nRyser: " << soln_ryser << std::endl; + csv_file.close(); + return 3; + } + } else { + // Only check Glynn vs Ryser when not running combinatorial + if (std::fabs(soln_glynn - soln_ryser) / + std::min(std::fabs(soln_glynn), std::fabs(soln_ryser)) > + TOLERANCE) { + std::cerr << std::scientific << "Bad permanent values:" + << "\nGlynn: " << soln_glynn + << "\nRyser: " << soln_ryser << std::endl; + csv_file.close(); + return 3; + } } } } @@ -130,11 +167,17 @@ int generate_tuning_data(std::ofstream &csv_file, const size_t m, const size_t n double mean_glynn = 0; double mean_ryser = 0; for (size_t i = 0; i != NUM_TRIALS; ++i) { - mean_combn += time_combn[i]; + if (run_combinatorial) { + mean_combn += time_combn[i]; + } mean_glynn += time_glynn[i]; mean_ryser += time_ryser[i]; } - mean_combn /= NUM_TRIALS; + if (run_combinatorial) { + mean_combn /= NUM_TRIALS; + } else { + mean_combn = std::numeric_limits::infinity(); + } mean_glynn /= NUM_TRIALS; mean_ryser /= NUM_TRIALS; diff --git a/tests/test_permanent.py b/tests/test_permanent.py index ee53532..e99012f 100644 --- a/tests/test_permanent.py +++ b/tests/test_permanent.py @@ -76,242 +76,42 @@ def test_int_dim_diff_gt_20_raises(fn, arg): np.array([[chr(x) for x in range(65, 69)] for _ in range(4)]), ], ) + @pytest.mark.parametrize("fn", FNS) def test_invalid_type_raises(fn, arg): r"""Ensure that matrices with an invalid dtype raise a ValueError.""" with npt.assert_raises(TypeError): fn(arg) +@pytest.mark.parametrize( + "args", + [ + (np.arange(1, 5, dtype=float).reshape(2, 2), 10), + (np.arange(1, 10, dtype=float).reshape(3, 3), 450), + (np.arange(1, 17, dtype=float).reshape(4, 4), 55456), + (np.arange(1, 50, dtype=float).reshape(7, 7), 5373548250000), + # Rectangular matrices + (np.arange(1, 7, dtype=float).reshape(2, 3), 58), + (np.arange(1, 9, dtype=float).reshape(2, 4), 190), + (np.arange(1, 15, dtype=float).reshape(2, 7), 1820), + (np.arange(1, 36, dtype=float).reshape(5, 7), 1521238320), + (np.arange(1, 43, dtype=float).reshape(6, 7), 117681979920), + # Special matrices + (np.ones((10, 10), dtype=float), factorial(10)), + (np.ones((12, 12), dtype=float), factorial(12)), + (np.identity(10, dtype=float), 1), + (np.identity(5, dtype=float), 1), + (np.diag(np.diag(np.ones((3, 7), dtype=float))), 1), + ], +) -# # Test square matrices -# def test_2by2_comb(): -# matrix = np.arange(1, 5, dtype=np.double).reshape(2, 2) -# assert abs((permanent.combinatoric(matrix) - 10) / 10) <= REL_ERROR -# npt.assert_allclose( -# -# -# def test_2by2_glynn(): -# matrix = np.arange(1, 5, dtype=np.double).reshape(2, 2) -# assert abs((permanent.glynn(matrix) - 10) / 10) <= REL_ERROR -# -# -# def test_2by2_ryser(): -# matrix = np.arange(1, 5, dtype=np.double).reshape(2, 2) -# assert abs((permanent.ryser(matrix) - 10) / 10) <= REL_ERROR -# -# -# def test_3by3_comb(): -# matrix = np.arange(1, 10, dtype=np.double).reshape(3, 3) -# assert abs((permanent.combinatoric(matrix) - 450) / 450) <= REL_ERROR -# -# -# def test_3by3_glynn(): -# matrix = np.arange(1, 10, dtype=np.double).reshape(3, 3) -# assert abs((permanent.glynn(matrix) - 450) / 450) <= REL_ERROR -# -# -# def test_3by3_ryser(): -# matrix = np.arange(1, 10, dtype=np.double).reshape(3, 3) -# assert abs((permanent.ryser(matrix) - 450) / 450) <= REL_ERROR -# -# -# def test_4by4_comb(): -# matrix = np.arange(1, 17, dtype=np.double).reshape(4, 4) -# assert abs((permanent.combinatoric(matrix) - 55456) / 55456) <= REL_ERROR -# -# -# def test_4by4_glynn(): -# matrix = np.arange(1, 17, dtype=np.double).reshape(4, 4) -# assert abs((permanent.glynn(matrix) - 55456) / 55456) <= REL_ERROR -# -# -# def test_4by4_ryser(): -# matrix = np.arange(1, 17, dtype=np.double).reshape(4, 4) -# assert abs((permanent.ryser(matrix) - 55456) / 55456) <= REL_ERROR -# -# -# def test_7by7_comb(): -# matrix = np.arange(1, 50, dtype=np.double).reshape(7, 7) -# assert ( -# abs((permanent.combinatoric(matrix) - 5373548250000) / 5373548250000) -# <= REL_ERROR -# ) -# -# -# def test_7by7_glynn(): -# matrix = np.arange(1, 50, dtype=np.double).reshape(7, 7) -# assert abs((permanent.glynn(matrix) - 5373548250000) / 5373548250000) <= REL_ERROR -# -# -# def test_7by7_ryser(): -# matrix = np.arange(1, 50, dtype=np.double).reshape(7, 7) -# assert abs((permanent.ryser(matrix) - 5373548250000) / 5373548250000) <= REL_ERROR -# -# -# ## Test rectangular matrices -# def test_2by3_comb(): -# matrix = np.arange(1, 7, dtype=np.double).reshape(2, 3) -# assert abs((permanent.combinatoric(matrix) - 58) / 58) <= REL_ERROR -# -# -# def test_2by3_glynn(): -# matrix = np.arange(1, 7, dtype=np.double).reshape(2, 3) -# assert abs((permanent.glynn(matrix) - 58) / 58) <= REL_ERROR -# -# -# def test_2by3_ryser(): -# matrix = np.arange(1, 7, dtype=np.double).reshape(2, 3) -# assert abs((permanent.ryser(matrix) - 58) / 58) <= REL_ERROR -# -# -# def test_2by4_comb(): -# matrix = np.arange(1, 9, dtype=np.double).reshape(2, 4) -# assert abs((permanent.combinatoric(matrix) - 190) / 190) <= REL_ERROR -# -# -# def test_2by4_glynn(): -# matrix = np.arange(1, 9, dtype=np.double).reshape(2, 4) -# assert abs((permanent.glynn(matrix) - 190) / 190) <= REL_ERROR -# -# -# def test_2by4_ryser(): -# matrix = np.arange(1, 9, dtype=np.double).reshape(2, 4) -# assert abs((permanent.ryser(matrix) - 190) / 190) <= REL_ERROR -# -# -# def test_2by7_comb(): -# matrix = np.arange(1, 15, dtype=np.double).reshape(2, 7) -# assert abs((permanent.combinatoric(matrix) - 1820) / 1820) <= REL_ERROR -# -# -# def test_2by7_glynn(): -# matrix = np.arange(1, 15, dtype=np.double).reshape(2, 7) -# assert abs((permanent.glynn(matrix) - 1820) / 1820) <= REL_ERROR -# -# -# def test_2by7_ryser(): -# matrix = np.arange(1, 15, dtype=np.double).reshape(2, 7) -# assert abs((permanent.ryser(matrix) - 1820) / 1820) <= REL_ERROR -# -# -# def test_5by7_comb(): -# matrix = np.arange(1, 36, dtype=np.double).reshape(5, 7) -# assert abs((permanent.combinatoric(matrix) - 1521238320) / 1521238320) <= REL_ERROR -# -# -# def test_5by7_glynn(): -# matrix = np.arange(1, 36, dtype=np.double).reshape(5, 7) -# assert abs((permanent.glynn(matrix) - 1521238320) / 1521238320) <= REL_ERROR -# -# -# def test_5by7_ryser(): -# matrix = np.arange(1, 36, dtype=np.double).reshape(5, 7) -# assert abs((permanent.ryser(matrix) - 1521238320) / 1521238320) <= REL_ERROR -# -# -# def test_6by7_comb(): -# matrix = np.arange(1, 43, dtype=np.double).reshape(6, 7) -# assert ( -# abs((permanent.combinatoric(matrix) - 117681979920) / 117681979920) <= REL_ERROR -# ) -# -# -# def test_6by7_glynn(): -# matrix = np.arange(1, 43, dtype=np.double).reshape(6, 7) -# assert abs((permanent.glynn(matrix) - 117681979920) / 117681979920) <= REL_ERROR -# -# -# def test_6by7_ryser(): -# matrix = np.arange(1, 43, dtype=np.double).reshape(6, 7) -# assert abs((permanent.ryser(matrix) - 117681979920) / 117681979920) <= REL_ERROR -# -# -# def test_ones_comb(): -# matrix = np.ones((10, 10), dtype=np.double) -# assert ( -# abs((permanent.combinatoric(matrix) - factorial(10)) / factorial(10)) -# <= REL_ERROR -# ) -# -# -# def test_ones_ryser(): -# matrix = np.ones((10, 10), dtype=np.double) -# assert ( -# abs((permanent.ryser(matrix) - factorial(10)) / factorial(10)) -# <= REL_ERROR -# ) -# -# -# def test_ones_glynn(): -# matrix = np.ones((10, 10), dtype=np.double) -# assert ( -# abs((permanent.glynn(matrix) - factorial(10)) / factorial(10)) -# <= REL_ERROR -# ) -# -# -# def test_ones_comb_big(): -# matrix = np.ones((12, 12), dtype=np.double) -# assert ( -# abs((permanent.combinatoric(matrix) - factorial(12)) / factorial(12)) -# <= REL_ERROR -# ) -# -# -# def test_ones_ryser_big(): -# matrix = np.ones((12, 12), dtype=np.double) -# assert ( -# abs((permanent.ryser(matrix) - factorial(12)) / factorial(12)) -# <= REL_ERROR -# ) -# -# -# def test_ones_glynn_big(): -# matrix = np.ones((12, 12), dtype=np.double) -# assert ( -# abs((permanent.glynn(matrix) - factorial(12)) / factorial(12)) -# <= REL_ERROR -# ) -# -# -# def test_identity_comb(): -# matrix = np.identity(10, dtype=np.double) -# assert abs((permanent.combinatoric(matrix) - 1) / 1) <= REL_ERROR -# -# -# def test_identity_ryser(): -# matrix = np.identity(10, dtype=np.double) -# assert abs((permanent.ryser(matrix) - 1) / 1) <= REL_ERROR -# -# -# def test_identity_glynn(): -# matrix = np.identity(10, dtype=np.double) -# assert abs((permanent.glynn(matrix) - 1) / 1) <= REL_ERROR -# -# -# def test_identity_ryser_odd(): -# matrix = np.identity(5, dtype=np.double) -# assert abs((permanent.ryser(matrix) - 1) / 1) <= REL_ERROR -# -# -# def test_identity_glynn_odd(): -# matrix = np.identity(5, dtype=np.double) -# assert abs((permanent.glynn(matrix) - 1) / 1) <= REL_ERROR -# -# -# def test_identity_comb_diag(): -# matrix = np.ones((3, 7), dtype=np.double) -# diag_matrix = np.diag(np.diag(matrix)) -# assert abs((permanent.combinatoric(diag_matrix) - 1) / 1) <= REL_ERROR -# -# -# def test_identity_ryser_diag(): -# matrix = np.ones((3, 7), dtype=np.double) -# diag_matrix = np.diag(np.diag(matrix)) -# assert abs((permanent.ryser(diag_matrix) - 1) / 1) <= REL_ERROR -# -# -# def test_identity_glynn_diag(): -# matrix = np.ones((3, 7), dtype=np.double) -# diag_matrix = np.diag(np.diag(matrix)) -# assert abs((permanent.glynn(diag_matrix) - 1) / 1) <= REL_ERROR +@pytest.mark.parametrize("fn", FNS) +def test_compute_permanent_extended(fn, args): + """Test permanent computation for various matrix sizes and types.""" + npt.assert_allclose( + fn(args[0]), + args[1], + atol=ATOL, + rtol=RTOL, + err_msg=f"Failed for matrix shape {args[0].shape}, expected {args[1]}" + ) diff --git a/tools/generate_tuning_header.py b/tools/generate_tuning_header.py index aa82556..fdfd7b5 100644 --- a/tools/generate_tuning_header.py +++ b/tools/generate_tuning_header.py @@ -1,82 +1,79 @@ import sys - import numpy as np - import pandas as pd - from sklearn import svm from sklearn.metrics import accuracy_score - CSV_FILE = sys.argv[1] - HEADER_FILE = sys.argv[2] - # Read output from tuning program tuning = pd.read_csv(CSV_FILE, usecols=["M/N", "N", "Fastest"]) +# Split data into two parts based on N=13 +tuning_small = tuning[tuning["N"] <= 13].copy() +tuning_large = tuning[tuning["N"] > 13].copy() + # Update label columns for ML -tuning.rename(columns={"N": "x", "M/N": "y", "Fastest": "target"}, inplace=True) - -# Find Ryser limit and update to dual class for SVM -# Locate rows where column value matches specified value -matching_row = tuning.loc[tuning["target"] == 0, "x"] - -if not matching_row.empty: - # Pull out specific column value from the last matching row - ryser_limit = matching_row.iloc[-1].copy() -else: - ryser_limit = 0 - -update_target = tuning["target"] == 0 -tuning.loc[update_target, "target"] = -1 - -# Update classes to -1/1, Combn = 1 -update_glynn = tuning["target"] == 2 -tuning.loc[update_glynn, "target"] = -1 - -features = tuning[["x", "y"]] -label = tuning["target"] -value_counts = tuning["target"].value_counts() - -# Create train/test split -size = tuning.shape[0] -test_size = int(np.round(size * 0.1, 0)) - -# Split dataset into training and testing sets -x_train = features[:-test_size].values -y_train = label[:-test_size].values -x_test = features[-test_size:].values -y_test = label[-test_size:].values - -### Train a linear model - add a hard margin -linear_model = svm.SVC(kernel="linear", C=100.0) -linear_model.fit(x_train, y_train) - -# Create grid to evaluate model -xx = np.linspace(-1, max(features["x"]) + 1, len(x_train)) -yy = np.linspace(0, max(features["y"]) + 1, len(y_train)) -YY, XX = np.meshgrid(yy, xx) -xy = np.vstack([XX.ravel(), YY.ravel()]).T -train_size = len(features[:-test_size]["x"]) - -# Get the separating hyperplane -Z = linear_model.decision_function(xy).reshape(XX.shape) - -# Check the accuracy of the model -predictions_linear = linear_model.predict(x_test) -accuracy_linear = accuracy_score(y_test, predictions_linear) - -# Get the coefficients and bias -coefficients = linear_model.coef_[0] -bias = linear_model.intercept_[0] - -# Write a header file with constants defined as macros -param_1 = coefficients[0] -param_2 = coefficients[1] -param_3 = bias -param_4 = ryser_limit +tuning_small.rename(columns={"N": "x", "M/N": "y", "Fastest": "target"}, inplace=True) +tuning_large.rename(columns={"N": "x", "M/N": "y", "Fastest": "target"}, inplace=True) + +# Find Combinatoric Square limit +matching_row = tuning["Fastest"] == 1 +combn_limit = tuning.loc[matching_row, "N"].iloc[-1] if matching_row.any() else 0 + +# Update classes for first SVM (N <= 13, all three algorithms) +update_ryser = tuning_small["target"] == 0 +tuning_small.loc[update_ryser, "target"] = -1 +update_glynn = tuning_small["target"] == 2 +tuning_small.loc[update_glynn, "target"] = -1 + +# Update classes for second SVM (N > 13, Glynn vs Ryser) +tuning_large["target"] = np.where(tuning_large["target"] == 0, -1, 1) + +# Train first SVM (N <= 13) +features_small = tuning_small[["x", "y"]] +label_small = tuning_small["target"] +size_small = tuning_small.shape[0] +test_size_small = int(np.round(size_small * 0.1, 0)) + +x_train_small = features_small[:-test_size_small].values +y_train_small = label_small[:-test_size_small].values +x_test_small = features_small[-test_size_small:].values +y_test_small = label_small[-test_size_small:].values + +linear_model_small = svm.SVC(kernel="linear", C=100.0) +linear_model_small.fit(x_train_small, y_train_small) + +# Train second SVM (N > 13) +features_large = tuning_large[["x", "y"]] +label_large = tuning_large["target"] +size_large = tuning_large.shape[0] +test_size_large = int(np.round(size_large * 0.1, 0)) + +x_train_large = features_large[:-test_size_large].values +y_train_large = label_large[:-test_size_large].values +x_test_large = features_large[-test_size_large:].values +y_test_large = label_large[-test_size_large:].values + +linear_model_large = svm.SVC(kernel="linear", C=100.0) +linear_model_large.fit(x_train_large, y_train_large) + +# Get coefficients and bias for both models +coefficients_small = linear_model_small.coef_[0] +bias_small = linear_model_small.intercept_[0] +coefficients_large = linear_model_large.coef_[0] +bias_large = linear_model_large.intercept_[0] + +# Write header file with all parameters +param_1 = coefficients_small[0] # First hyperplane x coefficient +param_2 = coefficients_small[1] # First hyperplane y coefficient +param_3 = bias_small # First hyperplane bias +param_4 = combn_limit # Combinatoric square limit +param_5 = coefficients_large[0] # Second hyperplane x coefficient +param_6 = coefficients_large[1] # Second hyperplane y coefficient +param_7 = bias_large # Second hyperplane bias +param_8 = 13.0 # Combinatorial limit try: with open(HEADER_FILE, "w") as file_ptr: @@ -95,6 +92,10 @@ static constexpr double PARAM_2 = {param_2:+16.9e}; static constexpr double PARAM_3 = {param_3:+16.9e}; static constexpr double PARAM_4 = {param_4:+16.9e}; + static constexpr double PARAM_5 = {param_5:+16.9e}; + static constexpr double PARAM_6 = {param_6:+16.9e}; + static constexpr double PARAM_7 = {param_7:+16.9e}; + static constexpr double PARAM_8 = {param_8:+16.9e}; }}; template @@ -109,6 +110,18 @@ template static constexpr double PARAM_4 = _tuning_params_t::PARAM_4; +template +static constexpr double PARAM_5 = _tuning_params_t::PARAM_5; + +template +static constexpr double PARAM_6 = _tuning_params_t::PARAM_6; + +template +static constexpr double PARAM_7 = _tuning_params_t::PARAM_7; + +template +static constexpr double PARAM_8 = _tuning_params_t::PARAM_8; + }} // namespace permanent #endif // permanent_tuning_h_ @@ -121,4 +134,4 @@ except Exception as e: print("Error occurred:", e) - exit(1) + exit(1) \ No newline at end of file