matching.cc

#ifndef SCIL_MATCHING_CC
#define SCIL_MATCHING_CC

#define LEPS 1e-8 //must be larger than 1e-10, because of the cut tree alg.

/*namespace {
extern "C" {
   void compute_tree (int n, int m, int* os, int* ot, double* oc, int* cts, int* ctt, double* ctc);
   }
}*/

using namespace boost;
using namespace SCIL;
using namespace std;

template <typename Graph>
MATCHING<Graph>::MATCHING(Graph& G_, map<vertex_descriptor, int>& NC_,
                          var_map<edge_descriptor>& VM_, bool perfect_)
   : G(G_), NC(NC_), VM(VM_), perfect(perfect_) {
   feps = 1.0e-4;
   BGL_FORALL_EDGES_T(e, G, Graph)
      if(VM[e].type() != Vartype_Binary) {
         cerr << "matching.cc: All variables in VM need to be";
         cerr << "Vartype_Binary" << endl;
         exit(1);
      }
   os = new int[num_edges(G) + num_vertices(G)];
   ot = new int[num_edges(G) + num_vertices(G)];
   oc =  new double [num_edges(G) + num_vertices(G)];
   cts = new int[num_vertices(G)];
   ctt = new int[num_vertices(G)];
   ctc = new double[num_vertices(G)];
   edg = new edge_descriptor[num_edges(G)];
}

template <typename Graph>
MATCHING<Graph>::~MATCHING() {
   delete[]os;
   delete[]ot;
   delete[]oc;
   delete[]cts;
   delete[]ctt;
   delete[]ctc;
   delete[]edg;
}

template <typename Graph>
void MATCHING<Graph>::init(subproblem& S) {
   if(S.configuration("MATCHING_Debug_Out")=="true")
      cerr << "MATCHING::init\n";

   typedef map<vertex_descriptor, int> cn_t;
   cn_t cn;
   associative_property_map<cn_t> pm_cn(cn);
   int comps = connected_components(G, pm_cn);
   if( comps != 1 ) {
      if(S.configuration("MATCHING_Debug_Out")=="true")
         cerr << "G is not connected! " << comps << endl;
      //so what?
      //exit(1);
   }

   //Check whether Graph is undirected
   typedef typename GraphTraits::directed_category directed_category;
   if (!(is_same<directed_category,undirected_tag>::value
         || is_base_and_derived<undirected_tag, directed_category>::value)) {
      cerr << "kconnected.cc: The graph has to be undirected." << endl;
      exit(1);
   };

   //Graph has to provides edges() and vertices()
   BOOST_CONCEPT_ASSERT((VertexAndEdgeListGraphConcept<Graph>));
   //Graph has to provide out_edges()
   BOOST_CONCEPT_ASSERT((IncidenceGraphConcept<Graph>));
   //Graph has to provide remove_edge()
   BOOST_CONCEPT_ASSERT((MutableGraphConcept<Graph>));

   int cap = 0;
   int i = 1;
   BGL_FORALL_VERTICES_T(v, G, Graph) {
      nid[v] = i++;
      cap += NC[v];

      if( NC[v] < 0 ) {
         if(S.configuration("MATCHING_Debug_Out")=="true")
            cerr << "MATCHING::init   Node capacities must be positive!" << endl;
         exit(1);
      }

      //add the capacity equalities
      row r;
      typename GraphTraits::out_edge_iterator oe_it, oe_it_end;
      for(tie(oe_it, oe_it_end) = out_edges(v, G); oe_it != oe_it_end; oe_it++)
         r += VM[*oe_it];
      if( perfect )
         S.add_basic_constraint(r==NC[v], Static);
      else
         S.add_basic_constraint(r<=NC[v], Static);
   }

   if( perfect && cap % 2 != 0 ) {
      if(S.configuration("MATCHING_Debug_Out")=="true")
         cerr << "MATCHING::init   Sum of node capacities must be even for perfect matchings!" << endl;
      exit(1);
   }

   i = 0;
   BGL_FORALL_EDGES_T(e, G, Graph) {
      os[i] = nid[source(e, G)];
      ot[i] = nid[target(e, G)];
      edg[i] = e;
      i++;
   }

   return;
}

