dnatag.cc

#include<scil/scil.h>
#include<boost/foreach.hpp>
#include<boost/graph/iteration_macros.hpp>
#include<boost/graph/adjacency_list.hpp>

using namespace SCIL;
using namespace boost;

#define P 100
#define T 50
#define Pr 0.02

typedef boost::adjacency_list<vecS,vecS,bidirectionalS> Graph;
typedef boost::graph_traits<Graph>::vertex_descriptor vertex_descriptor;
typedef boost::graph_traits<Graph>::edge_descriptor edge_descriptor;

class PrimerVariable : public var_obj {
  cons* C;
  std::list<int> PV;

  public:

  PrimerVariable(bool* B_, cons* C_) : var_obj(1.,0.,1.,Vartype_Binary) {
    C=new cons[P];
    for(int i=0; i<P; i++) {
      C[i]=C_[i];
      if(B_[i]) PV.push_back(i);
    };
  };

  PrimerVariable(cons* C_) : var_obj(100.,0.,1.,Vartype_Binary) {
    C=new cons[P];
    for(int i=0; i<P; i++) {
      C[i]=C_[i];
      PV.push_back(i);
    };
  };

  ~PrimerVariable() {
    delete[] C;
  };

  void non_zero_entries(column& Co) {
    cout<<PV.size()<<" HIER\n";
    BOOST_FOREACH(int i, PV) {
      Co+=(column) C[i];
    }
  }

};

double Solve_Pricing(Graph& G, double* C, bool* B, vector<int> N) { /*node_array<int> N*/
  ILP_Problem ILP(Optsense_Min);
  
  //ILP.set_silent();

  var_map<vertex_descriptor> VN;
  row r;
  BGL_FORALL_VERTICES(u,G,Graph) {
    if(N[u]==-1) {
      VN[u]=ILP.add_binary_variable(0.);
      r+=VN[u];
    } else { 
      VN[u]=ILP.add_binary_variable((double)-C[N[u]]);
      r-=VN[u];
    }
  }
  ILP.add_basic_constraint(r==0);  

  var_map<edge_descriptor> VE;
  BGL_FORALL_EDGES(e,G,Graph) {
    VE[e]=ILP.add_variable(0.,0.,0.,Vartype_Integer);
    ILP.add_basic_constraint((row) VE[e]>=VN[source(e,G)]+VN[target(e,G)]-1);
  }

  BGL_FORALL_VERTICES(u,G,Graph) if(in_degree(u,G)+out_degree(u,G) != 0) {
    row r1;
    BGL_FORALL_INEDGES (u,e,G,Graph) r1+=(row) VE[e];
    BGL_FORALL_OUTEDGES(u,e,G,Graph) r1+=(row) VE[e];
    ILP.add_basic_constraint(r1<=1);
  }

  ILP.optimize();

  double D=0;
  BGL_FORALL_VERTICES(u, G, Graph) if(N[u]!=-1) {
    if(ILP.get_solution(VN[u])==1) {
      B[N[u]]=1;
      D+=C[N[u]];
    } else B[N[u]]=0;
  };
  cout<<D<<" SOLIUTION\n";
  return D;
};

class PrimerPricer : public sym_constraint {
  private:
    Graph& G;
    cons* C;
    vector<int>& PN;
    bool no_more_pricing;

  public:

  PrimerPricer(Graph& G_, cons* C_, vector<int>& PN_) : G(G_), C(C_), PN(PN_) {
    no_more_pricing=false;
  };

  status feasible(subproblem& s) {
    return feasible_solution;
  };

  status LP_pricing(subproblem& s) {
    if(no_more_pricing) return no_variable_found;
    cout<<"HIER\n";
    double Co[P];
    bool B[P];
    for(int i=0; i<P; i++) Co[i]=s.value(C[i])+0.001;
    double d=Solve_Pricing(G, Co, B, PN);
    //for(int i=0; i<P; i++) cout<<B[i]<<" "; cout<<endl;
    cout<<d<<"----------------- ADD AN VARIABLE"<<endl;
    if(d>1.05) { s.add_variable(new PrimerVariable(B, C));
      return variable_found;
    }
    no_more_pricing=true;
    return no_variable_found;
  };
};

int main() {
  // Make Input

  vertex_descriptor PN[P];
  vertex_descriptor TN[T];

  Graph G;
  srand ( time(NULL) );    /* initialize random number generator */ 
  double p;


  for(int i=0; i<P; i++) PN[i]=add_vertex(G);

  for(int j=0; j<T; j++) {
    TN[j]=add_vertex(G);
    for(int i=0; i<P; i++) {
      p = (double)rand() / RAND_MAX;

      if(p<=Pr) add_edge(PN[i],TN[j],G);
    }
  };
  

  // Optimize

  vector<int> PNN(num_vertices(G));
  for(int i=0; i<P; i++) PNN[PN[i]]=i;
  for(int j=0; j<T; j++) PNN[TN[j]]=-1;

  ILP_Problem ILP(Optsense_Min); //, true);

  cons CS[P];
  for(int i=0; i<P; i++) CS[i]=ILP.add_basic_constraint(Greater,1);

  ILP.add_variable(new PrimerVariable(CS));

  ILP.add_sym_constraint(new PrimerPricer(G, CS, PNN));
    
  ILP.optimize();

  cout<<"OPTIMIZED\n";

  return 0;

}

---