Débit maximal et problème d'affectation linéaire

Voici un problème : votre entreprise affecte des sous-traitants pour exécuter les contrats. Vous parcourez vos listes et décidez quels sous-traitants sont disponibles pour un engagement d'un mois et vous parcourez vos contrats disponibles pour voir lesquels d'entre eux sont pour des tâches d'un mois. Étant donné que vous savez avec quelle efficacité chaque sous-traitant peut exécuter chaque contrat, comment affectez-vous les sous-traitants pour maximiser l'efficacité globale pour ce mois ?

Ceci est un exemple du problème d'affectation, et le problème peut être résolu avec l'algorithme hongrois classique.

L'algorithme hongrois (également connu sous le nom d'algorithme de Kuhn-Munkres) est un algorithme de temps polynomial qui maximise la correspondance des poids dans un graphe bipartite pondéré. Ici, les sous-traitants et les contrats peuvent être modélisés sous la forme d'un graphe bipartite, avec leur efficacité en tant que poids des arêtes entre le sous-traitant et les nœuds de contrat.

Dans cet article, vous découvrirez une implémentation de l'algorithme hongrois qui utilise l'algorithme d'Edmonds-Karp pour résoudre le problème d'affectation linéaire. Vous apprendrez également en quoi l'algorithme d'Edmonds-Karp est une légère modification de la méthode Ford-Fulkerson et en quoi cette modification est importante.

Le problème du débit maximal

Le problème de débit maximal lui-même peut être décrit de manière informelle comme le problème du déplacement d'un fluide ou d'un gaz à travers un réseau de tuyaux d'une source unique à un terminal unique. Cela se fait en supposant que la pression dans le réseau est suffisante pour garantir que le fluide ou le gaz ne peut pas s'attarder dans n'importe quelle longueur de tuyau ou de raccord de tuyau (les endroits où différentes longueurs de tuyau se rencontrent).

En apportant certaines modifications au graphique, le problème d'affectation peut être transformé en un problème de flux maximum.

Préliminaires

Les idées nécessaires pour résoudre ces problèmes surgissent dans de nombreuses disciplines mathématiques et d'ingénierie, souvent des concepts similaires sont connus sous des noms différents et exprimés de différentes manières (par exemple, les matrices d'adjacence et les listes d'adjacence). Étant donné que ces idées sont assez ésotériques, des choix sont faits concernant la manière dont ces concepts seront généralement définis pour un contexte donné.

Cet article ne supposera aucune connaissance préalable au-delà d'une petite introduction à la théorie des ensembles.

Les implémentations de cet article représentent les problèmes sous forme de graphes orientés (digraphe).

DiGraphs

Un digraphe a deux attributs, setOfNodes et setOfArcs . Ces deux attributs sont des ensembles (collections non ordonnées). Dans les blocs de code de cet article, j'utilise en fait le frozenset de Python, mais ce détail n'est pas particulièrement important.

 DiGraph = namedtuple('DiGraph', ['setOfNodes','setOfArcs'])

(Remarque : ce code, ainsi que tous les autres codes de cet article, sont écrits en Python 3.6.)

Nœuds

Un nœud n est composé de deux attributs :

• n.uid : Un identifiant unique.

  Cela signifie que pour deux nœuds x et y ,

 x.uid != y.uid

• n.datum : Cela représente n'importe quel objet de données.

 Node = namedtuple('Node', ['uid','datum'])

Arcs

Un arc a est composé de trois attributs :

• a.fromNode : Il s'agit d'un nœud , tel que défini ci-dessus.

• a.toNode : Il s'agit d'un nœud , tel que défini ci-dessus.

• a.datum : Cela représente n'importe quel objet de données.

 Arc = namedtuple('Arc', ['fromNode','toNode','datum'])

L'ensemble des arcs dans le digraphe représente une relation binaire sur les nœuds du digraphe . L'existence de l'arc a implique qu'une relation existe entre a.fromNode et a.toNode .

Dans un graphe orienté (par opposition à un graphe non orienté), l'existence d'une relation entre a.fromNode et a.toNode n'implique pas qu'il existe une relation similaire entre a.toNode et a.fromNode .

En effet, dans un graphe non orienté, la relation exprimée n'est pas nécessairement symétrique.

DiGraphs

Les nœuds et les arcs peuvent être utilisés pour définir un digraphe , mais pour plus de commodité, dans les algorithmes ci-dessous, un digraphe sera représenté en utilisant comme dictionnaire.

Voici une méthode qui peut convertir la représentation graphique ci-dessus en une représentation de dictionnaire similaire à celle-ci :

 def digraph_to_dict(G): G_as_dict = dict([]) for a in G.setOfArcs: if(a.fromNode not in G.setOfNodes): err_msg = 'There is no Node {a.fromNode.uid!s} to match the Arc from {a.fromNode.uid!s} to {a.toNode.uid!s}'.format(**locals()) pdg(G) raise KeyError(err_msg) if(a.toNode not in G.setOfNodes): err_msg = 'There is no Node {a.toNode.uid!s} to match the Arc from {a.fromNode.uid!s} to {a.toNode.uid!s}'.format(**locals()) pdg(G) raise KeyError(err_msg) G_as_dict[a.fromNode] = (G_as_dict[a.fromNode].union(frozenset([a]))) if (a.fromNode in G_as_dict) else frozenset([a]) for a in G.setOfArcs: if(a.fromNode not in G.setOfNodes): err_msg = 'There is no Node {a.fromNode.uid!s} to match the Arc from {a.fromNode.uid!s} to {a.toNode.uid!s}'.format(**locals()) raise KeyError(err_msg) if a.toNode not in G_as_dict: G_as_dict[a.toNode] = frozenset([]) return G_as_dict

Et en voici un autre qui le convertit en dictionnaire de dictionnaires, autre opération qui sera utile :

 def digraph_to_double_dict(G): G_as_dict = dict([]) for a in G.setOfArcs: if(a.fromNode not in G.setOfNodes): err_msg = 'There is no Node {a.fromNode.uid!s} to match the Arc from {a.fromNode.uid!s} to {a.toNode.uid!s}'.format(**locals()) raise KeyError(err_msg) if(a.toNode not in G.setOfNodes): err_msg = 'There is no Node {a.toNode.uid!s} to match the Arc from {a.fromNode.uid!s} to {a.toNode.uid!s}'.format(**locals()) raise KeyError(err_msg) if(a.fromNode not in G_as_dict): G_as_dict[a.fromNode] = dict({a.toNode : frozenset([a])}) else: if(a.toNode not in G_as_dict[a.fromNode]): G_as_dict[a.fromNode][a.toNode] = frozenset([a]) else: G_as_dict[a.fromNode][a.toNode] = G_as_dict[a.fromNode][a.toNode].union(frozenset([a])) for a in G.setOfArcs: if(a.fromNode not in G.setOfNodes): err_msg = 'There is no Node {a.fromNode.uid!s} to match the Arc from {a.fromNode.uid!s} to {a.toNode.uid!s}'.format(**locals()) raise KeyError(err_msg) if a.toNode not in G_as_dict: G_as_dict[a.toNode] = dict({}) return G_as_dict

Lorsque l'article parle d'un digraphe représenté par un dictionnaire, il mentionnera G_as_dict pour s'y référer.

Parfois, il est utile de récupérer un nœud à partir d'un digraphe G en passant par son uid (identifiant unique):

 def find_node_by_uid(find_uid, G): nodes = {n for n in G.setOfNodes if n.uid == find_uid} if(len(nodes) != 1): err_msg = 'Node with uid {find_uid!s} is not unique.'.format(**locals()) raise KeyError(err_msg) return nodes.pop()

Dans la définition des graphes, les gens utilisent parfois les termes nœud et sommet pour désigner le même concept ; il en est de même des termes arc et arête.

Deux représentations de graphes populaires en Python sont celle-ci qui utilise des dictionnaires et une autre qui utilise des objets pour représenter des graphes. La représentation dans cet article se situe quelque part entre ces deux représentations couramment utilisées.

Ceci est ma représentation digraphique . Il y en a beaucoup comme ça, mais celui-ci est le mien.

Promenades et Sentiers

Soit S_Arcs une séquence finie (collection ordonnée) d' arcs dans un digraphe G tel que si a est un arc quelconque dans S_Arcs à l'exception du dernier, et b suit a dans la séquence, alors il doit y avoir un nœud n = a.fromNode dans G.setOfNodes tel que a.toNode = b.fromNode .

En commençant par le premier arc dans S_Arcs , et en continuant jusqu'au dernier arc dans S_Arcs , collectez (dans l'ordre rencontré) tous les nœuds n tels que définis ci-dessus entre chacun des deux arcs consécutifs dans S_Arcs . Étiquetez la collection ordonnée de nœuds collectés lors de cette opération S_{Nodes} .

 def path_arcs_to_nodes(s_arcs): s_nodes = list([]) arc_it = iter(s_arcs) step_count = 0 last = None try: at_end = False last = a1 = next(arc_it) while (not at_end): s_nodes += [a1.fromNode] last = a2 = next(arc_it) step_count += 1 if(a1.toNode != a2.fromNode): err_msg = "Error at index {step_count!s} of Arc sequence.".format(**locals()) raise ValueError(err_msg) a1 = a2 except StopIteration as e: at_end = True if(last is not None): s_nodes += [last.toNode] return s_nodes

• Si un nœud apparaît plus d'une fois dans la séquence S_Nodes alors appelez S_Arcs a Walk on digraph G .

• Sinon, appelez S_Arcs un chemin de list(S_Nodes)[0] à list(S_Nodes)[-1] sur le digraphe G .

Nœud source

Appelez le nœud s un nœud source dans le digraphe G si s est dans G.setOfNodes et G.setOfArcs ne contient pas d' arc a tel que a.toNode = s .

 def source_nodes(G): to_nodes = frozenset({a.toNode for a in G.setOfArcs}) sources = G.setOfNodes.difference(to_nodes) return sources

Nœud terminal

