quadref.cc

#include <scil/constraints/quadref.h>
#include <boost/unordered_map.hpp>
#include <scil/core/monomial.h>
#include <sstream>

using namespace SCIL;
using namespace std;

//inserts a pair of variables with corresponding linearization variable into the hashmap
void QuadRef::hashPair(var& x, var& y, var& l){
   int index1, index2;
   if(x.var_pointer()->index() < y.var_pointer()->index()){
      index1 = x.var_pointer()->index();
      index2 = y.var_pointer()->index();
   } else {
      index1 = y.var_pointer()->index();
      index2 = x.var_pointer()->index();
   }
   std::stringstream ss;
   ss<<index1<<'#'<<index2;
   string newKey = ss.str();
   hashmap.insert(make_pair(newKey, l));
}

//adds a constraint to the internal constraint list
void QuadRef::addCons(cons_obj * c){
   cons_list.push_back(c);
}

//reformulates every constraint which has a reformulation flag given by the user
void QuadRef::reformulate_forced_constraints(){
   std::list<cons_obj*>::reverse_iterator rit;
   for(rit=cons_list.rbegin(); rit != cons_list.rend(); ++rit){
      cons_obj * cons = NULL;
      switch ((*rit)->get_qrStatus()){
         case SQK2:
            cons = reformulate_SQK2(*rit, true);
            break;
         case SQK3:
            reformulate_SQK3(*rit, true);
            break;
         case SQK2_SQK3:
            cons = reformulate_SQK2(*rit, true);
            reformulate_SQK3(*rit, true);
            break;
         default:
            break;
      }
      if(cons != NULL){
         IP.add_basic_constraint(cons, (*rit)->get_Act());
      }
   }
}

//checks every constraint in cons_list for reformulation without adding new variables
void QuadRef::check_and_reformulate(){
   std::list<cons_obj*>::reverse_iterator rit;
   for(rit=cons_list.rbegin(); rit != cons_list.rend(); ++rit){
      cons_obj * cons = NULL;
      if ((*rit)->get_qrStatus() == CHECK){
         switch (IP.get_qrType()){
            case SQK2:
               cons = reformulate_SQK2(*rit, false);
               break;
            case SQK3:
               reformulate_SQK3(*rit, false);
               break;
            case SQK2_SQK3:
               cons = reformulate_SQK2(*rit, false);
               reformulate_SQK3(*rit, false);
               break;
            default:
               break;
         }
         if (cons != NULL)
            IP.add_basic_constraint(cons, (*rit)->get_Act());
      }
   }
}

//constructs SQK2 constraint and returns it
cons_obj* QuadRef::reformulate_SQK2(cons_obj* c, bool force) {

   int index1, index2;
   string hashkey;
   boost::unordered_map<string, var>::iterator it;
   
   row r, q;
   c->non_zero_entries(r);
   r.normalize();
  
   //do not reformulate trivial constraints
   if(r.size() == 1 && (c->rhs() == 1 || c->rhs() == 0))
      return NULL;

   //construct quadratic constraint
   bool cancel = false;
   row_entry::list_iterator ri,rj;
   foreach(ri, r) {
      foreach(rj, r) {
         if(ri->Var == rj->Var){
            //linear case
            q += ri->Var * (ri->coeff * ri->coeff);
         } else {
            if((ri->Var).var_pointer()->index() < (rj->Var).var_pointer()->index()){
               index1 = (ri->Var).var_pointer()->index();
               index2 = (rj->Var).var_pointer()->index();
            } else {
               index1 = (rj->Var).var_pointer()->index();
               index2 = (ri->Var).var_pointer()->index();
            }
            std::stringstream ss;
            ss<<index1<<'#'<<index2;
            hashkey = ss.str();
            //search for already existing linearization variables
            it = hashmap.find(hashkey);
            if(force && it == hashmap.end()){
               //create new linearization variable
               q += IP.add_polynomial(ri->Var * rj->Var * ri->coeff * rj->coeff, false);
            } 
            else if(it != hashmap.end()){
               //use existing linearization variable
               q += (it->second) * (ri->coeff * rj->coeff);
            } 
            else{
               //if the reformulation is not forced and if no lin. variable is found cancel the reformulation
               cancel = true;
               break;
            }
         }       
      }
      if(cancel)
         break;
   }

   if(cancel)
      return NULL;
   
   double rhs = c->rhs() * c->rhs();

   q.normalize();

   //correction terms for coeffs < 0
   double d = 0.;
   foreach(ri, r) {
      if( ri->coeff < 0. )
         d += ri->coeff;
   }
   q += -2. * d * r;
   rhs += -2. * d * c->rhs();
   q.normalize();

   cons_obj* quad_cons;
   if( c->sense() == Less)
      quad_cons = (q<=rhs);
   else if( c->sense() == Greater )
      quad_cons = (q>=rhs);
   else
      quad_cons = (q==rhs);

   //set qrStatus for the reformulated constraint to NONE to avoid second reformulation
   quad_cons->set_qrStatus(NONE);

   return quad_cons;
}

