Maksymalny przepływ i problem przypisania liniowego

Oto problem: Twoja firma wyznacza kontrahentów do realizacji umów. Przeglądasz swoje dyżury i decydujesz, którzy kontrahenci są dostępni na jednomiesięczne zlecenie, a także przeglądasz dostępne kontrakty, aby zobaczyć, którzy z nich są na miesięczną pracę. Biorąc pod uwagę, że wiesz, jak skutecznie każdy wykonawca może wywiązać się z każdego kontraktu, jak przydzielić wykonawców, aby zmaksymalizować ogólną efektywność w danym miesiącu?

Jest to przykład problemu przypisania, a problem można rozwiązać za pomocą klasycznego algorytmu węgierskiego.

Algorytm węgierski (znany również jako algorytm Kuhna-Munkresa) to wielomianowy algorytm czasu, który maksymalizuje dopasowanie wag w ważonym grafie dwudzielnym. W tym przypadku wykonawców i kontrakty można zamodelować jako graf dwudzielny, a ich efektywność stanowi wagi krawędzi między wykonawcą a węzłami kontraktu.

W tym artykule poznasz implementację węgierskiego algorytmu, który wykorzystuje algorytm Edmondsa-Karpa do rozwiązania problemu przypisania liniowego. Dowiesz się również, w jaki sposób algorytm Edmondsa-Karpa jest niewielką modyfikacją metody Forda-Fulkersona i jak ważna jest ta modyfikacja.

Problem maksymalnego przepływu

Sam problem maksymalnego przepływu można opisać nieformalnie jako problem przenoszenia pewnej ilości płynu lub gazu przez sieć rur z jednego źródła do jednego terminala. Odbywa się to przy założeniu, że ciśnienie w sieci jest wystarczające, aby ciecz lub gaz nie zalegały na żadnej długości rury lub łącznika rurowego (w miejscach, w których spotykają się różne długości rur).

Dokonując pewnych zmian na wykresie, problem przypisania może zostać przekształcony w problem maksymalnego przepływu.

Czynności wstępne

Pomysły potrzebne do rozwiązania tych problemów pojawiają się w wielu dyscyplinach matematycznych i inżynierskich, często podobne koncepcje są znane pod różnymi nazwami i wyrażane na różne sposoby (np. macierze sąsiedztwa i listy sąsiedztwa). Ponieważ te idee są dość ezoteryczne, dokonuje się wyborów dotyczących tego, jak ogólnie te pojęcia zostaną zdefiniowane dla danego otoczenia.

Ten artykuł nie zakłada żadnej wcześniejszej wiedzy poza niewielką wstępną teorią mnogości.

Implementacje w tym poście przedstawiają problemy jako grafy skierowane (digraf).

Wykresy

Digraph ma dwa atrybuty, setOfNodes i setOfArcs . Oba te atrybuty są zestawami (kolekcje nieuporządkowane). W blokach kodu w tym poście faktycznie używam zamrożonego zestawu Pythona, ale ten szczegół nie jest szczególnie ważny.

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

(Uwaga: ten i cały inny kod w tym artykule są napisane w Pythonie 3.6.)

Węzły

Węzeł n składa się z dwóch atrybutów:

• n.uid : Unikalny identyfikator.

  Oznacza to, że dla dowolnych dwóch węzłów x i y ,

 x.uid != y.uid

• n.datum : Reprezentuje dowolny obiekt danych.

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

Łuki

Łuk a składa się z trzech atrybutów:

• a.fromNode : jest to węzeł zdefiniowany powyżej.

• a.toNode : To jest węzeł , jak zdefiniowano powyżej.

• a.datum : reprezentuje dowolny obiekt danych.

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

Zbiór łuków w dwuwykresie reprezentuje relację binarną na węzłach dwuwykresu . Istnienie arc a implikuje, że istnieje związek między a.fromNode i a.toNode .

W grafie skierowanym (w przeciwieństwie do grafu nieskierowanego) istnienie relacji między a.fromNode i a.toNode nie oznacza, że ​​istnieje podobna relacja między a.toNode i a.fromNode .

Dzieje się tak, ponieważ w grafie nieskierowanym wyrażona zależność niekoniecznie jest symetryczna.

Wykresy

Węzły i łuki mogą być użyte do zdefiniowania digrafu , ale dla wygody w poniższych algorytmach digraf będzie reprezentowany jako słownik.

Oto metoda, która może przekonwertować powyższą reprezentację wykresu na reprezentację słownikową podobną do tej:

 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

A oto kolejna, która przekształca to w słownik słowników, kolejna operacja, która będzie przydatna:

 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

Kiedy artykuł mówi o digrafie reprezentowanym przez słownik, będzie się do niego odwoływać G_as_dict .

Czasami pomocne jest pobranie węzła z wykresu G przez jego uid (unikalny identyfikator):

 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()

Definiując grafy, ludzie czasami używają terminów węzeł i wierzchołek w odniesieniu do tego samego pojęcia; to samo dotyczy terminów arc i edge.

Dwie popularne reprezentacje grafów w Pythonie to ta, która używa słowników i druga, która używa obiektów do reprezentowania grafów. Reprezentacja w tym poście znajduje się gdzieś pomiędzy tymi dwoma powszechnie używanymi reprezentacjami.

To jest moja reprezentacja dwugrafowa . Takich jest wiele, ale ten jest mój.

Spacery i ścieżki

Niech S_Arcs będzie skończoną sekwencją (uporządkowanym zbiorem) łuków na wykresie G tak, że jeśli a jest dowolnym łukiem w S_Arcs z wyjątkiem ostatniego, a b następuje po a w sekwencji, to musi istnieć węzeł n = a.fromNode w G.setOfNodes tak, że a.toNode = b.fromNode .

Zaczynając od pierwszego łuku w S_Arcs i kontynuując aż do ostatniego łuku w S_Arcs , zbierz (w napotkanej kolejności) wszystkie węzły n zdefiniowane powyżej między każdymi dwoma kolejnymi łukami w S_Arcs . Oznacz etykietą uporządkowaną kolekcję węzłów zebraną podczas tej operacji 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

• Jeśli którykolwiek węzeł pojawia się więcej niż jeden raz w sekwencji S_Nodes , wywołaj S_Arcs spacer po wykresie G .

• W przeciwnym razie wywołaj S_Arcs ścieżkę od list(S_Nodes)[0] do list(S_Nodes)[-1] na wykresie G .

Węzeł źródłowy

Call node s węzeł źródłowy w digrafie G , jeśli s jest w G.setOfNodes i G.setOfArcs nie zawiera łuku a takiego, że 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

Węzeł końcowy

