Portata massima e problema di assegnazione lineare

Ecco un problema: la tua azienda assegna appaltatori per adempiere ai contratti. Esamini i tuoi elenchi e decidi quali appaltatori sono disponibili per un impegno di un mese e guardi i tuoi contratti disponibili per vedere quali di loro sono per attività di un mese. Dato che sai con quanta efficacia ogni appaltatore può adempiere a ciascun contratto, come assegni gli appaltatori per massimizzare l'efficacia complessiva per quel mese?

Questo è un esempio del problema di assegnazione e il problema può essere risolto con il classico algoritmo ungherese.

L'algoritmo ungherese (noto anche come algoritmo di Kuhn-Munkres) è un algoritmo di tempo polinomiale che massimizza la corrispondenza del peso in un grafo bipartito ponderato. Qui, gli appaltatori e i contratti possono essere modellati come un grafo bipartito, con la loro efficacia come i pesi degli spigoli tra l'appaltatore e i nodi del contratto.

In questo articolo imparerai un'implementazione dell'algoritmo ungherese che utilizza l'algoritmo di Edmonds-Karp per risolvere il problema dell'assegnazione lineare. Imparerai anche come l'algoritmo Edmonds-Karp sia una leggera modifica del metodo Ford-Fulkerson e come questa modifica sia importante.

Il problema del flusso massimo

Lo stesso problema di flusso massimo può essere descritto informalmente come il problema di spostare un fluido o gas attraverso una rete di tubi da un'unica sorgente a un unico terminale. Ciò viene fatto partendo dal presupposto che la pressione nella rete sia sufficiente per garantire che il fluido o il gas non possano rimanere in nessuna lunghezza di tubo o raccordo (quei punti in cui si incontrano diverse lunghezze di tubo).

Apportando alcune modifiche al grafico, il problema di assegnazione può essere trasformato in un problema di flusso massimo.

Preliminari

Le idee necessarie per risolvere questi problemi sorgono in molte discipline matematiche e ingegneristiche, spesso concetti simili sono conosciuti con nomi diversi ed espressi in modi diversi (es. matrici di adiacenza e liste di adiacenza). Dal momento che queste idee sono piuttosto esoteriche, vengono fatte delle scelte riguardo a come generalmente questi concetti saranno definiti per un dato ambiente.

Questo articolo non presuppone alcuna conoscenza precedente oltre a una piccola teoria degli insiemi introduttiva.

Le implementazioni in questo post rappresentano i problemi come grafici diretti (digraph).

Grafici

Un digraph ha due attributi, setOfNodes e setOfArcs . Entrambi questi attributi sono insiemi (raccolte non ordinate). Nei blocchi di codice di questo post, sto effettivamente usando il frozenset di Python, ma quel dettaglio non è particolarmente importante.

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

(Nota: questo, e tutto il resto del codice in questo articolo, sono scritti in Python 3.6.)

Nodi

Un nodo n è composto da due attributi:

• n.uid : un identificatore univoco.

  Ciò significa che per due nodi y x

 x.uid != y.uid

• n.datum : Rappresenta qualsiasi oggetto dati.

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

Archi

Un arco a è composto da tre attributi:

• a.fromNode : questo è un nodo , come definito sopra.

• a.toNode : Questo è un nodo , come definito sopra.

• a.datum : rappresenta qualsiasi oggetto dati.

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

L'insieme degli archi nel digrafo rappresenta una relazione binaria sui nodi nel digrafo . L'esistenza dell'arco a implica che esista una relazione tra a.fromNode e a.toNode .

In un grafo diretto (al contrario di un grafo non orientato), l'esistenza di una relazione tra a.fromNode e a.toNode non implica che esista una relazione simile tra a.toNode e a.fromNode .

Questo perché in un grafo non orientato, la relazione espressa non è necessariamente simmetrica.

Grafici

Nodi e archi possono essere utilizzati per definire un digraph , ma per comodità, negli algoritmi seguenti, un digraph verrà rappresentato utilizzando come dizionario.

Ecco un metodo che può convertire la rappresentazione del grafico sopra in una rappresentazione del dizionario simile a questa:

 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

Ed eccone un altro che lo converte in un dizionario di dizionari, altra operazione che sarà 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

Quando l'articolo parla di un digrafo rappresentato da un dizionario, menzionerà G_as_dict per riferirsi ad esso.

A volte è utile recuperare un nodo da un digrafo G da esso attraverso il suo uid (identificatore univoco):

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

Nella definizione dei grafi, le persone a volte usano i termini nodo e vertice per riferirsi allo stesso concetto; lo stesso vale per i termini arco e spigolo.

Due rappresentazioni di grafici popolari in Python sono questa che usa dizionari e un'altra che usa oggetti per rappresentare grafici. La rappresentazione in questo post è una via di mezzo tra queste due rappresentazioni comunemente usate.

Questa è la mia rappresentazione del digramma . Ce ne sono molti simili, ma questo è mio.

Passeggiate e Percorsi

Sia S_Arcs una sequenza finita (raccolta ordinata) di archi in un digrafo G tale che se a è un qualsiasi arco in S_Arcs tranne l'ultimo, e b segue a nella sequenza, allora deve esserci un nodo n = a.fromNode in G.setOfNodes tale che a.toNode = b.fromNode .

A partire dal primo arco in S_Arcs e proseguendo fino all'ultimo arco in S_Arcs , raccogli (nell'ordine incontrato) tutti i nodi n come definito sopra tra ogni due archi consecutivi in S_Arcs . Etichetta la raccolta ordinata di nodi raccolti durante questa operazione 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

• Se un nodo appare più di una volta nella sequenza S_Nodes , chiama S_Arcs a Walk on digraph G .

• Altrimenti, chiama S_Arcs un percorso da list(S_Nodes)[0] a list(S_Nodes)[-1] sul digrafo G .

Nodo di origine

Chiama nodo s un nodo sorgente nel digrafo G se s è in G.setOfNodes e G.setOfArcs non contiene alcun arco a tale che 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

Nodo terminale

