mw_matching0.t

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

    pot[u] = leda_max(pot[u], (w[e]/2));
    pot[v] = leda_max(pot[v], (w[e]/2));
  }
  
  int free = G.number_of_nodes();
  mate.init(G, nil);

  forall_edges(e, G) {
    u = source(e);
    v = target(e);
    if ((pot[u] + pot[v] == w[e]) &&
         (mate[u] == nil) && (mate[v] == nil)) {
      mate[v] = u;
      mate[u] = v;
  #ifdef MW_MATCHING_DEBUG
      if (debug) cout << "(" << index(u) << ", " << index(v) << ")" << endl;
  #endif
      free -= 2;
    }    
  }
  
  if (perfect) {
    forall_nodes(u, G) {
      if (!mate[u]) {
        NT slack = MAX_VALUE(NT);
        forall_inout_edges(e, u) {
          v = opposite(u, e);
          slack = leda_min(pot[u] + pot[v] - w[e], slack);
        }
        pot[u] -= slack;
      }
    }
  }
#ifdef _CHECK
  cout << "\t\tCHECK FEASIBILITY..." << std::flush;
  check_feasibility(G, w, pot, mate, perfect);
  cout << "OK!" << endl;
#endif
#ifdef _INFO
  printf("ready! %d free vertices. TIME: %3.2f sec.\n", free, used_time(t));
#endif

  return free;
}




// ----- #include "fractional.t"


extern bool debug;

void alternate_cycle(node u, node_array<node> &mate, 
                     node_array<node> &pred) {

  node cur1 = pred[u];  
  while (cur1 != u) {
    mate[pred[cur1]] = cur1;
    mate[cur1] = pred[cur1];
    node h = pred[cur1];
    pred[cur1] = nil;
    cur1 = h;
    h = pred[cur1];
    pred[cur1] = nil;
    cur1 = h;
  }
  pred[u] = nil; 
}

void alternate_path(node u, node_array<int> &label,
                    node_array<node> &mate, node_array<node> &pred) {

  node cur = u;
  node pre = nil, nxt;
  
  while (cur) {
    if (label[cur] == even) {
      nxt = mate[cur];
      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) << "  " << std::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, bool &is_opt) {
#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;
          }
        }
      }
    }
  }

  is_opt = true;  
  forall_nodes(u, G)
    if (pred[u]) {
      alternate_cycle(u, mate, pred);
      free++; is_opt = false;
    } 
#ifdef _CHECK
  cout << "\t\tCHECK FEASIBILITY..." << std::flush;
  check_feasibility(G, w, pot, mate, perfect);
  cout << "OK!" << endl;
#endif
#ifdef _INFO
  printf("\tready! %d free vertices. is_opt = %d. TIME: %3.2f sec.\n", free, is_opt, used_time(t));

  double W = 0;
  forall_edges(e, G) {
    u = source(e);
    v = target(e);
  
    if (mate[u] && (mate[u] == v) &&
        mate[v] && (mate[v] == u))
    W += w[e];

    // M.push(e);
  }

  printf("\tw(M) = %.1f\n", W);
#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;

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

  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)

// ----- #include "SST.t"


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);