mwb_matching.t

A Library of Efficient Data Types and Algorithms,封装了常用的ADT及其相关算法的软件包

          { dist[b] = db; pred[b] = e;
            PQ.decrease_p(b,db); 
          }
        }

 }
    }

  }

 
 while (!RA.empty() )
 { node a = RA.pop();
   pred[a] = nil;
   NT pot_change = Delta - dist[a];
   if (pot_change <= 0 ) continue;
   pot[a] = pot[a] - pot_change;
 }

 while (!RB.empty() )
 { node b = RB.pop();
   pred[b] = nil;
   if (PQ.member(b)) PQ.del(b);
   NT pot_change = Delta - dist[b];
   if (pot_change <= 0 ) continue;
   pot[b] = pot[b] + pot_change;
 }

}


template <class NT>
__temp_func_inline
void mwbm_heuristic(graph& G, const list<node>& A, 
                    const edge_array<NT>& c, node_array<NT>& pot,
                    node_array<bool>& free)
{ 
  node a, b;  edge e, e2, eb; 
  node_array<edge> sec_edge(G,nil);

  forall( a, A )                      
  { NT max2 = 0; NT max = 0; eb = e2 = nil;
    // compute edges with largest and second largest slack
    forall_adj_edges( e, a )           
    { NT we = c[e] - pot[target(e)];     
      if ( we >= max2 ) 
      { if( we >= max )
        { max2 = max;  e2 = eb;
          max = we;    eb = e;
        }
        else 
        { max2 = we;   e2 = e;
        }
      }
    }

    if( eb ) 
    { b = target(eb);
      if( free[b] )                    
      { // match eb and change pot[] to make slack of e2 zero
        sec_edge[a] = e2;                
        pot[a] = max2;
        pot[b] = max-max2;
        G.rev_edge(eb);
        free[a] = free[b] = false;
      }
      else                              
      { // try to augment matching along 
        // path of length 3 given by sec_edge[]
        pot[a] = max;                  
        e2 = G.first_adj_edge(b);
        e = sec_edge[target(e2)];
        if( e && G.outdeg(target(e)) == 0 ) 
        { free[a] = free[G.target(e)] = false;
          G.rev_edge(e); G.rev_edge(e2); G.rev_edge(eb);
        }
      }
    }
    else pot[a] = 0;
  }
}
  

static int which_heuristic = 2;

template <class NT>
__temp_func_inline
list<edge> MAX_WEIGHT_BIPARTITE_MATCHING_T(graph& G,
                          const list<node>& A, const list<node>& B,   
                          const edge_array<NT>& c, node_array<NT>& pot)
{ node a,b,v; edge e;

  list<edge> result;
  forall_nodes(v,G) pot[v] = 0;

  if (G.number_of_edges() == 0 ) return result;
  
  // check that all edges are directed from A to B
  forall(b,B) assert(G.outdeg(b) == 0);

  node_array<bool> free(G,true);

  node_array<edge> pred(G,nil);
  node_array<NT>   dist(G,0);  
  node_pq<NT>      PQ(G);

  switch (which_heuristic) 
  { case 0: { // naive heuristic 
             NT C = 0;
             forall_edges(e,G) if (c[e] > C) C = c[e];
             forall(a,A)       pot[a] = C;
             break;
           }

    case 1: { // simple heuristic
             
             forall(a,A) 
             { edge e_max = nil; NT C_max = 0;
               forall_adj_edges(e,a)
                 if (c[e] >  C_max) { e_max = e; C_max = c[e]; }
               pot[a] = C_max;  
               if ( e_max != nil && free[b = G.target(e_max)] )
               { G.rev_edge(e_max);
                 free[a] = free[b] = false;
               }
             }
                 
             break;
           }

    default: { // refined heuristic
             mwbm_heuristic( G, A, c, pot, free);
             break;
           }
  }

  forall(a,A)
    if (free[a]) augment(G,a,c,pot,free,pred,dist,PQ);

  forall(b,B)
    { forall_out_edges(e,b) result.append(e); }

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

  return result;
}  



