ilp_problem.cc

#include<scil/variable.h>
#include<scil/global.h>
#include<scil/var_obj.h>
#include<scil/cons_obj.h>
#include<scil/aba_variable.h>
#include<scil/aba_constraint.h>
#include<scil/subproblem.h>
#include<scil/ilp_problem.h>
#include<scil/row.h>
#include<scil/primal_heur.h>
#include<scil/constraints/quadref.h>
#include<scil/core/polynomial.h>
#include<scil/core/monomial_inst.h>
#include<scil/core/bool_inst.h>
#include<scil/core/boolfunction.h>
#include <tr1/tuple>
#include<iostream>

using namespace SCIL;

ABA_SUB* ILP_Problem::firstSub() {
  myroot=new subproblem(this);

  std::list<Coefficient>::const_iterator co;
  foreach( co, coeff_list ){
    myroot->set_coefficient( tr1::get<0>(*co), tr1::get<1>(*co), tr1::get<2>(*co) );
  }

  std::list<sym_constraint*>::const_iterator c;
  foreach( c, sym_constraints ) { 
    myroot->add_sym_constraint(*c);
  }
  return myroot;
};


const
std::string ILP_Problem::configuration(const std::string& s) const {
  std::map<std::string,std::string>::const_iterator it = ConfMap.find( s );

  if( it==ConfMap.end() )  
    return "undef";
  else
    return it->second;
};

void ILP_Problem::configuration(const std::string& s1, const std::string& s2) {
   ConfMap[s1] = s2;
   return;
};

void ILP_Problem::read_configuration_file(const char* s) {
  ifstream ifs(s);
  std::cerr << "reading configuration file " << s << endl;

  while(ifs.good()) {
     std::string s1, s2;
    if(ifs.peek()=='#') std::getline(ifs,s2);
    else {
       if( !ifs.good() )
          return;
      ifs >> s1;
      ifs >> s2;
      if( s1.length() > 0 && s2.length() > 0 )
         configuration(s1, s2);
    }
  }
};


ILP_Problem::ILP_Problem (Optsense optSense, bool price, int _numcon, const char* problemName) : 
  ABA_MASTER (problemName, true, price, (optSense == Optsense_Min ? ABA_OPTSENSE::Min : ABA_OPTSENSE::Max)),numcon(_numcon) //, bran(NULL)
{
  ConfMap["undef"] = "undef";

  os = optSense;
  ABA_BUFFER < ABA_VARIABLE * >variables (this, 0);
  ABA_BUFFER < ABA_CONSTRAINT * >temp (this, 0);
  ABA_BUFFER < ABA_CONSTRAINT * >constraints (this, 0);
  initializePools (constraints, temp, variables, 
                   minimal_size, minimal_size, true);
  conPool ()->increase (minimal_size);
  cutPool ()->increase (minimal_size);
  nvars = 0;
  nconss = 0;
  primal_heur = NULL;
  implic = NULL;
  objint = 1;
  spb = false;
  monomialinstance = NULL;
  boolinstance = NULL;
  mon_splitstrat = NULL;
  bool_splitstrat = NULL;
  qr = new QuadRef(*this);
  qr_on_separate = true;
  qrType = SQK2_SQK3;
};



ILP_Problem::~ILP_Problem () {
   if( primal_heur != NULL )
      delete primal_heur;
   if( implic != NULL ) {
      delete implic;
      implic = NULL;
   }
   for (std::list<sym_constraint*>::iterator sc = sym_constraints.begin(); sc != sym_constraints.end(); sc++)
      delete *sc;
   std::list<boolfunction*>::iterator it;
   for (it = bf_list.begin(); it != bf_list.end(); it++)
         delete *it;
   for (it = bf_toDelete.begin(); it != bf_toDelete.end(); it++)
      delete *it;
};


void ILP_Problem::initializeOptimization ()
{   
  if(objint) objInteger(1);
  if(numcon >= 0) {
    maxConAdd(numcon);
    maxConBuffered(numcon);
  }
  if( spb )
     primalBound(pb);
};