Wywołaj węzeł t jako węzeł końcowy w digrafie G , jeśli t znajduje się w G.setOfNodes i G.setOfArcs nie zawiera łuku a takiego, że 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

Cięcia i st Cięcia

cut połączonego digrafu G jest podzbiorem łuków z G.setOfArcs , który dzieli zbiór węzłów G.setOfNodes w G . G jest połączony, jeśli każdy węzeł n w G.setOfNodes i ma co najmniej jeden łuk a w G.setOfArcs taki, że n = a.fromNode lub n = a.toNode , ale a.fromNode != a.toNode .

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

Powyższa definicja odnosi się do podzbioru łuków , ale może również określać podział węzłów G.setOfNodes .

Dla funkcji predecessors_of i successors_of n jest węzłem w zbiorze G.setOfNodes digrafu G G a cut jest kawałkiem 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

Niech cut będzie cięciem dwuznaku 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 jest cięciem digrafu G if: (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 nazywamy cięciem xy if (x in get_first_part(cut, G)) and (y in get_second_part(cut, G) ) and (x != y) . Kiedy węzeł x w cut xy jest węzłem źródłowym, a węzeł y w wycięciu xy jest węzłem końcowym , to cięcie to nazywa się cięciem st .

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

Sieci przepływu

Możesz użyć dwuznaku G do przedstawienia sieci przepływu.

Przypisz każdy węzeł n , gdzie n jest w G.setOfNodes n.datum , który jest FlowNodeDatum :

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

Przypisz każdemu łukowi a , gdzie a znajduje się w G.setOfArcs , a a.datum to FlowArcDatum .

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

flowNodeDatum.flowIn i flowNodeDatum.flowOut są dodatnimi liczbami rzeczywistymi.

flowArcDatum.capacity i flowArcDatum.flow są również dodatnimi liczbami rzeczywistymi.

Dla każdego węzła n w 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})

Digraph G reprezentuje teraz sieć przepływu .

Przepływ G odnosi się do a.flow dla wszystkich łuków a w G .

Możliwe przepływy

Niech digraf G reprezentuje sieć przepływów .

Sieć przepływowa reprezentowana przez G ma możliwe przepływy , jeżeli:

1. Dla każdego węzła n w G.setOfNodes z wyjątkiem węzłów źródłowych i końcowych : n.datum.flowIn = n.datum.flowOut .

2. Dla każdego łuku a w G.setOfNodes : a.datum.capacity <= a.datum.flow .

Warunek 1 nazywany jest ograniczeniem konserwacyjnym.

Warunek 2 nazywany jest ograniczeniem pojemności.

Wydajność cięcia

Zdolność cięcia st cut stCut z węzłem źródłowym s i węzłem końcowym t sieci przepływowej reprezentowanej przez wykres G wynosi:

 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

Minimalna wydajność cięcia

Niech stCut = stCut(s,t,cut) będzie st cięciem sieci przepływu reprezentowanej przez wykres G .

stCut to minimalna redukcja przepustowości sieci przepływu reprezentowana przez G , jeśli nie ma innego stcut stCutCompetitor w tej sieci przepływu, tak że:

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

Usuwanie przepływów

Chciałbym odnieść się do digrafu , który byłby wynikiem pobrania digrafu G i usunięcia wszystkich danych przepływu ze wszystkich węzłów w G.setOfNodes , a także wszystkich łuków w 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)

Problem z maksymalnym przepływem

Sieć przepływu reprezentowana jako digraf G , węzeł źródłowy s w G.setOfNodes i węzeł końcowy t w G.setOfNodes , G może reprezentować problem maksymalnego przepływu, jeżeli:

 (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)

Oznacz tę reprezentację:

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

Gdzie sourceNodeUid = s.uid , terminalNodeUid = t.uid , a maxFlowProblemStateUid to identyfikator wystąpienia problemu.

Maksymalne rozwiązanie przepływu

Niech maxFlowProblem reprezentuje problem maksymalnego przepływu . Rozwiązanie maxFlowProblem może być reprezentowane przez sieć przepływu reprezentowaną jako dwuwykres H .

Digraph H jest możliwym rozwiązaniem problemu maksymalnego przepływu na wejściu python maxFlowProblem , jeśli:

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

2. H jest siecią przepływową i ma możliwe przepływy .

Jeśli oprócz 1 i 2:

1. Nie może istnieć inna sieć przepływu reprezentowana przez digraf K taka, że strip_flows(G) == strip_flows(K) i find_node_by_uid(t.uid,G).flowIn < find_node_by_uid(t.uid,K).flowIn .

Wtedy H jest również optymalnym rozwiązaniem dla maxFlowProblem .

Innymi słowy, możliwe do zrealizowania rozwiązanie przepływu maksymalnego można przedstawić za pomocą wykresu , który:

1. Jest identyczny z wykresem G odpowiedniego problemu maksymalnego przepływu z wyjątkiem tego, że n.datum.flowIn , n.datum.flowOut i a.datum.flow dowolnego z węzłów i łuków mogą być różne.

2. Reprezentuje sieć przepływową z możliwymi przepływami .

I może reprezentować optymalne rozwiązanie maksymalnego przepływu, jeśli dodatkowo:

1. flowIn dla węzła odpowiadającego węzłowi końcowemu w problemie maksymalnego przepływu jest tak duży, jak to możliwe (gdy warunki 1 i 2 są nadal spełnione).

Jeśli wykres H reprezentuje możliwe rozwiązanie maksymalnego przepływu : find_node_by_uid(s.uid,H).flowOut = find_node_by_uid(t.uid,H).flowIn wynika z twierdzenia o maksymalnym przepływie, min cut (omówione poniżej). Nieformalnie, ponieważ zakłada się, że H ma wykonalne przepływy , oznacza to, że przepływ t może być ani „stworzony” (z wyjątkiem węzłów źródłowych ) ani „zniszczony” (z wyjątkiem węzła końcowego s ) podczas przekraczania dowolnego (innego) węzła ( ograniczenia konserwacyjne ).

Ponieważ problem maksymalnego przepływu zawiera tylko jeden węzeł źródłowy s i jeden węzeł końcowy t , cały przepływ „utworzony” w s musi zostać „zniszczony” w t lub sieć przepływów nie ma wykonalnych przepływów ( ograniczenie zachowania zostałoby naruszone ).

Niech dwuwykres H reprezentuje możliwe do zrealizowania maksymalne rozwiązanie przepływu ; powyższa wartość nazywana jest st. Flow Value of H .

Pozwolić:

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

Oznacza to, że mfps jest następcą stanu maxFlowProblem , co oznacza po prostu, że mfps jest dokładnie taki jak maxFlowProblem z wyjątkiem tego, że wartości a.flow dla łuków a w maxFlowProblem.setOfArcs mogą być inne niż a.flow dla łuków a w 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

