diff --git a/.github/ACKNOWLEDGMENTS.md b/.github/ACKNOWLEDGMENTS.md index 083acdd56..a706c5e88 100644 --- a/.github/ACKNOWLEDGMENTS.md +++ b/.github/ACKNOWLEDGMENTS.md @@ -48,4 +48,12 @@ * [Ali Asadi](https://github.com/maliasadi) (Western University) - :thread: Commander of threads -* [Sebastián Duque](https://github.com/sduquemesa) (Xanadu) - 🎧 Quantum sound explorer \ No newline at end of file +* [Sebastián Duque](https://github.com/sduquemesa) (Xanadu) - 🎧 Quantum sound explorer + +* [Jiaqi Zhao](https://github.com/JQZ1111) (Polytechnique Montréal) - Photon squeezer + +* [Dominic Leclerc](https://github.com/dleclerc33) (Polytechnique de Montréal) :fleur_de_lis: King of vaccum + +* [Benjamin Lanthier](https://github.com/benjaminlanthier) (Polytechnique de Montréal) :mage: Warlock of eigenvalues + +* [Brandon Turcotte](https://github.com/brandonpolymtl) (Polytechnique de Montréal) - :star2: The Quantum-plators diff --git a/.github/CHANGELOG.md b/.github/CHANGELOG.md index 47085736e..1b2396304 100644 --- a/.github/CHANGELOG.md +++ b/.github/CHANGELOG.md @@ -4,6 +4,9 @@ ### Improvements +* Permanent algorithms are implemented in numba, hence reducing the C++ dependencies of The Walrus. [#300](https://github.com/XanaduAI/thewalrus/pull/300) + + ### Bug fixes ### Breaking changes @@ -12,6 +15,8 @@ This release contains contributions from (in alphabetical order): +Benjamin Lanthier, Dominic Leclerc, Nicolas Quesada, Brandon Turcotte, Jiaqi Zhao + --- # Version 0.17.0 diff --git a/include/libwalrus.hpp b/include/libwalrus.hpp index 2d82eabf2..abe586574 100644 --- a/include/libwalrus.hpp +++ b/include/libwalrus.hpp @@ -16,7 +16,6 @@ #include #include #include -#include /** * @namespace libwalrus diff --git a/include/permanent.hpp b/include/permanent.hpp deleted file mode 100644 index 32490b9d3..000000000 --- a/include/permanent.hpp +++ /dev/null @@ -1,598 +0,0 @@ -// Copyright 2019 Xanadu Quantum Technologies Inc. - -// Licensed under the Apache License, Version 2.0 (the "License"); -// you may not use this file except in compliance with the License. -// You may obtain a copy of the License at - -// http://www.apache.org/licenses/LICENSE-2.0 - -// Unless required by applicable law or agreed to in writing, software -// distributed under the License is distributed on an "AS IS" BASIS, -// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -// See the License for the specific language governing permissions and -// limitations under the License. -#pragma once -#include -#include -#include -#include "fsum.hpp" - -typedef unsigned long long int ullint; -typedef long long int llint; -//typedef long double qp; - -#if (defined(__GNUC__) || defined(__GNUG__)) && !(defined(__clang__) || defined(__INTEL_COMPILER)) -typedef __float128 qp; -//#include -#else -typedef long double qp; -#endif - -/** - * Gray code generator. - * - * @param n - * @return the corresponding Gray code - */ -static inline llint igray(llint n) -{ - /* Right Shift the number by 1 - taking xor with original number */ - return n ^ (n >> 1); -} - -/** - * Returns left most set bit. - * - * @param n - * @return the left most set bit - */ -static inline int left_most_set_bit(llint n) -{ - if (n == 0) - return 0; - - int msb = 0; - - while (n != 0) { - n = n / 2; - msb++; - } - return msb; -} - - -/** - * Decimal to binary conversion - * - * @param \f$k\f$ - * @param \f$n\f$ - * @return \f$n\f$ bit binary representation of integer \f$k\f$ - */ -static inline std::vector dec2bin(llint &k, int &n) -{ - llint kk = k; - int i = n; - std::vector mat(n, 0); - - while ( kk > 0 && i > 0 ) { - mat[i-1] = kk % 2; - kk = kk/2; - i = i-1; - } - return mat; -} - -/** - * Get the next ordering index - * - * @param \f$l\f$ - * @param \f$k\f$ - * @return the \f$k+1\f$-th ordering index with updating \f$l\f$ from init index \f$k\f$ - */ -static inline size_t next_perm_ordering_index(std::vector &l, size_t k) -{ - l[0] = 0; - l[k] = l[k+1]; - l[k+1] = k+1; - return l[0]; -} - -namespace libwalrus { - -/** - * Returns the permanent of a matrix. - * - * \rst - * - * Returns the permanent of a matrix using Ryser's algorithm - * with Gray code ordering, that has a time-complexity of \f$O(n 2^n)\f$ - * - * \endrst - * - * - * @tparam T type of the matrix data - * @param mat a flattened vector of size \f$n^2\f$, representing an - * \f$n\times n\f$ row-ordered symmetric matrix. - * @return permanent of the input matrix - */ -template -inline T permanent(std::vector &mat) -{ - int n = std::sqrt(static_cast(mat.size())); - llint x = static_cast(pow(2,n) - 1) ; - -#ifdef _OPENMP - int nthreads = omp_get_max_threads(); - omp_set_num_threads(nthreads); -#else - int nthreads = 1; -#endif - - std::vector tot(nthreads, static_cast(0)); - - std::vector threadbound_low(nthreads); - std::vector threadbound_hi(nthreads); - - for (int i=0; i < nthreads; i++) { - threadbound_low[i] = i*x/nthreads; - threadbound_hi[i] = (i+1)*x/nthreads; - } - threadbound_hi[nthreads-1] = x; - - #pragma omp parallel for shared(tot) - for (int ii = 0; ii < nthreads; ii++) { - T permtmp = static_cast(0); - int cntr = 0; - std::vector chitmp(n, 0); - std::vector tmp(n, static_cast(0)); - - for (llint k = threadbound_low[ii]; k < threadbound_hi[ii]; k++) { - T rowsumprod = static_cast(1); - llint kg2 = igray(k+1); - llint sgntmp = kg2 - igray(k); - llint sig = sgntmp/std::abs(sgntmp); - int pos = n - left_most_set_bit(sgntmp); - - if ( k == threadbound_low[ii] ) { - chitmp = dec2bin(kg2, n); - - // loop rows of the k-th submatrix - for ( int j = 0; j < n; j++) { - T localsum = static_cast(0); - for (int id = 0; id < n; id++) { - localsum += static_cast(chitmp[id]) * mat[id*n+j]; - } - tmp[j] += localsum; - - // update product of row sums - rowsumprod *= tmp[j]; - } - - cntr = static_cast(std::accumulate(chitmp.begin(), chitmp.end(), 0)); - } - else { - cntr += sig; - - for(int j = 0; j < n; j++ ) { - if (sig < 0) - tmp[j] -= mat[pos * n + j]; - else - tmp[j] += mat[pos * n + j]; - - rowsumprod *= tmp[j]; - } - } - - if ( (static_cast(n)-cntr) % 2 == 0) - permtmp += rowsumprod; - else - permtmp -= rowsumprod; - - } - tot[ii] = permtmp; - } - - return static_cast(std::accumulate(tot.begin(), tot.end(), static_cast(0))); -} - -/** - * Returns the permanent of a matrix using fsum. - * - * \rst - * - * Returns the permanent of a matrix using Ryser's algorithm - * with Gray code ordering. - * - * \endrst - * - * - * @tparam T type of the matrix data - * @param mat a flattened vector of size \f$n^2\f$, representing an - * \f$n\times n\f$ row-ordered symmetric matrix. - * @return permanent of the input matrix - */ -template -inline double perm_fsum(std::vector &mat) -{ - int n = std::sqrt(static_cast(mat.size())); - llint x = static_cast(pow(2,n) - 1) ; - -#ifdef _OPENMP - int nthreads = omp_get_max_threads(); - omp_set_num_threads(nthreads); -#else - int nthreads = 1; -#endif - - std::vector tot(nthreads, static_cast(0)); - - std::vector threadbound_low(nthreads); - std::vector threadbound_hi(nthreads); - - for (int i=0; i < nthreads; i++) { - threadbound_low[i] = i*x/nthreads; - threadbound_hi[i] = (i+1)*x/nthreads; - } - threadbound_hi[nthreads-1] = x; - - #pragma omp parallel for shared(tot) - for (int ii = 0; ii < nthreads; ii++) { - - fsum::sc_partials permtmp; // = 0; - - int cntr = 0; - std::vector chitmp(n, 0); - std::vector tmp(n, static_cast(0)); - - for (llint k = threadbound_low[ii]; k < threadbound_hi[ii]; k++) { - T rowsumprod = static_cast(1); - llint kg2 = igray(k+1); - llint sgntmp = kg2 - igray(k); - llint sig = sgntmp/std::abs(sgntmp); - int pos = 0; - - pos = n - left_most_set_bit(sgntmp); - - if ( k == threadbound_low[ii] ) { - chitmp = dec2bin(kg2, n); - - for ( int j = 0; j < n; j++) { - fsum::sc_partials localsum; //= static_cast(0); - for (int id = 0; id < n; id++) { - localsum += static_cast(chitmp[id]) * mat[id*n+j]; - } - tmp[j] += localsum; - rowsumprod *= tmp[j]; - } - - cntr = static_cast(std::accumulate(chitmp.begin(), chitmp.end(), 0)); - } - else { - cntr += sig; - - for(int j = 0; j < n; j++ ) { - if (sig < 0) - tmp[j] -= mat[pos * n + j]; - else - tmp[j] += mat[pos * n + j]; - - rowsumprod *= tmp[j]; - } - } - - if ( (static_cast(n)-cntr) % 2 == 0) - permtmp += rowsumprod; - else - permtmp -= rowsumprod; - - } - tot[ii] = permtmp; - } - return static_cast(std::accumulate(tot.begin(), tot.end(), static_cast(0))); -} - -/** - * \rst - * - * Returns the permanent of a matrix using the Ryser's algo - * with Gray code ordering - * - * \endrst - * - * - * This is a wrapper around the templated function `libwalrus::permanent` - * for Python integration. It accepts and returns complex double numeric types, - * and returns sensible values for empty and non-even matrices. - * - * In addition, this wrapper function automatically casts all matrices - * to type `complex`, allowing for greater precision than - * supported by Python and NumPy. - * - * - * @param mat vector representing the flattened matrix - * @return the permanent - */ -std::complex permanent_quad(std::vector> &mat) -{ - std::vector> matq(mat.begin(), mat.end()); - std::complex perm = permanent(matq); - return static_cast>(perm); -} - -/** - * \rst - * - * Returns the permanent of a matrix using Ryser's algo - * with Gray code ordering - * - * \endrst - * - * - * This is a wrapper around the templated function `libwalrus::permanent` - * for Python integration. It accepts and returns double numeric types, - * and returns sensible values for empty and non-even matrices. - * - * In addition, this wrapper function automatically casts all matrices - * to type `long double`, allowing for greater precision than supported - * by Python and NumPy. - * - * - * @param mat vector representing the flattened matrix - * @return the permanent - */ -double permanent_quad(std::vector &mat) -{ - std::vector matq(mat.begin(), mat.end()); - qp perm = permanent(matq); - return static_cast(perm); -} - -/** - * \rst - * - * Returns the permanent of a matrix using Ryser's algo - * with Gray code ordering with fsum - * - * \endrst - * - * - * This is a wrapper around the templated function `libwalrus::perm_fsum` - * for Python integration. It accepts and returns double numeric types, - * and returns sensible values for empty and non-even matrices. - * - * - * @param mat vector representing the flattened matrix - * @return the permanent - */ -double permanent_fsum(std::vector &mat) -{ - std::vector matq(mat.begin(), mat.end()); - double perm = perm_fsum(matq); - return static_cast(perm); -} - - -/** - * Returns the permanent of a matrix (nthreads=1) - * - * \rst - * - * Returns the permanent of a matrix using the BBFG algorithm. - * This algorithm was given by Glynn (2010) with the time-complexity - * of \f$O(n 2^n)\f$ using Gray code ordering. - * - * \endrst - * - * - * @tparam T type of the matrix data - * @param mat a flattened vector of size \f$n^2\f$, representing an - * \f$n\times n\f$ row-ordered symmetric matrix. - * @return permanent of the input matrix - */ -template -inline T perm_BBFG_serial(std::vector &mat) -{ - const size_t n = static_cast(std::sqrt(static_cast(mat.size()))); - const double p = pow(2, n-1); - const ullint x = static_cast(p); - if (p != x) { - std::cerr << "overflow to inf" << std::endl; - exit(EXIT_FAILURE); - } - - constexpr T p1 = static_cast(1.0); - constexpr T p2 = static_cast(2.0); - constexpr T n2 = static_cast(-2.0); - constexpr T zero = static_cast(0); - - size_t i, j; - std::vector colsum(n, zero); - // init colsum - for (i=0; i < n; ++i) { - for (j=0; j < n; ++j) { - colsum[i] += mat[j*n + i]; - } - } - - T mulcolsum, coeff; - T total = zero; - ullint k, og=0, ng, gd; - int sgn=1, gdi; - - // over all 2^{n-1} permutations of delta - for (k=1; k < x+1; ++k) { - mulcolsum = std::accumulate( - colsum.begin(), - colsum.end(), - p1, - std::multiplies()); - total += sgn > 0 ? mulcolsum : -mulcolsum; - - // updating gray order - ng = igray(k); - gd = og ^ ng; - - coeff = og > ng ? p2 : n2; - gdi = left_most_set_bit(gd); - for (j=0; j < n; ++j) { - colsum[j] += coeff * mat[gdi * n + j]; - } - - sgn = -sgn; - og = ng; - } - - return total / static_cast(x); -} - -/** - * Returns the permanent of a matrix (nthreads=1) (2nd version) - * - * \rst - * - * Returns the permanent of a matrix using the BBFG algorithm. - * This algorithm was given by Glynn (2010) with the time-complexity - * of \f$O(n 2^n)\f$ using Gray code ordering. - * - * \endrst - * - * - * @tparam T type of the matrix data - * @param mat a flattened vector of size \f$n^2\f$, representing an - * \f$n\times n\f$ row-ordered symmetric matrix. - * @return permanent of the input matrix - */ -template -inline T perm_BBFG_serial2(std::vector &mat) -{ - const size_t n = static_cast(std::sqrt(static_cast(mat.size()))); - const long double p = pow(2, n-1); - const T x = static_cast(p); - if (p == HUGE_VAL || p != x) { - std::cerr << "overflow to inf" << std::endl; - exit(EXIT_FAILURE); - } - - constexpr T p2 = static_cast(2.0); - constexpr T p1 = static_cast(1.0); - constexpr T zero = static_cast(0); - - std::vector grays(n); - std::iota(grays.begin(), grays.end(), 0); - std::vector coeffs(n, p2); - std::vector colsum(n, zero); - T mulcolsum, total = p1; - size_t i, j, k=0; - int sgn=1; - - // init colsum - for (i=0; i < n; ++i) { - for (j=0; j < n; ++j) { - colsum[i] += mat[j*n + i]; - } - total *= colsum[i]; - } - - while (k < n-1) { - for (j=0; j < n; ++j) { - colsum[j] -= coeffs[k] * mat[k*n+j]; - } - - mulcolsum = std::accumulate( - colsum.begin(), - colsum.end(), - p1, - std::multiplies()); - - coeffs[k] = -coeffs[k]; - sgn = -sgn; - total += sgn > 0 ? mulcolsum : -mulcolsum; - k = next_perm_ordering_index(grays, k); - } - - return total / x; -} - -/** - * Returns the permanent of a matrix - * - * \rst - * - * Returns the permanent of a matrix using the BBFG algorithm. - * This algorithm was given by Glynn (2010) with the time-complexity - * of \f$O(n 2^n)\f$ using Gray code ordering. - * - * \endrst - * - * This is a wrapper function for computing permanent of a matrix - * based on Balasubramanian-Bax-Franklin-Glynn (BBFG) formula. - * - * - * @tparam T type of the matrix data - * @param mat a flattened vector of size \f$n^2\f$, representing an - * \f$n\times n\f$ row-ordered symmetric matrix. - * @return permanent of the input matrix - */ -template -inline T perm_BBFG(std::vector &mat) { - return perm_BBFG_serial2(mat); -} - -/** - * \rst - * - * Returns the permanent of a matrix using the BBFG formula. - * This algorithm was given by Glynn (2010) with the time-complexity - * of \f$O(n 2^n)\f$ using Gray code ordering. - * - * \endrst - * - * This is a wrapper around the templated function `libwalrus::perm_BBFG` - * for Python integration. It accepts and returns double numeric types, - * and returns sensible values for empty and non-even matrices. - * - * In addition, this wrapper function automatically casts all matrices - * to type `long double`, allowing for greater precision than supported - * by Python and NumPy. - * - * - * @param mat vector representing the flattened matrix - * @return the permanent - */ -double perm_BBFG_qp(std::vector &mat) -{ - std::vector matqp(mat.begin(), mat.end()); - qp perm = perm_BBFG(matqp); - return static_cast(perm); -} - -/** - * \rst - * - * Returns the permanent of a matrix using the BBFG formula. - * This algorithm was given by Glynn (2010) with the time-complexity - * of \f$O(n 2^n)\f$ using Gray code ordering. - * - * \endrst - * - * This is a wrapper around the templated function `libwalrus::perm_BBFG` - * for Python integration. It accepts and returns complex double numeric types, - * and returns sensible values for empty and non-even matrices. - * - * In addition, this wrapper function automatically casts all matrices - * to type `complex`, allowing for greater precision than - * supported by Python and NumPy. - * - * - * @param mat vector representing the flattened matrix - * @return the permanent - */ -std::complex perm_BBFG_cmplx(std::vector> &mat) -{ - std::vector> matx(mat.begin(), mat.end()); - std::complex perm = perm_BBFG(matx); - return static_cast>(perm); -} - -} diff --git a/tests/libwalrus_unittest.cpp b/tests/libwalrus_unittest.cpp index d64859343..cbe8f2309 100644 --- a/tests/libwalrus_unittest.cpp +++ b/tests/libwalrus_unittest.cpp @@ -20,156 +20,6 @@ const double tol = 1.0e-10f; const double tol2 = 1.0e-7f; -namespace permanent { - -TEST(PermanentRealFsum, CompleteGraph) { - std::vector mat2(4, 1.0); - std::vector mat3(9, 1.0); - std::vector mat4(16, 1.0); - - EXPECT_NEAR(2, libwalrus::permanent_fsum(mat2), tol); - EXPECT_NEAR(6, libwalrus::permanent_fsum(mat3), tol); - EXPECT_NEAR(24, libwalrus::permanent_fsum(mat4), tol); -} - -TEST(PermanentFsum, Random) { - std::vector mat(9, 1.0); - - std::default_random_engine generator; - generator.seed(20); - std::normal_distribution distribution(0.0, 1.0); - - for (int i = 0; i < 3; i++) { - for (int j = 0; j < 3; j++) { - double randnum = distribution(generator); - mat[i * 3 + j] = randnum; - } - } - - double expected = mat[2] * mat[4] * mat[6] + mat[1] * mat[5] * mat[6] + - mat[2] * mat[3] * mat[7] + mat[0] * mat[5] * mat[7] + - mat[1] * mat[3] * mat[8] + mat[0] * mat[4] * mat[8]; - - EXPECT_NEAR(expected, libwalrus::permanent_fsum(mat), tol); -} - -TEST(PermanentReal, CompleteGraph) { - std::vector mat2(4, 1.0); - std::vector mat3(9, 1.0); - std::vector mat4(16, 1.0); - - EXPECT_NEAR(2, libwalrus::permanent_quad(mat2), tol); - EXPECT_NEAR(6, libwalrus::permanent_quad(mat3), tol); - EXPECT_NEAR(24, libwalrus::permanent_quad(mat4), tol); -} - -TEST(PermanentReal, Random) { - std::vector mat(9, 1.0); - - std::default_random_engine generator; - generator.seed(20); - std::normal_distribution distribution(0.0, 1.0); - - for (int i = 0; i < 3; i++) { - for (int j = 0; j < 3; j++) { - double randnum = distribution(generator); - mat[i * 3 + j] = randnum; - } - } - - double expected = mat[2] * mat[4] * mat[6] + mat[1] * mat[5] * mat[6] + - mat[2] * mat[3] * mat[7] + mat[0] * mat[5] * mat[7] + - mat[1] * mat[3] * mat[8] + mat[0] * mat[4] * mat[8]; - - EXPECT_NEAR(expected, libwalrus::permanent_quad(mat), tol); -} - -TEST(PermanentComplex, Random) { - std::vector> mat(9, 1.0); - - std::default_random_engine generator; - generator.seed(20); - std::normal_distribution distribution(0.0, 1.0); - - for (int i = 0; i < 3; i++) { - for (int j = 0; j < 3; j++) { - double randnum1 = distribution(generator); - double randnum2 = distribution(generator); - mat[i * 3 + j] = std::complex(randnum1, randnum2); - } - } - - std::complex expected = - mat[2] * mat[4] * mat[6] + mat[1] * mat[5] * mat[6] + - mat[2] * mat[3] * mat[7] + mat[0] * mat[5] * mat[7] + - mat[1] * mat[3] * mat[8] + mat[0] * mat[4] * mat[8]; - - std::complex perm = libwalrus::permanent_quad(mat); - - EXPECT_NEAR(std::real(expected), std::real(perm), tol); - EXPECT_NEAR(std::imag(expected), std::imag(perm), tol); -} - -TEST(PermanentRealBBFG, CompleteGraph) { - std::vector mat2(4, 1.0); - std::vector mat3(9, 1.0); - std::vector mat4(16, 1.0); - - EXPECT_NEAR(2, libwalrus::perm_BBFG_qp(mat2), tol); - EXPECT_NEAR(6, libwalrus::perm_BBFG_qp(mat3), tol); - EXPECT_NEAR(24, libwalrus::perm_BBFG_qp(mat4), tol); -} - -TEST(PermanentRealBBFG, Random) { - std::vector mat(9, 1.0); - - std::default_random_engine generator; - generator.seed(20); - std::normal_distribution distribution(0.0, 1.0); - - for (int i = 0; i < 3; i++) { - for (int j = 0; j < 3; j++) { - double randnum = distribution(generator); - mat[i * 3 + j] = randnum; - } - } - - double expected = mat[2] * mat[4] * mat[6] + mat[1] * mat[5] * mat[6] + - mat[2] * mat[3] * mat[7] + mat[0] * mat[5] * mat[7] + - mat[1] * mat[3] * mat[8] + mat[0] * mat[4] * mat[8]; - - EXPECT_NEAR(expected, libwalrus::perm_BBFG_qp(mat), tol); -} - - -TEST(PermanentComplexBBFG, Random) { - std::vector> mat(9, 1.0); - - std::default_random_engine generator; - generator.seed(20); - std::normal_distribution distribution(0.0, 1.0); - - for (int i = 0; i < 3; i++) { - for (int j = 0; j < 3; j++) { - double randnum1 = distribution(generator); - double randnum2 = distribution(generator); - mat[i * 3 + j] = std::complex(randnum1, randnum2); - } - } - - std::complex expected = - mat[2] * mat[4] * mat[6] + mat[1] * mat[5] * mat[6] + - mat[2] * mat[3] * mat[7] + mat[0] * mat[5] * mat[7] + - mat[1] * mat[3] * mat[8] + mat[0] * mat[4] * mat[8]; - - std::complex perm = libwalrus::perm_BBFG_cmplx(mat); - - EXPECT_NEAR(std::real(expected), std::real(perm), tol); - EXPECT_NEAR(std::imag(expected), std::imag(perm), tol); -} - -} // namespace permanent - namespace recursive_real { // Unit tests for the real recursive_hafnian function diff --git a/thewalrus/__init__.py b/thewalrus/__init__.py index 0a9491ef1..c4aa14722 100644 --- a/thewalrus/__init__.py +++ b/thewalrus/__init__.py @@ -124,14 +124,8 @@ interferometer, grad_hermite_multidimensional, ) -from ._permanent import ( - perm, - perm_complex, - perm_real, - permanent_repeated, - perm_BBFG_real, - perm_BBFG_complex, -) +from ._permanent import perm, permanent_repeated + from ._torontonian import ( tor, threshold_detection_prob_displacement, diff --git a/thewalrus/_permanent.py b/thewalrus/_permanent.py index 1b341f1c2..aa42e2031 100644 --- a/thewalrus/_permanent.py +++ b/thewalrus/_permanent.py @@ -11,33 +11,49 @@ # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. # See the License for the specific language governing permissions and # limitations under the License. -""" -Permanent Python interface +r""" +Permanent Algorithms +============================= + +.. currentmodule:: thewalrus._permanent + +This submodule provides access to tools for finding the permanent of a matrix. The algorithms implemented +here was first derived in +* Ryser, Herbert John (1963). + Combinatorial Mathematics, The Carus Mathematical Monographs, Vol. 14, Mathematical Association of America. +* Glynn, David G. + (2010), "The permanent of a square matrix", European Journal of Combinatorics, 31 (7): 1887–1891. + `_ + +Summary +------- +.. autosummary:: + + perm + perm_ryser + perm_bbfg + permanent_repeated + +Code details +------------ """ import numpy as np - +from numba import jit from ._hafnian import hafnian_repeated -from .libwalrus import perm_complex, perm_real, perm_BBFG_real, perm_BBFG_complex -def perm(A, quad=True, fsum=False, method="bbfg"): - """Returns the permanent of a matrix using the `method` formula +def perm(A, method="bbfg"): + """Returns the permanent of a matrix using various methods. - For more direct control, you may wish to call :func:`perm_real` - or :func:`perm_complex` directly. Args: - A (array): a square array. - quad (bool): If ``True``, the input matrix is cast to a ``long double`` - matrix internally for a quadruple precision hafnian computation. - fsum (bool): Whether to use the ``fsum`` method for higher accuracy summation. - Note that if ``fsum`` is true, double precision will be used, and the - ``quad`` keyword argument will be ignored. - method (string): "ryser" calls the associated methods to - `Ryser formula `_, - and "bbfg" calls the associated methods to - `BBFG formula `_. - + A (array[float or complex]): a square array. + method (string): Set this to ``"ryser"`` to use the + `Ryser formula + `_, + or ``"bbfg"`` to use the + `BBFG formula + `_. Returns: np.float64 or np.complex128: the permanent of matrix A. """ @@ -68,15 +84,88 @@ def perm(A, quad=True, fsum=False, method="bbfg"): isRyser = bool(method != "bbfg") - if A.dtype == np.complex: - if np.any(np.iscomplex(A)): - return perm_complex(A, quad=quad) if isRyser else perm_BBFG_complex(A) - return ( - perm_real(np.float64(A.real), quad=quad, fsum=fsum) - if isRyser - else perm_BBFG_real(np.float64(A.real)) - ) - return perm_real(A, quad=quad, fsum=fsum) if isRyser else perm_BBFG_real(A) + return perm_ryser(A) if isRyser else perm_bbfg(A) + + +@jit(nopython=True) +def perm_ryser(M): # pragma: no cover + """ + Returns the permanent of a matrix using the Ryser formula in Gray ordering. + + The code is an re-implementation from a Python 2 code found in + `Permanent code golf + `_ + using Numba. + + Args: + M (array) : a square array. + + Returns: + np.float64 or np.complex128: the permanent of matrix M. + """ + n = len(M) + # row_comb keeps the sum of previous subsets. + # Every iteration, it removes a term and/or adds a new term + # to give the term to add for the next subset + row_comb = np.zeros((n), dtype=M.dtype) + total = 0 + old_grey = 0 + sign = +1 + binary_power_dict = [2 ** i for i in range(n)] + num_loops = 2 ** n + for k in range(0, num_loops): + bin_index = (k + 1) % num_loops + reduced = np.prod(row_comb) + total += sign * reduced + new_grey = bin_index ^ (bin_index // 2) + grey_diff = old_grey ^ new_grey + grey_diff_index = binary_power_dict.index(grey_diff) + new_vector = M[grey_diff_index] + direction = (old_grey > new_grey) - (old_grey < new_grey) + for i in range(n): + row_comb[i] += new_vector[i] * direction + sign = -sign + old_grey = new_grey + return total + + +@jit(nopython=True) +def perm_bbfg(M): # pragma: no cover + """ + Returns the permanent of a matrix using the bbfg formula in Gray ordering + + The code is a re-implementation from a Python 2 code found in + `Permanent code golf + `_ + using Numba. + + Args: + M (array) : a square array. + + Returns: + np.float64 or np.complex128: the permanent of a matrix M. + """ + + n = len(M) + row_comb = np.sum(M, 0) + total = 0 + old_gray = 0 + sign = +1 + binary_power_dict = [2 ** i for i in range(n)] + num_loops = 2 ** (n - 1) + for bin_index in range(1, num_loops + 1): + reduced = np.prod(row_comb) + total += sign * reduced + new_gray = bin_index ^ (bin_index // 2) + gray_diff = old_gray ^ new_gray + gray_diff_index = binary_power_dict.index(gray_diff) + new_vector = M[gray_diff_index] + direction = 2 * ((old_gray > new_gray) - (old_gray < new_gray)) + for i in range(n): + row_comb[i] += new_vector[i] * direction + sign = -sign + old_gray = new_gray + return total / num_loops def permanent_repeated(A, rpt): diff --git a/thewalrus/libwalrus.pyx b/thewalrus/libwalrus.pyx index d8321ea16..e0070308d 100644 --- a/thewalrus/libwalrus.pyx +++ b/thewalrus/libwalrus.pyx @@ -27,18 +27,10 @@ cdef extern from "../include/libwalrus.hpp" namespace "libwalrus": T hafnian[T](vector[T] &mat) T hafnian_recursive[T](vector[T] &mat) T loop_hafnian[T](vector[T] &mat) - T permanent[T](vector[T] &mat) - T perm_BBFG[T](vector[T] &mat) T hafnian_rpt[T](vector[T] &mat, vector[int] &nud) T loop_hafnian_rpt[T](vector[T] &mat, vector[T] &mu, vector[int] &nud) - double permanent_quad(vector[double] &mat) - double complex permanent_quad(vector[double complex] &mat) - double perm_fsum[T](vector[T] &mat) - double permanent_fsum(vector[double] &mat) - double perm_BBFG_qp(vector[double] &mat) - double complex perm_BBFG_cmplx(vector[double complex] &mat) double hafnian_recursive_quad(vector[double] &mat) double complex hafnian_recursive_quad(vector[double complex] &mat) @@ -241,98 +233,3 @@ def haf_real(double[:, :] A, bint loop=False, bint recursive=True, quad=True, bi -# ============================================================================== -# Permanent - -def perm_complex(double complex[:, :] A, quad=True): - """Returns the hafnian of a complex matrix A via the C++ libwalrus library. - - Args: - A (array): a np.float, square array - quad (bool): If ``True``, the input matrix is cast to a ``long double complex`` - matrix internally for a quadruple precision hafnian computation. - - Returns: - np.complex128: the hafnian of matrix A - """ - cdef int i, j, n = A.shape[0] - cdef vector[double complex] mat - - for i in range(n): - for j in range(n): - mat.push_back(A[i, j]) - - # Exposes a c function to python - if quad: - return permanent_quad(mat) - - return permanent(mat) - -def perm_real(double [:, :] A, quad=True, fsum=False): - """Returns the hafnian of a real matrix A via the C++ libwalrus library. - - Args: - A (array): a np.float64, square array - quad (bool): If ``True``, the input matrix is cast to a ``long double`` - matrix internally for a quadruple precision hafnian computation. - fsum (bool): If ``True``, ``fsum`` method is used for summation. - - - Returns: - np.float64: the hafnian of matrix A - """ - cdef int i, j, n = A.shape[0] - cdef vector[double] mat - - for i in range(n): - for j in range(n): - mat.push_back(A[i, j]) - - if fsum: - return permanent_fsum(mat) - - # Exposes a c function to python - if quad: - return permanent_quad(mat) - - return permanent(mat) - -def perm_BBFG_complex(double complex[:, :] A): - """Returns the hafnian of a complex matrix A via the C++ libwalrus library - using Balasubramanian-Bax-Franklin-Glynn formula. - - Args: - A (array): a np.float, square array - - Returns: - np.complex128: the hafnian of matrix A - """ - cdef int i, j, n = A.shape[0] - cdef vector[double complex] mat - - for i in range(n): - for j in range(n): - mat.push_back(A[i, j]) - - # Exposes a c function to python - return perm_BBFG_cmplx(mat) - -def perm_BBFG_real(double [:, :] A): - """Returns the hafnian of a real matrix A via the C++ libwalrus library - using Balasubramanian-Bax-Franklin-Glynn formula. - - Args: - A (array): a np.float64, square array - - Returns: - np.float64: the hafnian of matrix A - """ - cdef int i, j, n = A.shape[0] - cdef vector[double] mat - - for i in range(n): - for j in range(n): - mat.push_back(A[i, j]) - - # Exposes a c function to python - return perm_BBFG_qp(mat) diff --git a/thewalrus/tests/test_permanent.py b/thewalrus/tests/test_permanent.py index c553868a8..cc2fee075 100644 --- a/thewalrus/tests/test_permanent.py +++ b/thewalrus/tests/test_permanent.py @@ -18,14 +18,12 @@ import numpy as np from scipy.special import factorial as fac -from thewalrus import ( - perm, - perm_real, - perm_complex, - permanent_repeated, - perm_BBFG_real, - perm_BBFG_complex, -) +from thewalrus import perm, permanent_repeated + +perm_real = perm +perm_complex = perm +perm_BBFG_real = lambda x: perm(x, method="bbfg") +perm_BBFG_complex = lambda x: perm(x, method="bbfg") class TestPermanentWrapper: