cp_knapsack.cc

#include <algorithm>
#include <boost/foreach.hpp>
#include <scil/constraints/cp_knapsack.h>

namespace SCIL {

LCI::LCI(std::vector<int>coeff,int r, std::map<int,var>&  map)
   : cons_obj(Less, r), item_index_map(map) {
   coeff_array=coeff ;
   rhs = r ;
}

LCI::~LCI() {};

/*  
   double LCI::coeff(var_obj* v) {
   int index ;
   index = item_index_map[ v ] ;
   return coeff_array[ index ] ;
}
*/
void LCI::non_zero_entries(row& r) {   
   for(int i=0; i < int(coeff_array.size()); i++) {
      if(coeff_array[i]!=0)
         r+=item_index_map[i]*coeff_array[i];
   }
}

int CP_KNAPSACK::find_largest_index(std::vector<int> vec, int val) {
   std::vector<int>::iterator index;
   index = std::find(vec.begin(), vec.end(), val);
   //Check whether val was found in vec
   if (index != vec.end())
   {
      //Increase index to index of last element with a value equal to val
      while ( *(index+1) == *index ) index++;
      return index - vec.begin();
   }
   return -1;
}

CP_KNAPSACK::CP_KNAPSACK(cons_obj* KCt, double tolerance) : KC(KCt) {
   violation_tolerance = tolerance ;
}

void CP_KNAPSACK::init(subproblem& S) {
   if (S.configuration("CP_KNAPSACK_Debug_Out")=="true")
      cerr << "CP_KNAPSACK::init\n";
   S.add_basic_constraint(KC);

   number_variables = S.nVar() ; 
   righthandside = (int) KC->rhs() ;
   coeffs.resize(number_variables);
   liftcoeffs.resize(number_variables);
  
   // substitute negative coeffs  //
  
   for (int index=0; index <= number_variables-1; index++ )
   {
      item_index_map[index]=S.get_var(index); // TODO von mir!  
      int coeff =(int) KC->coeff( ((var) item_index_map[index]).var_pointer() )  ;  
       
      if ( coeff < 0 ) {
         coeffs[index] = - coeff ; 
         righthandside -= coeff ;
      }
      else {
         coeffs[index] = coeff ;
      }   
   }
  
   //  insert indices of coeffs a_j, with   //
   // (a_j <= righthandside) && ( a_j != 0) //
   // in set N.                             //
  
  
   for (int index=0; index <= number_variables-1; index++ ) {
      if ( (coeffs[index] <=  righthandside) && ( coeffs[index] != 0 ) )
            {
         N.insert( index ) ;
          
         if ( KC->coeff( ((var) item_index_map[index]).var_pointer() ) < 0 ) 
         {
            negative_coeff_indices.push_front(index) ; 
         }
      }
      
   }
  
  
  if ( ( N.size() > 0 ) ) no_sep = 0 ; 
  else no_sep = 1 ; 
  
}

void CP_KNAPSACK::update_Ti(int x)
{
   Ti += liftcoeffs [ x ] ;    
}

void CP_KNAPSACK::update_gammasum(int x)
{
   gammasum += liftcoeffs[x];            
}

int CP_KNAPSACK::init_righthandsum()
{
   rightsum = 0 ;
   BOOST_FOREACH(int x,C2)
   {
   
      rightsum += coeffs [ x ] ;
   }
   return rightsum ;
}

int CP_KNAPSACK::update_righthandsum(int x)
{
   rightsum -= coeffs[x] ;
   return rightsum ;
}  

//Comparison function to sort input by nonincreasing values
bool nonincreasing (double i, double j) { return (i>j); }

bool CP_KNAPSACK::compute_cover()
{
  
   // COVER GENEARTION //                                                  

   std::vector<int> cover_candidates; // Indices of cover candidates
  
   // write all indices of variables with fractional value > 0 in       //
   // array cover_candidates                                            //
  
   int index1 = 0 ;  
   BOOST_FOREACH(index,N) // check all variables 
   { 
      if ( (fractional_point[index] > 0.00001) ) 
      {
         cover_candidates [ index1 ] = index ;
         index1++; 
      }
   }
  
   if ( index1 == 0 )   // STOP, all variables equal to zero
   {    
      return 1 ;
   }
   int number_cover_candidates = index1; //save size of array 'cover_candidates'

   // sort array by nonincreasing values of fractional point            //
   std::sort(cover_candidates.begin(), cover_candidates.end(),nonincreasing);

   // try to generate a cover, by summing up  the coefficients of the      //
   // variables in order of nonincreasing values of fractional points      //  
   // The resulting cover indices are stored in array 'sort_array'         //
  
   bool cover_found = 0 ; //true if cover was found
   int sum = 0 ;
   index1 = 0 ;
   std::set<int> C ;

   while ( (cover_found == 0) && ( index1 < number_cover_candidates ) )
   {
      sum += coeffs [ cover_candidates[ index1 ] ];
      C.insert( cover_candidates[ index1 ] ) ;
      sort_array[ index1 ] = cover_candidates[ index1 ] ;
      index1 ++ ;
      
      if ( sum > righthandside ) 
      {
         cover_found = 1 ;
      }
      
   }

   if ( cover_found == 0 ) // STOP, no cover found  
   {           
      return 1 ;
   }
 
 
   // MINIMAL COVER CONVERSION //                                                  
   int sort_array_size = index1-1 ; // save size of array 'sort_candidates'

   // sort cover indices array by nondecreasing values of coefficients  //
   std::sort(sort_array.begin(), sort_array.end());
  
   // try to eleminate the variables of the Cover in order of  //
   // nondecreasing value of coefficients                      //
   bool minimal_cover = 0 ; //true if minimal cover was found
   index1 = 0 ;
  
   while( (!minimal_cover) && (index1 < sort_array_size-1 ) )
   {
      sum -= coeffs [ (sort_array[ index1 ]) ] ;
 
      if ( sum > righthandside )
            {
         C.erase( sort_array [ index1 ] ) ; 
         index1 ++ ;
      }
      else
      {
         minimal_cover = 1 ;
      }
   }
  
   // COVER PARTITION // 

   index1 = 0 ; 
   sum = 0 ;
  
   BOOST_FOREACH (int x,C)
   {
      if ( ( fractional_point[ x ] == 1 ) )
      {
         sum += coeffs[x] ;
         C2.insert( x ) ; 
      } 
      else
      {
         sum += coeffs[x] ;
         C1.insert( x ) ;
         index1 ++ ;
      }
      
   }
  
  
   //      //

   if ( ( C1.size() <= 1 ) && ( C2.size() > 1 ) ) 
   {
      
      C1.insert(C2.begin(), C2.end()) ;
      C2.clear() ;
   }
  
   if ( C1.size() <= 1 ) // STOP, cardinality C_1 is too small to generate LCI
   {
      return 1 ;
   }

   //Calculates N - (C1 + C2) and stores result in N_C 
   std::set<int> N_C ;
   N_C.insert(C1.begin(), C1.end());
   N_C.insert(C2.begin(), C2.end());
   std::vector<int> difference;
   std::set_difference(N.begin(), N.end(), N_C.begin(), N_C.end(), 
                       difference.begin());
   N_C.insert(difference.begin(), difference.end());
  
   BOOST_FOREACH(int x , N_C )
   {
      if ( fractional_point[ x ] >  0.00001 )  
      { 
         F.insert(x) ;
      }
      else
            {
         R.insert(x) ;
      }
   }
  
  
   // LIFTINGORDER     // 
  
   // lift variables in R in order of           //
   // nonincreasing value of fractional point   //
   std::vector<int>::iterator indx1 = liftsequence_on_C2.begin();

   if ( C2.size() != 0 )
   {    
      BOOST_FOREACH(int x, C2)
      {
         *indx1 = x ;
         indx1++ ;
      }
      std::sort(liftsequence_on_C2.begin(), indx1-1);
   }
 
  
   indx1 = liftsequence_on_R.begin() ;
 
   if ( R.size() != 0 )
   {      
      BOOST_FOREACH(int x, R)
      {
         *indx1 = x ;
         indx1 ++ ;
      }
      std::sort(liftsequence_on_R.begin(), indx1-1);
   }

   return 0 ;
}

void CP_KNAPSACK::lift_on_F( std::vector<int>* W )
{
   double contribute = -1 ;
   int  lift_index = -1 ;
   int actual_lift_coeff = -1 ;
   int lift_coeff = -1 ;
  
   BOOST_FOREACH(int x, to_lift_in_F)
   {
      int RHS = righthandside - (rightsum+coeffs[x]) ;
      
      if ( RHS >= 0 )
      {           
         int z = find_largest_index(*W, RHS);
         actual_lift_coeff = r - 1 + gammasum - z ;      
         if ( actual_lift_coeff * fractional_point[ x ] > contribute )
         {
            lift_coeff = actual_lift_coeff ;
            contribute = lift_coeff * fractional_point[ x ] ; 
            lift_index = x ; 
         }
          
      }
      
   }
  
   if ( lift_index != -1 )
   {
      
      liftcoeffs[ lift_index ] = lift_coeff ;
      to_lift_in_F.erase( lift_index ) ;
      update_Ti(lift_index) ;
      previous_lift_index = lift_index ;
      lifted_in_F = 1 ;
   }
  
}

void CP_KNAPSACK::lift_on_C2 (std:: vector<int>* W )
{
  
   int x = C2.size() - to_lift_in_C2.size() ;
   int lift_index = liftsequence_on_C2[ x ] ;
   int RHS = righthandside - update_righthandsum( lift_index ) ;
   int z = find_largest_index(*W, RHS);
   int lift_coeff = z - r + 1 - gammasum ; 
   liftcoeffs[ lift_index ] = lift_coeff ;
   update_gammasum(lift_index) ;
   to_lift_in_C2.erase( lift_index ) ;
   update_Ti(lift_index) ;
   previous_lift_index = lift_index ;
   lifted_in_C2 = 1 ;
  
}
  
void CP_KNAPSACK::lift_on_R ( std::vector<int>* W )
{
  
   int x = R.size() - to_lift_in_R.size() ;
   int lift_index = liftsequence_on_R[ x ] ;
   int RHS = righthandside - (rightsum+coeffs[lift_index]) ;
   int z = find_largest_index(*W, RHS);
   int lift_coeff = r - 1 + gammasum - z ; 
   liftcoeffs[ lift_index ] = lift_coeff ;
   to_lift_in_R.erase( lift_index ) ;
   update_Ti(lift_index) ;
   previous_lift_index = lift_index ;
}

void CP_KNAPSACK::compute_first_row()
{
 
   r = C1.size() ;
   std::vector<int>::iterator indx1 = sort_array.begin();
  
   BOOST_FOREACH(int x, C1)
   {
      *indx1 = coeffs[ x ] ;
      indx1++ ; 
   }
   std::sort(sort_array.begin(), indx1-1);
  
   W1 = new std::vector<int>() ;
   (*W1)[0] = 0 ;

   std::vector<int>::iterator it;
   for( it = W1->begin(); it != W1->end(); it++)
      *(it+1) = *it + sort_array[it - W1->begin()];

   W2 = W1 ;
}
  
void CP_KNAPSACK::compute_next_row()
{
   int T = Ti ;
   W2 = new std::vector<int>();
   int c = liftcoeffs[ previous_lift_index ] ;
   int a = coeffs [  previous_lift_index ] ; 
   int var1 = (*W1)[ index ] ;
   for( int index=0; index <= int((*W1).size()-1) ; index++)
   {  
      int var2;
      if ( index-c > 0 )
      {
         var2 = (*W1)[ index - c] ;
      }
      else 
      {
         var2 = 0 ;
      }
      
      var2 = var2 + a ;
      
      if ( var1 < var2 )
            {
         (*W2)[ index ] = var1 ;
      }
      else
      {
         (*W2)[ index ] = var2 ;
          
      }
      
   }
  
   for( int index = (*W1).size() ; index<=T; index++)
   {
      if ( index-c > 0 )
      {
         var1 = (*W1)[ index - c] ;
          
      }
      else 
      {
         var1 = 0 ;
      }
      
      (*W2)[ index ] = var1 + a ;
   }
   delete W1;
   W1 = W2 ;
}

void CP_KNAPSACK::setup_values(subproblem& S)
{
   if (S.configuration("CP_KNAPSACK_Debug_Out")=="true")
      cerr << "CP_KNAPSACK::setup_values\n";
   // write array containing the current fractional solution according to the substitution 
   for (int index=0; index <= number_variables-1; index++ )
   {
      int coeff = (int) KC->coeff(((var) item_index_map[index]).var_pointer() ) ;
      
      if ( coeff >= 0 )
      {
         fractional_point[index] = S.value( item_index_map[ index ] ) ; 
      }
      else
      {
         fractional_point[index] = 1 - S.value( item_index_map[ index ] ) ; 
      }
      
      
      
   }
   F.clear() ;
   R.clear() ;
   C1.clear() ;
   C2.clear() ;
   to_lift_in_F.clear() ;
   to_lift_in_C2.clear() ;
   to_lift_in_R.clear() ;
   BOOST_FOREACH(int &member, liftcoeffs)
      member = 0;   
}

CP_KNAPSACK::status CP_KNAPSACK::standard_separation(subproblem& S) 
{
   if (S.configuration("CP_KNAPSACK_Debug_Out")=="true")
      cerr << "CP_KNAPSACK::seperate\n";  
   if ( no_sep == 1)
   {
      return no_constraint_found ;
   }
  
   setup_values(S) ;

  
   if ( compute_cover() == 1  ) 
   {
      return no_constraint_found ;
   }
  
   BOOST_FOREACH (int x,C1)
   {
      liftcoeffs[ x ] = 1 ;   
   }

   init_righthandsum();
   Ti=C1.size() ;
   gammasum = 0 ;
  
   compute_first_row() ;
  
  
   to_lift_in_F = F ;
   to_lift_in_C2 = C2 ;
   to_lift_in_R = R ;
  
   bool TESTDONE = 0 ; //True if test is done
  
  
   while ( (to_lift_in_R.size() != 0 ) || (to_lift_in_C2.size() != 0) 
            || ( to_lift_in_F.size() != 0) )
   {
      lifted_in_F = 0 ;
      lifted_in_C2 = 0 ;
      

      if ( to_lift_in_F.size() != 0 )

      {   
         lift_on_F(W2) ;
      }
 

      if ( ( to_lift_in_F.size() == 0 ) && ( !TESTDONE ) )
      {
         
         double C1_size = C1.size()-1 ;

          
         
         double g = 0 ;
         
         BOOST_FOREACH (int x, F)
         {
            if (  KC->coeff ( ((var) item_index_map[x]).var_pointer() )  < 0 )
            {
                 
               g -= liftcoeffs[x]* S.value( item_index_map[ x ] )  ;
               C1_size -=  liftcoeffs[x] ;
            }
            else
            {
               g += liftcoeffs[x]* S.value( item_index_map[ x ] )  ;
                 
            } 
         }
         
         BOOST_FOREACH (int x, C1)
         {  
            if ( KC->coeff ( ((var) item_index_map[x]).var_pointer() ) < 0 )
            {
               g -= S.value( item_index_map[ x ] )  ;
               C1_size -=  liftcoeffs[x] ;
            }
            else
            {
               g += S.value( item_index_map[ x ] )   ;
               
            }
         }
         
         if ( g-violation_tolerance <= C1_size )
         {           
            return no_constraint_found ; 
         }
         
         TESTDONE = 1 ;   
      }

   

      if ( (!lifted_in_F ) && ( to_lift_in_C2.size() != 0 ) )

      {
         lift_on_C2(W2) ;
      }
 
      if  ( (!lifted_in_C2) && (!lifted_in_F) && ( to_lift_in_C2.size() == 0) 
             && ( to_lift_in_R.size() != 0 ) && ( to_lift_in_F.size() == 0 ) )  
      {
         lift_on_R(W2) ;
      }

      if ( (to_lift_in_R.size() != 0 ) || (to_lift_in_C2.size() != 0) 
            || ( to_lift_in_F.size() != 0) )
      {
         compute_next_row() ;
      }

      //delete W2 ;
    
   } 
  
  
   int RHS = C1.size()-1 + gammasum  ;
  
   BOOST_FOREACH ( index, negative_coeff_indices)
   {         
      liftcoeffs [index] = -(liftcoeffs [index]) ;
      RHS += liftcoeffs[index] ;
   }
  
   cons_obj *new_constraint = new LCI(liftcoeffs,RHS,item_index_map) ;
   S.add_basic_constraint( new_constraint, Dynamic ) ; 
   cout<<"cp knapsack constraint\n";
   return constraint_found ; // one LCI added 
}
void CP_KNAPSACK::info() {
   cout<<"Configurations:\n";
   cout<<"CP_KNAPSACK_Debug_Out [true|false]\n";
};

}

---