Oto wizualizacja mfps wraz z skojarzonym z nim maxFlowProb . Każdy łuk a na obrazie ma etykietę, te etykiety to a.datum.flowFrom / a.datum.flowTo , każdy węzeł n na obrazie ma etykietę, a te etykiety to n.uid {n.datum.flowIn / a.datum.flowOut} .

St Cut Flow

Niech mfps będzie reprezentować MaxFlowProblemState , a stCut reprezentować wycięcie mfps.G . Przepływ cięcia stCut jest zdefiniowany:

 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

st cut flow to suma przepływów z partycji zawierającej węzeł źródłowy do partycji zawierającej węzeł końcowy pomniejszona o sumę przepływów z partycji zawierającej węzeł końcowy do partycji zawierającej węzeł źródłowy .

Maksymalny przepływ, minimalne cięcie

Niech maxFlowProblem reprezentuje problem maksymalnego przepływu i niech rozwiązanie maxFlowProblem będzie reprezentowane przez sieć przepływów reprezentowaną jako digraf H .

Niech minStCut będzie minimalną przepustowością sieci przepływowej reprezentowanej przez maxFlowProblem.G .

Ponieważ w problemie maksymalnego przepływu przepływ zaczyna się tylko w jednym węźle źródłowym i kończy się w jednym węźle końcowym , a ze względu na ograniczenia przepustowości i ograniczenia dotyczące zachowania wiemy, że cały przepływ wchodzący w skład maxFlowProblem.terminalNodeUid musi przecinać dowolne , w szczególności musi przecinać minStCut . To znaczy:

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

Rozwiązywanie problemu maksymalnego przepływu

Podstawowym pomysłem na rozwiązanie problemu maksymalnego przepływu maxFlowProblem jest rozpoczęcie od rozwiązania maksymalnego przepływu reprezentowanego przez wykres H . Takim punktem wyjścia może być H = strip_flows(maxFlowProblem.G) . Zadanie polega więc na użyciu H i pewnej zachłannej modyfikacji wartości a.datum.flow niektórych łuków a w H.setOfArcs w celu wytworzenia innego rozwiązania maksymalnego przepływu reprezentowanego przez wykres K , tak że K nie może nadal reprezentować sieci przepływów z możliwymi przepływami i get_flow_value(H, maxFlowProblem) < get_flow_value(K, maxFlowProblem) . Dopóki ten proces trwa, jakość ( get_flow_value(K, maxFlowProblem) ) ostatnio napotkanego rozwiązania maksymalnego przepływu ( K ) jest lepsza niż jakiegokolwiek innego znalezionego rozwiązania maksymalnego przepływu . Jeśli proces osiągnie punkt, w którym wie, że żadna inna poprawa nie jest możliwa, proces może się zakończyć i zwróci optymalny maksymalny przepływ roztworu .

Powyższy opis jest ogólny i pomija wiele dowodów, takich jak to, czy taki proces jest możliwy lub jak długo może potrwać, podam jeszcze kilka szczegółów i algorytm.

Twierdzenie o maksymalnym przepływie i minimalnym odcięciu

Z książki Flows in Networks autorstwa Forda i Fulkersona, twierdzenie o maksymalnym przepływie, twierdzeniu o minimalnym odcięciu (Twierdzenie 5.1) to:

Dla dowolnej sieci maksymalna wartość przepływu od s do t jest równa minimalnej przepustowości wszystkich cięć oddzielających s i t .

Korzystając z definicji w tym poście, przekłada się to na:

Rozwiązanie maxFlowProblem reprezentowanego przez sieć przepływową reprezentowaną jako wykres H jest optymalne, jeśli:

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

Podoba mi się ten dowód twierdzenia, a Wikipedia ma inny.

Do udowodnienia poprawności i kompletności metody Forda-Fulkersona służy twierdzenie o maksymalnym przepływie i minimalnym wycięciu .

Podam również dowód twierdzenia w sekcji po augmentacji ścieżek .

Metoda Forda-Fulkersona i algorytm Edmondsa-Karpa

CLRS definiuje metodę Forda-Fulkersona w ten sposób (sekcja 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

Wykres resztowy

Wykres rezydualny sieci przepływowej reprezentowany jako dwuznak G może być przedstawiony jako dwuwykres 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) zwraca sumę a.datum.capacity dla wszystkich łuków w podzbiorze G.setOfArcs , gdzie łuk a jest w podzbiorze, jeśli a.fromNode = n i a.toNode = u .

• agg_n_to_u_cap(n,u,G_as_dict) zwraca sumę a.datum.flow dla wszystkich łuków w podzbiorze G.setOfArcs , gdzie łuk a jest w podzbiorze, jeśli a.fromNode = n i a.toNode = u .

Krótko mówiąc, wykres reszt G_f reprezentuje pewne działania, które można wykonać na wykresie G .

Każda para węzłów n,u w G.setOfNodes sieci przepływów reprezentowanych przez wykres G może generować 0, 1 lub 2 łuki w grafie resztkowym G_f G .

1. Para n,u nie generuje żadnych łuków w G_f , jeśli nie ma łuku a w G.setOfArcs takiego, że a.fromNode = n i a.toNode = u .

2. Para n,u generuje łuk a w G_f.setOfArcs , gdzie a reprezentuje łuk oznaczony jako łuk przepływu wypychającego od n do 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. Para n,u generuje łuk a w G_f.setOfArcs , gdzie a reprezentuje łuk oznaczony jako łuk przepływu ciągnącego od n do u a = Arc(n,u,datum=ResidualNode(n_to_u_cap_sum - n_to_u_flow_sum)) jeśli n_to_u_flow_sum > 0.0 .

• Każdy łuk przepływu wypychanego w G_f.setOfArcs reprezentuje akcję polegającą na dodaniu sumy x <= n_to_u_cap_sum - n_to_u_flow_sum przepływu do łuków w podzbiorze G.setOfArcs , gdzie łuk a jest w podzbiorze, jeśli a.fromNode = n i a.toNode = u .

• Każdy łuk przepływu ciągnącego w G_f.setOfArcs reprezentuje akcję odjęcia sumy przepływu x <= n_to_u_flow_sum do łuków w podzbiorze G.setOfArcs , gdzie łuk a jest w podzbiorze, jeśli a.fromNode = n i a.toNode = u .

Wykonywanie pojedynczego działania push lub pull z G_f na odpowiednich łukach w G może wygenerować sieć przepływu bez wykonalnych przepływów , ponieważ ograniczenia przepustowości lub ograniczenia zachowania mogą zostać naruszone w wygenerowanej sieci przepływu .

Oto wizualizacja wykresu rezydualnego poprzedniej przykładowej wizualizacji rozwiązania maksymalnego przepływu , w wizualizacji każdy łuk a reprezentuje a.residualCapacity .

Ścieżka rozszerzająca

Niech maxFlowProblem będzie problemem z maksymalnym przepływem i niech G_f = get_residual_graph_of(G) będzie grafem rezydualnym maxFlowProblem.G .

Ścieżka augmentingPath dla maxFlowProblem to dowolna ścieżka od find_node_by_uid(maxFlowProblem.sourceNode,G_f) do find_node_by_uid(maxFlowProblem.terminalNode,G_f) .

Okazuje się, że ścieżka augmentingPath może być zastosowana do rozwiązania maksymalnego przepływu reprezentowanego przez wykres H , generując inne rozwiązanie maksymalnego przepływu reprezentowane przez wykres K gdzie get_flow_value get_flow_value(H, maxFlowProblem) < get_flow_value(K, maxFlowProblem) jeśli H nie jest optymalne .

Oto jak:

 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

W powyższym, TOL to pewna wartość tolerancji dla zaokrąglania wartości przepływu w sieci. Ma to na celu uniknięcie kaskadowej niedokładności obliczeń zmiennoprzecinkowych. Na przykład użyłem TOL = 10 , aby oznaczać zaokrąglenie do 10 cyfr znaczących.

Niech K = augment(augmentingPath, H) , wtedy K reprezentuje możliwe rozwiązanie maksymalnego przepływu dla maxFlowProblem . Aby stwierdzenie było prawdziwe, sieć przepływów reprezentowana przez K musi mieć możliwe przepływy (nie naruszać ograniczenia przepustowości ani ograniczenia zachowania ).

Oto dlaczego: W powyższej metodzie każdy węzeł dodany do nowej sieci przepływu reprezentowanej przez wykres K jest albo dokładną kopią węzła z wykresu H , albo węzła n , który ma ten sam numer dodany do swojego n.datum.flowIn co jego n.datum.flowOut . Oznacza to, że warunek konserwatorski jest spełniony w K , o ile był spełniony w H . Warunek zachowania jest spełniony, ponieważ wyraźnie sprawdzamy, czy każdy nowy łuk a w sieci ma a.datum.flow <= a.datum.capacity ; tak długo, jak łuki ze zbioru H.setOfArcs , które zostały skopiowane bez modyfikacji do K.setOfArcs , nie naruszają ograniczenia pojemności , to K nie narusza ograniczenia pojemności .

Prawdą jest również, że get_flow_value(H, maxFlowProblem) < get_flow_value(K, maxFlowProblem) jeśli H nie jest optymalne .

Oto dlaczego: Aby ścieżka augmentingPath istniała w dwugrafowej reprezentacji grafu rezydualnego G_f z maksymalnym przepływem maxFlowProblem , to ostatni łuk a na augmentingPath musi być łukiem „push” i musi mieć a.toNode == maxFlowProblem.terminalNode . Ścieżka rozszerzająca jest zdefiniowana jako taka, która kończy się w węźle końcowym problemu maksymalnego przepływu, dla którego jest ścieżką rozszerzającą . Z definicji grafu resztowego jasno wynika, że ​​ostatni łuk na ścieżce rozszerzającej na tym grafie resztkowym musi być łukiem „pchania”, ponieważ każdy łuk „ciągnący” b na ścieżce rozszerzającej będzie miał b.toNode == maxFlowProblem.terminalNode i b.fromNode != maxFlowProblem.terminalNode z definicji ścieżki . Dodatkowo, z definicji path , jasno wynika, że ​​węzeł końcowy jest modyfikowany tylko raz metodą augment . W ten sposób augment modyfikuje maxFlowProblem.terminalNode.flowIn dokładnie raz i zwiększa wartość maxFlowProblem.terminalNode.flowIn , ponieważ ostatni łuk w augmentingPath musi być łukiem, który powoduje modyfikację w maxFlowProblem.terminalNode.flowIn podczas augment . Zgodnie z definicją augment , która ma zastosowanie do łuków „push”, wartość maxFlowProblem.terminalNode.flowIn może być tylko zwiększona, a nie zmniejszona.

Niektóre dowody z Sedgewick i Wayne

Książka Algorithms, czwarte wydanie autorstwa Roberta Sedgewicha i Kevina Wayne'a zawiera wspaniałe i krótkie dowody (strony 892-894), które będą przydatne. Odtworzę je tutaj, choć użyję języka pasującego do poprzednich definicji. Moje etykiety próbne są takie same jak w książce Sedgewicka.

Twierdzenie E: Dla dowolnego digrafu H reprezentującego możliwe rozwiązanie problemu maksymalnego przepływu maxFlowProblem , dla dowolnego stCut get_stcut_flow(stCut,H,maxFlowProblem) = get_flow_value(H, maxFlowProblem) .

Dowód: Niech stCut=stCut(maxFlowProblem.sourceNode,maxFlowProblem.terminalNode,set([a for a in H.setOfArcs if a.toNode == maxFlowProblem.terminalNode])) . Twierdzenie E obowiązuje dla stCut bezpośrednio z definicji wartości przepływu st . Załóżmy, że tam chcemy przenieść jakiś węzeł n z partycji s ( get_first_part(stCut.cut, G) ) i do partycji t (get_second_part(stCut.cut, G)) , w tym celu musimy zmienić stCut.cut , który może zmienić stcut_flow = get_stcut_flow(stCut,H,maxFlowProblem) i unieważnić propozycję E . Zobaczmy jednak, jak zmieni się wartość stcut_flow po wprowadzeniu tej zmiany. węzeł n jest w równowadze, co oznacza, że ​​suma przepływu do węzła n jest równa sumie przepływu z niego (jest to konieczne, aby H stanowiło możliwe rozwiązanie ). Zauważże cały przepływ będący częścią węzła wejściowego stcut_flow n wchodzi do niego z partycji s (przepływ wchodzący do węzła n z partycji t bezpośrednio lub pośrednio nie zostałby uwzględniony w wartości stcut_flow ponieważ zmierza w złym kierunku na podstawie definicji). Dodatkowo cały przepływ wychodzący z n ostatecznie (bezpośrednio lub pośrednio) wpłynie do węzła końcowego (udowodniono wcześniej). Kiedy przenosimy węzeł n do partycji t, cały przepływ wchodzący z partycji n musi zostać dodany do nowej wartości stcut_flow ; jednak cały przepływ wychodzący z n musi zostać odjęty od nowej wartości stcut_flow ; część przepływu kierująca się bezpośrednio do partycji t jest odejmowana, ponieważ ten przepływ jest teraz wewnętrzny dla nowej partycji t i nie jest liczony jako stcut_flow . Część przepływu z n nagłówków do węzłów w partycji s musi również zostać odjęta od stcut_flow : Po przeniesieniu n do partycji t, przepływy te zostaną skierowane z partycji t do partycji s i tak nie należy uwzględniać w stcut_flow , ponieważ te przepływy są usuwane, dopływ do partycji s i muszą być pomniejszone o sumę tych przepływów oraz odpływ z partycji s do partycji t (gdzie wszystkie przepływy z st musi się skończyć) należy zmniejszyć o tę samą kwotę. Ponieważ węzeł n był w równowadze przed procesem, aktualizacja doda tę samą wartość do nowej wartości stcut_flow , co odjęta, pozostawiając w ten sposób twierdzenie E prawdziwe po aktualizacji. Trafność zdania E wynika zatem z indukcji wielkości t-partycji.

Oto kilka przykładowych sieci przepływów, które pomogą zwizualizować mniej oczywiste przypadki, w których zachodzi propozycja E ; na obrazie czerwone obszary oznaczają partycję s, niebieskie obszary reprezentują partycję t, a zielone łuki oznaczają cięcie st . Na drugim obrazie przepływ między węzłem A a węzłem B wzrasta, podczas gdy przepływ do węzła końcowego t nie ulega zmianie.:

Wniosek: Żadna wartość przepływu st cut nie może przekroczyć wydajności żadnego st cut .

Twierdzenie F. (przepływ maksymalny, twierdzenie o minimalnym przekroju): Niech f będzie przepływem st . Następujące 3 warunki są równoważne:

1. Istnieje odcięcie , którego przepustowość jest równa wartości przepływu f .

2. f jest przepływem maksymalnym .

3. Nie ma ścieżki rozszerzającej w odniesieniu do f .

Warunek 1 implikuje warunek 2 w następstwie. Warunek 2 implikuje warunek 3, ponieważ istnienie ścieżki zwiększającej implikuje istnienie przepływu o większych wartościach, sprzecznych z maksymalizacją f . Warunek 3 implikuje warunek 1: Niech C_s będzie zbiorem wszystkich węzłów , do których można dotrzeć z s ze ścieżką rozszerzającą w grafie rezydualnym . Niech C_t będzie pozostałymi łukami , wtedy t musi być w C_t (z naszego założenia). Łuki przechodzące z C_s do C_t tworzą następnie cięcie st, które zawiera tylko łuki a gdzie a.datum.capacity = a.datum.flow lub a.datum.flow = 0.0 . Gdyby było inaczej, węzły połączone łukiem z pozostałą pojemnością resztkową do C_s byłyby w zbiorze C_s , ponieważ istniałaby wtedy ścieżka zwiększająca się od s do takiego węzła . Przepływ przez przecięcie jest równe przepustowości przecięcia (ponieważ łuki od C_s do C_t mają przepływ równy przepustowości), a także wartości przepływu (według propozycji E ).

To twierdzenie o maksymalnym przepływie i minimalnym odcięciu implikuje wcześniejsze stwierdzenie z Flows in Networks.

Wniosek (właściwość integralności): Gdy pojemności są liczbami całkowitymi, istnieje maksymalny przepływ o wartości całkowitej i algorytm Forda-Fulkersona znajduje go.

Dowód: Każda ścieżka zwiększająca zwiększa przepływ o dodatnią liczbę całkowitą, minimum niewykorzystanych pojemności w łukach „pchających” i przepływy w łukach „ciągnących”, z których wszystkie są zawsze dodatnimi liczbami całkowitymi.

Uzasadnia to opis metody Forda-Fulkersona z CLRS . Metoda polega na ciągłym znajdowaniu ścieżek rozszerzających i stosowaniu augment do najnowszego maxFlowSolution i wymyślaniu lepszych rozwiązań, aż do momentu, gdy nie będzie już ścieżek rozszerzających, co oznacza, że ​​najnowsze rozwiązanie maksymalnego przepływu jest optymalne .

Od Forda-Fulkersona do Edmonds-Karpa

Pozostałe pytania dotyczące rozwiązywania problemów maksymalnego przepływu to:

1. Jak należy budować ścieżki rozszerzające ?

2. Czy metoda zakończy się, jeśli użyjemy liczb rzeczywistych, a nie liczb całkowitych?

3. Ile czasu zajmie rozwiązanie (jeśli tak)?

Algorytm Edmondsa-Karpa określa, że ​​każda ścieżka rozszerzająca jest konstruowana przez pierwsze przeszukiwanie wszerz (BFS) grafu resztkowego ; okazuje się, że decyzja z punktu 1 powyżej również wymusi zakończenie działania algorytmu (punkt 2) i pozwoli na wyznaczenie asymptotycznej złożoności czasowej i przestrzennej.

Po pierwsze, oto implementacja 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

Użyłem deque z modułu kolekcji Pythona.

Aby odpowiedzieć na pytanie 2 powyżej, sparafrazuję inny dowód z Sedgewicka i Wayne'a: ​​Propozycja G. Liczba ścieżek rozszerzających potrzebnych w algorytmie Edmondsa-Karpa z N węzłami i łukami A wynosi co najwyżej NA/2 . Dowód: Każda ścieżka rozszerzania ma łuk wąskiego gardła — łuk , który jest usuwany z wykresu resztkowego, ponieważ odpowiada albo łukowi „pchającemu”, który zostaje wypełniony do pełnej pojemności, albo łukowi „ciągnącemu”, przez który przepływ staje się równy 0. Za każdym razem łuk staje się łukiem wąskim gardłem , długość każdej rozszerzającej się ścieżki przez niego musi wzrosnąć o współczynnik 2. Dzieje się tak dlatego, że każdy węzeł na ścieżce może pojawić się tylko raz lub wcale (z definicji ścieżki ), ponieważ ścieżki są eksplorowane od najkrótszej ścieżki do najdłuższej, co oznacza, że ​​co najmniej jeden węzeł musi zostać odwiedzony przez następną ścieżkę przechodzącą przez dany węzeł wąskiego gardła, co oznacza dodatkowe 2 łuki na ścieżce, zanim dotrzemy do węzła . Ponieważ ścieżka rozszerzająca ma długość co najwyżej N , każdy łuk może znajdować się na co najwyżej N/2 ścieżek rozszerzających , a całkowita liczba ścieżek rozszerzających wynosi co najwyżej NA/2 .

Algorytm Edmondsa-Karpa jest wykonywany w O(NA^2) . Jeśli podczas algorytmu zbadane zostaną co najwyżej ścieżki NA/2 , a eksploracja każdej ścieżki za pomocą BFS to N+A , wówczas najważniejszym terminem iloczynu i stąd złożoność asymptotyczna jest O(NA^2) .

Niech mfp będzie 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

Powyższa wersja jest nieefektywna i ma gorszą złożoność niż O(NA^2) , ponieważ za każdym razem konstruuje nowe rozwiązanie maksymalnego przepływu i nowy wykres resztkowy (zamiast modyfikować istniejące digrafy w miarę postępu algorytmu). Aby uzyskać prawdziwe rozwiązanie O(NA^2) , algorytm musi zachować zarówno wykres reprezentujący stan problemu maksymalnego przepływu, jak i powiązany z nim wykres rezydualny . Dlatego algorytm musi unikać niepotrzebnego iterowania po łukach i węzłach oraz aktualizować ich wartości i powiązane wartości na wykresie rezydualnym tylko w razie potrzeby.

Aby napisać szybszy algorytm Edmonds Karp , przepisałem kilka fragmentów kodu z powyższego. Mam nadzieję, że przejrzenie kodu, który wygenerował nowy digraf , pomogło w zrozumieniu o co chodzi. W szybkim algorytmie używam kilku nowych sztuczek i struktur danych Pythona, których nie chcę szczegółowo omawiać. Powiem, że a.fromNode i a.toNode są teraz traktowane jako stringi i uids do węzłów . Dla tego kodu niech mfps będzie 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

Oto wizualizacja tego, jak ten algorytm rozwiązuje przykładową sieć przepływu z góry. Wizualizacja pokazuje kroki, tak jak są one odzwierciedlone na wykresie G reprezentującym najbardziej aktualną sieć przepływów i jak są odzwierciedlone na wykresie rezydualnym tej sieci przepływów. Ścieżki rozszerzające na wykresie resztowym są pokazane jako ścieżki na czerwono, a wykres reprezentujący problem zbioru węzłów i łuków , na które ma wpływ dana ścieżka rozszerzająca, jest podświetlony na zielono. W każdym przypadku zaznaczę części wykresu, które zostaną zmienione (na czerwono lub zielono), a następnie pokażę wykres po zmianach (tylko na czarno).

Oto kolejna wizualizacja tego, jak ten algorytm rozwiązuje inną przykładową sieć przepływu . Zauważ, że ten używa liczb rzeczywistych i zawiera wiele łuków z tymi samymi wartościami fromNode i toNode .

**Zauważ również, że ponieważ łuki z resztowym punktem odniesienia mogą być częścią Ścieżki rozszerzającej, węzły dotknięte w DiGraph of the Flowed Network _mogą nie znajdować się na ścieżce w G! .

Wykresy dwustronne

Załóżmy, że mamy dwuwykres G , G jest dwudzielny, jeśli można podzielić węzły w G.setOfNodes na dwa zestawy ( part_1 i part_2 ) tak, że dla dowolnego łuku a G.setOfArcs nie może być prawdą , że a.fromNode w part_1 i a.toNode w part_1 . Nie może być również prawdą , że a.fromNode w part_2 i a.toNode w part_2 .

Innymi słowy, G jest dwudzielny , jeśli można go podzielić na dwa zestawy węzłów , tak że każdy łuk musi łączyć węzeł w jednym zestawie z węzłem w drugim zestawie.

Testowanie Dwustronne

Załóżmy, że mamy dwuznak G , chcemy sprawdzić, czy jest dwuczęściowy . Możemy to zrobić w O(|G.setOfNodes|+|G.setOfArcs|) przez zachłanne pokolorowanie wykresu na dwa kolory.

Najpierw musimy wygenerować nowy digraf H . Ten wykres będzie miał ten sam zestaw węzłów co G , ale będzie miał więcej łuków niż G . Każdy łuk a w G utworzy 2 łuki w H ; pierwszy łuk będzie identyczny z a , a drugi łuk odwraca kierunek 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)

Dopasowania i maksymalne dopasowania

Załóżmy, że mamy digraf G , a matching jest podzbiorem łuków z G.setOfArcs . matching jest dopasowaniem, jeśli dla dowolnych dwóch łuków a i b w matching : len(frozenset( {a.fromNode} ).union( {a.toNode} ).union( {b.fromNode} ).union( {b.toNode} )) == 4 . Innymi słowy, żadne dwa łuki w dopasowanych elementach nie mają wspólnego węzła .

Dopasowanie , jest maksymalnym matching , jeśli nie ma innego dopasowania alt_matching w G , takiego jak len(matching) < len(alt_matching) . Innymi słowy, matching jest maksymalnym dopasowaniem , jeśli jest to największy zestaw łuków z G.setOfArcs , który nadal spełnia definicję dopasowania (dodanie dowolnego łuku, który nie jest jeszcze w dopasowaniu, spowoduje przerwanie definicji dopasowania ).

Maksymalne dopasowanie matching idealnym dopasowaniem , jeśli dla każdego węzła n w G.setOfArcs istnieje łuk a w matching , gdzie a.fromNode == n or a.toNode == n .

Maksymalne dopasowanie dwustronne

Maksymalne dopasowanie dwudzielne to maksymalne dopasowanie na dwudzielnym wykresie G .

Biorąc pod uwagę, że G jest dwudzielne , problem znalezienia maksymalnego dopasowania dwudzielnego można przekształcić w problem maksymalnego przepływu, który można rozwiązać za pomocą algorytmu Edmondsa-Karpa , a następnie maksymalne dopasowanie dwudzielne można odzyskać z rozwiązania problemu maksymalnego przepływu .

Niech bipartition będzie dwudzielnością G .

Aby to zrobić, muszę wygenerować nowy digraph ( H ) z kilkoma nowymi węzłami ( H.setOfNodes ) i kilkoma nowymi łukami ( H.setOfArcs ). H.setOfNodes zawiera wszystkie węzły w G.setOfNodes i jeszcze dwa węzły , s ( węzeł źródłowy ) i t ( węzeł końcowy ).

H.setOfArcs będzie zawierać jeden łuk dla każdego G.setOfArcs . Jeśli łuk a znajduje się w G.setOfArcs , a.fromNode w bipartition.firstPart , a a.toNode w bipartition.secondPart , dodaj a w H (dodając FlowArcDatum(1,0) ).

Jeśli a.fromNode znajduje się w bipartition.secondPart , a a.toNode w bipartition.firstPart , należy uwzględnić Arc(a.toNode,a.fromNode,FlowArcDatum(1,0)) w H.setOfArcs .

Definicja grafu dwudzielnego zapewnia, że ​​żaden łuk nie łączy węzłów , w których oba węzły znajdują się w tej samej partycji. H.setOfArcs zawiera również łuk od węzła do każdego węzła s bipartition.firstPart . Wreszcie H.setOfArcs zawiera łuk każdego węzła w bipartition.secondPart do węzła t . a.datum.capacity = 1 dla wszystkich a w H.setOfArcs .

Najpierw podziel węzły w G.setOfNodes na dwa rozłączne zbiory ( part1 i part2 ) tak, aby żaden łuk w G.setOfArcs nie był skierowany z jednego zbioru do tego samego zbioru (ten podział jest możliwy, ponieważ G jest dwudzielny ). Następnie dodaj wszystkie łuki w G.setOfArcs , które są skierowane z part1 do part2 do H.setOfArcs . Następnie utwórz pojedynczy węzeł źródłowy s i pojedynczy węzeł końcowy t i utwórz więcej łuków

Następnie skonstruuj 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

Minimalna osłona węzła

Pokrycie węzła w digrafie G to zbiór węzłów ( cover ) z G.setOfNodes taki, że dla a łuku G.setOfArcs musi to być prawda: (a.fromNode in cover) or (a.toNode in cover) .

Minimalne pokrycie węzłów to najmniejszy możliwy zestaw węzłów na grafie, który nadal jest pokryciem węzłów . Twierdzenie Koniga mówi, że w grafie dwudzielnym rozmiar maksymalnego dopasowania na tym grafie jest równy rozmiarowi minimalnego pokrycia węzła i sugeruje, w jaki sposób pokrycie węzła może zostać odzyskane z maksymalnego dopasowania :

Załóżmy, że mamy bipartycję bipartition i maksymalne matching . Zdefiniuj nowy digraf H , H.setOfNodes=G.setOfNodes , łuki w H.setOfArcs są sumą dwóch zbiorów.

Pierwszy zestaw to łuki a w matching , z tą zmianą, że jeśli a.fromNode in bipartition.firstPart i a.toNode in bipartition.secondPart , a.fromNode i a.toNode zostaną zamienione w utworzonym łuku , dają takim łukom a.datum.inMatching=True atrybut wskazujący, że pochodzą one z łuków w dopasowaniu .

Drugi zestaw a łuki NIE w matching , z tą zmianą, że jeśli a.fromNode in bipartition.secondPart i a.toNode in bipartition.firstPart , to a.fromNode i a.toNode zostaną zamienione w utworzonym łuku (nadaj takim łukom a.datum.inMatching=False atrybut).

Następnie uruchom wyszukiwanie wgłębne (DFS), zaczynając od każdego węzła n w bipartition.firstPart , który nie jest ani n == a.fromNode ani n == a.toNodes dla dowolnego matching łuku a . Podczas DFS niektóre węzły są odwiedzane, a inne nie (przechowuj te informacje w polu n.datum.visited ). Minimalne pokrycie węzłów to suma węzłów {a.fromNode for a in H.setOfArcs if ( (a.datum.inMatching) and (a.fromNode.datum.visited) ) } oraz węzłów {a.fromNode for a in H.setOfArcs if (a.datum.inMatching) and (not a.toNode.datum.visited)} .

Można wykazać, że prowadzi to od maksymalnego dopasowania do minimalnego pokrycia węzła za pomocą dowodu przez sprzeczność, weź pewien łuk a , który rzekomo nie został pokryty i rozważ wszystkie cztery przypadki dotyczące tego, czy a.fromNode i a.toNode należą (czy jako toNode , czy fromNode ) matching dowolnego łuku w dopasowywaniu . Każdy przypadek prowadzi do sprzeczności ze względu na kolejność odwiedzania węzłów przez system DFS oraz fakt, że matching jest maksymalnym dopasowaniem .

Załóżmy, że mamy funkcję wykonującą te kroki i zwracającą zestaw węzłów obejmujący minimalne pokrycie węzła , gdy mamy digraf G i maksymalne 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

Problem przypisania liniowego

Problem przypisania liniowego polega na znalezieniu maksymalnego dopasowania wagi w ważonym dwudzielnym wykresie.

Problemy takie jak ten na samym początku tego postu można wyrazić jako liniowy problem przypisania . Mając zbiór pracowników, zbiór zadań oraz funkcję wskazującą opłacalność przydziału jednego pracownika do jednego zadania, chcemy zmaksymalizować sumę wszystkich przydziałów, które wykonujemy; jest to problem przypisania liniowego .

Załóżmy, że liczba zadań i pracowników jest równa, choć pokażę, że to założenie jest łatwe do usunięcia. W implementacji reprezentuję wagi łuków z atrybutem a.datum.weight dla łuku a .

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

Algorytm Kuhna-Munkresa

Algorytm Kuhna-Munkresa rozwiązuje problem przypisania liniowego . Dobra implementacja może zająć czas O(N^{4}) , (gdzie N jest liczbą węzłów na wykresie reprezentującym problem). Łatwiejsza do wyjaśnienia implementacja przyjmuje O(N^{5}) (dla wersji, która regeneruje DiGraphs ) i O(N^{4}) for (dla wersji, która obsługuje DiGraphs ). Jest to podobne do dwóch różnych implementacji algorytmu Edmondsa-Karpa .

W tym opisie pracuję tylko z pełnymi grafami dwudzielnymi (takimi, w których połowa węzłów znajduje się w jednej części dwudzielności , a druga połowa w drugiej części). W przypadku pracownika motywacja zadaniowa oznacza, że ​​jest tyle pracowników, ile zadań.

Wydaje się, że jest to istotny warunek (co jeśli te zestawy nie są równe!), ale łatwo jest rozwiązać ten problem; Mówię o tym, jak to zrobić w ostatniej sekcji.

Opisana tutaj wersja algorytmu wykorzystuje użyteczną koncepcję łuków o zerowej wadze . Niestety, ta koncepcja ma sens tylko wtedy, gdy rozwiązujemy minimalizację (jeśli zamiast maksymalizować zyski z naszych przydziałów zadań pracowniczych, zamiast tego minimalizowaliśmy koszty takich zadań).

Na szczęście łatwo jest przekształcić problem maksymalnego liniowego przypisania w problem minimalnego przypisania liniowego , ustawiając każdy łuk a weights na Ma.datum.weight , gdzie M=max({a.datum.weight for a in G.setOfArcs}) . Rozwiązanie pierwotnego problemu maksymalizacji będzie identyczne z rozwiązaniem problemu minimalizacji po zmianie wag łuków . Więc w pozostałej części załóżmy, że wprowadzamy tę zmianę.