Appelez le nœud t un nœud terminal dans le digraphe G si t est dans G.setOfNodes et G.setOfArcs ne contient pas d' arc a tel que a.fromNode = t .

 def terminal_nodes(G): from_nodes = frozenset({a.fromNode for a in G.setOfArcs}) terminals = G.setOfNodes.difference(from_nodes) return terminals

Coupes et st Coupes

Une cut d'un digraphe connexe G est un sous-ensemble d' arcs de G.setOfArcs qui partitionne l'ensemble des nœuds G.setOfNodes dans G . G est connexe si chaque nœud n dans G.setOfNodes et a au moins un arc a dans G.setOfArcs tel que n = a.fromNode ou n = a.toNode , mais a.fromNode != a.toNode .

 Cut = namedtuple('Cut', ['setOfCutArcs'])

La définition ci-dessus fait référence à un sous-ensemble d' arcs , mais elle peut également définir une partition des nœuds de G.setOfNodes .

Pour les fonctions predecessors_of et successors_of , n est un nœud dans l'ensemble G.setOfNodes du digraphe G , et cut est une coupe de G :

 def cut_predecessors_of(n, cut, G): allowed_arcs = G.setOfArcs.difference(frozenset(cut.setOfCutArcs)) predecessors = frozenset({}) previous_count = len(predecessors) reach_fringe = frozenset({n}) keep_going = True while( keep_going ): reachable_from = frozenset({a.fromNode for a in allowed_arcs if (a.toNode in reach_fringe)}) reach_fringe = reachable_from predecessors = predecessors.union(reach_fringe) current_count = len(predecessors) keep_going = current_count - previous_count > 0 previous_count = current_count return predecessors def cut_successors_of(n, cut, G): allowed_arcs = G.setOfArcs.difference(frozenset(cut.setOfCutArcs)) successors = frozenset({}) previous_count = len(successors) reach_fringe = frozenset({n}) keep_going = True while( keep_going ): reachable_from = frozenset({a.toNode for a in allowed_arcs if (a.fromNode in reach_fringe)}) reach_fringe = reachable_from successors = successors.union(reach_fringe) current_count = len(successors) keep_going = current_count - previous_count > 0 previous_count = current_count return successors

Soit cut une coupe du digraphe G .

 def get_first_part(cut, G): firstPartFringe = frozenset({a.fromNode for a in cut.setOfCutArcs}) firstPart = firstPartFringe for n in firstPart: preds = cut_predecessors_of(n,cut,G) firstPart = firstPart.union(preds) return firstPart def get_second_part(cut, G): secondPartFringe = frozenset({a.toNode for a in cut.setOfCutArcs}) secondPart = secondPartFringe for n in secondPart: succs= cut_successors_of(n,cut,G) secondPart = secondPart.union(succs) return secondPart

cut est une coupe du digraphe G si : (get_first_part(cut, G).union(get_second_part(cut, G)) == G.setOfNodes) and (len(get_first_part(cut, G).intersect(get_second_part(cut, G))) == 0) cut est appelé une coupe xy si (x in get_first_part(cut, G)) and (y in get_second_part(cut, G) ) and (x != y) . Lorsque le nœud x dans une cut xy est un nœud source et que le nœud y dans la coupe xy est un nœud terminal , alors cette coupe est appelée une ère coupe .

 STCut = namedtuple('STCut', ['s','t','cut'])

Réseaux de flux

Vous pouvez utiliser un digramme G pour représenter un réseau de flux.

Attribuez à chaque nœud n , où n est dans G.setOfNodes un n.datum qui est un FlowNodeDatum :

 FlowNodeDatum = namedtuple('FlowNodeDatum', ['flowIn','flowOut'])

Attribuez à chaque arc a , où a est dans G.setOfArcs et a.datum qui est un FlowArcDatum .

 FlowArcDatum = namedtuple('FlowArcDatum', ['capacity','flow'])

flowNodeDatum.flowIn et flowNodeDatum.flowOut sont des nombres réels positifs.

flowArcDatum.capacity et flowArcDatum.flow sont également des nombres réels positifs.

Pour chaque nœud node n dans G.setOfNodes :

 n.datum.flowIn = sum({a.datum.flow for a in G.Arcs if a.toNode == n}) n.datum.flowOut = sum({a.datum.flow for a in G.Arcs if a.fromNode == n})

Le digramme G représente maintenant un réseau de flux .

Le flot de G fait référence au a.flow pour tous les arcs a dans G .

Flux réalisables

Soit le digraphe G représentant un réseau de flux .

Le réseau de flux représenté par G a des flux réalisables si :

1. Pour chaque nœud n dans G.setOfNodes à l'exception des nœuds sources et des nœuds terminaux : n.datum.flowIn = n.datum.flowOut .

2. Pour chaque arc a dans G.setOfNodes : a.datum.capacity <= a.datum.flow .

La condition 1 est appelée une contrainte de conservation.

La condition 2 est appelée contrainte de capacité.

Capacité de coupe

La capacité de coupe d'une st coupe stCut de nœud source s et de nœud terminal t d'un réseau de flux représenté par un digraphe G est :

 def cut_capacity(stCut, G): part_1 = get_first_part(stCut.cut,G) part_2 = get_second_part(stCut.cut,G) s_part = part_1 if stCut.s in part_1 else part_2 t_part = part_1 if stCut.t in part_1 else part_2 cut_capacity = sum({a.datum.capacity for a in stCut.cut.setOfCutArcs if ( (a.fromNode in s_part) and (a.toNode in t_part) )}) return cut_capacity

Capacité minimale de coupe

Soit stCut = stCut(s,t,cut) une ère coupe d'un réseau de flux représenté par un digraphe G .

stCut est la coupure de capacité minimale du réseau de flux représenté par G s'il n'y a pas d'autre stCutCompetitor de st cut dans ce réseau de flux tel que :

 cut_capacity(stCut, G) < cut_capacity(stCutCompetitor, G)

Éliminer les flux

Je voudrais faire référence au digraphe qui serait le résultat de la prise d'un digraphe G et de la suppression de toutes les données de flux de tous les nœuds de G.setOfNodes ainsi que de tous les arcs de G.setOfArcs .

 def strip_flows(G): new_nodes= frozenset( (Node(n.uid, FlowNodeDatum(0.0,0.0)) for n in G.setOfNodes) ) new_arcs = frozenset([]) for a in G.setOfArcs: new_fromNode = Node(a.fromNode.uid, FlowNodeDatum(0.0,0.0)) new_toNode = Node(a.toNode.uid, FlowNodeDatum(0.0,0.0)) new_arc = Arc(new_fromNode, new_toNode, FlowArcDatum(a.datum.capacity, 0.0)) new_arcs = new_arcs.union(frozenset([new_arc])) return DiGraph(new_nodes, new_arcs)

Problème de débit maximal

Un réseau de flux représenté par un digraphe G , un nœud source s dans G.setOfNodes et un nœud terminal t dans G.setOfNodes , G peut représenter un problème de flux maximal si :

 (len(list(source_nodes(G))) == 1) and (next(iter(source_nodes(G))) == s) and (len(list(terminal_nodes(G))) == 1) and (next(iter(terminal_nodes(G))) == t)

Étiquetez cette représentation :

 MaxFlowProblemState = namedtuple('MaxFlowProblemState', ['G','sourceNodeUid','terminalNodeUid','maxFlowProblemStateUid'])

Où sourceNodeUid = s.uid , terminalNodeUid = t.uid et maxFlowProblemStateUid est un identifiant pour l'instance du problème.

Solution de débit maximal

Soit maxFlowProblem représenter un problème de débit maximum . La solution à maxFlowProblem peut être représentée par un réseau de flux représenté par un digraphe H .

Le digraphe H est une solution réalisable au problème de débit maximal sur l'entrée python maxFlowProblem si :

1. strip_flows(maxFlowProblem.G) == strip_flows(H) .

2. H est un réseau de flux et a des flux réalisables .

Si en plus de 1 et 2 :

1. Il ne peut y avoir aucun autre réseau de flux représenté par le digraphe K tel que strip_flows(G) == strip_flows(K) et find_node_by_uid(t.uid,G).flowIn < find_node_by_uid(t.uid,K).flowIn .

Alors H est aussi une solution optimale à maxFlowProblem .

En d'autres termes, une solution de débit maximal réalisable peut être représentée par un digraphe , qui :

1. Est identique au digraphe G du problème de flux maximal correspondant à l'exception que le n.datum.flowIn , n.datum.flowOut et le a.datum.flow de l'un des nœuds et des arcs peuvent être différents.

2. Représente un réseau de flux qui a des flux réalisables .

Et, il peut représenter une solution optimale de débit maximum si en plus :

1. Le flowIn pour le nœud correspondant au nœud terminal dans le problème de flux maximal est aussi grand que possible (lorsque les conditions 1 et 2 sont toujours satisfaites).

Si le digraphe H représente une solution de débit maximum réalisable : find_node_by_uid(s.uid,H).flowOut = find_node_by_uid(t.uid,H).flowIn cela découle du débit max, théorème de coupe min (discuté ci-dessous). De manière informelle, puisque H est supposé avoir des flux réalisables , cela signifie que le flux ne peut être ni "créé" (sauf au nœud source s ) ni "détruit" (sauf au nœud terminal t ) tout en traversant un (autre) nœud ( contraintes de conservation ).

Puisqu'un problème de flux maximal ne contient qu'un seul nœud source s et un seul nœud terminal t , tous les flux "créés" en s doivent être "détruits" en t ou le réseau de flux n'a pas de flux réalisables (la contrainte de conservation aurait été violée ).

Soit le digraphe H une solution de débit maximal réalisable ; la valeur ci-dessus est appelée la st Flow Value de H .

Laisser:

 mfps=MaxFlowProblemState(H, maxFlowProblem.sourceNodeUid, maxFlowProblem.terminalNodeUid, maxFlowProblem.maxFlowProblemStateUid)