template <typename Graph>
int MATCHING<Graph>::exactBlossomSeparation(subproblem& S) {
   if(S.configuration("MATCHING_Debug_Out")=="true")
      cerr << "MATCHING::exactBlossomSeparation\n";

   int numfound = 0;
   Graph H;
   map<edge_descriptor, double> ctc;
   map<vertex_descriptor, vertex_descriptor> p;
   map<vertex_descriptor, double> fl;
   double x;
   BGL_FORALL_EDGES_T(b, G, Graph) {
      x = S.value(VM[b]);
      ctc[b] = (x < .5) ? x: 1. - x;
   }
   computeCutTree(G, ctc, p, fl);

   //for each cut induced by an edge in the cut tree determine the cut edges in G 
   BGL_FORALL_VERTICES_T(v, G, Graph) {
      if( p[v] == v )
         continue;
      row s;
      int rt = 0;

      //temporarily delete one edge of the cut tree
      vertex_descriptor ov = p[v];
      p[v] = v;

      //determine the cut nodes
      typedef map<vertex_descriptor, set<vertex_descriptor>* > vertex_partition;
      vertex_partition vpartition;
      //Initialize the vertex_partition
      BGL_FORALL_VERTICES_T(w, G, Graph) {
         vpartition[w] = new set<vertex_descriptor>();
         vpartition[w]->insert(w);
      }

      //Union sets of vertices w with sets of vertices p[w]
      BGL_FORALL_VERTICES_T(w, G, Graph) {
         set<vertex_descriptor> *tmp = vpartition[p[w]];
         vpartition[w]->insert(tmp->begin(), tmp->end());
         typename set<vertex_descriptor>::iterator it;
         //Update the pointers of the vertices in vpartition[p[w]]
         for(it = tmp->begin(); it != tmp->end(); it++)
            vpartition[*it] = vpartition[w];
      }

      vertex_descriptor t = *(vertices(G).first);

      int bW = 0;
      BGL_FORALL_VERTICES_T(w, G, Graph) {
               if( *(vpartition[w]) == *(vpartition[t]) )
                  bW += NC[v]%2; 
               //bW += NC[w]; 
      }

      list<edge_descriptor> F; //the cut edges with 1 - x^* < x^*
      list<edge_descriptor> delta; //the cut edges that are not in F
      BGL_FORALL_EDGES_T(f, G, Graph) {
         if( !(*(vpartition[source(f, G)]) == *(vpartition[target(f, G)])) ) {
            if( S.value(VM[f]) > .5 )
               F.push_back(f);
            else
               delta.push_back(f);
         }
      }

      // if b(W) + |F| is odd, we take F and construct the blossom 
      // inequality, if lhs < 1.
      // Otherwise we construct F' (symmetric difference) and proceed. 

      edge_descriptor f;
      if( ((bW + F.size()) % 2) == 0 ) {
         edge_descriptor ff;
         double m;
         double minm = 2.;
         bool fffound = false;
         BOOST_FOREACH(f, delta) {
            m = fabs((2. * S.value(VM[f])) - 1.);
            if( m <= minm ) {
               minm = m;
               ff = f;
            }
         }
         BOOST_FOREACH(f, F) {
            m = fabs((2. * S.value(VM[f])) - 1.);
            if( m <= minm ) {
               minm = m;
               ff = f;
               fffound = true;
            }
         }
         // the symmetric diff. of F and {ff} is 
         // F \ ff, if ff \in F, F \cup ff else 
         if( fffound ) {
            F.remove(ff);
            delta.push_back(ff);
         }
         else {
            F.push_back(ff);
            delta.remove(ff);
         }
      }

      //construct the violated blossom inequality
      double rhs = 1. - double(F.size()); 
      row r;
      BOOST_FOREACH(f, delta) {
         r += VM[f];
      }
      BOOST_FOREACH(f, F) {
         r -= VM[f];
      }
     
      if( r.size() > 0 ) {
         if(S.configuration("MATCHING_Debug_Out")=="true")
            cerr << "constraint found" << endl;
         if( !perfect ) {
            rhs += double(rt);
            r += s;
         }

         //r.normalize(true);
         cons_obj* c = (r>=rhs);
         if( c->violation(S) > feps ) {
            if(S.configuration("MATCHING_Debug_Out")=="true")
               cerr << r << " >= " << rhs << endl;
            S.add_basic_constraint(r>=rhs); 
            numfound++;
         }
         delete c;
      }
      p[v] = ov;
   }
   if(S.configuration("MATCHING_Debug_Out")=="true")
      cerr << "MAT: " << numfound << endl;
   return numfound;
}

template <typename Graph>
typename MATCHING<Graph>::status MATCHING<Graph>::feasible(solution& S) {
   if(S.originating_subproblem()->configuration("MATCHING_Debug_Out")=="true")
      cerr << "MATCHING::feasible\n";
/*   double b;

   BGL_FORALL_VERTICES_T(v, G, Graph) {
      b = 0.;
      typename GraphTraits::out_edge_iterator oe_it, oe_it_end;
      for(tie(oe_it, oe_it_end) = out_edges(v, G); oe_it != oe_it_end; oe_it++) 
         b += S.value(VM[*oe_it]);
      if( perfect ) {
         if( fabs(b - double(NC[v])) > LEPS )
            return infeasible_solution;
      } else {
         if ( b + LEPS > double(NC[v]) )
            return infeasible_solution;
      }
   }*/
   return feasible_solution;
}