Chiama il nodo t un nodo terminale nel digrafo G se t è in G.setOfNodes e G.setOfArcs non contiene un arco a tale che 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

Tagli e Tagli

Un cut di un digrafo connesso G è un sottoinsieme di archi di G.setOfArcs che partiziona l'insieme di nodi G.setOfNodes in G . G è connesso se ogni nodo n in G.setOfNodes e ha almeno un arco a in G.setOfArcs tale che n = a.fromNode o n = a.toNode , ma a.fromNode != a.toNode .

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

La definizione sopra si riferisce ad un sottoinsieme di archi , ma può anche definire una partizione dei nodi di G.setOfNodes .

Per le funzioni predecessors_of e successors_of , n è un nodo nell'insieme G.setOfNodes del digraph G e cut è un taglio di 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

Sia cut un taglio del digrafo 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 è un taglio del digrafo G se: (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 è chiamato taglio xy se (x in get_first_part(cut, G)) and (y in get_second_part(cut, G) ) and (x != y) . Quando il nodo x in un cut xy è un nodo sorgente e il nodo y nel taglio xy è un nodo terminale , allora questo taglio è chiamato taglio st .

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

Reti di flusso

È possibile utilizzare un digrafo G per rappresentare una rete di flusso.

Assegna a ogni nodo n , dove n è in G.setOfNodes un n.datum che è un FlowNodeDatum :

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

Assegna a ciascun arco a , dove a è in G.setOfArcs e a.datum che è un FlowArcDatum .

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

flowNodeDatum.flowIn e flowNodeDatum.flowOut sono numeri reali positivi.

flowArcDatum.capacity e flowArcDatum.flow sono numeri reali positivi.

Per ogni nodo nodo n in 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})

Il digrafo G rappresenta ora una rete di flusso .

Il flusso di G si riferisce al a.flow a per tutti gli archi a in G .

Flussi fattibili

Sia il digrafo G rappresentare una rete di flusso .

La rete di flussi rappresentata da G ha flussi ammissibili se:

1. Per ogni nodo n in G.setOfNodes ad eccezione dei nodi sorgente e terminali : n.datum.flowIn = n.datum.flowOut .

2. Per ogni arco a in G.setOfNodes : a.datum.capacity <= a.datum.flow .

La condizione 1 è chiamata vincolo di conservazione.

La condizione 2 è chiamata vincolo di capacità.

Capacità di taglio

La capacità di taglio di un st cut stCut con nodo sorgente s e nodo terminale t di una rete di flusso rappresentata da un digrafo G è:

 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

Taglio di capacità minima

Sia stCut = stCut(s,t,cut) un taglio st di una rete di flusso rappresentata da un digrafo G .

stCut è il taglio di capacità minima della rete di flusso rappresentato da G se non sono presenti altri st cut stCutCompetitor in questa rete di flusso tale che:

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

Strappare i flussi via

Vorrei fare riferimento al digrafo che sarebbe il risultato di prendere un digrafo G e rimuovere tutti i dati di flusso da tutti i nodi in G.setOfNodes e anche tutti gli archi in 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)

Problema di flusso massimo

Una rete di flusso rappresentata come un digrafo G , un nodo sorgente s in G.setOfNodes e un nodo terminale t in G.setOfNodes , G può rappresentare un problema di flusso massimo se:

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

Etichetta questa rappresentazione:

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

Dove sourceNodeUid = s.uid , terminalNodeUid = t.uid e maxFlowProblemStateUid è un identificatore per l'istanza del problema.

Soluzione di flusso massimo

Lascia che maxFlowProblem rappresenti un problema di flusso massimo . La soluzione a maxFlowProblem può essere rappresentata da una rete di flusso rappresentata come un digrafo H .

Digraph H è una soluzione fattibile al problema del flusso massimo su input python maxFlowProblem se:

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

2. H è una rete di flusso e ha flussi ammissibili .

Se oltre a 1 e 2:

1. Non può esserci nessun'altra rete di flusso rappresentata dal digrafo K tale che strip_flows(G) == strip_flows(K) e find_node_by_uid(t.uid,G).flowIn < find_node_by_uid(t.uid,K).flowIn .

Allora H è anche una soluzione ottimale per maxFlowProblem .

In altre parole una soluzione fattibile di portata massima può essere rappresentata da un digrafo , che:

1. È identico al digramma G del corrispondente problema di flusso massimo con l'eccezione che n.datum.flowIn , n.datum.flowOut e a.datum.flow di qualsiasi nodo e arco possono essere diversi.

2. Rappresenta una rete di flusso che ha flussi fattibili .

Inoltre, può rappresentare una soluzione di flusso massimo ottimale se inoltre:

1. Il flowIn per il nodo corrispondente al nodo terminale nel problema del flusso massimo è il più grande possibile (quando le condizioni 1 e 2 sono ancora soddisfatte).

Se il digrafo H rappresenta una soluzione di flusso massimo ammissibile : find_node_by_uid(s.uid,H).flowOut = find_node_by_uid(t.uid,H).flowIn questo segue dal teorema max flow, min cut (discusso di seguito). Informalmente, poiché si presume che H abbia flussi ammissibili, ciò significa che il flusso non può essere né "creato" (tranne nel nodo sorgente s ) né "distrutto" (tranne nel nodo terminale t ) mentre attraversa qualsiasi (altro) nodo ( vincoli di conservazione ).

Poiché un problema di flusso massimo contiene solo un singolo nodo sorgente s e un singolo nodo terminale t , tutto il flusso "creato" in s deve essere "distrutto" in t altrimenti la rete di flusso non ha flussi ammissibili (il vincolo di conservazione sarebbe stato violato ).

Sia il digramma H rappresentare una soluzione di flusso massimo ammissibile ; il valore sopra è chiamato st Flow Value di H .

Permettere:

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

Ciò significa che mfps è uno stato successore di maxFlowProblem , il che significa semplicemente che mfps è esattamente come maxFlowProblem con l'eccezione che i valori di a.flow per archi a in maxFlowProblem.setOfArcs possono essere diversi da a.flow per archi a in 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

Ecco una visualizzazione di un mfps insieme al suo associato maxFlowProb . Ogni arco a nell'immagine ha un'etichetta, queste etichette sono a.datum.flowFrom / a.datum.flowTo , ogni nodo n nell'immagine ha un'etichetta e queste etichette sono n.uid {n.datum.flowIn / a.datum.flowOut} .

st Tagliare il flusso

Lascia che mfps rappresenti un MaxFlowProblemState e che stCut rappresenti un taglio di mfps.G . Il flusso di taglio di stCut è definito:

 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 è la somma dei flussi dalla partizione contenente il nodo sorgente alla partizione contenente il nodo terminale meno la somma dei flussi dalla partizione contenente il nodo terminale alla partizione contenente il nodo sorgente .

Flusso massimo, taglio minimo

Lascia che maxFlowProblem rappresenti un problema di flusso massimo e che la soluzione di maxFlowProblem sia rappresentata da una rete di flusso rappresentata dal digrafo H .

Sia minStCut il taglio di capacità minimo della rete di flusso rappresentato da maxFlowProblem.G .

Poiché nel problema del flusso massimo il flusso ha origine in un solo nodo sorgente e termina in un singolo nodo terminale e, a causa dei vincoli di capacità e di conservazione , sappiamo che tutto il flusso che entra in maxFlowProblem.terminalNodeUid deve attraversare qualsiasi st cut , in particolare deve attraversare minStCut . Questo significa:

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

Risolvere il problema del flusso massimo

L'idea di base per risolvere un problema di portata massima maxFlowProblem è di iniziare con una soluzione di portata massima rappresentata dal digrafo H . Tale punto di partenza può essere H = strip_flows(maxFlowProblem.G) . Il compito è quindi di utilizzare H e con qualche avida modifica dei valori a.datum.flow di alcuni archi a in H.setOfArcs per produrre un'altra soluzione di flusso massimo rappresentata dal digrafo K tale che K non possa ancora rappresentare una rete di flusso con flussi ammissibili e get_flow_value(H, maxFlowProblem) < get_flow_value(K, maxFlowProblem) . Finché questo processo continua, la qualità ( get_flow_value(K, maxFlowProblem) ) della soluzione di flusso massimo rilevata più di recente ( K ) è migliore di qualsiasi altra soluzione di flusso massimo trovata. Se il processo raggiunge un punto in cui sa che nessun altro miglioramento è possibile, il processo può terminare e restituirà la soluzione di flusso massimo ottimale.

La descrizione sopra è generale e salta molte prove, ad esempio se un tale processo è possibile o quanto tempo potrebbe richiedere, fornirò alcuni dettagli in più e l'algoritmo.

Il flusso massimo, teorema di taglio minimo

Dal libro Flows in Networks di Ford e Fulkerson, l'affermazione del teorema max flow, min cut (Teorema 5.1) è:

Per qualsiasi rete, il valore di flusso massimo da s a t è uguale alla capacità di taglio minima di tutti i tagli che separano s e t .

Utilizzando le definizioni in questo post, ciò si traduce in:

La soluzione a un maxFlowProblem rappresentato da una rete di flusso rappresentata dal digrafo H è ottimale se:

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

Mi piace questa dimostrazione del teorema e Wikipedia ne ha un altro.

Il teorema max flow, min cut viene utilizzato per dimostrare la correttezza e la completezza del metodo Ford-Fulkerson .

Darò anche una dimostrazione del teorema nella sezione dopo i cammini aumentanti .

Il metodo Ford-Fulkerson e l'algoritmo di Edmonds-Karp

CLRS definisce il metodo Ford-Fulkerson in questo modo (sezione 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

Grafico residuo

Il Grafico Residuo di una rete di flusso rappresentato dal digrafo G può essere rappresentato come un digrafo 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) restituisce la somma di a.datum.capacity per tutti gli archi nel sottoinsieme di G.setOfArcs dove arco a è nel sottoinsieme se a.fromNode = n e a.toNode = u .

• agg_n_to_u_cap(n,u,G_as_dict) restituisce la somma di a.datum.flow per tutti gli archi nel sottoinsieme di G.setOfArcs dove arco a è nel sottoinsieme se a.fromNode = n e a.toNode = u .

In breve, il grafico residuo G_f rappresenta alcune azioni che possono essere eseguite sul digrafo G .

Ciascuna coppia di nodi n,u in G.setOfNodes della rete di flusso rappresentata dal digrafo G può generare 0, 1 o 2 archi nel grafo residuo G_f di G .

1. La coppia n,u non genera alcun arco in G_f se non c'è un arco a in G.setOfArcs tale che a.fromNode = n e a.toNode = u .

2. La coppia n,u genera l' arco a in G_f.setOfArcs dove a rappresenta un arco etichettato come arco di flusso di spinta da n a 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 coppia n,u genera l' arco a in G_f.setOfArcs dove a rappresenta un arco etichettato come arco pull flow da n a u a = Arc(n,u,datum=ResidualNode(n_to_u_cap_sum - n_to_u_flow_sum)) se n_to_u_flow_sum > 0.0 .

• Ciascun arco di flusso di spinta in G_f.setOfArcs rappresenta l'azione di sommare un totale di x <= n_to_u_cap_sum - n_to_u_flow_sum flusso agli archi nel sottoinsieme di G.setOfArcs dove arco a è nel sottoinsieme se a.fromNode = n e a.toNode = u .

• Ciascun arco di flusso di attrazione in G_f.setOfArcs rappresenta l'azione di sottrarre un totale di x <= n_to_u_flow_sum flusso agli archi nel sottoinsieme di G.setOfArcs dove arco a è nel sottoinsieme se a.fromNode = n e a.toNode = u .

L'esecuzione di una singola azione push o pull da G_f sugli archi applicabili in G potrebbe generare una rete di flusso senza flussi ammissibili perché i vincoli di capacità o di conservazione potrebbero essere violati nella rete di flusso generata.

Ecco una visualizzazione del grafico residuo dell'esempio precedente visualizzazione di una soluzione di flusso massimo , nella visualizzazione ogni arco a rappresenta a.residualCapacity .

Percorso in aumento

Sia maxFlowProblem un problema di flusso massimo e sia G_f = get_residual_graph_of(G) il grafico residuo di maxFlowProblem.G .

Un percorso di augmentingPath per maxFlowProblem è qualsiasi percorso da find_node_by_uid(maxFlowProblem.sourceNode,G_f) a find_node_by_uid(maxFlowProblem.terminalNode,G_f) .

Si scopre che un percorso augmentingPath può essere applicato a una soluzione di flusso massimo rappresentata dal digrafo H generando un'altra soluzione di flusso massimo rappresentata dal digrafo K dove get_flow_value(H, maxFlowProblem) < get_flow_value(K, maxFlowProblem) se H non è ottimale .

Ecco come:

 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

In quanto sopra, TOL è un valore di tolleranza per arrotondare i valori di flusso nella rete. Questo per evitare imprecisioni a cascata dei calcoli in virgola mobile. Quindi, ad esempio, ho usato TOL = 10 per indicare arrotondare a 10 cifre significative.

Sia K = augment(augmentingPath, H) , quindi K rappresenta una soluzione di flusso massimo fattibile per maxFlowProblem . Perché l'affermazione sia vera, la rete di flusso rappresentata da K deve avere flussi ammissibili (non violare il vincolo di capacità o il vincolo di conservazione .

Ecco perché: nel metodo sopra, ogni nodo aggiunto alla nuova rete di flusso rappresentata dal digrafo K è o una copia esatta di un nodo dal digrafo H o un nodo n che ha avuto lo stesso numero aggiunto al suo n.datum.flowIn come suo n.datum.flowOut . Ciò significa che il vincolo di conservazione è soddisfatto in K fintanto che è stato soddisfatto in H . Il vincolo di conservazione è soddisfatto perché controlliamo esplicitamente che ogni nuovo arco a nella rete abbia a.datum.flow <= a.datum.capacity ; quindi, fintanto che gli archi dell'insieme H.setOfArcs che sono stati copiati non modificati in K.setOfArcs non violano il vincolo di capacità , allora K non viola il vincolo di capacità .

È anche vero che get_flow_value(H, maxFlowProblem) < get_flow_value(K, maxFlowProblem) se H non è ottimale .

Ecco perché: Affinché un percorso augmentingPath esista nella rappresentazione digrafica del grafo residuo G_f di un problema di flusso massimo maxFlowProblem , l'ultimo arco a su augmentingPath deve essere un arco 'push' e deve avere a.toNode == maxFlowProblem.terminalNode . Un percorso aumentante è definito come quello che termina al nodo terminale del problema di flusso massimo per il quale è un percorso aumentante . Dalla definizione del grafo residuo , è chiaro che l'ultimo arco in un percorso crescente su quel grafo residuo deve essere un arco 'push' perché qualsiasi arco 'pull' b nel percorso crescente avrà b.toNode == maxFlowProblem.terminalNode e b.fromNode != maxFlowProblem.terminalNode dalla definizione del percorso . Inoltre, dalla definizione di percorso , è chiaro che il nodo terminale viene modificato solo una volta dal metodo di augment . Quindi augment modifica maxFlowProblem.terminalNode.flowIn esattamente una volta e aumenta il valore di maxFlowProblem.terminalNode.flowIn perché l'ultimo arco in augmentingPath deve essere l' arco che causa la modifica in maxFlowProblem.terminalNode.flowIn durante l' augment . Dalla definizione di augment applicata agli archi "push", maxFlowProblem.terminalNode.flowIn può essere solo aumentato, non diminuito.

Alcune prove da Sedgewick e Wayne

Il libro Algorithms, quarta edizione di Robert Sedgewich e Kevin Wayne ha alcune meravigliose e brevi prove (pagine 892-894) che saranno utili. Li ricreerò qui, anche se userò un linguaggio che si adatta alle definizioni precedenti. Le mie etichette per le bozze sono le stesse del libro di Sedgewick.

Proposizione E: Per qualsiasi digrafo H che rappresenta una soluzione di flusso massimo fattibile a un problema di flusso massimo maxFlowProblem , per qualsiasi stCut get_stcut_flow(stCut,H,maxFlowProblem) = get_flow_value(H, maxFlowProblem) .

Dimostrazione: Let stCut=stCut(maxFlowProblem.sourceNode,maxFlowProblem.terminalNode,set([a for a in H.setOfArcs if a.toNode == maxFlowProblem.terminalNode])) . La proposizione E vale per stCut direttamente dalla definizione del valore di flusso st . Supponiamo di voler spostare qualche nodo n dalla partizione s ( get_first_part(stCut.cut, G) ) e nella partizione t (get_second_part(stCut.cut, G)) , per farlo dobbiamo cambiare stCut.cut , che potrebbe modificare stcut_flow = get_stcut_flow(stCut,H,maxFlowProblem) e invalidare la proposta E . Tuttavia, vediamo come cambierà il valore di stcut_flow mentre apportiamo questa modifica. il nodo n è in equilibrio, il che significa che la somma del flusso nel nodo n è uguale alla somma del flusso in uscita da esso (questo è necessario affinché H rappresenti una soluzione ammissibile ). Si noti che tutto il flusso che fa parte dello stcut_flow che entra nel nodo n lo entra dalla partizione s (il flusso che entra nel nodo n dalla partizione t direttamente o indirettamente non sarebbe stato contato nel valore stcut_flow perché si sta dirigendo nella direzione sbagliata in base alla definizione). Inoltre, tutto il flusso in uscita da n finirà (direttamente o indirettamente) nel nodo terminale (dimostrato in precedenza). Quando spostiamo il nodo n nella partizione t, tutto il flusso che entra in n dalla partizione s deve essere aggiunto al nuovo valore stcut_flow ; tuttavia, tutto il flusso in uscita da n deve essere sottratto dal nuovo valore di stcut_flow ; la parte dell'intestazione del flusso direttamente nella partizione a t viene sottratta perché questo flusso è ora interno alla nuova partizione a t e non viene contato come stcut_flow . Anche la parte del flusso da n che si dirige verso i nodi nella partizione s deve essere sottratta da stcut_flow : dopo che n è stato spostato nella partizione t, questi flussi saranno diretti dalla partizione t e nella partizione s e così non devono essere contabilizzati in stcut_flow , poiché questi flussi vengono rimossi l'afflusso nella partizione s e devono essere ridotti della somma di questi flussi e il flusso in uscita dalla partizione s nella partizione t (dove tutti i flussi da st deve finire) deve essere ridotto dello stesso importo. Poiché il nodo n era in equilibrio prima del processo, l'aggiornamento avrà aggiunto lo stesso valore al nuovo valore stcut_flow sottratto, lasciando così la proposizione E vera dopo l'aggiornamento. La validità della proposizione E segue quindi dall'induzione sulla dimensione della partizione t.