void ILP_Problem::initializeParameters ()
{
  //if(exists(".abacus")) {
    //readParameters (".abacus");
    //}
    //Parameters are read from $ABACUS_DIR/.abacus only!
};

var ILP_Problem::add_binary_variable(double obj, Activation a)
{
   return add_variable(obj, 0, 1, Vartype_Binary, a);
}

var ILP_Problem::add_variable (double obj, double lBound, double uBound,
                               Vartype t, Activation a)
{
  if (nvars == varPool ()->size ())
    varPool ()->increase (2 * nvars);
  
  var_obj * v = new var_obj (obj, lBound, uBound, t);
  return add_variable (v, a);
};


var ILP_Problem::add_variable (var_obj * v, Activation a)
{
  v->init (*this, nvars, a);
  varPool ()->insert (v->AVar ());
  nvars++;
  if(!isInteger(v->obj()) || v->vt != Vartype_Integer) objint = 0;
  return v;
};

void ILP_Problem::updateBestSolution (solution& Sol_)
{
  Sol=Sol_;
}

cons
ILP_Problem::add_basic_constraint (cons_obj * c, Activation a, quadRefStatus qrStatus)
{
   c->init (*this, nconss, a);
   if(qrStatus != CHECK)
      c->set_qrStatus(qrStatus);
   qr->addCons(c);
   nconss++;
   if (conPool ()->number () == conPool ()->size ())
      conPool ()->increase (2 * conPool ()->size ());
   conPool ()->insert (c->Acons ());
   return c;
};


cons ILP_Problem::add_basic_constraint(cons_sense s, double r, Activation a) {
  /*{\Mop adds a primitive basic constraint to the root of the 
    BCP-tree. |s| specifies the sense (one of |Less|, |Equal|, or 
    |Greater|) and |r| the right-hand-side of the constraint.}*/
  cons_obj* I=new cons_obj();
    I->set_sense(s);
    I->set_rhs(r);
    add_basic_constraint(I, a);
    return I;
}

void ILP_Problem::add_sym_constraint(sym_constraint* c) {
  /*{\Mop adds the symbolic constraint |c| to the model.}*/
  sym_constraints.push_back(c);
  return;
}

void ILP_Problem::add_pol_constraint(pol_constraint c) {
   polynomial pol = c.p;
   std::list<monomial> mons = pol.get_p();
   row r;
   std::list<monomial>::iterator pos;
   std::list<monomial>::const_iterator mo;
   foreach(mo, mons) {
      pos = lower_bound(monomials.begin(), monomials.end(), *mo);
      //construct constraint using linearization variable and monomial coeff
      if(pos != monomials.end() && *pos == *mo )
         r += pos->get_L() * mo->get_coeff();
      else
         r += add_polynomial(*mo, false);
   }
   if( c.s == Less )
      add_basic_constraint(r<=c.rhs);
   else if( c.s == Greater )
      add_basic_constraint(r>=c.rhs);
   else
      add_basic_constraint(r==c.rhs);
   return;
}

