Flujo Máximo y el Problema de Asignación Lineal

Aquí hay un problema: su empresa asigna contratistas para cumplir con los contratos. Revisa sus listas y decide qué contratistas están disponibles para un compromiso de un mes y revisa sus contratos disponibles para ver cuáles de ellos son para tareas de un mes. Dado que sabe con qué eficacia cada contratista puede cumplir con cada contrato, ¿cómo asigna contratistas para maximizar la eficacia general de ese mes?

Este es un ejemplo del problema de asignación, y el problema se puede resolver con el algoritmo húngaro clásico.

El algoritmo húngaro (también conocido como algoritmo de Kuhn-Munkres) es un algoritmo de tiempo polinomial que maximiza la coincidencia de peso en un gráfico bipartito ponderado. Aquí, los contratistas y los contratos se pueden modelar como un gráfico bipartito, con su efectividad como los pesos de los bordes entre el contratista y los nodos del contrato.

En este artículo, aprenderá sobre una implementación del algoritmo húngaro que utiliza el algoritmo de Edmonds-Karp para resolver el problema de asignación lineal. También aprenderá cómo el algoritmo de Edmonds-Karp es una ligera modificación del método Ford-Fulkerson y cómo esta modificación es importante.

El problema del flujo máximo

El problema del flujo máximo en sí puede describirse de manera informal como el problema de mover algún fluido o gas a través de una red de tuberías desde una sola fuente a una sola terminal. Esto se hace asumiendo que la presión en la red es suficiente para garantizar que el fluido o el gas no permanezcan en ningún tramo de tubería o accesorio de tubería (aquellos lugares donde se encuentran diferentes longitudes de tubería).

Al hacer ciertos cambios en el gráfico, el problema de asignación se puede convertir en un problema de flujo máximo.

Preliminares

Las ideas necesarias para resolver estos problemas surgen en muchas disciplinas matemáticas y de ingeniería, a menudo conceptos similares se conocen con diferentes nombres y se expresan de diferentes maneras (por ejemplo, matrices de adyacencia y listas de adyacencia). Dado que estas ideas son bastante esotéricas, se eligen con respecto a qué tan generalmente se definirán estos conceptos para cualquier entorno dado.

Este artículo no asumirá ningún conocimiento previo más allá de una pequeña introducción a la teoría de conjuntos.

Las implementaciones en esta publicación representan los problemas como gráficos dirigidos (dígrafo).

dígrafos

Un dígrafo tiene dos atributos, setOfNodes y setOfArcs . Ambos atributos son conjuntos (colecciones desordenadas). En los bloques de código de esta publicación, en realidad estoy usando el conjunto congelado de Python, pero ese detalle no es particularmente importante.

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

(Nota: este y todos los demás códigos de este artículo están escritos en Python 3.6).

Nodos

Un nodo n se compone de dos atributos:

• n.uid : Un identificador único.

  Esto significa que para cualesquiera dos nodos x e y ,

 x.uid != y.uid

• n.datum : Esto representa cualquier objeto de datos.

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

arcos

Un arco a se compone de tres atributos:

• a.fromNode : Este es un nodo , como se define arriba.

• a.toNode : Este es un nodo , como se define arriba.

• a.datum : Esto representa cualquier objeto de datos.

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

El conjunto de arcos en el dígrafo representa una relación binaria en los nodos del dígrafo . La existencia de arc a implica que existe una relación entre a.fromNode y a.toNode .

En un gráfico dirigido (a diferencia de un gráfico no dirigido), la existencia de una relación entre a.fromNode y a.toNode no implica que exista una relación similar entre a.toNode y a.fromNode .

Esto se debe a que en un gráfico no dirigido, la relación que se expresa no es necesariamente simétrica.

dígrafos

Los nodos y los arcos se pueden usar para definir un dígrafo , pero por conveniencia, en los algoritmos a continuación, un dígrafo se representará como un diccionario.

Aquí hay un método que puede convertir la representación gráfica anterior en una representación de diccionario similar a esta:

 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

Y aquí hay otro que lo convierte en un diccionario de diccionarios, otra operación que te será útil:

 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

Cuando el artículo habla de un dígrafo representado por un diccionario, mencionará G_as_dict para referirse a él.

A veces es útil obtener un nodo de un dígrafo G a través de su uid (identificador único):

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

Al definir gráficos, las personas a veces usan los términos nodo y vértice para referirse al mismo concepto; lo mismo ocurre con los términos arco y arista.

Dos representaciones gráficas populares en Python son esta que usa diccionarios y otra que usa objetos para representar gráficos. La representación en esta publicación está en algún lugar entre estas dos representaciones de uso común.

Esta es mi representación digráfica . Hay muchos como este, pero este es el mío.

Paseos y Senderos

Sea S_Arcs una secuencia finita (colección ordenada) de arcos en un dígrafo G tal que si a es cualquier arco en S_Arcs excepto el último, y b sigue a en la secuencia, entonces debe haber un nodo n = a.fromNode en G.setOfNodes tal que a.toNode = b.fromNode .

Comenzando desde el primer arco en S_Arcs y continuando hasta el último arco en S_Arcs , recopile (en el orden en que se encuentran) todos los nodos n como se definió anteriormente entre cada dos arcos consecutivos en S_Arcs . Etiquete la colección ordenada de nodos recopilados durante esta operación S_{Nodes} .

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

• Si algún nodo aparece más de una vez en la secuencia S_Nodes , llame a S_Arcs a Walk en el dígrafo G

• De lo contrario, llama a S_Arcs una ruta desde list(S_Nodes)[0] a list(S_Nodes)[-1] en el dígrafo G .

Nodo de origen

Llame al nodo s un nodo de origen en el dígrafo G si s está en G.setOfNodes y G.setOfArcs no contiene a arco tal que a.toNode = s .

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

Nodo terminal

Llame al nodo t un nodo terminal en el dígrafo G si t está en G.setOfNodes y G.setOfArcs no contiene a arco tal que a.fromNode = t .

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

Cortes y cortes de st

Un corte cut de un dígrafo conectado G es un subconjunto de arcos de G.setOfArcs que divide el conjunto de nodos G.setOfNodes en G . G está conectado si cada nodo n en G.setOfNodes y tiene al menos un arco a en G.setOfArcs tal que n = a.fromNode o n = a.toNode , pero a.fromNode != a.toNode .

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

La definición anterior se refiere a un subconjunto de arcos , pero también puede definir una partición de los nodos de G.setOfNodes .

Para las funciones predecessors_of y successors_of , n es un nodo en el conjunto G.setOfNodes del dígrafo G , y cut es un corte de G :

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

Sea cut un corte del dígrafo 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 es un corte del dígrafo G si: (get_first_part(cut, G).union(get_second_part(cut, G)) == G.setOfNodes) and (len(get_first_part(cut, G).intersect(get_second_part(cut, G))) == 0) cut se llama corte xy si (x in get_first_part(cut, G)) and (y in get_second_part(cut, G) ) and (x != y) . Cuando el nodo x en un cut xy es un nodo fuente y el nodo y en el corte xy es un nodo terminal , entonces este corte se denomina corte st .

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

Redes de flujo

Puede usar un dígrafo G para representar una red de flujo.

Asigne a cada nodo n , donde n está en G.setOfNodes un n.datum que es un FlowNodeDatum :

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

Asigne a cada arco a , donde a está en G.setOfArcs y a.datum que es un FlowArcDatum .

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

flowNodeDatum.flowIn y flowNodeDatum.flowOut son números reales positivos.

flowArcDatum.capacity y flowArcDatum.flow también son números reales positivos.

Para cada nodo nodo n en 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})

El dígrafo G ahora representa una red de flujo .

El flujo de G se refiere al a.flow a para todos los arcos a en G

Flujos factibles

Sea el dígrafo G una red de flujo .

La red de flujo representada por G tiene flujos factibles si:

1. Para cada nodo n en G.setOfNodes excepto para los nodos de origen y los nodos terminales : n.datum.flowIn = n.datum.flowOut .

2. Para cada arco a en G.setOfNodes : a.datum.capacity <= a.datum.flow .

La condición 1 se llama restricción de conservación.

La condición 2 se denomina restricción de capacidad.

Capacidad de corte

La capacidad de corte de un st corte stCut con nodo fuente s y nodo terminal t de una red de flujo representada por un dígrafo G es:

 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

Corte de capacidad mínima

Sea stCut = stCut(s,t,cut) un st cut de una red de flujo representada por un dígrafo G .

stCut es el corte de capacidad mínima de la red de flujo representada por G si no hay otro st cut stCutCompetitor en esta red de flujo tal que:

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

Eliminando los flujos

Me gustaría referirme al dígrafo que sería el resultado de tomar un dígrafo G y eliminar todos los datos de flujo de todos los nodos en G.setOfNodes y también todos los arcos en 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 de flujo máximo

Una red de flujo representada como un dígrafo G , un nodo fuente s en G.setOfNodes y un nodo terminal t en G.setOfNodes , G puede representar un problema de flujo máximo si:

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

Etiqueta esta representación:

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

Donde sourceNodeUid = s.uid , terminalNodeUid = t.uid y maxFlowProblemStateUid es un identificador para la instancia del problema.

Solución de flujo máximo

Sea maxFlowProblem un problema de flujo máximo . La solución a maxFlowProblem se puede representar mediante una red de flujo representada como un dígrafo H .

Digraph H es una solución factible para el problema de flujo máximo en la entrada python maxFlowProblem si:

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

2. H es una red de flujo y tiene flujos factibles .

Si además de 1 y 2:

1. No puede haber otra red de flujo representada por el dígrafo K tal que strip_flows(G) == strip_flows(K) y find_node_by_uid(t.uid,G).flowIn < find_node_by_uid(t.uid,K).flowIn .

Entonces H también es una solución óptima para maxFlowProblem .

En otras palabras, una solución de flujo máximo factible se puede representar mediante un dígrafo , que:

1. Es idéntico al dígrafo G del problema de flujo máximo correspondiente con la excepción de que el n.datum.flowIn , n.datum.flowOut y el a.datum.flow de cualquiera de los nodos y arcos pueden ser diferentes.

2. Representa una red de flujo que tiene flujos factibles .

Y, puede representar una solución óptima de flujo máximo si además:

1. El flowIn de entrada para el nodo correspondiente al nodo terminal en el problema de flujo máximo es tan grande como sea posible (cuando las condiciones 1 y 2 aún se cumplen).

Si el dígrafo H representa una solución de flujo máximo factible : find_node_by_uid(s.uid,H).flowOut = find_node_by_uid(t.uid,H).flowIn esto se deduce del flujo máximo, el teorema de corte mínimo (discutido a continuación). De manera informal, dado que se supone que H tiene flujos factibles, esto significa que el flujo no se puede 'crear' (excepto en el nodo de origen s ) ni 'destruir' (excepto en el nodo terminal t ) al cruzar cualquier (otro) nodo ( restricciones de conservación ).

Dado que un problema de flujo máximo contiene solo un nodo fuente único s y un nodo terminal único t , todo el flujo 'creado' en s debe ser 'destruido' en t o la red de flujo no tiene flujos factibles (se habría violado la restricción de conservación ).

Deje que el dígrafo H represente una solución de flujo máximo factible ; el valor anterior se denomina valor de flujo st de H .

Dejar:

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

Esto significa que mfps es un estado sucesor de maxFlowProblem , lo que significa que mfps es exactamente como maxFlowProblem con la excepción de que los valores de a.flow para arcos a en maxFlowProblem.setOfArcs pueden ser diferentes de a.flow para arcos a en 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

