Get tests/exceptions working

theochem · Oct 8, 2024 · 84ae74f · 84ae74f
1 parent f291ab0
commit 84ae74f
Show file tree

Hide file tree

Showing 6 changed files with 490 additions and 252 deletions.
diff --git a/include/permanent/common.h b/include/permanent/common.h
@@ -21,15 +21,13 @@ template <typename Type, typename IntType = void, typename = void>
 struct _result_t;
 
 template <typename Type>
-struct _result_t<Type, void,
-                 std::enable_if_t<std::is_scalar_v<Type> && !std::is_integral_v<Type>, void>>
+struct _result_t<Type, void, std::enable_if_t<std::is_scalar_v<Type>, void>>
 {
   typedef double type;
 };
 
 template <typename Type>
-struct _result_t<std::complex<Type>, void,
-                 std::enable_if_t<std::is_scalar_v<Type> && !std::is_integral_v<Type>, void>>
+struct _result_t<std::complex<Type>, void, std::enable_if_t<std::is_scalar_v<Type>, void>>
 {
   typedef std::complex<double> type;
 };

diff --git a/include/permanent/glynn.h b/include/permanent/glynn.h
@@ -14,6 +14,159 @@ result_t<T, I> glynn_square(const size_t m, const size_t n, const T *ptr)
 {
   (void)n;
 
+  /* Fill delta array (all +1 to start), and permutation array with [0...m]. */
+
+  sgn_t delta[64];
+  size_t perm[64 + 1];
+  for (size_t i = 0; i < m; ++i) {
+    perm[i] = i;
+    delta[i] = 1;
+  }
+
+  /* Handle first permutation. */
+
+  // result_t<T, I> vec[64];
+  result_t<T, I> out = 1;
+  for (size_t j = 0; j < m; ++j) {
+    result_t<T, I> sum = 0;
+    for (size_t i = 0; i < m; ++i) {
+      sum += ptr[i * m + j] * delta[i];
+    }
+    out *= sum;
+    // vec[j] = sum;
+  }
+  // for (size_t j = 0; j < m; ++j) {
+  //     out *= vec[j];
+  // }
+
+  /* Iterate from the second to the final permutation. */
+
+  std::size_t bound = m - 1;
+  std::size_t pos = 0;
+  sgn_t sign = 1;
+  while (pos != bound) {
+    /* Update sign and delta. */
+
+    sign *= -1;
+    delta[bound - pos] *= -1;
+
+    /* Compute term. */
+
+    result_t<T, I> prod = 1;
+    for (size_t j = 0; j < m; ++j) {
+      result_t<T, I> sum = 0;
+      for (size_t i = 0; i < m; ++i) {
+        sum += ptr[i * m + j] * delta[i];
+      }
+      // vec[j] = sum;
+      prod *= sum;
+    }
+
+    /* Multiply by the product of the vectors in delta. */
+
+    // for (size_t i = 0; i < m; ++i) {
+    //     prod *= vec[i];
+    // }
+    out += sign * prod;
+
+    /* Go to next permutation. */
+
+    perm[0] = 0;
+    perm[pos] = perm[pos + 1];
+    ++pos;
+    perm[pos] = pos;
+    pos = perm[0];
+  }
+
+  /* Divide by external factor and return permanent. */
+
+  return out / (1UL << bound);
+}
+
+template <typename T, typename I = void>
+result_t<T, I> glynn_rectangular(const size_t m, const size_t n, const T *ptr)
+{
+  /* Fill delta array (all +1 to start), and permutation array with [0...n].
+   */
+
+  sgn_t delta[64];
+  size_t perm[64 + 1];
+  for (size_t i = 0; i < n; ++i) {
+    delta[i] = 1;
+    perm[i] = i;
+  }
+
+  /* Handle first permutation. */
+
+  result_t<T, I> out = 1;
+  for (size_t j = 0; j < n; ++j) {
+    result_t<T, I> sum = 0;
+    for (size_t i = 0; i < m; ++i) {
+      sum += ptr[i * n + j] * delta[i];
+    }
+    for (size_t k = m; k < n; ++k) {
+      sum += delta[k];
+    }
+    // vec[j] = sum;
+    out *= sum;
+  }
+  // for (i = 0; i < n; ++i) {
+  //     out *= vec[i];
+  // }
+
+  /* Iterate from the second to the final permutation. */
+
+  size_t pos = 0;
+  size_t bound = n - 1;
+  sgn_t sign = 1;
+  while (pos != bound) {
+    /* Update sign and delta. */
+
+    sign *= -1;
+    delta[bound - pos] *= -1;
+
+    /* Compute term. */
+
+    result_t<T, I> prod = 1;
+    for (size_t j = 0; j < n; ++j) {
+      result_t<T, I> sum = 0;
+      for (size_t i = 0; i < m; ++i) {
+        sum += ptr[i * n + j] * delta[i];
+      }
+
+      for (size_t k = m; k < n; ++k) {
+        sum += delta[k];
+      }
+      prod *= sum;
+    }
+
+    /* Multiply by the product of the vectors in delta. */
+
+    // T prod = 1.0;
+    // for (i = 0; i < n; ++i) {
+    //     prod *= vec[i];
+    // }
+    out += sign * prod;
+
+    /* Go to next permutation. */
+
+    perm[0] = 0;
+    perm[pos] = perm[pos + 1];
+    ++pos;
+    perm[pos] = pos;
+    pos = perm[0];
+  }
+
+  /* Divide by external factor and return permanent. */
+
+  return (out / (1UL << bound)) / factorial<result_t<T, I>>(n - m);
+}
+
+template <typename T, typename I = void>
+result_t<T, I> glynn_square1(const size_t m, const size_t n, const T *ptr)
+{
+  (void)n;
+
   // Fill delta array ([+1...+1]) and permutation array ([0...m])
   size_t perm[64 + 1];
   sgn_t delta[64];
@@ -72,7 +225,7 @@ result_t<T, I> glynn_square(const size_t m, const size_t n, const T *ptr)
 }
 
 template <typename T, typename I = void>
