atour.cc

#ifndef ATOUR_CONSTRAINT_CC
#define ATOUR_CONSTRAINT_CC


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

#define VALONE 1000000

template <typename Graph>
class ATOUR<Graph>::dir_degree_equality : public cons_obj {
private:
  const Graph& G;
  vertex_descriptor v;
  double deg;
  var_map<edge_descriptor>& VM;
  bool in;
  
public:

   dir_degree_equality(subproblem& S_,
                  const Graph& G_, vertex_descriptor v_, var_map<edge_descriptor>& VM_,
                   double d, bool i_)
       /*{create The constraint $\sum_{e \in \delta(v)} x_e = d$, where $\delta(v)$ are the 
          edges which are adjacent to v and for which a variable was defined.}*/
     : cons_obj(Equal, d), G(G_), VM(VM_) { 
     v=v_; deg=d, in=i_; 
   };

   virtual
      ~dir_degree_equality()
      {};
  
   vertex_descriptor gnode() { return v; };
  
   virtual void non_zero_entries(row& r) {
      if(in) {
         typename GraphTraits::out_edge_iterator it, it_end;
         for (tie(it, it_end) = out_edges(v,G); it != it_end; it++) {
                  r+=VM[*it];
         }
      } else {
         typename GraphTraits::in_edge_iterator it, it_end;
         for (tie(it, it_end) = in_edges(v,G); it != it_end; it++) {
            r+=VM[*it];
         }
      };
   };
};

template <typename Graph>
class ATOUR<Graph>::dir_cutset_inequality : public cons_obj {
private:
   Graph& G;
   var_map<edge_descriptor>& VM;
   map<vertex_descriptor, bool> cut;
   double rhs;
  
public:
  
   dir_cutset_inequality(Graph& Gt, 
      var_map<edge_descriptor>& VM_, 
      list<vertex_descriptor>& L) 
      : cons_obj(Less, L.size()-1), G(Gt), VM(VM_) {
      typename GraphTraits::vertex_iterator v_it, v_it_end;
      for(tie(v_it, v_it_end) = vertices(G); v_it != v_it_end; v_it++)
         cut[*v_it] = false;
      typename list<vertex_descriptor>::iterator it;
      for (it=L.begin(); it != L.end(); it++)
         cut[*it]=true;
   }

   void non_zero_entries(row& r) { 
      typename GraphTraits::edge_iterator it, it_end;
      for(tie(it, it_end) = edges(G); it != it_end; it++) {
         if(coeff(*it)==1) r+=VM[*it];
      }
   }

   virtual
   ~dir_cutset_inequality() {};

   double coeff(edge_descriptor e) { 
      if((cut[source(e, G)]) && (cut[target(e, G)])) return 1;
      return(0);
   };
};

template <typename Graph>
ATOUR<Graph>::ATOUR(Graph& G_, var_map<edge_descriptor>& VM_)
   : G(G_), VM(VM_) {
   feps = 1.0e-4;
   //Graph has to provide in_edges() and out_edges()
   BOOST_CONCEPT_ASSERT((BidirectionalGraphConcept<Graph>));
   //Graph has to provide vertices() and edges()
   BOOST_CONCEPT_ASSERT((VertexAndEdgeListGraphConcept<Graph>));
   typename GraphTraits::edge_iterator it, it_end;
   for (tie(it, it_end) = edges(G); it != it_end; it++)
      if(VM[*it].type() != Vartype_Binary) {
         cerr << "atour.cc: All variables in VM need to be Vartype_Binary" << endl;
         exit(1);
      }
};

template <typename Graph>
void ATOUR<Graph>::init(subproblem& S) {
   if (S.configuration("ATOUR_Debug_Out")=="true")
      cerr << "ATOUR::init\n";
   typename GraphTraits::vertex_iterator it, it_end;
   for (tie(it, it_end) = vertices(G); it != it_end; it++) {
      S.add_basic_constraint(new dir_degree_equality(S,G,*it,VM,1, true));
      S.add_basic_constraint(new dir_degree_equality(S,G,*it,VM,1, false));
   }
}