Aquí hay una visualización de un mfps junto con su maxFlowProb asociado. Cada arco a en la imagen tiene una etiqueta, estas etiquetas son a.datum.flowFrom / a.datum.flowTo , cada nodo n en la imagen tiene una etiqueta, y estas etiquetas son n.uid {n.datum.flowIn / a.datum.flowOut} .

Flujo de corte st

Deje que mfps represente un MaxFlowProblemState y deje que stCut represente un corte de mfps.G El flujo de corte de stCut se define:

 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

El primer flujo de corte es la suma de los flujos desde la partición que contiene el nodo de origen hasta la partición que contiene el nodo terminal menos la suma de los flujos desde la partición que contiene el nodo terminal hasta la partición que contiene el nodo de origen .

Flujo máximo, corte mínimo

Deje que maxFlowProblem represente un problema de flujo máximo y deje que la solución de maxFlowProblem se represente mediante una red de flujo representada como dígrafo H .

Sea minStCut el corte de capacidad mínima de la red de flujo representada por maxFlowProblem.G .

Debido a que en el problema de flujo máximo, el flujo se origina en un solo nodo de origen y termina en un solo nodo terminal y, debido a las restricciones de capacidad y las restricciones de conservación , sabemos que todo el flujo que ingresa a maxFlowProblem.terminalNodeUid debe cruzar cualquier corte de st . , en particular debe cruzar minStCut . Esto significa:

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

Resolviendo el Problema de Flujo Máximo

La idea básica para resolver un problema de flujo máximo maxFlowProblem es comenzar con una solución de flujo máximo representada por el dígrafo H . Tal punto de partida puede ser H = strip_flows(maxFlowProblem.G) . La tarea es entonces usar H y mediante alguna modificación codiciosa de los valores de a.datum.flow de algunos arcos a en H.setOfArcs para producir otra solución de flujo máximo representada por el dígrafo K tal que K aún no puede representar una red de flujo con flujos factibles y get_flow_value(H, maxFlowProblem) < get_flow_value(K, maxFlowProblem) . Mientras este proceso continúe, la calidad ( get_flow_value(K, maxFlowProblem) ) de la solución de flujo máximo encontrada más recientemente ( K ) es mejor que cualquier otra solución de flujo máximo que se haya encontrado. Si el proceso llega a un punto en el que sabe que no es posible ninguna otra mejora, el proceso puede terminar y devolverá la solución óptima de flujo máximo .

La descripción anterior es general y omite muchas pruebas, como si dicho proceso es posible o cuánto tiempo puede llevar. Daré algunos detalles más y el algoritmo.

El flujo máximo, teorema de corte mínimo

Del libro Flows in Networks de Ford y Fulkerson, la declaración del flujo máximo, teorema de corte mínimo (Teorema 5.1) es:

Para cualquier red, el valor de flujo máximo de s a t es igual a la capacidad de corte mínima de todos los cortes que separan s y t .

Usando las definiciones en esta publicación, eso se traduce en:

La solución a un maxFlowProblem representado por una red de flujo representada como dígrafo H es óptima si:

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

Me gusta esta prueba del teorema y Wikipedia tiene otra.

El teorema de flujo máximo, corte mínimo se utiliza para probar la corrección y la integridad del método de Ford-Fulkerson .

También daré una prueba del teorema en la sección después de aumentar caminos .

El método de Ford-Fulkerson y el algoritmo de Edmonds-Karp

CLRS define el método Ford-Fulkerson así (sección 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

Gráfico de residuos

El gráfico residual de una red de flujo representada como el dígrafo G se puede representar como un dígrafo 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) devuelve la suma de a.datum.capacity para todos los arcos en el subconjunto de G.setOfArcs donde arc a está en el subconjunto si a.fromNode = n y a.toNode = u .

• agg_n_to_u_cap(n,u,G_as_dict) devuelve la suma de a.datum.flow para todos los arcos en el subconjunto de G.setOfArcs donde arc a está en el subconjunto si a.fromNode = n y a.toNode = u .

Brevemente, el gráfico residual G_f representa ciertas acciones que se pueden realizar en el dígrafo G

Cada par de nodos n,u en G.setOfNodes de la red de flujo representada por el dígrafo G puede generar 0, 1 o 2 arcos en el gráfico residual G_f de G .

1. El par n,u no genera ningún arco en G_f si no hay un arco a en G.setOfArcs tal que a.fromNode = n y a.toNode = u .

2. El par n,u genera el arco a en G_f.setOfArcs donde a representa un arco etiquetado como arco de flujo de empuje de 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. El par n,u genera el arco a en G_f.setOfArcs donde a representa un arco etiquetado como arco de flujo de extracción de n a u a = Arc(n,u,datum=ResidualNode(n_to_u_cap_sum - n_to_u_flow_sum)) si n_to_u_flow_sum > 0.0 .

• Cada arco de flujo de inserción en G_f.setOfArcs representa la acción de agregar un total de x <= n_to_u_cap_sum - n_to_u_flow_sum flujo a los arcos en el subconjunto de G.setOfArcs donde arc a está en el subconjunto si a.fromNode = n y a.toNode = u

• Cada arco de flujo de extracción en G_f.setOfArcs representa la acción de restar un total de flujo x <= n_to_u_flow_sum a arcos en el subconjunto de G.setOfArcs donde arc a está en el subconjunto si a.fromNode = n y a.toNode = u .

Realizar una acción individual de empujar o jalar desde G_f en los arcos aplicables en G podría generar una red de flujo sin flujos factibles porque las restricciones de capacidad o las restricciones de conservación podrían violarse en la red de flujo generada.

Aquí hay una visualización del gráfico residual de la visualización de ejemplo anterior de una solución de flujo máximo , en la visualización, cada arco a representa una a.residualCapacity .