//compare Cunha and Lucena 'A Hybrid Relax and Cut
//Branch-And-Cut Algorithm for the DCMST-Problem
template <typename Graph>
int MATCHING<Graph>::heuristicBlossomSeparation(subproblem& S) {
   if(S.configuration("MATCHING_Debug_Out")=="true")
      cerr << "MATCHING::heuristicBlossomSeparation\n";

   int numfound = 0;

   list<edge_descriptor> edgesH;
   list<edge_descriptor> edgesNotH;
   map<vertex_descriptor, vertex_descriptor> verticesHtoG;
   map<vertex_descriptor, vertex_descriptor> verticesGtoH;
   map<edge_descriptor, edge_descriptor> edgesHtoG;
   Graph H;
   BGL_FORALL_VERTICES_T(v, G, Graph) {
      vertex_descriptor newNodeInH = add_vertex(H);
      verticesHtoG[newNodeInH] = v;
      verticesGtoH[v] = newNodeInH;
   }
   BGL_FORALL_EDGES_T(e, G, Graph) {
      edge_descriptor newEdgeInH = (add_edge(verticesGtoH[source(e, G)], 
                                             verticesGtoH[target(e, G)],
                                             H)).first;
      edgesHtoG[newEdgeInH] = e;
   }

   list<edge_descriptor> toBeDeleted;
   BGL_FORALL_EDGES_T(e, H, Graph) {
      if( (S.value(VM[edgesHtoG[e]]) < LEPS)
          || (S.value(VM[edgesHtoG[e]]) > 1. - LEPS) ) {
         edgesNotH.push_back(edgesHtoG[e]);
         toBeDeleted.push_back(e);
      } else
         edgesH.push_back(edgesHtoG[e]);
   }
   edge_descriptor e;
   BOOST_FOREACH(e, toBeDeleted)
      remove_edge(e, H);

   typedef map<vertex_descriptor, int> compnumH_t;
   compnumH_t compnumH;
   associative_property_map<compnumH_t> pm_compnumH(compnumH);
   int comps = connected_components(H, pm_compnumH);

   if( comps < 2 )
      return 0;

   compnumH_t compnum;
   BGL_FORALL_VERTICES_T(v, H, Graph)
      compnum[verticesHtoG[v]] = compnumH[v];
   
   for( int i = 0; i < comps; i++ ) {
      row r;
      double rhs = 1.;
      edge_descriptor e;
      BOOST_FOREACH(e, edgesNotH) {
               vertex_descriptor s = source(e, G);
               vertex_descriptor t = target(e, G);
               if( ((compnum[s] == i) || (compnum[t] == i)) && compnum[s] != compnum[t] ) {
                  if( S.value(VM[e]) > 1. - LEPS ) {
                     r -= VM[e];
                     rhs -= 1.;
                  }
                  else 
                     r += VM[e];
               }
      }

      int b = 0;
      int c = 0;
      double sl = 0.;
      BGL_FORALL_VERTICES_T(v, G, Graph) {
         if( compnum[v] == i ) {
            rhs -= NC[v];
            b += NC[v];
            typename GraphTraits::out_edge_iterator it, it_end;
            for(tie(it, it_end) = out_edges(v, G); it != it_end; it++) {
               r -= VM[*it];
               sl += S.value(VM[*it]);
               if( S.value(VM[*it]) > 1. - LEPS )
                      c++;
                  }
               }
      }


      if( ((b - c) % 2) && (double(b) - sl < 1. - LEPS) ) {
         //FIXME
         r.normalize();
         if( r.size() > 0 ) {
            S.add_basic_constraint(r>=rhs);
            numfound++;
         }
      }
   }
   if(S.configuration("MATCHING_Debug_Out")=="true")
      cerr << "h: numfound:" << numfound << endl;
   return numfound;
}

template <typename Graph>
typename MATCHING<Graph>::status MATCHING<Graph>::standard_separation(subproblem& S) {
   if(S.configuration("MATCHING_Debug_Out")=="true")
      cerr << "MATCHING::standard_separation\n";

   if( heuristicBlossomSeparation(S) )
      return constraint_found; 
   if( perfect && exactBlossomSeparation(S) )
      return constraint_found;
   else
      return no_constraint_found;
}

template <typename Graph>
typename MATCHING<Graph>::status MATCHING<Graph>::fast_separation(subproblem& S) {
   if(S.configuration("MATCHING_Debug_Out")=="true")
      cerr << "MATCHING::fast_separation\n";

   if( heuristicBlossomSeparation(S) )
      return constraint_found;
   else
      return no_constraint_found;
}


template <typename Graph>
void MATCHING<Graph>::info() {
   cout<<"Configurations:\n";
   cout<<"MATCHING_Debug_Out [true|false]\n";
};
#endif

---