cp_knapsack.h

#line 2 "cp_knapsack.lw"
#include <LEDA/map.h>
#include <LEDA/set.h>
#include <LEDA/basic.h>
#include <LEDA/d_array.h>
#include <LEDA/window.h>
#include <LEDA/file.h>
#include <LEDA/array.h>

class cmp_by_fractional: public leda_cmp_base<int>
{
 private:

map<int,double>* fractional_point ;

 public:

 cmp_by_fractional(map<int,double>* fract)
   {
     fractional_point = fract ;
   }
 int operator()(const int& a, const int& b) const

    {

      if ( (*fractional_point) [ a ] < (*fractional_point) [ b ] )
        { return 1 ;}
      if ( (*fractional_point) [ a ] > (*fractional_point) [ b ] )
        { return -1 ;}
      return 0 ; 
 
    }

} ;     

class cmp_by_coefficient: public leda_cmp_base<int>
{
 private:

array<int>* coeffs ;

 public:

 cmp_by_coefficient(array<int>* coeff)
   {
     coeffs = coeff ;
   }
 int operator()(const int& a, const int& b) const

   {

     if ( (*coeffs) [ a ] > (*coeffs) [ b ] )
       { return 1 ;}
     if ( (*coeffs) [ a ] < (*coeffs) [ b ] )
       { return -1 ;}
     return 0 ;
   }
} ;             

class cmp_decrease: public leda_cmp_base<double>
{

 public:

 cmp_decrease()
   {}
 int operator()(const double& a, const double& b) const

   {

     if ( a  > b  )
       { return -1 ;}
     if ( a  < b  )
       { return 1 ;}
     return 0 ;
   }
} ;           

class cmp_increase_by_coefficient: public leda_cmp_base<int>
{
 private:

array<int>* coeffs ;

 public:

 cmp_increase_by_coefficient(array<int>* coeff)
   {
     coeffs = coeff ;
   }
 int operator()(const int& a, const int& b) const

   {

     if ( (*coeffs) [ a ] > (*coeffs) [ b ] )
       { return 1 ;}
     if ( (*coeffs) [ a ] < (*coeffs) [ b ] )
       { return -1 ;}
     return 0 ;
   }
} ;             

class  cmp_nondecrease_by_reduced_costs: public leda_cmp_base<int>
{
 private:

map<int,double>* red_costs ;

 public:

  cmp_nondecrease_by_reduced_costs(map<int,double>* redcosts)
   {
     red_costs = redcosts ;
   }
 int operator()(const int& a, const int& b) const

   {

     if ( (*red_costs) [ a ] > (*red_costs) [ b ] )
       { return 1 ;}
     if ( (*red_costs) [ a ] < (*red_costs) [ b ] )
       { return -1 ;}
     return 0 ;
   }
} ;       



class LCI : public cons_obj {

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

  ~LCI()
    {};

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


class cp_knapsack : public sym_constraint {
 private:

  int my_count;
  cons_obj* KC ;
  map<int,var> item_index_map;
  cons_obj* new_constraint ;
  double violation_tolerance ;
  bool lifted_in_C2 ;
  bool TESTDONE ;
  bool lifted_in_F ;
  bool cover_found ;
  bool minimal_cover ;
  bool negative_coeffs ;
  bool test ;
  bool violated ;
  double contribute ; 
  int id ;
  int constraint_count ;
  int var1, var2 ;
  int index ; 
  int z ;
  int a ;
  int  c ;
  int sum ;
  int x ;
  double RHS ;
  double C1_size ;
  int rightsum ;
  int gammasum ;
  int Ti ;
  int number_variables ;
  int number_cover_candidates ; 
  int sort_array_size ;
  int index1 ;
  int index2 ;
  int coeff ;
  double  g ;
  int r ;
  int T ;
  int previous_T ;
  int lift_index ;
  int previous_lift_index ; 
  int lift_coeff ;
  int actual_lift_coeff ;
  int righthandside ;
  bool no_sep ;
  double fractional ;
  list<int> negative_coeff_indizes ;
  array<int>  liftcoeffs ;
  map<int,double> red_costs ;
  array<int>  coeffs ;
  array<int> cover_candidates ;
  map<int,double>  fractional_point ;
  //array<double>  red_costs ;
  int zippl ;
  array<bool> coeffs_sign ;
  array<int> sort_array ;
  array<int> liftsequence_on_R ;
  array<int> liftsequence_on_C2 ;
  array<int>* W1 ;
  array<int>* W2 ;
  set<int> F ;
  set<int> R ;
  set<int> C1 ;
  set<int> C2 ;
  set<int> to_lift_in_F ;
  set<int> to_lift_in_C2 ;
  set<int> to_lift_in_R ;
  set<int> N_C ;
  set<int> C ;
  set<int> L ;
  set<int> N ;
 
 public:
 
  cp_knapsack(cons_obj* KCt, double tolerance=0.001) : KC(KCt) {
     violation_tolerance = tolerance ;
  }