Ecco alcuni esempi di reti di flusso per aiutare a visualizzare i casi meno ovvi in ​​cui vale la proposizione E ; nell'immagine, le aree rosse indicano la partizione a s, le aree blu rappresentano la partizione a t e gli archi verdi indicano un taglio a st . Nella seconda immagine, il flusso tra il nodo A e il nodo B aumenta mentre il flusso nel nodo terminale t non cambia.:

Corollario: Nessun valore di flusso di st cut può superare la capacità di qualsiasi st cut .

Proposizione F. (flusso massimo, teorema di taglio minimo): Sia f un flusso st . Le seguenti 3 condizioni sono equivalenti:

1. Esiste un primo taglio la cui portata è uguale al valore della portata f .

2. f è un flusso massimo .

3. Non esiste un cammino crescente rispetto a f .

La condizione 1 implica la condizione 2 per corollario. La condizione 2 implica la condizione 3 perché l'esistenza di un cammino crescente implica l'esistenza di un flusso con valori maggiori, contraddicendo la massimalità di f . La condizione 3 implica la condizione 1: Sia C_s l'insieme di tutti i nodi che possono essere raggiunti da s con un cammino crescente nel grafo residuo . Sia C_t gli archi rimanenti, quindi t deve essere in C_t (secondo la nostra ipotesi). Gli archi che attraversano da C_s a C_t formano quindi un primo taglio che contiene solo archi a dove a.datum.capacity = a.datum.flow o a.datum.flow = 0.0 . Se così fosse, allora i nodi collegati da un arco con capacità residua residua a C_s sarebbero nell'insieme C_s poiché ci sarebbe quindi un percorso crescente da s a tale nodo . Il flusso attraverso il primo taglio è uguale alla capacità del primo taglio (poiché gli archi da C_s a C_t hanno flusso uguale alla capacità) e anche al valore del primo flusso (per la proposizione E ).

Questa affermazione del teorema max flow, min cut implica la precedente affermazione di Flows in Networks.

Corollario (proprietà di integralità): quando le capacità sono interi, esiste un flusso massimo a valore intero e l'algoritmo Ford-Fulkerson lo trova.

Dimostrazione: ogni percorso crescente aumenta il flusso di un intero positivo, il minimo delle capacità inutilizzate negli archi 'push' e i flussi negli archi 'pull', che sono sempre numeri interi positivi.

Ciò giustifica la descrizione del metodo Ford-Fulkerson di CLRS . Il metodo consiste nel continuare a trovare percorsi aumentanti augment applicare l'aumento all'ultima maxFlowSolution soluzioni migliori, fino a quando non c'è più percorso aumentante, il che significa che l'ultima soluzione di flusso massimo è ottimale .

Da Ford-Fulkerson a Edmonds-Karp

Le restanti domande relative alla risoluzione dei problemi di flusso massimo sono:

1. Come dovrebbero essere costruiti i percorsi aumentanti ?

2. Il metodo terminerà se utilizziamo numeri reali e non interi?

3. Quanto tempo ci vorrà per terminare (se lo fa)?

L' algoritmo di Edmonds-Karp specifica che ogni cammino aumentante è costruito da una ricerca in ampiezza (BFS) del grafo residuo ; si scopre che questa decisione del punto 1 sopra forzerà anche l'algoritmo a terminare (punto 2) e consentirà di determinare la complessità asintotica del tempo e dello spazio.

Innanzitutto, ecco un'implementazione 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

Ho usato una deque dal modulo delle raccolte python.

Per rispondere alla domanda 2 sopra, parafraserò un'altra dimostrazione di Sedgewick e Wayne: Proposizione G. Il numero di cammini aumentanti necessari nell'algoritmo Edmonds-Karp con N nodi e A archi è al massimo NA/2 . Dimostrazione: ogni percorso aumentante ha un arco di collo di bottiglia , un arco che viene cancellato dal grafo residuo perché corrisponde o a un arco "push" che viene riempito fino alla capacità o a un arco "pull" attraverso il quale il flusso diventa 0. Ogni volta un arc diventa un arco collo di bottiglia, la lunghezza di qualsiasi percorso crescente che lo attraversa deve aumentare di un fattore 2. Questo perché ogni nodo in un percorso può apparire solo una volta o per niente (dalla definizione di percorso ) poiché i percorsi sono in corso esplorato dal percorso più breve al più lungo ciò significa che almeno un altro nodo deve essere visitato dal percorso successivo che passa attraverso il particolare nodo del collo di bottiglia, il che significa altri 2 archi sul percorso prima di arrivare al nodo . Poiché il percorso crescente è di lunghezza al massimo N , ogni arco può essere su al massimo N/2 cammini aumentanti , e il numero totale di cammini aumentanti è al massimo NA/2 .

L' algoritmo Edmonds-Karp viene eseguito in O(NA^2) . Se durante l'algoritmo verranno esplorati al massimo NA/2 percorsi ed esplorando ogni percorso con BFS è N+A allora il termine più significativo del prodotto e quindi la complessità asintotica è O(NA^2) .