Algorytm Kuhna-Munkresa rozwiązuje minimalne dopasowanie wagi w ważonym grafie dwudzielnym przez sekwencję maksymalnych dopasowań w nieważonych grafach dwudzielnych . Jeśli a znajdziemy idealne dopasowanie w reprezentacji dwuwykresu liniowego problemu przypisania i jeśli waga każdego łuku w dopasowaniu wynosi zero, to znaleźliśmy dopasowanie o minimalnej wadze, ponieważ to dopasowanie sugeruje, że wszystkie węzły w dwugrafie zostały dopasowany łukiem o najniższym możliwym koszcie (żaden koszt nie może być niższy niż 0, na podstawie wcześniejszych definicji).

Do dopasowywania nie można dodać żadnych innych łuków (ponieważ wszystkie węzły są już dopasowane) i żadne łuki nie powinny być usuwane z dopasowywania , ponieważ każdy możliwy łuk zastępczy będzie miał co najmniej taką samą wagę.

Jeśli znajdziemy maksymalne dopasowanie podwykresu G , który zawiera tylko łuki o zerowej wadze i nie jest to idealne dopasowanie , nie mamy pełnego rozwiązania (ponieważ dopasowanie nie jest idealne ). Możemy jednak utworzyć nowy digraf H , zmieniając wagi łuków w G.setOfArcs w taki sposób, że pojawiają się nowe łuki o zerowej wadze, a optymalne rozwiązanie H jest takie samo, jak optymalne rozwiązanie G . Ponieważ gwarantujemy, że w każdej iteracji powstaje co najmniej jeden łuk o zerowej wadze , gwarantujemy, że osiągniemy idealne dopasowanie w nie więcej niż |G.setOfNodes|^{2}=N^{2} takich iteracji.

Załóżmy, że w bipartition bipartition , bipartition.firstPart zawiera węzły reprezentujące procesy robocze, a bipartition.secondPart reprezentuje węzły reprezentujące zadania.

Algorytm rozpoczyna się od wygenerowania nowego wykresu H . H.setOfNodes = G.setOfNodes . Niektóre łuki w H są generowane z węzłów n w bipartition.firstPart . Każdy taki węzeł n generuje łuk b w H.setOfArcs dla każdego łuku a w bipartition.G.setOfArcs gdzie a.fromNode = n lub a.toNode = n , b=Arc(a.fromNode, a.toNode, a.datum.weight - z) gdzie z=min(x.datum.weight for x in G.setOfArcs if ( (x.fromNode == n) or (x.toNode == n) )) .

Więcej łuków w H jest generowanych z węzłów n w bipartition.secondPart . Każdy taki węzeł n generuje łuk b w H.setOfArcs dla każdego łuku a w bipartition.G.setOfArcs gdzie a.fromNode = n lub a.toNode = n , b=Arc(a.fromNode, a.toNode, ArcWeightDatum(a.datum.weight - z)) gdzie z=min(x.datum.weight for x in G.setOfArcs if ( (x.fromNode == n) or (x.toNode == n) )) .

KMA: Następnie utwórz nowy wykres K składający się tylko z łuków o zerowej wadze z H , oraz węzłów padających na te łuki . Utwórz podwójną partycję na węzłach w K , a następnie użyj solve_mbm( bipartition ) bipartition aby uzyskać maksymalne dopasowanie ( matching ) na K . Jeśli matching jest idealnym dopasowaniem w H ( łuki w matching występują na wszystkich węzłach w H.setOfNodes ) H.setOfNodes wówczas matching jest optymalnym rozwiązaniem problemu przypisania liniowego .

W przeciwnym razie, jeśli matching nie jest idealne , wygeneruj minimalne pokrycie węzła K przy użyciu node_cover = get_min_node_cover(matching, bipartition(K)) . Następnie zdefiniuj z=min({a.datum.weight for a in H.setOfArcs if a not in node_cover}) . Zdefiniuj 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)} != w poprzednim wyrażeniu działa jako operator XOR .Następnie arcs = arcs1.union(arcs2.union(arcs3)) , H=DiGraph(nodes,arcs) do etykieta KMA Algorytm jest kontynuowany aż do uzyskania idealnego dopasowania.To dopasowanie jest również rozwiązaniem problemu liniowego przypisania .

 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

Ta implementacja jest O(N^{5}) , ponieważ generuje nowe maksymalne matching w każdej iteracji; podobnie jak w poprzednich dwóch implementacjach Edmondsa-Karpa , ten algorytm można zmodyfikować tak, aby śledził dopasowanie i inteligentnie dostosowywał je do każdej iteracji. Kiedy to zostanie zrobione, złożoność staje się O(N^{4}) . Bardziej zaawansowana i nowsza wersja tego algorytmu (wymagająca bardziej zaawansowanych struktur danych) może działać w O(N^{3}) . Szczegóły zarówno prostszej implementacji powyżej, jak i bardziej zaawansowanej implementacji można znaleźć w tym poście, który zmotywował ten wpis na blogu.

Żadna z operacji na wagach łuku nie modyfikuje końcowego przypisania zwracanego przez algorytm. Oto dlaczego: Ponieważ nasze grafy wejściowe są zawsze kompletnymi grafami dwudzielnymi , rozwiązanie musi odwzorować każdy węzeł w jednej partycji na inny węzeł w drugiej partycji, za pomocą łuku między tymi dwoma węzłami . Zauważ, że operacje wykonywane na wagach łuków nigdy nie zmieniają kolejności (uporządkowanej według wagi) łuków padających na konkretny węzeł .

Tak więc, gdy algorytm kończy się w idealnym pełnym dwudzielnym dopasowaniu , każdemu węzłowi przypisywany jest łuk o zerowej wadze , ponieważ względna kolejność łuków z tego węzła nie zmieniła się podczas algorytmu, a ponieważ łuk o zerowej wadze jest najtańszym możliwym łukiem i idealne pełne dopasowanie dwudzielne gwarantuje, że dla każdego węzła istnieje jeden taki łuk . Oznacza to, że wygenerowane rozwiązanie jest rzeczywiście takie samo jak rozwiązanie z oryginalnego problemu przypisania liniowego bez jakiejkolwiek modyfikacji wag łuków .

Niezrównoważone zadania

Wygląda na to, że algorytm jest dość ograniczony, ponieważ zgodnie z opisem działa tylko na pełnych grafach dwudzielnych (tych, w których połowa węzłów znajduje się w jednej części dwudzielności , a druga połowa w drugiej części). W przypadku pracownika motywacja zadaniowa oznacza, że ​​jest tyle pracowników, ile zadań (wydaje się dość ograniczające).

Istnieje jednak łatwa transformacja, która usuwa to ograniczenie. Załóżmy, że jest mniej pracowników niż zadań, dodajemy kilka fikcyjnych pracowników (wystarczająco, aby wynikowy wykres był kompletnym grafem dwudzielnym ). 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.

Dalsza lektura na blogu Toptal Engineering:

• Wykresy analizy danych za pomocą Python/NetworkX