  void init(subproblem& S) {
    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!  
         coeff =(int) KC->coeff (((var) item_index_map[index]).var_pointer() )  ;  
         
         if ( coeff < 0 )
           {
             coeffs[index] = - coeff ; 
             righthandside -= coeff ;
           }
         else
           {
             coeffs[index] = coeff ;
           }      
       }
     
     //  insert indizes 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_indizes.push(index) ; 
               }
           }
         
       }
        

     if ( ( N.size() > 0 ) )
      {
        cover_candidates.resize( N.size() ) ;
        sort_array.resize( N.size() ) ;
        liftsequence_on_C2.resize( N.size() ) ;
        liftsequence_on_R.resize( N.size() ) ;
        no_sep = 0 ; 
      }

    else
      {
        no_sep = 1 ; 
      }
     
   }
  
 
 public: 

  // This function adds the value of the actually computed lifting coefficient    //
  // to the current value of T_i.                                                 //
  void update_Ti(int x)
  {
    Ti += liftcoeffs [ x ] ;    
  }

  // This function adds the value of the actually computed lifting coefficient    //
  // of a variable from the set C_2 to the actual value of gammasum.              //
  // gammasum represents the value of the sum \sum_{j \in C_2 \cap l_i} \gamma_j  //
  void update_gammasum(int x)
  {
    gammasum += liftcoeffs[x];            
  }
  // This function initializes the value rightsum with \sum_{j \in C_2} a_j       //
  // rightsum is used to compute the value b- \sum_{j \in L'_i} a_j               //
  int init_righthandsum()
  {
    rightsum = 0 ;
    forall (x,C2)
      {
        
        rightsum += coeffs [ x ] ;
      }
    return rightsum ;
  }
  // This function substracts the value of coeff[x] from the value rightsum       //
  // and returns the value                                                        //
  int  update_righthandsum(int x)
  {
    rightsum -= coeffs[x] ;
    return rightsum ;
  }  


