mw_matching.t

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

    }
    label = l;
  }

  bool improve_connection(c_pq_item it, NT x, node u) {

    if (!it) return false;
    NT old_min = min_prio();
    if (decrease_p(it, x)) inf(it)->best_adj = u;
    return old_min != min_prio();    
  }

  bool delete_connection(c_pq_item it) {

    if (!it) return false;
    NT old_min = min_prio();
    del_item(it);
    inf(it)->best_adj = nil;
    return old_min != min_prio();
  }

  void append_subblossom(blossom<NT> *CUR, NT Delta,
                         list<blossom<NT>*> &T, node_slist &Q) {

    if (CUR->label == odd) 
      CUR->status_change(even, Delta, T, Q);

    if (!CUR->trivial())
      CUR->pot += -2*CUR->offset + 2*Delta;

    T.del(CUR->item_in_T);
    CUR->item_in_T = nil;

    concat(*CUR);
    CUR->split_item = last_item();
    subblossom_p.append(CUR);   
  }

  void expand(NT Delta) {
      
    blossom<NT> *CUR;
    forall(CUR, subblossom_p) {
      split_at_item(CUR->split_item, *CUR, *this);
      CUR->offset = offset;
      CUR->label = label;

      if (!CUR->trivial() && label == odd)
        CUR->pot += 2*offset + 2*Delta;
      else if (!CUR->trivial()) {
        assert(!CUR->item_in_T);
        assert(CUR->label == even || CUR->label == unlabeled);
        CUR->pot += 2*offset;
      }
    }
  }

  int restore_matching(blossom<NT> *BASE, blossom<NT> *DISC) {

    
    while (subblossom_p.tail() != BASE) {
      subblossom_p.append(subblossom_p.pop());
      shrink_path.append(shrink_path.pop());
      shrink_path.append(shrink_path.pop());
    }
    
    BASE->mate = mate;
    BASE->base = base;

    node ba, m;
    int dist = 0, pos = 1;
    list_item p_it   = shrink_path.first();
    list_item sub_it = subblossom_p.first();
    blossom<NT> *CUR = subblossom_p.inf(sub_it), *ADJ;

    while (CUR != BASE) {
      if (CUR == DISC) dist = pos;
      sub_it = subblossom_p.succ(sub_it); pos++;
      ADJ    = subblossom_p.inf(sub_it);
      if (ADJ == DISC) dist = pos; 

      p_it = shrink_path.succ(p_it); 
      p_it = shrink_path.succ(p_it); ba = shrink_path.inf(p_it);
      p_it = shrink_path.succ(p_it); m = shrink_path.inf(p_it);

      CUR->base = ba; CUR->mate = m;
      ADJ->base = m; ADJ->mate = ba;
         
      sub_it = subblossom_p.succ(sub_it); pos++;
      CUR    = subblossom_p.inf(sub_it);
      p_it   = shrink_path.succ(p_it);
    }
    if (CUR == DISC) dist = pos; 
    return dist;
  }

  LEDA_MEMORY(blossom<NT>);
};

template<class NT> class vertex {

public:

  NT   pot;  
  node my_node;  
  node best_adj;

  vertex(NT d, node u) {
    pot = d;
    my_node = u;
    best_adj = nil;
  }

  LEDA_MEMORY(vertex<NT>);
};

#define _BLOSSOM_OF(this_node) \
(this_node ? blossom_of(item_of[this_node],NT(0)) : nil)

template<class NT>
inline NT compute_potential(blossom<NT> *CUR, NT Delta, c_pq_item it) {

  int a = (it == nil ? -2 : 1);
  int sigma = 0;
  if (CUR->item_in_T) sigma = (CUR->label == even ? -1 : 1);
  NT stored = (it == nil ? CUR->pot : CUR->pot_of(it));

  return stored + a * CUR->offset + a * sigma * Delta;
}

// sn: still some compilers do not support default arguments for 
//     template functions