Ruta de aumento

Sea maxFlowProblem un problema de flujo máximo y sea G_f = get_residual_graph_of(G) el gráfico residual de maxFlowProblem.G .

Una ruta de augmentingPath AugmentingPath para maxFlowProblem es cualquier ruta desde find_node_by_uid(maxFlowProblem.sourceNode,G_f) hasta find_node_by_uid(maxFlowProblem.terminalNode,G_f) .

Resulta que se puede aplicar una ruta de augmentingPath a una solución de flujo máximo representada por el dígrafo H generando otra solución de flujo máximo representada por el dígrafo K donde get_flow_value(H, maxFlowProblem) < get_flow_value(K, maxFlowProblem) si H no es óptimo .

Así es cómo:

 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

En lo anterior, TOL es un valor de tolerancia para redondear los valores de flujo en la red. Esto es para evitar la imprecisión en cascada de los cálculos de coma flotante. Entonces, por ejemplo, usé TOL = 10 para significar redondear a 10 dígitos significativos.

Sea K = augment(augmentingPath, H) , entonces K representa una solución factible de flujo máximo para maxFlowProblem . Para que la declaración sea verdadera, la red de flujo representada por K debe tener flujos factibles (no violar la restricción de capacidad o la restricción de conservación) .

He aquí por qué: en el método anterior, cada nodo agregado a la nueva red de flujo representada por el dígrafo K es una copia exacta de un nodo del dígrafo H o un nodo n al que se le ha agregado el mismo número a su n.datum.flowIn como es n.datum.flowOut . Esto significa que la restricción de conservación se cumple en K siempre que se satisfaga en H . La restricción de conservación se cumple porque verificamos explícitamente que cualquier nuevo arco a en la red tenga a.datum.flow <= a.datum.capacity ; por lo tanto, siempre que los arcos del conjunto H.setOfArcs que se copiaron sin modificar en K.setOfArcs no violen la restricción de capacidad , entonces K no viola la restricción de capacidad .

También es cierto que get_flow_value(H, maxFlowProblem) < get_flow_value(K, maxFlowProblem) si H no es óptimo .

Este es el motivo: para que exista una ruta de augmentingPath en la representación de dígrafo del gráfico residual G_f de un problema de flujo máximo maxFlowProblem , entonces el último arco a en la ruta de augmentingPath debe ser un arco de "empuje" y debe tener a.toNode == maxFlowProblem.terminalNode . Una ruta de aumento se define como aquella que termina en el nodo terminal del problema de flujo máximo para el cual es una ruta de aumento. A partir de la definición del gráfico residual , está claro que el último arco en una ruta de aumento en ese gráfico residual debe ser un arco de "empuje" porque cualquier arco b de "tracción" en la ruta de aumento tendrá b.toNode == maxFlowProblem.terminalNode y b.fromNode != maxFlowProblem.terminalNode de la definición de ruta . Además, a partir de la definición de ruta , está claro que el nodo terminal solo se modifica una vez mediante el método de augment . Por lo tanto, el maxFlowProblem.terminalNode.flowIn augment una vez y aumenta el valor de maxFlowProblem.terminalNode.flowIn porque el último arco en la ruta de augmentingPath debe ser el arco que provoca la modificación en maxFlowProblem.terminalNode.flowIn durante el augment . A partir de la definición de augment tal como se aplica a los arcos 'push', maxFlowProblem.terminalNode.flowIn solo se puede aumentar, no disminuir.

Algunas pruebas de Sedgewick y Wayne

El libro Algorithms, cuarta edición de Robert Sedgewich y Kevin Wayne tiene algunas pruebas cortas y maravillosas (páginas 892-894) que serán útiles. Los recrearé aquí, aunque usaré un lenguaje que encaje con las definiciones anteriores. Mis etiquetas para las pruebas son las mismas que en el libro de Sedgewick.

Proposición E: para cualquier dígrafo H que represente una solución factible de flujo máximo para un problema de flujo máximo maxFlowProblem , para cualquier stCut get_stcut_flow(stCut,H,maxFlowProblem) = get_flow_value(H, maxFlowProblem) .

Prueba: Let stCut=stCut(maxFlowProblem.sourceNode,maxFlowProblem.terminalNode,set([a for a in H.setOfArcs if a.toNode == maxFlowProblem.terminalNode])) . La proposición E se cumple para stCut directamente de la definición de st flow value . Supongamos que allí deseamos mover algún nodo n de la partición s ( get_first_part(stCut.cut, G) ) y dentro de la partición t (get_second_part(stCut.cut, G)) , para hacerlo necesitamos cambiar stCut.cut , que podría cambiar stcut_flow = get_stcut_flow(stCut,H,maxFlowProblem) e invalidar la proposición E . Sin embargo, veamos cómo cambiará el valor de stcut_flow a medida que hagamos este cambio. el nodo n está en equilibrio, lo que significa que la suma del flujo hacia el nodo n es igual a la suma del flujo que sale (esto es necesario para que H represente una solución factible ). Tenga en cuenta que todo el flujo que forma parte del nodo de entrada n de stcut_flow ingresa desde la partición s (el flujo que ingresa al nodo n desde la partición t, ya sea directa o indirectamente, no se habría contado en el valor de stcut_flow porque se dirige en la dirección incorrecta según la definición). Además, todo el flujo que sale de n eventualmente (directa o indirectamente) fluirá hacia el nodo terminal (probado anteriormente). Cuando movemos el nodo n a la partición t, todo el flujo que ingresa a n desde la partición s debe agregarse al nuevo valor de stcut_flow ; sin embargo, todo el flujo que sale de n debe restarse del nuevo valor de stcut_flow ; la parte del flujo que se dirige directamente a la partición t se resta porque este flujo ahora es interno a la nueva partición t y no se cuenta como stcut_flow . La parte del flujo de n que se dirige a los nodos en la partición s también se debe restar de stcut_flow : después de que n se mueva a la partición t, estos flujos se dirigirán desde la partición t a la partición s y así no debe tenerse en cuenta en stcut_flow , ya que estos flujos se eliminan del flujo de entrada a la partición s y deben reducirse por la suma de estos flujos, y el flujo de salida de la partición s a la partición t (donde todos los flujos de st debe terminar) debe reducirse en la misma cantidad. Como el nodo n estaba en equilibrio antes del proceso, la actualización habrá agregado el mismo valor al nuevo valor de stcut_flow que restó, dejando así la proposición E verdadera después de la actualización. La validez de la proposición E se sigue entonces de la inducción sobre el tamaño de la partición t.

