_mwb_matching.c

数据类型和算法库LEDA 数据类型和算法库LEDA

/*******************************************************************************
+
+  LEDA  3.0
+
+
+  _mwb_matching.c
+
+
+  Copyright (c) 1992  by  Max-Planck-Institut fuer Informatik
+  Im Stadtwald, 6600 Saarbruecken, FRG     
+  All rights reserved.
+ 
*******************************************************************************/


/*******************************************************************************
*                                                                              *
*  MAX_WEIGHT_BIPARTITE_MATCHING  (maximum weight bipartite matching)          *
*                                                                              *
*******************************************************************************/

#include <LEDA/graph_alg.h>
#include <LEDA/prio.h>


#if defined(REAL_NUMBERS)
typedef double num_type;
typedef node_array<double> num_node_array;
typedef edge_array<double> num_edge_array;
#else
typedef int num_type;
typedef node_array<int> num_node_array;
typedef edge_array<int> num_edge_array;
#endif


static void dijkstra(const graph& G, const list<node>&       free,
                                     const num_node_array&   u,
                                     const num_edge_array&   weight,
                                                num_type     maxval,
                                           num_node_array&   dist,
                                           node_array<edge>& pred)
{ priority_queue<node,num_type> PQ;
  node_array<pq_item> I(G);
  node v;
  edge e;

  forall(v,free)
  { dist[v] = 0;
    I[v] = PQ.insert(v,0);
   }

  while (! PQ.empty())
  { v = PQ.del_min();
    num_type dv = dist[v];
    num_type uv = u[v];
    forall_adj_edges(e,v)
    { node w = G.target(e);
      num_type dw = dist[w];
      num_type  c = dv + uv + u[w] - weight[e];

      if (c < dv) c = dv;  /* rounding errors: u[v]+u[w]-weight[e]
                                               might become negative
                           */

      if (c < dw)
      { if (dw == maxval)
           I[w] = PQ.insert(w,c);
        else
           PQ.decrease_inf(I[w],c);
        dist[w] = c;
        pred[w] = e;
       }
     }
   }
}


static int heuristic(graph& G, const list<node>& A,
                               const list<node>& B,
                               const num_edge_array& weight,
                                     num_type maxval, 
                                     num_type minval,
                                     num_node_array& u)
{ node v;
  edge e;

  num_type L1 = minval;
  num_type L2 = maxval;

  /* right side (B):  u_j = max c_ij,  L2  = min u_j
                             i                j
  */

  forall(v,B) u[v] = 0;

  forall_edges(e,G)
  { v = target(e);
    num_type we = weight[e];
    if (we > u[v]) u[v] = we;
   }

  forall(v,B)
  { num_type uv = u[v];
    if (L2 > uv) L2 = uv;
   }


  /* left side (A):  u_i = max (c_ij - u_j),  L1  = max -u_i
                            j                        i
  */

  forall_edges(e,G)
  { v = source(e);
    num_type c = weight[e] - u[target(e)];
    if (c > u[v]) u[v] = c;
   }

  forall(v,A)
  if (G.outdeg(v) > 0)
  { num_type uv = -u[v];
    if (L1 < uv) L1 = uv;
   }


  num_type L = (L1 > L2) ? L2 : L1;   /* L = min(L1,L2) */

  forall(v,A)
  { u[v] += L;
    /* if (u[v] < 0) u[v] = 0;    required in non-maxcard case */
   }

  forall(v,B)  u[v] -= L;


  /* construct subgraph G1 of edges (i,j) with reduced cost p_ij =0 */


  GRAPH<node,edge> G1;

  node_array<node> copy_node(G,nil);

  forall_nodes(v,G) copy_node[v] = G1.new_node(v);

  forall_edges(e,G)
  { node v = source(e);
    node w = target(e);
    if (u[v] + u[w] == weight[e])    /* p_ij = 0 */
      G1.new_edge(copy_node[v],copy_node[w],e);
   }

  list<edge> EL = MAX_CARD_BIPARTITE_MATCHING(G1);

  forall(e,EL) G.rev_edge(G1.inf(e));

  return EL.length();

}




list<edge> MAX_WEIGHT_BIPARTITE_MATCHING(graph& G,
                                         const list<node>& A,
                                         const list<node>& B,
                                         const num_edge_array& weight)
{

node v;
edge e;
edge p;

num_type MAX;
num_type MIN;

Max_Value(MAX);
Min_Value(MIN);

num_node_array  dist(G,MAX);
num_node_array  u(G,MIN);
node_array<edge>    pred(G,nil);

list<node> free;
node_array<list_item> free_item(G,nil);

/* without heuristic

forall_nodes(v,G) u[v] = 0;
forall(v,A)
   forall_adj_edges(e,v) if (weight[e] > u[v]) u[v] = weight[e];
*/

heuristic(G,A,B,weight,MAX,MIN,u);

forall(v,A)
   if (G.indeg(v) == 0)
    { free_item[v] = free.append(v);
      dist[v] = 0;
     }

while (!free.empty())
{

  dijkstra(G,free,u,weight,MAX,dist,pred);

  num_type  d_min = MAX;
  node v_min;

  forall(v,B)
    if (outdeg(v)==0 && dist[v] < d_min)
     { v_min = v;
       d_min = dist[v];
      }


/* if it is not required that the matching is of max. cardinality

  num_type min_A = MAX;
  node a;
  forall(v,A)
    if (dist[v] < MAX && (dist[v]+u[v]) < min_A)
    { a = v;
      min_A = dist[v] + u[v];
     }

  if (min_A < d_min)
  { v_min = a;
    d_min = min_A;
   }
*/

  if (d_min == MAX) break;

  e = pred[v_min];
  p = pred[source(e)];

  while (p)
  { G.rev_edge(e);
    e = p;
    p = pred[source(p)];
   }

   free.del(free_item[source(e)]);

   G.rev_edge(e);


  forall(v,A)
  { num_type dv = dist[v];
    if (d_min > dv) u[v] -= d_min-dv;
    dist[v] = MAX;
   }

  forall(v,B)
  { num_type dv = dist[v];
    if (d_min > dv) u[v] += d_min-dv;
    dist[v] = MAX;
   }


}

list<edge> result;
forall(v,B) forall_adj_edges(e,v) result.append(e);

forall(e,result) G.rev_edge(e);

return result;

}