template<class NT>
inline NT compute_potential(blossom<NT> *CUR, NT Delta) {
    return compute_potential(CUR,Delta, (c_pq_item)nil);
}



template<class NT> 
__temp_func_inline
void alternate_path(blossom<NT>* RESP, node_array<c_pq_item> &item_of) {

  if (!RESP) return;
  blossom<NT> *CUR = RESP;
  node pred = RESP->base, disc = nil, mate;

  while (CUR) {
    if (CUR->label == even) {
      mate = CUR->mate;
      CUR->mate = disc;
      CUR->base = pred;
      CUR = _BLOSSOM_OF(mate);
    }
    else {  // CUR->label == odd
      pred = CUR->pred;
      disc = CUR->disc;
      CUR->mate = pred;
      CUR->base = disc;
      CUR = _BLOSSOM_OF(pred);    
    }
  }
}

template <class NT> 
__temp_func_inline
int restore_matching(blossom<NT>* RESP, blossom<NT>* BASE, blossom<NT>* DISC,
                     list<blossom<NT>*> &subblossom_p, list<node> &shrink_path) {

  BASE->mate = RESP->mate;
  BASE->base = RESP->base;

  while (subblossom_p.tail() != BASE) {
    subblossom_p.append(subblossom_p.pop());
    shrink_path.append(shrink_path.pop());
    shrink_path.append(shrink_path.pop());
  }

  int dist, pos = 1;
  node base, mate;

  list_item p_it   = shrink_path.first();
  list_item sub_it = subblossom_p.first();
  blossom<NT>* CUR = subblossom_p.inf(sub_it), *ADJ;

  while (CUR != BASE) {

    if (CUR == DISC) dist = pos;

    pos++;
    sub_it = subblossom_p.succ(sub_it);
    ADJ    = subblossom_p.inf(sub_it);

    if (ADJ == DISC) dist = pos; 
  
    p_it = shrink_path.succ(p_it);
    p_it = shrink_path.succ(p_it); base = shrink_path.inf(p_it);
    p_it = shrink_path.succ(p_it); mate = shrink_path.inf(p_it);

    CUR->base = base; CUR->mate = mate;
    ADJ->base = mate; ADJ->mate = base;
    
    pos++;
    sub_it = subblossom_p.succ(sub_it);
    CUR    = subblossom_p.inf(sub_it);
    p_it   = shrink_path.succ(p_it);
  }
  if (CUR == DISC) dist = pos;
  return dist;
}

template<class NT> 
__temp_func_inline
void unpack_blossom(blossom<NT> *RESP, const node_array<c_pq_item> &item_of,
                    node_array<NT> &pot, node_array<int> &b, 
                    array<two_tuple<NT, int> > &BT,
                    int &k, int parent, NT Delta) {

  if (RESP->trivial()) {
    
    node cur = RESP->node_of(RESP->first_item());
    if (b[cur] != -1) return;
    pot[cur] = compute_potential(RESP, Delta, item_of[cur]);
    b[cur] = parent;      
    #ifdef MW_MATCHING_DEBUG
      if (debug) cout << "trivial blossom " << index(cur) << ": dv:" << pot[cur] 
                      << " immediate super blossom: " << b[cur] << endl;
    #endif
  }
  else {
    
    if (k > BT.high()) BT.resize(2*k+1);
    BT[k].first()  = compute_potential(RESP, Delta);
    BT[k].second() = parent;
    int idx = k++;
    #ifdef MW_MATCHING_DEBUG
      if (debug) {
        cout << "non-trivial blossom "; RESP->print_vertices();
        cout << ": dv: " << BT[idx].first() << ", index: " << idx << " and parent blossom " 
             << BT[idx].second() << endl;
      }
    #endif

    RESP->expand(Delta);
    blossom<NT> *BASE = _BLOSSOM_OF(RESP->base);
    blossom<NT> *DISC = nil;
    RESP->restore_matching(BASE, DISC);

    blossom<NT>* CUR;
    forall(CUR, RESP->subblossom_p) 
      unpack_blossom(CUR, item_of, pot, b, BT, k, idx, Delta);

    // mod. 12/09/00
    RESP->destroy();
  }
}