Aquí hay algunos ejemplos de redes de flujo para ayudar a visualizar los casos menos obvios donde se cumple la proposición E ; en la imagen, las áreas rojas indican la partición s, las áreas azules representan la partición t y los arcos verdes indican un corte st . En la segunda imagen, el flujo entre el nodo A y el nodo B aumenta mientras que el flujo hacia el nodo terminal t no cambia:

Corolario: Ningún valor de flujo de corte de st puede exceder la capacidad de cualquier corte de st .

Proposición F. (flujo máximo, teorema de corte mínimo): Sea f un flujo st . Las siguientes 3 condiciones son equivalentes:

1. Existe un st corte cuya capacidad es igual al valor del caudal f .

2. f es un flujo máximo .

3. No hay camino de aumento con respecto a f .

La condición 1 implica la condición 2 por el corolario. La condición 2 implica la condición 3 porque la existencia de un camino creciente implica la existencia de un flujo con valores mayores, contradiciendo la maximalidad de f . La condición 3 implica la condición 1: Sea C_s el conjunto de todos los nodos a los que se puede llegar desde s con un camino creciente en el grafo residual . Sean C_t los arcos restantes, entonces t debe estar en C_t (según nuestra suposición). Los arcos que cruzan de C_s a C_t luego forman un corte st que contiene solo arcos a donde a.datum.capacity = a.datum.flow o a.datum.flow = 0.0 . Si esto fuera de otro modo, los nodos conectados por un arco con capacidad residual remanente a C_s estarían en el conjunto C_s ya que entonces habría un camino creciente desde s hasta dicho nodo . El flujo a través del st corte es igual a la capacidad del st corte (ya que los arcos de C_s a C_t tienen un flujo igual a la capacidad) y también al valor del st flujo (por la proposición E ).

Esta afirmación del teorema de flujo máximo, corte mínimo implica la afirmación anterior de Flows in Networks.

Corolario (propiedad de integralidad): cuando las capacidades son números enteros, existe un flujo máximo de valor entero y el algoritmo de Ford-Fulkerson lo encuentra.

Prueba: cada ruta de aumento aumenta el flujo en un número entero positivo, el mínimo de las capacidades no utilizadas en los arcos de 'empuje' y los flujos en los arcos de 'tirón', todos los cuales son siempre números enteros positivos.

Esto justifica la descripción del método Ford-Fulkerson de CLRS . El método consiste en seguir encontrando rutas de aumento y aplicar el maxFlowSolution augment encontrar mejores soluciones, hasta que no haya más rutas de aumento, lo que significa que la última solución de flujo máximo es óptima .

De Ford-Fulkerson a Edmonds-Karp

Las preguntas restantes con respecto a la solución de problemas de flujo máximo son:

1. ¿Cómo deben construirse los caminos de aumento ?

2. ¿Terminará el método si usamos números reales y no enteros?

3. ¿Cuánto tiempo llevará terminar (si lo hace)?

El algoritmo de Edmonds-Karp especifica que cada ruta de aumento se construye mediante una búsqueda en amplitud (BFS) del gráfico residual ; resulta que esta decisión del punto 1 anterior también forzará la terminación del algoritmo (punto 2) y permitirá determinar la complejidad asintótica de tiempo y espacio.

Primero, aquí hay una implementación de 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

Usé un deque del módulo de colecciones de python.

Para responder a la pregunta 2 anterior, parafrasearé otra prueba de Sedgewick y Wayne: Proposición G. La cantidad de rutas de aumento necesarias en el algoritmo de Edmonds-Karp con N nodos y A arcos es como máximo NA/2 . Prueba: cada ruta de aumento tiene un arco de cuello de botella , un arco que se elimina del gráfico residual porque corresponde a un arco de "empuje" que se llena al máximo o a un arco de "tirón" a través del cual el flujo se convierte en 0. Cada vez que un arco se convierte en un arco de cuello de botella, la longitud de cualquier ruta de aumento a través de él debe aumentar en un factor de 2. Esto se debe a que cada nodo en una ruta puede aparecer solo una vez o no aparecer en absoluto (según la definición de ruta ) ya que las rutas están siendo explorado del camino más corto al más largo, lo que significa que al menos un nodo más debe ser visitado por el siguiente camino que pasa por el nodo de cuello de botella en particular, lo que significa 2 arcos adicionales en el camino antes de que lleguemos al nodo . Dado que el trayecto de aumento tiene una longitud máxima de N , cada arco puede estar en un máximo de N/2 trayectos de aumento, y el número total de trayectos de aumento es como máximo NA/2 .

El algoritmo de Edmonds-Karp se ejecuta en O(NA^2) . Si se explorarán como máximo NA/2 rutas durante el algoritmo y explorar cada ruta con BFS es N+A , entonces el término más significativo del producto y, por lo tanto, la complejidad asintótica es O(NA^2) .

Sea mfp un maxFlowProblemState .

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

