mw_matching0.t

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

}

template<class NT> 
inline c_pq_item new_blossom(NT d, node b, blossom<NT>* &B) {
  B = new blossom<NT>(b);
  vertex<NT> *v = new vertex<NT>(d, b);
  return B->init(MAX_VALUE(NT), v);
} 

template<class NT> 
class blossom : public virtual concat_pq<NT, vertex<NT>*> {
  
public:
  
  LABEL label;
  NT    pot;
  NT    offset;

  node base, mate;
  node disc, pred;

  list<node>         shrink_path;
  list<blossom<NT>*> subblossom_p; 

  c_pq_item          split_item;

  pq_item item_in_pq;

  list_item item_in_T;
  list_item item_in_O;

  double marker1, marker2;
  
  const NT min_prio() const 
  { return (find_min() ? prio(find_min()) : MAX_VALUE(NT)); }

  NT   pot_of  (c_pq_item it) const { return inf(it)->pot;      } 
  node node_of (c_pq_item it) const { return inf(it)->my_node;  }
  node best_adj(c_pq_item it) const { return inf(it)->best_adj; }

  bool trivial() const { return size() == 1; }

  void print_vertices() const {
    vertex<NT> *Cur;
    forall(Cur, *this) cout << index(Cur->my_node) << "  ";
  }

  blossom(node ba = nil) : concat_pq<NT, vertex<NT>*>() {

    set_owner(leda_cast(this));
    label = even;
    pot = offset = 0; 
    base = ba; mate = nil;
    disc = pred = nil;
    marker1 = marker2 = 0;
    item_in_T  = item_in_O = nil; 
    item_in_pq = nil;
    split_item = nil;
  }

  ~blossom() { clear(); }

  void destroy() {
    c_pq_item it;
    forall_items(it, *this)
      delete this->inf(it);
    this->clear();
    delete this;
  }

  void status_change(LABEL l, NT Delta, list<blossom<NT>*> &T, node_slist &Q) {

    if (l == unlabeled) {
      
      assert((label != l) && item_in_T);
      offset += (label == odd ? Delta : -Delta);
      T.del(item_in_T);
      item_in_T = nil;
    }
    else if (l == odd) {
      
      assert((label != l) && !item_in_T);
      offset -= Delta;
      item_in_T = T.append(this);    
    }
    else if (l == even) {
      
      assert((label != l) || !item_in_T);
      if (label == odd) offset += 2*Delta;
      else {  // non-tree blossom
        offset += Delta;
        item_in_T = T.append(this);
      }  
      
      c_pq_item it;
      forall_items(it, *this) {
        inf(it)->pot += offset;
        Q.append(node_of(it));
      }
      if (!trivial()) pot -= 2*offset;
      offset = 0;
    }
    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;
  }

  // GSchange
  // 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 = clean_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);

  bool is_opt;
 
  #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, is_opt); break; }      
  }
  #ifdef _STAT
  heur_t = used_time(heur_t);
  #endif

  // GSchange
#ifdef _FORMER
  if (free == 0) {  
#else
  if (is_opt) {
#endif
    
    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];