template<class NT>
__temp_func_inline
void print_pqs(const NT delta1, node resp_d1,
               const NT delta2a, edge resp_d2a, 
               const p_queue<NT, blossom<NT>*> &delta2b,
               const p_queue<NT, edge> &delta3,
               const p_queue<NT, blossom<NT>*> &delta4,
               const NT Delta, bool perfect) {

  pq_item it;

  if (!perfect) 
    cout << "delta1: " << delta1-Delta << " (" << index(resp_d1) << ")" << endl;
  cout << "delta2a: ";
  if (delta2a != MAX_VALUE(NT)) 
    cout << delta2a-Delta << " (" << index(source(resp_d2a)) << ", " 
         << index(target(resp_d2a)) << ")" << endl;
  else cout << delta2a << endl;

  cout << "delta2b: " << endl;
  forall_items(it, delta2b) {
    cout << delta2b.prio(it) - Delta << "\t"; 
    delta2b.inf(it)->print_vertices(); cout << endl;
  }

 cout << "delta3: " << endl;
  forall_items(it, delta3)
    cout << delta3.prio(it) - Delta << "\t"
         << "(" << index(source(delta3.inf(it))) << ", " 
         << index(target(delta3.inf(it))) << ")" << endl;

  cout << "delta4: " << endl;
  forall_items(it, delta4) {
    cout << delta4.prio(it) - Delta << "\t";
    delta4.inf(it)->print_vertices(); cout << endl;
  }
}

