Graphs

makefile

@@ -7,7 +7,7 @@ setup = setup.py
 default : $(target) test
 
 # Compile
-$(target): $(srcdir)/permanent.c
+$(target): $(srcdir)/* 
 	python $(setup) build_ext --inplace 
 
 test : 

run-tests.py

@@ -4,6 +4,7 @@ import multiprocessing as mp
 import numpy as np
 import lib
 import time
+from matplotlib import pyplot as plt
 
 def perm_ryser(a):
     ''' the permanent calculated using the ryser formula. much faster than the naive approach '''
@@ -16,23 +17,30 @@ def perm_ryser(a):
     terms=map(get_term, indeces) 
     return np.sum(terms)*((-1)**n)
 
-dimension=9
-real=np.random.uniform(-1, 1, dimension*dimension).reshape((dimension, dimension))
-imag=np.random.uniform(-1, 1, dimension*dimension).reshape((dimension, dimension))
-submatrix=real+1j*imag
+maxtime=1
+dimensions=range(1,11)
 
-#print lib.permanent(submatrix), perm_ryser(submatrix)
+for (function, label) in zip((lib.permanent, perm_ryser), ("C", "Python")):
+    counts=[]
+    for dimension in dimensions:
+        print dimension
+        real=np.random.uniform(-1, 1, dimension*dimension).reshape((dimension, dimension))
+        imag=np.random.uniform(-1, 1, dimension*dimension).reshape((dimension, dimension))
+        submatrix=real+1j*imag
 
-#sys.exit(0)
-#t=time.clock()
-#for i in range(1000):
-    #perm_ryser(submatrix)
-#t1=time.clock()-t
+        t=time.clock()
+        n=0
+        while time.clock()-t < maxtime:
+            for i in range(5):
+                function(submatrix)
+            n+=5
+        counts.append(n)
 
-t=time.clock()
-for i in range(1000):
-    lib.permanent(submatrix)
-t2=time.clock()-t
-print t2
+    plt.plot(dimensions, counts, '.-', label=label)
 
-#print t1/t2
+plt.ylabel('Number of permanents per second')
+plt.xlabel('Dimension')
+plt.xlim(min(dimensions), max(dimensions))
+plt.legend()
+plt.semilogy()
+plt.savefig('out.pdf')

src/bithacks.h

@@ -8,7 +8,7 @@ int countbits(unsigned int n)
     q = (q & 0x0F0F0F0F0F0F0F0F) + ((q & 0xF0F0F0F0F0F0F0F0) >> 4);
     q = (q & 0x00FF00FF00FF00FF) + ((q & 0xFF00FF00FF00FF00) >> 8);
     q = (q & 0x0000FFFF0000FFFF) + ((q & 0xFFFF0000FFFF0000) >> 16);
-    q = (q & 0x00000000FFFFFFFF) + ((q & 0xFFFFFFFF00000000) >> 32); // This last & isq't strictly qecessary.
+    q = (q & 0x00000000FFFFFFFF) + ((q & 0xFFFFFFFF00000000) >> 32); // This last & isn't strictly qecessary.
     return q;
 }
 

src/npy_util.h

@@ -40,8 +40,8 @@ void complex_inc(npy_complex128 *a, npy_complex128 b) {
 
 // Multipy a number by another one
 void complex_multiply(npy_complex128 *a, npy_complex128 b) { 
-    npy_complex128 c = {.real=a->real, .imag=a->imag};
-    a->real  = c.real*b.real-c.imag*b.imag;
-    a->imag  = c.real*b.imag+c.imag*b.real;
+    double c = a->real;
+    a->real  = a->real*b.real-a->imag*b.imag;
+    a->imag  = c*b.imag+a->imag*b.real;
 }
 

src/permanent.c

@@ -4,7 +4,7 @@
 #include <numpy/arrayobject.h>
 #include "npy_util.h"
 #include "bithacks.h"
-#include "permanents.h"
+#include "ryser.h"
 
 // Forward function declaration 
 static PyObject *permanent(PyObject *self, PyObject *args);    
@@ -28,6 +28,6 @@ static PyObject *permanent(PyObject *self, PyObject *args) {
   if (!PyArg_ParseTuple(args, "O!", &PyArray_Type, &submatrix)) {return NULL;}
 
   // Compute the permanent
-  npy_complex128 p = perm_ryser(submatrix);
+  npy_complex128 p = ryser(submatrix);
   return PyComplex_FromDoubles(p.real, p.imag);
 }

src/permanents.h

@@ -1,22 +0,0 @@
-// Ryser's algorithm 
-static npy_complex128 perm_ryser(PyArrayObject *submatrix) {
-    int n = (int) PyArray_DIM(submatrix, 0);
-    npy_complex128 rowsum, rowsumprod;
-    npy_complex128 perm = complex_zero;
-    int exp = 1 << n; 
-    int i, y, z;
-    for (i=0; i<exp; ++i) {
-        rowsumprod = complex_one;
-        for (y=0; y<n; ++y) {               
-            rowsum = complex_zero;
-            for (z=0; z<n; ++z) { 
-                if ((i & (1 << z)) != 0) { complex_inc(&rowsum, SM(z, y)); }
-            }
-            complex_multiply(&rowsumprod, rowsum);
-        }
-        complex_inc(&perm, complex_float_prod(rowsumprod, bitparity(i)));
-    }
-    if (n%2 == 1) {perm=complex_float_prod(perm, -1);}
-    return perm;
-}
-