Lascia che mfp sia 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 versione sopra è inefficiente e ha una complessità peggiore di O(NA^2) poiché costruisce una nuova soluzione di flusso massimo e ogni volta un nuovo grafo residuo (piuttosto che modificare i digrafi esistenti man mano che l'algoritmo avanza). Per ottenere una vera soluzione O(NA^2) l'algoritmo deve mantenere sia il digrafo che rappresenta lo stato del problema di flusso massimo sia il grafo residuo associato. Quindi l'algoritmo deve evitare di iterare inutilmente su archi e nodi e aggiornare i loro valori e valori associati nel grafico residuo solo se necessario.

Per scrivere un algoritmo Edmonds Karp più veloce, ho riscritto diversi pezzi di codice da quanto sopra. Spero che l'analisi del codice che ha generato un nuovo digrafo sia stato utile per capire cosa sta succedendo. Nell'algoritmo veloce, utilizzo alcuni nuovi trucchi e strutture dati Python che non voglio approfondire in dettaglio. Dirò che a.fromNode e a.toNode sono ora trattati come stringhe e uid per nodes . Per questo codice, lascia che mfps sia 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

Ecco una visualizzazione di come questo algoritmo risolve la rete di flusso di esempio dall'alto. La visualizzazione mostra i passaggi così come sono riflessi nel digrafo G che rappresenta la rete di flusso più aggiornata e come si riflettono nel grafico residuo di tale rete di flusso. I cammini aumentanti nel grafo residuo sono mostrati come cammini rossi, e il digrafo che rappresenta il problema l'insieme dei nodi e degli archi interessati da un dato cammino aumentante è evidenziato in verde. In ogni caso, evidenzierò le parti del grafico che verranno modificate (in rosso o in verde) e poi mostrerò il grafico dopo le modifiche (solo in nero).

Ecco un'altra visualizzazione di come questo algoritmo risolve una diversa rete di flusso di esempio. Si noti che questo utilizza numeri reali e contiene più archi con gli stessi valori fromNode e toNode .

**Si noti inoltre che poiché gli archi con un ResidualDatum 'pull' possono far parte del percorso di aumento, i nodi interessati nel DiGraph della rete volata _potrebbero non trovarsi su un percorso in G! .

Grafici bipartiti

Supponiamo di avere un digrafo G , G è bipartito se è possibile partizionare i nodi in G.setOfNodes in due insiemi ( part_1 e part_2 ) tali che per qualsiasi arco a in G.setOfArcs non può essere vero che a.fromNode in part_1 e a.toNode nella part_1 . Inoltre non può essere vero che a.fromNode in part_2 e a.toNode in part_2 .

In altre parole G è bipartito se può essere partizionato in due insiemi di nodi in modo tale che ogni arco debba connettere un nodo in un insieme a un nodo nell'altro insieme.

Test bipartito

Supponiamo di avere un digrafo G , vogliamo verificare se è bipartito . Possiamo farlo in O(|G.setOfNodes|+|G.setOfArcs|) colorando avidamente il grafico in due colori.

Innanzitutto, dobbiamo generare un nuovo digrafo H . Questo grafico avrà lo stesso insieme di nodi di G , ma avrà più archi di G Ogni arco a in G creerà 2 archi in H ; il primo arco sarà identico a a e il secondo arco inverte il direttore di 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)

Abbinamenti e abbinamenti massimi

