Implement C++ part of 2-cut optimization

Makefile

 default: groupConnectivity.so clean_obj
 
-groupConnectivity.so: groupConnectivity.pyx group-connectivity.h setup.py compileTimeOptions.h generateCompileTimeOptions.sh rings.h fast-array.h parmap.py
+groupConnectivity.so: groupConnectivity.pyx group-connectivity.h setup.py compileTimeOptions.h generateCompileTimeOptions.sh rings.h fast-array.h twoCuts.h parmap.py
 	./generateCompileTimeOptions.sh > options.h
 	python setup.py build_ext
 	cp build/lib*/groupConnectivity.so .

group-connectivity.h

 #ifndef __GROUP_CONNECTIVITY_H__
 #define __GROUP_CONNECTIVITY_H__
 
+#include <stdlib.h>
+#include <stdint.h>
 #include <vector>
 #include "options.h"
 #include "rings.h"
 #include "fast-array.h"
 
-
 typedef int EdgeId;
 typedef int DirectedEdgeId;
 typedef std::pair<DirectedEdgeId, DirectedEdgeId> TwoCut;
@@ -39,6 +40,8 @@ struct Tester : public AbstractTester {
   typedef typename Ring::T T;
   typedef ::Mapping<Ring> Mapping;
 
+# include "twoCuts.h"
+
   size_t numForb;
   int edges;
 
@@ -47,11 +50,14 @@ struct Tester : public AbstractTester {
 # if EXPLICIT_NORMAL_EDGES
   std::vector< EdgeId > normalEdges;
 # endif
+# if USE_TWO_CUTS
+  std::vector< TwoCutInt > twoCuts;
+# endif
 
   virtual void init(
     int edges_,
     std::vector<EdgeId> spanningTree,
-    std::vector<TwoCut> twoCuts,
+    std::vector<TwoCut> twoCuts_,
     std::vector< std::vector<DirectedEdgeId> > elementaryCycles
   ) {
     edges = edges_;
@@ -64,6 +70,10 @@ struct Tester : public AbstractTester {
     ClassesType C(num_classes, 0);
     C.swap(classes);
 
+#   if EXPLICIT_NORMAL_EDGES
+    std::vector<bool> edgeList(edges, 1);
+#   endif
+
     for (auto e : spanningTree)
       classEdges.push_back(e - 1);
 
@@ -79,8 +89,25 @@ struct Tester : public AbstractTester {
       }
     }
 
+#   if USE_TWO_CUTS
+    for (auto cut : twoCuts_) {
+      if (cut.first < 0) {
+        cut.first *= -1;
+        cut.second *= -1;
+      }
+
+      TwoCutInt intCut;
+      intCut.a = cut.first - 1;
+      intCut.flip = (cut.second < 0);
+      intCut.b = (intCut.flip ? - cut.second : cut.second) - 1;
+
+      twoCuts.push_back(intCut);
+      edgeList[intCut.a] = false;
+      edgeList[intCut.b] = false;
+    }
+#   endif
+
 #   if EXPLICIT_NORMAL_EDGES
-    std::vector<bool> edgeList(edges, 1);
     for (size_t i = 0; i < edges; i++)
       if (edgeList[i]) normalEdges.push_back(i);
 #   endif
@@ -103,6 +130,18 @@ struct Tester : public AbstractTester {
       }
     }
 
+#   if USE_TWO_CUTS
+    for (const auto& cut : twoCuts) {
+      int val = cut.getValue(forb) + 1;
+      if (val >= cut.max_value) {
+        cut.encode(forb, 0);
+      } else {
+        cut.encode(forb, val);
+        return true;
+      }
+    }
+#   endif
+
     return false;
   }
 #else // !USE_NEXT_FORB
@@ -127,6 +166,17 @@ struct Tester : public AbstractTester {
         map.combine(map[i.first], i.second);
     }
 
+#   if USE_TWO_CUTS
+    for (const auto &c : twoCuts) {
+      T m_a = map[c.a];
+      T m_b = map[c.b];
+      if (map[c.a] > map[c.b]) {
+        map.assign(c.a, m_b);
+        map.assign(c.b, m_a);
+      }
+    }
+#   endif
+
     return map;
   }
 
@@ -138,6 +188,13 @@ struct Tester : public AbstractTester {
       index += map[i];
     }
 
+#   if USE_TWO_CUTS
+    for (const auto& c : twoCuts) {
+      index *= c.max_value;
+      index += c.getValue(map);
+    }
+#   endif
+
     return index;
   }
 
@@ -184,5 +241,10 @@ struct Tester : public AbstractTester {
   virtual ~Tester() {}
 };
 
+#if USE_TWO_CUTS
+template < typename R >
+const typename Tester<R>::TwoCutInt::Tables Tester<R>::TwoCutInt::S;
+#endif
+
 #endif
 

twoCuts.h

+// Internals of Tester class
+#if USE_TWO_CUTS
+
+struct TwoCutInt {
+  EdgeId a, b;
+  bool flip;
+
+  static const int max_value = ((Ring::size - 1) * (Ring::size - 2)) / 2;
+
+  struct Tables {
+    char decodeMatrix[Ring::size][Ring::size];
+    struct { char a, b; } encodeTable[max_value];
+
+    Tables() {
+      encodeTable[0] = {1, 2};
+
+      for (int i = 1; i < max_value; i++) {
+        bool carry = (encodeTable[i-1].b >= Ring::size);
+        encodeTable[i].a = encodeTable[i-1].a + carry;
+        encodeTable[i].b = carry ? encodeTable[i].a + 1 : encodeTable[i-1].b + 1;
+      }
+
+      memset(decodeMatrix, -1, sizeof(decodeMatrix));
+
+      for (int i = 0; i < max_value; i++)
+        decodeMatrix[encodeTable[i].a][encodeTable[i].b] = i;
+    }
+  };
+
+  static const Tables S;
+
+  inline int getValue(const Mapping &map) const {
+    return S.decodeMatrix[map[a]][flip ? Ring::negate(map[b]) : map[b]];
+  }
+
+  inline void encode(Mapping& map, int value) const {
+    auto vals = S.encodeTable[value];
+    map.assign(a, vals.a);
+    map.assign(b, flip ? Ring::negate(vals.b) : vals.b);
+  }
+};
+
+#endif
+