Cela signifie que mfps est un état successeur de maxFlowProblem , ce qui signifie simplement que mfps est exactement comme maxFlowProblem à l'exception que les valeurs de a.flow pour les arcs a dans maxFlowProblem.setOfArcs peuvent être différentes de a.flow pour les arcs a dans mfps.setOfArcs .

 def get_mfps_flow(mfps): flow_from_s = find_node_by_uid(mfps.sourceNodeUid,mfps.G).datum.flowOut flow_to_t = find_node_by_uid(mfps.terminalNodeUid,mfps.G).datum.flowIn if( (flow_to_t - flow_from_s) > 0): raise Exception('Infeasible st flow') return flow_to_t

Voici une visualisation d'un mfps avec son maxFlowProb associé. Chaque arc a dans l'image a une étiquette, ces étiquettes sont a.datum.flowFrom / a.datum.flowTo , chaque nœud n dans l'image a une étiquette, et ces étiquettes sont n.uid {n.datum.flowIn / a.datum.flowOut} .

premier flux de coupe

Laissez mfps représenter un MaxFlowProblemState et laissez stCut représenter une coupe de mfps.G . Le flux de coupe de stCut est défini :

 def get_stcut_flow(stCut,mfps): s = find_node_by_uid(mfps.sourceNodeUid, mfps.G) t = find_node_by_uid(mfps.terminalNodeUid, mfps.G) part_1 = get_first_part(stCut.cut,mfps.G) part_2 = get_second_part(stCut.cut,mfps.G) s_part = part_1 if s in part_1 else part_2 t_part = part_1 if t in part_1 else part_2 s_t_sum = sum(a.datum.flow for a in mfps.G.setOfArcs if (a in stCut.cut.setOfCutArcs) and (a.fromNode in s_part) ) t_s_sum = sum(a.datum.flow for a in mfps.G.setOfArcs if (a in stCut.cut.setOfCutArcs) and (a.fromNode in t_part) ) cut_flow = s_t_sum - t_s_sum return cut_flow

Le premier débit de coupure est la somme des flux de la partition contenant le nœud source à la partition contenant le nœud terminal moins la somme des flux de la partition contenant le nœud terminal à la partition contenant le nœud source .

Débit maximum, coupe minimum

Supposons que maxFlowProblem représente un problème de flux maximal et que la solution de maxFlowProblem soit représentée par un réseau de flux représenté par le digraphe H .

Soit minStCut la coupure de capacité minimale du réseau de flux représenté par maxFlowProblem.G .

Parce que dans le problème de flux maximal, le flux provient d'un seul nœud source et se termine à un seul nœud terminal et, en raison des contraintes de capacité et des contraintes de conservation , nous savons que l'ensemble du flux entrant dans maxFlowProblem.terminalNodeUid doit traverser n'importe quelle ère coupe , en particulier il doit croiser minStCut . Ça signifie:

 get_flow_value(H, maxFlowProblem) <= cut_capacity(minStCut, maxFlowProblem.G)

Résolution du problème de débit maximal

L'idée de base pour résoudre un problème de débit maximum maxFlowProblem est de commencer avec une solution de débit maximum représentée par le digraphe H . Un tel point de départ peut être H = strip_flows(maxFlowProblem.G) . La tâche consiste alors à utiliser H et par une modification gloutonne des valeurs a.datum.flow de certains arcs a dans H.setOfArcs pour produire une autre solution de flux maximum représentée par le digraphe K tel que K ne puisse toujours pas représenter un réseau de flux avec des flux réalisables et get_flow_value(H, maxFlowProblem) < get_flow_value(K, maxFlowProblem) . Tant que ce processus se poursuit, la qualité ( get_flow_value(K, maxFlowProblem) ) de la dernière solution de débit maximum rencontrée ( K ) est meilleure que toute autre solution de débit maximum trouvée. Si le processus atteint un point où il sait qu'aucune autre amélioration n'est possible, le processus peut se terminer et il renverra la solution de débit maximal optimale.

La description ci-dessus est générale et ignore de nombreuses preuves telles que si un tel processus est possible ou combien de temps cela peut prendre, je donnerai quelques détails supplémentaires et l'algorithme.

Le débit maximal, théorème de coupe minimale

D'après le livre Flows in Networks de Ford et Fulkerson, l'énoncé du théorème max flow, min cut (théorème 5.1) est :

Pour tout réseau, la valeur maximale du débit de s à t est égale à la capacité de coupe minimale de toutes les coupes séparant s et t .

En utilisant les définitions de cet article, cela se traduit par :

La solution à un maxFlowProblem représenté par un réseau de flux représenté par le digraphe H est optimale si :

 get_flow_value(H, maxFlowProblem) == cut_capacity(minStCut, maxFlowProblem.G)

J'aime cette preuve du théorème et Wikipedia en a une autre.

Le théorème débit max, coupe min est utilisé pour prouver l'exactitude et l'exhaustivité de la méthode de Ford-Fulkerson .

Je donnerai également une preuve du théorème dans la section après l' augmentation des chemins .

La méthode Ford-Fulkerson et l'algorithme Edmonds-Karp

CLRS définit la méthode Ford-Fulkerson comme suit (section 26.2) :

 FORD-FULKERSON-METHOD(G, s, t): initialize flow f to 0 while there exists an augmenting path p in the residual graph G_f augment flow f along

Graphique résiduel

Le graphe résiduel d'un réseau de flux représenté par le digraphe G peut être représenté par un digraphe G_f :

 ResidualDatum = namedtuple('ResidualDatum', ['residualCapacity','action']) def agg_n_to_u_cap(n,u,G_as_dict): arcs_out = G_as_dict[n] return sum([a.datum.capacity for a in arcs_out if( (a.fromNode == n) and (a.toNode == u) ) ]) def agg_n_to_u_flow(n,u,G_as_dict): arcs_out = G_as_dict[n] return sum([a.datum.flow for a in arcs_out if( (a.fromNode == n) and (a.toNode == u) ) ]) def get_residual_graph_of(G): G_as_dict = digraph_to_dict(G) residual_nodes = G.setOfNodes residual_arcs = frozenset([]) for n in G_as_dict: arcs_from = G_as_dict[n] nodes_to = frozenset([find_node_by_uid(a.toNode.uid,G) for a in arcs_from]) for u in nodes_to: n_to_u_cap_sum = agg_n_to_u_cap(n,u,G_as_dict) n_to_u_flow_sum = agg_n_to_u_flow(n,u,G_as_dict) if(n_to_u_cap_sum > n_to_u_flow_sum): flow = round(n_to_u_cap_sum - n_to_u_flow_sum, TOL) residual_arcs = residual_arcs.union( frozenset([Arc(n,u,datum=ResidualDatum(flow, 'push'))]) ) if(n_to_u_flow_sum > 0.0): flow = round(n_to_u_flow_sum, TOL) residual_arcs = residual_arcs.union( frozenset([Arc(u,n,datum=ResidualDatum(flow, 'pull'))]) ) return DiGraph(residual_nodes, residual_arcs)

• agg_n_to_u_cap(n,u,G_as_dict) renvoie la somme de a.datum.capacity pour tous les arcs du sous-ensemble de G.setOfArcs où l'arc a est dans le sous-ensemble si a.fromNode = n et a.toNode = u .

• agg_n_to_u_cap(n,u,G_as_dict) renvoie la somme de a.datum.flow pour tous les arcs du sous-ensemble de G.setOfArcs où l'arc a est dans le sous-ensemble si a.fromNode = n et a.toNode = u .

Brièvement, le graphe résiduel G_f représente certaines actions qui peuvent être effectuées sur le digraphe G .

Chaque paire de nœuds n,u dans G.setOfNodes du réseau de flux représenté par le digraphe G peut générer 0, 1 ou 2 arcs dans le graphe résiduel G_f de G .

1. Le couple n,u ne génère aucun G_f s'il n'y a pas d' arc a dans G.setOfArcs tel que a.fromNode = n et a.toNode = u .

2. La paire n,u génère l' arc a dans G_f.setOfArcs où a représente un arc étiqueté arc de flux poussé de n vers u a = Arc(n,u,datum=ResidualNode(n_to_u_cap_sum - n_to_u_flow_sum)) if n_to_u_cap_sum > n_to_u_flow_sum .

3. La paire n,u génère l' arc a dans G_f.setOfArcs où a représente un arc étiqueté arc de flux tiré de n à u a = Arc(n,u,datum=ResidualNode(n_to_u_cap_sum - n_to_u_flow_sum)) if n_to_u_flow_sum > 0.0 .

• Chaque arc de flux de poussée dans G_f.setOfArcs représente l'action d'ajouter un total de flux x <= n_to_u_cap_sum - n_to_u_flow_sum aux arcs dans le sous-ensemble de G.setOfArcs où l'arc a est dans le sous-ensemble si a.fromNode = n et a.toNode = u .

• Chaque arc de flux d'extraction dans G_f.setOfArcs représente l'action de soustraire un total de flux x <= n_to_u_flow_sum aux arcs dans le sous-ensemble de G.setOfArcs où l'arc a est dans le sous-ensemble si a.fromNode = n et a.toNode = u .

L'exécution d'une action push ou pull individuelle à partir de G_f sur les arcs applicables dans G peut générer un réseau de flux sans flux réalisables car les contraintes de capacité ou les contraintes de conservation peuvent être violées dans le réseau de flux généré.

Voici une visualisation du graphique résiduel de l'exemple précédent de visualisation d'une solution de débit maximum , dans la visualisation chaque arc a représente a.residualCapacity .

Chemin d'augmentation

Soit maxFlowProblem un problème de flux max , et soit G_f = get_residual_graph_of(G) le graphe résiduel de maxFlowProblem.G .

Un chemin d' augmentingPath pour maxFlowProblem est n'importe quel chemin de find_node_by_uid(maxFlowProblem.sourceNode,G_f) à find_node_by_uid(maxFlowProblem.terminalNode,G_f) .

Il s'avère qu'un chemin d' augmentingPath peut être appliqué à une solution de débit max représentée par le digraphe H générant une autre solution de débit max représentée par le digraphe K où get_flow_value(H, maxFlowProblem) < get_flow_value(K, maxFlowProblem) si H n'est pas optimal .

