mw_matching.t

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

      mate[cur] = pre;
      cur = nxt;      
    }
    else {
      pre = cur;
      mate[cur] = pred[cur];
      nxt = pred[cur];
      pred[cur] = nil;
      cur = nxt;      
    }  
  }
}

template<class NT>
__temp_func_inline
void destroy_tree(node_slist &T, node_array<int> &label, node_array<NT> &pot, 
                  node_array<node> &mate, node_array<node> &pred, NT Delta) {

  node v;
  while (!T.empty()) {
    v = T.pop();
    if (label[v] == even) pot[v] -= Delta;
    else {
      pot[v] += Delta;
      if (mate[v]) pred[v] = nil;  // only for vertices not on cycle
    }      
    label[v] = unlabeled;
  }
}

void seek_lca(node u, node v, node &lca, 
              node_array<node> &mate, node_array<node> &pred,
              node_array<double> &P1, node_array<double> &P2, 
              double &lock) {

  node cur1 = u, cur2 = v;
  P1[cur1] = P2[cur2] = ++lock;

  while ((P1[cur2] != lock) && (P2[cur1] != lock) && 
         (mate[cur1] || mate[cur2])) {
    
    if (mate[cur1]) {      
      cur1 = pred[mate[cur1]];
      P1[cur1] = lock;
    }
    if (mate[cur2]) {
      cur2 = pred[mate[cur2]];
      P2[cur2] = lock;
    }
  }
  if (P1[cur2] == lock)  // cur2 is lca
    lca = cur2;  
  else if (P1[cur1] == lock)  // cur1 is lca
    lca = cur1;
  else lca = nil;
}

void construct_cycle(node u, node v, node lca, 
                     node_array<node> &mate, node_array<node> &pred) {

  node cur1 = u, cur2 = v; 
  while (cur1 != lca) {
    // set pred data to reversed tree path; delete mate entries
    node h = mate[cur1];
    mate[cur1] = nil;
    cur1 = pred[h];
    pred[h] = mate[h];
    mate[h] = nil;
    pred[cur1] = h;
  }
  while (cur2 != lca) {
    // set pred data to tree path; delete mate entries
    node h = mate[cur2];
    pred[cur2] = mate[cur2];
    mate[cur2] = nil;
    mate[h] = nil;
    cur2 = pred[h];
  }  
  pred[u] = v;
#ifdef _CHECK
  if (debug) cout << "check cycle..." << index(cur1) << "  " << flush; 
  // check correctness of odd cycle
  cur1 = pred[lca];
  int cnt = 1;
  while (cur1 != lca) {
    cur1 = pred[cur1];
    assert(mate[cur1] == nil);
    cnt++;
  }
  assert(cnt % 2);
#endif
}

