From f250d95023f28526d468c9b7df347ddb58fdc2fd Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Radek=20Hu=C5=A1ek?= <husek@iuuk.mff.cuni.cz>
Date: Mon, 24 Jan 2022 13:57:47 +0100
Subject: [PATCH] add experiments/redice-cycle

---
 experiments/reduce-cycle.py | 163 ++++++++++++++++++++++++++++++++++++
 1 file changed, 163 insertions(+)
 create mode 100755 experiments/reduce-cycle.py

diff --git a/experiments/reduce-cycle.py b/experiments/reduce-cycle.py
new file mode 100755
index 0000000..4023f00
--- /dev/null
+++ b/experiments/reduce-cycle.py
@@ -0,0 +1,163 @@
+#!/usr/bin/python
+
+"""
+Application of Theorem 6.12 to triangles, 4-cycles (Theorem 6.13)
+and 5.cycles (Theorem 6.14).
+"""
+
+import sys, os
+sys.path.append(os.path.dirname(__file__) + "/..")
+
+from graph_tools.base import *
+from graph_tools.parameters import CircuitDoubleCover, VertexCount
+from graph_tools.misc import sun
+
+from sage.all import MixedIntegerLinearProgram
+from itertools import permutations
+
+verbose = False
+
+def print_v(fmt, *args):
+  if verbose:
+    print(fmt % args)
+
+def solve_gadget(g, alt_gadgets, c2, param=CircuitDoubleCover):
+  N = g.size()
+  print_v("N: %i", N)
+
+  g_n = VertexCount.finalize(g.eval_gadget(VertexCount))
+
+  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(param.enumerate_boundaries(N))
+  print_v("boundaries: %s", BOUNDARIES)
+
+  def fake_gadget(b):
+    return FakeGadget(N, [ BoundaryValue(b, 1) ])
+
+  def join(a, b):
+    return Gadget.join([ a, b ], [ ((1, i), (2, i)) for i in range(1,N+1) ], [])
+
+  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))
+
+    const = sum( a * V[b] for a, b in zip(mv, BOUNDARIES) ) >= c2**((n - g_n) / 2)
+    P.add_constraint(const)
+    print_v("constraint: %s", const)
+    for i, a in enumerate(mv):
+      dual_constraints[i] += a * V_dual[var]
+    dual_objective += c2**((n - g_n) / 2) * V_dual[var]
+
+  for i, ag in enumerate(alt_gadgets): gadget_to_constraint(ag, i)
+  print_v("dual objective: %s", dual_objective)
+  P_dual.set_objective(dual_objective)
+
+  mv = gadget_to_multiplicity_vector(g)
+  objective = sum( a * V[b] for a, b in zip(mv, BOUNDARIES) )
+  print_v("objective: %s", objective)
+  P.set_objective(objective)
+
+  for b, const in zip(mv, dual_constraints):
+    print_v("dual constraint: %s", const <= b)
+    P_dual.add_constraint(const <= b)
+
+  print_v("%s", P)
+  print_v("%s", P_dual)
+
+  ret = ( P.solve(), P_dual.solve() )
+
+  #P_dual.show()
+  for i, v in sorted(P_dual.get_values(V_dual).items()):
+    print_v('w_%s = %s', i, v)
+
+  return ret
+
+
+def friend_boundary(g, param):
+  N = g.size();
+  BOUNDARIES = list(param.enumerate_boundaries(N))
+
+  def test_boundary(b):
+    fg = FakeGadget(N, [ BoundaryValue(b, 1) ])
+    joins = [ ((1, i), (2, i)) for i in range(1,N+1) ]
+    return (b, Gadget.join([ fg, g ], joins, []).eval(param))
+
+  return [ test_boundary(b) for b in BOUNDARIES ]
+
+
+def rot(l, x):
+  return l[x:] + l[:x]
+
+
+G4 = [
+  Gadget.join([CUBIC_VERTEX] * 2, [((1,1), (2,1))], [ (1,2), (1,3), (2,2), (2,3) ]),
+  Gadget.join([CUBIC_VERTEX] * 2, [((1,1), (2,1))], [ (1,2), (2,2), (1,3), (2,3) ]),
+  Gadget.join([CUBIC_VERTEX] * 2, [((1,1), (2,1))], [ (1,2), (2,2), (2,3), (1,3) ]),
+]
+
+G4f = [
+  Gadget.join([FREE_EDGE]*2, [], [ (1,1), (1,2), (2,1), (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) ]),
+]
+
+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) ])
+]
+
+G5f = [
+  Gadget.join([ CUBIC_VERTEX, FREE_EDGE ], [], list(out))
+  for out in permutations([ (1,1), (1,2), (1,3), (2,1), (2,2) ])
+]
+
+G5f_smart = G5fp + [
+  Gadget.join([ CUBIC_VERTEX, FREE_EDGE ], [], rot([ (1,1), (1,2), (2,1), (1,3), (2,2) ], i))
+  for i in range(5)
+]
+
+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)
+]
+
+
+if __name__ == "__main__":
+  import sys
+  verbose = "-v" in sys.argv
+
+  print("3-cycle replaced by vertex: %s" % (
+    solve_gadget(sun(3), [ CUBIC_VERTEX ], 2),))
+
+  print("4-cycle replaced by all 3 matchings (c = sqrt(%s)): %s" % (
+    2, solve_gadget(sun(4), G4f, 2),))
+  print("4-cycle replaced by noncrossing matchings (c = sqrt(%s)): %s" % (
+    2, solve_gadget(sun(4), G4f[:2], 2),))
+
+  print("5-cycle replaced by a non-crossing cubic vertex and an edge (c = sqrt(%s)): %s" % (
+    (5/2)**0.5, solve_gadget(sun(5), G5fp, (5/2)**0.5),))
+  print("5-cycle replaced by a tree (c = sqrt(%s)): %s" % (
+    3.75**0.5, solve_gadget(sun(5), G5t, 3.75**0.5),))
+  print("5-cycle replaced by a cubic vertex and an edge (c = sqrt(%s)): %s" % (
+    3.75**0.5, solve_gadget(sun(5), G5f, 3.75**0.5),))
+  print("5-cycle replaced by a cubic vertex and an edge (smart) (c = sqrt(%s)): %s" % (
+    3.75**0.5, solve_gadget(sun(5), G5f_smart, 3.75**0.5),))
+
-- 
GitLab