-result_t<T, I> glynn_rectangular(const size_t m, const size_t n, const T *ptr)
+result_t<T, I> glynn_rectangular1(const size_t m, const size_t n, const T *ptr)
 {
   // Fill delta array ([+1...+1]) and permutation array ([0...m])
   size_t perm[64 + 1];

diff --git a/include/permanent/tuning.default.h b/include/permanent/tuning.default.h
@@ -8,10 +8,10 @@ namespace permanent {
 template <typename Type, typename IntType = void>
 struct _tuning_params_t
 {
-  static constexpr double PARAM_1 = +0.000000000e+00;
-  static constexpr double PARAM_2 = +0.000000000e+00;
-  static constexpr double PARAM_3 = +0.000000000e+00;
-  static constexpr double PARAM_4 = +0.000000000e+00;
+  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;
 };
 
 template <typename Type, typename IntType = void>

diff --git a/src/main.cc b/src/main.cc
@@ -259,11 +259,11 @@ static PyObject *py_opt(PyObject *module, PyObject *object)
 
   size_t m = PyArray_DIMS(matrix)[0];
   size_t n = PyArray_DIMS(matrix)[1];
-  if (m > n) {
-    return py_opt(module, PyArray_Transpose(matrix, nullptr));
-  } else if (m > 64 || n > 64) {
+  if (m > 64 || n > 64) {
     PyErr_SetString(PyExc_ValueError, "Argument array must have <=64 rows and columns");
     return nullptr;
+  } else if (m > n) {
+    return py_opt(module, PyArray_Transpose(matrix, nullptr));
   }
 
   PERMANENT_DISPATCH_ARRAY_TYPE(permanent::opt, matrix)
@@ -282,7 +282,10 @@ static PyObject *py_combinatoric(PyObject *module, PyObject *object)
 
   size_t m = PyArray_DIMS(matrix)[0];
   size_t n = PyArray_DIMS(matrix)[1];
-  if (m > n) {
+  if (m > 64 || n > 64) {
+    PyErr_SetString(PyExc_ValueError, "Argument array must have <=64 rows and columns");
+    return nullptr;
+  } else if (m > n) {
     return py_combinatoric(module, PyArray_Transpose(matrix, nullptr));
   }
 
@@ -302,7 +305,10 @@ static PyObject *py_glynn(PyObject *module, PyObject *object)
 
   size_t m = PyArray_DIMS(matrix)[0];
   size_t n = PyArray_DIMS(matrix)[1];
-  if (m > n) {
+  if (m > 64 || n > 64) {
+    PyErr_SetString(PyExc_ValueError, "Argument array must have <=64 rows and columns");
+    return nullptr;
+  } else if (m > n) {
     return py_glynn(module, PyArray_Transpose(matrix, nullptr));
   }
 
@@ -320,7 +326,10 @@ static PyObject *py_ryser(PyObject *module, PyObject *object)
 
   size_t m = PyArray_DIMS(matrix)[0];
   size_t n = PyArray_DIMS(matrix)[1];
-  if (m > n) {
+  if (m > 64 || n > 64) {
+    PyErr_SetString(PyExc_ValueError, "Argument array must have <=64 rows and columns");
+    return nullptr;
+  } else if (m > n) {
     return py_ryser(module, PyArray_Transpose(matrix, nullptr));
   }
 

diff --git a/src/tuning.cc b/src/tuning.cc
@@ -54,19 +54,19 @@ int generate_tuning_data(std::ofstream &csv_file, const size_t m, const size_t n
   if (m == n) {
     for (size_t i = 0; i != NUM_TRIALS; ++i) {
       begin = std::clock();
-      soln_combn = permanent::combinatoric<T, I>(m, n, array);
+      soln_combn = permanent::combinatoric_square<T, I>(m, n, array);
       ensure(soln_combn);
       end = std::clock();
       time_combn[i] = static_cast<double>(end - begin);
 
       begin = std::clock();
-      soln_glynn = permanent::glynn<T, I>(m, n, array);
+      soln_glynn = permanent::glynn_square<T, I>(m, n, array);
       ensure(soln_glynn);
       end = std::clock();
       time_glynn[i] = static_cast<double>(end - begin);
 
       begin = std::clock();
-      soln_ryser = permanent::ryser<T, I>(m, n, array);
+      soln_ryser = permanent::ryser_square<T, I>(m, n, array);
       ensure(soln_ryser);
       end = std::clock();
       time_ryser[i] = static_cast<double>(end - begin);