/*******************************************************************************
+
+  LEDA 3.5
+
+  _mc_matching.c
+
+  This file is part of the LEDA research version (LEDA-R) that can be 
+  used free of charge in academic research and teaching. Any commercial
+  use of this software requires a license which is distributed by the
+  LEDA Software GmbH, Postfach 151101, 66041 Saarbruecken, FRG
+  (fax +49 681 31104).
+
+  Copyright (c) 1991-1997  by  Max-Planck-Institut fuer Informatik
+  Im Stadtwald, 66123 Saarbruecken, Germany     
+  All rights reserved.
+ 
*******************************************************************************/
/*******************************************************************************
*
*  MAX_CARD_MATCHING  (maximum cardinality matching)
*
*  An implementation of Edmonds' algorithm 
*
*
*  J. Edmonds:  Paths, trees, and flowers 
*               Canad. J. Math., Vol. 17, 1965, 449-467
*
*  R.E. Tarjan: Data Structures and Network Algorithms,
*               CBMS-NFS Regional Conference Series in Applied Mathematics,
*               Vol. 44, 1983
*
*
*  running time: O(n*m*alpha(n))
*
*
*  S. Naeher (October 1991)
*
*******************************************************************************/


#include <LEDA/graph_alg.h>
#include <LEDA/node_partition.h>


enum LABEL {ODD, EVEN, UNREACHED};

static int cmp_degree(const node& v,const node& w)  // compare indeg of v and w
{ return indeg(v) - indeg(w); }


static int greedy( graph &G, node_array<node>& mate )   // greedy heuristic
{ node v;
  int count = 0;
  forall_nodes(v,G)
    if (mate[v] == nil)
    { edge e;
      forall_adj_edges(e,v)
      { node w = target(e);
        if (v != w && mate[w] == nil) 
        { mate[v] = w;
          mate[w] = v;
          count++;
          break;
         }
       }
     }
   return count;
 }



static void heuristic( graph &G, node_array<node>& mate )
{
 // Heuristic (Markus Paul):
 // finds almost all augmenting paths of length <=3 with two passes over
 // the adjacency lists
 // ("almost": discovery of a blossom {v,w,x,v} leads to a skip of the
 // edge (x,v), even if the base v stays unmatched - it's not worth while
 // to fix this problem)
 // if all adjacent nodes w of v are matched, try to find an other 
 // partner for mate[x], and match v and w on success

  node u,v;
  bool found;
  node_array<int> all_matched(G,false);

  forall_nodes( v, G )
  {
    edge e = G.first_adj_edge(v);

    if( mate[v]==nil )
    {
        while(e != nil && mate[target(e)]!=nil) e = G.adj_succ(e);

        if( e != nil ) 
          { mate[v] = target(e);  
            mate[target(e)] = v ;  
           }
        else // second pass
          {
	    all_matched[v] = true;

            forall_adj_edges(e,v)
            {
              node w = target(e);
              node x = mate[w];
    
              if( ! all_matched[x] ) 
              {
                while((found=G.next_adj_node(u,x)) && (u==v || mate[u]!=nil));
    
                if( found ) 
                  { mate[u] = x;  mate[x] = u ;
                    mate[v] = w;  mate[w] = v ;
                    G.init_adj_iterator(x);
                    break;
    	           }
	        else
                   all_matched[x] = true;
               }
             }
           }
        } 
    }

}




static void find_path(list<node>& L, node_array<int>&  label,
                                     node_array<node>& pred,
                                     node_array<node>& mate,
                                     node_array<edge>& bridge,
                                     node x, node y)
{
  /* computes an even length alternating path from x to y begining with a 
     matching edge (Tarjan: Data Structures and Network Algorithms, page 122)
     Preconditions: 
     a) x and y are even or shrinked
     b) when x was made part of a blossom for the first time, y was a base 
        and predecessor of the base of that blossom
   */

 if (x==y) 
 { // [ x ]
   L.append(x);
   return;
  }

 if (label[x] == EVEN) 
 { // [ x --> mate[x] --> path(pred[mate[x]],y) ]
   find_path(L,label,pred,mate,bridge,pred[mate[x]],y);
   L.push(mate[x]);
   L.push(x);
   return;
 }

 if (label[x] == ODD) 
 { // [ x --> REV(path(source(bridge),mate[x])) --> path(target(bridge),y)) ]
   node v;
   list<node> L1;
   find_path(L,label,pred,mate,bridge,target(bridge[x]),y);
   find_path(L1,label,pred,mate,bridge,source(bridge[x]),mate[x]);
   forall(v,L1) L.push(v);
   L.push(x);
   return;
  }
}




