Correction du set covering problem

# -*- 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))