Voici comment:

 def augment(augmentingPath, H): augmentingPath = list(augmentingPath) H_as_dict = digraph_to_dict(H) new_nodes = frozenset({}) new_arcs = frozenset({}) visited_nodes = frozenset({}) visited_arcs = frozenset({}) bottleneck_residualCapacity = min( augmentingPath, key=lambda a: a.datum.residualCapacity ).datum.residualCapacity for x in augmentingPath: from_node_uid = x.fromNode.uid if x.datum.action == 'push' else x.toNode.uid to_node_uid = x.toNode.uid if x.datum.action == 'push' else x.fromNode.uid from_node = find_node_by_uid( from_node_uid, H ) to_node = find_node_by_uid( to_node_uid, H ) load = bottleneck_residualCapacity if x.datum.action == 'push' else -bottleneck_residualCapacity for a in H_as_dict[from_node]: if(a.toNode == to_node): test_sum = a.datum.flow + load new_flow = a.datum.flow new_from_node_flow_out = a.fromNode.datum.flowOut new_to_node_flow_in = a.toNode.datum.flowIn new_from_look = {n for n in new_nodes if (n.uid == a.fromNode.uid)} new_to_look = {n for n in new_nodes if (n.uid == a.toNode.uid) } prev_from_node = new_from_look.pop() if (len(new_from_look)>0) else a.fromNode prev_to_node = new_to_look.pop() if (len(new_to_look)>0) else a.toNode new_nodes = new_nodes.difference(frozenset({prev_from_node})) new_nodes = new_nodes.difference(frozenset({prev_to_node})) if(test_sum > a.datum.capacity): new_flow = a.datum.capacity new_from_node_flow_out = prev_from_node.datum.flowOut - a.datum.flow + a.datum.capacity new_to_node_flow_in = prev_to_node.datum.flowIn - a.datum.flow + a.datum.capacity load = test_sum - a.datum.capacity elif(test_sum < 0.0): new_flow = 0.0 new_from_node_flow_out = prev_from_node.datum.flowOut - a.datum.flow new_to_node_flow_in = prev_to_node.datum.flowIn - a.datum.flow load = test_sum else: new_flow = test_sum new_from_node_flow_out = prev_from_node.datum.flowOut - a.datum.flow + new_flow new_to_node_flow_in = prev_to_node.datum.flowIn - a.datum.flow + new_flow load = 0.0 new_from_node_flow_out = round(new_from_node_flow_out, TOL) new_to_node_flow_in = round(new_to_node_flow_in, TOL) new_flow = round(new_flow, TOL) new_from_node = Node(prev_from_node.uid, FlowNodeDatum(prev_from_node.datum.flowIn, new_from_node_flow_out)) new_to_node = Node(prev_to_node.uid, FlowNodeDatum(new_to_node_flow_in, prev_to_node.datum.flowOut)) new_arc = Arc(new_from_node, new_to_node, FlowArcDatum(a.datum.capacity, new_flow)) visited_nodes = visited_nodes.union(frozenset({a.fromNode,a.toNode})) visited_arcs = visited_arcs.union(frozenset({a})) new_nodes = new_nodes.union(frozenset({new_from_node, new_to_node})) new_arcs = new_arcs.union(frozenset({new_arc})) not_visited_nodes = H.setOfNodes.difference(visited_nodes) not_visited_arcs = H.setOfArcs.difference(visited_arcs) full_new_nodes = new_nodes.union(not_visited_nodes) full_new_arcs = new_arcs.union(not_visited_arcs) G = DiGraph(full_new_nodes, full_new_arcs) full_new_arcs_update = frozenset([]) for a in full_new_arcs: new_from_node = a.fromNode new_to_node = a.toNode new_from_node = find_node_by_uid( a.fromNode.uid, G ) new_to_node = find_node_by_uid( a.toNode.uid, G ) full_new_arcs_update = full_new_arcs_update.union( {Arc(new_from_node, new_to_node, FlowArcDatum(a.datum.capacity, a.datum.flow))} ) G = DiGraph(full_new_nodes, full_new_arcs_update) return G

Dans ce qui précède, TOL est une valeur de tolérance pour arrondir les valeurs de débit dans le réseau. Cela permet d'éviter l'imprécision en cascade des calculs en virgule flottante. Ainsi, par exemple, j'ai utilisé TOL = 10 pour signifier arrondir à 10 chiffres significatifs.

Soit K = augment(augmentingPath, H) , alors K représente une solution de débit max faisable pour maxFlowProblem . Pour que la déclaration soit vraie, le réseau de flux représenté par K doit avoir des flux réalisables (ne pas violer la contrainte de capacité ou la contrainte de conservation .

Voici pourquoi : dans la méthode ci-dessus, chaque nœud ajouté au nouveau réseau de flux représenté par le digraphe K est soit une copie exacte d'un nœud du digraphe H , soit un nœud n auquel le même numéro a été ajouté à son n.datum.flowIn comme son n.datum.flowOut . Cela signifie que la contrainte de conservation est satisfaite dans K tant qu'elle est satisfaite dans H . La contrainte de conservation est satisfaite car on vérifie explicitement que tout nouvel arc a dans le réseau a a.datum.flow <= a.datum.capacity ; ainsi, tant que les arcs de l'ensemble H.setOfArcs qui ont été copiés sans modification dans K.setOfArcs ne violent pas la contrainte de capacité , alors K ne viole pas la contrainte de capacité .

Il est également vrai que get_flow_value(H, maxFlowProblem) < get_flow_value(K, maxFlowProblem) si H n'est pas optimal .

Voici pourquoi : pour qu'un chemin d' augmentingPath existe dans la représentation G_f du graphe résiduel G_f d'un problème de flux max maxFlowProblem , le dernier arc a sur augmentingPath doit être un arc "push" et il doit avoir a.toNode == maxFlowProblem.terminalNode . Un chemin d'augmentation est défini comme celui qui se termine au nœud terminal du problème de débit maximal pour lequel il s'agit d'un chemin d'augmentation . D'après la définition du graphe résiduel , il est clair que le dernier arc d'un chemin d'augmentation sur ce graphe résiduel doit être un arc "poussé" car tout arc "tiré" b dans le chemin d'augmentation aura b.toNode == maxFlowProblem.terminalNode et b.fromNode != maxFlowProblem.terminalNode à partir de la définition de path . De plus, à partir de la définition de path , il est clair que le nœud terminal n'est modifié qu'une seule fois par la méthode augment . Ainsi, augment modifie maxFlowProblem.terminalNode.flowIn exactement une fois et augmente la valeur de maxFlowProblem.terminalNode.flowIn car le dernier arc de l' augmentingPath doit être l' arc qui provoque la modification de maxFlowProblem.terminalNode.flowIn pendant augment . D'après la définition d' augment telle qu'elle s'applique aux arcs "push", le maxFlowProblem.terminalNode.flowIn ne peut être qu'augmenté, pas diminué.

Quelques preuves de Sedgewick et Wayne

Le livre Algorithms, quatrième édition de Robert Sedgewich et Kevin Wayne contient de merveilleuses et courtes preuves (pages 892-894) qui seront utiles. Je vais les recréer ici, même si j'utiliserai un langage correspondant aux définitions précédentes. Mes étiquettes pour les épreuves sont les mêmes que dans le livre de Sedgewick.

Proposition E : Pour tout digraphe H représentant une solution de débit maximum réalisable à un problème de débit maximum maxFlowProblem , pour tout stCut get_stcut_flow(stCut,H,maxFlowProblem) = get_flow_value(H, maxFlowProblem) .

Preuve : Soit stCut=stCut(maxFlowProblem.sourceNode,maxFlowProblem.terminalNode,set([a for a in H.setOfArcs if a.toNode == maxFlowProblem.terminalNode])) . La proposition E est valable pour stCut directement à partir de la définition de st valeur de flux . Supposons que nous souhaitions déplacer un nœud n de la s-partition ( get_first_part(stCut.cut, G) ) et dans la t-partition (get_second_part(stCut.cut, G)) , pour ce faire, nous devons changer stCut.cut , ce qui pourrait changer stcut_flow = get_stcut_flow(stCut,H,maxFlowProblem) et invalider la proposition E . Cependant, voyons comment la valeur de stcut_flow changera au fur et à mesure que nous apporterons ce changement. le nœud n est à l'équilibre, ce qui signifie que la somme des flux entrant dans le nœud n est égale à la somme des flux sortant de celui-ci (cela est nécessaire pour que H représente une solution réalisable ). Notez que tous les flux qui font partie du nœud n entrant dans stcut_flow entrent depuis la partition s (le flux entrant dans le nœud n depuis la partition t, directement ou indirectement, n'aurait pas été compté dans la valeur stcut_flow car il se dirige dans la mauvaise direction selon la définition). De plus, tout flux sortant de n finira par s'écouler (directement ou indirectement) dans le nœud terminal (prouvé précédemment). Lorsque nous déplaçons le nœud n dans la t-partition, tout le flux entrant dans n depuis la s-partition doit être ajouté à la nouvelle valeur stcut_flow ; cependant, tous les flux sortant de n doivent être soustraits de la nouvelle valeur stcut_flow ; la partie du flux se dirigeant directement vers la t-partition est soustraite car ce flux est désormais interne à la nouvelle t-partition et n'est pas compté comme stcut_flow . La partie du flux de n se dirigeant vers les nœuds de la s-partition doit également être soustraite de stcut_flow : Après que n soit déplacé dans la t-partition, ces flux seront dirigés de la t-partition vers la s-partition et ainsi ne doit pas être pris en compte dans stcut_flow , puisque ces flux sont supprimés de l'afflux dans la partition s et doivent être réduits par la somme de ces flux, et le flux sortant de la partition s vers la partition t (où tous les flux de st doit finir) doit être réduit du même montant. Comme le nœud n était à l'équilibre avant le processus, la mise à jour aura ajouté la même valeur à la nouvelle valeur stcut_flow qu'elle a soustraite, laissant ainsi la proposition E vraie après la mise à jour. La validité de la proposition E découle alors de l'induction sur la taille de la t-partition.

Voici quelques exemples de réseaux de flux pour aider à visualiser les cas moins évidents où la proposition E est valable ; dans l'image, les zones rouges indiquent la partition s, les zones bleues représentent la partition t et les arcs verts indiquent une ère coupe . Dans la deuxième image, le flux entre le nœud A et le nœud B augmente alors que le flux vers le nœud terminal t ne change pas :

