improve test_cycle.py

test_cycle.py

 from graph_tools.base import *
-from graph_tools.parameters import CircuitDoubleCover, VertexCount
+from graph_tools.parameters import CircuitDoubleCover, CycleDoubleCover, VertexCount
 from graph_tools.misc import sun
 
 from sage.all import MixedIntegerLinearProgram
 from itertools import permutations
 
-def solve_gadget(g, alt_gadgets):
+verbose = False
+
+def solve_gadget(g, alt_gadgets, param=CycleDoubleCover(5)):
   N = g.size()
 
   P = MixedIntegerLinearProgram(maximization=False)
   V = P.new_variable(real=True, nonnegative=True)
+  P_dual = MixedIntegerLinearProgram(maximization=True)
+  V_dual = P_dual.new_variable(real=True, nonnegative=True)
 
-  BOUNDARIES = list(CircuitDoubleCover.enumerate_boundaries(N))
+  BOUNDARIES = list(param.enumerate_boundaries(N))
 
   def fake_gadget(b):
     return FakeGadget(N, [ Boundary(b, 1) ])
@@ -19,38 +23,54 @@ def solve_gadget(g, alt_gadgets):
   def join(a, b):
     return Gadget.join([ a, b ], [ ((1, i), (2, i)) for i in range(1,N+1) ], [])
 
-  def gadget_to_poly(gadget):
-    return sum(
-      V[b] * float(join(fake_gadget(b), gadget).eval(CircuitDoubleCover))
-      for b in BOUNDARIES
-    )
+  def gadget_to_multiplicity_vector(gadget):
+    return [
+      join(fake_gadget(b), gadget).eval(param) for b in BOUNDARIES
+    ]
+
+  dual_objective = 0
+  dual_constraints = [0] * len(BOUNDARIES)
+  def gadget_to_constraint(gadget, var):
+    nonlocal dual_objective
+
+    mv = gadget_to_multiplicity_vector(gadget)
+    n = VertexCount.finalize(gadget.eval_gadget(VertexCount))
 