row ILP_Problem::add_polynomial(polynomial p, bool update) {
   p.normalize();
   row r;    
   std::list<monomial> mons = p.get_p();
   std::list<monomial>::iterator mon;
   std::list<monomial>::iterator pos;
   std::vector<var>::iterator varit;
   std::vector<var> varlist;
   foreach(mon, mons) { 
      varlist = mon->get_A();
      foreach(varit, varlist){
         if(!(varit->type()==Vartype_Binary)){
            cerr<<"All variables in a polynomial have to be binary!"<<endl;
            exit(Fatal);
         }
      }
      if( mon->size() == 1 ) {
         //monomial is linear
         if(update)
            set_obj_coefficient((mon->get_A())[0], get_obj_coefficient((mon->get_A())[0]) + mon->get_coeff());
         r += mon->get_coeff() * (mon->get_A())[0];
         continue;
      }
      if( monomials.size() > 0 ) {
         //search for the monomial in global monomial list
         pos = lower_bound(monomials.begin(), monomials.end(), *mon);
         if(pos!=monomials.end() && *mon == *pos) {
            r+=mon->get_coeff() * pos->get_L();
            if(update){
               //if monomial exists the coefficient is updated
               pos->set_coeff(pos->get_coeff() + mon->get_coeff());
               set_obj_coefficient(pos->get_L(), pos->get_coeff());
            } 
         }else {
            //add monomial to global list and create linearization variable
            double coeff = mon->get_coeff();
            if(!update)
               mon->set_coeff(0);
            mon->linearize(*this);
            r+=coeff * mon->get_L();
            monomials.insert(pos, *mon);
            //if the monomial consists of two variables add it to QuadRef
            if(mon->size() == 2){
               qr->hashPair((mon->get_A())[0], (mon->get_A())[1], (mon->get_L()));
            }
         }
      }else{
         //add monomial to global list and create linearization variable
         double coeff = mon->get_coeff();
         if(!update)
            mon->set_coeff(0);
         mon->linearize(*this);
         r+=coeff * mon->get_L();
         monomials.push_back(*mon);
         //if the monomial consists of two variables add it to QuadRef
         if(mon->size() == 2){
            qr->hashPair((mon->get_A())[0], (mon->get_A())[1], (mon->get_L()));
         }

      }
   }
   return r;
}

row ILP_Problem::add_polynomial(monomial m, bool update) {
   return add_polynomial(polynomial(m), update);
}

var ILP_Problem::add_boolfunction(boolfunction* bf, double c) {
   var linVar;
   std::list<boolfunction*>::iterator pos;
   if(bf->get_op() == BASIC){
      if(bf->is_negated()){
         linVar = add_binary_variable(c);
         add_basic_constraint(linVar + bf->get_v() == 1, Static);
      } else {
         set_obj_coefficient(bf->get_v(), get_obj_coefficient(bf->get_v()) + c);
         linVar = bf->get_v();
      }
      bf_toDelete.push_back(bf);
   } else {
      bf->simplify();
      pos = lower_bound(bf_list.begin(), bf_list.end(), bf, bfcompare_less);
      if(pos!=bf_list.end() && bfcompare_equal(*pos, bf)) {
         linVar = (*pos)->get_rLinear().begin()->Var;
         set_obj_coefficient(linVar, get_obj_coefficient(linVar) + c);
         bf_toDelete.push_back(bf);
      } else {
         linVar = bf->linearize(*this);
         set_obj_coefficient(linVar, c);
         bf_list.insert(pos, bf);
      }
   }
   return linVar;
}
   
void ILP_Problem::add_bool_constraint(boolfunction* bf, bool_cons_type bct){
   std::list<boolfunction*> bfl;
   bfl.push_back(bf);
   add_bool_constraint(bfl, bct);
}

void ILP_Problem::add_bool_constraint(std::list<boolfunction*>& bfl, bool_cons_type bct){
   var v;
   row r;
   boolfunction* bf;
   switch(bct) {
      case C0:
         for(std::list<boolfunction*>::iterator it = bfl.begin(); it != bfl.end(); it++){ 
            v = add_boolfunction(*it, .0);
            add_basic_constraint(row(v) == 0, Static);
         }
         break;
      case C1:
         for(std::list<boolfunction*>::iterator it = bfl.begin(); it != bfl.end(); it++){ 
            v = add_boolfunction(*it, .0);
            add_basic_constraint(row(v) == 1, Static);
         }
         break;
      case CS:
      case CE:
         std::list<boolfunction*>::iterator it = bfl.begin();
         bf = (*it)->copy();
         r = 0;
         r += add_boolfunction(*it, .0);
         it++;
         for(;it != bfl.end(); it++){
            r += add_boolfunction(*it, .0);
            bf->add(OR, *it);
         }
         v = add_boolfunction(bf, .0);
         add_basic_constraint(r == v, Static);
         if(bct == CE)
            add_basic_constraint(row(v) == 1, Static);
         break;
   }
}

