#!/usr/bin/env python3
# Christoph Dürr, Mathieu Mari, Ulrike Schmidt-Kraepelin 2022
# verifies that on the Petersen graph the ratio is 5/3

# here is how we numbered the edges of the graph
# https://www.lip6.fr/Christoph.Durr/pricing/naming-convention.png

"""
Here is how we number (vertices) and edges.

                    (3)
                 .   :   .
              .      8      .
           3         :         2
        .           (8)           .
     .              . .             .
  .                .   .              .
(4)....           13   11         .... (2)
 .    9...        .     .       ..7    .
 .        (9) .......12......(7)       .
  .            14 .     .   10         .
   .             .       .             .
   .             .  . .   .            .
    4           . .    .   .           1
     .          (5)      (6)          .
     .        .             .         .
      .     5                 6      .
       .  .                     .   .
        (0)............0...........(1)
"""

from pulp import *

glpk = GLPK()

def rotate(S):
	""" S is a set of edges, and we rotate then counter clock wise
	"""
	return [5 * (e//5) + ((e+1) % 5) for e in S]

# alternating cycles to avoid
C = [[[5,6,7,9],[1,4,10,14]],
	 [[0,7,9],[1,4,12]],
	 [[2,9,10,11],[3,7,13,14]],
	 [[0,10,14],[5,6,12]],
	 [[1,3,5,11],[2,4,6,13]]]

# complete C with rotations
for r in range(4 * 5):
		S, T = C[-5]
		C.append([rotate(S),rotate(T)])

# graph[i] = set of edges adjacent to vertex i
graph = [[0,4,5], [5,10,13]]
# complete with rotations
for r in range(2 * 4):
	graph.append(rotate(graph[-2]))

# adj[i] = set of edges sharing an endpoint with edge i
adj = [set() for i in range(15)]
for S in graph:
	for i in S:
		for j in S:
			adj[i].add(j)

def can_extend(S, i):
	"""is S union {i} a matching ?
	""" 
	for j in adj[i]:
			if j in S:
				return False
	return True	

# M[k] = list of all matchings of size k, M[0] we don't care 
M = [None, {frozenset([i]) for i in range(15)}]
for k in (2,3,4):
	M.append(set())
	for S in M[-2]:
		for i in range(15):
			if can_extend(S, i):
				Si = S | {i}
				M[-1].add(Si)

prob = LpProblem("Petersen")
# x[i]=1 iff edge i is kept
x = [LpVariable("x%i"%i, cat='Binary') for i in range(15)]
# the set described by x is keepable
for S,T in C:
	# S,T should not be alternating for x (S included in x and T disjoint from x)
	diff = lpSum(x[i] for i in S) - lpSum(x[i] for i in T)
	k = len(T)
	prob +=  diff >= 1 - k
	prob +=  diff <= k - 1
# x contains a matching of size 4
y4 = []
for S in M[4]:
	# yS => S is contained in x
	name = "y" + "_".join(str(i) for i in S)
	yS = LpVariable(name, cat='Binary') 
	xS = lpSum(x[i] for i in S)
	y4.append(yS)
	prob += xS >= 4 * yS 
# x must contain at least one matching of size 4
prob += lpSum(y4) >= 1
# minmax matching is at least 4
for S in M[3]:
	# matching S={a,b,c} is in x => zS
	name = "z" + "_".join(str(i) for i in S)
	zS = LpVariable(name, cat='Binary') 
	xS = lpSum(x[i] for i in S)
	prob += xS <= 2 + zS
	# if zS, then there is an edge d such that S+d is a matching and it is in x
	prob += lpSum(x[d] for d in range(15) if can_extend(S, d)) >= zS

prob.solve(glpk)