template <class NT>
__temp_func_inline
void check_invariants_max_weight_assignment(const graph& G, 
const edge_array<NT>& c,
        const node_array<NT>& pot, const list<node>& A,
        const list<node>& B, const node_array<edge>& pred,node a0)
{ node a,b,v; edge e;
  forall(b,B)
  { assert(G.outdeg(b) <= 1);  assert(pred[b] == nil);
    assert(pot[b] >= 0);
    assert(G.outdeg(b) == 1 || pot[b] == 0);
    e = G.first_adj_edge(b);
    if (e) 
    { a = G.target(e); 
      assert(pot[a] + pot[b] == c[e]);
    }
  }
  a = a0;
  while (a)
  { assert(G.indeg(a) <= 1);  assert(pred[a] == nil);
 
    forall_adj_edges(e,a)
    { b = G.target(e);
      assert(pot[a] + pot[b] >= c[e]);
    }
    a = G.pred_node(a);
  }
}



LEDA_END_NAMESPACE
#include <LEDA/stack.h>
LEDA_BEGIN_NAMESPACE

template <class NT>
__temp_func_inline
bool max_weight_assignment_augment(graph& G, 
                    node a, const edge_array<NT>& c,
                    node_array<NT>& pot, node_array<bool>& free, 
                    node_array<edge>& pred, node_array<NT>& dist, 
                    node_pq<NT>& PQ)
{ 
  dist[a] = 0;

  //node best_node_in_A = a; 
  NT   minA           = pot[a];
  NT   Delta;

  // mod. 12/2000 by GS
  NT   upper_bound = 0;
  bool upper_bound_is_defined = false;  // upper_bound = +\infty 

  stack<node> RA;  RA.push(a);
  stack<node> RB;

  node a1 = a; edge e;
  
  forall_adj_edges(e,a1)
  { node b = G.target(e);
    NT db = dist[a1] + (pot[a1] + pot[b] - c[e]); 

    // mod. 12/2000 by GS
    if (PRUNE && upper_bound_is_defined && db >= upper_bound) continue;
    if (free[b]) { upper_bound = db; upper_bound_is_defined = true; }

    if ( pred[b] == nil ) 
    { dist[b] = db; pred[b] = e; RB.push(b);  
      PQ.insert(b,db); 
    }
    else
    if ( db < dist[b] )
    { dist[b] = db; pred[b] = e;
      PQ.decrease_p(b,db); 
    }
  }



  while ( true )
  { node b; NT db;
    if (PQ.empty()) { return false; }
    else { b = PQ.del_min(); db = dist[b]; }

    if ( free[b] )
    { Delta = db;
      
      augment_path_to(G,b,pred);
      free[a] = free[b] = false;
      break;
 
    }
    else
    {
     e = G.first_adj_edge(b); 
     node a1 = G.target(e);
     pred[a1] = e; RA.push(a1);
     dist[a1] = db; 
           
     if (db + pot[a1] < minA)
     { //best_node_in_A = a1;
       minA = db + pot[a1];
     }

     
     forall_adj_edges(e,a1)
     { node b = G.target(e);
       NT db = dist[a1] + (pot[a1] + pot[b] - c[e]); 

       // mod. 12/2000 by GS
       if (PRUNE && upper_bound_is_defined && db >= upper_bound) continue;
       if (free[b]) { upper_bound = db; upper_bound_is_defined = true; }

       if ( pred[b] == nil ) 
       { dist[b] = db; pred[b] = e; RB.push(b);  
         PQ.insert(b,db); 
       }
       else
       if ( db < dist[b] )
       { dist[b] = db; pred[b] = e;
         PQ.decrease_p(b,db); 
       }
     }

 }
  }

 
 while (!RA.empty() )
 { node a = RA.pop();
   pred[a] = nil;
   NT pot_change = Delta - dist[a];
   if (pot_change <= 0 ) continue;
   pot[a] = pot[a] - pot_change;
 }

 while (!RB.empty() )
 { node b = RB.pop();
   pred[b] = nil;
   if (PQ.member(b)) PQ.del(b);
   NT pot_change = Delta - dist[b];
   if (pot_change <= 0 ) continue;
   pot[b] = pot[b] + pot_change;
 }


 return true;
}