void ILP_Problem::extendSolution(solution& S) {
   for( std::list<monomial>::iterator it = monomials.begin(); it != monomials.end(); it++ ) {
      std::vector<var> a = it->get_A();
      //avoid monomials with only one variable
      if(a.size()<2)
         continue;
      S.set_value(it->get_L(), 1.);
      for( unsigned int i = 0; i < a.size(); i++ ) {
         if( S.value(a[i]) < 1e-6 ) {
            S.set_value(it->get_L(), 0.);
            break;
         }
      }
   }
   for( std::list<boolfunction*>::iterator bfit = bf_list.begin(); bfit != bf_list.end(); bfit++) {
      if ((*bfit)->get_op() == BASIC || (*bfit)->get_rLinear().size() > 1 || (*bfit)->get_rLinear().begin()->coeff != 1)
         continue;
      if ((*bfit)->evaluate(S))
         S.set_value((*bfit)->get_rLinear().begin()->Var,1.);
      else
         S.set_value((*bfit)->get_rLinear().begin()->Var,0.);
   }
}

void ILP_Problem::set_coefficient(var v, cons i, double d) {
  /*{\Mop sets the coefficient for variable $v$ and ineqality $i$
    to $d$ in the LP-Matrix.}*/
  coeff_list.push_back(Coefficient(v,i,d));
};

void ILP_Problem::set_obj_coefficient(var v, double d) {
   (v.Avar_pointer())->obj_ = d;
   (v.var_pointer())->obj_ = d;
   return;
}

double ILP_Problem::get_obj_coefficient(var v) const {
   return (v.Avar_pointer())->obj_;
}

void ILP_Problem::set_primal_heuristic(primal_heuristic* P) { // 
  /*{\Mop adds a primal heuristic to the model.}*/
  primal_heur=P;
};

void ILP_Problem::set_implicator(implicator* P) { // 
   implic = P;
}

void ILP_Problem::set_monomial_split_strategy(mon_split_strategy* spst) {
   mon_splitstrat=spst;
}

void ILP_Problem::set_bool_split_strategy(bool_split_strategy* spst) {
   bool_splitstrat=spst;
}

void ILP_Problem::set_qr_on_separate(bool qr){
   qr_on_separate = qr;
}

bool ILP_Problem::get_qr_on_separate() {
   return qr_on_separate;
}

void ILP_Problem::set_qrType(quadRefStatus qt) {
   qrType = qt;
}

quadRefStatus ILP_Problem::get_qrType(){
   return qrType;
}

#if 0
void ILP_Problem::set_branching(Branching* B) {
   bran = B;
   cerr << "Branching set to " << bran << endl;
   return;
}
#endif

void ILP_Problem::optimize() {
  /*{\Mop Computes the optimal solution for the current model.}*/

   add_sym_constraint(qr);
   //new monomials may be added by this method
   qr->reformulate_forced_constraints();

   if(monomials.size() > 0){
      //decompose all monomials to reduce polynomial instance to quadratic instance
      monomialinstance = new monomial_inst(monomials, *this, mon_splitstrat);
      add_sym_constraint(monomialinstance);
      qr->check_and_reformulate();
   }

   if(bf_list.size() > 0) {
      boolinstance = new bool_inst(bf_list, *this, bool_splitstrat);
      add_sym_constraint(boolinstance);
   }

   ABA_MASTER::optimize();
};


double ILP_Problem::get_solution(var v) {
  /*{\Mop returns the value of $v$ in the best found solution.}*/
  return Sol.value(v);
}

double ILP_Problem::get_solution(row& r) {
  row_entry::list_constiterator re;
  double d=0;
  foreach(re, r) d+=re->coeff*get_solution(re->Var);
  return d;
}

solution ILP_Problem::get_solution() {
   return Sol;
}

double ILP_Problem::get_optimum() {
   //FIXME
   return primalBound();
  /*
  var v;
  double d=0.;
  forall_variables(v, *myroot) d+=get_solution(v)*v.obj();
  return d;
  */
}

---