Supponiamo di avere un digrafo G e che la matching sia un sottoinsieme di archi da G.setOfArcs . matching è una corrispondenza se per due archi qualsiasi a in matching : len( b len(frozenset( {a.fromNode} ).union( {a.toNode} ).union( {b.fromNode} ).union( {b.toNode} )) == 4 . In altre parole, non esistono due archi in una corrispondenza che condividono un nodo .

Matching matching , è un abbinamento massimo se non ci sono altre corrispondenze alt_matching in G tali che len(matching) < len(alt_matching) . In altre parole, la matching è una corrispondenza massima se è il più grande insieme di archi di G.setOfArcs che soddisfa ancora la definizione di corrispondenza (l'aggiunta di qualsiasi arco non già nella corrispondenza interromperà la definizione della corrispondenza ).

Una corrispondenza di matching massima è una corrispondenza perfetta se ogni per il nodo n in G.setOfArcs esiste un arco a nella matching dove a.fromNode == n or a.toNode == n .

Massima corrispondenza bipartita

Un massimo abbinamento bipartito è un massimo abbinamento su un digrafo G che è bipartito .

Dato che G è bipartito , il problema di trovare un massimo abbinamento bipartito può essere trasformato in un problema di massimo flusso risolvibile con l'algoritmo di Edmonds-Karp e quindi il massimo abbinamento bipartito può essere recuperato dalla soluzione del problema di massimo flusso .

Sia bipartition una bipartizione di G .

Per fare ciò, ho bisogno di generare un nuovo digraph ( H ) con alcuni nuovi nodi ( H.setOfNodes ) e alcuni nuovi archi ( H.setOfArcs ). H.setOfNodes contiene tutti i nodi in G.setOfNodes e altri due nodess , s (un nodo di origine ) t (un nodo terminale ).

H.setOfArcs conterrà un arco per ogni G.setOfArcs . Se un arco a è in G.setOfArcs e a.fromNode è in bipartition.firstPart e a.toNode è in bipartition.secondPart , includere a in H (aggiungendo un FlowArcDatum(1,0) ).

Se a.fromNode è in bipartition.secondPart e a.toNode è in bipartition.firstPart , includi Arc(a.toNode,a.fromNode,FlowArcDatum(1,0)) in H.setOfArcs .

La definizione di un grafo bipartito assicura che nessun arco connetta i nodi in cui entrambi i nodi si trovano nella stessa partizione. H.setOfArcs contiene anche un arco dal nodo s a ciascun nodo in bipartition.firstPart . Infine, H.setOfArcs contiene un arco per ogni nodo in bipartition.secondPart al nodo t . a.datum.capacity = 1 per tutti a in H.setOfArcs .

Prima partiziona i nodi in G.setOfNodes i due insiemi disgiunti ( part1 e part2 ) in modo tale che nessun arco in G.setOfArcs sia diretto da un insieme allo stesso insieme (questa partizione è possibile perché G è bipartito ). Quindi, aggiungi tutti gli archi in G.setOfArcs che sono diretti da part1 a part2 in H.setOfArcs . Quindi crea un singolo nodo sorgente s e un singolo nodo terminale t e crea altri archi

Quindi, costruisci 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

Copertura del nodo minima

Una copertura di nodo in un digrafo G è un insieme di nodi ( cover ) da G.setOfNodes tale che per qualsiasi arco a di G.setOfArcs questo deve essere vero: (a.fromNode in cover) or (a.toNode in cover) .

Una copertura minima di nodi è l'insieme più piccolo possibile di nodi nel grafico che è ancora una copertura di nodi . Il teorema di Konig afferma che in un grafo bipartito , la dimensione della corrispondenza massima su quel grafo è uguale alla dimensione della copertura minima del nodo e suggerisce come la copertura del nodo può essere recuperata da una corrispondenza massima :

Supponiamo di avere la bipartizione bipartition e la corrispondenza massima matching . Definisci un nuovo digrafo H , H.setOfNodes=G.setOfNodes , gli archi in H.setOfArcs sono l'unione di due insiemi.

Il primo set è archi a in matching , con la modifica che se a.fromNode in bipartition.firstPart e a.toNode in bipartition.secondPart quindi a.fromNode e a.toNode vengono scambiati nell'arco creato, danno a tali archi un a.datum.inMatching=True per indicare che sono stati derivati ​​da archi in una corrispondenza .

Il secondo set è archi a NOT in matching , con la modifica che se a.fromNode in bipartition.secondPart e a.toNode in bipartition.firstPart quindi a.fromNode e a.toNode vengono scambiati nell'arco creato (dare a tali archi a a.datum.inMatching=False ).

Quindi, esegui una ricerca in profondità (DFS) a partire da ciascun nodo n in bipartition.firstPart che non è né n == a.fromNode né n == a.toNodes per qualsiasi arco a in matching . Durante il DFS, alcuni nodi vengono visitati e altri no (memorizzare queste informazioni in un campo n.datum.visited ). La copertura minima del nodo è l'unione dei nodi {a.fromNode for a in H.setOfArcs if ( (a.datum.inMatching) and (a.fromNode.datum.visited) ) } e dei nodi {a.fromNode for a in H.setOfArcs if (a.datum.inMatching) and (not a.toNode.datum.visited)} .

Si può dimostrare che ciò porta da una corrispondenza massima a un nodo minimo coperto da una dimostrazione per contraddizione, prendi un arco a che presumibilmente non era coperto e considera tutti e quattro i casi relativi all'appartenenza di a.fromNode e a.toNode (sia come toNode che fromNode ) a qualsiasi arco in corrispondenza di matching . Ogni caso porta a una contraddizione dovuta all'ordine in cui DFS visita i nodi e al fatto che la matching è una corrispondenza massima .

Supponiamo di avere una funzione per eseguire questi passaggi e restituire l'insieme di nodi che comprende la copertura minima del nodo quando viene fornito il digrafo G e la corrispondenza massima di 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

Il problema dell'assegnazione lineare

Il problema dell'assegnazione lineare consiste nel trovare un abbinamento di peso massimo in un grafo bipartito pesato.

Problemi come quello all'inizio di questo post possono essere espressi come un problema di assegnazione lineare . Dato un insieme di lavoratori, un insieme di compiti e una funzione che indica la redditività dell'assegnazione di un lavoratore a un compito, vogliamo massimizzare la somma di tutti gli incarichi che svolgiamo; questo è un problema di assegnazione lineare .

Assumiamo che il numero di compiti e lavoratori sia uguale, anche se mostrerò che questo presupposto è facile da rimuovere. Nell'implementazione, rappresento i pesi dell'arco con un attributo a.datum.weight per un arco a .

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

Algoritmo di Kuhn-Munkres

L'algoritmo di Kuhn-Munkres risolve il problema dell'assegnazione lineare . Una buona implementazione può richiedere tempo O(N^{4}) , (dove N è il numero di nodi nel digrafo che rappresentano il problema). Un'implementazione più facile da spiegare prende O(N^{5}) (per una versione che rigenera DiGraphs ) e O(N^{4}) for (per una versione che mantiene DiGraphs ). Questo è simile alle due diverse implementazioni dell'algoritmo Edmonds-Karp .

Per questa descrizione, sto lavorando solo con grafici bipartiti completi (quelli in cui metà dei nodi si trovano in una parte della bipartizione e l'altra metà nella seconda parte). Nel lavoratore, motivazione del compito significa che ci sono tanti lavoratori quanti sono i compiti.

Questa sembra una condizione significativa (e se questi set non fossero uguali!), ma è facile risolvere questo problema; Parlo di come farlo nell'ultima sezione.

La versione dell'algoritmo qui descritta utilizza l'utile concetto di archi a peso zero . Sfortunatamente, questo concetto ha senso solo quando stiamo risolvendo una minimizzazione (se invece di massimizzare i profitti dei nostri incarichi di lavoro stessimo invece riducendo al minimo il costo di tali incarichi).

Fortunatamente, è facile trasformare un problema di assegnazione lineare massima in un problema di assegnazione lineare minima impostando ciascuno dei pesi dell'arco a su Ma.datum.weight dove M=max({a.datum.weight for a in G.setOfArcs}) . La soluzione del problema di massimizzazione originale sarà identica alla soluzione del problema di minimizzazione dopo che i pesi dell'arco sono stati modificati. Quindi, per il resto, supponiamo di apportare questa modifica.

L' algoritmo di Kuhn-Munkres risolve la corrispondenza del peso minimo in un grafo bipartito ponderato mediante una sequenza di corrispondenze massime nei grafi bipartiti non ponderati. Se a troviamo un abbinamento perfetto sulla rappresentazione in digrafo del problema di assegnazione lineare , e se il peso di ogni arco nell'abbinamento è zero, allora abbiamo trovato l' abbinamento di peso minimo poiché questo abbinamento suggerisce che tutti i nodi nel digrafo sono stati abbinato a un arco con il costo più basso possibile (nessun costo può essere inferiore a 0, in base alle definizioni precedenti).

Nessun altro arco può essere aggiunto alla corrispondenza (perché tutti i nodi sono già abbinati) e nessun arco deve essere rimosso dalla corrispondenza poiché qualsiasi arco di sostituzione possibile avrà almeno un valore di peso altrettanto grande.

Se troviamo una corrispondenza massima del sottografo di G che contiene solo archi di peso zero e non è una corrispondenza perfetta , non abbiamo una soluzione completa (poiché la corrispondenza non è perfetta ). Tuttavia, possiamo produrre un nuovo digrafo H modificando i pesi degli archi in G.setOfArcs in modo che appaiano nuovi archi di peso 0 e la soluzione ottima di H sia la stessa della soluzione ottima di G . Poiché garantiamo che ad ogni iterazione venga prodotto almeno un arco di peso zero , garantiamo che arriveremo a una corrispondenza perfetta in non più di |G.setOfNodes|^{2}=N^{2} di tali iterazioni.

Si supponga che in bipartition bipartition , bipartition.firstPart contenga nodi che rappresentano i lavoratori e bipartition.secondPart rappresenti nodi che rappresentano le attività.

L'algoritmo inizia generando un nuovo digrafo H . H.setOfNodes = G.setOfNodes . Alcuni archi in H sono generati dai nodi n in bipartition.firstPart . Ciascuno di questi nodi n genera un arco b in H.setOfArcs per ogni arco a in bipartition.G.setOfArcs dove a.fromNode = n o a.toNode = n , b=Arc(a.fromNode, a.toNode, a.datum.weight - z) dove z=min(x.datum.weight for x in G.setOfArcs if ( (x.fromNode == n) or (x.toNode == n) )) .

Più archi in H vengono generati dai nodi n in bipartition.secondPart . Ciascuno di questi nodi n genera un arco b in H.setOfArcs per ogni arco a in bipartition.G.setOfArcs dove a.fromNode = n o a.toNode = n , b=Arc(a.fromNode, a.toNode, ArcWeightDatum(a.datum.weight - z)) dove z=min(x.datum.weight for x in G.setOfArcs if ( (x.fromNode == n) or (x.toNode == n) )) .

KMA: Quindi, forma un nuovo digrafo K composto solo dagli archi di peso zero da H e dai nodi incidenti su quegli archi . Forma una bipartition sui nodi in K , quindi usa solve_mbm( bipartition ) per ottenere una corrispondenza massima ( matching ) su K . Se la matching è una corrispondenza perfetta in H (gli archi nella matching sono incidenti su tutti i nodi in H.setOfNodes ), la matching è una soluzione ottimale al problema di assegnazione lineare .

Altrimenti, se la matching non è perfetta , genera la copertura minima del nodo di K usando node_cover = get_min_node_cover(matching, bipartition(K)) . Quindi, definisci z=min({a.datum.weight for a in H.setOfArcs if a not in node_cover}) . Definisci 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)} . Il simbolo != nell'espressione precedente funge da operatore XOR. Quindi arcs = arcs1.union(arcs2.union(arcs3)) . Successivamente, H=DiGraph(nodes,arcs) . Torna a l'etichetta KMA.L'algoritmo continua fino a quando non viene prodotto un abbinamento perfetto.Questo abbinamento è anche la soluzione al problema di assegnazione lineare .

 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

Questa implementazione è O(N^{5}) perché genera una nuova corrispondenza massima di matching ad ogni iterazione; simile alle due precedenti implementazioni di Edmonds-Karp, questo algoritmo può essere modificato in modo che tenga traccia della corrispondenza e lo adatti in modo intelligente a ogni iterazione. Al termine, la complessità diventa O(N^{4}) . Una versione più avanzata e più recente di questo algoritmo (che richiede alcune strutture dati più avanzate) può essere eseguita in O(N^{3}) . I dettagli sia dell'implementazione più semplice di cui sopra che dell'implementazione più avanzata possono essere trovati in questo post che ha motivato questo post sul blog.

Nessuna delle operazioni sui pesi degli archi modifica l'assegnazione finale restituita dall'algoritmo. Ecco perché: poiché i nostri grafici di input sono sempre grafici bipartiti completi , una soluzione deve mappare ogni nodo in una partizione a un altro nodo nella seconda partizione, tramite l' arco tra questi due nodi . Si noti che le operazioni eseguite sui pesi degli archi non cambiano mai l'ordine (ordinato per peso) degli archi incidenti su un nodo particolare.

Pertanto, quando l'algoritmo termina con una perfetta corrispondenza bipartita completa , a ciascun nodo viene assegnato un arco di peso zero , poiché l'ordine relativo degli archi da quel nodo non è cambiato durante l'algoritmo e poiché un arco di peso zero è l'arco più economico possibile e la perfetta corrispondenza bipartita completa garantisce l'esistenza di un tale arco per ciascun nodo . Ciò significa che la soluzione generata è effettivamente la stessa della soluzione del problema di assegnazione lineare originale senza alcuna modifica dei pesi dell'arco .

Compiti sbilanciati

Sembra che l'algoritmo sia piuttosto limitato poiché, come descritto, opera solo su grafi bipartiti completi (quelli in cui metà dei nodi si trovano in una parte della bipartizione e l'altra metà nella seconda parte). Nel lavoratore, la motivazione del compito significa che ci sono tanti lavoratori quanti sono i compiti (sembra piuttosto limitante).

Tuttavia, esiste una facile trasformazione che rimuove questa restrizione. Supponiamo che ci siano meno lavoratori rispetto alle attività, aggiungiamo alcuni lavoratori fittizi (abbastanza per rendere il grafo risultante un grafo bipartito completo ). 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.

Ulteriori letture sul blog di Toptal Engineering:

• Graph Data Science con Python/NetworkX