Corollaire : Aucune valeur de débit de st cut ne peut dépasser la capacité de n'importe quel st cut .

Proposition F. (max flow, min cut theorem) : Soit f un er flux . Les 3 conditions suivantes sont équivalentes :

1. Il existe une ère coupe dont la capacité est égale à la valeur du débit f .

2. f est un débit max .

3. Il n'y a pas de chemin augmentant par rapport à f .

La condition 1 implique la condition 2 par le corollaire. La condition 2 implique la condition 3 car l'existence d'un chemin augmentant implique l'existence d'un flux avec des valeurs plus grandes, contredisant la maximalité de f . La condition 3 implique la condition 1 : Soit C_s l'ensemble de tous les nœuds accessibles depuis s avec un chemin croissant dans le graphe résiduel . Soit C_t les arcs restants , alors t doit être dans C_t (par notre hypothèse). Les arcs passant de C_s à C_t forment alors une ère coupe qui ne contient que des arcs a où soit a.datum.capacity = a.datum.flow soit a.datum.flow = 0.0 . S'il en était autrement, alors les nœuds connectés par un arc avec une capacité résiduelle restante à C_s seraient dans l'ensemble C_s puisqu'il y aurait alors un chemin augmentant de s à un tel nœud . Le flux à travers la ère coupe est égal à la capacité de la ère coupe (puisque les arcs de C_s à C_t ont un flux égal à la capacité) et aussi à la valeur du er flux (par la proposition E ).

Cette déclaration du débit max, théorème de coupe min implique la déclaration précédente de Flows in Networks.