list<edge> MAX_CARD_MATCHING(graph& G, int heur)
{ 

  // input: directed graph G, we make it undirected (symmetric) by inserting
  //        reversal edges


  list<edge> R;

  int    strue = 0;
  bool   done  = false;

  node_array<int>  label(G,UNREACHED);
  node_array<node> pred(G,nil);
  node_array<int>  path1(G,0);
  node_array<int>  path2(G,0);
  node_array<edge> bridge(G,nil);

  node_array<node> mate(G,nil);


  // make graph symmetric

  Make_Bidirected(G,R);

  G.make_map();


  switch (heur) {

  case 1: { greedy(G,mate);
            break;
          }

  case 2: { G.sort_nodes(cmp_degree);  // sort nodes by degree
            heuristic(G,mate);
            break;
           }
  }


  while (! done)  /* main loop */
  {
    node_list Q;
    node v;

    node_partition base(G);    // now base(v) = v for all nodes v

    done = true;
   
    forall_nodes(v,G)
    { pred[v] = nil;
  
      if (mate[v] == nil)
      { label[v] = EVEN; 
        Q.append(v); 
       }
      else label[v] = UNREACHED;
    }
  
    while (!Q.empty())    // search for augmenting path
    {
      node v = Q.pop();
      edge e;

      forall_adj_edges(e,v)
      {
        node w = G.target(e);

        if (v == w) continue;    // ignore self-loops

        if (base(v) == base(w) || (label[w] == ODD && base(w) == w)) 
        continue;   // do nothing

        if (label[w] == UNREACHED)
        { label[w] = ODD;
          pred[w] = v;
          label[mate[w]] = EVEN;
          Q.append(mate[w]);
         }
        else  // base(v) != base(w) && (label[w] == EVEN || base(w) !=w)
        { 
          node hv = base(v);
          node hw = base(w);
  
          strue++;
          path1[hv] = path2[hw] = strue;
  
          while ((path1[hw] != strue && path2[hv] != strue) &&
                 (mate[hv] != nil || mate[hw] != nil) )
          {
            if (mate[hv] != nil)
            { hv = base(pred[mate[hv]]);
              path1[hv] = strue;
             }
  
            if (mate[hw] != nil)
            { hw = base(pred[mate[hw]]);
              path2[hw] = strue;
             }
           }
  
          if (path1[hw] == strue || path2[hv] == strue)  // Shrink Blossom
          { node x;
            node y;
            node b = (path1[hw] == strue) ? hw : hv;    // Base

#if defined(REPORT_BLOSSOMS)
  cout << "SHRINK BLOSSOM\n";
  cout << "bridge = "; 
  G.print_edge(e);
  newline;
  cout << "base   = "; 
  G.print_node(b);
  newline;
  cout << "path1  = ";
#endif

            x = base(v);
            while (x != b)
            { 

#if defined(REPORT_BLOSSOMS)
  G.print_node(x); 
  cout << " ";
#endif
              base.union_blocks(x,b);
              base.make_rep(b);

              x = mate[x];

#if defined(REPORT_BLOSSOMS)
  G.print_node(x); 
  cout << " ";
#endif
              y = base(pred[x]);
              base.union_blocks(x,b);
              base.make_rep(b);

              Q.append(x);

              bridge[x] = e;
              x = y;
             }

#if defined(REPORT_BLOSSOMS)
  G.print_node(b);
  newline;
  cout << "path2  = ";
#endif

            x = base(w);
            while (x != b)
            { 

#if defined(REPORT_BLOSSOMS)
  G.print_node(x); 
  cout << " ";
#endif
              base.union_blocks(x,b);
              base.make_rep(b);

              x = mate[x];

#if defined(REPORT_BLOSSOMS)
  G.print_node(x); 
  cout << " ";
#endif
              y = base(pred[x]);

              base.union_blocks(x,b);
              base.make_rep(b);

              Q.append(x);

              bridge[x] = G.reversal(e);

              x = y;
             }

#if defined(REPORT_BLOSSOMS)
  G.print_node(b);
  newline;
  newline;
#endif
          }
          else  // augment path
          {
            list<node> P0,P1;

            find_path(P0,label,pred,mate,bridge,v,hv);
            P0.push(w);

            find_path(P1,label,pred,mate,bridge,w,hw);
            P1.pop();

            while(! P0.empty())
             { node a = P0.pop();
               node b = P0.pop();
               mate[a] = b;
               mate[b] = a;
              }

            while(! P1.empty())
             { node a = P1.pop();
               node b = P1.pop();
               mate[a] = b;
               mate[b] = a;
              }

             Q.clear();                /* stop search (while Q not empty) */
             done = false;             /* continue main loop              */
             G.init_adj_iterator(v);
             break;                    /* forall_adj_edges(e,v) ...       */
           }

        } // base(v) != base(w) && (label[w] == EVEN || base(w) !=w)
  
      } // for all adjacent edges
  
    } // while Q not empty

  } // while not done


 // restore graph (only original edges in result!)

 edge e;
 forall(e,R) G.del_edge(e);  


 list<edge> result;

/*
  node v;
  forall_nodes(v,G)
    if (mate[v] != nil)
    { edge e;
      forall_adj_edges(e,v)
        if (mate[target(e)] == v) 
        { result.append(e);
          break;
         }
     }
*/

 forall_edges(e,G) 
 { node v = source(e);
   node w = target(e);
   if ( v != w  &&  mate[v] == w ) 
   { result.append(e);
     mate[v] = v;
     mate[w] = w;
    }
  }

 return result;

}

