tour.cc

#ifndef SCIL_TOUR_CC
#define SCIL_TOUR_CC

using namespace SCIL;
using namespace std;

#define VALONE 1000000

template <typename Graph>
class TOUR<Graph>::degree_equality : public cons_obj
{
   private:
      const Graph& G;
      vertex_descriptor v;
      double deg;
      var_map<edge_descriptor>& VM;

   public:

   // The constructor saves the references to the graph and to the var_map<edge>
   // Furthermore we save the right hand side and the given node of the constraint.
   // We initialize the sense and the right hand side with equal and the given
   // right hand side.
   degree_equality(const Graph& G_, vertex_descriptor v_, var_map<edge_descriptor>& VM_,
                  double d)
      : cons_obj(Equal, d), G(G_), VM(VM_) { 
      v=v_; deg=d; 
   };
  
   virtual
   ~degree_equality()
   {};
  
   vertex_descriptor gnode() { return v; };

   // To generate the row, we can simply iterate over all edges adjacent to
   // the node of the degree constraint and add the appropriate variables.
   virtual void non_zero_entries(row& r) {
      edge_descriptor e;
      BGL_FORALL_OUTEDGES_T(v, e, G, Graph) {
         r+=VM[e];
      }
   };
}; // END_OF_CLASS_DEGREE_EQUALITY

template <typename Graph>
class TOUR<Graph>::cutset_inequality : public cons_obj {
   private:
      Graph& G;
      var_map<edge_descriptor>& VM;
      map<vertex_descriptor, bool> cut;
      double rhs;
  
   public:
  
      cutset_inequality(Graph& Gt, var_map<edge_descriptor>& VM_, list<vertex_descriptor>& L) 
         : cons_obj(Greater, 2), G(Gt), VM(VM_)
      {
         vertex_descriptor v;
         BGL_FORALL_VERTICES_T(v, G, Graph)
            cut[v] = false;
         BOOST_FOREACH(v, L)
            cut[v] = true;
      };

      void non_zero_entries(row& r) { 
         edge_descriptor e;
         BGL_FORALL_EDGES_T(e, G, Graph) {
            if(coeff(e)==1) r+=VM[e];
         }
      }