template<class NT> 
__temp_func_inline
int jump_start(const graph &G, const edge_array<NT> &w, 
               node_array<NT> &pot, node_array<node> &mate, 
               bool perfect) {
#ifdef _INFO
  cout << "\tSOLVING half integral matching problem..." << endl;
  float t = used_time();
#endif
  
  
  edge e; 
  node u, v, r;

  node_array<int>  label(G);
  node_array<node> pred(G, nil);

  NT               delta1;
  NT               delta2a;
  node_pq<NT>      delta2b(G); 

  node             resp_d1  = 0;
  edge             resp_d2a = 0;
  node_array<edge> resp_d2b(G);

  NT               Delta = 0;

  slist<edge>        tight;   
  slist<node>        Q;       
  node_slist         T(G);

  double             lock = 0;
  node_array<double> P1(G, 0);
  node_array<double> P2(G, 0); 
  
  int free = greedy_matching(G, w, pot, mate, perfect);
  if (free == 0) return free;
  forall_nodes(u, G) 
    label[u] = (mate[u] ? unlabeled : even);

  forall_nodes(r, G) {
    if (mate[r] || pred[r]) continue;
#ifdef MW_MATCHING_DEBUG
    if (debug) cout << "Grow search tree from " << index(r) << endl;
#endif

    
    delta1 = delta2a = MAX_VALUE(NT);
    delta2b.clear(); 
    Q.clear(); tight.clear(); 

    pot[r] += Delta;
    T.append(r); Q.append(r);

    bool terminate = false;
    while (!terminate) {

      
      while (!Q.empty()) {
        u = Q.pop();    
        NT pot_u = pot[u] - Delta;

        if (!perfect) {
          
          if (pot_u < leda_min(delta1, delta2a) - Delta) {
            delta1 = pot_u + Delta;
            resp_d1 = u;
            if (delta1 == Delta) break;
          }
        }
        forall_inout_edges(e, u) {
          v = opposite(u, e);
          if (label[v] == odd) continue;
          NT pot_v = pot[v] - (((label[v] == even) && T.member(v)) ? Delta : 0);
          NT pi = pot_u + pot_v - w[e];
                  
          if (pi == 0) { 
            
            if ((label[v] == unlabeled) && mate[v]) tight.append(e);
            else {
              tight.clear(); Q.clear();
              tight.append(e);
              break;
            }
          }
          else {  
            
            #if !defined(_NO_PRUNING)
            if (label[v] == even && T.member(v)) {
              if (pi/2 + Delta >= leda_min(delta1, delta2a)) continue;
            }
            else if (pi + Delta >= leda_min(delta1, delta2a)) continue;
            #endif
            if ((label[v] == unlabeled) && mate[v]) {
              
              if (delta2b.member(v)) {
                if (pi < delta2b.prio(v) - Delta) {
                  delta2b.decrease_p(v, pi + Delta);
                  resp_d2b[v] = e;
                }
              }
              else {
                delta2b.insert(v, pi + Delta);
                resp_d2b[v] = e;
              }
            }
            else {
              
              if ((label[v] == even) && T.member(v)) pi /= 2;

              if (pi < delta2a - Delta) {
                delta2a = pi + Delta;
                resp_d2a = e;
              }
            }
          }   
        }    
      }

      if (delta1 == Delta) { 
        
        alternate_path(resp_d1, label, mate, pred);
        destroy_tree(T, label, pot, mate, pred, Delta);
        label[resp_d1] = even;
        terminate = true; 
      }
      else if (!tight.empty()) { 
        
        while (!tight.empty()) {
          e = tight.pop();
          u = (T.member(source(e)) ? source(e) : target(e));
          v = opposite(u, e);

          if (label[v] == odd || label[u] == odd) continue;

          if (label[v] == unlabeled) {
            if (mate[v]) {  
              
              label[v] = odd;
              pred[v] = u;
              pot[v] -= Delta;
              T.append(v);
              delta2b.del(v);
              resp_d2b[v] = nil;

              node m = mate[v]; 
              label[m] = even;
              pot[m] += Delta;
              T.append(m);
              Q.append(m);
              delta2b.del(m);
              resp_d2b[m] = nil;
            }
            else {  // v on half-valued cycle
              
              alternate_cycle(v, mate, pred);
              alternate_path(u, label, mate, pred);
              mate[u] = v;
              mate[v] = u;
              destroy_tree(T, label, pot, mate, pred, Delta);
              free--;
              terminate = true;
              break;
            }
          }
          else {  // label[v] == even
            if (T.member(v)) { 
              
              node lca;
              seek_lca(u, v, lca, mate, pred, P1, P2, lock);
              alternate_path(lca, label, mate, pred);
              construct_cycle(u, v, lca, mate, pred);
              destroy_tree(T, label, pot, mate, pred, Delta);
              free--;
              terminate = true; 
              break;
            }
            else { 
              
              alternate_path(u, label, mate, pred);
              mate[u] = v;
              mate[v] = u;
              label[v] = unlabeled;
              destroy_tree(T, label, pot, mate, pred, Delta);
              free -= 2;
              terminate = true; 
              break;
            }
          }           
        } 
      }
      else {
        
        NT cand2b = (delta2b.empty() ? \
                     MAX_VALUE(NT) : delta2b.prio(delta2b.find_min()));

        NT delta = leda_min(delta1, leda_min(delta2a, cand2b));

        if (delta == MAX_VALUE(NT) && perfect) {
          mate.init(G, nil);
          return 0; 
        }
        #ifdef MW_MATCHING_DEBUG  
        if (debug) {
          cout << "\tPERFORMED DUAL ADJUSTMENT by " << delta-Delta << endl;
          cout << "\tVALUE OF DELTA: " << delta << endl;
        }
        #endif
        Delta = delta;  // corresponds to Delta += (delta - Delta)
        
        if (delta1 == Delta)
          continue;
        else if (delta2a == Delta) {
          tight.append(resp_d2a);
          resp_d2a = nil;
        }
        else {
          while (!delta2b.empty() && 
                 (delta2b.prio(delta2b.find_min()) == Delta)) {
            u = delta2b.del_min();
            tight.append(resp_d2b[u]);
            resp_d2b[u] = nil;
          }
        }
      }
    }
  }
  
  forall_nodes(u, G)
    if (pred[u]) {
      alternate_cycle(u, mate, pred);
      free++;
    } 
#ifdef _CHECK
  cout << "\t\tCHECK FEASIBILITY..." << flush;
  check_feasibility(G, w, pot, mate, perfect);
  cout << "OK!" << endl;
#endif
#ifdef _INFO
  printf("\tready! %d free vertices. TIME: %3.2f sec.\n", free, used_time(t));
#endif

  return free;
}

template<class NT> 
__temp_func_inline
list<edge> fractional_matching(const graph &G, 
                               const edge_array<NT> &w, 
                               node_array<NT> &pot,
                               bool perfect) {

  node_array<node> mate;

  jump_start(G, w, pot, mate, perfect);

  list<edge> M;

  node_array<int> b(G, -1);
  array<two_tuple<NT, int> > BT;
  int k = 0;

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

#ifdef _CHECK
  if (!perfect) 
    CHECK_MAX_WEIGHT_MATCHING(G, w, M, pot, BT, b);
  else     
    CHECK_MAX_WEIGHT_PERFECT_MATCHING(G, w, M, pot, BT, b);
#endif 

  return M;
}




#if defined(_SST_APPROACH)



#ifdef _STAT
int   adj_c, 
      alternate_c, grow_c, shrink_c, expand_c, augment_c, destroy_c,
      scan_c;             
float alternate_t, grow_t, shrink_t, expand_t, augment_t, destroy_t, 
      scan_t;
float heur_t, init_t, alg_t, extract_t; 
#endif


template<class NT> class vertex; 
template<class NT> class blossom;


template<class NT> 
inline blossom<NT>* blossom_of(c_pq_item it, NT) {
  return (it ? (blossom<NT>*)(get_owner(it)) : nil);
}

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;