bool compute_cover()
{

  // COVER GENEARTION //                                                  

  // write all indices of variables with fractional value > 0 in       //
  // array cover_candidates                                            //

  index1 = 0 ;  
  forall(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 ;
    }
  number_cover_candidates = index1 ; // save size of array 'cover_candidates'

  // sort array by nonincreasing values of fractional point            //
  cmp_by_fractional nonincrease_by_fractional( &fractional_point ) ;
  cover_candidates.sort(nonincrease_by_fractional,0 , number_cover_candidates-1) ;

  // try to generate a cover, by summing up  the coefficients of the      //
  // variables in order of nonincreasing values of fractional point       //  
  // The resulting cover indizes are stored in array 'sort_array'         //
  
  cover_found = 0 ;
  sum = 0 ;
  index1 = 0 ;

  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 //                                                  
  sort_array_size = index1-1 ; // save size of array 'sort_candidates'

  // sort cover indizes array by nondecreasing values of coefficients  //
  cmp_by_coefficient nondecrease_by_coeff( &coeffs ) ;
  sort_array.sort(nondecrease_by_coeff,0,sort_array_size) ;
  
  // try to eleminate the variables of the Cover in order of  //
  // nondecreasing value of coefficients                      //
  minimal_cover = 0 ;
  index1 = 0 ;
  
  while( (!minimal_cover) && (index1 < sort_array_size-1 ) )
    {
      sum -= coeffs [ (sort_array[ index1 ]) ] ;
 
     if ( sum > righthandside )
        {
          C.del( sort_array [ index1 ] ) ; 
          index1 ++ ;
        }
      else
        {
          minimal_cover = 1 ;
        }
    }
  
  // COVER PARTITION // 

  index1 = 0 ; 
  sum = 0 ;
  
  forall (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 += C2 ;
      C2.clear() ;
    }
  
  if ( C1.size() <= 1 ) // STOP, cardinality C_1 is too small to generate LCI
    {
      return 1 ;
    }
 
  N_C = N - ( C1 + C2 ) ;
  
  forall(x , N_C )
    {
      if ( fractional_point[ x ] >  0.00001 )  
        {
          F.insert(x) ;
      }
      else
        {
          R.insert(x) ;
        }
    }
  
  
  // LIFTINGORDER     // 
  
  cmp_nondecrease_by_reduced_costs reduced( &red_costs ) ;
  
  // lift variables in R in order of          //
  // nonincreasing value of fractional point   //
  index1 = 0;

  if ( C2.size() != 0 )
    {    
      forall( x, C2 )
        {
          liftsequence_on_C2[ index1 ] = x ;
         index1++ ;
        }
      liftsequence_on_C2.sort( reduced,0,index1-1) ;  
      //liftsequence_on_C2.sort( nonincrease_by_fractional,0,index1-1) ;
    }
 
  
 index1 = 0 ;
 
 if ( R.size() != 0 )
   {      
     forall( x, R )
       {
         liftsequence_on_R[ index1 ] = x ;
         index1 ++ ;
       }
     liftsequence_on_R.sort( reduced,0,index1-1) ;
     //liftsequence_on_R.sort( nonincrease_by_fractional,0,index1-1) ;
   }
 

   return 0 ;

}
  




  void lift_on_F( array<int>* W )
  {
    contribute = -1 ;
    lift_index = -1 ;
    actual_lift_coeff = -1 ;
    lift_coeff = -1 ;
    
    forall( x, to_lift_in_F )
      {
       RHS = righthandside - (rightsum+coeffs[x]) ; ;
     
       if ( RHS >= 0 )
         {        
           z = (*W).binary_locate( 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.del( lift_index ) ;
        L.insert( lift_index ) ;
        update_Ti(lift_index) ;
        previous_lift_index = lift_index ;
        lifted_in_F = 1 ;
      }
    
  }

  
  void lift_on_C2 ( array<int>* W )
  {
  
    x = C2.size() - to_lift_in_C2.size() ;
    lift_index = liftsequence_on_C2[ x ] ;
    RHS = righthandside - update_righthandsum( lift_index ) ;
    z = (*W).binary_locate( RHS ) ;
    lift_coeff = z - r +1 - gammasum ; 
    liftcoeffs[ lift_index ] = lift_coeff ;
    update_gammasum(lift_index) ;
    to_lift_in_C2.del( lift_index ) ;
    L.insert( lift_index ) ;
    update_Ti(lift_index) ;
    previous_lift_index = lift_index ;
    lifted_in_C2 = 1 ;
    
 
  }
  
  
  void lift_on_R ( array<int>* W )
  {
   
    x = R.size() - to_lift_in_R.size() ;
    lift_index = liftsequence_on_R[ x ] ;
    RHS = righthandside - (rightsum+coeffs[lift_index]) ;
    z = (*W).binary_locate( RHS )   ;
    lift_coeff = r - 1 + gammasum - z ; 
    liftcoeffs[ lift_index ] = lift_coeff ;
    to_lift_in_R.del( lift_index ) ;
    L.insert( lift_index ) ;
    update_Ti(lift_index) ;
    previous_lift_index = lift_index ;
    
  }
  
  void compute_first_row()
    
  {
    r = C1.size() ;
    index1 = 0;

    forall( x, C1 )
      {
        sort_array[ index1 ] = coeffs[ x ] ;
        index1 ++ ; 
      }
    sort_array.sort(0,index1-1) ;



    W1 = new array<int>(0,r) ;
    (*W1)[0] = 0 ;
    
   for (int t=0; t <= r-1 ;  t++ )
     {     
       (*W1)[ t+1 ] = (*W1)[ t ] + sort_array[ t ] ;
     }
   
   W2 = W1 ;
  }
  
  void compute_next_row()
  {
    T = Ti ;
    W2 = new  array<int>(0,T) ;
    c = liftcoeffs[ previous_lift_index ] ;
    a = coeffs [  previous_lift_index ] ; 
    for( int index=0; index <= (*W1).size()-1 ; index++)
      {  
        var1 = (*W1)[ index ] ;    
        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 setup_values(subproblem& S)
  {
    
    // write array containing the current fractional solution according the substitution 
   for (int index=0; index <= number_variables-1; index++ )
     {
         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 ] ) ; 
           }
         
          
          
     }
   red_costs.clear() ;
   forall(x,N)
     {
       red_costs[x] = S.red_cost( item_index_map[ x ] ) ; 
     }
   
   test = 0 ; 
   F.clear() ;
   R.clear() ;
   C1.clear() ;
   C2.clear() ;
   to_lift_in_F.clear() ;
   to_lift_in_C2.clear() ;
   to_lift_in_R.clear() ;
   N_C.clear() ;
   C.clear() ;
   liftcoeffs.init(0) ;
   
  }


status separate(subproblem& S) 
{

  if ( no_sep == 1)
    {
      return no_constraint_found ;
    }

  setup_values(S) ;

  
  if ( compute_cover() == 1  ) 
    {
      return no_constraint_found ;
    }
  
  forall (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 ;
  
  TESTDONE = 0 ;
  
  
  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 ) )
       {
         
         C1_size = C1.size()-1 ;

          
         
         g = 0 ;
         
         forall (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 ] )  ;
                 
               } 
           }
         
         forall (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 ;
    
    } 
  
  
  RHS = C1.size()-1 + gammasum  ;
  
  forall ( index, negative_coeff_indizes)
    {         
      liftcoeffs [index] = -(liftcoeffs [index]) ;
      RHS += liftcoeffs[index] ;
    }
  
  new_constraint = new LCI(liftcoeffs,RHS,item_index_map) ;
  S.add_basic_constraint( new_constraint, Dynamic ) ; 
  my_count ++ ;
  cout<<"cp knapsack constraint\n";
  return constraint_found ; // one LCI added 
 
 
   
}
};

---