   virtual
      ~cutset_inequality() {};

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

template <typename Graph>
TOUR<Graph>::TOUR(Graph& G_, var_map<edge_descriptor>& VM_)
   : G(G_), VM(VM_) { 
   edge_descriptor e;
   feps = 1.0e-4;
   BGL_FORALL_EDGES_T(e, G, Graph)
      if(VM[e].type() != Vartype_Binary) {
         cerr << "tour.cc: All variables in VM need to be Vartype_Binary" << endl;
         exit(1);
      }
};

template <typename Graph>
void TOUR<Graph>::init(subproblem& S) {
   if(S.configuration("TOUR_Debug_Out")=="true")
      cerr << "TOUR::init\n";
   vertex_descriptor v;
   BGL_FORALL_VERTICES_T(v, G, Graph) {
      S.add_basic_constraint(new degree_equality(G,v,VM,2));
   }
};
  

template <typename Graph>
typename TOUR<Graph>::status TOUR<Graph>::standard_separation(subproblem& S) {
   if(S.configuration("TOUR_Debug_Out")=="true")
      cerr << "TOUR::standard_separation\n";
   return separation(S, false);
}

template <typename Graph>
typename TOUR<Graph>::status TOUR<Graph>::fast_separation(subproblem& S) {
   if(S.configuration("TOUR_Debug_Out")=="true")
      cerr << "TOUR::fast_separation\n";
   return separation(S, true);
}



template <typename Graph>
typename TOUR<Graph>::status TOUR<Graph>::separation(subproblem& S, bool fast_separation) {
   // Preparation : Read the current LP-values of all edges
   //               We multiply the values with VALONE and round to the next integer
   int i=0;
   map<edge_descriptor, int> val;
   edge_descriptor e;
   BGL_FORALL_EDGES_T(e, G, Graph) {
      val[e]=(int) (VALONE*S.value(VM[e]));
   }

   // Heuristic : We create a copy H of G, where all edges with LP-value less than
   //             1 are deleted from the graph. We compute the connected components
   //             of this graph and check the subtour elemination constraints
   //             corresponding to the components for violation
   Graph H;
   vertex_descriptor u;
   map<vertex_descriptor, vertex_descriptor> HtoG;
   map<vertex_descriptor, vertex_descriptor> GtoH;
  
   // Copy all nodes of G
   int index = 0;
   map<vertex_descriptor, int> vertex_index;
   BGL_FORALL_VERTICES_T(u, G, Graph) { 
      GtoH[u]=add_vertex(H);
      HtoG[GtoH[u]]=u;
      vertex_index[GtoH[u]] = index;
      index++;
   }
  
      // Copy the edges with value close to 1
      BGL_FORALL_EDGES_T(e, G, Graph) if(val[e]>=VALONE-5)
         add_edge(GtoH[source(e, G)], GtoH[target(e, G)], H);
  
      // Compute the connected components
      vector<int> CC(num_vertices(H));
      int k = connected_components(H, &CC[0]);
      list<vertex_descriptor> CCL[k];
      BGL_FORALL_VERTICES_T(u,H, Graph) CCL[CC[vertex_index[u]]].push_back(HtoG[u]);

      // Create the corresponding subtour elemination constraints and check for
      // violation
      if(k>1) {
         for(int j=0;j<k;j++) {
            cons_obj* c=new cutset_inequality(G,VM,CCL[j]);
            if (c->violation(S) > feps) { 
               i++;
               S.add_basic_constraint(c, Dynamic); 
            } else delete c;
         }
      }

      // If the heuristic found a violated constraint, we stop the separation
      if(i>0) { 
         cout<<i<<" Heuristic\n"; 
         return constraint_found; 
      }

   // If only the fast heuristic should be applied, we stop the separation  
   if(fast_separation) 
      return no_constraint_found;

   // exact separation : We compute the minimal cut in G and check the
   //                    corresponding subtour elimination constraint for
   //                    violation

   MC<Graph> mc;
   MIN_CUT<Graph> *min_cut = new MIN_CUT<Graph>(G, val);        // calc min-cut of uG
   min_cut->run();
   mc = min_cut->minCut;
   delete min_cut;
   list<vertex_descriptor> cut; // min-cut of G
   BOOST_FOREACH( vertex_descriptor v, mc.nodes )
     cut.push_back(v);


   cons_obj* c=new cutset_inequality(G,VM,cut);
   if (c->violation(S) > feps) {
      i++; S.add_basic_constraint(c, Dynamic);
   } else delete c;

   // return the correct value
   if(i>0) {
      return constraint_found;
   }

   BGL_FORALL_VERTICES_T(u, H, Graph)
      clear_vertex(u, H);

   // Copy the edges with fractional value
   BGL_FORALL_EDGES_T(e, G, Graph) 
      if((val[e]>=5) && (val[e]<=VALONE-5))
         add_edge(GtoH[source(e, G)], GtoH[target(e, G)], H);
  
   // Compute the connected components
   {
      vector<int> CC(num_vertices(H));
      int k=connected_components(H,&CC[0]);
      list<vertex_descriptor> CCL[k];
      BGL_FORALL_VERTICES_T(u,H, Graph) CCL[CC[vertex_index[u]]].push_back(HtoG[u]);

      // Create the corresponding comb constraints and check for
      // violation
      vertex_descriptor v;
      map<vertex_descriptor, bool> R;
      if(k>1) {
         BGL_FORALL_VERTICES_T(v, G, Graph)
            R[v] = false;
         for(int j=0;j<k;j++) {
            row r;
            BOOST_FOREACH(v, CCL[j]) R[v]=true;
            double d=0;
            int t=0;
            BGL_FORALL_EDGES_T(e,G, Graph) {
               if((R[source(e, G)]) && (R[target(e, G)])) { 
                  d+=S.value(VM[e]);
                  r+=VM[e];
               }
               if((R[source(e, G)]!=R[target(e, G)]) && (val[e]>VALONE-10)) {
                  d+=S.value(VM[e]);
                  r+=VM[e];
                  t++;
               };
            };
            if((t % 2==1) && (d>CCL[j].size()+(t-1)/2+0.01)) {
               i++;
               S.add_basic_constraint(r<=CCL[j].size()+(t-1)/2);
               cout<<"comb\n";
            };
            BOOST_FOREACH(v, CCL[j]) R[v]=false;
         }
      }
   }
   if(i>0) {
      return constraint_found;
   }

   return no_constraint_found;
}; // END_OF_SEPARATE


template <typename Graph>
typename TOUR<Graph>::status TOUR<Graph>::feasible(solution& S) {
   if(S.originating_subproblem()->configuration("TOUR_Debug_Out")=="true")
      cerr << "TOUR::feasible\n";
   // We construct a copy H of G, where we delete all edges with value different
   // from one. We check all degree equalities. If all degree equalities are
   // satisfied, we know that H is a Hamiltonian cycle if H is connected

   // Create H
   map<vertex_descriptor, vertex_descriptor> GtoH;
   Graph H;
   vertex_descriptor u;
   BGL_FORALL_VERTICES_T(u, G, Graph)
      GtoH[u]=add_vertex(H);
   edge_descriptor e;
   BGL_FORALL_EDGES_T(e,G, Graph) {
      if(S.value(VM[e])>0.5) add_edge(GtoH[source(e, G)], GtoH[target(e, G)], H);
   }

   // check degree equalities
   BGL_FORALL_VERTICES_T(u,H, Graph) if (out_degree(u, H)!=2) return infeasible_solution;

   // check connectedness
   vector<int> comp(num_vertices(H));
   int cc = connected_components(H,&comp[0]);
   if (cc == 1) return feasible_solution;
   return infeasible_solution;
}; // END_OF_FEASIBLE

template<typename Graph>
void TOUR<Graph>::info() {
   cout<<"Method: Cutting Planes\n\n";

   cout<<"Configurations:\n";
   cout<<"TOUR_heuristic_only [true|false]\n";
   cout<<"TOUR_Debug_Out [true|false]\n";
};

#undef VALONE

#endif

---