Corollaire (propriété d'intégralité) : lorsque les capacités sont des nombres entiers, il existe un flux maximal de valeurs entières, et l'algorithme de Ford-Fulkerson le trouve.

Preuve : Chaque chemin augmentant augmente le débit d'un entier positif, le minimum des capacités inutilisées dans les arcs 'push' et les flux dans les arcs 'pull', qui sont toujours des entiers positifs.

Cela justifie la description de la méthode Ford-Fulkerson du CLRS . La méthode consiste à continuer à trouver des chemins d'augmentation et à appliquer l' augment à la dernière maxFlowSolution proposant de meilleures solutions, jusqu'à ce qu'il n'y ait plus de chemin d'augmentation, ce qui signifie que la dernière solution de débit maximal est optimale .

De Ford-Fulkerson à Edmonds-Karp

Les questions restantes concernant la résolution des problèmes de débit maximum sont :

1. Comment construire les chemins augmentants ?

2. La méthode se terminera-t-elle si nous utilisons des nombres réels et non des entiers ?

3. Combien de temps faudra-t-il pour résilier (si c'est le cas) ?

L' algorithme d'Edmonds-Karp spécifie que chaque chemin augmentant est construit par une recherche en largeur (BFS) du graphe résiduel ; il s'avère que cette décision du point 1 ci-dessus va aussi forcer l'algorithme à se terminer (point 2) et permet de déterminer la complexité asymptotique en temps et en espace.

Tout d'abord, voici une implémentation BFS :

 def bfs(sourceNode, terminalNode, G_f): G_f_as_dict = digraph_to_dict(G_f) parent_arcs = dict([]) visited = frozenset([]) deq = deque([sourceNode]) while len(deq) > 0: curr = deq.popleft() if curr == terminalNode: break for a in G_f_as_dict.get(curr): if (a.toNode not in visited): visited = visited.union(frozenset([a.toNode])) parent_arcs[a.toNode] = a deq.extend([a.toNode]) path = list([]) curr = terminalNode while(curr != sourceNode): if (curr not in parent_arcs): err_msg = 'No augmenting path from {} to {}.'.format(sourceNode.uid, terminalNode.uid) raise StopIteration(err_msg) path.extend([parent_arcs[curr]]) curr = parent_arcs[curr].fromNode path.reverse() test = deque([path[0].fromNode]) for a in path: if(test[-1] != a.fromNode): err_msg = 'Bad path: {}'.format(path) raise ValueError(err_msg) test.extend([a.toNode]) return path

J'ai utilisé un deque du module python collections.

Pour répondre à la question 2 ci-dessus, je vais paraphraser une autre preuve de Sedgewick et Wayne : Proposition G. Le nombre de chemins d'augmentation nécessaires dans l'algorithme d' Edmonds-Karp avec N nœuds et A arcs est au plus NA/2 . Preuve : Chaque chemin d'augmentation a un arc de goulot d'étranglement - un arc qui est supprimé du graphe résiduel car il correspond soit à un arc "poussant" qui se remplit à pleine capacité, soit à un arc "tirant" à travers lequel le flux devient 0. Chaque fois qu'un arc devient un arc de goulot d'étranglement, la longueur de tout chemin d'augmentation à travers celui-ci doit augmenter d'un facteur 2. En effet, chaque nœud d'un chemin peut n'apparaître qu'une seule fois ou pas du tout (d'après la définition de path ) puisque les chemins sont en train d'être exploré du chemin le plus court au plus long, cela signifie qu'au moins un nœud supplémentaire doit être visité par le chemin suivant qui passe par le nœud goulot d'étranglement particulier, ce qui signifie 2 arcs supplémentaires sur le chemin avant d'arriver au nœud . Puisque le chemin d'augmentation est de longueur au plus N , chaque arc peut être sur au plus N/2 chemins d'augmentation et le nombre total de chemins d'augmentation est au plus NA/2 .

L'algorithme d' Edmonds-Karp s'exécute en O(NA^2) . Si au plus NA/2 chemins seront explorés pendant l'algorithme et que l'exploration de chaque chemin avec BFS est N+A alors le terme le plus significatif du produit et donc la complexité asymptotique est O(NA^2) .

Soit mfp un maxFlowProblemState .

 def edmonds_karp(mfp): sid, tid = mfp.sourceNodeUid, mfp.terminalNodeUid no_more_paths = False loop_count = 0 while(not no_more_paths): residual_digraph = get_residual_graph_of(mfp.G) try: asource = find_node_by_uid(mfp.sourceNodeUid, residual_digraph) aterm = find_node_by_uid(mfp.terminalNodeUid, residual_digraph) apath = bfs(asource, aterm, residual_digraph) paint_mfp_path(mfp, loop_count, apath) G = augment(apath, mfp.G) s = find_node_by_uid(sid, G) t = find_node_by_uid(tid, G) mfp = MaxFlowProblemState(G, s.uid, t.uid, mfp.maxFlowProblemStateUid) loop_count += 1 except StopIteration as e: no_more_paths = True return mfp

La version ci-dessus est inefficace et a une complexité pire que O(NA^2) car elle construit une nouvelle solution de débit maximum et un nouveau graphique résiduel à chaque fois (plutôt que de modifier les digraphes existants au fur et à mesure que l'algorithme avance). Pour arriver à une vraie solution O(NA^2) , l'algorithme doit maintenir à la fois le digraphe représentant l' état du problème de débit maximum et son graphe résiduel associé. Ainsi, l'algorithme doit éviter d'itérer inutilement sur les arcs et les nœuds et mettre à jour leurs valeurs et les valeurs associées dans le graphe résiduel uniquement si nécessaire.

Pour écrire un algorithme d' Edmonds Karp plus rapide, j'ai réécrit plusieurs morceaux de code de ce qui précède. J'espère que parcourir le code qui a généré un nouveau digraphe a été utile pour comprendre ce qui se passe. Dans l'algorithme rapide, j'utilise de nouvelles astuces et des structures de données Python que je ne veux pas détailler. Je dirai que a.fromNode et a.toNode sont maintenant traités comme des chaînes et des uids aux nœuds . Pour ce code, laissez mfps être un maxFlowProblemState

 import uuid def fast_get_mfps_flow(mfps): flow_from_s = {n for n in mfps.G.setOfNodes if n.uid == mfps.sourceNodeUid}.pop().datum.flowOut flow_to_t = {n for n in mfps.G.setOfNodes if n.uid == mfps.terminalNodeUid}.pop().datum.flowIn if( (flow_to_t - flow_from_s) > 0.00): raise Exception('Infeasible st flow') return flow_to_t def fast_e_k_preprocess(G): G = strip_flows(G) get = dict({}) get['nodes'] = dict({}) get['node_to_node_capacity'] = dict({}) get['node_to_node_flow'] = dict({}) get['arcs'] = dict({}) get['residual_arcs'] = dict({}) for a in G.setOfArcs: if(a.fromNode not in G.setOfNodes): err_msg = 'There is no Node {a.fromNode.uid!s} to match the Arc from {a.fromNode.uid!s} to {a.toNode.uid!s}'.format(**locals()) raise KeyError(err_msg) if(a.toNode not in G.setOfNodes): err_msg = 'There is no Node {a.toNode.uid!s} to match the Arc from {a.fromNode.uid!s} to {a.toNode.uid!s}'.format(**locals()) raise KeyError(err_msg) get['nodes'][a.fromNode.uid] = a.fromNode get['nodes'][a.toNode.uid] = a.toNode lark = Arc(a.fromNode.uid, a.toNode.uid, FlowArcDatumWithUid(a.datum.capacity, a.datum.flow, uuid.uuid4())) if(a.fromNode.uid not in get['arcs']): get['arcs'][a.fromNode.uid] = dict({a.toNode.uid : dict({lark.datum.uid : lark})}) else: if(a.toNode.uid not in get['arcs'][a.fromNode.uid]): get['arcs'][a.fromNode.uid][a.toNode.uid] = dict({lark.datum.uid : lark}) else: get['arcs'][a.fromNode.uid][a.toNode.uid][lark.datum.uid] = lark for a in G.setOfArcs: if a.toNode.uid not in get['arcs']: get['arcs'][a.toNode.uid] = dict({}) for n in get['nodes']: get['residual_arcs'][n] = dict() get['node_to_node_capacity'][n] = dict() get['node_to_node_flow'][n] = dict() for u in get['nodes']: n_to_u_cap_sum = sum(a.datum.capacity for a in G.setOfArcs if (a.fromNode.uid == n) and (a.toNode.uid == u) ) n_to_u_flow_sum = sum(a.datum.flow for a in G.setOfArcs if (a.fromNode.uid == n) and (a.toNode.uid == u) ) if(n_to_u_cap_sum > n_to_u_flow_sum): flow = round(n_to_u_cap_sum - n_to_u_flow_sum, TOL) get['residual_arcs'][n][u] = Arc(n,u,ResidualDatum(flow, 'push')) if(n_to_u_flow_sum > 0.0): flow = round(n_to_u_flow_sum, TOL) get['residual_arcs'][u][n] = Arc(u,n,ResidualDatum(flow, 'pull')) get['node_to_node_capacity'][n][u] = n_to_u_cap_sum get['node_to_node_flow'][n][u] = n_to_u_flow_sum return get def fast_bfs(sid, tid, get): parent_of = dict([]) visited = frozenset([]) deq = coll.deque([sid]) while len(deq) > 0: n = deq.popleft() if n == tid: break for u in get['residual_arcs'][n]: if (u not in visited): visited = visited.union(frozenset({u})) parent_of[u] = get['residual_arcs'][n][u] deq.extend([u]) path = list([]) curr = tid while(curr != sid): if (curr not in parent_of): err_msg = 'No augmenting path from {} to {}.'.format(sid, curr) raise StopIteration(err_msg) path.extend([parent_of[curr]]) curr = parent_of[curr].fromNode path.reverse() return path def fast_edmonds_karp(mfps): sid, tid = mfps.sourceNodeUid, mfps.terminalNodeUid get = fast_e_k_preprocess(mfps.G) no_more_paths, loop_count = False, 0 while(not no_more_paths): try: apath = fast_bfs(sid, tid, get) get = fast_augment(apath, get) loop_count += 1 except StopIteration as e: no_more_paths = True nodes = frozenset(get['nodes'].values()) arcs = frozenset({}) for from_node in get['arcs']: for to_node in get['arcs'][from_node]: for arc in get['arcs'][from_node][to_node]: arcs |= frozenset({get['arcs'][from_node][to_node][arc]}) G = DiGraph(nodes, arcs) mfps = MaxFlowProblemState(G, sourceNodeUid=sid, terminalNodeUid=tid, maxFlowProblemStateUid=mfps.maxFlowProblemStateUid) return mfps

Voici une visualisation de la façon dont cet algorithme résout l'exemple de réseau de flux d'en haut. La visualisation montre les étapes telles qu'elles sont reflétées dans le digramme G représentant le réseau de flux le plus à jour et telles qu'elles sont reflétées dans le graphique résiduel de ce réseau de flux. Les chemins d'augmentation dans le graphique résiduel sont représentés par des chemins rouges, et le digraphe représentant le problème de l'ensemble de nœuds et d' arcs affectés par un chemin d'augmentation donné est surligné en vert. Dans chaque cas, je mettrai en surbrillance les parties du graphique qui seront modifiées (en rouge ou en vert), puis je montrerai le graphique après les modifications (juste en noir).

Voici une autre visualisation de la façon dont cet algorithme résout un autre exemple de réseau de flux . Notez que celui-ci utilise des nombres réels et contient plusieurs arcs avec les mêmes valeurs fromNode et toNode .

** Notez également que parce que les arcs avec un ResidualDatum "pull" peuvent faire partie du chemin d'augmentation, les nœuds affectés dans le DiGraph du réseau Flown _ peuvent ne pas être sur un chemin dans G! .

Graphiques bipartis

Supposons que nous ayons un digraphe G , G est biparti s'il est possible de partitionner les nœuds de G.setOfNodes en deux ensembles ( part_1 et part_2 ) tels que pour tout arc a dans G.setOfArcs , il ne peut pas être vrai que a.fromNode dans part_1 et a.toNode dans part_1 . Il ne peut pas non plus être vrai que a.fromNode dans part_2 et a.toNode dans part_2 .

En d'autres termes, G est biparti s'il peut être partitionné en deux ensembles de nœuds de sorte que chaque arc doit connecter un nœud d'un ensemble à un nœud de l'autre ensemble.

Test bipartite

Supposons que nous ayons un digraphe G , nous voulons tester s'il est biparti . Nous pouvons le faire en O(|G.setOfNodes|+|G.setOfArcs|) en colorant le graphe en deux couleurs.

Premièrement, nous devons générer un nouveau digraphe H . Ce graphe aura aura le même ensemble de nœuds que G , mais il aura plus d' arcs que G . Chaque arc a dans G créera 2 arcs dans H ; le premier arc sera identique à a , et le second arc inversera le directeur de a ( b = Arc(a.toNode,a.fromNode,a.datum) ).

 Bipartition = coll.namedtuple('Bipartition',['firstPart', 'secondPart', 'G']) def bipartition(G): nodes = frozenset(G.setOfNodes arcs = frozenset(G.setOfArcs) arcs = arcs.union( frozenset( {Arc(a.toNode,a.fromNode,a.datum) for a in G.setOfArcs} ) ) H = DiGraph(nodes, arcs) H_as_dict = digraph_to_dict(H) color = dict([]) some_node = next(iter(H.setOfNodes)) deq = coll.deque([some_node]) color[some_node] = -1 while len(deq) > 0: curr = deq.popleft() for a in H_as_dict.get(curr): if (a.toNode not in color): color[a.toNode] = -1*color[curr] deq.extend([a.toNode]) elif(color[curr] == color[a.toNode]): print(curr,a.toNode) raise Exception('Not Bipartite.') firstPart = frozenset( {n for n in color if color[n] == -1 } ) secondPart = frozenset( {n for n in color if color[n] == 1 } ) if( firstPart.union(secondPart) != G.setOfNodes ): raise Exception('Not Bipartite.') return Bipartition(firstPart, secondPart, G)

Correspondances et correspondances maximales

Supposons que nous ayons un digraphe G et que matching soit un sous-ensemble d' arcs de G.setOfArcs . matching est une correspondance si pour deux arcs a et b dans matching : len(frozenset( {a.fromNode} ).union( {a.toNode} ).union( {b.fromNode} ).union( {b.toNode} )) == 4 . En d'autres termes, deux arcs dans une correspondance ne partagent pas un nœud .

Matching matching , est une correspondance maximale s'il n'y a pas d'autre alt_matching correspondant dans G tel que len(matching) < len(alt_matching) . En d'autres termes, matching est une correspondance maximale si c'est le plus grand ensemble d' arcs de G.setOfArcs qui satisfait toujours la définition de la correspondance (l'ajout de tout arc qui n'est pas déjà dans la correspondance cassera la définition de la correspondance ).

Une matching maximale est une correspondance parfaite si pour chaque nœud n dans G.setOfArcs il existe un arc a dans la matching où a.fromNode == n or a.toNode == n .

Correspondance bipartite maximale

Un couplage bipartite maximum est un couplage maximum sur un digraphe G qui est biparti .

Étant donné que G est bipartite , le problème de trouver une correspondance bipartite maximale peut être transformé en un problème de débit maximal résoluble avec l'algorithme d' Edmonds-Karp , puis la correspondance bipartite maximale peut être récupérée à partir de la solution du problème de débit maximal .

Soit bipartition une bipartition de G .

Pour ce faire, je dois générer un nouveau digraphe ( H ) avec de nouveaux nœuds ( H.setOfNodes ) et de nouveaux arcs ( H.setOfArcs ). H.setOfNodes contient tous les nœuds de G.setOfNodes et deux autres nœuds , s (un nœud source ) et t (un nœud terminal ).

H.setOfArcs contiendra un arc pour chaque G.setOfArcs . Si un arc a est dans G.setOfArcs et a.fromNode est dans bipartition.firstPart et a.toNode est dans bipartition.secondPart alors incluez a dans H (en ajoutant un FlowArcDatum(1,0) ).

Si a.fromNode est dans bipartition.secondPart et a.toNode est dans bipartition.firstPart , alors incluez Arc(a.toNode,a.fromNode,FlowArcDatum(1,0)) dans H.setOfArcs .

La définition d'un graphe bipartite garantit qu'aucun arc ne relie les nœuds où les deux nœuds se trouvent dans la même partition. H.setOfArcs contient également un arc du nœud s à chaque nœud dans bipartition.firstPart . Enfin, H.setOfArcs contient un arc de chaque nœud dans bipartition.secondPart au nœud t . a.datum.capacity = 1 pour tout a dans H.setOfArcs .

Commencez par partitionner les nœuds dans G.setOfNodes en deux ensembles disjoints ( part1 et part2 ) de sorte qu'aucun arc dans G.setOfArcs ne soit dirigé d'un ensemble vers le même ensemble (cette partition est possible car G est bipartite ). Ensuite, ajoutez tous les arcs dans G.setOfArcs qui sont dirigés de part1 à part2 dans H.setOfArcs . Créez ensuite un seul nœud source s et un seul nœud terminal t et créez d'autres arcs

Ensuite, construisez un maxFlowProblemState .

 def solve_mbm( bipartition ): s = Node(uuid.uuid4(), FlowNodeDatum(0,0)) t = Node(uuid.uuid4(), FlowNodeDatum(0,0)) translate = {} arcs = frozenset([]) for a in bipartition.G.setOfArcs: if ( (a.fromNode in bipartition.firstPart) and (a.toNode in bipartition.secondPart) ): fark = Arc(a.fromNode,a.toNode,FlowArcDatum(1,0)) arcs = arcs.union({fark}) translate[frozenset({a.fromNode.uid,a.toNode.uid})] = a elif ( (a.toNode in bipartition.firstPart) and (a.fromNode in bipartition.secondPart) ): bark = Arc(a.toNode,a.fromNode,FlowArcDatum(1,0)) arcs = arcs.union({bark}) translate[frozenset({a.fromNode.uid,a.toNode.uid})] = a arcs1 = frozenset( {Arc(s,n,FlowArcDatum(1,0)) for n in bipartition.firstPart } ) arcs2 = frozenset( {Arc(n,t,FlowArcDatum(1,0)) for n in bipartition.secondPart } ) arcs = arcs.union(arcs1.union(arcs2)) nodes = frozenset( {Node(n.uid,FlowNodeDatum(0,0)) for n in bipartition.G.setOfNodes} ).union({s}).union({t}) G = DiGraph(nodes, arcs) mfp = MaxFlowProblemState(G, s.uid, t.uid, 'bipartite') result = edmonds_karp(mfp) lookup_set = {a for a in result.G.setOfArcs if (a.datum.flow > 0) and (a.fromNode.uid != s.uid) and (a.toNode.uid != t.uid)} matching = {translate[frozenset({a.fromNode.uid,a.toNode.uid})] for a in lookup_set} return matching

Couverture minimale des nœuds

Une couverture de nœuds dans un digraphe G est un ensemble de nœuds ( cover ) de G.setOfNodes tel que pour tout arc a de G.setOfArcs cela doit être vrai : (a.fromNode in cover) or (a.toNode in cover) .

Une couverture minimale de nœuds est le plus petit ensemble possible de nœuds dans le graphe qui est toujours une couverture de nœuds . Le théorème de Konig stipule que dans un graphe bipartite , la taille de la correspondance maximale sur ce graphe est égale à la taille de la couverture de nœuds minimale , et il suggère comment la couverture de nœuds peut être récupérée à partir d'une correspondance maximale :

Supposons que nous ayons la bipartition bipartition et la correspondance maximale matching . Définissez un nouveau digraphe H , H.setOfNodes=G.setOfNodes , les arcs dans H.setOfArcs sont l'union de deux ensembles.

Le premier ensemble est arcs a in matching , avec le changement que si a.fromNode in bipartition.firstPart et a.toNode in bipartition.secondPart alors a.fromNode et a.toNode sont échangés dans l' arc créé donnent à ces arcs un a.datum.inMatching=True attribut pour indiquer qu'ils ont été dérivés d' arcs dans un fichier .

Le deuxième ensemble est arcs a NOT in matching , avec le changement que si a.fromNode in bipartition.secondPart et a.toNode in bipartition.firstPart alors a.fromNode et a.toNode sont échangés dans l' arc créé (donnez à ces arcs un a.datum.inMatching=False attribut).

Ensuite, exécutez une première recherche en profondeur (DFS) à partir de chaque nœud n dans bipartition.firstPart qui n'est ni n == a.fromNode ni n == a.toNodes pour tout arc a dans matching . Pendant le DFS, certains nœuds sont visités et d'autres non (stockez ces informations dans un champ n.datum.visited ). La couverture minimale des nœuds est l'union des nœuds {a.fromNode for a in H.setOfArcs if ( (a.datum.inMatching) and (a.fromNode.datum.visited) ) } et des nœuds {a.fromNode for a in H.setOfArcs if (a.datum.inMatching) and (not a.toNode.datum.visited)} .

On peut montrer que cela conduit d'une correspondance maximale à une couverture minimale de nœuds par une preuve par contradiction, prenez un arc a qui n'était pas censé être couvert et considérez les quatre cas concernant l'appartenance de a.fromNode et a.toNode (que ce soit en tant que toNode ou fromNode ) à n'importe quel arc dans matching matching . Chaque cas conduit à une contradiction due à l'ordre dans lequel DFS visite les nœuds et au fait que la matching est une correspondance maximale .

Supposons que nous ayons une fonction pour exécuter ces étapes et renvoyer l'ensemble de nœuds comprenant la couverture de nœuds minimale lorsqu'on lui donne le digraphe G , et la correspondance maximale matching :

 ArcMatchingDatum = coll.namedtuple('ArcMatchingDatum', ['inMatching' ]) NodeMatchingDatum = coll.namedtuple('NodeMatchingDatum', ['visited']) def dfs_helper(snodes, G): sniter, do_stop = iter(snodes), False visited, visited_uids = set(), set() while(not do_stop): try: stack = [ next(sniter) ] while len(stack) > 0: curr = stack.pop() if curr.uid not in visited_uids: visited = visited.union( frozenset( {Node(curr.uid, NodeMatchingDatum(False))} ) ) visited_uids = visited.union(frozenset({curr.uid})) succ = frozenset({a.toNode for a in G.setOfArcs if a.fromNode == curr}) stack.extend( succ.difference( frozenset(visited) ) ) except StopIteration as e: stack, do_stop = [], True return visited def get_min_node_cover(matching, bipartition): nodes = frozenset( { Node(n.uid, NodeMatchingDatum(False)) for n in bipartition.G.setOfNodes} ) G = DiGraph(nodes, None) charcs = frozenset( {a for a in matching if ( (a.fromNode in bipartition.firstPart) and (a.toNode in bipartition.secondPart) )} ) arcs0 = frozenset( { Arc(find_node_by_uid(a.toNode.uid,G), find_node_by_uid(a.fromNode.uid,G), ArcMatchingDatum(True) ) for a in charcs } ) arcs1 = frozenset( { Arc(find_node_by_uid(a.fromNode.uid,G), find_node_by_uid(a.toNode.uid,G), ArcMatchingDatum(True) ) for a in matching.difference(charcs) } ) not_matching = bipartition.G.setOfArcs.difference( matching ) charcs = frozenset( {a for a in not_matching if ( (a.fromNode in bipartition.secondPart) and (a.toNode in bipartition.firstPart) )} ) arcs2 = frozenset( { Arc(find_node_by_uid(a.toNode.uid,G), find_node_by_uid(a.fromNode.uid,G), ArcMatchingDatum(False) ) for a in charcs } ) arcs3 = frozenset( { Arc(find_node_by_uid(a.fromNode.uid,G), find_node_by_uid(a.toNode.uid,G), ArcMatchingDatum(False) ) for a in not_matching.difference(charcs) } ) arcs = arcs0.union(arcs1.union(arcs2.union(arcs3))) G = DiGraph(nodes, arcs) bip = Bipartition({find_node_by_uid(n.uid,G) for n in bipartition.firstPart},{find_node_by_uid(n.uid,G) for n in bipartition.secondPart},G) match_from_nodes = frozenset({find_node_by_uid(a.fromNode.uid,G) for a in matching}) match_to_nodes = frozenset({find_node_by_uid(a.toNode.uid,G) for a in matching}) snodes = bip.firstPart.difference(match_from_nodes).difference(match_to_nodes) visited_nodes = dfs_helper(snodes, bip.G) not_visited_nodes = bip.G.setOfNodes.difference(visited_nodes) H = DiGraph(visited_nodes.union(not_visited_nodes), arcs) cover1 = frozenset( {a.fromNode for a in H.setOfArcs if ( (a.datum.inMatching) and (a.fromNode.datum.visited) ) } ) cover2 = frozenset( {a.fromNode for a in H.setOfArcs if ( (a.datum.inMatching) and (not a.toNode.datum.visited) ) } ) min_cover_nodes = cover1.union(cover2) true_min_cover_nodes = frozenset({find_node_by_uid(n.uid, bipartition.G) for n in min_cover_nodes}) return min_cover_nodes

Le problème d'affectation linéaire

Le problème d'affectation linéaire consiste à trouver un appariement de poids maximum dans un graphe biparti pondéré.

Des problèmes comme celui au tout début de cet article peuvent être exprimés comme un problème d'affectation linéaire . Étant donné un ensemble de travailleurs, un ensemble de tâches et une fonction indiquant la rentabilité d'une affectation d'un travailleur à une tâche, nous voulons maximiser la somme de toutes les affectations que nous effectuons ; c'est un problème d'affectation linéaire .

Supposons que le nombre de tâches et de travailleurs soient égaux, bien que je montre que cette hypothèse est facile à supprimer. Dans l'implémentation, je représente les poids d'arc avec un attribut a.datum.weight pour un arc a .

 WeightArcDatum = namedtuple('WeightArcDatum', [weight])

Algorithme de Kuhn-Munkres

L'algorithme de Kuhn-Munkres résout le problème d'affectation linéaire . Une bonne implémentation peut prendre O(N^{4}) temps, (où N est le nombre de nœuds dans le digraphe représentant le problème). Une implémentation plus facile à expliquer prend O(N^{5}) (pour une version qui régénère DiGraphs ) et O(N^{4}) for (pour une version qui maintient DiGraphs ). Ceci est similaire aux deux implémentations différentes de l'algorithme d' Edmonds-Karp .

Pour cette description, je ne travaille qu'avec des graphes bipartis complets (ceux où la moitié des nœuds sont dans une partie de la bipartition et l'autre moitié dans la seconde partie). Chez le travailleur, la motivation à la tâche signifie qu'il y a autant de travailleurs que de tâches.

Cela semble être une condition importante (et si ces ensembles ne sont pas égaux !) mais il est facile de résoudre ce problème ; Je parle de la façon de le faire dans la dernière section.

La version de l'algorithme décrite ici utilise le concept utile des arcs de poids nul . Malheureusement, ce concept n'a de sens que lorsque nous résolvons une minimisation (si, plutôt que de maximiser les bénéfices de nos affectations de tâches, nous minimisions plutôt le coût de ces affectations).

Heureusement, il est facile de transformer un problème d'affectation linéaire maximum en un problème d'affectation linéaire minimum en définissant chacun des poids de l' arc a sur Ma.datum.weight où M=max({a.datum.weight for a in G.setOfArcs}) . La solution au problème de maximisation d'origine sera identique au problème de minimisation de la solution après la modification des poids d' arc . Donc, pour le reste, supposons que nous apportions ce changement.

L' algorithme de Kuhn-Munkres résout l'appariement de poids minimum dans un graphe biparti pondéré par une séquence d' appariements maximum dans des graphes bipartis non pondérés. Si a nous trouvons une correspondance parfaite sur la représentation digraphique du problème d'affectation linéaire , et si le poids de chaque arc dans la correspondance est nul, alors nous avons trouvé la correspondance de poids minimum puisque cette correspondance suggère que tous les nœuds du digraphe ont été correspondant à un arc avec le coût le plus bas possible (aucun coût ne peut être inférieur à 0, sur la base des définitions précédentes).

Aucun autre arc ne peut être ajouté à la correspondance (car tous les nœuds sont déjà mis en correspondance) et aucun arc ne doit être supprimé de la correspondance car tout arc de remplacement possible aura une valeur de poids au moins aussi élevée.

Si nous trouvons une correspondance maximale du sous-graphe de G qui ne contient que des arcs de poids nul , et que ce n'est pas une correspondance parfaite , nous n'avons pas de solution complète (puisque la correspondance n'est pas parfaite ). Cependant, nous pouvons produire un nouveau digraphe H en modifiant les poids des arcs dans G.setOfArcs de manière à ce que de nouveaux arcs de poids 0 apparaissent et que la solution optimale de H soit la même que la solution optimale de G . Puisque nous garantissons qu'au moins un arc de poids nul est produit à chaque itération, nous garantissons que nous arriverons à une correspondance parfaite en pas plus de |G.setOfNodes|^{2}=N^{2} de telles itérations.

Supposons que dans bipartition bipartition , bipartition.firstPart contient des nœuds représentant des travailleurs et bipartition.secondPart représente des nœuds représentant des tâches.

L'algorithme commence par générer un nouveau digraphe H . H.setOfNodes = G.setOfNodes . Certains arcs dans H sont générés à partir des nœuds n dans bipartition.firstPart . Chacun de ces nœuds n génère un arc b dans H.setOfArcs pour chaque arc a dans bipartition.G.setOfArcs où a.fromNode = n or a.toNode = n , b=Arc(a.fromNode, a.toNode, a.datum.weight - z) où z=min(x.datum.weight for x in G.setOfArcs if ( (x.fromNode == n) or (x.toNode == n) )) .

Plus d' arcs dans H sont générés à partir des nœuds n dans bipartition.secondPart . Chacun de ces nœuds n génère un arc b dans H.setOfArcs pour chaque arc a dans bipartition.G.setOfArcs où a.fromNode = n ou a.toNode = n , b=Arc(a.fromNode, a.toNode, ArcWeightDatum(a.datum.weight - z)) où z=min(x.datum.weight for x in G.setOfArcs if ( (x.fromNode == n) or (x.toNode == n) )) .

KMA : Ensuite, formez un nouveau digraphe K composé uniquement des arcs de poids nul de H , et des nœuds incidents sur ces arcs . Formez une bipartition sur les nœuds de K , puis utilisez solve_mbm( bipartition ) pour obtenir un matching maximum ( matching ) sur K . Si matching est une correspondance parfaite dans H (les arcs de matching sont incidents sur tous les nœuds de H.setOfNodes ), alors la matching est une solution optimale au problème d'affectation linéaire .

Sinon, si matching n'est pas parfaite , générez la couverture de nœud minimale de K en utilisant node_cover = get_min_node_cover(matching, bipartition(K)) . Ensuite, définissez z=min({a.datum.weight for a in H.setOfArcs if a not in node_cover}) . Définir nodes = H.setOfNodes , arcs1 = {Arc(a.fromNode,a.toNode,ArcWeightDatum(a.datum.weigh-z)) for a in H.setOfArcs if ( (a.fromNode not in node_cover) and (a.toNode not in node_cover)} , arcs2 = {Arc(a.fromNode,a.toNode,ArcWeightDatum(a.datum.weigh)) for a in H.setOfArcs if ( (a.fromNode not in node_cover) != (a.toNode not in node_cover)} , arcs3 = {Arc(a.fromNode,a.toNode,ArcWeightDatum(a.datum.weigh+z)) for a in H.setOfArcs if ( (a.fromNode in node_cover) and (a.toNode in node_cover)} . Le symbole != dans l'expression précédente agit comme un opérateur XOR. Alors arcs = arcs1.union(arcs2.union(arcs3)) . Ensuite, H=DiGraph(nodes,arcs) . Revenir à l'étiquette KMA L'algorithme continue jusqu'à ce qu'un appariement parfait soit produit Cet appariement est aussi la solution du problème d'affectation linéaire .

 def kuhn_munkres( bipartition ): nodes = bipartition.G.setOfNodes arcs = frozenset({}) for n in bipartition.firstPart: z = min( {x.datum.weight for x in bipartition.G.setOfArcs if ( (x.fromNode == n) or (x.toNode == n) )} ) arcs = arcs.union( frozenset({Arc(a.fromNode, a.toNode, ArcWeightDatum(a.datum.weight - z)) }) ) for n in bipartition.secondPart: z = min( {x.datum.weight for x in bipartition.G.setOfArcs if ( (x.fromNode == n) or (x.toNode == n) )} ) arcs = arcs.union( frozenset({Arc(a.fromNode, a.toNode, ArcWeightDatum(a.datum.weight - z)) }) ) H = DiGraph(nodes, arcs) assignment, value = dict({}), 0 not_done = True while( not_done ): zwarcs = frozenset( {a for a in H.setOfArcs if a.datum.weight == 0} ) znodes = frozenset( {n.fromNode for n in zwarcs} ).union( frozenset( {n.toNode for n in zwarcs} ) ) K = DiGraph(znodes, zwarcs) k_bipartition = bipartition(K) matching = solve_mbm( k_bipartition ) mnodes = frozenset({a.fromNode for a in matching}).union(frozenset({a.toNode for a in matching})) if( len(mnodes) == len(H.setOfNodes) ): for a in matching: assignment[ a.fromNode.uid ] = a.toNode.uid value = sum({a.datum.weight for a in matching}) not_done = False else: node_cover = get_min_node_cover(matching, bipartition(K)) z = min( frozenset( {a.datum.weight for a in H.setOfArcs if a not in node_cover} ) ) nodes = H.setOfNodes arcs1 = frozenset( {Arc(a.fromNode,a.toNode,ArcWeightDatum(a.datum.weigh-z)) for a in H.setOfArcs if ( (a.fromNode not in node_cover) and (a.toNode not in node_cover)} ) arcs2 = frozenset( {Arc(a.fromNode,a.toNode,ArcWeightDatum(a.datum.weigh)) for a in H.setOfArcs if ( (a.fromNode not in node_cover) != (a.toNode not in node_cover)} ) arcs3 = frozenset( {Arc(a.fromNode,a.toNode,ArcWeightDatum(a.datum.weigh+z)) for a in H.setOfArcs if ( (a.fromNode in node_cover) and (a.toNode in node_cover)} ) arcs = arcs1.union(arcs2.union(arcs3)) H = DiGraph(nodes,arcs) return value, assignment

Cette implémentation est O(N^{5}) car elle génère une nouvelle matching maximale à chaque itération ; similaire aux deux précédentes implémentations d' Edmonds-Karp , cet algorithme peut être modifié afin qu'il garde une trace de l'appariement et l'adapte intelligemment à chaque itération. Lorsque cela est fait, la complexité devient O(N^{4}) . Une version plus avancée et plus récente de cet algorithme (nécessitant des structures de données plus avancées) peut s'exécuter en O(N^{3}) . Les détails de l'implémentation plus simple ci-dessus et de l'implémentation plus avancée peuvent être trouvés dans cet article qui a motivé cet article de blog.

Aucune des opérations sur les poids d' arc ne modifie l'affectation finale renvoyée par l'algorithme. Voici pourquoi : étant donné que nos graphes d'entrée sont toujours des graphes bipartis complets , une solution doit mapper chaque nœud d'une partition à un autre nœud de la seconde partition, via l' arc entre ces deux nœuds . Notez que les opérations effectuées sur les poids des arcs ne changent jamais l'ordre (ordonné par poids) des arcs incidents sur un nœud particulier.

Ainsi, lorsque l'algorithme se termine par une correspondance bipartite parfaite et complète , chaque nœud se voit attribuer un arc de poids nul , puisque l'ordre relatif des arcs de ce nœud n'a pas changé au cours de l'algorithme, et qu'un arc de poids nul est l'arc le moins cher possible et l' appariement bipartite parfait et complet garantit qu'un tel arc existe pour chaque nœud . Cela signifie que la solution générée est en effet la même que la solution du problème d'affectation linéaire d'origine sans aucune modification des poids des arcs .

Affectations déséquilibrées

Il semble que l'algorithme soit assez limité puisque tel que décrit, il ne fonctionne que sur des graphes bipartis complets (ceux où la moitié des nœuds sont dans une partie de la bipartition et l'autre moitié dans la seconde partie). Chez le travailleur, la motivation à la tâche signifie qu'il y a autant de travailleurs que de tâches (semble assez limitatif).

Cependant, il existe une transformation facile qui supprime cette restriction. Supposons qu'il y ait moins de travailleurs que de tâches, nous ajoutons quelques travailleurs fictifs (assez pour faire du graphe résultant un graphe biparti complet ). Each dummy worker has an arc directed from the worker to each of the tasks. Each such arc has weight 0 (placing it in a matching gives no added profit). After this change the graph is a complete bipartite graph which we can solve for. Any task assigned a dummy worker is not initiated.

Suppose that there are more tasks than workers. We add some dummy tasks (enough to make the resulting graph a complete bipartite graph ). Each dummy task has an arc directed from each worker to the dummy task. Each such arc has a weight of 0 (placing it in a matching gives no added profit). After this change the graph is a complete bipartite graph which we can solve for. Any worker assigned to dummy task is not employed during the period.

A Linear Assignment Example

Finally, let's do an example with the code I've been using. I'm going to modify the example problem from here. We have 3 tasks: we need to clean the bathroom , sweep the floor , and wash the windows .

The workers available to use are Alice , Bob , Charlie , and Diane . Each of the workers gives us the wage they require per task. Here are the wages per worker:

If we want to pay the least amount of money, but still get all the tasks done, who should do what task? Start by introducing a dummy task to make the digraph representing the problem bipartite.

Supposing that the problem is encoded in a digraph , then kuhn_munkres( bipartition(G) ) will solve the problem and return the assignment. It's easy to verify that the optimal (lowest cost) assignment will cost $5.

Here's a visualization of the solution the code above generates:

That is it. You now know everything you need to know about the linear assignment problem.

You can find all of the code from this article on GitHub.

Lectures complémentaires sur le blog Toptal Engineering :

• Graphique Data Science avec Python/NetworkX