template <typename Graph>
typename ATOUR<Graph>::status ATOUR<Graph>::heuristic_separation(subproblem& S) {
   typedef adjacency_list< vecS, vecS, undirectedS > uGraph;
   typedef graph_traits<uGraph>::edge_descriptor uEdge;
   typedef graph_traits<uGraph>::vertex_descriptor uVertex_descriptor;
   
   // Preparation
   int i=0;
   map<edge_descriptor, int> val;
   typename GraphTraits::edge_iterator e_it, e_it_end;
   for(tie(e_it, e_it_end) = edges(G); e_it != e_it_end; e_it++){
      val[*e_it] = max ((int) (VALONE*S.value(VM[*e_it])), 0);
   }
  
   // Heuristic
   uGraph H;
   map<uVertex_descriptor, vertex_descriptor> HtoG;
   map<vertex_descriptor, uVertex_descriptor> GtoH;
  
   typename GraphTraits::vertex_iterator v_it, v_it_end;
   for(tie(v_it, v_it_end) = vertices(G); v_it != v_it_end; v_it++) {
      GtoH[*v_it] = add_vertex(H);  
      HtoG[GtoH[*v_it]] = *v_it;
   }

   for(tie(e_it, e_it_end) = edges(G); e_it != e_it_end; e_it++) {
      if(val[*e_it]>=VALONE-5) add_edge( GtoH[source(*e_it, G)], GtoH[target(*e_it, G)], H);
   }

   typedef map<uVertex_descriptor, int> CC_t;
   CC_t CC;
   associative_property_map<CC_t> pm_CC(CC);

   int k = connected_components(H, pm_CC);

   typename graph_traits<uGraph>::vertex_iterator uv_it, uv_it_end;
   if(k>1) {
      for(int j=0;j<k;j++) {
         list<vertex_descriptor> CCL_one_element;
         for(tie(uv_it, uv_it_end) = vertices(H); uv_it != uv_it_end; uv_it++)
            if( pm_CC[*uv_it] == j )
               CCL_one_element.push_back(HtoG[*uv_it]);
         cons_obj* c=new dir_cutset_inequality(G,VM,CCL_one_element);
         if (c->violation(S) > feps) { 
                  i++;
                  S.add_basic_constraint(new dir_cutset_inequality(G,VM,CCL_one_element), Dynamic); 
         }
         delete c;
      }
   }

   // Heuristic succeeded?
   if(i>0) { 
      return constraint_found; 
   }
   return no_constraint_found;
}

template <typename Graph>
typename ATOUR<Graph>::status ATOUR<Graph>::exact_separation(subproblem& S){
   // Create undirected copy of Graph G 
   // This is needed to perform min-cut
   // I think this workaround should soon be removed 

   typedef adjacency_list< vecS, vecS, undirectedS > uGraph;
   typedef graph_traits<uGraph>::edge_descriptor uEdge;
   typedef graph_traits<uGraph>::vertex_descriptor uVertex_descriptor;

   uGraph uG(num_vertices(G));

   std::map<vertex_descriptor, unsigned int> GtoUG;
   vector <vertex_descriptor> UGtoG(num_vertices(G));  
   int n=0;
   vertex_descriptor v;
   map<uVertex_descriptor, int> vertex_index;
   BGL_FORALL_VERTICES_T( v,G,Graph ){
      UGtoG[n] = v;
      GtoUG[v] = n;
      vertex_index[GtoUG[v]] = n;
      n++;
   }

   edge_descriptor e;
   uEdge ue;
   map<uEdge, int> valUG;
   BGL_FORALL_EDGES_T( e,G,Graph ){
      if(!edge(GtoUG[source(e,G)], GtoUG[target(e,G)],uG).second){
      ue = add_edge( GtoUG[source(e,G)], GtoUG[target(e,G)], uG ).first;      
      valUG[ue] = max ((int) (VALONE*S.value(VM[e])), 0);
      }else
         valUG[edge(GtoUG[source(e,G)], GtoUG[target(e,G)],uG).first] += max ((int) (VALONE*S.value(VM[e])), 0);
   }
 
 
   // exact separation   
   MC<uGraph> mc;
   MIN_CUT<uGraph> *min_cut = new MIN_CUT<uGraph>(uG, valUG);   // calc min-cut of uG
   min_cut->run();
   mc = min_cut->minCut;
   delete min_cut;
   std::list<vertex_descriptor> cut; // min-cut of G
   BOOST_FOREACH( long unsigned int ui, mc.nodes )
     cut.push_back(UGtoG[ui]);

   int i = 0;
   cons_obj* c = new dir_cutset_inequality(G, VM, cut);
   if (c->violation(S) > feps) {
       S.add_basic_constraint(new dir_cutset_inequality(G, VM, cut), Dynamic);
       i++;
   }
   delete c;
   if(i>0) {
      //cout<<"Exact CUT found inequality\n";
      return constraint_found;
   }
   else
      return no_constraint_found;
}

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

   return heuristic_separation(S);
}

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

   if(heuristic_separation(S) == no_constraint_found)
      return exact_separation(S);
   else
      return constraint_found;

}

template <typename Graph>
typename ATOUR<Graph>::status ATOUR<Graph>::feasible(solution& S) {
   if (S.originating_subproblem()->configuration("ATOUR_Debug_Out")=="true")
      cerr << "ATOUR::feasible\n";
   map<vertex_descriptor, vertex_descriptor> GtoH;
   Graph H;
   typename GraphTraits::vertex_iterator it, it_end;
   for (tie(it, it_end) = vertices(G); it != it_end; it++)
      GtoH[*it]=add_vertex(H);
   typename GraphTraits::edge_iterator e_it, e_it_end;
   for (tie(e_it, e_it_end) = edges(G); e_it != e_it_end; e_it++)
      if(S.value(VM[*e_it])>0.5)
         add_edge(GtoH[source(*e_it, G)], GtoH[target(*e_it,G)], H);
   vector<int> component(num_vertices(H));
   //Test if graph H is connected. If H has just 1 connected component, it is.
   if (connected_components(H, &component[0]) == 1) 
      return feasible_solution;
   else return infeasible_solution;
}

template <typename Graph>
void ATOUR<Graph>::info() {
   cout<<"Configurations:\n";
   cout<<"ATOUR_Debug_Out [true|false]\n";
};

#undef VALONE

#endif

---