spantree.cc

#ifndef SPAN_TREE_CC
#define SPAN_TREE_CC

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

using namespace SCIL;
using namespace LEDA;

#define VALONE 10000

class SpanTree::st_cutset_inequality : public cons_obj {
private:
  graph& G;
  scil_map<edge>& VM;
  node_array<bool> cut;
  double rhs;
  
public:
  
  st_cutset_inequality(graph& Gt, scil_map<edge>& VM_, list<node>& L) 
  : cons_obj(Greater, 1), 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
  ~st_cutset_inequality() {};

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


SpanTree::SpanTree(graph& G_, scil_map<edge>& VM_)
  : G(G_), VM(VM_)
{ };

void SpanTree::init(subproblem& S) {
  row r;
  edge e;
  forall_edges(e, G) { r+=VM(e); }
  S.add_basic_constraint(r == G.number_of_nodes()-1);
  
  node v;
  GtoH.init(G);
  GtoH1.init(G);
  GtoH2.init(G);
  forall_nodes(v, G) {
    GtoH[v]=H.new_node();
  };

  forall_edges(e, G) {
    GtoH1[e]=H.new_edge(GtoH[source(e)], GtoH[target(e)]);
    GtoH2[e]=H.new_edge(GtoH[target(e)], GtoH[source(e)]);
  };

  son=H.new_node();
  tan=H.new_node();
  forall_nodes(v, G) {
    H.new_edge(son, GtoH[v]);
    H.new_edge(GtoH[v], tan);
  };
  
  Cap.init(H);
}

void SpanTree::min_cut(node u, edge_array<int>& C, edge_array<int>& F, 
                       node_array<bool>& R, int& k) {
  R[u]=true;
  k++;
  edge e;
  forall_out_edges(e, u) if((F[e]!=C[e]) && (!R[target(e)])) {
    min_cut(target(e), C, F, R, k);
  };
  forall_in_edges(e, u) if((F[e]!=0) && (!R[source(e)])) {
    min_cut(source(e), C, F, R, k);
  };
};

SpanTree::status SpanTree::separate(subproblem& S) {

  cout<<"Separation SpanTree\n";

  status sta=no_constraint_found;
  // Preperation
  int j, k;
  node_array<int> Con(H);
  edge e,f;
  forall_edges(e, G) { 
    j=(int) (VALONE*S.value(VM(e)));
    if(j<0) j=0;
    Cap[GtoH1[e]]=j;
    Cap[GtoH2[e]]=j;
    Con[GtoH[source(e)]]+=j;
    Con[GtoH[target(e)]]+=j;
  };
  
  forall_out_edges(e, son) {
    Cap[e]=leda_max(Con[target(e)]-2*VALONE,0);
  };
  forall_in_edges(e, tan) {
    Cap[e]=leda_max(2*VALONE-Con[source(e)],0);
  };
  
  if(sta==constraint_found) return sta;            
  
  edge_array<int> Flow(H);
  forall_out_edges(e, son) {
    k=Cap[e];
    Cap[e]=1000*VALONE;
    j=MAX_FLOW(H, son, tan, Cap, Flow);
    // Make Constraint
    node_array<bool> R(H, false);
    j=0;
    min_cut(son, Cap, Flow, R, j);
    row r;
    double d=0;
    forall_edges(f, G) {
      if((R[GtoH[source(f)]]) && (R[GtoH[target(f)]])) {
        r+=VM(f);
        d+=S.value(VM(f));
      };
      }
    if(d>((double) j)-1.99) {
      S.add_basic_constraint(r<=j-2, Dynamic);
      sta=constraint_found;
    };
    Cap[e]=k;
  };
  return sta;
}

SpanTree::status SpanTree::feasible(solution& S) {
  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)]);
  }
  node_array<int> C(H);
  int i=COMPONENTS(H, C);
  if (i==1) return feasible_solution;
  cout<<"integer feasible, but infeasible\n";
  return infeasible_solution;
}

#undef VALONE

#endif

---