-  def gadget_to_constraint(gadget):
-    poly = gadget_to_poly(gadget)
-    n = VertexCount.finalize(gadget.eval(VertexCount))
+    P.add_constraint(sum( a * V[b] for a, b in zip(mv, BOUNDARIES) ) >= 2**(n // 2))
+    for i, a in enumerate(mv):
+      dual_constraints[i] += a * V_dual[var]
+    dual_objective += 2**(n // 2) * V_dual[var]
 
-    P.add_constraint(poly >= 2**(n // 2))
+  for i, ag in enumerate(alt_gadgets): gadget_to_constraint(ag, i)
+  P_dual.set_objective(dual_objective)
 
-  for ag in alt_gadgets: gadget_to_constraint(ag)
+  mv = gadget_to_multiplicity_vector(g)
+  P.set_objective(sum( a * V[b] for a, b in zip(mv, BOUNDARIES) ))
 
-  P.set_objective(gadget_to_poly(g))
+  for b, const in zip(mv, dual_constraints):
+    P_dual.add_constraint(const <= b)
 
-  if False:
-    P.show()
-    for i, v in sorted(P.get_values(V).items()):
+  ret = (
+    P.solve(), P_dual.solve(),
+    2**(VertexCount.finalize(g.eval_gadget(VertexCount))//2)
+  )
+
+  if verbose:
+    #P_dual.show()
+    for i, v in sorted(P_dual.get_values(V_dual).items()):
       print('w_%s = %s' % (i, v))
 
-  return (floor(P.solve()), 2**(g.eval(UnderlayingGraph)[0].value.num_verts()//2))
+  return ret
 
 
-def friend_boundary(g):
+def friend_boundary(g, param):
   N = g.size();
-  BOUNDARIES = list(CircuitDoubleCover.enumerate_boundaries(N))
+  BOUNDARIES = list(param.enumerate_boundaries(N))
 
   def test_boundary(b):
     fg = FakeGadget(N, [ Boundary(b, 1) ])
     joins = [ ((1, i), (2, i)) for i in range(1,N+1) ]
-    return (b, Gadget.join([ fg, g ], joins, []).eval(CircuitDoubleCover))
+    return (b, Gadget.join([ fg, g ], joins, []).eval(param))
 
   return [ test_boundary(b) for b in BOUNDARIES ]
 
@@ -67,34 +87,57 @@ G4 = [
 
 G4f = [
   Gadget.join([FREE_EDGE]*2, [], [ (1,1), (1,2), (2,1), (2,2) ]),
-  Gadget.join([FREE_EDGE]*2, [], [ (1,1), (2,1), (1,2), (2,2) ]),
   Gadget.join([FREE_EDGE]*2, [], [ (1,1), (2,1), (2,2), (1,2) ]),
+  Gadget.join([FREE_EDGE]*2, [], [ (1,1), (2,1), (1,2), (2,2) ]),
 ]
 
-G5 = [
-  Gadget.join([CUBIC_VERTEX] * 3, [((1,1), (2,1)), ((2,2), (3,2))], out)
+G5fp = [
+  Gadget.join([ CUBIC_VERTEX, FREE_EDGE ], [], rot([ (1,1), (1,2), (1,3), (2,1), (2,2) ], i))
+  for i in range(5)
+]
+
+G5t = [
+  Gadget.join([CUBIC_VERTEX] * 3, [((1,1), (2,1)), ((2,2), (3,2))], list(out))
   for out in permutations([ (1,2), (1,3), (2,3), (3,1), (3,3) ])
-] + [
-  Gadget.join([ CUBIC_VERTEX, FREE_EDGE ], [], out)
+]
+
+G5f = [
+  Gadget.join([ CUBIC_VERTEX, FREE_EDGE ], [], list(out))
   for out in permutations([ (1,1), (1,2), (1,3), (2,1), (2,2) ])
 ]
 
 G5s = [
   Gadget.join([CUBIC_VERTEX] * 3, [((1,1), (2,1)), ((2,2), (3,2))], rot([ (1,2), (1,3), (2,3), (3,1), (3,3) ], i))
   for i in range(5)
-] + [
-  Gadget.join([ CUBIC_VERTEX, FREE_EDGE ], [], rot([ (1,1), (1,2), (1,3), (2,1), (2,2) ], i))
-  for i in range(5)
 ]
 
 G6 = [
-  Gadget.join([CUBIC_VERTEX]*4, [((1,1), (2,1)), ((1,2), (3,1)), ((1,3), (4,1))], out)
-  for out in [[ (2,2), (2,3), (3,2), (3,3), (4,2), (4,3) ]]
+  Gadget.join([FREE_EDGE]*3, [], rot([ (1,1), (1,2), (2,1), (2,2), (3,1), (3,2) ], i))
+  for i in range(6)
 ]
 
 if __name__ == "__main__":
-  # print(solve_gadget(G4[0], G4f))
-  l = sorted(friend_boundary(sun(5)), key=lambda x: x[1])
-  print(len(l))
-  for x in l: print(x)
+  import sys
+  verbose = "-v" in sys.argv
+
+#  print("6-cycle replaced by matchings: %s" % (
+#    solve_gadget(sun(6), G6),))
+
+  print("3-cycle replaced by vertex: %s" % (
+    solve_gadget(sun(3), [ CUBIC_VERTEX ]),))
+  print("4-cycle replaced by all 3 matchings: %s" % (
+    solve_gadget(sun(4), G4f),))
+  print("4-cycle replaced by noncrossing matchings: %s" % (
+    solve_gadget(sun(4), G4f[:2]),))
+  print("5-cycle replaced by a non-crossing cubic vertex and an edge: %s" % (
+    solve_gadget(sun(5), G5fp),))
+  print("5-cycle replaced by a tree: %s" % (
+    solve_gadget(sun(5), G5t),))
+  print("5-cycle replaced by a cubic vertex and an edge: %s" % (
+    solve_gadget(sun(5), G5f),))
+  print("5-cycle replaced by a non-crossing tree: %s" % (
+    solve_gadget(sun(5), G5s),))
+#  l = sorted(friend_boundary(sun(5)), key=lambda x: x[1])
+#  print(len(l))
+#  for x in l: print(x)