La versión anterior es ineficiente y tiene peor complejidad que O(NA^2) ya que construye una nueva solución de flujo máximo y un nuevo gráfico residual cada vez (en lugar de modificar los dígrafos existentes a medida que avanza el algoritmo). Para llegar a una verdadera solución O(NA^2) , el algoritmo debe mantener tanto el dígrafo que representa el estado del problema de flujo máximo como su gráfico residual asociado. Por lo tanto, el algoritmo debe evitar iterar sobre arcos y nodos innecesariamente y actualizar sus valores y valores asociados en el gráfico residual solo cuando sea necesario.

Para escribir un algoritmo de Edmonds Karp más rápido, reescribí varios fragmentos de código del anterior. Espero que revisar el código que generó un nuevo dígrafo haya sido útil para comprender lo que está sucediendo. En el algoritmo rápido, uso algunos trucos nuevos y estructuras de datos de Python que no quiero analizar en detalle. Diré que a.fromNode y a.toNode ahora se tratan como cadenas y uids a nodos . Para este código, deje que mfps sea un maxFlowProblemState

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

Aquí hay una visualización de cómo este algoritmo resuelve la red de flujo de ejemplo de arriba. La visualización muestra los pasos tal como se reflejan en el dígrafo G que representa la red de flujo más actualizada y tal como se reflejan en el gráfico residual de esa red de flujo. Las rutas de aumento en el gráfico residual se muestran como rutas rojas, y el dígrafo que representa el problema, el conjunto de nodos y arcos afectados por una ruta de aumento dada, se resalta en verde. En cada caso, resaltaré las partes del gráfico que se cambiarán (en rojo o verde) y luego mostraré el gráfico después de los cambios (solo en negro).

Aquí hay otra visualización de cómo este algoritmo resuelve una red de flujo de ejemplo diferente. Tenga en cuenta que este usa números reales y contiene varios arcos con los mismos valores fromNode y toNode .

**También tenga en cuenta que debido a que los arcos con un ResidualDatum 'pull' pueden ser parte de la ruta de aumento, los nodos afectados en el DiGraph de la red volada _pueden no estar en una ruta en G! .

Gráficos bipartitos

Supongamos que tenemos un dígrafo G , G es bipartito si es posible dividir los nodos en G.setOfNodes en dos conjuntos ( part_1 y part_2 ) de modo que para cualquier arco a en G.setOfArcs no puede ser cierto que a.fromNode en part_1 y a.toNode en part_1 . Tampoco puede ser cierto que a.fromNode en part_2 y a.toNode en part_2 .

En otras palabras, G es bipartito si se puede dividir en dos conjuntos de nodos de modo que cada arco debe conectar un nodo de un conjunto con un nodo del otro conjunto.

Prueba bipartita

Supongamos que tenemos un dígrafo G , queremos probar si es bipartito . Podemos hacer esto en O(|G.setOfNodes|+|G.setOfArcs|) coloreando con avidez el gráfico en dos colores.

Primero, necesitamos generar un nuevo dígrafo H . Este gráfico tendrá el mismo conjunto de nodos que G , pero tendrá más arcos que G . Cada arco a en G creará 2 arcos en H ; el primer arco será idéntico a a , y el segundo arco invierte el director de a ( b = Arc(a.toNode,a.fromNode,a.datum) ).

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

Coincidencias y Máximo de Coincidencias

Supongamos que tenemos un dígrafo G y matching es un subconjunto de arcos de G.setOfArcs . matching es una coincidencia si para cualquiera de los dos arcos a y b matching : len(frozenset( {a.fromNode} ).union( {a.toNode} ).union( {b.fromNode} ).union( {b.toNode} )) == 4 . En otras palabras, no hay dos arcos en una coincidencia que compartan un nodo .

Matching matching , es una coincidencia máxima si no hay otra coincidencia alt_matching en G tal que len(matching) < len(alt_matching) . En otras palabras, matching es una coincidencia máxima si es el conjunto de arcos más grande de G.setOfArcs que todavía satisface la definición de coincidencia (la adición de cualquier arco que no esté ya en la coincidencia romperá la definición de coincidencia ).

Una matching máxima es una coincidencia perfecta si cada nodo n en G.setOfArcs existe un arco a en matching donde a.fromNode == n or a.toNode == n .

Coincidencia bipartita máxima

Una coincidencia bipartita máxima es una coincidencia máxima en un dígrafo G que es bipartito .

Dado que G es bipartita , el problema de encontrar una coincidencia bipartita máxima se puede transformar en un problema de flujo máximo que se pueda resolver con el algoritmo de Edmonds-Karp y luego se puede recuperar la coincidencia bipartita máxima de la solución al problema de flujo máximo .

Sea bipartition una bipartición de G .

Para hacer esto, necesito generar un nuevo dígrafo ( H ) con algunos nodos nuevos ( H.setOfNodes ) y algunos arcos nuevos ( H.setOfArcs ). H.setOfNodes contiene todos los nodos en G.setOfNodes y dos nodos más, s (un nodo fuente ) y t (un nodo terminal ).

H.setOfArcs contendrá un arco para cada G.setOfArcs . Si un arco a está en G.setOfArcs y a.fromNode está en bipartition.firstPart y a.toNode está en bipartition.secondPart , entonces incluya a en H (agregando un FlowArcDatum(1,0) ).

Si a.fromNode está en bipartition.secondPart y a.toNode está en bipartition.firstPart , incluya Arc(a.toNode,a.fromNode,FlowArcDatum(1,0)) en H.setOfArcs .

La definición de un gráfico bipartito asegura que ningún arco conecte ningún nodo donde ambos nodos estén en la misma partición. H.setOfArcs también contiene un arco desde el nodo s hasta cada nodo en bipartition.firstPart . Finalmente, H.setOfArcs contiene un arco de cada nodo en bipartition.secondPart al nodo t . a.datum.capacity = 1 para todo a en H.setOfArcs .