template <class NT>
__temp_func_inline
list<edge> MAX_WEIGHT_ASSIGNMENT_T(graph& G, 
                          const list<node>& A, 
                          const list<node>& B,  
                          const edge_array<NT>& c,  
                          node_array<NT>& pot)
{ node a,b; edge e;
  
  // check that all edges are directed from A to B
  forall(b,B) assert(G.outdeg(b) == 0);

  node_array<bool> free(G,true);

  node_array<edge> pred(G,nil);
  node_array<NT>   dist(G,0);  
  node_pq<NT>      PQ(G);

  forall(b,B) pot[b] = 0;

  switch (which_heuristic) 
  { case 0: { // no heuristic 
             NT C = 0;
             forall_edges(e,G) if (c[e] > C) C = c[e];
             forall(a,A)       pot[a] = C;
             break;
           }

    case 1: { // simple heuristic
             
             forall(a,A) 
             { edge e_max = nil; NT C_max = 0;
               forall_adj_edges(e,a)
                 if (c[e] >  C_max) { e_max = e; C_max = c[e]; }
               pot[a] = C_max;  
               if ( e_max != nil && free[b = G.target(e_max)] )
               { G.rev_edge(e_max);
                 free[a] = free[b] = false;
               }
             }
                 
             break;
           }

    default: { // refined heuristic
             mwbm_heuristic( G, A, c, pot, free);
             break;
           }
  }

  list<edge> result;

  forall(a,A)
  { if (free[a] && 
      !max_weight_assignment_augment(G,a,c,pot,free,pred,dist,PQ))
    { // graph has no perfect matching
      forall(b,B) 
       forall_out_edges(e,b) G.rev_edge(e);
      list<edge> result;
      return result;  // return empty list
    }
  }

  // TODO

  //forall(b,B) assert(G.outdeg(b) == 1);
  //forall(a,A) assert(G.indeg(a) == 1);

  forall(b,B)
    forall_adj_edges(e,b) result.append(e);

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

  //forall(b,B) assert(G.outdeg(b) == 0);

  return result;
}  



template <class NT>
__temp_func_inline
list<edge> MIN_WEIGHT_ASSIGNMENT_T(graph& G,                         
                        const list<node>& A, const list<node>& B, 
                        const edge_array<NT>& c, node_array<NT>& pot)
{ edge_array<NT> w(G);
  edge e;
  forall_edges(e,G) w[e] = - c[e];

  list<edge> M = MAX_WEIGHT_ASSIGNMENT_T(G,A,B,w,pot);
  node v;
  forall_nodes(v,G) pot[v] = -pot[v];
  return M;
}


template <class NT>
__temp_func_inline
list<edge> MWMCB_MATCHING_T(graph& G, 
                          const list<node>& A, const list<node>& B,
                          const edge_array<NT>& c, node_array<NT>& pot)
{ NT C = 0;
  edge e;
  forall_edges(e,G)
  { if (c[e] > C) C = c[e];
    if (-c[e] > C) C = -c[e];
  }
  int k = leda_max(A.size(),B.size());
  
  C = 1 + 2*k*C;
  edge_array<NT> c_L(G);
  forall_edges(e,G) c_L[e] = c[e] + C;
  list<edge> M = MAX_WEIGHT_BIPARTITE_MATCHING_T(G,A,B,c_L,pot);
#ifdef TEST
  leda_assert(CHECK_MWBM_T(G,c_L,M,pot),"check in MWMCB_MATCHING_T failed");
#endif
  return M;
}
                                    

template <class NT>
__temp_func_inline
list<edge> MWMCB_MATCHING_T(graph& G, 
                          const list<node>& A, 
                          const list<node>& B,
                          const edge_array<NT>& c)
{ node_array<NT> pot(G);
  list<edge> M = MWMCB_MATCHING_T(G, 
                           A,B,c,pot);
  return M;
}

#if LEDA_ROOT_INCL_ID == 450466
#undef LEDA_ROOT_INCL_ID
#include <LEDA/POSTAMBLE.h>
#endif



LEDA_END_NAMESPACE