_mwb_matching.c
上传用户:gzelex
上传日期:2007-01-07
资源大小:707k
文件大小:7k
开发平台:

MultiPlatform

  1. /*******************************************************************************
  2. +
  3. +  LEDA-R  3.2.3
  4. +
  5. +  _mwb_matching.c
  6. +
  7. +  Copyright (c) 1995  by  Max-Planck-Institut fuer Informatik
  8. +  Im Stadtwald, 66123 Saarbruecken, Germany     
  9. +  All rights reserved.
  10. *******************************************************************************/
  11. /*******************************************************************************
  12. *                                                                              *
  13. *  MAX_WEIGHT_BIPARTITE_MATCHING  (maximum weight bipartite matching)          *
  14. *        *
  15. *  modified 01/94, M.Paul        *
  16. *                                                                              *
  17. *******************************************************************************/
  18. #include <LEDA/graph_alg.h>
  19. #include <LEDA/node_pq.h>
  20. static list<edge> compute_MWBM( graph& G,
  21.                           const list<node>& A, const list<node>& B,
  22.                           const edge_array<num_type>& weight );
  23. list<edge> MAX_WEIGHT_BIPARTITE_MATCHING( graph& G,
  24.                  const list<node>& A, const list<node>& B,
  25.                  const edge_array<num_type>& weight )
  26. {
  27.   list<edge> mwm;
  28.   if( ! (A.empty() || B.empty()) ) 
  29.   {
  30.     /* change the graph such that A becomes the node set with
  31.        smaller cardinality (if necessary)
  32.      */
  33.     if( B.size() < A.size() ) 
  34.     { G.rev();
  35.       mwm = compute_MWBM( G, B, A, weight );
  36.       G.rev();
  37.      }
  38.     else
  39.       mwm = compute_MWBM( G, A, B, weight );
  40.   }
  41.   return mwm;
  42. }
  43. /* return the list of possible endpoints of augmenting paths which
  44.    only need a minimal change of the dual function u
  45.  */
  46. static list<node> dijkstra( const graph& G, 
  47.     const node_list& free,
  48.                             const node_array<num_type>& u,
  49.                             const edge_array<num_type>& weight,
  50.                                   num_type maxval,
  51.                                   node_array<num_type>& dist,
  52.                                   node_array<edge>& pred,
  53.                             const node_array<int>& side)
  54.   node_pq<num_type> PQ(G);
  55.   num_type a_min = maxval, // minimal distance to an endpoint of 
  56.       b_min = maxval, // an augmenting path in A resp. B
  57.            dv, dw, c, uv;
  58.   list<node> alist, blist; // the lists of minimal endpoints
  59.   node v, w;
  60.   edge e;
  61.   forall(v,free) {  // initialize distances and PQ
  62.     dist[v] = 0;
  63.     PQ.insert(v,0);
  64.     if( u[v]<=a_min ) { // compute initial nodelist for A
  65.       if( u[v]<a_min ) {
  66.         alist.clear(); 
  67.         a_min = u[v]; 
  68.       }
  69.       alist.append(v);
  70.     }
  71.   }
  72.   while( !PQ.empty() ) { 
  73.     v = PQ.del_min(); // dist[v] is the final distance of v
  74.     dv = dist[v];
  75.     uv = u[v];
  76.     /* if v is in A, update the nodelist for A
  77.      */
  78.     if( !side[v] && uv+dv<=a_min ) {
  79.       if( uv+dv<a_min ) {
  80.         alist.clear(); 
  81.         a_min = uv+dv;
  82.       }
  83.       alist.append(v);
  84.     }
  85.     /* stop if the actual distance becomes greater than the 
  86.        distance to an endpoint of an augmenting path already seen
  87.      */
  88.     if( a_min<dv || b_min<dv ) 
  89.       break;
  90.       
  91.     /* if v is in B, update the nodelist for B
  92.      */
  93.     if( side[v] && !G.outdeg(v) ) {
  94.       b_min = dv;
  95.       blist.append(v);
  96.     }
  97.     /* perform one step of Dijkstra's algorithm
  98.      */
  99.     forall_adj_edges( e, v ) { 
  100.       w = G.target(e);
  101.       dw = dist[w];
  102.       c = dv + uv + u[w] - weight[e];
  103.       if (c < dv) c = dv;  /* rounding errors: u[v]+u[w]-weight[e]
  104.                                                might become negative */
  105.       if( c < dw ) { 
  106.         if( dw == maxval )
  107.            PQ.insert(w,c);
  108.         else
  109.            PQ.decrease_p(w,c);
  110.         dist[w] = c;
  111.         pred[w] = e;
  112.       }
  113.     }
  114.   }
  115.   /* return the appropriate list of nodes
  116.    */
  117.   if( a_min==b_min )
  118.     alist.conc( blist );
  119.   
  120.   return a_min<=b_min ? alist : blist;
  121. }
  122. /* augment matching in one pass along paths of length 1 and 3 
  123. */
  124. static int mwbm_heuristic( graph& G, 
  125.            const list<node>& A,
  126.            const edge_array<num_type>& weight, 
  127.            node_array<num_type>& u)
  128.   node v, w;
  129.   edge e, e2, eb;
  130.   int n = 0;
  131.   num_type max, max2, we;
  132.   node_array<edge> sec_edge(G,0);
  133.   forall( v, A ) { // for all nodes in A ...
  134.     max2 = 0; max = 0; eb = 0;
  135.     forall_adj_edges( e, v ) { // compute edges with largest
  136.       we = weight[e]-u[target(e)]; // and second largest slack
  137.       if( we >= max2 ) {
  138.         if( we >= max ) {
  139.   max2 = max;  e2 = eb;
  140.           max = we;  eb = e;
  141.         }
  142. else {
  143.   max2 = we;  e2 = e;
  144. }
  145.       }
  146.     }
  147.     if( eb ) {
  148.       w = target(eb);
  149.       if( !G.outdeg(w) ) { // match eb and change u[] such
  150.         sec_edge[v] = e2; // that e2 also gets slack 0
  151.         u[v] = max2;
  152.         u[w] = max-max2;
  153.         G.rev_edge(eb);
  154.         n++;
  155.       }
  156.       else { // try to augment matching along
  157.         u[v] = max; // path of length 3 given by sec_edge[]
  158. e2 = G.first_adj_edge(w);
  159. e = sec_edge[target(e2)];
  160. if( e && !G.outdeg(target(e)) ) {
  161.   G.rev_edge(e); G.rev_edge(e2); G.rev_edge(eb);
  162.   n++;
  163. }
  164.       }
  165.     }
  166.   }
  167.   return n;
  168. }
  169. /* compute a maximum weight matching in G
  170.  */
  171. static list<edge> compute_MWBM( graph& G, 
  172. const list<node>& A,
  173. const list<node>& B,
  174. const edge_array<num_type>& weight )
  175. {
  176.   num_type MAX, MIN;
  177.   Max_Value(MAX); Min_Value(MIN);
  178.   
  179.   num_type d_min;
  180.   node v, v_min;
  181.   edge e;
  182.   node_array<edge>     pred(G,nil);
  183.   node_array<num_type> dist(G,MAX), u(G,0);
  184.   node_array<int>      side(G,1), used(G,0);
  185.   int mark=0; // nodes v with used[v]<mark are unused
  186.   node_list free;
  187.   list<node> vlist;
  188.   mwbm_heuristic(G,A,weight,u);
  189.   forall(v,A) {
  190.      if( G.indeg(v)==0 && G.outdeg(v) ) {
  191.        free.append(v);
  192.        dist[v] = 0;
  193.       }
  194.       side[v]=0;
  195.   }
  196.   
  197.   while( !free.empty() ) {
  198.     mark++; // mark all nodes as unused
  199.     /* compute the list of possible endpoints
  200.      */
  201.     vlist = dijkstra(G,free,u,weight,MAX,dist,pred,side);
  202.     forall( v_min, vlist ) {
  203.       v=v_min; 
  204.       /* if v_min is not the first node of the list, check if the
  205.          augmenting path to v_min contains an used node
  206.        */
  207.       if( v_min!=vlist.head() ) {
  208.         e=pred[v];
  209.         while( e && used[v]<mark ) { 
  210.           v = source(e);
  211.           e = pred[v];
  212.         }
  213.       }
  214.       /* if all nodes are unused, augment the matching along the path
  215.        */
  216.       if( used[v]<mark ) {
  217.         v = v_min; e = pred[v];
  218.         while( e ) { 
  219.           used[v]=mark;
  220.           v = source(e);
  221.           G.rev_edge(e);
  222.           e = pred[v];
  223.         }
  224.         free.del(v);
  225.         used[v]=mark;
  226.       }
  227.     }
  228.     /* change the dual function
  229.      */
  230.     v_min = vlist.head();
  231.     d_min = side[v_min] ? dist[v_min] : u[v_min]+dist[v_min];
  232.     
  233.     forall(v,A) { 
  234.       if (d_min > dist[v]) 
  235.         u[v] -= d_min-dist[v];
  236.       dist[v] = MAX;
  237.     }
  238.   
  239.     forall(v,B) { 
  240.       if (d_min > dist[v]) 
  241.         u[v] += d_min-dist[v];
  242.       dist[v] = MAX;
  243.     }
  244.   }
  245.   /* compute the result
  246.    */
  247.   list<edge> result;
  248.   forall(v,B) 
  249.     forall_adj_edges(e,v) result.append(e);
  250.   forall(e,result) G.rev_edge(e);
  251.   
  252.   return result;
  253. }