Primero divida los nodos en G.setOfNodes los dos conjuntos disjuntos ( part1 y part2 ) de modo que ningún arco en G.setOfArcs se dirija de un conjunto al mismo conjunto (esta partición es posible porque G es bipartito ). A continuación, agregue todos los arcos en G.setOfArcs que se dirigen desde part1 a part2 en H.setOfArcs . Luego cree un solo nodo de origen s y un solo nodo terminal t y cree algunos arcos más

Luego, construya 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

Cobertura de nodo mínimo

Una cubierta de nodo en un dígrafo G es un conjunto de nodos ( cover ) de G.setOfNodes tal que para cualquier arco a de G.setOfArcs esto debe ser cierto: (a.fromNode in cover) or (a.toNode in cover) .

Una cobertura mínima de nodos es el conjunto más pequeño posible de nodos en el gráfico que sigue siendo una cobertura de nodos . El teorema de Konig establece que en un gráfico bipartito , el tamaño de la coincidencia máxima en ese gráfico es igual al tamaño de la cobertura de nodo mínima , y ​​sugiere cómo se puede recuperar la cobertura de nodo de una coincidencia máxima :

Supongamos que tenemos la bipartición bipartition y la máxima coincidencia matching . Defina un nuevo dígrafo H , H.setOfNodes=G.setOfNodes , los arcos en H.setOfArcs son la unión de dos conjuntos.

El primer conjunto son los arcos a en matching , con el cambio de que si a.fromNode in bipartition.firstPart y a.toNode in bipartition.secondPart , entonces a.fromNode y a.toNode se intercambian en el arco creado, dan a dichos arcos un a.datum.inMatching=True atributo para indicar que se derivaron de arcos en una coincidencia .

El segundo conjunto son arcos a NOT en matching , con el cambio de que si a.fromNode in bipartition.secondPart y a.toNode in bipartition.firstPart entonces a.fromNode y a.toNode se intercambian en el arco creado (da a tales arcos a.datum.inMatching=False ).

A continuación, ejecute una búsqueda en profundidad (DFS) a partir de cada nodo n en bipartition.firstPart que no es ni n == a.fromNode ni n == a.toNodes para ningún arco a en matching . Durante el DFS, se visitan algunos nodos y otros no (almacene esta información en un campo n.datum.visited ). La cobertura mínima de nodos es la unión de los nodos {a.fromNode for a in H.setOfArcs if ( (a.datum.inMatching) and (a.fromNode.datum.visited) ) } y los nodos {a.fromNode for a in H.setOfArcs if (a.datum.inMatching) and (not a.toNode.datum.visited)} .

Se puede demostrar que esto lleva de una coincidencia máxima a una cobertura de nodo mínima mediante una prueba por contradicción, tome un arco a que supuestamente no estaba cubierto y considere los cuatro casos con respecto a si a.fromNode y a.toNode pertenecen (ya sea como toNode o fromNode ) a cualquier arco en coincidencia matching . Cada caso conduce a una contradicción debido al orden en que DFS visita los nodos y al hecho de que la matching es una coincidencia máxima .

Supongamos que tenemos una función para ejecutar estos pasos y devolver el conjunto de nodos que comprende la cobertura mínima de nodos cuando se le da el dígrafo G y la matching máxima :

 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

El problema de asignación lineal

El problema de asignación lineal consiste en encontrar una coincidencia de peso máximo en un gráfico bipartito ponderado.

Problemas como el que se encuentra al comienzo de esta publicación se pueden expresar como un problema de asignación lineal . Dado un conjunto de trabajadores, un conjunto de tareas y una función que indica la rentabilidad de la asignación de un trabajador a una tarea, queremos maximizar la suma de todas las asignaciones que hacemos; este es un problema de asignación lineal .

Suponga que el número de tareas y trabajadores es igual, aunque mostraré que esta suposición es fácil de eliminar. En la implementación, represento pesos de arco con un atributo a.datum.weight para un arco a .

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

Algoritmo de Kuhn-Munkres

El Algoritmo de Kuhn-Munkres resuelve el problema de asignación lineal . Una buena implementación puede llevar O(N^{4}) tiempo (donde N es el número de nodos en el dígrafo que representa el problema). Una implementación que es más fácil de explicar toma O(N^{5}) (para una versión que regenera DiGraphs ) y O(N^{4}) para (para una versión que mantiene DiGraphs ). Esto es similar a las dos implementaciones diferentes del algoritmo de Edmonds-Karp .

Para esta descripción, solo estoy trabajando con grafos bipartitos completos (aquellos donde la mitad de los nodos están en una parte de la bipartición y la otra mitad en la segunda parte). En el trabajador, la motivación por tareas esto significa que hay tantos trabajadores como tareas.

Esto parece una condición importante (¡qué pasa si estos conjuntos no son iguales!) pero es fácil solucionar este problema; Hablo de cómo hacerlo en la última sección.

La versión del algoritmo descrita aquí utiliza el concepto útil de arcos de peso cero . Desafortunadamente, este concepto solo tiene sentido cuando estamos resolviendo una minimización (si en lugar de maximizar las ganancias de nuestras asignaciones de tareas de trabajadores, en cambio, minimicáramos el costo de tales asignaciones).

Afortunadamente, es fácil convertir un problema de asignación lineal máxima en un problema de asignación lineal mínima configurando cada arco a pesos a Ma.datum.weight donde M=max({a.datum.weight for a in G.setOfArcs}) . La solución al problema de maximización original será idéntica a la solución del problema de minimización después de cambiar los pesos del arco . Entonces, para el resto, supongamos que hacemos este cambio.

