# -*- coding: utf-8 -*-
"""
Created on Wed Nov 4 10:21:04 2015
@author: Fabrice et Jill-Jênn
"""
import copy
from collections import deque
import numpy
from scipy.optimize import linprog
class SCPProblem:
def __init__(self, filename = None):
self.nodeNames = {}
self.nodeDistances = {}
self.nodeCoordinates = {}
self.nodeCount = 0
if filename:
self.read(filename)
def read(self, filename):
with open(filename, 'r') as f:
self.nodeCount = 0
for line in f:
if line[0] == 'c':
continue
elif line[0] == 'p':
continue
elif line[0] == 'v':
c, u, x, y = line.split(' ', 4)
if not str(u) in self.nodeNames:
self.nodeNames[u] = self.nodeCount
self.nodeCoordinates[u] = (float(x), float(y))
self.nodeCount += 1
elif line[0] == 'a':
c, u, v, d = line.split(' ', 4)
if not u in self.nodeDistances:
self.nodeDistances[u] = {}
if not v in self.nodeDistances:
self.nodeDistances[v] = {}
self.nodeDistances[u][v] = float(d)
self.nodeDistances[v][u] = float(d)
else:
print("Skipped: {0}".format(line))
continue
def get_remaining_variables(a):
if not a:
return []
return sorted(list(a.values())[0].keys())
def display_constraints(a):
variables = get_remaining_variables(a)
print(' ' + ''.join(variables))
for key in a:
print(key + ''.join('1' if a[key][var] else ' ' for var in variables))
def allPairsDistances(problem):
'''Calcul des plus courtes distances entre toutes les paires de noeuds.'''
V = problem.nodeNames.keys()
# print(V)
weight = {u: {v: float('inf') for v in V} for u in V}
for u in V:
weight[u][u] = 0
for u in problem.nodeDistances:
for v in problem.nodeDistances[u]:
weight[u][v] = problem.nodeDistances[u][v]
# weight[u][v] : poids de l'arête de u à v si elle existe, 0 si u = v, infini sinon
for k in V:
for u in V:
for v in V:
weight[u][v] = min(weight[u][v],
weight[u][k] + weight[k][v])
# weight[u][v] : distance de u à v
return weight
def solveInit(problem, d):
'''Initialisation du problème.'''
weight = allPairsDistances(problem)
V = problem.nodeNames.keys()
a = {}
X = {}
for u in V:
X[u] = None # X[u] = None pour chaque variable non encore affectée
a[u] = {}
for v in V:
# a[u][v] = True if weight[u][v] <= d else False
a[u][v] = weight[u][v] <= d
return a, X
def solveRedColonne(a, X):
'''Réduction du problème par colonne.'''
n = len(X)
variables = get_remaining_variables(a)
# print(variables)
# display_constraints(a)
removeCols = set()
for j in variables:
for k in variables:
if all(a[i][j] <= a[i][k] for i in a) and \
any(a[i][j] < a[i][k] for i in a):
removeCols.add(j)
# print('Colonnes à supprimer', removeCols)
for j in removeCols:
for i in a:
if j in a[i]:
del a[i][j]
X[j] = 0
return a, X
def solveRedLigne1(a, X):
'''Réduction du problème par ligne, première technique.'''
# display_constraints(a)
n = len(X)
vars_at_1 = set()
for i in a:
s = 0
for j in a[i]:
s += a[i][j]
if a[i][j] == True:
pos_of_last_seen_1 = j
if s == 1:
vars_at_1.add(pos_of_last_seen_1)
# if sum(a[i]) == 1: # a.count(1) == 1 dans le cas matrice
# j = a[i].index(1)
# print('Variables trouvées', vars_at_1)
for j in vars_at_1:
X[j] = 1
# print('Variables', X)
removeRows = set()
for i in a:
for j in vars_at_1:
if a[i][j] == True:
removeRows.add(i)
del a[i][j] # On supprime la colonne même si la variable est affectée à False
# print('Lignes à supprimer', removeRows)
for i in removeRows:
del a[i]
return a, X
def solveRedLigne2(a, X):
'''Réduction du problème par ligne, seconde technique.'''
n = len(X)
removeRows = set()
for m in a:
for n in a:
if (all(a[m][j] <= a[n][j] for j in a[m]) and
any(a[m][j] < a[n][j] for j in a[m])):
removeRows.add(n)
# print('Lignes à supprimer', removeRows)
for i in removeRows:
del a[i]
return a, X
def solveIsSolution(n):
'''Teste si le noeud n est solution.'''
a, X = n
# ???
return len(a) == 0
def solveGetObjective(n):
'''Fonction objectif.'''
a, X = n
s = 0
for j in X:
if X[j] is not None:
s += X[j]
return s
def solveSplit(n):
'''Séparation du noeud n. Renvoie la liste des fils.'''
a, X = n
for j in X:
if X[j] is None:
# print('Variable de split', j)
# aj0 : on supprime la colonne j
aj0 = copy.deepcopy(a)
Xj0 = copy.deepcopy(X)
Xj0[j] = 0
removeRows = set()
for i in aj0:
del aj0[i][j]
if len(aj0[i]) == 0:
removeRows.add(i)
# aj0, Xj0 = solveRedColonne(aj0, Xj0)
# aj0, Xj0 = solveRedLigne1(aj0, Xj0)
# aj0, Xj0 = solveRedLigne2(aj0, Xj0)
for i in removeRows:
del aj0[i]
# aj1 : on supprime toutes les lignes i ayant a[i][j] == 1
aj1 = copy.deepcopy(a)
Xj1 = copy.deepcopy(X)
Xj1[j] = 1
removeRows = set()
for i in aj1:
if aj1[i][j] == True:
removeRows.add(i)
del aj1[i][j]
if len(aj1[i]) == 0:
removeRows.add(i)
for i in removeRows:
del aj1[i]
# aj1, Xj1 = solveRedColonne(aj1, Xj1)
# aj1, Xj1 = solveRedLigne1(aj1, Xj1)
# aj1, Xj1 = solveRedLigne2(aj1, Xj1)
for i in aj1:
assert len(aj1[i]) > 0
return [(aj0, Xj0), (aj1, Xj1)]
def solveEvaluate(n):
'''Évaluation du noeud n.'''
# Ô simplexe, renvoie-moi la valeur de l'instance
# Reconstruire l'instance en passant à l'opposé
# Pour tout i, ∑_j -a[i][j] <= -1
a, X = n
# display_constraints(a)
# print(X)
nb_variables = 0
for j in X:
if X[j] is None:
nb_variables += 1
instance_A = []
instance_b = []
for i in a:
instance_A.append([-a[i][j] for j in a[i]])
instance_b.append(-1)
# print(X)
# print(nb_variables, instance_A, instance_b)
# print(linprog([1] * nb_variables, instance_A, instance_b, bounds=[(0, 1)] * nb_variables).x)
result = linprog([1] * nb_variables, instance_A, instance_b, bounds=[(0, 1)] * nb_variables)
if result.status != 0:
return float('inf') # No solution found
return solveGetObjective(n) + sum(result.x)
def solveBandB(a, X):
'''Branch-and-Bound solver.'''
# Le code le plus important de ce TP
todo = []
todo.append((a, X))
best = float('inf'), (a, X)
while todo:
n = todo.pop()
if solveIsSolution(n):
obj = solveGetObjective(n)
if obj < best[0]:
best = obj, n
else:
children = solveSplit(n)
for child in children:
borne_min_sur_obj = solveEvaluate(n)
# print(borne_min_sur_obj)
if borne_min_sur_obj > best[0]:
# Même Dieu n'a pas trouvé mieux
continue
else:
todo.append(child)
return best
def solve(problem, dCouverture):
X = {}
o = 0
print('\n\nSolving for dCouv = ', dCouverture)
a, X = solveInit(problem, dCouverture)
# a est une matrice de booléens
# display_constraints(a)
# réduction du problème par colonne
a, X = solveRedColonne(a, X)
# réduction du problème par ligne (1)
a, X = solveRedLigne1(a, X)
# réduction du problème par ligne (2)
a, X = solveRedLigne2(a, X)
o = sum((1 if X[i] else 0) for i in X)
# Si le problème n'est pas encore résolu
# il faut passer par Branch and Bound
if (len(a)):
print('Problem is not yet solved.\nSwitching to Branch and Bound')
o, X = solveBandB(a, X)
print('Problem is done')
return o, X
"""pb = SCPProblem('scp-ex1.gr')
for i in range(3, 20):
print('For dCouv = '+str(i), solve(pb, i))"""
pb = SCPProblem('france-1.gr')
print(solve(pb, 130))
