tour.cc

#ifndef SCIL_TOUR_CC
#define SCIL_TOUR_CC

#include<scil/scil.h>
#include<LEDA/graph.h>
#include<LEDA/graph_alg.h>
#include<LEDA/set.h>
#include<scil/constraints/tour.h>

using namespace SCIL;
using namespace LEDA;

#define VALONE 1000000

class TOUR::degree_equality : public cons_obj
{
  private:
  const graph& G;
  node v;
  double deg;
  scil_map<edge>& VM;

  public:

  // The constructor saves the references to the graph and to the scil_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_, node v_, scil_map<edge>& VM_,
                  double d)
  : cons_obj(Equal, d), G(G_), VM(VM_) { 
    v=v_; deg=d; 
  };
  
  virtual
  ~degree_equality()
  {};
  
  node 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 e;
    forall_inout_edges(e, v) {
      r+=VM(e);
    }
  };
}; // END_OF_CLASS_DEGREE_EQUALITY

class TOUR::cutset_inequality : public cons_obj {
  private:
  graph& G;
  scil_map<edge>& VM;
  node_array<bool> cut;
  double rhs;
  
  public:
  
  cutset_inequality(graph& Gt, scil_map<edge>& VM_, list<node>& L) 
  : cons_obj(Greater, 2), G(Gt), VM(VM_) {
    node v;
    cut.init(G, false);
    forall(v,L) cut[v]=true;
  }

  void non_zero_entries(row& r) { 
    edge e;
    forall_edges(e,G) {
      if(coeff(e)==1) r+=VM(e);
    }
  }

  virtual
  ~cutset_inequality() {};

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

TOUR::TOUR(graph& G_, scil_map<edge>& VM_)
  /*{\Mcreate Ensures that the edges in VM form a Hamiltonian cycle.\\
    For a demo, see: gtsp.}*/
  : G(G_), VM(VM_) { 
};

void TOUR::init(subproblem& S) {
  node v;
  forall_nodes(v, G) {
    S.add_basic_constraint(new degree_equality(G,v,VM,2));
  }
};
  
TOUR::status TOUR::separate(subproblem& S) {
  // Preperation : Read the current LP-values of all edges
  //               We multiply the value with VALONE and round to the next integer
  int i=0;
  edge_array<int> val(G);
  edge e;
  forall_edges(e, G) 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;
  node u;
  node_array<node> HtoG(H, G.number_of_nodes(), u);
  node_array<node> GtoH(G);
  
  // Copy all nodes of G
  forall_nodes(u, G) { GtoH[u]=H.new_node(); HtoG[GtoH[u]]=u; }
  
  // Copy the edges with value close to 1
  forall_edges(e, G) if(val[e]>=VALONE-5)
    H.new_edge(GtoH[source(e)], GtoH[target(e)]);
  
  // Compute the connected components
  node_array<int> CC(H);
  int k=COMPONENTS(H,CC);
  list<node> CCL[k];
  forall_nodes(u,H) CCL[CC[u]].append(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->violated(S)) { 
        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(S.configuration("TOUR_heuristic_only")=="true") 
  //  return no_constraint_found;

  // exact separation : We compute the minimal cut in G and check the
  //                    corresponding subtour elimination constraint for
  //                    violation
  list<node> mc=MIN_CUT(G, val);
  cons_obj* c=new cutset_inequality(G,VM,mc);
  if (c->violated(S)) {
    i++; S.add_basic_constraint(c, Dynamic);
  } else delete c;

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

  return no_constraint_found;

  H.del_all_edges();

  // Copy the edges with fractional value
  forall_edges(e, G) if((val[e]>=5) && (val[e]<=VALONE-5))
    H.new_edge(GtoH[source(e)], GtoH[target(e)]);
  
  // Compute the connected components
  {
    node_array<int> CC(H);
    int k=COMPONENTS(H,CC);
    list<node> CCL[k];
    forall_nodes(u,H) CCL[CC[u]].append(HtoG[u]);

    // Create the corresponding comb constraints and check for
    // violation
    node v;
    node_array<bool> R(G, false);
    if(k>1) {
      for(int j=0;j<k;j++) {
        row r;
        forall(v, CCL[j]) R[v]=true;
        double d=0;
        int t=0;
        forall_edges(e,G) {
          if((R[source(e)]) && (R[target(e)])) { 
            d+=S.value(VM(e));
            r+=VM(e);
          }
          if((R[source(e)]!=R[target(e)]) && (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";
        };
        forall(v, CCL[j]) R[v]=false;
      }
    }
  }
  if(i>0) {
    return constraint_found;
  }

  return no_constraint_found;
}; // END_OF_SEPARATE


TOUR::status TOUR::feasible(solution& S) {
  // 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
  node_array<node> GtoH(G);
  graph H;
  node u;
  forall_nodes(u, G)
    GtoH[u]=H.new_node();
  edge e;
  forall_edges(e,G) {
    if(S.value(VM(e))>0.5) H.new_edge(GtoH[source(e)], GtoH[target(e)]);
  }

  // check degree equalities
  forall_nodes(u,H) if (H.degree(u)!=2) return infeasible_solution;

  // check connectedness
  if (Is_Connected(H)) return feasible_solution;
  return infeasible_solution;
}; // END_OF_FEASIBLE

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

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

#undef VALONE

#endif

---