El algoritmo de Kuhn-Munkres resuelve la coincidencia de peso mínimo en un gráfico bipartito ponderado mediante una secuencia de coincidencias máximas en gráficos bipartitos no ponderados. Si a encontramos una coincidencia perfecta en la representación del dígrafo del problema de asignación lineal , y si el peso de cada arco en la coincidencia es cero, entonces hemos encontrado la coincidencia de peso mínimo ya que esta coincidencia sugiere que todos los nodos en el dígrafo han sido emparejado por un arco con el costo más bajo posible (ningún costo puede ser menor que 0, según definiciones anteriores).

No se pueden agregar otros arcos a la coincidencia (porque todos los nodos ya están emparejados) y no se debe eliminar ningún arco de la coincidencia porque cualquier posible arco de reemplazo tendrá al menos un valor de peso igual de grande.

Si encontramos una coincidencia máxima del subgrafo de G que contiene solo arcos de peso cero , y no es una coincidencia perfecta , no tenemos una solución completa (ya que la coincidencia no es perfecta ). Sin embargo, podemos producir un nuevo dígrafo H cambiando los pesos de los arcos en G.setOfArcs de manera que aparezcan nuevos arcos de peso 0 y la solución óptima de H sea la misma que la solución óptima de G . Como garantizamos que se produce al menos un arco de peso cero en cada iteración, garantizamos que llegaremos a una coincidencia perfecta en no más de |G.setOfNodes|^{2}=N^{2} iteraciones.

Supongamos que en bipartition bipartition , bipartition.firstPart contiene nodos que representan trabajadores y bipartition.secondPart representa nodos que representan tareas.

El algoritmo comienza generando un nuevo dígrafo H . H.setOfNodes = G.setOfNodes . Algunos arcos en H se generan a partir de nodos n en bipartition.firstPart . Cada nodo n genera un arco b en H.setOfArcs para cada arco a en bipartition.G.setOfArcs donde a.fromNode = n o a.toNode = n , b=Arc(a.fromNode, a.toNode, a.datum.weight - z) donde z=min(x.datum.weight for x in G.setOfArcs if ( (x.fromNode == n) or (x.toNode == n) )) .

Se generan más arcos en H a partir de los nodos n en bipartition.secondPart . Cada nodo n genera un arco b en H.setOfArcs para cada arco a en bipartition.G.setOfArcs donde a.fromNode = n o a.toNode = n , b=Arc(a.fromNode, a.toNode, ArcWeightDatum(a.datum.weight - z)) donde z=min(x.datum.weight for x in G.setOfArcs if ( (x.fromNode == n) or (x.toNode == n) )) .

KMA: A continuación, forme un nuevo dígrafo K compuesto solo por los arcos de peso cero de H y los nodos incidentes en esos arcos . Forme una bipartition en los nodos en K , luego use solve_mbm( bipartition ) para obtener una coincidencia máxima ( matching ) en K . Si matching es un emparejamiento perfecto en H (los arcos en el matching inciden en todos los nodos en H.setOfNodes ), entonces el matching es una solución óptima para el problema de asignación lineal .

De lo contrario, si matching no es perfecta , genere la cobertura de nodo mínima de K usando node_cover = get_min_node_cover(matching, bipartition(K)) . A continuación, defina z=min({a.datum.weight for a in H.setOfArcs if a not in node_cover}) . Defina 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)} . El símbolo != en la expresión anterior actúa como un operador XOR. Entonces arcs = arcs1.union(arcs2.union(arcs3)) . A continuación, H=DiGraph(nodes,arcs) . Vuelva a la etiqueta KMA . El algoritmo continúa hasta que se produce una coincidencia perfecta . Esta coincidencia es también la solución al problema de asignación lineal .

 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

Esta implementación es O(N^{5}) porque genera una nueva matching máxima en cada iteración; Al igual que las dos implementaciones anteriores de Edmonds-Karp, este algoritmo se puede modificar para que realice un seguimiento de la coincidencia y lo adapte de manera inteligente a cada iteración. Cuando se hace esto, la complejidad se convierte en O(N^{4}) . Una versión más avanzada y más reciente de este algoritmo (que requiere algunas estructuras de datos más avanzadas) puede ejecutarse en O(N^{3}) . Los detalles tanto de la implementación más simple anterior como de la implementación más avanzada se pueden encontrar en esta publicación que motivó esta publicación de blog.

Ninguna de las operaciones sobre los pesos del arco modifica la asignación final devuelta por el algoritmo. He aquí por qué: dado que nuestros gráficos de entrada son siempre gráficos bipartitos completos, una solución debe asignar cada nodo en una partición a otro nodo en la segunda partición, a través del arco entre estos dos nodos . Tenga en cuenta que las operaciones realizadas en los pesos de los arcos nunca cambian el orden (ordenado por peso) de los arcos que inciden en un nodo en particular.

Por lo tanto, cuando el algoritmo termina en una coincidencia bipartita completa perfecta , a cada nodo se le asigna un arco de peso cero , ya que el orden relativo de los arcos de ese nodo no ha cambiado durante el algoritmo, y dado que un arco de peso cero es el arco más barato posible y el emparejamiento bipartito completo perfecto garantiza que existe uno de esos arcos para cada nodo . Esto significa que la solución generada es de hecho la misma que la solución del problema de asignación lineal original sin ninguna modificación de los pesos de los arcos .

Asignaciones desequilibradas

Parece que el algoritmo es bastante limitado ya que, como se describe, opera solo en grafos bipartitos completos (aquellos en los que la mitad de los nodos están en una parte de la bipartición y la otra mitad en la segunda parte). En el trabajador, la motivación por tareas esto significa que hay tantos trabajadores como tareas (parece bastante limitante).

Sin embargo, existe una transformación fácil que elimina esta restricción. Supongamos que hay menos trabajadores que tareas, agregamos algunos trabajadores ficticios (suficientes para hacer que el gráfico resultante sea un gráfico 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.

Lecturas adicionales en el blog de ingeniería de Toptal:

• Grafique la ciencia de datos con Python/NetworkX