//constructs SQK3 constraints and adds them to the internal list
void QuadRef::reformulate_SQK3(cons_obj* c, bool force) {
   int index1, index2;
   string hashkey;
   boost::unordered_map<string, var>::iterator it;

   row r, s;

   c->non_zero_entries(r);
   r.normalize();

   //do not reformulate trivial constraints
   if(r.size() == 1 && (c->rhs() == 1 || c->rhs() == 0))
      return;

   //construct quadratic constraints
   bool cancel;
   row_entry::list_iterator ri,rj;
   foreach(ri, r) {
      row t;
      cancel = false;
      foreach(rj, r) {
         if(ri->Var == rj->Var){
            //linear case
            t += ri->Var * (ri->coeff - c->rhs());
         }
         else {
            if((ri->Var).var_pointer()->index() < (rj->Var).var_pointer()->index()){
               index1 = (ri->Var).var_pointer()->index();
               index2 = (rj->Var).var_pointer()->index();
            } else {
               index1 = (rj->Var).var_pointer()->index();
               index2 = (ri->Var).var_pointer()->index();
            }
            std::stringstream ss;
            ss<<index1<<'#'<<index2;
            hashkey = ss.str();
            //search for already existing linearization variables
            it = hashmap.find(hashkey);
            if(force && it == hashmap.end()) {
               //create new linearization variable
               t += IP.add_polynomial(ri->Var * rj->Var * rj->coeff, false);
            }
            else if(it != hashmap.end()){
               //use existing linearization variable
               t += (it->second) * rj->coeff;
            }
            else{
               //if the reformulation is not forced and if no lin. variable is found cancel the reformulation
               cancel = true;
               break;
            }
         }
      }
      if(cancel)
         continue;
      if( ri->Var == r.begin()->Var) {
         s -= t;
         s += r;
         s.normalize();
         //add constraint to internal list (for later separation)
         if( c->sense() == Less)
            lc.push_back(s<=c->rhs());
         else if( c->sense() == Greater )
            lc.push_back(s>=c->rhs());
         else
            lc.push_back(s==c->rhs());
      }
      //add constraint to internal list (for later separation)
      if( c->sense() == Less )
         lc.push_back(t<=0.);
      else if( c->sense() == Greater )
         lc.push_back(t>=0.);
      else
         lc.push_back(t==0.);
   }
}

//separates SQK3 constraints
QuadRef::status QuadRef::standard_separation(subproblem& S) {
   std::list<cons_obj*>::iterator c;
   int num = 0;
   cons_obj* cons;
   foreach(c, lc) {
      if( (*c)->violation(S) > 1.0e-4 ) {
         row r;
         (*c)->non_zero_entries(r);
         if( (*c)->sense() == Less )
            cons = r<=(*c)->rhs();
         else if( (*c)->sense() == Greater )
            cons = r>=(*c)->rhs();
         else
            cons = r==(*c)->rhs();
         //sets qrStatus to NONE to avoid second reformulation
         cons->set_qrStatus(NONE);
         S.add_basic_constraint(cons);
         num++;
      }
   }
   if( num ) {
      return constraint_found;
   }

   return no_constraint_found;
}

---