template<class NT> 
__temp_func_inline
list<edge> MWM_SST(const graph &G, const edge_array<NT> &w, 
                   node_array<NT> &pot, array<two_tuple<NT, int> > &BT, 
                   node_array<int> &b, int heur/* =1*/, bool perfect/* =false*/)
{
#ifdef _STAT
  adj_c = alternate_c = grow_c = shrink_c = expand_c = augment_c = scan_c = destroy_c = 0;
  heur_t = init_t = alg_t = extract_t = alternate_t = grow_t = shrink_t = \
expand_t = augment_t = scan_t = destroy_t = 0;
  init_t = used_time();
#endif

  
  NT                        delta1; 
  NT                        delta2a; 
  p_queue<NT, blossom<NT>*> delta2b;
  p_queue<NT, edge>         delta3;         
  p_queue<NT, blossom<NT>*> delta4; 

  node resp_d1;
  edge resp_d2a;

  NT Delta = 0;

  node_slist Q(G);

  list<blossom<NT>*> T;
  list<blossom<NT>*> O;

  node_array<c_pq_item> item_of(G);

  edge        e;
  node        resp, opst, cur, adj, u, v, r;
  blossom<NT> *RESP, *OPST, *CUR, *ADJ, *R;

  list<edge>  M;

  const NT INFTY = MAX_VALUE(NT);

  double lock = 0;
  
  int free = G.number_of_nodes();
  node_array<node> mate(G, nil);

  #ifdef _STAT
  heur_t = used_time();
  #endif
  switch(heur) {
    case 0:  { 
               forall_nodes(u, G) {
                 if (degree(u) == 0) { pot[u] = 0; continue; }
                 NT max = -INFTY;
                 forall_inout_edges(e, u) max = leda_max(w[e], max);
                 pot[u] = max/2;
               } break; }
    case 1:  { 
               free = greedy_matching(G, w, pot, mate, perfect); break; }
    default: { 
               free = jump_start(G, w, pot, mate, perfect); break; }      
  }
  #ifdef _STAT
  heur_t = used_time(heur_t);
  #endif

  if (free == 0) {
    
    forall_edges(e, G) {
      u = source(e);
      v = target(e);
      
      if (mate[u] && (mate[u] == v) &&
          mate[v] && (mate[v] == u))
        M.push(e);
    }
    return M;
  }

  forall_nodes(u, G) {
    item_of[u] = new_blossom(pot[u], u, CUR);

    if (mate[u]) {
      CUR->mate = mate[u];    
      CUR->label = unlabeled;
    }
  }
#ifdef _STAT
  init_t = used_time(init_t);
  alg_t = used_time();
#endif

  forall_nodes(r, G) {
#ifdef MW_MATCHING_DEBUG
    if (debug) cout << "GROW SEARCH TREE rooted at " << index(r) << endl;
#endif
    R = _BLOSSOM_OF(r);
    if (R->mate) continue;

    
    delta1 = delta2a = INFTY;
    delta2b.clear(); delta3.clear(); delta4.clear();
    Q.clear();

    R->status_change(even, Delta, T, Q);

    bool terminate = false;
    while (!terminate) {

      
      #ifdef _STAT
      scan_c++;
      float scan_h = used_time();
      #endif

      NT cur_pot, adj_pot, actual_p, stored_p;

      while (!Q.empty()) {
        cur     = Q.pop();
        CUR     = _BLOSSOM_OF(cur);
        cur_pot = compute_potential(CUR, Delta, item_of[cur]);

        if (!perfect) {
          
          if (cur_pot < delta1 - Delta) {
            delta1 = cur_pot + Delta;
            resp_d1 = cur;
            // if (delta1 == Delta) break; 
          }
        }
        forall_inout_edges(e, cur) {
          adj = opposite(cur, e);     
          ADJ = _BLOSSOM_OF(adj);

          
          // do not consider edges within a blossom
          if (CUR == ADJ) continue;

          // do not consider tree edges
          if ((ADJ->label == odd) &&
              ((ADJ->base == adj && ADJ->mate == cur) ||
               (ADJ->disc == adj && ADJ->pred == cur))) continue;

          adj_pot  = compute_potential(ADJ, Delta, item_of[adj]);
          actual_p = cur_pot + adj_pot - w[e];

      #ifdef MW_MATCHING_DEBUG
        if (debug) {
          cout << "\tSCANNING edge (" << index(cur) << ", " << index(adj) << "): " << 
            cur_pot << " + " << adj_pot << " - " << w[e] << " = " << actual_p << endl;
        }
      #endif
          
          #if !defined(_NO_PRUNING)
          if ((ADJ->label == even) && ADJ->item_in_T) {
            if (actual_p/2 + Delta > leda_min(delta1, delta2a)) continue;
          }
          else if (actual_p + Delta > leda_min(delta1, delta2a)) continue;
          #endif

          if ((ADJ->label == even) && !ADJ->item_in_T) {
            
            if (actual_p < delta2a - Delta) {
              delta2a  = actual_p + Delta;
              resp_d2a = e;
              if (delta2a == Delta) break;    
            } 
          }
          else if (ADJ->label == unlabeled) {  
            
            stored_p = actual_p - ADJ->offset + Delta;
            if (ADJ->improve_connection(item_of[adj], stored_p, cur))
              if (ADJ->item_in_pq)
                delta2b.decrease_p(ADJ->item_in_pq, actual_p + Delta);
              else {
                ADJ->item_in_pq = delta2b.insert(actual_p + Delta, ADJ);
                ADJ->item_in_O = O.append(ADJ);
              } 
          }
          else if (ADJ->label == even)  // ADJ is even tree blossom
            delta3.insert(actual_p/2 + Delta, e);
          else if (ADJ->label == odd) {  
            stored_p = actual_p - ADJ->offset;
            ADJ->improve_connection(item_of[adj], stored_p, cur);
          }
        }
      }

      #ifdef _STAT
      scan_t += used_time(scan_h);
      #endif         
#ifdef MW_MATCHING_DEBUG
    if (debug) {
      map<blossom<NT>*, bool> printed(false);
      blossom<NT>* B;
      forall_nodes(v, G) {
        B = _BLOSSOM_OF(v);
        LABEL label = B->label;
        NT dv;

        dv = compute_potential(B, Delta,item_of[v]);