cons_knapsack.c

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*                  This file is part of the program and library             */
/*         SCIP --- Solving Constraint Integer Programs                      */
/*                                                                           */
/*  Copyright 2002-2022 Zuse Institute Berlin                                */
/*                                                                           */
/*  Licensed under the Apache License, Version 2.0 (the "License");          */
/*  you may not use this file except in compliance with the License.         */
/*  You may obtain a copy of the License at                                  */
/*                                                                           */
/*      http://www.apache.org/licenses/LICENSE-2.0                           */
/*                                                                           */
/*  Unless required by applicable law or agreed to in writing, software      */
/*  distributed under the License is distributed on an "AS IS" BASIS,        */
/*  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. */
/*  See the License for the specific language governing permissions and      */
/*  limitations under the License.                                           */
/*                                                                           */
/*  You should have received a copy of the Apache-2.0 license                */
/*  along with SCIP; see the file LICENSE. If not visit scipopt.org.         */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

/**@file   cons_knapsack.c
 * @ingroup DEFPLUGINS_CONS
 * @brief  Constraint handler for knapsack constraints of the form  \f$a^T x \le b\f$, x binary and \f$a \ge 0\f$.
 * @author Tobias Achterberg
 * @author Xin Liu
 * @author Kati Wolter
 * @author Michael Winkler
 * @author Tobias Fischer
 */

/*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/

#include "blockmemshell/memory.h"
#include "scip/cons_knapsack.h"
#include "scip/cons_linear.h"
#include "scip/cons_logicor.h"
#include "scip/cons_setppc.h"
#include "scip/pub_cons.h"
#include "scip/pub_event.h"
#include "scip/pub_implics.h"
#include "scip/pub_lp.h"
#include "scip/pub_message.h"
#include "scip/pub_misc.h"
#include "scip/pub_misc_select.h"
#include "scip/pub_misc_sort.h"
#include "scip/pub_sepa.h"
#include "scip/pub_var.h"
#include "scip/scip_branch.h"
#include "scip/scip_conflict.h"
#include "scip/scip_cons.h"
#include "scip/scip_copy.h"
#include "scip/scip_cut.h"
#include "scip/scip_event.h"
#include "scip/scip_general.h"
#include "scip/scip_lp.h"
#include "scip/scip_mem.h"
#include "scip/scip_message.h"
#include "scip/scip_nlp.h"
#include "scip/scip_numerics.h"
#include "scip/scip_param.h"
#include "scip/scip_prob.h"
#include "scip/scip_probing.h"
#include "scip/scip_sol.h"
#include "scip/scip_solvingstats.h"
#include "scip/scip_tree.h"
#include "scip/scip_var.h"
#include <ctype.h>
#include <string.h>

#ifdef WITH_CARDINALITY_UPGRADE
#include "scip/cons_cardinality.h"
#endif

/* constraint handler properties */
#define CONSHDLR_NAME          "knapsack"
#define CONSHDLR_DESC          "knapsack constraint of the form  a^T x <= b, x binary and a >= 0"
#define CONSHDLR_SEPAPRIORITY   +600000 /**< priority of the constraint handler for separation */
#define CONSHDLR_ENFOPRIORITY   -600000 /**< priority of the constraint handler for constraint enforcing */
#define CONSHDLR_CHECKPRIORITY  -600000 /**< priority of the constraint handler for checking feasibility */
#define CONSHDLR_SEPAFREQ             0 /**< frequency for separating cuts; zero means to separate only in the root node */
#define CONSHDLR_PROPFREQ             1 /**< frequency for propagating domains; zero means only preprocessing propagation */
#define CONSHDLR_EAGERFREQ          100 /**< frequency for using all instead of only the useful constraints in separation,
                                         *   propagation and enforcement, -1 for no eager evaluations, 0 for first only */
#define CONSHDLR_MAXPREROUNDS        -1 /**< maximal number of presolving rounds the constraint handler participates in (-1: no limit) */
#define CONSHDLR_DELAYSEPA        FALSE /**< should separation method be delayed, if other separators found cuts? */
#define CONSHDLR_DELAYPROP        FALSE /**< should propagation method be delayed, if other propagators found reductions? */
#define CONSHDLR_NEEDSCONS         TRUE /**< should the constraint handler be skipped, if no constraints are available? */

#define CONSHDLR_PRESOLTIMING            SCIP_PRESOLTIMING_ALWAYS
#define CONSHDLR_PROP_TIMING             SCIP_PROPTIMING_BEFORELP

#define EVENTHDLR_NAME         "knapsack"
#define EVENTHDLR_DESC         "bound change event handler for knapsack constraints"
#define EVENTTYPE_KNAPSACK SCIP_EVENTTYPE_LBCHANGED \
                         | SCIP_EVENTTYPE_UBTIGHTENED \
                         | SCIP_EVENTTYPE_VARFIXED \
                         | SCIP_EVENTTYPE_VARDELETED \
                         | SCIP_EVENTTYPE_IMPLADDED /**< variable events that should be caught by the event handler */

#define LINCONSUPGD_PRIORITY    +100000 /**< priority of the constraint handler for upgrading of linear constraints */

#define MAX_USECLIQUES_SIZE        1000 /**< maximal number of items in knapsack where clique information is used */
#define MAX_ZEROITEMS_SIZE        10000 /**< maximal number of items to store in the zero list in preprocessing */

#define KNAPSACKRELAX_MAXDELTA        0.1 /**< maximal allowed rounding distance for scaling in knapsack relaxation */
#define KNAPSACKRELAX_MAXDNOM      1000LL /**< maximal allowed denominator in knapsack rational relaxation */
#define KNAPSACKRELAX_MAXSCALE     1000.0 /**< maximal allowed scaling factor in knapsack rational relaxation */

#define DEFAULT_SEPACARDFREQ          1 /**< multiplier on separation frequency, how often knapsack cuts are separated */
#define DEFAULT_MAXROUNDS             5 /**< maximal number of separation rounds per node (-1: unlimited) */
#define DEFAULT_MAXROUNDSROOT        -1 /**< maximal number of separation rounds in the root node (-1: unlimited) */
#define DEFAULT_MAXSEPACUTS          50 /**< maximal number of cuts separated per separation round */
#define DEFAULT_MAXSEPACUTSROOT     200 /**< maximal number of cuts separated per separation round in the root node */
#define DEFAULT_MAXCARDBOUNDDIST    0.0 /**< maximal relative distance from current node's dual bound to primal bound compared
                                         *   to best node's dual bound for separating knapsack cuts */
#define DEFAULT_DISAGGREGATION     TRUE /**< should disaggregation of knapsack constraints be allowed in preprocessing? */
#define DEFAULT_SIMPLIFYINEQUALITIES TRUE/**< should presolving try to simplify knapsacks */
#define DEFAULT_NEGATEDCLIQUE      TRUE /**< should negated clique information be used in solving process */

#define MAXABSVBCOEF               1e+5 /**< maximal absolute coefficient in variable bounds used for knapsack relaxation */
#define USESUPADDLIFT             FALSE /**< should lifted minimal cover inequalities using superadditive up-lifting be separated in addition */

#define DEFAULT_PRESOLUSEHASHING   TRUE /**< should hash table be used for detecting redundant constraints in advance */
#define HASHSIZE_KNAPSACKCONS       500 /**< minimal size of hash table in linear constraint tables */

#define DEFAULT_PRESOLPAIRWISE     TRUE /**< should pairwise constraint comparison be performed in presolving? */
#define NMINCOMPARISONS          200000 /**< number for minimal pairwise presolving comparisons */
#define MINGAINPERNMINCOMPARISONS 1e-06 /**< minimal gain per minimal pairwise presolving comparisons to repeat pairwise 
                                         *   comparison round */
#define DEFAULT_DUALPRESOLVING     TRUE /**< should dual presolving steps be performed? */
#define DEFAULT_DETECTCUTOFFBOUND  TRUE /**< should presolving try to detect constraints parallel to the objective
                                         *   function defining an upper bound and prevent these constraints from
                                         *   entering the LP */
#define DEFAULT_DETECTLOWERBOUND TRUE   /**< should presolving try to detect constraints parallel to the objective
                                         *   function defining a lower bound and prevent these constraints from
                                         *   entering the LP */
#define DEFAULT_CLIQUEEXTRACTFACTOR 0.5 /**< lower clique size limit for greedy clique extraction algorithm (relative to largest clique) */
#define MAXCOVERSIZEITERLEWI       1000 /**< maximal size for which LEWI are iteratively separated by reducing the feasible set */

#define DEFAULT_USEGUBS           FALSE /**< should GUB information be used for separation? */
#define GUBCONSGROWVALUE              6 /**< memory growing value for GUB constraint array */
#define GUBSPLITGNC1GUBS          FALSE /**< should GNC1 GUB conss without F vars be split into GOC1 and GR GUB conss? */
#define DEFAULT_CLQPARTUPDATEFAC   1.5  /**< factor on the growth of global cliques to decide when to update a previous
                                         *   (negated) clique partition (used only if updatecliquepartitions is set to TRUE) */
#define DEFAULT_UPDATECLIQUEPARTITIONS FALSE /**< should clique partition information be updated when old partition seems outdated? */
#define MAXNCLIQUEVARSCOMP 1000000      /**< limit on number of pairwise comparisons in clique partitioning algorithm */
#ifdef WITH_CARDINALITY_UPGRADE
#define DEFAULT_UPGDCARDINALITY   FALSE /**< if TRUE then try to update knapsack constraints to cardinality constraints */
#endif

/* @todo maybe use event SCIP_EVENTTYPE_VARUNLOCKED to decide for another dual-presolving run on a constraint */

/*
 * Data structures
 */

/** constraint handler data */
struct SCIP_ConshdlrData
{
   int*                  ints1;              /**< cleared memory array, all entries are set to zero in initpre, if you use this
                                              *   you have to clear it at the end, exists only in presolving stage */
   int*                  ints2;              /**< cleared memory array, all entries are set to zero in initpre, if you use this
                                              *   you have to clear it at the end, exists only in presolving stage */
   SCIP_Longint*         longints1;          /**< cleared memory array, all entries are set to zero in initpre, if you use this
                                              *   you have to clear it at the end, exists only in presolving stage */
   SCIP_Longint*         longints2;          /**< cleared memory array, all entries are set to zero in initpre, if you use this
                                              *   you have to clear it at the end, exists only in presolving stage */
   SCIP_Bool*            bools1;             /**< cleared memory array, all entries are set to zero in initpre, if you use this
                                              *   you have to clear it at the end, exists only in presolving stage */
   SCIP_Bool*            bools2;             /**< cleared memory array, all entries are set to zero in initpre, if you use this
                                              *   you have to clear it at the end, exists only in presolving stage */
   SCIP_Bool*            bools3;             /**< cleared memory array, all entries are set to zero in initpre, if you use this
                                              *   you have to clear it at the end, exists only in presolving stage */
   SCIP_Bool*            bools4;             /**< cleared memory array, all entries are set to zero in initpre, if you use this
                                              *   you have to clear it at the end, exists only in presolving stage */
   SCIP_Real*            reals1;             /**< cleared memory array, all entries are set to zero in consinit, if you use this
                                              *   you have to clear it at the end */
   int                   ints1size;          /**< size of ints1 array */
   int                   ints2size;          /**< size of ints2 array */
   int                   longints1size;      /**< size of longints1 array */
   int                   longints2size;      /**< size of longints2 array */
   int                   bools1size;         /**< size of bools1 array */
   int                   bools2size;         /**< size of bools2 array */
   int                   bools3size;         /**< size of bools3 array */
   int                   bools4size;         /**< size of bools4 array */
   int                   reals1size;         /**< size of reals1 array */
   SCIP_EVENTHDLR*       eventhdlr;          /**< event handler for bound change events */
   SCIP_Real             maxcardbounddist;   /**< maximal relative distance from current node's dual bound to primal bound compared
                                              *   to best node's dual bound for separating knapsack cuts */
   int                   sepacardfreq;       /**< multiplier on separation frequency, how often knapsack cuts are separated */
   int                   maxrounds;          /**< maximal number of separation rounds per node (-1: unlimited) */
   int                   maxroundsroot;      /**< maximal number of separation rounds in the root node (-1: unlimited) */
   int                   maxsepacuts;        /**< maximal number of cuts separated per separation round */
   int                   maxsepacutsroot;    /**< maximal number of cuts separated per separation round in the root node */
   SCIP_Bool             disaggregation;     /**< should disaggregation of knapsack constraints be allowed in preprocessing? */
   SCIP_Bool             simplifyinequalities;/**< should presolving try to cancel down or delete coefficients in inequalities */
   SCIP_Bool             negatedclique;      /**< should negated clique information be used in solving process */
   SCIP_Bool             presolpairwise;     /**< should pairwise constraint comparison be performed in presolving? */
   SCIP_Bool             presolusehashing;   /**< should hash table be used for detecting redundant constraints in advance */
   SCIP_Bool             dualpresolving;     /**< should dual presolving steps be performed? */
   SCIP_Bool             usegubs;            /**< should GUB information be used for separation? */
   SCIP_Bool             detectcutoffbound;  /**< should presolving try to detect constraints parallel to the objective
                                              *   function defining an upper bound and prevent these constraints from
                                              *   entering the LP */
   SCIP_Bool             detectlowerbound;   /**< should presolving try to detect constraints parallel to the objective
                                              *   function defining a lower bound and prevent these constraints from
                                              *   entering the LP */
   SCIP_Bool             updatecliquepartitions; /**< should clique partition information be updated when old partition seems outdated? */
   SCIP_Real             cliqueextractfactor;/**< lower clique size limit for greedy clique extraction algorithm (relative to largest clique) */
   SCIP_Real             clqpartupdatefac;   /**< factor on the growth of global cliques to decide when to update a previous
                                              *   (negated) clique partition (used only if updatecliquepartitions is set to TRUE) */
#ifdef WITH_CARDINALITY_UPGRADE
   SCIP_Bool             upgdcardinality;    /**< if TRUE then try to update knapsack constraints to cardinality constraints */
   SCIP_Bool             upgradedcard;       /**< whether we have already upgraded knapsack constraints to cardinality constraints */
#endif
};


/** constraint data for knapsack constraints */
struct SCIP_ConsData
{
   SCIP_VAR**            vars;               /**< variables in knapsack constraint */
   SCIP_Longint*         weights;            /**< weights of variables in knapsack constraint */
   SCIP_EVENTDATA**      eventdata;          /**< event data for bound change events of the variables */
   int*                  cliquepartition;    /**< clique indices of the clique partition */
   int*                  negcliquepartition; /**< clique indices of the negated clique partition */
   SCIP_ROW*             row;                /**< corresponding LP row */
   SCIP_NLROW*           nlrow;              /**< corresponding NLP row */
   int                   nvars;              /**< number of variables in knapsack constraint */
   int                   varssize;           /**< size of vars, weights, and eventdata arrays */
   int                   ncliques;           /**< number of cliques in the clique partition */
   int                   nnegcliques;        /**< number of cliques in the negated clique partition */
   int                   ncliqueslastnegpart;/**< number of global cliques the last time a negated clique partition was computed */
   int                   ncliqueslastpart;   /**< number of global cliques the last time a clique partition was computed */
   SCIP_Longint          capacity;           /**< capacity of knapsack */
   SCIP_Longint          weightsum;          /**< sum of all weights */
   SCIP_Longint          onesweightsum;      /**< sum of weights of variables fixed to one */
   unsigned int          presolvedtiming:5;  /**< max level in which the knapsack constraint is already presolved */
   unsigned int          sorted:1;           /**< are the knapsack items sorted by weight? */
   unsigned int          cliquepartitioned:1;/**< is the clique partition valid? */
   unsigned int          negcliquepartitioned:1;/**< is the negated clique partition valid? */
   unsigned int          merged:1;           /**< are the constraint's equal variables already merged? */
   unsigned int          cliquesadded:1;     /**< were the cliques of the knapsack already added to clique table? */
   unsigned int          varsdeleted:1;      /**< were variables deleted after last cleanup? */
   unsigned int          existmultaggr:1;    /**< does this constraint contain multi-aggregations */
};

/** event data for bound changes events */
struct SCIP_EventData
{
   SCIP_CONS*            cons;               /**< knapsack constraint to process the bound change for */
   SCIP_Longint          weight;             /**< weight of variable */
   int                   filterpos;          /**< position of event in variable's event filter */
};


/** data structure to combine two sorting key values */
struct sortkeypair
{
   SCIP_Real             key1;               /**< first sort key value */
   SCIP_Real             key2;               /**< second sort key value */
};
typedef struct sortkeypair SORTKEYPAIR;

/** status of GUB constraint */
enum GUBVarstatus
{
   GUBVARSTATUS_UNINITIAL        = -1,       /** unintitialized variable status */
   GUBVARSTATUS_CAPACITYEXCEEDED =  0,       /** variable with weight exceeding the knapsack capacity */
   GUBVARSTATUS_BELONGSTOSET_R   =  1,       /** variable in noncovervars R */
   GUBVARSTATUS_BELONGSTOSET_F   =  2,       /** variable in noncovervars F */
   GUBVARSTATUS_BELONGSTOSET_C2  =  3,       /** variable in covervars C2 */
   GUBVARSTATUS_BELONGSTOSET_C1  =  4        /** variable in covervars C1 */
};
typedef enum GUBVarstatus GUBVARSTATUS;

/** status of variable in GUB constraint */
enum GUBConsstatus
{
   GUBCONSSTATUS_UNINITIAL         = -1,     /** unintitialized GUB constraint status */
   GUBCONSSTATUS_BELONGSTOSET_GR   =  0,     /** all GUB variables are in noncovervars R */
   GUBCONSSTATUS_BELONGSTOSET_GF   =  1,     /** all GUB variables are in noncovervars F (and noncovervars R) */
   GUBCONSSTATUS_BELONGSTOSET_GC2  =  2,     /** all GUB variables are in covervars C2 */
   GUBCONSSTATUS_BELONGSTOSET_GNC1 =  3,     /** some GUB variables are in covervars C1, others in noncovervars R or F */
   GUBCONSSTATUS_BELONGSTOSET_GOC1 =  4      /** all GUB variables are in covervars C1 */
};
typedef enum GUBConsstatus GUBCONSSTATUS;

/** data structure of GUB constraints */
struct SCIP_GUBCons
{
   int*                  gubvars;            /**< indices of GUB variables in knapsack constraint */
   GUBVARSTATUS*         gubvarsstatus;      /**< status of GUB variables */
   int                   ngubvars;           /**< number of GUB variables */
   int                   gubvarssize;        /**< size of gubvars array */
};
typedef struct SCIP_GUBCons SCIP_GUBCONS;

/** data structure of a set of GUB constraints */
struct SCIP_GUBSet
{
   SCIP_GUBCONS**        gubconss;           /**< GUB constraints in GUB set */
   GUBCONSSTATUS*        gubconsstatus;      /**< status of GUB constraints */
   int                   ngubconss;          /**< number of GUB constraints */
   int                   nvars;              /**< number of variables in knapsack constraint */
   int*                  gubconssidx;        /**< index of GUB constraint (in gubconss array) of each knapsack variable */
   int*                  gubvarsidx;         /**< index in GUB constraint (in gubvars array) of each knapsack variable  */
};
typedef struct SCIP_GUBSet SCIP_GUBSET;

/*
 * Local methods
 */

/** comparison method for two sorting key pairs */
static
SCIP_DECL_SORTPTRCOMP(compSortkeypairs)
{
   SORTKEYPAIR* sortkeypair1 = (SORTKEYPAIR*)elem1;
   SORTKEYPAIR* sortkeypair2 = (SORTKEYPAIR*)elem2;

   if( sortkeypair1->key1 < sortkeypair2->key1 )
      return -1;
   else if( sortkeypair1->key1 > sortkeypair2->key1 )
      return +1;
   else if( sortkeypair1->key2 < sortkeypair2->key2 )
      return -1;
   else if( sortkeypair1->key2 > sortkeypair2->key2 )
      return +1;
   else 
      return 0;
}

/** creates event data */
static
SCIP_RETCODE eventdataCreate(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_EVENTDATA**      eventdata,          /**< pointer to store event data */
   SCIP_CONS*            cons,               /**< constraint */
   SCIP_Longint          weight              /**< weight of variable */
   )
{
   assert(eventdata != NULL);

   SCIP_CALL( SCIPallocBlockMemory(scip, eventdata) );
   (*eventdata)->cons = cons;
   (*eventdata)->weight = weight;

   return SCIP_OKAY;
}  

/** frees event data */
static
SCIP_RETCODE eventdataFree(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_EVENTDATA**      eventdata           /**< pointer to event data */
   )
{
   assert(eventdata != NULL);

   SCIPfreeBlockMemory(scip, eventdata);

   return SCIP_OKAY;
}

/** sorts items in knapsack with nonincreasing weights */
static
void sortItems(
   SCIP_CONSDATA*        consdata            /**< constraint data */
   )
{
   assert(consdata != NULL);
   assert(consdata->nvars == 0 || consdata->vars != NULL);
   assert(consdata->nvars == 0 || consdata->weights != NULL);
   assert(consdata->nvars == 0 || consdata->eventdata != NULL);
   assert(consdata->nvars == 0 || (consdata->cliquepartition != NULL && consdata->negcliquepartition != NULL));

   if( !consdata->sorted )
   {
      int pos;
      int lastcliquenum;
      int v;

      /* sort of five joint arrays of Long/pointer/pointer/ints/ints,
       * sorted by first array in non-increasing order via sort template */
      SCIPsortDownLongPtrPtrIntInt(
         consdata->weights,
         (void**)consdata->vars,
         (void**)consdata->eventdata,
         consdata->cliquepartition,
         consdata->negcliquepartition,
         consdata->nvars);

      v = consdata->nvars - 1;
      /* sort all items with same weight according to their variable index, used for hash value for fast pairwise comparison of all constraints */
      while( v >= 0 )
      {
         int w = v - 1;

         while( w >= 0 && consdata->weights[v] == consdata->weights[w] )
            --w;

         if( v - w > 1 )
         {
            /* sort all corresponding parts of arrays for which the weights are equal by using the variable index */
            SCIPsortPtrPtrIntInt(
               (void**)(&(consdata->vars[w+1])),
               (void**)(&(consdata->eventdata[w+1])),
               &(consdata->cliquepartition[w+1]),
               &(consdata->negcliquepartition[w+1]),
               SCIPvarComp,
               v - w);
         }
         v = w;
      }

      /* we need to make sure that our clique numbers of our normal clique will be in increasing order without gaps */
      if( consdata->cliquepartitioned )
      {
         lastcliquenum = 0;

         for( pos = 0; pos < consdata->nvars; ++pos )
         {
            /* if the clique number in the normal clique at position pos is greater than the last found clique number the
             * partition is invalid */
            if( consdata->cliquepartition[pos] > lastcliquenum )
            {
               consdata->cliquepartitioned = FALSE;
               break;
            }
            else if( consdata->cliquepartition[pos] == lastcliquenum )
               ++lastcliquenum;
         }
      }
      /* we need to make sure that our clique numbers of our negated clique will be in increasing order without gaps */
      if( consdata->negcliquepartitioned )
      {
         lastcliquenum = 0;

         for( pos = 0; pos < consdata->nvars; ++pos )
         {
            /* if the clique number in the negated clique at position pos is greater than the last found clique number the
             * partition is invalid */
            if( consdata->negcliquepartition[pos] > lastcliquenum )
            {
               consdata->negcliquepartitioned = FALSE;
               break;
            }
            else if( consdata->negcliquepartition[pos] == lastcliquenum )
               ++lastcliquenum;
         }
      }

      consdata->sorted = TRUE;
   }
#ifndef NDEBUG
   {
      /* check if the weight array is sorted in a non-increasing way */
      int i;
      for( i = 0; i < consdata->nvars-1; ++i )
         assert(consdata->weights[i] >= consdata->weights[i+1]);
   }
#endif
}

/** calculates a partition of the variables into cliques */
static
SCIP_RETCODE calcCliquepartition(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONSHDLRDATA*    conshdlrdata,       /**< knapsack constraint handler data */
   SCIP_CONSDATA*        consdata,           /**< constraint data */
   SCIP_Bool             normalclique,       /**< Should normal cliquepartition be created? */
   SCIP_Bool             negatedclique       /**< Should negated cliquepartition be created? */
   )
{
   SCIP_Bool ispartitionoutdated;
   SCIP_Bool isnegpartitionoutdated;
   assert(consdata != NULL);
   assert(consdata->nvars == 0 || (consdata->cliquepartition != NULL && consdata->negcliquepartition != NULL));

   /* rerun eventually if number of global cliques increased considerably since last partition */
   ispartitionoutdated = (conshdlrdata->updatecliquepartitions && consdata->ncliques > 1
         && SCIPgetNCliques(scip) >= (int)(conshdlrdata->clqpartupdatefac * consdata->ncliqueslastpart));

   if( normalclique && ( !consdata->cliquepartitioned || ispartitionoutdated ) )
   {
      SCIP_CALL( SCIPcalcCliquePartition(scip, consdata->vars, consdata->nvars, consdata->cliquepartition, &consdata->ncliques) );
      consdata->cliquepartitioned = TRUE;
      consdata->ncliqueslastpart = SCIPgetNCliques(scip);
   }

   /* rerun eventually if number of global cliques increased considerably since last negated partition */
   isnegpartitionoutdated = (conshdlrdata->updatecliquepartitions && consdata->nnegcliques > 1
         && SCIPgetNCliques(scip) >= (int)(conshdlrdata->clqpartupdatefac * consdata->ncliqueslastnegpart));

   if( negatedclique && (!consdata->negcliquepartitioned || isnegpartitionoutdated) )
   {
      SCIP_CALL( SCIPcalcNegatedCliquePartition(scip, consdata->vars, consdata->nvars, consdata->negcliquepartition, &consdata->nnegcliques) );
      consdata->negcliquepartitioned = TRUE;
      consdata->ncliqueslastnegpart = SCIPgetNCliques(scip);
   }
   assert(!consdata->cliquepartitioned || consdata->ncliques <= consdata->nvars);
   assert(!consdata->negcliquepartitioned || consdata->nnegcliques <= consdata->nvars);

   return SCIP_OKAY;
}

/** installs rounding locks for the given variable in the given knapsack constraint */
static
SCIP_RETCODE lockRounding(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONS*            cons,               /**< knapsack constraint */
   SCIP_VAR*             var                 /**< variable of constraint entry */
   )
{
   SCIP_CALL( SCIPlockVarCons(scip, var, cons, FALSE, TRUE) );

   return SCIP_OKAY;
}

/** removes rounding locks for the given variable in the given knapsack constraint */
static
SCIP_RETCODE unlockRounding(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONS*            cons,               /**< knapsack constraint */
   SCIP_VAR*             var                 /**< variable of constraint entry */
   )
{
   SCIP_CALL( SCIPunlockVarCons(scip, var, cons, FALSE, TRUE) );

   return SCIP_OKAY;
}

/** catches bound change events for variables in knapsack */
static
SCIP_RETCODE catchEvents(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONS*            cons,               /**< constraint */
   SCIP_CONSDATA*        consdata,           /**< constraint data */
   SCIP_EVENTHDLR*       eventhdlr           /**< event handler to call for the event processing */
   )
{
   int i;

   assert(cons != NULL);
   assert(consdata != NULL);
   assert(consdata->nvars == 0 || consdata->vars != NULL);
   assert(consdata->nvars == 0 || consdata->weights != NULL);
   assert(consdata->nvars == 0 || consdata->eventdata != NULL);

   for( i = 0; i < consdata->nvars; i++)
   {
      SCIP_CALL( eventdataCreate(scip, &consdata->eventdata[i], cons, consdata->weights[i]) );
      SCIP_CALL( SCIPcatchVarEvent(scip, consdata->vars[i], EVENTTYPE_KNAPSACK,
            eventhdlr, consdata->eventdata[i], &consdata->eventdata[i]->filterpos) );
   }

   return SCIP_OKAY;
}

/** drops bound change events for variables in knapsack */
static
SCIP_RETCODE dropEvents(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONSDATA*        consdata,           /**< constraint data */
   SCIP_EVENTHDLR*       eventhdlr           /**< event handler to call for the event processing */
   )
{
   int i;

   assert(consdata != NULL);
   assert(consdata->nvars == 0 || consdata->vars != NULL);
   assert(consdata->nvars == 0 || consdata->weights != NULL);
   assert(consdata->nvars == 0 || consdata->eventdata != NULL);

   for( i = 0; i < consdata->nvars; i++)
   {
      SCIP_CALL( SCIPdropVarEvent(scip, consdata->vars[i], EVENTTYPE_KNAPSACK,
            eventhdlr, consdata->eventdata[i], consdata->eventdata[i]->filterpos) );
      SCIP_CALL( eventdataFree(scip, &consdata->eventdata[i]) );
   }

   return SCIP_OKAY;
}

/** ensures, that vars and vals arrays can store at least num entries */
static
SCIP_RETCODE consdataEnsureVarsSize(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONSDATA*        consdata,           /**< knapsack constraint data */
   int                   num,                /**< minimum number of entries to store */
   SCIP_Bool             transformed         /**< is constraint from transformed problem? */
   )
{
   assert(consdata != NULL);
   assert(consdata->nvars <= consdata->varssize);

   if( num > consdata->varssize )
   {
      int newsize;

      newsize = SCIPcalcMemGrowSize(scip, num);
      SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->vars, consdata->varssize, newsize) );
      SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->weights, consdata->varssize, newsize) );
      if( transformed )
      {
         SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->eventdata, consdata->varssize, newsize) );
         SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->cliquepartition, consdata->varssize, newsize) );
         SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &consdata->negcliquepartition, consdata->varssize, newsize) );
      }
      else
      {
         assert(consdata->eventdata == NULL);
         assert(consdata->cliquepartition == NULL);
         assert(consdata->negcliquepartition == NULL);
      }
      consdata->varssize = newsize;
   }
   assert(num <= consdata->varssize);

   return SCIP_OKAY;
}

/** updates all weight sums for fixed and unfixed variables */
static
void updateWeightSums(
   SCIP_CONSDATA*        consdata,           /**< knapsack constraint data */
   SCIP_VAR*             var,                /**< variable for this weight */
   SCIP_Longint          weightdelta         /**< difference between the old and the new weight of the variable */
   )
{
   assert(consdata != NULL);
   assert(var != NULL);

   consdata->weightsum += weightdelta;

   if( SCIPvarGetLbLocal(var) > 0.5 )
      consdata->onesweightsum += weightdelta;

   assert(consdata->weightsum >= 0);
   assert(consdata->onesweightsum >= 0);
}

/** creates knapsack constraint data */
static
SCIP_RETCODE consdataCreate(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONSDATA**       consdata,           /**< pointer to store constraint data */
   int                   nvars,              /**< number of variables in knapsack */
   SCIP_VAR**            vars,               /**< variables of knapsack */
   SCIP_Longint*         weights,            /**< weights of knapsack items */
   SCIP_Longint          capacity            /**< capacity of knapsack */
   )
{
   int v;
   SCIP_Longint constant;

   assert(consdata != NULL);

   SCIP_CALL( SCIPallocBlockMemory(scip, consdata) );

   constant = 0L;
   (*consdata)->vars = NULL;
   (*consdata)->weights = NULL;
   (*consdata)->nvars = 0;
   if( nvars > 0 )
   {
      SCIP_VAR** varsbuffer;
      SCIP_Longint* weightsbuffer;
      int k;

      SCIP_CALL( SCIPallocBufferArray(scip, &varsbuffer, nvars) );
      SCIP_CALL( SCIPallocBufferArray(scip, &weightsbuffer, nvars) );

      k = 0;
      for( v = 0; v < nvars; ++v )
      {
         assert(vars[v] != NULL);
         assert(SCIPvarIsBinary(vars[v]));

         /* all weight have to be non negative */
         assert( weights[v] >= 0 );

         if( weights[v] > 0 )
         {
            /* treat fixed variables as constants if problem compression is enabled */
            if( SCIPisConsCompressionEnabled(scip) && SCIPvarGetLbGlobal(vars[v]) > SCIPvarGetUbGlobal(vars[v]) - 0.5 )
            {
               /* only if the variable is fixed to 1, we add its weight to the constant */
               if( SCIPvarGetUbGlobal(vars[v]) > 0.5 )
                  constant += weights[v];
            }
            else
            {
               varsbuffer[k] = vars[v];
               weightsbuffer[k] = weights[v];
               ++k;
            }
         }
      }
      assert(k >= 0);
      assert(constant >= 0);

      (*consdata)->nvars = k;

      /* copy the active variables and weights into the constraint data structure */
      if( k > 0 )
      {
         SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*consdata)->vars, varsbuffer, k) );
         SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*consdata)->weights, weightsbuffer, k) );
      }

      /* free buffer storage */
      SCIPfreeBufferArray(scip, &weightsbuffer);
      SCIPfreeBufferArray(scip, &varsbuffer);
   }

   (*consdata)->varssize = (*consdata)->nvars;
   (*consdata)->capacity = capacity - constant;
   (*consdata)->eventdata = NULL;
   (*consdata)->cliquepartition = NULL;
   (*consdata)->negcliquepartition = NULL;
   (*consdata)->row = NULL;
   (*consdata)->nlrow = NULL;
   (*consdata)->weightsum = 0;
   (*consdata)->onesweightsum = 0;
   (*consdata)->ncliques = 0;
   (*consdata)->nnegcliques = 0;
   (*consdata)->presolvedtiming = 0;
   (*consdata)->sorted = FALSE;
   (*consdata)->cliquepartitioned = FALSE;
   (*consdata)->negcliquepartitioned = FALSE;
   (*consdata)->ncliqueslastpart = -1;
   (*consdata)->ncliqueslastnegpart = -1;
   (*consdata)->merged = FALSE;
   (*consdata)->cliquesadded = FALSE;
   (*consdata)->varsdeleted = FALSE;
   (*consdata)->existmultaggr = FALSE;

   /* get transformed variables, if we are in the transformed problem */
   if( SCIPisTransformed(scip) )
   {
      SCIP_CALL( SCIPgetTransformedVars(scip, (*consdata)->nvars, (*consdata)->vars, (*consdata)->vars) );

      for( v = 0; v < (*consdata)->nvars; v++ )
      {
         SCIP_VAR* var = SCIPvarGetProbvar((*consdata)->vars[v]);
         assert(var != NULL);
         (*consdata)->existmultaggr = (*consdata)->existmultaggr || (SCIPvarGetStatus(var) == SCIP_VARSTATUS_MULTAGGR);
      }

      /* allocate memory for additional data structures */
      SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*consdata)->eventdata, (*consdata)->nvars) );
      SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*consdata)->cliquepartition, (*consdata)->nvars) );
      SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*consdata)->negcliquepartition, (*consdata)->nvars) );
   }

   /* calculate sum of weights and capture variables */
   for( v = 0; v < (*consdata)->nvars; ++v )
   {
      /* calculate sum of weights */
      updateWeightSums(*consdata, (*consdata)->vars[v], (*consdata)->weights[v]);

      /* capture variables */
      SCIP_CALL( SCIPcaptureVar(scip, (*consdata)->vars[v]) );
   }
   return SCIP_OKAY;
}

/** frees knapsack constraint data */
static
SCIP_RETCODE consdataFree(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONSDATA**       consdata,           /**< pointer to the constraint data */
   SCIP_EVENTHDLR*       eventhdlr           /**< event handler to call for the event processing */
   )
{
   assert(consdata != NULL);
   assert(*consdata != NULL);

   if( (*consdata)->row != NULL )
   {
      SCIP_CALL( SCIPreleaseRow(scip, &(*consdata)->row) );
   }
   if( (*consdata)->nlrow != NULL )
   {
      SCIP_CALL( SCIPreleaseNlRow(scip, &(*consdata)->nlrow) );
   }
   if( (*consdata)->eventdata != NULL )
   {
      SCIP_CALL( dropEvents(scip, *consdata, eventhdlr) );
      SCIPfreeBlockMemoryArray(scip, &(*consdata)->eventdata, (*consdata)->varssize);
   }
   if( (*consdata)->negcliquepartition != NULL )
   {
      SCIPfreeBlockMemoryArray(scip, &(*consdata)->negcliquepartition, (*consdata)->varssize);
   }
   if( (*consdata)->cliquepartition != NULL )
   {
      SCIPfreeBlockMemoryArray(scip, &(*consdata)->cliquepartition, (*consdata)->varssize);
   }
   if( (*consdata)->vars != NULL )
   {
      int v;

      /* release variables */
      for( v = 0; v < (*consdata)->nvars; v++ )
      {
         assert((*consdata)->vars[v] != NULL);
         SCIP_CALL( SCIPreleaseVar(scip, &((*consdata)->vars[v])) );
      }

      assert( (*consdata)->weights != NULL );
      assert( (*consdata)->varssize > 0 );
      SCIPfreeBlockMemoryArray(scip, &(*consdata)->vars, (*consdata)->varssize);
      SCIPfreeBlockMemoryArray(scip, &(*consdata)->weights, (*consdata)->varssize);
   }

   SCIPfreeBlockMemory(scip, consdata);

   return SCIP_OKAY;
}

/** changes a single weight in knapsack constraint data */
static
void consdataChgWeight(
   SCIP_CONSDATA*        consdata,           /**< knapsack constraint data */
   int                   item,               /**< item number */
   SCIP_Longint          newweight           /**< new weight of item */
   )
{
   SCIP_Longint oldweight;
   SCIP_Longint weightdiff;

   assert(consdata != NULL);
   assert(0 <= item && item < consdata->nvars);

   oldweight = consdata->weights[item];
   weightdiff = newweight - oldweight;
   consdata->weights[item] = newweight;

   /* update weight sums for all and fixed variables */
   updateWeightSums(consdata, consdata->vars[item], weightdiff);

   if( consdata->eventdata != NULL )
   {
      assert(consdata->eventdata[item] != NULL);
      assert(consdata->eventdata[item]->weight == oldweight);
      consdata->eventdata[item]->weight = newweight;
   }

   consdata->presolvedtiming = 0;
   consdata->sorted = FALSE;

   /* recalculate cliques extraction after a weight was increased */
   if( oldweight < newweight )
   {
      consdata->cliquesadded = FALSE;
   }
}

/** creates LP row corresponding to knapsack constraint */
static 
SCIP_RETCODE createRelaxation(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONS*            cons                /**< knapsack constraint */
   )
{
   SCIP_CONSDATA* consdata;
   int i;

   consdata = SCIPconsGetData(cons);
   assert(consdata != NULL);
   assert(consdata->row == NULL);

   SCIP_CALL( SCIPcreateEmptyRowCons(scip, &consdata->row, cons, SCIPconsGetName(cons),
         -SCIPinfinity(scip), (SCIP_Real)consdata->capacity,
         SCIPconsIsLocal(cons), SCIPconsIsModifiable(cons), SCIPconsIsRemovable(cons)) );

   SCIP_CALL( SCIPcacheRowExtensions(scip, consdata->row) );
   for( i = 0; i < consdata->nvars; ++i )
   {
      SCIP_CALL( SCIPaddVarToRow(scip, consdata->row, consdata->vars[i], (SCIP_Real)consdata->weights[i]) );
   }
   SCIP_CALL( SCIPflushRowExtensions(scip, consdata->row) );

   return SCIP_OKAY;
}

/** adds linear relaxation of knapsack constraint to the LP */
static
SCIP_RETCODE addRelaxation(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONS*            cons,               /**< knapsack constraint */
   SCIP_Bool*            cutoff              /**< whether a cutoff has been detected */
   )
{
   SCIP_CONSDATA* consdata;

   assert( cutoff != NULL );
   *cutoff = FALSE;

   consdata = SCIPconsGetData(cons);
   assert(consdata != NULL);

   if( consdata->row == NULL )
   {
      SCIP_CALL( createRelaxation(scip, cons) );
   }
   assert(consdata->row != NULL);

   /* insert LP row as cut */
   if( !SCIProwIsInLP(consdata->row) )
   {
      SCIPdebugMsg(scip, "adding relaxation of knapsack constraint <%s> (capacity %" SCIP_LONGINT_FORMAT "): ",
         SCIPconsGetName(cons), consdata->capacity);
      SCIPdebug( SCIP_CALL(SCIPprintRow(scip, consdata->row, NULL)) );
      SCIP_CALL( SCIPaddRow(scip, consdata->row, FALSE, cutoff) );
   }

   return SCIP_OKAY;
}

/** adds knapsack constraint as row to the NLP, if not added yet */
static
SCIP_RETCODE addNlrow(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONS*            cons                /**< knapsack constraint */
   )
{
   SCIP_CONSDATA* consdata;

   assert(SCIPisNLPConstructed(scip));

   /* skip deactivated, redundant, or local linear constraints (the NLP does not allow for local rows at the moment) */
   if( !SCIPconsIsActive(cons) || !SCIPconsIsChecked(cons) || SCIPconsIsLocal(cons) )
      return SCIP_OKAY;

   consdata = SCIPconsGetData(cons);
   assert(consdata != NULL);

   if( consdata->nlrow == NULL )
   {
      SCIP_Real* coefs;
      int i;

      SCIP_CALL( SCIPallocBufferArray(scip, &coefs, consdata->nvars) );
      for( i = 0; i < consdata->nvars; ++i )
         coefs[i] = (SCIP_Real)consdata->weights[i];  /*lint !e613*/

      SCIP_CALL( SCIPcreateNlRow(scip, &consdata->nlrow, SCIPconsGetName(cons), 0.0,
            consdata->nvars, consdata->vars, coefs, NULL,
            -SCIPinfinity(scip), (SCIP_Real)consdata->capacity, SCIP_EXPRCURV_LINEAR) );

      assert(consdata->nlrow != NULL);

      SCIPfreeBufferArray(scip, &coefs);
   }

   if( !SCIPnlrowIsInNLP(consdata->nlrow) )
   {
      SCIP_CALL( SCIPaddNlRow(scip, consdata->nlrow) );
   }

   return SCIP_OKAY;
}

/** checks knapsack constraint for feasibility of given solution: returns TRUE iff constraint is feasible */
static
SCIP_RETCODE checkCons(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONS*            cons,               /**< constraint to check */
   SCIP_SOL*             sol,                /**< solution to check, NULL for current solution */
   SCIP_Bool             checklprows,        /**< Do constraints represented by rows in the current LP have to be checked? */
   SCIP_Bool             printreason,        /**< Should the reason for the violation be printed? */
   SCIP_Bool*            violated            /**< pointer to store whether the constraint is violated */
   )
{
   SCIP_CONSDATA* consdata;

   assert(violated != NULL);

   consdata = SCIPconsGetData(cons);
   assert(consdata != NULL);

   SCIPdebugMsg(scip, "checking knapsack constraint <%s> for feasibility of solution %p (lprows=%u)\n",
      SCIPconsGetName(cons), (void*)sol, checklprows);

   *violated = FALSE;

   if( checklprows || consdata->row == NULL || !SCIProwIsInLP(consdata->row) )
   {
      SCIP_Real sum;
      SCIP_Longint integralsum;
      SCIP_Bool ishuge;
      SCIP_Real absviol;
      SCIP_Real relviol;
      int v;

      /* increase age of constraint; age is reset to zero, if a violation was found only in case we are in
       * enforcement
       */
      if( sol == NULL )
      {
         SCIP_CALL( SCIPincConsAge(scip, cons) );
      }

      sum = 0.0;
      integralsum = 0;
      /* we perform a more exact comparison if the capacity does not exceed the huge value */
      if( SCIPisHugeValue(scip, (SCIP_Real) consdata->capacity) )
      {
         ishuge = TRUE;

         /* sum over all weight times the corresponding solution value */
         for( v = consdata->nvars - 1; v >= 0; --v )
         {
            assert(SCIPvarIsBinary(consdata->vars[v]));
            sum += consdata->weights[v] * SCIPgetSolVal(scip, sol, consdata->vars[v]);
         }
      }
      else
      {
         ishuge = FALSE;

         /* sum over all weight for which the variable has a solution value of 1 in feastol */
         for( v = consdata->nvars - 1; v >= 0; --v )
         {
            assert(SCIPvarIsBinary(consdata->vars[v]));

            if( SCIPgetSolVal(scip, sol, consdata->vars[v]) > 0.5 )
               integralsum += consdata->weights[v];
         }
      }

      /* calculate constraint violation and update it in solution */
      absviol = ishuge ? sum : (SCIP_Real)integralsum;
      absviol -= consdata->capacity;
      relviol = SCIPrelDiff(absviol + consdata->capacity, (SCIP_Real)consdata->capacity);
      if( sol != NULL )
         SCIPupdateSolLPConsViolation(scip, sol, absviol, relviol);

      if( SCIPisFeasPositive(scip, absviol) )
      {
         *violated = TRUE;

         /* only reset constraint age if we are in enforcement */
         if( sol == NULL )
         {
            SCIP_CALL( SCIPresetConsAge(scip, cons) );
         }

         if( printreason )
         {
            SCIP_CALL( SCIPprintCons(scip, cons, NULL) );

            SCIPinfoMessage(scip, NULL, ";\n");
            SCIPinfoMessage(scip, NULL, "violation: the capacity is violated by %.15g\n", absviol);
         }
      }
   }

   return SCIP_OKAY;
}

/* IDX computes the integer index for the optimal solution array */
#define IDX(j,d) ((j)*(intcap)+(d))

/** solves knapsack problem in maximization form exactly using dynamic programming;
 *  if needed, one can provide arrays to store all selected items and all not selected items
 *
 * @note in case you provide the solitems or nonsolitems array you also have to provide the counter part, as well
 *
 * @note the algorithm will first compute a greedy solution and terminate
 *       if the greedy solution is proven to be optimal.
 *       The dynamic programming algorithm runs with a time and space complexity
 *       of O(nitems * capacity).
 *
 * @todo If only the objective is relevant, it is easy to change the code to use only one slice with O(capacity) space.
 *       There are recursive methods (see the book by Kellerer et al.) that require O(capacity) space, but it remains
 *       to be checked whether they are faster and whether they can reconstruct the solution.
 *       Dembo and Hammer (see Kellerer et al. Section 5.1.3, page 126) found a method that relies on a fast probing method.
 *       This fixes additional elements to 0 or 1 similar to a reduced cost fixing.
 *       This could be implemented, however, it would be technically a bit cumbersome,
 *       since one needs the greedy solution and the LP-value for this.
 *       This is currently only available after the redundant items have already been sorted out.
 */
SCIP_RETCODE SCIPsolveKnapsackExactly(
   SCIP*                 scip,               /**< SCIP data structure */
   int                   nitems,             /**< number of available items */
   SCIP_Longint*         weights,            /**< item weights */
   SCIP_Real*            profits,            /**< item profits */
   SCIP_Longint          capacity,           /**< capacity of knapsack */
   int*                  items,              /**< item numbers */
   int*                  solitems,           /**< array to store items in solution, or NULL */
   int*                  nonsolitems,        /**< array to store items not in solution, or NULL */
   int*                  nsolitems,          /**< pointer to store number of items in solution, or NULL */
   int*                  nnonsolitems,       /**< pointer to store number of items not in solution, or NULL */
   SCIP_Real*            solval,             /**< pointer to store optimal solution value, or NULL */
   SCIP_Bool*            success             /**< pointer to store if an error occured during solving
                                              *   (normally a memory problem) */
   )
{
   SCIP_RETCODE retcode;
   SCIP_Real* tempsort;
   SCIP_Real* optvalues;
   int intcap;
   int d;
   int j;
   int greedymedianpos;
   SCIP_Longint weightsum;
   int* myitems;
   SCIP_Longint* myweights;
   SCIP_Real* realweights;
   int* allcurrminweight;
   SCIP_Real* myprofits;
   int nmyitems;
   SCIP_Longint gcd;
   SCIP_Longint minweight;
   SCIP_Longint maxweight;
   int currminweight;
   SCIP_Longint greedysolweight;
   SCIP_Real greedysolvalue;
   SCIP_Real greedyupperbound;
   SCIP_Bool eqweights;
   SCIP_Bool intprofits;

   assert(weights != NULL);
   assert(profits != NULL);
   assert(capacity >= 0);
   assert(items != NULL);
   assert(nitems >= 0);
   assert(success != NULL);

   *success = TRUE;

#ifndef NDEBUG
   for( j = nitems - 1; j >= 0; --j )
      assert(weights[j] >= 0);
#endif

   SCIPdebugMsg(scip, "Solving knapsack exactly.\n");

   /* initializing solution value */
   if( solval != NULL )
      *solval = 0.0;

   /* init solution information */
   if( solitems != NULL )
   {
      assert(items != NULL);
      assert(nsolitems != NULL);
      assert(nonsolitems != NULL);
      assert(nnonsolitems != NULL);

      *nnonsolitems = 0;
      *nsolitems = 0;
   }

   /* allocate temporary memory */
   SCIP_CALL( SCIPallocBufferArray(scip, &myweights, nitems) );
   SCIP_CALL( SCIPallocBufferArray(scip, &myprofits, nitems) );
   SCIP_CALL( SCIPallocBufferArray(scip, &myitems, nitems) );
   nmyitems = 0;
   weightsum = 0;
   minweight = SCIP_LONGINT_MAX;
   maxweight = 0;

   /* remove unnecessary items */
   for( j = 0; j < nitems; ++j )
   {
      assert(0 <= weights[j] && weights[j] < SCIP_LONGINT_MAX);

      /* item does not fit */
      if( weights[j] > capacity )
      {
         if( solitems != NULL )
            nonsolitems[(*nnonsolitems)++] = items[j]; /*lint !e413*/
      }
      /* item is not profitable */
      else if( profits[j] <= 0.0 )
      {
         if( solitems != NULL )
            nonsolitems[(*nnonsolitems)++] = items[j]; /*lint !e413*/
      }
      /* item always fits */
      else if( weights[j] == 0 )
      {
         if( solitems != NULL )
            solitems[(*nsolitems)++] = items[j]; /*lint !e413*/

         if( solval != NULL )
            *solval += profits[j];
      }
      /* all important items */
      else
      {
         myweights[nmyitems] = weights[j];
         myprofits[nmyitems] = profits[j];
         myitems[nmyitems] = items[j];

         /* remember smallest item */
         if( myweights[nmyitems] < minweight )
            minweight = myweights[nmyitems];

         /* remember bigest item */
         if( myweights[nmyitems] > maxweight )
            maxweight = myweights[nmyitems];

         weightsum += myweights[nmyitems];
         ++nmyitems;
      }
   }

   intprofits = TRUE;
   /* check if all profits are integer to strengthen the upper bound on the greedy solution */
   for( j = 0; j < nmyitems && intprofits; ++j )
      intprofits = intprofits && SCIPisIntegral(scip, myprofits[j]);

   /* if no item is left then goto end */
   if( nmyitems == 0 )
   {
      SCIPdebugMsg(scip, "After preprocessing no items are left.\n");

      goto TERMINATE;
   }

   /* if all items fit, we also do not need to do the expensive stuff later on */
   if( weightsum > 0 && weightsum <= capacity )
   {
      SCIPdebugMsg(scip, "After preprocessing all items fit into knapsack.\n");

      for( j = nmyitems - 1; j >= 0; --j )
      {
         if( solitems != NULL )
            solitems[(*nsolitems)++] = myitems[j]; /*lint !e413*/

         if( solval != NULL )
            *solval += myprofits[j];
      }

      goto TERMINATE;
   }

   assert(0 < minweight && minweight <= capacity );
   assert(0 < maxweight && maxweight <= capacity);

   /* make weights relatively prime */
   eqweights = TRUE;
   if( maxweight > 1 )
   {
      /* determine greatest common divisor */
      gcd = myweights[nmyitems - 1];
      for( j = nmyitems - 2; j >= 0 && gcd >= 2; --j )
         gcd = SCIPcalcGreComDiv(gcd, myweights[j]);

      SCIPdebugMsg(scip, "Gcd is %" SCIP_LONGINT_FORMAT ".\n", gcd);

      /* divide by greatest common divisor */
      if( gcd > 1 )
      {
         for( j = nmyitems - 1; j >= 0; --j )
         {
            myweights[j] /= gcd;
            eqweights = eqweights && (myweights[j] == 1);
         }
         capacity /= gcd;
         minweight /= gcd;
      }
      else
         eqweights = FALSE;
   }
   assert(minweight <= capacity);

   /* if only one item fits, then take the best */
   if( minweight > capacity / 2 )
   {
      int p;

      SCIPdebugMsg(scip, "Only one item fits into knapsack, so take the best.\n");

      p = nmyitems - 1;

      /* find best item */
      for( j = nmyitems - 2; j >= 0; --j )
      {
         if( myprofits[j] > myprofits[p] )
            p = j;
      }

      /* update solution information */
      if( solitems != NULL )
      {
         assert(nsolitems != NULL && nonsolitems != NULL && nnonsolitems != NULL);

         solitems[(*nsolitems)++] = myitems[p];
         for( j = nmyitems - 1; j >= 0; --j )
         {
            if( j != p )
               nonsolitems[(*nnonsolitems)++] = myitems[j];
         }
      }
      /* update solution value */
      if( solval != NULL )
         *solval += myprofits[p];

      goto TERMINATE;
   }

   /* if all items have the same weight, then take the best */
   if( eqweights )
   {
      SCIP_Real addval = 0.0;

      SCIPdebugMsg(scip, "All weights are equal, so take the best.\n");

      SCIPsortDownRealIntLong(myprofits, myitems, myweights, nmyitems);

      /* update solution information */
      if( solitems != NULL || solval != NULL )
      {
         SCIP_Longint i;

         /* if all items would fit we had handled this case before */
         assert((SCIP_Longint) nmyitems > capacity);
         assert(nsolitems != NULL && nonsolitems != NULL && nnonsolitems != NULL);

         /* take the first best items into the solution */
         for( i = capacity - 1; i >= 0; --i )
         {
            if( solitems != NULL )
               solitems[(*nsolitems)++] = myitems[i];
            addval += myprofits[i];
         }

         if( solitems != NULL )
         {
            /* the rest are not in the solution */
            for( i = nmyitems - 1; i >= capacity; --i )
               nonsolitems[(*nnonsolitems)++] = myitems[i];
         }
      }
      /* update solution value */
      if( solval != NULL )
      {
         assert(addval > 0.0);
         *solval += addval;
      }

      goto TERMINATE;
   }

   SCIPdebugMsg(scip, "Determine greedy solution.\n");

   /* sort myitems (plus corresponding arrays myweights and myprofits) such that
    * p_1/w_1 >= p_2/w_2 >= ... >= p_n/w_n, this is only used for the greedy solution
    */
   SCIP_CALL( SCIPallocBufferArray(scip, &tempsort, nmyitems) );
   SCIP_CALL( SCIPallocBufferArray(scip, &realweights, nmyitems) );

   for( j = 0; j < nmyitems; ++j )
   {
      tempsort[j] = myprofits[j]/((SCIP_Real) myweights[j]);
      realweights[j] = (SCIP_Real)myweights[j];
   }

   SCIPselectWeightedDownRealLongRealInt(tempsort, myweights, myprofits, myitems, realweights,
      (SCIP_Real)capacity, nmyitems, &greedymedianpos);

   SCIPfreeBufferArray(scip, &realweights);
   SCIPfreeBufferArray(scip, &tempsort);

   /* initialize values for greedy solution information */
   greedysolweight = 0;
   greedysolvalue = 0.0;

   /* determine greedy solution */
   for( j = 0; j < greedymedianpos; ++j )
   {
      assert(myweights[j] <= capacity);

      /* update greedy solution weight and value */
      greedysolweight += myweights[j];
      greedysolvalue += myprofits[j];
   }

   assert(0 < greedysolweight && greedysolweight <= capacity);
   assert(greedysolvalue > 0.0);

   /* If the greedy solution is optimal by comparing to the LP solution, we take this solution. This happens if:
    * - the greedy solution reaches the capacity, because then the LP solution is integral;
    * - the greedy solution has an objective that is at least the LP value rounded down in case that all profits are integer, too. */
   greedyupperbound = greedysolvalue + myprofits[j] * (SCIP_Real) (capacity - greedysolweight)/((SCIP_Real) myweights[j]);
   if( intprofits )
      greedyupperbound = SCIPfloor(scip, greedyupperbound);
   if( greedysolweight == capacity || SCIPisGE(scip, greedysolvalue, greedyupperbound) )
   {
      SCIPdebugMsg(scip, "Greedy solution is optimal.\n");

      /* update solution information */
      if( solitems != NULL )
      {
         int l;

         assert(nsolitems != NULL && nonsolitems != NULL && nnonsolitems != NULL);

         /* collect items */
         for( l = 0; l < j; ++l )
            solitems[(*nsolitems)++] = myitems[l];
         for ( ; l < nmyitems; ++l )
            nonsolitems[(*nnonsolitems)++] = myitems[l];
      }
      /* update solution value */
      if( solval != NULL )
      {
         assert(greedysolvalue > 0.0);
         *solval += greedysolvalue;
      }

      goto TERMINATE;
   }

   /* in the following table we do not need the first minweight columns */
   capacity -= (minweight - 1);

   /* we can only handle integers */
   if( capacity >= INT_MAX )
   {
      SCIPdebugMsg(scip, "Capacity is to big, so we cannot handle it here.\n");

      *success = FALSE;
      goto TERMINATE;
   }
   assert(capacity < INT_MAX);

   intcap = (int)capacity;
   assert(intcap >= 0);
   assert(nmyitems > 0);
   assert(sizeof(size_t) >= sizeof(int)); /*lint !e506*/ /* no following conversion should be messed up */

   /* this condition checks whether we will try to allocate a correct number of bytes and do not have an overflow, while
    * computing the size for the allocation
    */
   if( intcap < 0 || (intcap > 0 && (((size_t)nmyitems) > (SIZE_MAX / (size_t)intcap / sizeof(*optvalues)) || ((size_t)nmyitems) * ((size_t)intcap) * sizeof(*optvalues) > ((size_t)INT_MAX) )) ) /*lint !e571*/
   {
      SCIPdebugMsg(scip, "Too much memory (%lu) would be consumed.\n", (unsigned long) (((size_t)nmyitems) * ((size_t)intcap) * sizeof(*optvalues))); /*lint !e571*/

      *success = FALSE;
      goto TERMINATE;
   }

   /* allocate temporary memory and check for memory exceedance */
   retcode = SCIPallocBufferArray(scip, &optvalues, nmyitems * intcap);
   if( retcode == SCIP_NOMEMORY )
   {
      SCIPdebugMsg(scip, "Did not get enough memory.\n");

      *success = FALSE;
      goto TERMINATE;
   }
   else
   {
      SCIP_CALL( retcode );
   }

   SCIPdebugMsg(scip, "Start real exact algorithm.\n");

   /* we memorize at each step the current minimal weight to later on know which value in our optvalues matrix is valid;
    * each value entries of the j-th row of optvalues is valid if the index is >= allcurrminweight[j], otherwise it is
    * invalid; a second possibility would be to clear the whole optvalues, which should be more expensive than storing
    * 'nmyitem' values
    */
   SCIP_CALL( SCIPallocBufferArray(scip, &allcurrminweight, nmyitems) );
   assert(myweights[0] - minweight < INT_MAX);
   currminweight = (int) (myweights[0] - minweight);
   allcurrminweight[0] = currminweight;

   /* fills first row of dynamic programming table with optimal values */
   for( d = currminweight; d < intcap; ++d )
      optvalues[d] = myprofits[0];

   /* fills dynamic programming table with optimal values */
   for( j = 1; j < nmyitems; ++j )
   {
      int intweight;

      /* compute important part of weight, which will be represented in the table */
      intweight = (int)(myweights[j] - minweight);
      assert(0 <= intweight && intweight < intcap);

      /* copy all nonzeros from row above */
      for( d = currminweight; d < intweight && d < intcap; ++d )
         optvalues[IDX(j,d)] = optvalues[IDX(j-1,d)];

      /* update corresponding row */
      for( d = intweight; d < intcap; ++d )
      {
         /* if index d < current minweight then optvalues[IDX(j-1,d)] is not initialized, i.e. should be 0 */
         if( d < currminweight )
            optvalues[IDX(j,d)] = myprofits[j];
         else
         {
            SCIP_Real sumprofit;

            if( d - myweights[j] < currminweight )
               sumprofit = myprofits[j];
            else
               sumprofit = optvalues[IDX(j-1,(int)(d-myweights[j]))] + myprofits[j];

            optvalues[IDX(j,d)] = MAX(sumprofit, optvalues[IDX(j-1,d)]);
         }
      }

      /* update currminweight */
      if( intweight < currminweight )
         currminweight = intweight;

      allcurrminweight[j] = currminweight;
   }

   /* update optimal solution by following the table */
   if( solitems != NULL )
   {
      assert(nsolitems != NULL && nonsolitems != NULL && nnonsolitems != NULL);
      d = intcap - 1;

      SCIPdebugMsg(scip, "Fill the solution vector after solving exactly.\n");

      /* insert all items in (non-) solution vector */
      for( j = nmyitems - 1; j > 0; --j )
      {
         /* if the following condition holds this means all remaining items does not fit anymore */
         if( d < allcurrminweight[j] )
         {
            /* we cannot have exceeded our capacity */
            assert((SCIP_Longint) d >= -minweight);
            break;
         }

         /* collect solution items; the first condition means that no further item can fit anymore, but this does */
         if( d < allcurrminweight[j-1] || optvalues[IDX(j,d)] > optvalues[IDX(j-1,d)] )
         {
            solitems[(*nsolitems)++] = myitems[j];

            /* check that we do not have an underflow */
            assert(myweights[j] <= (INT_MAX + (SCIP_Longint) d));
            d = (int)(d - myweights[j]);
         }
         /* collect non-solution items */
         else
            nonsolitems[(*nnonsolitems)++] = myitems[j];
      }

      /* insert remaining items */
      if( d >= allcurrminweight[j] )
      {
         assert(j == 0);
         solitems[(*nsolitems)++] = myitems[j];
      }
      else
      {
         assert(j >= 0);
         assert(d < allcurrminweight[j]);

         for( ; j >= 0; --j )
            nonsolitems[(*nnonsolitems)++] = myitems[j];
      }

      assert(*nsolitems + *nnonsolitems == nitems);
   }

   /* update solution value */
   if( solval != NULL )
      *solval += optvalues[IDX(nmyitems-1,intcap-1)];
   SCIPfreeBufferArray(scip, &allcurrminweight);

   /* free all temporary memory */
   SCIPfreeBufferArray(scip, &optvalues);

 TERMINATE:
   SCIPfreeBufferArray(scip, &myitems);
   SCIPfreeBufferArray(scip, &myprofits);
   SCIPfreeBufferArray(scip, &myweights);

   return SCIP_OKAY;
}

/** solves knapsack problem in maximization form approximately by solving the LP-relaxation of the problem using Dantzig's
 *  method and rounding down the solution; if needed, one can provide arrays to store all selected items and all not
 *  selected items
 */
SCIP_RETCODE SCIPsolveKnapsackApproximately(
   SCIP*                 scip,               /**< SCIP data structure */
   int                   nitems,             /**< number of available items */
   SCIP_Longint*         weights,            /**< item weights */
   SCIP_Real*            profits,            /**< item profits */
   SCIP_Longint          capacity,           /**< capacity of knapsack */
   int*                  items,              /**< item numbers */
   int*                  solitems,           /**< array to store items in solution, or NULL */
   int*                  nonsolitems,        /**< array to store items not in solution, or NULL */
   int*                  nsolitems,          /**< pointer to store number of items in solution, or NULL */
   int*                  nnonsolitems,       /**< pointer to store number of items not in solution, or NULL */
   SCIP_Real*            solval              /**< pointer to store optimal solution value, or NULL */
   )
{
   SCIP_Real* tempsort;
   SCIP_Longint solitemsweight;
   SCIP_Real* realweights;
   int j;
   int criticalindex;

   assert(weights != NULL);
   assert(profits != NULL);
   assert(capacity >= 0);
   assert(items != NULL);
   assert(nitems >= 0);

   if( solitems != NULL )
   {
      *nsolitems = 0;
      *nnonsolitems = 0;
   }
   if( solval != NULL )
      *solval = 0.0;

   /* initialize data for median search */
   SCIP_CALL( SCIPallocBufferArray(scip, &tempsort, nitems) );
   SCIP_CALL( SCIPallocBufferArray(scip, &realweights, nitems) );
   for( j = nitems - 1; j >= 0; --j )
   {
      tempsort[j] = profits[j]/((SCIP_Real) weights[j]);
      realweights[j] = (SCIP_Real)weights[j];
   }

   /* partially sort indices such that all elements that are larger than the break item appear first */
   SCIPselectWeightedDownRealLongRealInt(tempsort, weights, profits, items, realweights, (SCIP_Real)capacity, nitems, &criticalindex);

   /* selects items as long as they fit into the knapsack */
   solitemsweight = 0;
   for( j = 0; j < nitems && solitemsweight + weights[j] <= capacity; ++j )
   {
      if( solitems != NULL )
         solitems[(*nsolitems)++] = items[j];

      if( solval != NULL )
         (*solval) += profits[j];
      solitemsweight += weights[j];
   }
   if ( solitems != NULL )
   {
      for( ; j < nitems; j++ )
         nonsolitems[(*nnonsolitems)++] = items[j];
   }

   SCIPfreeBufferArray(scip, &realweights);
   SCIPfreeBufferArray(scip, &tempsort);

   return SCIP_OKAY;
}

#ifdef SCIP_DEBUG
/** prints all nontrivial GUB constraints and their LP solution values */
static
void GUBsetPrint(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_GUBSET*          gubset,             /**< GUB set data structure */
   SCIP_VAR**            vars,               /**< variables in knapsack constraint */
   SCIP_Real*            solvals             /**< solution values of variables in knapsack constraint; or NULL */
   )
{
   int nnontrivialgubconss;
   int c;

   nnontrivialgubconss = 0;

   SCIPdebugMsg(scip, "   Nontrivial GUBs of current GUB set:\n");

   /* print out all nontrivial GUB constraints, i.e., with more than one variable */
   for( c = 0; c < gubset->ngubconss; c++ )
   {
      SCIP_Real gubsolval;

      assert(gubset->gubconss[c]->ngubvars >= 0);

      /* nontrivial GUB */
      if( gubset->gubconss[c]->ngubvars > 1 )
      {
         int v;

         gubsolval = 0.0;
         SCIPdebugMsg(scip, "   GUB<%d>:\n", c);

         /* print GUB var */
         for( v = 0; v < gubset->gubconss[c]->ngubvars; v++ )
         {
            int currentvar;

            currentvar = gubset->gubconss[c]->gubvars[v];
            if( solvals != NULL )
            {
               gubsolval += solvals[currentvar];
               SCIPdebugMsg(scip, "      +<%s>(%4.2f)\n", SCIPvarGetName(vars[currentvar]), solvals[currentvar]);
            }
            else
            {
               SCIPdebugMsg(scip, "      +<%s>\n", SCIPvarGetName(vars[currentvar]));
            }
         }

         /* check whether LP solution satisfies the GUB constraint */
         if( solvals != NULL )
         {
            SCIPdebugMsg(scip, "      =%4.2f <= 1 %s\n", gubsolval,
               SCIPisFeasGT(scip, gubsolval, 1.0) ? "--> violated" : "");
         }
         else
         {
            SCIPdebugMsg(scip, "      <= 1 %s\n", SCIPisFeasGT(scip, gubsolval, 1.0) ? "--> violated" : "");
         }
         nnontrivialgubconss++;
      }
   }

   SCIPdebugMsg(scip, "   --> %d/%d nontrivial GUBs\n", nnontrivialgubconss, gubset->ngubconss);
}
#endif

/** creates an empty GUB constraint */
static
SCIP_RETCODE GUBconsCreate(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_GUBCONS**        gubcons             /**< pointer to store GUB constraint data */
   )
{
   assert(scip != NULL);
   assert(gubcons != NULL);

   /* allocate memory for GUB constraint data structures */
   SCIP_CALL( SCIPallocBuffer(scip, gubcons) );
   (*gubcons)->gubvarssize = GUBCONSGROWVALUE;
   SCIP_CALL( SCIPallocBufferArray(scip, &(*gubcons)->gubvars, (*gubcons)->gubvarssize) );
   SCIP_CALL( SCIPallocBufferArray(scip, &(*gubcons)->gubvarsstatus, (*gubcons)->gubvarssize) );

   (*gubcons)->ngubvars = 0;

   return SCIP_OKAY;
}

/** frees GUB constraint */
static
void GUBconsFree(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_GUBCONS**        gubcons             /**< pointer to GUB constraint data structure */
   )
{
   assert(scip != NULL);
   assert(gubcons != NULL);
   assert((*gubcons)->gubvars != NULL);
   assert((*gubcons)->gubvarsstatus != NULL);

   /* free allocated memory */
   SCIPfreeBufferArray(scip, &(*gubcons)->gubvarsstatus);
   SCIPfreeBufferArray(scip, &(*gubcons)->gubvars);
   SCIPfreeBuffer(scip, gubcons);
}

/** adds variable to given GUB constraint */
static
SCIP_RETCODE GUBconsAddVar(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_GUBCONS*         gubcons,            /**< GUB constraint data */
   int                   var                 /**< index of given variable in knapsack constraint */
   )
{
   assert(scip != NULL);
   assert(gubcons != NULL);
   assert(gubcons->ngubvars >= 0 && gubcons->ngubvars < gubcons->gubvarssize);
   assert(gubcons->gubvars != NULL);
   assert(gubcons->gubvarsstatus != NULL);
   assert(var >= 0);

   /* add variable to GUB constraint */
   gubcons->gubvars[gubcons->ngubvars] = var;
   gubcons->gubvarsstatus[gubcons->ngubvars] = GUBVARSTATUS_UNINITIAL;
   gubcons->ngubvars++;

   /* increase space allocated to GUB constraint if the number of variables reaches the size */
   if( gubcons->ngubvars == gubcons->gubvarssize )
   {
      int newlen;

      newlen = gubcons->gubvarssize + GUBCONSGROWVALUE;
      SCIP_CALL( SCIPreallocBufferArray(scip, &gubcons->gubvars, newlen) );
      SCIP_CALL( SCIPreallocBufferArray(scip, &gubcons->gubvarsstatus, newlen) );

      gubcons->gubvarssize = newlen;
   }

   return SCIP_OKAY;
}

/** deletes variable from its current GUB constraint */
static
SCIP_RETCODE GUBconsDelVar(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_GUBCONS*         gubcons,            /**< GUB constraint data */
   int                   var,                /**< index of given variable in knapsack constraint */
   int                   gubvarsidx          /**< index of the variable in its current GUB constraint */
   )
{
   assert(scip != NULL);
   assert(gubcons != NULL);
   assert(var >= 0);
   assert(gubvarsidx >= 0 && gubvarsidx < gubcons->ngubvars);
   assert(gubcons->ngubvars >= gubvarsidx+1);
   assert(gubcons->gubvars[gubvarsidx] == var);

   /* delete variable from GUB by swapping it replacing in by the last variable in the GUB constraint */
   gubcons->gubvars[gubvarsidx] = gubcons->gubvars[gubcons->ngubvars-1];
   gubcons->gubvarsstatus[gubvarsidx] = gubcons->gubvarsstatus[gubcons->ngubvars-1];
   gubcons->ngubvars--;

   /* decrease space allocated for the GUB constraint, if the last GUBCONSGROWVALUE+1 array entries are now empty */
   if( gubcons->ngubvars < gubcons->gubvarssize - GUBCONSGROWVALUE && gubcons->ngubvars > 0 )
   {
      int newlen;

      newlen = gubcons->gubvarssize - GUBCONSGROWVALUE;

      SCIP_CALL( SCIPreallocBufferArray(scip, &gubcons->gubvars, newlen) );
      SCIP_CALL( SCIPreallocBufferArray(scip, &gubcons->gubvarsstatus, newlen) );

      gubcons->gubvarssize = newlen;
   }

   return SCIP_OKAY;
}

/** moves variable from current GUB constraint to a different existing (nonempty) GUB constraint */
static
SCIP_RETCODE GUBsetMoveVar(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_GUBSET*          gubset,             /**< GUB set data structure */
   SCIP_VAR**            vars,               /**< variables in knapsack constraint */
   int                   var,                /**< index of given variable in knapsack constraint */
   int                   oldgubcons,         /**< index of old GUB constraint of given variable */
   int                   newgubcons          /**< index of new GUB constraint of given variable */
   )
{
   int oldgubvaridx;
   int replacevar;
   int j;

   assert(scip != NULL);
   assert(gubset != NULL);
   assert(var >= 0);
   assert(oldgubcons >= 0 && oldgubcons < gubset->ngubconss);
   assert(newgubcons >= 0 && newgubcons < gubset->ngubconss);
   assert(oldgubcons != newgubcons);
   assert(gubset->gubconssidx[var] == oldgubcons);
   assert(gubset->gubconss[oldgubcons]->ngubvars > 0);
   assert(gubset->gubconss[newgubcons]->ngubvars >= 0);

   SCIPdebugMsg(scip, "   moving variable<%s> from GUB<%d> to GUB<%d>\n", SCIPvarGetName(vars[var]), oldgubcons, newgubcons);

   oldgubvaridx = gubset->gubvarsidx[var];

   /* delete variable from old GUB constraint by replacing it by the last variable of the GUB constraint */
   SCIP_CALL( GUBconsDelVar(scip, gubset->gubconss[oldgubcons], var, oldgubvaridx) );

   /* in GUB set, update stored index of variable in old GUB constraint for the variable used for replacement;
    * replacement variable is given by old position of the deleted variable
    */
   replacevar = gubset->gubconss[oldgubcons]->gubvars[oldgubvaridx];
   assert(gubset->gubvarsidx[replacevar] == gubset->gubconss[oldgubcons]->ngubvars);
   gubset->gubvarsidx[replacevar] = oldgubvaridx;

   /* add variable to the end of new GUB constraint */
   SCIP_CALL( GUBconsAddVar(scip, gubset->gubconss[newgubcons], var) );
   assert(gubset->gubconss[newgubcons]->gubvars[gubset->gubconss[newgubcons]->ngubvars-1] == var);

   /* in GUB set, update stored index of GUB of moved variable and stored index of variable in this GUB constraint */
   gubset->gubconssidx[var] = newgubcons;
   gubset->gubvarsidx[var] = gubset->gubconss[newgubcons]->ngubvars-1;

   /* delete old GUB constraint if it became empty */
   if( gubset->gubconss[oldgubcons]->ngubvars == 0 )
   {
      SCIPdebugMsg(scip, "deleting empty GUB cons<%d> from current GUB set\n", oldgubcons);
#ifdef SCIP_DEBUG
      GUBsetPrint(scip, gubset, vars, NULL);
#endif

      /* free old GUB constraint */
      GUBconsFree(scip, &gubset->gubconss[oldgubcons]);

      /* if empty GUB was not the last one in GUB set data structure, replace it by last GUB constraint */
      if( oldgubcons != gubset->ngubconss-1 )
      {
         gubset->gubconss[oldgubcons] = gubset->gubconss[gubset->ngubconss-1];
         gubset->gubconsstatus[oldgubcons] = gubset->gubconsstatus[gubset->ngubconss-1];

         /* in GUB set, update stored index of GUB constraint for all variable of the GUB constraint used for replacement;
          * replacement GUB is given by old position of the deleted GUB
          */
         for( j = 0; j < gubset->gubconss[oldgubcons]->ngubvars; j++ )
         {
            assert(gubset->gubconssidx[gubset->gubconss[oldgubcons]->gubvars[j]] == gubset->ngubconss-1);
            gubset->gubconssidx[gubset->gubconss[oldgubcons]->gubvars[j]] = oldgubcons;
         }
      }

      /* update number of GUB constraints */
      gubset->ngubconss--;

      /* variable should be at given new position, unless new GUB constraint replaced empty old GUB constraint
       * (because it was at the end of the GUB constraint array)
       */
      assert(gubset->gubconssidx[var] == newgubcons
         || (newgubcons == gubset->ngubconss && gubset->gubconssidx[var] == oldgubcons));
   }
#ifndef NDEBUG
   else
      assert(gubset->gubconssidx[var] == newgubcons);
#endif

   return SCIP_OKAY;
}

/** swaps two variables in the same GUB constraint */
static
void GUBsetSwapVars(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_GUBSET*          gubset,             /**< GUB set data structure */
   int                   var1,               /**< first variable to be swapped */
   int                   var2                /**< second variable to be swapped */
   )
{
   int gubcons;
   int var1idx;
   GUBVARSTATUS var1status;
   int var2idx;
   GUBVARSTATUS var2status;

   assert(scip != NULL);
   assert(gubset != NULL);

   gubcons = gubset->gubconssidx[var1];
   assert(gubcons == gubset->gubconssidx[var2]);

   /* nothing to be done if both variables are the same */
   if( var1 == var2 )
      return;

   /* swap index and status of variables in GUB constraint */
   var1idx = gubset->gubvarsidx[var1];
   var1status = gubset->gubconss[gubcons]->gubvarsstatus[var1idx];
   var2idx = gubset->gubvarsidx[var2];
   var2status = gubset->gubconss[gubcons]->gubvarsstatus[var2idx];

   gubset->gubvarsidx[var1] = var2idx;
   gubset->gubconss[gubcons]->gubvars[var1idx] = var2;
   gubset->gubconss[gubcons]->gubvarsstatus[var1idx] = var2status;

   gubset->gubvarsidx[var2] = var1idx;
   gubset->gubconss[gubcons]->gubvars[var2idx] = var1;
   gubset->gubconss[gubcons]->gubvarsstatus[var2idx] = var1status;
}

/** initializes partition of knapsack variables into nonoverlapping trivial GUB constraints (GUB with one variable) */
static
SCIP_RETCODE GUBsetCreate(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_GUBSET**         gubset,             /**< pointer to store GUB set data structure */
   int                   nvars,              /**< number of variables in the knapsack constraint */
   SCIP_Longint*         weights,            /**< weights of variables in knapsack constraint */
   SCIP_Longint          capacity            /**< capacity of knapsack */
   )
{
   int i;

   assert(scip != NULL);
   assert(gubset != NULL);
   assert(nvars > 0);
   assert(weights != NULL);
   assert(capacity >= 0);

   /* allocate memory for GUB set data structures */
   SCIP_CALL( SCIPallocBuffer(scip, gubset) );
   SCIP_CALL( SCIPallocBufferArray(scip, &(*gubset)->gubconss, nvars) );
   SCIP_CALL( SCIPallocBufferArray(scip, &(*gubset)->gubconsstatus, nvars) );
   SCIP_CALL( SCIPallocBufferArray(scip, &(*gubset)->gubconssidx, nvars) );
   SCIP_CALL( SCIPallocBufferArray(scip, &(*gubset)->gubvarsidx, nvars) );
   (*gubset)->ngubconss = nvars;
   (*gubset)->nvars = nvars;

   /* initialize the set of GUB constraints */
   for( i = 0; i < nvars; i++ )
   {
      /* assign each variable to a new (trivial) GUB constraint */
      SCIP_CALL( GUBconsCreate(scip, &(*gubset)->gubconss[i]) );
      SCIP_CALL( GUBconsAddVar(scip, (*gubset)->gubconss[i], i) );

      /* set status of GUB constraint to initial */
      (*gubset)->gubconsstatus[i] = GUBCONSSTATUS_UNINITIAL;

      (*gubset)->gubconssidx[i] = i;
      (*gubset)->gubvarsidx[i] = 0;
      assert((*gubset)->gubconss[i]->ngubvars == 1);

      /* already updated status of variable in GUB constraint if it exceeds the capacity of the knapsack */
      if( weights[i] > capacity )
         (*gubset)->gubconss[(*gubset)->gubconssidx[i]]->gubvarsstatus[(*gubset)->gubvarsidx[i]] = GUBVARSTATUS_CAPACITYEXCEEDED;
   }

   return SCIP_OKAY;
}

/** frees GUB set data structure */
static
void GUBsetFree(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_GUBSET**         gubset              /**< pointer to GUB set data structure */
   )
{
   int i;

   assert(scip != NULL);
   assert(gubset != NULL);
   assert((*gubset)->gubconss != NULL);
   assert((*gubset)->gubconsstatus != NULL);
   assert((*gubset)->gubconssidx != NULL);
   assert((*gubset)->gubvarsidx != NULL);

   /* free all GUB constraints */
   for( i = (*gubset)->ngubconss-1; i >= 0; --i )
   {
      assert((*gubset)->gubconss[i] != NULL);
      GUBconsFree(scip, &(*gubset)->gubconss[i]);
   }

   /* free allocated memory */
   SCIPfreeBufferArray( scip, &(*gubset)->gubvarsidx );
   SCIPfreeBufferArray( scip, &(*gubset)->gubconssidx );
   SCIPfreeBufferArray( scip, &(*gubset)->gubconsstatus );
   SCIPfreeBufferArray( scip, &(*gubset)->gubconss );
   SCIPfreeBuffer(scip, gubset);
}

#ifndef NDEBUG
/** checks whether GUB set data structure is consistent */
static
SCIP_RETCODE GUBsetCheck(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_GUBSET*          gubset,             /**< GUB set data structure */
   SCIP_VAR**            vars                /**< variables in the knapsack constraint */
   )
{
   int i;
   int gubconsidx;
   int gubvaridx;
   SCIP_VAR* var1;
   SCIP_VAR* var2;
   SCIP_Bool var1negated;
   SCIP_Bool var2negated;

   assert(scip != NULL);
   assert(gubset != NULL);

   SCIPdebugMsg(scip, "   GUB set consistency check:\n");

   /* checks for all knapsack vars consistency of stored index of associated gubcons and corresponding index in gubvars */
   for( i = 0; i < gubset->nvars; i++ )
   {
      gubconsidx = gubset->gubconssidx[i];
      gubvaridx = gubset->gubvarsidx[i];

      if( gubset->gubconss[gubconsidx]->gubvars[gubvaridx] != i )
      {
         SCIPdebugMsg(scip, "   var<%d> should be in GUB<%d> at position<%d>, but stored is var<%d> instead\n", i,
            gubconsidx, gubvaridx, gubset->gubconss[gubconsidx]->gubvars[gubvaridx] );
      }
      assert(gubset->gubconss[gubconsidx]->gubvars[gubvaridx] == i);
   }

   /* checks for each GUB whether all pairs of its variables have a common clique */
   for( i = 0; i < gubset->ngubconss; i++ )
   {
      int j;

      for( j = 0; j < gubset->gubconss[i]->ngubvars; j++ )
      {
         int k;

         /* get corresponding active problem variable */
         var1 = vars[gubset->gubconss[i]->gubvars[j]];
         var1negated = FALSE;
         SCIP_CALL( SCIPvarGetProbvarBinary(&var1, &var1negated) );

         for( k = j+1; k < gubset->gubconss[i]->ngubvars; k++ )
         {
            /* get corresponding active problem variable */
            var2 = vars[gubset->gubconss[i]->gubvars[k]];
            var2negated = FALSE;
            SCIP_CALL( SCIPvarGetProbvarBinary(&var2, &var2negated) );

            if( !SCIPvarsHaveCommonClique(var1, !var1negated, var2, !var2negated, TRUE) )
            {
               SCIPdebugMsg(scip, "   GUB<%d>: var<%d,%s> and var<%d,%s> do not share a clique\n", i, j,
                  SCIPvarGetName(vars[gubset->gubconss[i]->gubvars[j]]), k,
                  SCIPvarGetName(vars[gubset->gubconss[i]->gubvars[k]]));
               SCIPdebugMsg(scip, "   GUB<%d>: var<%d,%s> and var<%d,%s> do not share a clique\n", i, j,
                  SCIPvarGetName(var1), k,
                  SCIPvarGetName(var2));
            }

            /* @todo: in case we used also negated cliques for the GUB partition, this assert has to be changed */
            assert(SCIPvarsHaveCommonClique(var1, !var1negated, var2, !var2negated, TRUE));
         }
      }
   }
   SCIPdebugMsg(scip, "   --> successful\n");

   return SCIP_OKAY;
}
#endif

/** calculates a partition of the given set of binary variables into cliques;
 *  afterwards the output array contains one value for each variable, such that two variables got the same value iff they
 *  were assigned to the same clique;
 *  the first variable is always assigned to clique 0, and a variable can only be assigned to clique i if at least one of
 *  the preceding variables was assigned to clique i-1;
 *  note: in contrast to SCIPcalcCliquePartition(), variables with LP value 1 are put into trivial cliques (with one
 *  variable) and for the remaining variables, a partition with a small number of cliques is constructed
 */

static
SCIP_RETCODE GUBsetCalcCliquePartition(
   SCIP*const            scip,               /**< SCIP data structure */
   SCIP_VAR**const       vars,               /**< binary variables in the clique from which at most one can be set to 1 */
   int const             nvars,              /**< number of variables in the clique */
   int*const             cliquepartition,    /**< array of length nvars to store the clique partition */
   int*const             ncliques,           /**< pointer to store number of cliques actually contained in the partition */
   SCIP_Real*            solvals             /**< solution values of all given binary variables */
   )
{
   SCIP_VAR** tmpvars;
   SCIP_VAR** cliquevars;
   SCIP_Bool* cliquevalues;
   SCIP_Bool* tmpvalues;
   int* varseq;
   int* sortkeys;
   int ncliquevars;
   int maxncliquevarscomp;
   int nignorevars;
   int nvarsused;
   int i;

   assert(scip != NULL);
   assert(nvars == 0 || vars != NULL);
   assert(nvars == 0 || cliquepartition != NULL);
   assert(ncliques != NULL);

   if( nvars == 0 )
   {
      *ncliques = 0;
      return SCIP_OKAY;
   }

   /* allocate temporary memory for storing the variables of the current clique */
   SCIP_CALL( SCIPallocBufferArray(scip, &cliquevars, nvars) );
   SCIP_CALL( SCIPallocBufferArray(scip, &cliquevalues, nvars) );
   SCIP_CALL( SCIPallocBufferArray(scip, &tmpvalues, nvars) );
   SCIP_CALL( SCIPduplicateBufferArray(scip, &tmpvars, vars, nvars) );
   SCIP_CALL( SCIPallocBufferArray(scip, &varseq, nvars) );
   SCIP_CALL( SCIPallocBufferArray(scip, &sortkeys, nvars) );

   /* initialize the cliquepartition array with -1 */
   /* initialize the tmpvalues array */
   for( i = nvars - 1; i >= 0; --i )
   {
      tmpvalues[i] = TRUE;
      cliquepartition[i] = -1;
   }

   /* get corresponding active problem variables */
   SCIP_CALL( SCIPvarsGetProbvarBinary(&tmpvars, &tmpvalues, nvars) );

   /* ignore variables with LP value 1 (will be assigned to trivial GUBs at the end) and sort remaining variables
    * by nondecreasing number of cliques the variables are in
    */
   nignorevars = 0;
   nvarsused = 0;
   for( i = 0; i < nvars; i++ )
   {
      if( SCIPisFeasEQ(scip, solvals[i], 1.0) )
      {
         /* variables with LP value 1 are put to the end of varseq array and will not be sorted */
         varseq[nvars-1-nignorevars] = i;
         nignorevars++;
      }
      else
      {
         /* remaining variables are put to the front of varseq array and will be sorted by their number of cliques */
         varseq[nvarsused] = i;
         sortkeys[nvarsused] = SCIPvarGetNCliques(tmpvars[i], tmpvalues[i]);
         nvarsused++;
      }
   }
   assert(nvarsused + nignorevars == nvars);

   /* sort variables with LP value less than 1 by nondecreasing order of the number of cliques they are in */
   SCIPsortIntInt(sortkeys, varseq, nvarsused);

   maxncliquevarscomp = MIN(nvars*nvars, MAXNCLIQUEVARSCOMP);

   /* calculate the clique partition */
   *ncliques = 0;
   for( i = 0; i < nvars; ++i )
   {
      if( cliquepartition[varseq[i]] == -1 )
      {
         int j;

         /* variable starts a new clique */
         cliquepartition[varseq[i]] = *ncliques;
         cliquevars[0] = tmpvars[varseq[i]];
         cliquevalues[0] = tmpvalues[varseq[i]];
         ncliquevars = 1;

         /* if variable is not active (multi-aggregated or fixed), it cannot be in any clique and
          * if the variable has LP value 1 we do not want it to be in nontrivial cliques
          */
         if( SCIPvarIsActive(tmpvars[varseq[i]]) && i < nvarsused )
         {
            /* greedily fill up the clique */
            for( j = i + 1; j < nvarsused; ++j )
            {
               /* if variable is not active (multi-aggregated or fixed), it cannot be in any clique */
               if( cliquepartition[varseq[j]] == -1 && SCIPvarIsActive(tmpvars[varseq[j]]) )
               {
                  int k;

                  /* check if every variable in the actual clique is in clique with the new variable */
                  for( k = ncliquevars - 1; k >= 0; --k )
                  {
                     if( !SCIPvarsHaveCommonClique(tmpvars[varseq[j]], tmpvalues[varseq[j]], cliquevars[k],
                           cliquevalues[k], TRUE) )
                        break;
                  }

                  if( k == -1 )
                  {
                     /* put the variable into the same clique */
                     cliquepartition[varseq[j]] = cliquepartition[varseq[i]];
                     cliquevars[ncliquevars] = tmpvars[varseq[j]];
                     cliquevalues[ncliquevars] = tmpvalues[varseq[j]];
                     ++ncliquevars;
                  }
               }
            }
         }

         /* this clique is finished */
         ++(*ncliques);
      }
      assert(cliquepartition[varseq[i]] >= 0 && cliquepartition[varseq[i]] < i + 1);

      /* break if we reached the maximal number of comparisons */
      if( i * nvars > maxncliquevarscomp )
         break;
   }
   /* if we had too many variables fill up the cliquepartition and put each variable in a separate clique */
   for( ; i < nvars; ++i )
   {
      if( cliquepartition[varseq[i]] == -1 )
      {
         cliquepartition[varseq[i]] = *ncliques;
         ++(*ncliques);
      }
   }

   /* free temporary memory */
   SCIPfreeBufferArray(scip, &sortkeys);
   SCIPfreeBufferArray(scip, &varseq);
   SCIPfreeBufferArray(scip, &tmpvars);
   SCIPfreeBufferArray(scip, &tmpvalues);
   SCIPfreeBufferArray(scip, &cliquevalues);
   SCIPfreeBufferArray(scip, &cliquevars);

   return SCIP_OKAY;
}

/** constructs sophisticated partition of knapsack variables into non-overlapping GUBs; current partition uses trivial GUBs */
static
SCIP_RETCODE GUBsetGetCliquePartition(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_GUBSET*          gubset,             /**< GUB set data structure */
   SCIP_VAR**            vars,               /**< variables in the knapsack constraint */
   SCIP_Real*            solvals             /**< solution values of all knapsack variables */
   )
{
   int* cliquepartition;
   int* gubfirstvar;
   int ncliques;
   int currentgubconsidx;
   int newgubconsidx;
   int cliqueidx;
   int nvars;
   int i;

   assert(scip != NULL);
   assert(gubset != NULL);
   assert(vars != NULL);

   nvars = gubset->nvars;
   assert(nvars >= 0);

   /* allocate temporary memory for clique partition */
   SCIP_CALL( SCIPallocBufferArray(scip, &cliquepartition, nvars) );

   /* compute sophisticated clique partition */
   SCIP_CALL( GUBsetCalcCliquePartition(scip, vars, nvars, cliquepartition, &ncliques, solvals) );

   /* allocate temporary memory for GUB set data structure */
   SCIP_CALL( SCIPallocBufferArray(scip, &gubfirstvar, ncliques) );

   /* translate GUB partition into GUB set data structure */
   for( i = 0; i < ncliques; i++ )
   {
      /* initialize first variable for every GUB */
      gubfirstvar[i] = -1;
   }
   /* move every knapsack variable into GUB defined by clique partition */
   for( i = 0; i < nvars; i++ )
   {
      assert(cliquepartition[i] >= 0);

      cliqueidx = cliquepartition[i];
      currentgubconsidx = gubset->gubconssidx[i];
      assert(gubset->gubconss[currentgubconsidx]->ngubvars == 1 );

      /* variable is first element in GUB constraint defined by clique partition */
      if( gubfirstvar[cliqueidx] == -1 )
      {
         /* corresponding GUB constraint in GUB set data structure was already constructed (as initial trivial GUB);
          * note: no assert for gubconssidx, because it can changed due to deleting empty GUBs in GUBsetMoveVar()
          */
         assert(gubset->gubvarsidx[i] == 0);
         assert(gubset->gubconss[gubset->gubconssidx[i]]->gubvars[gubset->gubvarsidx[i]] == i);

         /* remember the first variable found for the current GUB */
         gubfirstvar[cliqueidx] = i;
      }
      /* variable is additional element of GUB constraint defined by clique partition */
      else
      {
         assert(gubfirstvar[cliqueidx] >= 0 && gubfirstvar[cliqueidx] < i);

         /* move variable to GUB constraint defined by clique partition; index of this GUB constraint is given by the
          * first variable of this GUB constraint
          */
         newgubconsidx = gubset->gubconssidx[gubfirstvar[cliqueidx]];
         assert(newgubconsidx != currentgubconsidx); /* because initially every variable is in a different GUB */
         SCIP_CALL( GUBsetMoveVar(scip, gubset, vars, i, currentgubconsidx, newgubconsidx) );

         assert(gubset->gubconss[gubset->gubconssidx[i]]->gubvars[gubset->gubvarsidx[i]] == i);
      }
   }

#ifdef SCIP_DEBUG
   /* prints GUB set data structure */
   GUBsetPrint(scip, gubset, vars, solvals);
#endif

#ifndef NDEBUG
   /* checks consistency of GUB set data structure */
   SCIP_CALL( GUBsetCheck(scip, gubset, vars) );
#endif

   /* free temporary memory */
   SCIPfreeBufferArray(scip, &gubfirstvar);
   SCIPfreeBufferArray(scip, &cliquepartition);

   return SCIP_OKAY;
}

/** gets a most violated cover C (\f$\sum_{j \in C} a_j > a_0\f$) for a given knapsack constraint \f$\sum_{j \in N} a_j x_j \leq a_0\f$
 *  taking into consideration the following fixing: \f$j \in C\f$, if \f$j \in N_1 = \{j \in N : x^*_j = 1\}\f$ and
 *  \f$j \in N \setminus C\f$, if \f$j \in N_0 = \{j \in N : x^*_j = 0\}\f$, if one exists.
 */
static
SCIP_RETCODE getCover(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_VAR**            vars,               /**< variables in knapsack constraint */
   int                   nvars,              /**< number of variables in knapsack constraint */
   SCIP_Longint*         weights,            /**< weights of variables in knapsack constraint */
   SCIP_Longint          capacity,           /**< capacity of knapsack */
   SCIP_Real*            solvals,            /**< solution values of all problem variables */
   int*                  covervars,          /**< pointer to store cover variables */
   int*                  noncovervars,       /**< pointer to store noncover variables */
   int*                  ncovervars,         /**< pointer to store number of cover variables */
   int*                  nnoncovervars,      /**< pointer to store number of noncover variables */
   SCIP_Longint*         coverweight,        /**< pointer to store weight of cover */
   SCIP_Bool*            found,              /**< pointer to store whether a cover was found */
   SCIP_Bool             modtransused,       /**< should modified transformed separation problem be used to find cover */
   int*                  ntightened,         /**< pointer to store number of variables with tightened upper bound */
   SCIP_Bool*            fractional          /**< pointer to store whether the LP sol for knapsack vars is fractional */
   )
{
   SCIP_Longint* transweights;
   SCIP_Real* transprofits;
   SCIP_Longint transcapacity;
   SCIP_Longint fixedonesweight;
   SCIP_Longint itemsweight;
   SCIP_Bool infeasible;
   int* fixedones;
   int* fixedzeros;
   int* items;
   int nfixedones;
   int nfixedzeros;
   int nitems;
   int j;

   assert(scip != NULL);
   assert(vars != NULL);
   assert(nvars > 0);
   assert(weights != NULL);
   assert(capacity >= 0);
   assert(solvals != NULL);
   assert(covervars != NULL);
   assert(noncovervars != NULL);
   assert(ncovervars != NULL);
   assert(nnoncovervars != NULL);
   assert(coverweight != NULL);
   assert(found != NULL);
   assert(ntightened != NULL);
   assert(fractional != NULL);

   SCIPdebugMsg(scip, "   get cover for knapsack constraint\n");

   /* allocates temporary memory */
   SCIP_CALL( SCIPallocBufferArray(scip, &transweights, nvars) );
   SCIP_CALL( SCIPallocBufferArray(scip, &transprofits, nvars) );
   SCIP_CALL( SCIPallocBufferArray(scip, &fixedones, nvars) );
   SCIP_CALL( SCIPallocBufferArray(scip, &fixedzeros, nvars) );
   SCIP_CALL( SCIPallocBufferArray(scip, &items, nvars) );

   *found = FALSE;
   *ncovervars = 0;
   *nnoncovervars = 0;
   *coverweight = 0;
   *fractional = TRUE;

   /* gets the following sets
    *  N_1 = {j in N : x*_j = 1} (fixedones),
    *  N_0 = {j in N : x*_j = 0} (fixedzeros) and
    *  N\(N_0 & N_1) (items),
    * where x*_j is the solution value of variable x_j
    */
   nfixedones = 0;
   nfixedzeros = 0;
   nitems = 0;
   fixedonesweight = 0;
   itemsweight = 0;
   *ntightened = 0;
   for( j = 0; j < nvars; j++ )
   {
      assert(SCIPvarIsBinary(vars[j]));

      /* tightens upper bound of x_j if weight of x_j is greater than capacity of knapsack */
      if( weights[j] > capacity )
      {
         SCIP_CALL( SCIPtightenVarUb(scip, vars[j], 0.0, FALSE, &infeasible, NULL) );
         assert(!infeasible);
         (*ntightened)++;
         continue;
      }

      /* variable x_j has solution value one */
      if( SCIPisFeasEQ(scip, solvals[j], 1.0) )
      {
         fixedones[nfixedones] = j;
         nfixedones++;
         fixedonesweight += weights[j];
      }
      /* variable x_j has solution value zero */
      else if( SCIPisFeasEQ(scip, solvals[j], 0.0) )
      {
         fixedzeros[nfixedzeros] = j;
         nfixedzeros++;
      }
      /* variable x_j has fractional solution value */
      else
      {
         assert( SCIPisFeasGT(scip, solvals[j], 0.0) && SCIPisFeasLT(scip, solvals[j], 1.0) );
         items[nitems] = j;
         nitems++;
         itemsweight += weights[j];
      }
   }
   assert(nfixedones + nfixedzeros + nitems == nvars - (*ntightened));

   /* sets whether the LP solution x* for the knapsack variables is fractional; if it is not fractional we stop
    * the separation routine
    */
   assert(nitems >= 0);
   if( nitems == 0 )
   {
      *fractional = FALSE;
      goto TERMINATE;
   }
   assert(*fractional);

   /* transforms the traditional separation problem (under consideration of the following fixing:
    * z_j = 1 for all j in N_1, z_j = 0 for all j in N_0)
    *
    *   min sum_{j in N\(N_0 & N_1)} (1 - x*_j) z_j
    *       sum_{j in N\(N_0 & N_1)} a_j z_j >= (a_0 + 1) - sum_{j in N_1} a_j
    *                                    z_j in {0,1}, j in N\(N_0 & N_1)
    *
    * to a knapsack problem in maximization form by complementing the variables
    *
    * sum_{j in N\(N_0 & N_1)} (1 - x*_j) -
    *   max sum_{j in N\(N_0 & N_1)} (1 - x*_j) z_j
    *       sum_{j in N\(N_0 & N_1)} a_j z_j <= sum_{j in N\N_0} a_j - (a_0 + 1)
    *                                    z_j in {0,1}, j in N\(N_0 & N_1)
    */

   /* gets weight and profit of variables in transformed knapsack problem */
   for( j = 0; j < nitems; j++ )
   {
      transweights[j] = weights[items[j]];
      transprofits[j] = 1.0 - solvals[items[j]];
   }
   /* gets capacity of transformed knapsack problem */
   transcapacity = fixedonesweight + itemsweight - capacity - 1;

   /* if capacity of transformed knapsack problem is less than zero, there is no cover
    * (when variables fixed to zero are not used)
    */
   if( transcapacity < 0 )
   {
      assert(!(*found));
      goto TERMINATE;
   }

   if( modtransused )
   {
      /* transforms the modified separation problem (under consideration of the following fixing:
       * z_j = 1 for all j in N_1, z_j = 0 for all j in N_0)
       *
       *   min sum_{j in N\(N_0 & N_1)} (1 - x*_j) a_j z_j
       *       sum_{j in N\(N_0 & N_1)} a_j z_j >= (a_0 + 1) - sum_{j in N_1} a_j
       *                                    z_j in {0,1}, j in N\(N_0 & N_1)
       *
       * to a knapsack problem in maximization form by complementing the variables
       *
       * sum_{j in N\(N_0 & N_1)} (1 - x*_j) a_j -
       *   max sum_{j in N\(N_0 & N_1)} (1 - x*_j) a_j z_j
       *       sum_{j in N\(N_0 & N_1)} a_j z_j <= sum_{j in N\N_0} a_j - (a_0 + 1)
       *                                    z_j in {0,1}, j in N\(N_0 & N_1)
       */

      /* gets weight and profit of variables in modified transformed knapsack problem */
      for( j = 0; j < nitems; j++ )
      {
         transprofits[j] *= weights[items[j]];
         assert(SCIPisFeasPositive(scip, transprofits[j]));
      }
   }

   /* solves (modified) transformed knapsack problem approximately by solving the LP-relaxation of the (modified)
    * transformed knapsack problem using Dantzig's method and rounding down the solution.
    * let z* be the solution, then
    *   j in C,          if z*_j = 0 and
    *   i in N\C,        if z*_j = 1.
    */
   SCIP_CALL( SCIPsolveKnapsackApproximately(scip, nitems, transweights, transprofits, transcapacity, items,
         noncovervars, covervars, nnoncovervars, ncovervars, NULL) );
   /*assert(checkSolveKnapsack(scip, nitems, transweights, transprofits, items, weights, solvals, modtransused));*/

   /* constructs cover C (sum_{j in C} a_j > a_0) */
   for( j = 0; j < *ncovervars; j++ )
   {
      (*coverweight) += weights[covervars[j]];
   }

   /* adds all variables from N_1 to C */
   for( j = 0; j < nfixedones; j++ )
   {
      covervars[*ncovervars] = fixedones[j];
      (*ncovervars)++;
      (*coverweight) += weights[fixedones[j]];
   }

   /* adds all variables from N_0 to N\C */
   for( j = 0; j < nfixedzeros; j++ )
   {
      noncovervars[*nnoncovervars] = fixedzeros[j];
      (*nnoncovervars)++;
   }
   assert((*ncovervars) + (*nnoncovervars) == nvars - (*ntightened));
   assert((*coverweight) > capacity);
   *found = TRUE;

 TERMINATE:
   /* frees temporary memory */
   SCIPfreeBufferArray(scip, &items);
   SCIPfreeBufferArray(scip, &fixedzeros);
   SCIPfreeBufferArray(scip, &fixedones);
   SCIPfreeBufferArray(scip, &transprofits);
   SCIPfreeBufferArray(scip, &transweights);

   SCIPdebugMsg(scip, "   get cover for knapsack constraint -- end\n");

   return SCIP_OKAY;
}

#ifndef NDEBUG
/** checks if minweightidx is set correctly
 */
static
SCIP_Bool checkMinweightidx(
   SCIP_Longint*         weights,            /**< weights of variables in knapsack constraint */
   SCIP_Longint          capacity,           /**< capacity of knapsack */
   int*                  covervars,          /**< pointer to store cover variables */
   int                   ncovervars,         /**< pointer to store number of cover variables */
   SCIP_Longint          coverweight,        /**< pointer to store weight of cover */
   int                   minweightidx,       /**< index of variable in cover variables with minimum weight */
   int                   j                   /**< current index in cover variables */
   )
{
   SCIP_Longint minweight;
   int i;

   assert(weights != NULL);
   assert(covervars != NULL);
   assert(ncovervars > 0);

   minweight = weights[covervars[minweightidx]];

   /* checks if all cover variables before index j have weight greater than minweight */
   for( i = 0; i < j; i++ )
   {
      assert(weights[covervars[i]] > minweight);
      if( weights[covervars[i]] <= minweight )
         return FALSE;
   }

   /* checks if all variables before index j cannot be removed, i.e. i cannot be the next minweightidx */
   for( i = 0; i < j; i++ )
   {
      assert(coverweight - weights[covervars[i]] <= capacity);
      if( coverweight - weights[covervars[i]] > capacity )
         return FALSE;
   }
   return TRUE;
}
#endif


/** gets partition \f$(C_1,C_2)\f$ of minimal cover \f$C\f$, i.e. \f$C_1 \cup C_2 = C\f$ and \f$C_1 \cap C_2 = \emptyset\f$,
 *  with \f$C_1\f$ not empty; chooses partition as follows \f$C_2 = \{ j \in C : x^*_j = 1 \}\f$ and \f$C_1 = C \setminus C_2\f$
 */
static
void getPartitionCovervars(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_Real*            solvals,            /**< solution values of all problem variables */
   int*                  covervars,          /**< cover variables */
   int                   ncovervars,         /**< number of cover variables */
   int*                  varsC1,             /**< pointer to store variables in C1 */
   int*                  varsC2,             /**< pointer to store variables in C2 */
   int*                  nvarsC1,            /**< pointer to store number of variables in C1 */
   int*                  nvarsC2             /**< pointer to store number of variables in C2 */
   )
{
   int j;

   assert(scip != NULL);
   assert(ncovervars >= 0);
   assert(solvals != NULL);
   assert(covervars != NULL);
   assert(varsC1 != NULL);
   assert(varsC2 != NULL);
   assert(nvarsC1 != NULL);
   assert(nvarsC2 != NULL);

   *nvarsC1 = 0;
   *nvarsC2 = 0;
   for( j = 0; j < ncovervars; j++ )
   {
      assert(SCIPisFeasGT(scip, solvals[covervars[j]], 0.0));

      /* variable has solution value one */
      if( SCIPisGE(scip, solvals[covervars[j]], 1.0) )
      {
         varsC2[*nvarsC2] = covervars[j];
         (*nvarsC2)++;
      }
      /* variable has solution value less than one */
      else
      {
         assert(SCIPisLT(scip, solvals[covervars[j]], 1.0));
         varsC1[*nvarsC1] = covervars[j];
         (*nvarsC1)++;
      }
   }
   assert((*nvarsC1) + (*nvarsC2) == ncovervars);
}

/** changes given partition (C_1,C_2) of minimal cover C, if |C1| = 1, by moving one and two (if possible) variables from
 *  C2 to C1 if |C1| = 1 and |C1| = 0, respectively.
 */
static
SCIP_RETCODE changePartitionCovervars(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_Longint*         weights,            /**< weights of variables in knapsack constraint */
   int*                  varsC1,             /**< pointer to store variables in C1 */
   int*                  varsC2,             /**< pointer to store variables in C2 */
   int*                  nvarsC1,            /**< pointer to store number of variables in C1 */
   int*                  nvarsC2             /**< pointer to store number of variables in C2 */
   )
{
   SCIP_Real* sortkeysC2;
   int j;

   assert(*nvarsC1 >= 0 && *nvarsC1 <= 1);
   assert(*nvarsC2 > 0);

   /* allocates temporary memory */
   SCIP_CALL( SCIPallocBufferArray(scip, &sortkeysC2, *nvarsC2) );

   /* sorts variables in C2 such that a_1 >= .... >= a_|C2| */
   for( j = 0; j < *nvarsC2; j++ )
      sortkeysC2[j] = (SCIP_Real) weights[varsC2[j]];
   SCIPsortDownRealInt(sortkeysC2, varsC2, *nvarsC2);

   /* adds one or two variable from C2 with smallest weight to C1 and removes them from C2 */
   assert(*nvarsC2 == 1 || weights[varsC2[(*nvarsC2)-1]] <= weights[varsC2[(*nvarsC2)-2]]);
   while( *nvarsC1 < 2 && *nvarsC2 > 0 )
   {
      varsC1[*nvarsC1] = varsC2[(*nvarsC2)-1];
      (*nvarsC1)++;
      (*nvarsC2)--;
   }

   /* frees temporary memory */
   SCIPfreeBufferArray(scip, &sortkeysC2);

   return SCIP_OKAY;
}

/** changes given partition (C_1,C_2) of feasible set C, if |C1| = 1, by moving one variable from C2 to C1 */
static
SCIP_RETCODE changePartitionFeasiblesetvars(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_Longint*         weights,            /**< weights of variables in knapsack constraint */
   int*                  varsC1,             /**< pointer to store variables in C1 */
   int*                  varsC2,             /**< pointer to store variables in C2 */
   int*                  nvarsC1,            /**< pointer to store number of variables in C1 */
   int*                  nvarsC2             /**< pointer to store number of variables in C2 */
   )
{
   SCIP_Real* sortkeysC2;
   int j;

   assert(*nvarsC1 >= 0 && *nvarsC1 <= 1);
   assert(*nvarsC2 > 0);

   /* allocates temporary memory */
   SCIP_CALL( SCIPallocBufferArray(scip, &sortkeysC2, *nvarsC2) );

   /* sorts variables in C2 such that a_1 >= .... >= a_|C2| */
   for( j = 0; j < *nvarsC2; j++ )
      sortkeysC2[j] = (SCIP_Real) weights[varsC2[j]];
   SCIPsortDownRealInt(sortkeysC2, varsC2, *nvarsC2);

   /* adds variable from C2 with smallest weight to C1 and removes it from C2 */
   assert(*nvarsC2 == 1 || weights[varsC2[(*nvarsC2)-1]] <= weights[varsC2[(*nvarsC2)-2]]);
   varsC1[*nvarsC1] = varsC2[(*nvarsC2)-1];
   (*nvarsC1)++;
   (*nvarsC2)--;

   /* frees temporary memory */
   SCIPfreeBufferArray(scip, &sortkeysC2);

   return SCIP_OKAY;
}


/** gets partition \f$(F,R)\f$ of \f$N \setminus C\f$ where \f$C\f$ is a minimal cover, i.e. \f$F \cup R = N \setminus C\f$
 *  and \f$F \cap R = \emptyset\f$; chooses partition as follows \f$R = \{ j \in N \setminus C : x^*_j = 0 \}\f$ and
 *  \f$F = (N \setminus C) \setminus F\f$
 */
static
void getPartitionNoncovervars(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_Real*            solvals,            /**< solution values of all problem variables */
   int*                  noncovervars,       /**< noncover variables */
   int                   nnoncovervars,      /**< number of noncover variables */
   int*                  varsF,              /**< pointer to store variables in F */
   int*                  varsR,              /**< pointer to store variables in R */
   int*                  nvarsF,             /**< pointer to store number of variables in F */
   int*                  nvarsR              /**< pointer to store number of variables in R */
   )
{
   int j;

   assert(scip != NULL);
   assert(nnoncovervars >= 0);
   assert(solvals != NULL);
   assert(noncovervars != NULL);
   assert(varsF != NULL);
   assert(varsR != NULL);
   assert(nvarsF != NULL);
   assert(nvarsR != NULL);

   *nvarsF = 0;
   *nvarsR = 0;

   for( j = 0; j < nnoncovervars; j++ )
   {
      /* variable has solution value zero */
      if( SCIPisFeasEQ(scip, solvals[noncovervars[j]], 0.0) )
      {
         varsR[*nvarsR] = noncovervars[j];
         (*nvarsR)++;
      }
      /* variable has solution value greater than zero */
      else
      {
         assert(SCIPisFeasGT(scip, solvals[noncovervars[j]], 0.0));
         varsF[*nvarsF] = noncovervars[j];
         (*nvarsF)++;
      }
   }
   assert((*nvarsF) + (*nvarsR) == nnoncovervars);
}

/** sorts variables in F, C_2, and R according to the second level lifting sequence that will be used in the sequential
 *  lifting procedure
 */
static
SCIP_RETCODE getLiftingSequence(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_Real*            solvals,            /**< solution values of all problem variables */
   SCIP_Longint*         weights,            /**< weights of variables in knapsack constraint */
   int*                  varsF,              /**< pointer to store variables in F */
   int*                  varsC2,             /**< pointer to store variables in C2 */
   int*                  varsR,              /**< pointer to store variables in R */
   int                   nvarsF,             /**< number of variables in F */
   int                   nvarsC2,            /**< number of variables in C2 */
   int                   nvarsR              /**< number of variables in R */
   )
{
   SORTKEYPAIR** sortkeypairsF;
   SORTKEYPAIR* sortkeypairsFstore;
   SCIP_Real* sortkeysC2;
   SCIP_Real* sortkeysR;
   int j;

   assert(scip != NULL);
   assert(solvals != NULL);
   assert(weights != NULL);
   assert(varsF != NULL);
   assert(varsC2 != NULL);
   assert(varsR != NULL);
   assert(nvarsF >= 0);
   assert(nvarsC2 >= 0);
   assert(nvarsR >= 0);

   /* allocates temporary memory */
   SCIP_CALL( SCIPallocBufferArray(scip, &sortkeypairsF, nvarsF) );
   SCIP_CALL( SCIPallocBufferArray(scip, &sortkeypairsFstore, nvarsF) );
   SCIP_CALL( SCIPallocBufferArray(scip, &sortkeysC2, nvarsC2) );
   SCIP_CALL( SCIPallocBufferArray(scip, &sortkeysR, nvarsR) );

   /* gets sorting key for variables in F corresponding to the following lifting sequence
    *  sequence 1: non-increasing absolute difference between x*_j and the value the variable is fixed to, i.e.
    *              x*_1 >= x*_2 >= ... >= x*_|F|
    * in case of equality uses
    *  sequence 4: non-increasing a_j, i.e. a_1 >= a_2 >= ... >= a_|C_2|
    */
   for( j = 0; j < nvarsF; j++ )
   {
      sortkeypairsF[j] = &(sortkeypairsFstore[j]);
      sortkeypairsF[j]->key1 = solvals[varsF[j]];
      sortkeypairsF[j]->key2 = (SCIP_Real) weights[varsF[j]];
   }

   /* gets sorting key for variables in C_2 corresponding to the following lifting sequence
    *  sequence 4: non-increasing a_j, i.e. a_1 >= a_2 >= ... >= a_|C_2|
    */
   for( j = 0; j < nvarsC2; j++ )
      sortkeysC2[j] = (SCIP_Real) weights[varsC2[j]];

   /* gets sorting key for variables in R corresponding to the following lifting sequence
    *  sequence 4: non-increasing a_j, i.e. a_1 >= a_2 >= ... >= a_|R|
    */
   for( j = 0; j < nvarsR; j++ )
      sortkeysR[j] = (SCIP_Real) weights[varsR[j]];

   /* sorts F, C2 and R */
   if( nvarsF > 0 )
   {
      SCIPsortDownPtrInt((void**)sortkeypairsF, varsF, compSortkeypairs, nvarsF);
   }
   if( nvarsC2 > 0 )
   {
      SCIPsortDownRealInt(sortkeysC2, varsC2, nvarsC2);
   }
   if( nvarsR > 0)
   {
      SCIPsortDownRealInt(sortkeysR, varsR, nvarsR);
   }

   /* frees temporary memory */
   SCIPfreeBufferArray(scip, &sortkeysR);
   SCIPfreeBufferArray(scip, &sortkeysC2);
   SCIPfreeBufferArray(scip, &sortkeypairsFstore);
   SCIPfreeBufferArray(scip, &sortkeypairsF);

   return SCIP_OKAY;
}

/** categorizes GUBs of knapsack GUB partion into GOC1, GNC1, GF, GC2, and GR and computes a lifting sequence of the GUBs
 *  for the sequential GUB wise lifting procedure
 */
static
SCIP_RETCODE getLiftingSequenceGUB(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_GUBSET*          gubset,             /**< GUB set data structure */
   SCIP_Real*            solvals,            /**< solution values of variables in knapsack constraint */
   SCIP_Longint*         weights,            /**< weights of variables in knapsack constraint */
   int*                  varsC1,             /**< variables in C1 */
   int*                  varsC2,             /**< variables in C2 */
   int*                  varsF,              /**< variables in F */
   int*                  varsR,              /**< variables in R */
   int                   nvarsC1,            /**< number of variables in C1 */
   int                   nvarsC2,            /**< number of variables in C2 */
   int                   nvarsF,             /**< number of variables in F */
   int                   nvarsR,             /**< number of variables in R */
   int*                  gubconsGC1,         /**< pointer to store GUBs in GC1(GNC1+GOC1) */
   int*                  gubconsGC2,         /**< pointer to store GUBs in GC2 */
   int*                  gubconsGFC1,        /**< pointer to store GUBs in GFC1(GNC1+GF) */
   int*                  gubconsGR,          /**< pointer to store GUBs in GR */
   int*                  ngubconsGC1,        /**< pointer to store number of GUBs in GC1(GNC1+GOC1) */
   int*                  ngubconsGC2,        /**< pointer to store number of GUBs in GC2 */
   int*                  ngubconsGFC1,       /**< pointer to store number of GUBs in GFC1(GNC1+GF) */
   int*                  ngubconsGR,         /**< pointer to store number of GUBs in GR */
   int*                  ngubconscapexceed,  /**< pointer to store number of GUBs with only capacity exceeding variables */
   int*                  maxgubvarssize      /**< pointer to store the maximal size of GUB constraints */
   )
{
   SORTKEYPAIR** sortkeypairsGFC1;
   SORTKEYPAIR* sortkeypairsGFC1store;
   SCIP_Real* sortkeysC1;
   SCIP_Real* sortkeysC2;
   SCIP_Real* sortkeysR;
   int* nC1varsingubcons;
   int var;
   int gubconsidx;
   int varidx;
   int ngubconss;
   int ngubconsGOC1;
   int targetvar;
   int nvarsprocessed;
   int i;
   int j;

#if GUBSPLITGNC1GUBS
   SCIP_Bool gubconswithF;
   int origngubconss;
   origngubconss = gubset->ngubconss;
#endif

   assert(scip != NULL);
   assert(gubset != NULL);
   assert(solvals != NULL);
   assert(weights != NULL);
   assert(varsC1 != NULL);
   assert(varsC2 != NULL);
   assert(varsF != NULL);
   assert(varsR != NULL);
   assert(nvarsC1 > 0);
   assert(nvarsC2 >= 0);
   assert(nvarsF >= 0);
   assert(nvarsR >= 0);
   assert(gubconsGC1 != NULL);
   assert(gubconsGC2 != NULL);
   assert(gubconsGFC1 != NULL);
   assert(gubconsGR != NULL);
   assert(ngubconsGC1 != NULL);
   assert(ngubconsGC2 != NULL);
   assert(ngubconsGFC1 != NULL);
   assert(ngubconsGR != NULL);
   assert(maxgubvarssize != NULL);

   ngubconss = gubset->ngubconss;
   nvarsprocessed = 0;
   ngubconsGOC1 = 0;

   /* GUBs are categorized into different types according to the variables in volved
    * - GOC1:  involves variables in C1 only           -- no C2, R, F
    * - GNC1:  involves variables in C1 and F (and R)  -- no C2
    * - GF:    involves variables in F  (and R) only   -- no C1, C2
    * - GC2:   involves variables in C2 only           -- no C1, R, F
    * - GR:    involves variables in R  only           -- no C1, C2, F
    * which requires splitting GUBs in case they include variable in F and R.
    *
    * afterwards all GUBs (except GOC1 GUBs, which we do not need to lift) are sorted by a two level lifting sequence.
    * - first  ordering level is: GFC1 (GNC1+GF), GC2, and GR.
    * - second ordering level is
    *    GFC1:   non-increasing number of variables in F and non-increasing max{x*_k : k in GFC1_j} in case of equality
    *    GC2:    non-increasing max{ a_k : k in GC2_j}; note that |GFC2_j| = 1
    *    GR:     non-increasing max{ a_k : k in GR_j}
    *
    * in additon, another GUB union, which is helpful for the lifting procedure, is formed
    * - GC1:   GUBs of category GOC1 and GNC1
    * with second ordering level non-decreasing min{ a_k : k in GC1_j };
    * note that min{ a_k : k in GC1_j } always comes from the first variable in the GUB
    */

   /* allocates temporary memory */
   SCIP_CALL( SCIPallocBufferArray(scip, &sortkeysC1, nvarsC1) );
   SCIP_CALL( SCIPallocBufferArray(scip, &sortkeysC2, nvarsC2) );
   SCIP_CALL( SCIPallocBufferArray(scip, &sortkeysR, nvarsR) );

   /* to get the GUB lifting sequence, we first sort all variables in F, C2, and R
    * - F:      non-increasing x*_j and non-increasing a_j in case of equality
    * - C2:     non-increasing a_j
    * - R:      non-increasing a_j
    * furthermore, sort C1 variables as needed for initializing the minweight table (non-increasing a_j).
    */

   /* gets sorting key for variables in C1 corresponding to the following ordering
    *  non-decreasing a_j, i.e. a_1 <= a_2 <= ... <= a_|C_1|
    */
   for( j = 0; j < nvarsC1; j++ )
   {
      /* gets sortkeys */
      sortkeysC1[j] = (SCIP_Real) weights[varsC1[j]];

      /* update status of variable in its gub constraint */
      gubconsidx = gubset->gubconssidx[varsC1[j]];
      varidx = gubset->gubvarsidx[varsC1[j]];
      gubset->gubconss[gubconsidx]->gubvarsstatus[varidx] = GUBVARSTATUS_BELONGSTOSET_C1;
   }

   /* gets sorting key for variables in F corresponding to the following ordering
    *  non-increasing x*_j, i.e., x*_1 >= x*_2 >= ... >= x*_|F|, and
    *  non-increasing a_j,  i.e., a_1  >= a_2  >= ... >= a_|F| in case of equality
    * and updates status of each variable in F in GUB set data structure
    */
   for( j = 0; j < nvarsF; j++ )
   {
      /* update status of variable in its gub constraint */
      gubconsidx = gubset->gubconssidx[varsF[j]];
      varidx = gubset->gubvarsidx[varsF[j]];
      gubset->gubconss[gubconsidx]->gubvarsstatus[varidx] = GUBVARSTATUS_BELONGSTOSET_F;
   }

   /* gets sorting key for variables in C2 corresponding to the following ordering
    *  non-increasing a_j,  i.e., a_1  >= a_2  >= ... >= a_|C2|
    * and updates status of each variable in F in GUB set data structure
    */
   for( j = 0; j < nvarsC2; j++ )
   {
      /* gets sortkeys */
      sortkeysC2[j] = (SCIP_Real) weights[varsC2[j]];

      /* update status of variable in its gub constraint */
      gubconsidx = gubset->gubconssidx[varsC2[j]];
      varidx = gubset->gubvarsidx[varsC2[j]];
      gubset->gubconss[gubconsidx]->gubvarsstatus[varidx] = GUBVARSTATUS_BELONGSTOSET_C2;
   }

   /* gets sorting key for variables in R corresponding to the following ordering
    *  non-increasing a_j,  i.e., a_1  >= a_2  >= ... >= a_|R|
    * and updates status of each variable in F in GUB set data structure
    */
   for( j = 0; j < nvarsR; j++ )
   {
      /* gets sortkeys */
      sortkeysR[j] = (SCIP_Real) weights[varsR[j]];

      /* update status of variable in its gub constraint */
      gubconsidx = gubset->gubconssidx[varsR[j]];
      varidx = gubset->gubvarsidx[varsR[j]];
      gubset->gubconss[gubconsidx]->gubvarsstatus[varidx] = GUBVARSTATUS_BELONGSTOSET_R;
   }

   /* sorts C1, F, C2 and R */
   assert(nvarsC1 > 0);
   SCIPsortRealInt(sortkeysC1, varsC1, nvarsC1);

   if( nvarsC2 > 0 )
   {
      SCIPsortDownRealInt(sortkeysC2, varsC2, nvarsC2);
   }
   if( nvarsR > 0)
   {
      SCIPsortDownRealInt(sortkeysR, varsR, nvarsR);
   }

   /* frees temporary memory */
   SCIPfreeBufferArray(scip, &sortkeysR);
   SCIPfreeBufferArray(scip, &sortkeysC2);
   SCIPfreeBufferArray(scip, &sortkeysC1);

   /* allocate and initialize temporary memory for sorting GUB constraints */
   SCIP_CALL( SCIPallocBufferArray(scip, &sortkeypairsGFC1, ngubconss) );
   SCIP_CALL( SCIPallocBufferArray(scip, &sortkeypairsGFC1store, ngubconss) );
   SCIP_CALL( SCIPallocBufferArray(scip, &nC1varsingubcons, ngubconss) );
   BMSclearMemoryArray(nC1varsingubcons, ngubconss);
   for( i = 0; i < ngubconss; i++)
   {
      sortkeypairsGFC1[i] = &(sortkeypairsGFC1store[i]);
      sortkeypairsGFC1[i]->key1 = 0.0;
      sortkeypairsGFC1[i]->key2 = 0.0;
   }
   *ngubconsGC1 = 0;
   *ngubconsGC2 = 0;
   *ngubconsGFC1 = 0;
   *ngubconsGR = 0;
   *ngubconscapexceed = 0;
   *maxgubvarssize = 0;

#ifndef NDEBUG
   for( i = 0; i < gubset->ngubconss; i++ )
      assert(gubset->gubconsstatus[i] == GUBCONSSTATUS_UNINITIAL);
#endif

   /* stores GUBs of group GC1 (GOC1+GNC1) and part of the GUBs of group GFC1 (GNC1 GUBs) and sorts variables in these GUBs
    * s.t. C1 variables come first (will automatically be sorted by non-decreasing weight).
    * gets sorting keys for GUBs of type GFC1 corresponding to the following ordering
    *    non-increasing number of variables in F, and
    *    non-increasing max{x*_k : k in GFC1_j} in case of equality
    */
   for( i = 0; i < nvarsC1; i++ )
   {
      int nvarsC1capexceed;

      nvarsC1capexceed = 0;

      var = varsC1[i];
      gubconsidx = gubset->gubconssidx[var];
      varidx = gubset->gubvarsidx[var];

      assert(gubconsidx >= 0 && gubconsidx < ngubconss);
      assert(gubset->gubconss[gubconsidx]->gubvarsstatus[varidx] == GUBVARSTATUS_BELONGSTOSET_C1);

      /* current C1 variable is put to the front of its GUB where C1 part is stored by non-decreasing weigth;
       * note that variables in C1 are already sorted by non-decreasing weigth
       */
      targetvar = gubset->gubconss[gubconsidx]->gubvars[nC1varsingubcons[gubconsidx]];
      GUBsetSwapVars(scip, gubset, var, targetvar);
      nC1varsingubcons[gubconsidx]++;

      /* the GUB was already handled (status set and stored in its group) by another variable of the GUB */
      if( gubset->gubconsstatus[gubconsidx] != GUBCONSSTATUS_UNINITIAL )
      {
         assert(gubset->gubconsstatus[gubconsidx] == GUBCONSSTATUS_BELONGSTOSET_GOC1
            || gubset->gubconsstatus[gubconsidx] == GUBCONSSTATUS_BELONGSTOSET_GNC1);
         continue;
      }

      /* determine the status of the current GUB constraint, GOC1 or GNC1; GUBs involving R variables are split into
       * GOC1/GNC1 and GF, if wanted. also update sorting key if GUB is of type GFC1 (GNC1)
       */
#if GUBSPLITGNC1GUBS
      gubconswithF = FALSE;
#endif
      for( j = 0; j < gubset->gubconss[gubconsidx]->ngubvars; j++ )
      {
         assert(gubset->gubconss[gubconsidx]->gubvarsstatus[j] != GUBVARSTATUS_BELONGSTOSET_C2);

         /* C1-variable: update number of C1/capacity exceeding variables */
         if( gubset->gubconss[gubconsidx]->gubvarsstatus[j] == GUBVARSTATUS_BELONGSTOSET_C1 )
         {
            nvarsC1capexceed++;
            nvarsprocessed++;
         }
         /* F-variable: update sort key (number of F variables in GUB) of corresponding GFC1-GUB */
         else if( gubset->gubconss[gubconsidx]->gubvarsstatus[j] == GUBVARSTATUS_BELONGSTOSET_F )
         {
#if GUBSPLITGNC1GUBS
            gubconswithF = TRUE;
#endif
            sortkeypairsGFC1[*ngubconsGFC1]->key1 += 1.0;

            if( solvals[gubset->gubconss[gubconsidx]->gubvars[j]] > sortkeypairsGFC1[*ngubconsGFC1]->key2 )
               sortkeypairsGFC1[*ngubconsGFC1]->key2 = solvals[gubset->gubconss[gubconsidx]->gubvars[j]];
         }
         else if( gubset->gubconss[gubconsidx]->gubvarsstatus[j] == GUBVARSTATUS_CAPACITYEXCEEDED )
         {
            nvarsC1capexceed++;
         }
         else
            assert(gubset->gubconss[gubconsidx]->gubvarsstatus[j] == GUBVARSTATUS_BELONGSTOSET_R);
      }

      /* update set of GC1 GUBs */
      gubconsGC1[*ngubconsGC1] = gubconsidx;
      (*ngubconsGC1)++;

      /* update maximum size of all GUB constraints */
      if( gubset->gubconss[gubconsidx]->gubvarssize > *maxgubvarssize )
         *maxgubvarssize = gubset->gubconss[gubconsidx]->gubvarssize;

      /* set status of GC1-GUB (GOC1 or GNC1) and update set of GFC1 GUBs */
      if( nvarsC1capexceed == gubset->gubconss[gubconsidx]->ngubvars )
      {
         gubset->gubconsstatus[gubconsidx] = GUBCONSSTATUS_BELONGSTOSET_GOC1;
         ngubconsGOC1++;
      }
      else
      {
#if GUBSPLITGNC1GUBS
         /* only variables in C1 and R -- no in F: GUB will be split into GR and GOC1 GUBs */
         if( !gubconswithF )
         {
            GUBVARSTATUS movevarstatus;

            assert(gubset->ngubconss < gubset->nvars);

            /* create a new GUB for GR part of splitting */
            SCIP_CALL( GUBconsCreate(scip, &gubset->gubconss[gubset->ngubconss]) );
            gubset->ngubconss++;
            ngubconss = gubset->ngubconss;

            /* fill GR with R variables in current GUB */
            for( j = gubset->gubconss[gubconsidx]->ngubvars-1; j >= 0; j-- )
            {
                movevarstatus = gubset->gubconss[gubconsidx]->gubvarsstatus[j];
                if( movevarstatus != GUBVARSTATUS_BELONGSTOSET_C1 )
                {
                   assert(movevarstatus == GUBVARSTATUS_BELONGSTOSET_R || movevarstatus == GUBVARSTATUS_CAPACITYEXCEEDED);
                   SCIP_CALL( GUBsetMoveVar(scip, gubset, vars, gubset->gubconss[gubconsidx]->gubvars[j],
                         gubconsidx, ngubconss-1) );
                   gubset->gubconss[ngubconss-1]->gubvarsstatus[gubset->gubconss[ngubconss-1]->ngubvars-1] =
                      movevarstatus;
                }
            }

            gubset->gubconsstatus[gubconsidx] = GUBCONSSTATUS_BELONGSTOSET_GOC1;
            ngubconsGOC1++;

            gubset->gubconsstatus[ngubconss-1] = GUBCONSSTATUS_BELONGSTOSET_GR;
            gubconsGR[*ngubconsGR] = ngubconss-1;
            (*ngubconsGR)++;
         }
         /* variables in C1, F, and maybe R: GNC1 GUB */
         else
         {
            assert(gubconswithF);

            gubset->gubconsstatus[gubconsidx] = GUBCONSSTATUS_BELONGSTOSET_GNC1;
            gubconsGFC1[*ngubconsGFC1] = gubconsidx;
            (*ngubconsGFC1)++;
         }
#else
         gubset->gubconsstatus[gubconsidx] = GUBCONSSTATUS_BELONGSTOSET_GNC1;
         gubconsGFC1[*ngubconsGFC1] = gubconsidx;
         (*ngubconsGFC1)++;
#endif
      }
   }

   /* stores GUBs of group GC2 (only trivial GUBs); sorting is not required because the C2 variables (which we loop over)
    * are already sorted correctly
    */
   for( i = 0; i < nvarsC2; i++ )
   {
      var = varsC2[i];
      gubconsidx = gubset->gubconssidx[var];
      varidx = gubset->gubvarsidx[var];

      assert(gubconsidx >= 0 && gubconsidx < ngubconss);
      assert(gubset->gubconss[gubconsidx]->ngubvars == 1);
      assert(varidx == 0);
      assert(gubset->gubconss[gubconsidx]->gubvarsstatus[varidx] == GUBVARSTATUS_BELONGSTOSET_C2);
      assert(gubset->gubconsstatus[gubconsidx] == GUBCONSSTATUS_UNINITIAL);

      /* set status of GC2 GUB */
      gubset->gubconsstatus[gubconsidx] = GUBCONSSTATUS_BELONGSTOSET_GC2;

      /* update group of GC2 GUBs */
      gubconsGC2[*ngubconsGC2] = gubconsidx;
      (*ngubconsGC2)++;

      /* update maximum size of all GUB constraints */
      if( gubset->gubconss[gubconsidx]->gubvarssize > *maxgubvarssize )
         *maxgubvarssize = gubset->gubconss[gubconsidx]->gubvarssize;

      nvarsprocessed++;
   }

   /* stores remaining part of the GUBs of group GFC1 (GF GUBs) and gets GUB sorting keys corresp. to following ordering
    *    non-increasing number of variables in F, and
    *    non-increasing max{x*_k : k in GFC1_j} in case of equality
    */
   for( i = 0; i < nvarsF; i++ )
   {
      var = varsF[i];
      gubconsidx = gubset->gubconssidx[var];
      varidx = gubset->gubvarsidx[var];

      assert(gubconsidx >= 0 && gubconsidx < ngubconss);
      assert(gubset->gubconss[gubconsidx]->gubvarsstatus[varidx] == GUBVARSTATUS_BELONGSTOSET_F);

      nvarsprocessed++;

      /* the GUB was already handled (status set and stored in its group) by another variable of the GUB */
      if( gubset->gubconsstatus[gubconsidx] != GUBCONSSTATUS_UNINITIAL )
      {
         assert(gubset->gubconsstatus[gubconsidx] == GUBCONSSTATUS_BELONGSTOSET_GF
              || gubset->gubconsstatus[gubconsidx] == GUBCONSSTATUS_BELONGSTOSET_GNC1);
         continue;
      }

      /* set status of GF GUB */
      gubset->gubconsstatus[gubconsidx] = GUBCONSSTATUS_BELONGSTOSET_GF;

      /* update sorting key of corresponding GFC1 GUB */
      for( j = 0; j < gubset->gubconss[gubconsidx]->ngubvars; j++ )
      {
         assert(gubset->gubconss[gubconsidx]->gubvarsstatus[j] != GUBVARSTATUS_BELONGSTOSET_C2
            && gubset->gubconss[gubconsidx]->gubvarsstatus[j] != GUBVARSTATUS_BELONGSTOSET_C1);

         /* F-variable: update sort key (number of F variables in GUB) of corresponding GFC1-GUB */
         if( gubset->gubconss[gubconsidx]->gubvarsstatus[j] == GUBVARSTATUS_BELONGSTOSET_F )
         {
            sortkeypairsGFC1[*ngubconsGFC1]->key1 += 1.0;

            if( solvals[gubset->gubconss[gubconsidx]->gubvars[j]] > sortkeypairsGFC1[*ngubconsGFC1]->key2 )
               sortkeypairsGFC1[*ngubconsGFC1]->key2 = solvals[gubset->gubconss[gubconsidx]->gubvars[j]];
         }
      }

      /* update set of GFC1 GUBs */
      gubconsGFC1[*ngubconsGFC1] = gubconsidx;
      (*ngubconsGFC1)++;

      /* update maximum size of all GUB constraints */
      if( gubset->gubconss[gubconsidx]->gubvarssize > *maxgubvarssize )
         *maxgubvarssize = gubset->gubconss[gubconsidx]->gubvarssize;
   }

   /* stores GUBs of group GR; sorting is not required because the R variables (which we loop over) are already sorted
    * correctly
    */
   for( i = 0; i < nvarsR; i++ )
   {
      var = varsR[i];
      gubconsidx = gubset->gubconssidx[var];
      varidx = gubset->gubvarsidx[var];

      assert(gubconsidx >= 0 && gubconsidx < ngubconss);
      assert(gubset->gubconss[gubconsidx]->gubvarsstatus[varidx] == GUBVARSTATUS_BELONGSTOSET_R);

      nvarsprocessed++;

      /* the GUB was already handled (status set and stored in its group) by another variable of the GUB */
      if( gubset->gubconsstatus[gubconsidx] != GUBCONSSTATUS_UNINITIAL )
      {
         assert(gubset->gubconsstatus[gubconsidx] == GUBCONSSTATUS_BELONGSTOSET_GR
              || gubset->gubconsstatus[gubconsidx] == GUBCONSSTATUS_BELONGSTOSET_GF
              || gubset->gubconsstatus[gubconsidx] == GUBCONSSTATUS_BELONGSTOSET_GNC1);
         continue;
      }

      /* set status of GR GUB */
      gubset->gubconsstatus[gubconsidx] = GUBCONSSTATUS_BELONGSTOSET_GR;

      /* update set of GR GUBs */
      gubconsGR[*ngubconsGR] = gubconsidx;
      (*ngubconsGR)++;

      /* update maximum size of all GUB constraints */
      if( gubset->gubconss[gubconsidx]->gubvarssize > *maxgubvarssize )
         *maxgubvarssize = gubset->gubconss[gubconsidx]->gubvarssize;
   }
   assert(nvarsprocessed == nvarsC1 + nvarsC2 + nvarsF + nvarsR);

   /* update number of GUBs with only capacity exceeding variables (will not be used for lifting) */
   (*ngubconscapexceed) =  ngubconss - (ngubconsGOC1 + (*ngubconsGC2) + (*ngubconsGFC1) + (*ngubconsGR));
   assert(*ngubconscapexceed >= 0);
#ifndef NDEBUG
   {
      int check;

      check = 0;

      /* remaining not handled GUBs should only contain capacity exceeding variables */
      for( i = 0; i < ngubconss; i++ )
      {
         if( gubset->gubconsstatus[i] ==  GUBCONSSTATUS_UNINITIAL )
            check++;
      }
      assert(check == *ngubconscapexceed);
   }
#endif

   /* sort GFCI GUBs according to computed sorting keys */
   if( (*ngubconsGFC1) > 0 )
   {
      SCIPsortDownPtrInt((void**)sortkeypairsGFC1, gubconsGFC1, compSortkeypairs, (*ngubconsGFC1));
   }

   /* free temporary memory */
#if GUBSPLITGNC1GUBS
   ngubconss = origngubconss;
#endif
   SCIPfreeBufferArray(scip, &nC1varsingubcons);
   SCIPfreeBufferArray(scip, &sortkeypairsGFC1store);
   SCIPfreeBufferArray(scip, &sortkeypairsGFC1);

   return SCIP_OKAY;
}

/** enlarges minweight table to at least the given length */
static
SCIP_RETCODE enlargeMinweights(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_Longint**        minweightsptr,      /**< pointer to minweights table */
   int*                  minweightslen,      /**< pointer to store number of entries in minweights table (incl. z=0) */
   int*                  minweightssize,     /**< pointer to current size of minweights table */
   int                   newlen              /**< new length of minweights table */
   )
{
   int j;

   assert(minweightsptr != NULL);
   assert(*minweightsptr != NULL);
   assert(minweightslen != NULL);
   assert(*minweightslen >= 0);
   assert(minweightssize != NULL);
   assert(*minweightssize >= 0);

   if( newlen > *minweightssize )
   {
      int newsize;

      /* reallocate table memory */
      newsize = SCIPcalcMemGrowSize(scip, newlen);
      SCIP_CALL( SCIPreallocBufferArray(scip, minweightsptr, newsize) );
      *minweightssize = newsize;
   }
   assert(newlen <= *minweightssize);

   /* initialize new elements */
   for( j = *minweightslen; j < newlen; ++j )
      (*minweightsptr)[j] = SCIP_LONGINT_MAX;
   *minweightslen = newlen;

   return SCIP_OKAY;
}

/** lifts given inequality
 *    sum_{j in M_1} x_j <= alpha_0
 *  valid for
 *    S^0 = { x in {0,1}^|M_1| : sum_{j in M_1} a_j x_j <= a_0 - sum_{j in M_2} a_j }
 *  to a valid inequality
 *    sum_{j in M_1} x_j + sum_{j in F} alpha_j x_j + sum_{j in M_2} alpha_j x_j + sum_{j in R} alpha_j x_j
 *    <= alpha_0 + sum_{j in M_2} alpha_j
 *  for
 *    S = { x in {0,1}^|N| : sum_{j in N} a_j x_j <= a_0 };
 *  uses sequential up-lifting for the variables in F, sequential down-lifting for the variable in M_2, and
 *  sequential up-lifting for the variables in R; procedure can be used to strengthen minimal cover inequalities and
 *  extended weight inequalities.
 */
static
SCIP_RETCODE sequentialUpAndDownLifting(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_VAR**            vars,               /**< variables in knapsack constraint */
   int                   nvars,              /**< number of variables in knapsack constraint */
   int                   ntightened,         /**< number of variables with tightened upper bound */
   SCIP_Longint*         weights,            /**< weights of variables in knapsack constraint */
   SCIP_Longint          capacity,           /**< capacity of knapsack */
   SCIP_Real*            solvals,            /**< solution values of all problem variables */
   int*                  varsM1,             /**< variables in M_1 */
   int*                  varsM2,             /**< variables in M_2 */
   int*                  varsF,              /**< variables in F */
   int*                  varsR,              /**< variables in R */
   int                   nvarsM1,            /**< number of variables in M_1 */
   int                   nvarsM2,            /**< number of variables in M_2 */
   int                   nvarsF,             /**< number of variables in F */
   int                   nvarsR,             /**< number of variables in R */
   int                   alpha0,             /**< rights hand side of given valid inequality */
   int*                  liftcoefs,          /**< pointer to store lifting coefficient of vars in knapsack constraint */
   SCIP_Real*            cutact,             /**< pointer to store activity of lifted valid inequality */
   int*                  liftrhs             /**< pointer to store right hand side of the lifted valid inequality */
   )
{
   SCIP_Longint* minweights;
   SCIP_Real* sortkeys;
   SCIP_Longint fixedonesweight;
   int minweightssize;
   int minweightslen;
   int j;
   int w;

   assert(scip != NULL);
   assert(vars != NULL);
   assert(nvars >= 0);
   assert(weights != NULL);
   assert(capacity >= 0);
   assert(solvals != NULL);
   assert(varsM1 != NULL);
   assert(varsM2 != NULL);
   assert(varsF != NULL);
   assert(varsR != NULL);
   assert(nvarsM1 >= 0 && nvarsM1 <= nvars - ntightened);
   assert(nvarsM2 >= 0 && nvarsM2 <= nvars - ntightened);
   assert(nvarsF >= 0 && nvarsF <= nvars - ntightened);
   assert(nvarsR >= 0 && nvarsR <= nvars - ntightened);
   assert(nvarsM1 + nvarsM2 + nvarsF + nvarsR == nvars  - ntightened);
   assert(alpha0 >= 0);
   assert(liftcoefs != NULL);
   assert(cutact != NULL);
   assert(liftrhs != NULL);

   /* allocates temporary memory */
   minweightssize = nvarsM1 + 1;
   SCIP_CALL( SCIPallocBufferArray(scip, &minweights, minweightssize) );
   SCIP_CALL( SCIPallocBufferArray(scip, &sortkeys, nvarsM1) );

   /* initializes data structures */
   BMSclearMemoryArray(liftcoefs, nvars);
   *cutact = 0.0;

   /* sets lifting coefficient of variables in M1, sorts variables in M1 such that a_1 <= a_2 <= ... <= a_|M1|
    * and calculates activity of the current valid inequality
    */
   for( j = 0; j < nvarsM1; j++ )
   {
      assert(liftcoefs[varsM1[j]] == 0);
      liftcoefs[varsM1[j]] = 1;
      sortkeys[j] = (SCIP_Real) (weights[varsM1[j]]);
      (*cutact) += solvals[varsM1[j]];
   }

   SCIPsortRealInt(sortkeys, varsM1, nvarsM1);

   /* initializes (i = 1) the minweight table, defined as: minweights_i[w] =
    *   min   sum_{j in M_1} a_j x_j + sum_{k=1}^{i-1} a_{j_k}     x_{j_k}
    *   s.t.  sum_{j in M_1}     x_j + sum_{k=1}^{i-1} alpha_{j_k} x_{j_k} >= w
    *                                    x_j in {0,1} for j in M_1 & {j_i,...,j_i-1},
    * for i = 1,...,t with t = |N\M1| and w = 0,...,|M1| + sum_{k=1}^{i-1} alpha_{j_k};
    */
   minweights[0] = 0;
   for( w = 1; w <= nvarsM1; w++ )
      minweights[w] = minweights[w-1] + weights[varsM1[w-1]];
   minweightslen = nvarsM1 + 1;

   /* gets sum of weights of variables fixed to one, i.e. sum of weights of variables in M_2 */
   fixedonesweight = 0;
   for( j = 0; j < nvarsM2; j++ )
      fixedonesweight += weights[varsM2[j]];
   assert(fixedonesweight >= 0);

   /* initializes right hand side of lifted valid inequality */
   *liftrhs = alpha0;

   /* sequentially up-lifts all variables in F: */
   for( j = 0; j < nvarsF; j++ )
   {
      SCIP_Longint weight;
      int liftvar;
      int liftcoef;
      int z;

      liftvar = varsF[j];
      weight = weights[liftvar];
      assert(liftvar >= 0 && liftvar < nvars);
      assert(SCIPisFeasGT(scip, solvals[liftvar], 0.0));
      assert(weight > 0);

      /* knapsack problem is infeasible:
       *   sets z = 0
       */
      if( capacity - fixedonesweight - weight < 0 )
      {
         z = 0;
      }
      /* knapsack problem is feasible:
       *   sets z = max { w : 0 <= w <= liftrhs, minweights_i[w] <= a_0 - fixedonesweight - a_{j_i}  } = liftrhs,
       *   if minweights_i[liftrhs] <= a_0 - fixedonesweight - a_{j_i}
       */
      else if( minweights[*liftrhs] <= capacity - fixedonesweight - weight )
      {
         z = *liftrhs;
      }
      /* knapsack problem is feasible:
       *   uses binary search to find z = max { w : 0 <= w <= liftrhs, minweights_i[w] <= a_0 - fixedonesweight - a_{j_i} }
       */
      else
      {
         int left;
         int right;
         int middle;

         assert((*liftrhs) + 1 >= minweightslen || minweights[(*liftrhs) + 1] > capacity - fixedonesweight - weight);
         left = 0;
         right = (*liftrhs) + 1;
         while( left < right - 1 )
         {
            middle = (left + right) / 2;
            assert(0 <= middle && middle < minweightslen);
            if( minweights[middle] <= capacity - fixedonesweight - weight )
               left = middle;
            else
               right = middle;
         }
         assert(left == right - 1);
         assert(0 <= left && left < minweightslen);
         assert(minweights[left] <= capacity - fixedonesweight - weight );
         assert(left == minweightslen - 1 || minweights[left+1] > capacity - fixedonesweight - weight);

         /* now z = left */
         z = left;
         assert(z <= *liftrhs);
      }

      /* calculates lifting coefficients alpha_{j_i} = liftrhs - z */
      liftcoef = (*liftrhs) - z;
      liftcoefs[liftvar] = liftcoef;
      assert(liftcoef >= 0 && liftcoef <= (*liftrhs) + 1);

      /* minweight table and activity of current valid inequality will not change, if alpha_{j_i} = 0 */
      if( liftcoef == 0 )
         continue;

      /* updates activity of current valid inequality */
      (*cutact) += liftcoef * solvals[liftvar];

      /* enlarges current minweight table:
       *  from minweightlen = |M1| + sum_{k=1}^{i-1} alpha_{j_k} + 1 entries
       *  to                  |M1| + sum_{k=1}^{i  } alpha_{j_k} + 1 entries
       * and sets minweights_i[w] = infinity for
       *  w = |M1| + sum_{k=1}^{i-1} alpha_{j_k} + 1 , ... , |M1| + sum_{k=1}^{i} alpha_{j_k}
       */
      SCIP_CALL( enlargeMinweights(scip, &minweights, &minweightslen, &minweightssize, minweightslen + liftcoef) );

      /* updates minweight table: minweight_i+1[w] =
       *   min{ minweights_i[w], a_{j_i}},                                 if w <  alpha_j_i
       *   min{ minweights_i[w], minweights_i[w - alpha_j_i] + a_j_i},     if w >= alpha_j_i
       */
      for( w = minweightslen - 1; w >= 0; w-- )
      {
         SCIP_Longint min;
         if( w < liftcoef )
         {
            min = MIN(minweights[w], weight);
            minweights[w] = min;
         }
         else
         {
            assert(w >= liftcoef);
            min = MIN(minweights[w], minweights[w - liftcoef] + weight);
            minweights[w] = min;
         }
      }
   }
   assert(minweights[0] == 0);

   /* sequentially down-lifts all variables in M_2: */
   for( j = 0; j < nvarsM2; j++ )
   {
      SCIP_Longint weight;
      int liftvar;
      int liftcoef;
      int left;
      int right;
      int middle;
      int z;

      liftvar = varsM2[j];
      weight = weights[liftvar];
      assert(SCIPisFeasEQ(scip, solvals[liftvar], 1.0));
      assert(liftvar >= 0 && liftvar < nvars);
      assert(weight > 0);

      /* uses binary search to find
       *   z = max { w : 0 <= w <= |M_1| + sum_{k=1}^{i-1} alpha_{j_k}, minweights_[w] <= a_0 - fixedonesweight + a_{j_i}}
       */
      left = 0;
      right = minweightslen;
      while( left < right - 1 )
      {
         middle = (left + right) / 2;
         assert(0 <= middle && middle < minweightslen);
         if( minweights[middle] <= capacity - fixedonesweight + weight )
            left = middle;
         else
            right = middle;
      }
      assert(left == right - 1);
      assert(0 <= left && left < minweightslen);
      assert(minweights[left] <= capacity - fixedonesweight + weight );
      assert(left == minweightslen - 1 || minweights[left+1] > capacity - fixedonesweight + weight);

      /* now z = left */
      z = left;
      assert(z >= *liftrhs);

      /* calculates lifting coefficients alpha_{j_i} = z - liftrhs */
      liftcoef = z - (*liftrhs);
      liftcoefs[liftvar] = liftcoef;
      assert(liftcoef >= 0);

      /* updates sum of weights of variables fixed to one */
      fixedonesweight -= weight;

      /* updates right-hand side of current valid inequality */
      (*liftrhs) += liftcoef;
      assert(*liftrhs >= alpha0);

      /* minweight table and activity of current valid inequality will not change, if alpha_{j_i} = 0 */
      if( liftcoef == 0 )
         continue;

      /* updates activity of current valid inequality */
      (*cutact) += liftcoef * solvals[liftvar];

      /* enlarges current minweight table:
       *  from minweightlen = |M1| + sum_{k=1}^{i-1} alpha_{j_k} + 1 entries
       *  to                  |M1| + sum_{k=1}^{i  } alpha_{j_k} + 1 entries
       * and sets minweights_i[w] = infinity for
       *  w = |M1| + sum_{k=1}^{i-1} alpha_{j_k} + 1 , ... , |M1| + sum_{k=1}^{i} alpha_{j_k}
       */
      SCIP_CALL( enlargeMinweights(scip, &minweights, &minweightslen, &minweightssize, minweightslen + liftcoef) );

      /* updates minweight table: minweight_i+1[w] =
       *   min{ minweights_i[w], a_{j_i}},                                 if w <  alpha_j_i
       *   min{ minweights_i[w], minweights_i[w - alpha_j_i] + a_j_i},     if w >= alpha_j_i
       */
      for( w = minweightslen - 1; w >= 0; w-- )
      {
         SCIP_Longint min;
         if( w < liftcoef )
         {
            min = MIN(minweights[w], weight);
            minweights[w] = min;
         }
         else
         {
            assert(w >= liftcoef);
            min = MIN(minweights[w], minweights[w - liftcoef] + weight);
            minweights[w] = min;
         }
      }
   }
   assert(fixedonesweight == 0);
   assert(*liftrhs >= alpha0);

   /* sequentially up-lifts all variables in R: */
   for( j = 0; j < nvarsR; j++ )
   {
      SCIP_Longint weight;
      int liftvar;
      int liftcoef;
      int z;

      liftvar = varsR[j];
      weight = weights[liftvar];
      assert(liftvar >= 0 && liftvar < nvars);
      assert(SCIPisFeasEQ(scip, solvals[liftvar], 0.0));
      assert(weight > 0);
      assert(capacity - weight >= 0);
      assert((*liftrhs) + 1 >= minweightslen || minweights[(*liftrhs) + 1] > capacity - weight);

      /* sets z = max { w : 0 <= w <= liftrhs, minweights_i[w] <= a_0 - a_{j_i} } = liftrhs,
       * if minweights_i[liftrhs] <= a_0 - a_{j_i}
       */
      if( minweights[*liftrhs] <= capacity - weight )
      {
         z = *liftrhs;
      }
      /* uses binary search to find z = max { w : 0 <= w <= liftrhs, minweights_i[w] <= a_0 - a_{j_i} }
       */
      else
      {
         int left;
         int right;
         int middle;

         left = 0;
         right = (*liftrhs) + 1;
         while( left < right - 1)
         {
            middle = (left + right) / 2;
            assert(0 <= middle && middle < minweightslen);
            if( minweights[middle] <= capacity - weight )
               left = middle;
            else
               right = middle;
         }
         assert(left == right - 1);
         assert(0 <= left && left < minweightslen);
         assert(minweights[left] <= capacity - weight );
         assert(left == minweightslen - 1 || minweights[left+1] > capacity - weight);

         /* now z = left */
         z = left;
         assert(z <= *liftrhs);
      }

      /* calculates lifting coefficients alpha_{j_i} = liftrhs - z */
      liftcoef = (*liftrhs) - z;
      liftcoefs[liftvar] = liftcoef;
      assert(liftcoef >= 0 && liftcoef <= *liftrhs);

      /* minweight table and activity of current valid inequality will not change, if alpha_{j_i} = 0 */
      if( liftcoef == 0 )
         continue;

      /* updates activity of current valid inequality */
      (*cutact) += liftcoef * solvals[liftvar];

      /* updates minweight table: minweight_i+1[w] =
       *   min{ minweight_i[w], a_{j_i}},                                if w <  alpha_j_i
       *   min{ minweight_i[w], minweight_i[w - alpha_j_i] + a_j_i},     if w >= alpha_j_i
       */
      for( w = *liftrhs; w >= 0; w-- )
      {
         SCIP_Longint min;
         if( w < liftcoef )
         {
            min = MIN(minweights[w], weight);
            minweights[w] = min;
         }
         else
         {
            assert(w >= liftcoef);
            min = MIN(minweights[w], minweights[w - liftcoef] + weight);
            minweights[w] = min;
         }
      }
   }

   /* frees temporary memory */
   SCIPfreeBufferArray(scip, &sortkeys);
   SCIPfreeBufferArray(scip, &minweights);

   return SCIP_OKAY;
}

/** adds two minweight values in a safe way, i.e,, ensures no overflow */
static
SCIP_Longint safeAddMinweightsGUB(
   SCIP_Longint          val1,               /**< first value to add */
   SCIP_Longint          val2                /**< second value to add */
   )
{
   assert(val1 >= 0);
   assert(val2 >= 0);

   if( val1 >= SCIP_LONGINT_MAX || val2 >= SCIP_LONGINT_MAX )
      return SCIP_LONGINT_MAX;
   else
   {
      assert(val1 <= SCIP_LONGINT_MAX - val2);
      return (val1 + val2);
   }
}

/** computes minweights table for lifting with GUBs by combining unfished and fished tables */
static
void computeMinweightsGUB(
   SCIP_Longint*         minweights,         /**< minweight table to compute */
   SCIP_Longint*         finished,           /**< given finished table */
   SCIP_Longint*         unfinished,         /**< given unfinished table */
   int                   minweightslen       /**< length of minweight, finished, and unfinished tables */
   )
{
   int w1;
   int w2;

   /* minweights_i[w] = min{finished_i[w1] + unfinished_i[w2] : w1>=0, w2>=0, w1+w2=w};
    * note that finished and unfished arrays sorted by non-decreasing weight
    */

   /* initialize minweight with w2 = 0 */
   w2 = 0;
   assert(unfinished[w2] == 0);
   for( w1 = 0; w1 < minweightslen; w1++ )
      minweights[w1] = finished[w1];

   /* consider w2 = 1, ..., minweightslen-1 */
   for( w2 = 1; w2 < minweightslen; w2++ )
   {
      if( unfinished[w2] >= SCIP_LONGINT_MAX )
         break;

      for( w1 = 0; w1 < minweightslen - w2; w1++ )
      {
         SCIP_Longint temp;

         temp = safeAddMinweightsGUB(finished[w1], unfinished[w2]);
         if( temp <= minweights[w1+w2] )
            minweights[w1+w2] = temp;
      }
   }
}

/** lifts given inequality
 *    sum_{j in C_1} x_j <= alpha_0
 *  valid for
 *    S^0 = { x in {0,1}^|C_1| : sum_{j in C_1} a_j x_j <= a_0 - sum_{j in C_2} a_j;
 *                               sum_{j in Q_i} x_j <= 1, forall i in I }
 *  to a valid inequality
 *    sum_{j in C_1} x_j + sum_{j in F} alpha_j x_j + sum_{j in C_2} alpha_j x_j + sum_{j in R} alpha_j x_j
 *    <= alpha_0 + sum_{j in C_2} alpha_j
 *  for
 *    S = { x in {0,1}^|N| : sum_{j in N} a_j x_j <= a_0; sum_{j in Q_i} x_j <= 1, forall i in I };
 *  uses sequential up-lifting   for the variables in GUB constraints in gubconsGFC1,
 *       sequential down-lifting for the variables in GUB constraints in gubconsGC2, and
 *       sequential up-lifting   for the variabels in GUB constraints in gubconsGR.
 */
static
SCIP_RETCODE sequentialUpAndDownLiftingGUB(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_GUBSET*          gubset,             /**< GUB set data structure */
   SCIP_VAR**            vars,               /**< variables in knapsack constraint */
   int                   ngubconscapexceed,  /**< number of GUBs with only capacity exceeding variables */
   SCIP_Longint*         weights,            /**< weights of variables in knapsack constraint */
   SCIP_Longint          capacity,           /**< capacity of knapsack */
   SCIP_Real*            solvals,            /**< solution values of all knapsack variables */
   int*                  gubconsGC1,         /**< GUBs in GC1(GNC1+GOC1) */
   int*                  gubconsGC2,         /**< GUBs in GC2 */
   int*                  gubconsGFC1,        /**< GUBs in GFC1(GNC1+GF) */
   int*                  gubconsGR,          /**< GUBs in GR */
   int                   ngubconsGC1,        /**< number of GUBs in GC1(GNC1+GOC1) */
   int                   ngubconsGC2,        /**< number of GUBs in GC2 */
   int                   ngubconsGFC1,       /**< number of GUBs in GFC1(GNC1+GF) */
   int                   ngubconsGR,         /**< number of GUBs in GR */
   int                   alpha0,             /**< rights hand side of given valid inequality */
   int*                  liftcoefs,          /**< pointer to store lifting coefficient of vars in knapsack constraint */
   SCIP_Real*            cutact,             /**< pointer to store activity of lifted valid inequality */
   int*                  liftrhs,            /**< pointer to store right hand side of the lifted valid inequality */
   int                   maxgubvarssize      /**< maximal size of GUB constraints */
   )
{
   SCIP_Longint* minweights;
   SCIP_Longint* finished;
   SCIP_Longint* unfinished;
   int* gubconsGOC1;
   int* gubconsGNC1;
   int* liftgubvars;
   SCIP_Longint fixedonesweight;
   SCIP_Longint weight;
   SCIP_Longint weightdiff1;
   SCIP_Longint weightdiff2;
   SCIP_Longint min;
   int minweightssize;
   int minweightslen;
   int nvars;
   int varidx;
   int liftgubconsidx;
   int liftvar;
   int sumliftcoef;
   int liftcoef;
   int ngubconsGOC1;
   int ngubconsGNC1;
   int left;
   int right;
   int middle;
   int nliftgubvars;
   int tmplen;
   int tmpsize;
   int j;
   int k;
   int w;
   int z;
#ifndef NDEBUG
   int ngubconss;
   int nliftgubC1;

   assert(gubset != NULL);
   ngubconss = gubset->ngubconss;
#else
   assert(gubset != NULL);
#endif

   nvars = gubset->nvars;

   assert(scip != NULL);
   assert(vars != NULL);
   assert(nvars >= 0);
   assert(weights != NULL);
   assert(capacity >= 0);
   assert(solvals != NULL);
   assert(gubconsGC1 != NULL);
   assert(gubconsGC2 != NULL);
   assert(gubconsGFC1 != NULL);
   assert(gubconsGR != NULL);
   assert(ngubconsGC1 >= 0 && ngubconsGC1 <= ngubconss - ngubconscapexceed);
   assert(ngubconsGC2 >= 0 && ngubconsGC2 <= ngubconss - ngubconscapexceed);
   assert(ngubconsGFC1 >= 0 && ngubconsGFC1 <= ngubconss - ngubconscapexceed);
   assert(ngubconsGR >= 0 && ngubconsGR <= ngubconss - ngubconscapexceed);
   assert(alpha0 >= 0);
   assert(liftcoefs != NULL);
   assert(cutact != NULL);
   assert(liftrhs != NULL);

   minweightssize = ngubconsGC1+1;

   /* allocates temporary memory */
   SCIP_CALL( SCIPallocBufferArray(scip, &liftgubvars, maxgubvarssize) );
   SCIP_CALL( SCIPallocBufferArray(scip, &gubconsGOC1, ngubconsGC1) );
   SCIP_CALL( SCIPallocBufferArray(scip, &gubconsGNC1, ngubconsGC1) );
   SCIP_CALL( SCIPallocBufferArray(scip, &minweights, minweightssize) );
   SCIP_CALL( SCIPallocBufferArray(scip, &finished, minweightssize) );
   SCIP_CALL( SCIPallocBufferArray(scip, &unfinished, minweightssize) );

   /* initializes data structures */
   BMSclearMemoryArray(liftcoefs, nvars);
   *cutact = 0.0;

   /* gets GOC1 and GNC1 GUBs, sets lifting coefficient of variables in C1 and calculates activity of the current
    * valid inequality
    */
   ngubconsGOC1 = 0;
   ngubconsGNC1 = 0;
   for( j = 0; j < ngubconsGC1; j++ )
   {
      if( gubset->gubconsstatus[gubconsGC1[j]] == GUBCONSSTATUS_BELONGSTOSET_GOC1 )
      {
         gubconsGOC1[ngubconsGOC1] = gubconsGC1[j];
         ngubconsGOC1++;
      }
      else
      {
         assert(gubset->gubconsstatus[gubconsGC1[j]] == GUBCONSSTATUS_BELONGSTOSET_GNC1);
         gubconsGNC1[ngubconsGNC1] = gubconsGC1[j];
         ngubconsGNC1++;
      }
      for( k = 0; k < gubset->gubconss[gubconsGC1[j]]->ngubvars
              && gubset->gubconss[gubconsGC1[j]]->gubvarsstatus[k] == GUBVARSTATUS_BELONGSTOSET_C1; k++ )
      {
         varidx = gubset->gubconss[gubconsGC1[j]]->gubvars[k];
         assert(varidx >= 0 && varidx < nvars);
         assert(liftcoefs[varidx] == 0);

         liftcoefs[varidx] = 1;
         (*cutact) += solvals[varidx];
      }
      assert(k >= 1);
   }
   assert(ngubconsGOC1 + ngubconsGFC1 + ngubconsGC2 + ngubconsGR == ngubconss - ngubconscapexceed);
   assert(ngubconsGOC1 + ngubconsGNC1 == ngubconsGC1);

   /* initialize the minweight tables, defined as: for i = 1,...,m with m = |I| and w = 0,...,|gubconsGC1|;
    * - finished_i[w] =
    *   min   sum_{k = 1,2,...,i-1} sum_{j in Q_k} a_j x_j
    *   s.t.  sum_{k = 1,2,...,i-1} sum_{j in Q_k} alpha_j x_j  >= w
    *                               sum_{j in Q_k} x_j <= 1
    *                               x_j in {0,1} forall j in Q_k forall k = 1,2,...,i-1,
    * - unfinished_i[w] =
    *   min   sum_{k = i+1,...,m} sum_{j in Q_k && j in C1} a_j x_j
    *   s.t.  sum_{k = i+1,...,m} sum_{j in Q_k && j in C1} x_j  >= w
    *                             sum_{j in Q_k} x_j <= 1
    *                             x_j in {0,1} forall j in Q_k forall k = 1,2,...,i-1,
    * - minweights_i[w] = min{finished_i[w1] + unfinished_i[w2] : w1>=0, w2>=0, w1+w2=w};
    */

   /* initialize finished table; note that variables in GOC1 GUBs (includes C1 and capacity exceeding variables)
    * are sorted s.t. C1 variables come first and are sorted by non-decreasing weight.
    * GUBs in the group GCI are sorted by non-decreasing min{ a_k : k in GC1_j } where min{ a_k : k in GC1_j } always
    * comes from the first variable in the GUB
    */
   assert(ngubconsGOC1 <= ngubconsGC1);
   finished[0] = 0;
   for( w = 1; w <= ngubconsGOC1; w++ )
   {
      liftgubconsidx = gubconsGOC1[w-1];

      assert(gubset->gubconsstatus[liftgubconsidx] == GUBCONSSTATUS_BELONGSTOSET_GOC1);
      assert(gubset->gubconss[liftgubconsidx]->gubvarsstatus[0] == GUBVARSTATUS_BELONGSTOSET_C1);

      varidx = gubset->gubconss[liftgubconsidx]->gubvars[0];

      assert(varidx >= 0 && varidx < nvars);
      assert(liftcoefs[varidx] == 1);

      min = weights[varidx];
      finished[w] = finished[w-1] + min;

#ifndef NDEBUG
      for( k = 1; k < gubset->gubconss[liftgubconsidx]->ngubvars
              && gubset->gubconss[liftgubconsidx]->gubvarsstatus[k] == GUBVARSTATUS_BELONGSTOSET_C1; k++ )
      {
         varidx = gubset->gubconss[liftgubconsidx]->gubvars[k];
         assert(varidx >= 0 && varidx < nvars);
         assert(liftcoefs[varidx] == 1);
         assert(weights[varidx] >= min);
      }
#endif
   }
   for( w = ngubconsGOC1+1; w <= ngubconsGC1; w++ )
      finished[w] = SCIP_LONGINT_MAX;

   /* initialize unfinished table; note that variables in GNC1 GUBs
    * are sorted s.t. C1 variables come first and are sorted by non-decreasing weight.
    * GUBs in the group GCI are sorted by non-decreasing min{ a_k : k in GC1_j } where min{ a_k : k in GC1_j } always
    * comes from the first variable in the GUB
    */
   assert(ngubconsGNC1 <= ngubconsGC1);
   unfinished[0] = 0;
   for( w = 1; w <= ngubconsGNC1; w++ )
   {
      liftgubconsidx = gubconsGNC1[w-1];

      assert(gubset->gubconsstatus[liftgubconsidx] == GUBCONSSTATUS_BELONGSTOSET_GNC1);
      assert(gubset->gubconss[liftgubconsidx]->gubvarsstatus[0] == GUBVARSTATUS_BELONGSTOSET_C1);

      varidx = gubset->gubconss[liftgubconsidx]->gubvars[0];

      assert(varidx >= 0 && varidx < nvars);
      assert(liftcoefs[varidx] == 1);

      min = weights[varidx];
      unfinished[w] = unfinished[w-1] + min;

#ifndef NDEBUG
      for( k = 1; k < gubset->gubconss[liftgubconsidx]->ngubvars
              && gubset->gubconss[liftgubconsidx]->gubvarsstatus[k] == GUBVARSTATUS_BELONGSTOSET_C1; k++ )
      {
         varidx = gubset->gubconss[liftgubconsidx]->gubvars[k];
         assert(varidx >= 0 && varidx < nvars);
         assert(liftcoefs[varidx] == 1);
         assert(weights[varidx] >= min );
      }
#endif
   }
   for( w = ngubconsGNC1 + 1; w <= ngubconsGC1; w++ )
      unfinished[w] = SCIP_LONGINT_MAX;

   /* initialize minweights table; note that variables in GC1 GUBs
    * are sorted s.t. C1 variables come first and are sorted by non-decreasing weight.
    * we can directly initialize minweights instead of computing it from finished and unfinished (which would be more time
    * consuming) because is it has to be build using weights from C1 only.
    */
   assert(ngubconsGOC1 + ngubconsGNC1 == ngubconsGC1);
   minweights[0] = 0;
   for( w = 1; w <= ngubconsGC1; w++ )
   {
      liftgubconsidx = gubconsGC1[w-1];

      assert(gubset->gubconsstatus[liftgubconsidx] == GUBCONSSTATUS_BELONGSTOSET_GOC1
          || gubset->gubconsstatus[liftgubconsidx] == GUBCONSSTATUS_BELONGSTOSET_GNC1);
      assert(gubset->gubconss[liftgubconsidx]->gubvarsstatus[0] == GUBVARSTATUS_BELONGSTOSET_C1);

      varidx = gubset->gubconss[liftgubconsidx]->gubvars[0];

      assert(varidx >= 0 && varidx < nvars);
      assert(liftcoefs[varidx] == 1);

      min = weights[varidx];
      minweights[w] = minweights[w-1] + min;

#ifndef NDEBUG
      for( k = 1; k < gubset->gubconss[liftgubconsidx]->ngubvars
              && gubset->gubconss[liftgubconsidx]->gubvarsstatus[k] == GUBVARSTATUS_BELONGSTOSET_C1; k++ )
      {
         varidx = gubset->gubconss[liftgubconsidx]->gubvars[k];
         assert(varidx >= 0 && varidx < nvars);
         assert(liftcoefs[varidx] == 1);
         assert(weights[varidx] >= min);
      }
#endif
   }
   minweightslen = ngubconsGC1 + 1;

   /* gets sum of weights of variables fixed to one, i.e. sum of weights of C2 variables GC2 GUBs */
   fixedonesweight = 0;
   for( j = 0; j < ngubconsGC2; j++ )
   {
      varidx = gubset->gubconss[gubconsGC2[j]]->gubvars[0];

      assert(gubset->gubconss[gubconsGC2[j]]->ngubvars == 1);
      assert(varidx >= 0 && varidx < nvars);
      assert(gubset->gubconss[gubconsGC2[j]]->gubvarsstatus[0] == GUBVARSTATUS_BELONGSTOSET_C2);

      fixedonesweight += weights[varidx];
   }
   assert(fixedonesweight >= 0);

   /* initializes right hand side of lifted valid inequality */
   *liftrhs = alpha0;

   /* sequentially up-lifts all variables in GFC1 GUBs */
   for( j = 0; j < ngubconsGFC1; j++ )
   {
      liftgubconsidx = gubconsGFC1[j];
      assert(liftgubconsidx >= 0 && liftgubconsidx < ngubconss);

      /* GNC1 GUB: update unfinished table (remove current GUB, i.e., remove min weight of C1 vars in GUB) and
       * compute minweight table via updated unfinished table and aleady upto date finished table;
       */
      k = 0;
      if( gubset->gubconsstatus[liftgubconsidx] == GUBCONSSTATUS_BELONGSTOSET_GNC1 )
      {
         assert(gubset->gubconsstatus[liftgubconsidx] == GUBCONSSTATUS_BELONGSTOSET_GNC1);
         assert(gubset->gubconss[liftgubconsidx]->gubvarsstatus[0] == GUBVARSTATUS_BELONGSTOSET_C1);
         assert(ngubconsGNC1 > 0);

         /* get number of C1 variables of current GNC1 GUB and put them into array of variables in GUB that
          * are considered for the lifting, i.e., not capacity exceeding
          */
         for( ; k < gubset->gubconss[liftgubconsidx]->ngubvars
                && gubset->gubconss[liftgubconsidx]->gubvarsstatus[k] == GUBVARSTATUS_BELONGSTOSET_C1; k++ )
            liftgubvars[k] = gubset->gubconss[liftgubconsidx]->gubvars[k];
         assert(k >= 1);

         /* update unfinished table by removing current GNC1 GUB, i.e, remove C1 variable with minimal weight
          * unfinished[w] = MAX{unfinished[w], unfinished[w+1] - weight}, "weight" is the minimal weight of current GUB
          */
         weight = weights[liftgubvars[0]];

         weightdiff2 = unfinished[ngubconsGNC1] - weight;
         unfinished[ngubconsGNC1] = SCIP_LONGINT_MAX;
         for( w = ngubconsGNC1-1; w >= 1; w-- )
         {
            weightdiff1 = weightdiff2;
            weightdiff2 = unfinished[w] - weight;

            if( unfinished[w] < weightdiff1 )
               unfinished[w] = weightdiff1;
            else
               break;
         }
         ngubconsGNC1--;

         /* computes minweights table by combining unfished and fished tables */
         computeMinweightsGUB(minweights, finished, unfinished, minweightslen);
         assert(minweights[0] == 0);
      }
      /* GF GUB: no update of unfinished table (and minweight table) required because GF GUBs have no C1 variables and
       * are therefore not in the unfinished table
       */
      else
         assert(gubset->gubconsstatus[liftgubconsidx] == GUBCONSSTATUS_BELONGSTOSET_GF);

#ifndef NDEBUG
      nliftgubC1 = k;
#endif
      nliftgubvars = k;
      sumliftcoef = 0;

      /* compute lifting coefficient of F and R variables in GNC1 and GF GUBs (C1 vars have already liftcoef 1) */
      for( ; k < gubset->gubconss[liftgubconsidx]->ngubvars; k++ )
      {
         if( gubset->gubconss[liftgubconsidx]->gubvarsstatus[k] == GUBVARSTATUS_BELONGSTOSET_F
             || gubset->gubconss[liftgubconsidx]->gubvarsstatus[k] == GUBVARSTATUS_BELONGSTOSET_R )
         {
            liftvar = gubset->gubconss[liftgubconsidx]->gubvars[k];
            weight = weights[liftvar];
            assert(weight > 0);
            assert(liftvar >= 0 && liftvar < nvars);
            assert(capacity - weight >= 0);

            /* put variable into array of variables in GUB that are considered for the lifting,
             *  i.e., not capacity exceeding
             */
            liftgubvars[nliftgubvars] = liftvar;
            nliftgubvars++;

            /* knapsack problem is infeasible:
             * sets z = 0
             */
            if( capacity - fixedonesweight - weight < 0 )
            {
               z = 0;
            }
            /* knapsack problem is feasible:
             *   sets z = max { w : 0 <= w <= liftrhs, minweights_i[w] <= a_0 - fixedonesweight - a_{j_i}  } = liftrhs,
             *   if minweights_i[liftrhs] <= a_0 - fixedonesweight - a_{j_i}
             */
            else if( minweights[*liftrhs] <= capacity - fixedonesweight - weight )
            {
               z = *liftrhs;
            }
            /* knapsack problem is feasible:
             *   binary search to find z = max {w : 0 <= w <= liftrhs, minweights_i[w] <= a_0 - fixedonesweight - a_{j_i}}
             */
            else
            {
               assert((*liftrhs) + 1 >= minweightslen || minweights[(*liftrhs) + 1] > capacity - fixedonesweight - weight);
               left = 0;
               right = (*liftrhs) + 1;
               while( left < right - 1 )
               {
                  middle = (left + right) / 2;
                  assert(0 <= middle && middle < minweightslen);
                  if( minweights[middle] <= capacity - fixedonesweight - weight )
                     left = middle;
                  else
                     right = middle;
               }
               assert(left == right - 1);
               assert(0 <= left && left < minweightslen);
               assert(minweights[left] <= capacity - fixedonesweight - weight);
               assert(left == minweightslen - 1 || minweights[left+1] > capacity - fixedonesweight - weight);

               /* now z = left */
               z = left;
               assert(z <= *liftrhs);
            }

            /* calculates lifting coefficients alpha_{j_i} = liftrhs - z */
            liftcoef = (*liftrhs) - z;
            liftcoefs[liftvar] = liftcoef;
            assert(liftcoef >= 0 && liftcoef <= (*liftrhs) + 1);

            /* updates activity of current valid inequality */
            (*cutact) += liftcoef * solvals[liftvar];

            /* updates sum of all lifting coefficients in GUB */
            sumliftcoef += liftcoefs[liftvar];
         }
         else
            assert(gubset->gubconss[liftgubconsidx]->gubvarsstatus[k] == GUBVARSTATUS_CAPACITYEXCEEDED);
      }
      /* at least one variable is in F or R (j = number of C1 variables in current GUB) */
      assert(nliftgubvars > nliftgubC1);

      /* activity of current valid inequality will not change if (sum of alpha_{j_i} in GUB) = 0
       * and finished and minweight table can be updated easily as only C1 variables need to be considered;
       * not needed for GF GUBs
       */
      if( sumliftcoef == 0 )
      {
         if( gubset->gubconsstatus[liftgubconsidx] == GUBCONSSTATUS_BELONGSTOSET_GNC1 )
         {
            weight = weights[liftgubvars[0]];
            /* update finished table and minweights table by applying special case of
             * finished[w] = MIN{finished[w], finished[w-1] + weight}, "weight" is the minimal weight of current GUB
             * minweights[w] = MIN{minweights[w], minweights[w-1] + weight}, "weight" is the minimal weight of current GUB
             */
            for( w = minweightslen-1; w >= 1; w-- )
            {
               SCIP_Longint tmpval;

               tmpval = safeAddMinweightsGUB(finished[w-1], weight);
               finished[w] = MIN(finished[w], tmpval);

               tmpval = safeAddMinweightsGUB(minweights[w-1], weight);
               minweights[w] = MIN(minweights[w], tmpval);
            }
         }
         else
            assert(gubset->gubconsstatus[liftgubconsidx] == GUBCONSSTATUS_BELONGSTOSET_GF);

         continue;
      }

      /* enlarges current minweights tables(finished, unfinished, minweights):
       *  from minweightlen = |gubconsGC1| + sum_{k=1,2,...,i-1}sum_{j in Q_k} alpha_j + 1 entries
       *  to                  |gubconsGC1| + sum_{k=1,2,...,i  }sum_{j in Q_k} alpha_j + 1 entries
       *  and sets minweights_i[w] = infinity for
       *  w = |gubconsGC1| + sum_{k=1,2,..,i-1}sum_{j in Q_k} alpha_j+1,..,|C1| + sum_{k=1,2,..,i}sum_{j in Q_k} alpha_j
       */
      tmplen = minweightslen; /* will be updated in enlargeMinweights() */
      tmpsize = minweightssize;
      SCIP_CALL( enlargeMinweights(scip, &unfinished, &tmplen, &tmpsize, tmplen + sumliftcoef) );
      tmplen = minweightslen;
      tmpsize = minweightssize;
      SCIP_CALL( enlargeMinweights(scip, &finished, &tmplen, &tmpsize, tmplen + sumliftcoef) );
      SCIP_CALL( enlargeMinweights(scip, &minweights, &minweightslen, &minweightssize, minweightslen + sumliftcoef) );

      /* update finished table and minweight table;
       * note that instead of computing minweight table from updated finished and updated unfinished table again
       * (for the lifting coefficient, we had to update unfinished table and compute minweight table), we here
       * only need to update the minweight table and the updated finished in the same way (i.e., computing for minweight
       * not needed because only finished table changed at this point and the change was "adding" one weight)
       *
       * update formular for minweight table is: minweight_i+1[w] =
       *   min{ minweights_i[w], min{ minweights_i[w - alpha_k]^{+} + a_k : k in GUB_j_i } }
       * formular for finished table has the same pattern.
       */
      for( w = minweightslen-1; w >= 0; w-- )
      {
         SCIP_Longint minminweight;
         SCIP_Longint minfinished;

         for( k = 0; k < nliftgubvars; k++ )
         {
            liftcoef = liftcoefs[liftgubvars[k]];
            weight = weights[liftgubvars[k]];

            if( w < liftcoef )
            {
               minfinished = MIN(finished[w], weight);
               minminweight = MIN(minweights[w], weight);

               finished[w] = minfinished;
               minweights[w] = minminweight;
            }
            else
            {
               SCIP_Longint tmpval;

               assert(w >= liftcoef);

               tmpval = safeAddMinweightsGUB(finished[w-liftcoef], weight);
               minfinished = MIN(finished[w], tmpval);

               tmpval = safeAddMinweightsGUB(minweights[w-liftcoef], weight);
               minminweight = MIN(minweights[w], tmpval);

               finished[w] = minfinished;
               minweights[w] = minminweight;
            }
         }
      }
      assert(minweights[0] == 0);
   }
   assert(ngubconsGNC1 == 0);

   /* note: now the unfinished table no longer exists, i.e., it is "0, MAX, MAX, ..." and minweight equals to finished;
    * therefore, only work with minweight table from here on
    */

   /* sequentially down-lifts C2 variables contained in trivial GC2 GUBs */
   for( j = 0; j < ngubconsGC2; j++ )
   {
      liftgubconsidx = gubconsGC2[j];

      assert(liftgubconsidx >=0 && liftgubconsidx < ngubconss);
      assert(gubset->gubconsstatus[liftgubconsidx] == GUBCONSSTATUS_BELONGSTOSET_GC2);
      assert(gubset->gubconss[liftgubconsidx]->ngubvars == 1);
      assert(gubset->gubconss[liftgubconsidx]->gubvarsstatus[0] == GUBVARSTATUS_BELONGSTOSET_C2);

      liftvar = gubset->gubconss[liftgubconsidx]->gubvars[0]; /* C2 GUBs contain only one variable */
      weight = weights[liftvar];

      assert(liftvar >= 0 && liftvar < nvars);
      assert(SCIPisFeasEQ(scip, solvals[liftvar], 1.0));
      assert(weight > 0);

      /* uses binary search to find
       *   z = max { w : 0 <= w <= |C_1| + sum_{k=1}^{i-1} alpha_{j_k}, minweights_[w] <= a_0 - fixedonesweight + a_{j_i}}
       */
      left = 0;
      right = minweightslen;
      while( left < right - 1 )
      {
         middle = (left + right) / 2;
         assert(0 <= middle && middle < minweightslen);
         if( minweights[middle] <= capacity - fixedonesweight + weight )
            left = middle;
         else
            right = middle;
      }
      assert(left == right - 1);
      assert(0 <= left && left < minweightslen);
      assert(minweights[left] <= capacity - fixedonesweight + weight);
      assert(left == minweightslen - 1 || minweights[left + 1] > capacity - fixedonesweight + weight);

      /* now z = left */
      z = left;
      assert(z >= *liftrhs);

      /* calculates lifting coefficients alpha_{j_i} = z - liftrhs */
      liftcoef = z - (*liftrhs);
      liftcoefs[liftvar] = liftcoef;
      assert(liftcoef >= 0);

      /* updates sum of weights of variables fixed to one */
      fixedonesweight -= weight;

      /* updates right-hand side of current valid inequality */
      (*liftrhs) += liftcoef;
      assert(*liftrhs >= alpha0);

      /* minweight table and activity of current valid inequality will not change, if alpha_{j_i} = 0 */
      if( liftcoef == 0 )
         continue;

      /* updates activity of current valid inequality */
      (*cutact) += liftcoef * solvals[liftvar];

      /* enlarges current minweight table:
       *  from minweightlen = |gubconsGC1| + sum_{k=1,2,...,i-1}sum_{j in Q_k} alpha_j + 1 entries
       *  to                  |gubconsGC1| + sum_{k=1,2,...,i  }sum_{j in Q_k} alpha_j + 1 entries
       * and sets minweights_i[w] = infinity for
       *  w = |C1| + sum_{k=1,2,...,i-1}sum_{j in Q_k} alpha_j + 1 , ... , |C1| + sum_{k=1,2,...,i}sum_{j in Q_k} alpha_j
       */
      SCIP_CALL( enlargeMinweights(scip, &minweights, &minweightslen, &minweightssize, minweightslen + liftcoef) );

      /* updates minweight table: minweight_i+1[w] =
       *  min{ minweights_i[w], a_{j_i}},                                 if w <  alpha_j_i
       *  min{ minweights_i[w], minweights_i[w - alpha_j_i] + a_j_i},     if w >= alpha_j_i
       */
      for( w = minweightslen - 1; w >= 0; w-- )
      {
         if( w < liftcoef )
         {
            min = MIN(minweights[w], weight);
            minweights[w] = min;
         }
         else
         {
            SCIP_Longint tmpval;

            assert(w >= liftcoef);

            tmpval = safeAddMinweightsGUB(minweights[w-liftcoef], weight);
            min = MIN(minweights[w], tmpval);
            minweights[w] = min;
         }
      }
   }
   assert(fixedonesweight == 0);
   assert(*liftrhs >= alpha0);

   /* sequentially up-lifts variables in GUB constraints in GR GUBs */
   for( j = 0; j < ngubconsGR; j++ )
   {
      liftgubconsidx = gubconsGR[j];

      assert(liftgubconsidx >=0 && liftgubconsidx < ngubconss);
      assert(gubset->gubconsstatus[liftgubconsidx] == GUBCONSSTATUS_BELONGSTOSET_GR);

      sumliftcoef = 0;
      nliftgubvars = 0;
      for( k = 0; k < gubset->gubconss[liftgubconsidx]->ngubvars; k++ )
      {
         if(gubset->gubconss[liftgubconsidx]->gubvarsstatus[k] == GUBVARSTATUS_BELONGSTOSET_R )
         {
            liftvar = gubset->gubconss[liftgubconsidx]->gubvars[k];
            weight = weights[liftvar];
            assert(weight > 0);
            assert(liftvar >= 0 && liftvar < nvars);
            assert(capacity - weight >= 0);
            assert((*liftrhs) + 1 >= minweightslen || minweights[(*liftrhs) + 1] > capacity - weight);

            /* put variable into array of variables in GUB that are considered for the lifting,
             *  i.e., not capacity exceeding
             */
            liftgubvars[nliftgubvars] = liftvar;
            nliftgubvars++;

            /* sets z = max { w : 0 <= w <= liftrhs, minweights_i[w] <= a_0 - a_{j_i} } = liftrhs,
             * if minweights_i[liftrhs] <= a_0 - a_{j_i}
             */
            if( minweights[*liftrhs] <= capacity - weight )
            {
               z = *liftrhs;
            }
            /* uses binary search to find z = max { w : 0 <= w <= liftrhs, minweights_i[w] <= a_0 - a_{j_i} }
             */
            else
            {
               left = 0;
               right = (*liftrhs) + 1;
               while( left < right - 1 )
               {
                  middle = (left + right) / 2;
                  assert(0 <= middle && middle < minweightslen);
                  if( minweights[middle] <= capacity - weight )
                     left = middle;
                  else
                     right = middle;
               }
               assert(left == right - 1);
               assert(0 <= left && left < minweightslen);
               assert(minweights[left] <= capacity - weight);
               assert(left == minweightslen - 1 || minweights[left + 1] > capacity - weight);

               /* now z = left */
               z = left;
               assert(z <= *liftrhs);
            }
            /* calculates lifting coefficients alpha_{j_i} = liftrhs - z */
            liftcoef = (*liftrhs) - z;
            liftcoefs[liftvar] = liftcoef;
            assert(liftcoef >= 0 && liftcoef <= (*liftrhs) + 1);

            /* updates activity of current valid inequality */
            (*cutact) += liftcoef * solvals[liftvar];

            /* updates sum of all lifting coefficients in GUB */
            sumliftcoef += liftcoefs[liftvar];
         }
         else
            assert(gubset->gubconss[liftgubconsidx]->gubvarsstatus[k] == GUBVARSTATUS_CAPACITYEXCEEDED);
      }
      assert(nliftgubvars >= 1); /* at least one variable is in R */

      /* minweight table and activity of current valid inequality will not change if (sum of alpha_{j_i} in GUB) = 0 */
      if( sumliftcoef == 0 )
         continue;

      /* updates minweight table: minweight_i+1[w] =
       *   min{ minweights_i[w], min{ minweights_i[w - alpha_k]^{+} + a_k : k in GUB_j_i } }
       */
      for( w = *liftrhs; w >= 0; w-- )
      {
         for( k = 0; k < nliftgubvars; k++ )
         {
            liftcoef = liftcoefs[liftgubvars[k]];
            weight = weights[liftgubvars[k]];

            if( w < liftcoef )
            {
               min = MIN(minweights[w], weight);
               minweights[w] = min;
            }
            else
            {
               SCIP_Longint tmpval;

               assert(w >= liftcoef);

               tmpval = safeAddMinweightsGUB(minweights[w-liftcoef], weight);
               min = MIN(minweights[w], tmpval);
               minweights[w] = min;
            }
         }
      }
      assert(minweights[0] == 0);
   }

   /* frees temporary memory */
   SCIPfreeBufferArray(scip, &minweights);
   SCIPfreeBufferArray(scip, &finished);
   SCIPfreeBufferArray(scip, &unfinished);
   SCIPfreeBufferArray(scip, &liftgubvars);
   SCIPfreeBufferArray(scip, &gubconsGOC1 );
   SCIPfreeBufferArray(scip, &gubconsGNC1);

   return SCIP_OKAY;
}

/** lifts given minimal cover inequality 
 *  \f[
 *    \sum_{j \in C} x_j \leq |C| - 1 
 *  \f]
 *  valid for 
 *  \f[
 *    S^0 = \{ x \in {0,1}^{|C|} : \sum_{j \in C} a_j x_j \leq a_0 \}
 *  \f]
 *  to a valid inequality 
 *  \f[
 *    \sum_{j \in C} x_j + \sum_{j \in N \setminus C} \alpha_j x_j \leq |C| - 1
 *  \f]
 *  for
 *  \f[ 
 *    S = \{ x \in {0,1}^{|N|} : \sum_{j \in N} a_j x_j \leq a_0 \}; 
 *  \f]
 *  uses superadditive up-lifting for the variables in \f$N \setminus C\f$.
 */ 
static
SCIP_RETCODE superadditiveUpLifting(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_VAR**            vars,               /**< variables in knapsack constraint */
   int                   nvars,              /**< number of variables in knapsack constraint */
   int                   ntightened,         /**< number of variables with tightened upper bound */
   SCIP_Longint*         weights,            /**< weights of variables in knapsack constraint */
   SCIP_Longint          capacity,           /**< capacity of knapsack */
   SCIP_Real*            solvals,            /**< solution values of all problem variables */
   int*                  covervars,          /**< cover variables */
   int*                  noncovervars,       /**< noncover variables */
   int                   ncovervars,         /**< number of cover variables */
   int                   nnoncovervars,      /**< number of noncover variables */
   SCIP_Longint          coverweight,        /**< weight of cover */
   SCIP_Real*            liftcoefs,          /**< pointer to store lifting coefficient of vars in knapsack constraint */
   SCIP_Real*            cutact              /**< pointer to store activity of lifted valid inequality */
   )
{
   SCIP_Longint* maxweightsums;
   SCIP_Longint* intervalends;
   SCIP_Longint* rhos;
   SCIP_Real* sortkeys;
   SCIP_Longint lambda;
   int j;
   int h;

   assert(scip != NULL);
   assert(vars != NULL);
   assert(nvars >= 0);
   assert(weights != NULL);
   assert(capacity >= 0);
   assert(solvals != NULL);
   assert(covervars != NULL);
   assert(noncovervars != NULL);
   assert(ncovervars > 0 && ncovervars <= nvars);
   assert(nnoncovervars >= 0 && nnoncovervars <= nvars - ntightened);
   assert(ncovervars + nnoncovervars == nvars - ntightened);
   assert(liftcoefs != NULL);
   assert(cutact != NULL);

   /* allocates temporary memory */
   SCIP_CALL( SCIPallocBufferArray(scip, &sortkeys, ncovervars) );
   SCIP_CALL( SCIPallocBufferArray(scip, &maxweightsums, ncovervars + 1) );
   SCIP_CALL( SCIPallocBufferArray(scip, &intervalends, ncovervars) );
   SCIP_CALL( SCIPallocBufferArray(scip, &rhos, ncovervars) );

   /* initializes data structures */
   BMSclearMemoryArray(liftcoefs, nvars);
   *cutact = 0.0;

   /* sets lifting coefficient of variables in C, sorts variables in C such that a_1 >= a_2 >= ... >= a_|C| 
    * and calculates activity of current valid inequality 
    */
   for( j = 0; j < ncovervars; j++ )
   {
      assert(liftcoefs[covervars[j]] == 0.0);
      liftcoefs[covervars[j]] = 1.0;
      sortkeys[j] = (SCIP_Real) weights[covervars[j]];
      (*cutact) += solvals[covervars[j]];
   }
   SCIPsortDownRealInt(sortkeys, covervars, ncovervars);

   /* calculates weight excess of cover C */
   lambda = coverweight - capacity;
   assert(lambda > 0);

   /* calculates A_h for h = 0,...,|C|, I_h for h = 1,...,|C| and rho_h for h = 1,...,|C| */
   maxweightsums[0] = 0;
   for( h = 1; h <= ncovervars; h++ )
   {
      maxweightsums[h] = maxweightsums[h-1] + weights[covervars[h-1]];
      intervalends[h-1] = maxweightsums[h] - lambda;
      rhos[h-1] = MAX(0, weights[covervars[h-1]] - weights[covervars[0]] + lambda);
   }

   /* sorts variables in N\C such that a_{j_1} <= a_{j_2} <= ... <= a_{j_t} */
   for( j = 0; j < nnoncovervars; j++ )
      sortkeys[j] = (SCIP_Real) (weights[noncovervars[j]]);
   SCIPsortRealInt(sortkeys, noncovervars, nnoncovervars);

   /* calculates lifting coefficient for all variables in N\C */
   h = 0;
   for( j = 0; j < nnoncovervars; j++ )
   {
      int liftvar;
      SCIP_Longint weight;
      SCIP_Real liftcoef;

      liftvar = noncovervars[j];
      weight = weights[liftvar];

      while( intervalends[h] < weight )
         h++;

      if( h == 0 )
         liftcoef = h;
      else 
      {
         if( weight <= intervalends[h-1] + rhos[h] )
         {
            SCIP_Real tmp1;
            SCIP_Real tmp2;
            tmp1 =  (SCIP_Real) (intervalends[h-1] + rhos[h] - weight);
            tmp2 =  (SCIP_Real) rhos[1];
            liftcoef = h - ( tmp1 / tmp2 );
         }
         else
            liftcoef = h;
      }      

      /* sets lifting coefficient */
      assert(liftcoefs[liftvar] == 0.0);
      liftcoefs[liftvar] = liftcoef;

      /* updates activity of current valid inequality */
      (*cutact) += liftcoef * solvals[liftvar];
   } 

   /* frees temporary memory */
   SCIPfreeBufferArray(scip, &rhos);
   SCIPfreeBufferArray(scip, &intervalends);
   SCIPfreeBufferArray(scip, &maxweightsums);
   SCIPfreeBufferArray(scip, &sortkeys);

   return SCIP_OKAY;
}


/** separates lifted minimal cover inequalities using sequential up- and down-lifting and GUB information, if wanted, for
 *  given knapsack problem
*/
static
SCIP_RETCODE separateSequLiftedMinimalCoverInequality(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONS*            cons,               /**< originating constraint of the knapsack problem, or NULL */
   SCIP_SEPA*            sepa,               /**< originating separator of the knapsack problem, or NULL */
   SCIP_VAR**            vars,               /**< variables in knapsack constraint */
   int                   nvars,              /**< number of variables in knapsack constraint */
   int                   ntightened,         /**< number of variables with tightened upper bound */
   SCIP_Longint*         weights,            /**< weights of variables in knapsack constraint */
   SCIP_Longint          capacity,           /**< capacity of knapsack */
   SCIP_Real*            solvals,            /**< solution values of all problem variables */
   int*                  mincovervars,       /**< mincover variables */
   int*                  nonmincovervars,    /**< nonmincover variables */
   int                   nmincovervars,      /**< number of mincover variables */
   int                   nnonmincovervars,   /**< number of nonmincover variables */
   SCIP_SOL*             sol,                /**< primal SCIP solution to separate, NULL for current LP solution */
   SCIP_GUBSET*          gubset,             /**< GUB set data structure, NULL if no GUB information should be used */
   SCIP_Bool*            cutoff,             /**< pointer to store whether a cutoff has been detected */
   int*                  ncuts               /**< pointer to add up the number of found cuts */
   )
{
   int* varsC1;
   int* varsC2;
   int* varsF;
   int* varsR;
   int nvarsC1;
   int nvarsC2;
   int nvarsF;
   int nvarsR;
   SCIP_Real cutact;
   int* liftcoefs;
   int liftrhs;

   assert( cutoff != NULL );
   *cutoff = FALSE;

   /* allocates temporary memory */
   SCIP_CALL( SCIPallocBufferArray(scip, &varsC1, nvars) );
   SCIP_CALL( SCIPallocBufferArray(scip, &varsC2, nvars) );
   SCIP_CALL( SCIPallocBufferArray(scip, &varsF, nvars) );
   SCIP_CALL( SCIPallocBufferArray(scip, &varsR, nvars) );
   SCIP_CALL( SCIPallocBufferArray(scip, &liftcoefs, nvars) );

   /* gets partition (C_1,C_2) of C, i.e. C_1 & C_2 = C and C_1 cap C_2 = emptyset, with C_1 not empty; chooses partition
    * as follows
    *   C_2 = { j in C : x*_j = 1 } and
    *   C_1 = C\C_2
    */
   getPartitionCovervars(scip, solvals, mincovervars, nmincovervars, varsC1, varsC2, &nvarsC1, &nvarsC2);
   assert(nvarsC1 + nvarsC2 == nmincovervars);
   assert(nmincovervars > 0);
   assert(nvarsC1 >= 0); /* nvarsC1 > 0 does not always hold, because relaxed knapsack conss may already be violated */

   /* changes partition (C_1,C_2) of minimal cover C, if |C1| = 1, by moving one variable from C2 to C1 */
   if( nvarsC1 < 2 && nvarsC2 > 0)
   {
      SCIP_CALL( changePartitionCovervars(scip, weights, varsC1, varsC2, &nvarsC1, &nvarsC2) );
      assert(nvarsC1 >= 1);
   }
   assert(nvarsC2 == 0 || nvarsC1 >= 1);

   /* gets partition (F,R) of N\C, i.e. F & R = N\C and F cap R = emptyset; chooses partition as follows
    *   R = { j in N\C : x*_j = 0 } and
    *   F = (N\C)\F
    */
   getPartitionNoncovervars(scip, solvals, nonmincovervars, nnonmincovervars, varsF, varsR, &nvarsF, &nvarsR);
   assert(nvarsF + nvarsR == nnonmincovervars);
   assert(nvarsC1 + nvarsC2 + nvarsF + nvarsR == nvars - ntightened);

   /* lift cuts without GUB information */
   if( gubset == NULL )
   {
      /* sorts variables in F, C_2, R according to the second level lifting sequence that will be used in the sequential
       * lifting procedure
       */
      SCIP_CALL( getLiftingSequence(scip, solvals, weights, varsF, varsC2, varsR, nvarsF, nvarsC2, nvarsR) );

      /* lifts minimal cover inequality sum_{j in C_1} x_j <= |C_1| - 1 valid for
       *
       *    S^0 = { x in {0,1}^|C_1| : sum_{j in C_1} a_j x_j <= a_0 - sum_{j in C_2} a_j }
       *
       * to a valid inequality sum_{j in C_1} x_j + sum_{j in N\C_1} alpha_j x_j <= |C_1| - 1 + sum_{j in C_2} alpha_j for
       *
       *      S = { x in {0,1}^|N|   : sum_{j in N}   a_j x_j <= a_0 },
       *
       * uses sequential up-lifting for the variables in F, sequential down-lifting for the variable in C_2 and sequential
       * up-lifting for the variables in R according to the second level lifting sequence
       */
      SCIP_CALL( sequentialUpAndDownLifting(scip, vars, nvars, ntightened, weights, capacity, solvals, varsC1, varsC2,
            varsF, varsR, nvarsC1, nvarsC2, nvarsF, nvarsR, nvarsC1 - 1, liftcoefs, &cutact, &liftrhs) );
   }
   /* lift cuts with GUB information */
   else
   {
      int* gubconsGC1;
      int* gubconsGC2;
      int* gubconsGFC1;
      int* gubconsGR;
      int ngubconsGC1;
      int ngubconsGC2;
      int ngubconsGFC1;
      int ngubconsGR;
      int ngubconss;
      int nconstightened;
      int maxgubvarssize;

      assert(nvars == gubset->nvars);

      ngubconsGC1 = 0;
      ngubconsGC2 = 0;
      ngubconsGFC1 = 0;
      ngubconsGR = 0;
      ngubconss = gubset->ngubconss;
      nconstightened = 0;
      maxgubvarssize = 0;

      /* allocates temporary memory */
      SCIP_CALL( SCIPallocBufferArray(scip, &gubconsGC1, ngubconss) );
      SCIP_CALL( SCIPallocBufferArray(scip, &gubconsGC2, ngubconss) );
      SCIP_CALL( SCIPallocBufferArray(scip, &gubconsGFC1, ngubconss) );
      SCIP_CALL( SCIPallocBufferArray(scip, &gubconsGR, ngubconss) );

      /* categorizies GUBs of knapsack GUB partion into GOC1, GNC1, GF, GC2, and GR and computes a lifting sequence of
       * the GUBs for the sequential GUB wise lifting procedure
       */
      SCIP_CALL( getLiftingSequenceGUB(scip, gubset, solvals, weights, varsC1, varsC2, varsF, varsR, nvarsC1,
            nvarsC2, nvarsF, nvarsR, gubconsGC1, gubconsGC2, gubconsGFC1, gubconsGR, &ngubconsGC1, &ngubconsGC2,
            &ngubconsGFC1, &ngubconsGR, &nconstightened, &maxgubvarssize) );

      /* lifts minimal cover inequality sum_{j in C_1} x_j <= |C_1| - 1 valid for
       *
       *    S^0 = { x in {0,1}^|C_1| : sum_{j in C_1} a_j x_j <= a_0 - sum_{j in C_2} a_j,
       *                               sum_{j in Q_i} x_j <= 1, forall i in I }
       *
       * to a valid inequality sum_{j in C_1} x_j + sum_{j in N\C_1} alpha_j x_j <= |C_1| - 1 + sum_{j in C_2} alpha_j for
       *
       *      S = { x in {0,1}^|N|   : sum_{j in N}   a_j x_j <= a_0, sum_{j in Q_i} x_j <= 1, forall i in I },
       *
       *  uses sequential up-lifting   for the variables in GUB constraints in gubconsGFC1,
       *       sequential down-lifting for the variables in GUB constraints in gubconsGC2, and
       *       sequential up-lifting   for the variabels in GUB constraints in gubconsGR.
       */
      SCIP_CALL( sequentialUpAndDownLiftingGUB(scip, gubset, vars, nconstightened, weights, capacity, solvals, gubconsGC1,
            gubconsGC2, gubconsGFC1, gubconsGR, ngubconsGC1, ngubconsGC2, ngubconsGFC1, ngubconsGR,
            MIN(nvarsC1 - 1, ngubconsGC1), liftcoefs, &cutact, &liftrhs, maxgubvarssize) );

      /* frees temporary memory */
      SCIPfreeBufferArray(scip, &gubconsGR);
      SCIPfreeBufferArray(scip, &gubconsGFC1);
      SCIPfreeBufferArray(scip, &gubconsGC2);
      SCIPfreeBufferArray(scip, &gubconsGC1);
   }

   /* checks, if lifting yielded a violated cut */
   if( SCIPisEfficacious(scip, (cutact - liftrhs)/sqrt((SCIP_Real)MAX(liftrhs, 1))) )
   {
      SCIP_ROW* row;
      char name[SCIP_MAXSTRLEN];
      int j;

      /* creates LP row */
      assert( cons == NULL || sepa == NULL );
      if ( cons != NULL )
      {
         (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "%s_mcseq%" SCIP_LONGINT_FORMAT "", SCIPconsGetName(cons), SCIPconshdlrGetNCutsFound(SCIPconsGetHdlr(cons)));
         SCIP_CALL( SCIPcreateEmptyRowCons(scip, &row, cons, name, -SCIPinfinity(scip), (SCIP_Real)liftrhs,
            cons != NULL ? SCIPconsIsLocal(cons) : FALSE, FALSE,
            cons != NULL ? SCIPconsIsRemovable(cons) : TRUE) );
      }
      else if ( sepa != NULL )
      {
         (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "%s_mcseq_%" SCIP_LONGINT_FORMAT "", SCIPsepaGetName(sepa), SCIPsepaGetNCutsFound(sepa));
         SCIP_CALL( SCIPcreateEmptyRowSepa(scip, &row, sepa, name, -SCIPinfinity(scip), (SCIP_Real)liftrhs, FALSE, FALSE, TRUE) );
      }
      else
      {
         (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "nn_mcseq_%d", *ncuts);
         SCIP_CALL( SCIPcreateEmptyRowUnspec(scip, &row, name, -SCIPinfinity(scip), (SCIP_Real)liftrhs, FALSE, FALSE, TRUE) );
      }

      /* adds all variables in the knapsack constraint with calculated lifting coefficient to the cut */
      SCIP_CALL( SCIPcacheRowExtensions(scip, row) );
      assert(nvarsC1 + nvarsC2 + nvarsF + nvarsR == nvars - ntightened);
      for( j = 0; j < nvarsC1; j++ )
      {
         SCIP_CALL( SCIPaddVarToRow(scip, row, vars[varsC1[j]], 1.0) );
      }
      for( j = 0; j < nvarsC2; j++ )
      {
         if( liftcoefs[varsC2[j]] > 0 )
         {
            SCIP_CALL( SCIPaddVarToRow(scip, row, vars[varsC2[j]], (SCIP_Real)liftcoefs[varsC2[j]]) );
         }
      }
      for( j = 0; j < nvarsF; j++ )
      {
         if( liftcoefs[varsF[j]] > 0 )
         {
            SCIP_CALL( SCIPaddVarToRow(scip, row, vars[varsF[j]], (SCIP_Real)liftcoefs[varsF[j]]) );
         }
      }
      for( j = 0; j < nvarsR; j++ )
      {
         if( liftcoefs[varsR[j]] > 0 )
         {
            SCIP_CALL( SCIPaddVarToRow(scip, row, vars[varsR[j]], (SCIP_Real)liftcoefs[varsR[j]]) );
         }
      }
      SCIP_CALL( SCIPflushRowExtensions(scip, row) );

      /* checks, if cut is violated enough */
      if( SCIPisCutEfficacious(scip, sol, row) )
      {
         if( cons != NULL )
         {
            SCIP_CALL( SCIPresetConsAge(scip, cons) );
         }
         SCIP_CALL( SCIPaddRow(scip, row, FALSE, cutoff) );
         (*ncuts)++;
      }
      SCIP_CALL( SCIPreleaseRow(scip, &row) );
   }

   /* frees temporary memory */
   SCIPfreeBufferArray(scip, &liftcoefs);
   SCIPfreeBufferArray(scip, &varsR);
   SCIPfreeBufferArray(scip, &varsF);
   SCIPfreeBufferArray(scip, &varsC2);
   SCIPfreeBufferArray(scip, &varsC1);

   return SCIP_OKAY;
}

/** separates lifted extended weight inequalities using sequential up- and down-lifting for given knapsack problem */
static
SCIP_RETCODE separateSequLiftedExtendedWeightInequality(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONS*            cons,               /**< constraint that originates the knapsack problem, or NULL */
   SCIP_SEPA*            sepa,               /**< originating separator of the knapsack problem, or NULL */
   SCIP_VAR**            vars,               /**< variables in knapsack constraint */
   int                   nvars,              /**< number of variables in knapsack constraint */
   int                   ntightened,         /**< number of variables with tightened upper bound */
   SCIP_Longint*         weights,            /**< weights of variables in knapsack constraint */
   SCIP_Longint          capacity,           /**< capacity of knapsack */
   SCIP_Real*            solvals,            /**< solution values of all problem variables */
   int*                  feassetvars,        /**< variables in feasible set */
   int*                  nonfeassetvars,     /**< variables not in feasible set */
   int                   nfeassetvars,       /**< number of variables in feasible set */
   int                   nnonfeassetvars,    /**< number of variables not in feasible set */
   SCIP_SOL*             sol,                /**< primal SCIP solution to separate, NULL for current LP solution */
   SCIP_Bool*            cutoff,             /**< whether a cutoff has been detected */
   int*                  ncuts               /**< pointer to add up the number of found cuts */
   )
{
   int* varsT1;
   int* varsT2;
   int* varsF;
   int* varsR;
   int* liftcoefs;
   SCIP_Real cutact;
   int nvarsT1;
   int nvarsT2;
   int nvarsF;
   int nvarsR;
   int liftrhs;
   int j;

   assert( cutoff != NULL );
   *cutoff = FALSE;

   /* allocates temporary memory */
   SCIP_CALL( SCIPallocBufferArray(scip, &varsT1, nvars) );
   SCIP_CALL( SCIPallocBufferArray(scip, &varsT2, nvars) );
   SCIP_CALL( SCIPallocBufferArray(scip, &varsF, nvars) );
   SCIP_CALL( SCIPallocBufferArray(scip, &varsR, nvars) );
   SCIP_CALL( SCIPallocBufferArray(scip, &liftcoefs, nvars) );

   /* gets partition (T_1,T_2) of T, i.e. T_1 & T_2 = T and T_1 cap T_2 = emptyset, with T_1 not empty; chooses partition
    * as follows
    *   T_2 = { j in T : x*_j = 1 } and
    *   T_1 = T\T_2
    */
   getPartitionCovervars(scip, solvals, feassetvars, nfeassetvars, varsT1, varsT2, &nvarsT1, &nvarsT2);
   assert(nvarsT1 + nvarsT2 == nfeassetvars);

   /* changes partition (T_1,T_2) of feasible set T, if |T1| = 0, by moving one variable from T2 to T1 */
   if( nvarsT1 == 0 && nvarsT2 > 0)
   {
      SCIP_CALL( changePartitionFeasiblesetvars(scip, weights, varsT1, varsT2, &nvarsT1, &nvarsT2) );
      assert(nvarsT1 == 1);
   }
   assert(nvarsT2 == 0 || nvarsT1 > 0);

   /* gets partition (F,R) of N\T, i.e. F & R = N\T and F cap R = emptyset; chooses partition as follows
    *   R = { j in N\T : x*_j = 0 } and
    *   F = (N\T)\F
    */
   getPartitionNoncovervars(scip, solvals, nonfeassetvars, nnonfeassetvars, varsF, varsR, &nvarsF, &nvarsR);
   assert(nvarsF + nvarsR == nnonfeassetvars);
   assert(nvarsT1 + nvarsT2 + nvarsF + nvarsR == nvars - ntightened);

   /* sorts variables in F, T_2, and R according to the second level lifting sequence that will be used in the sequential
    * lifting procedure (the variable removed last from the initial cover does not have to be lifted first, therefore it
    * is included in the sorting routine)
    */
   SCIP_CALL( getLiftingSequence(scip, solvals, weights, varsF, varsT2, varsR, nvarsF, nvarsT2, nvarsR) );

   /* lifts extended weight inequality sum_{j in T_1} x_j <= |T_1| valid for
    *
    *    S^0 = { x in {0,1}^|T_1| : sum_{j in T_1} a_j x_j <= a_0 - sum_{j in T_2} a_j }
    *
    * to a valid inequality sum_{j in T_1} x_j + sum_{j in N\T_1} alpha_j x_j <= |T_1| + sum_{j in T_2} alpha_j for
    *
    *      S = { x in {0,1}^|N|   : sum_{j in N}   a_j x_j <= a_0 },
    *
    * uses sequential up-lifting for the variables in F, sequential down-lifting for the variable in T_2 and sequential
    * up-lifting for the variabels in R according to the second level lifting sequence
    */
   SCIP_CALL( sequentialUpAndDownLifting(scip, vars, nvars, ntightened, weights, capacity, solvals, varsT1, varsT2, varsF, varsR,
         nvarsT1, nvarsT2, nvarsF, nvarsR, nvarsT1, liftcoefs, &cutact, &liftrhs) );

   /* checks, if lifting yielded a violated cut */
   if( SCIPisEfficacious(scip, (cutact - liftrhs)/sqrt((SCIP_Real)MAX(liftrhs, 1))) )
   {
      SCIP_ROW* row;
      char name[SCIP_MAXSTRLEN];

      /* creates LP row */
      assert( cons == NULL || sepa == NULL );
      if( cons != NULL )
      {
         (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "%s_ewseq%" SCIP_LONGINT_FORMAT "", SCIPconsGetName(cons), SCIPconshdlrGetNCutsFound(SCIPconsGetHdlr(cons)));
         SCIP_CALL( SCIPcreateEmptyRowConshdlr(scip, &row, SCIPconsGetHdlr(cons), name, -SCIPinfinity(scip), (SCIP_Real)liftrhs,
               cons != NULL ? SCIPconsIsLocal(cons) : FALSE, FALSE,
               cons != NULL ? SCIPconsIsRemovable(cons) : TRUE) );
      }
      else if ( sepa != NULL )
      {
         (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "%s_ewseq_%" SCIP_LONGINT_FORMAT "", SCIPsepaGetName(sepa), SCIPsepaGetNCutsFound(sepa));
         SCIP_CALL( SCIPcreateEmptyRowSepa(scip, &row, sepa, name, -SCIPinfinity(scip), (SCIP_Real)liftrhs, FALSE, FALSE, TRUE) );
      }
      else
      {
         (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "nn_ewseq_%d", *ncuts);
         SCIP_CALL( SCIPcreateEmptyRowUnspec(scip, &row, name, -SCIPinfinity(scip), (SCIP_Real)liftrhs, FALSE, FALSE, TRUE) );
      }

      /* adds all variables in the knapsack constraint with calculated lifting coefficient to the cut */
      SCIP_CALL( SCIPcacheRowExtensions(scip, row) );
      assert(nvarsT1 + nvarsT2 + nvarsF + nvarsR == nvars - ntightened);
      for( j = 0; j < nvarsT1; j++ )
      {
         SCIP_CALL( SCIPaddVarToRow(scip, row, vars[varsT1[j]], 1.0) );
      }
      for( j = 0; j < nvarsT2; j++ )
      {
         if( liftcoefs[varsT2[j]] > 0 )
         {
            SCIP_CALL( SCIPaddVarToRow(scip, row, vars[varsT2[j]], (SCIP_Real)liftcoefs[varsT2[j]]) );
         }
      }
      for( j = 0; j < nvarsF; j++ )
      {
         if( liftcoefs[varsF[j]] > 0 )
         {
            SCIP_CALL( SCIPaddVarToRow(scip, row, vars[varsF[j]], (SCIP_Real)liftcoefs[varsF[j]]) );
         }
      }
      for( j = 0; j < nvarsR; j++ )
      {
         if( liftcoefs[varsR[j]] > 0 )
         {
            SCIP_CALL( SCIPaddVarToRow(scip, row, vars[varsR[j]], (SCIP_Real)liftcoefs[varsR[j]]) );
         }
      }
      SCIP_CALL( SCIPflushRowExtensions(scip, row) );

      /* checks, if cut is violated enough */
      if( SCIPisCutEfficacious(scip, sol, row) )
      {
         if( cons != NULL )
         {
            SCIP_CALL( SCIPresetConsAge(scip, cons) );
         }
         SCIP_CALL( SCIPaddRow(scip, row, FALSE, cutoff) );
         (*ncuts)++;
      }
      SCIP_CALL( SCIPreleaseRow(scip, &row) );
   }

   /* frees temporary memory */
   SCIPfreeBufferArray(scip, &liftcoefs);
   SCIPfreeBufferArray(scip, &varsR);
   SCIPfreeBufferArray(scip, &varsF);
   SCIPfreeBufferArray(scip, &varsT2);
   SCIPfreeBufferArray(scip, &varsT1);

   return SCIP_OKAY;
}

/** separates lifted minimal cover inequalities using superadditive up-lifting for given knapsack problem */
static
SCIP_RETCODE separateSupLiftedMinimalCoverInequality(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONS*            cons,               /**< constraint that originates the knapsack problem, or NULL */
   SCIP_SEPA*            sepa,               /**< originating separator of the knapsack problem, or NULL */
   SCIP_VAR**            vars,               /**< variables in knapsack constraint */
   int                   nvars,              /**< number of variables in knapsack constraint */
   int                   ntightened,         /**< number of variables with tightened upper bound */
   SCIP_Longint*         weights,            /**< weights of variables in knapsack constraint */
   SCIP_Longint          capacity,           /**< capacity of knapsack */
   SCIP_Real*            solvals,            /**< solution values of all problem variables */
   int*                  mincovervars,       /**< mincover variables */
   int*                  nonmincovervars,    /**< nonmincover variables */
   int                   nmincovervars,      /**< number of mincover variables */
   int                   nnonmincovervars,   /**< number of nonmincover variables */
   SCIP_Longint          mincoverweight,     /**< weight of minimal cover */
   SCIP_SOL*             sol,                /**< primal SCIP solution to separate, NULL for current LP solution */
   SCIP_Bool*            cutoff,             /**< whether a cutoff has been detected */
   int*                  ncuts               /**< pointer to add up the number of found cuts */
   )
{
   SCIP_Real* realliftcoefs;
   SCIP_Real cutact;
   int liftrhs;

   assert( cutoff != NULL );
   *cutoff = FALSE;
   cutact = 0.0;

   /* allocates temporary memory */
   SCIP_CALL( SCIPallocBufferArray(scip, &realliftcoefs, nvars) );

   /* lifts minimal cover inequality sum_{j in C} x_j <= |C| - 1 valid for
    *
    *    S^0 = { x in {0,1}^|C| : sum_{j in C} a_j x_j <= a_0 }
    *
    * to a valid inequality sum_{j in C} x_j + sum_{j in N\C} alpha_j x_j <= |C| - 1 for
    *
    *      S = { x in {0,1}^|N|   : sum_{j in N}   a_j x_j <= a_0 },
    *
    * uses superadditive up-lifting for the variables in N\C.
    */
   SCIP_CALL( superadditiveUpLifting(scip, vars, nvars, ntightened, weights, capacity, solvals, mincovervars,
         nonmincovervars, nmincovervars, nnonmincovervars, mincoverweight, realliftcoefs, &cutact) );
   liftrhs = nmincovervars - 1;

   /* checks, if lifting yielded a violated cut */
   if( SCIPisEfficacious(scip, (cutact - liftrhs)/sqrt((SCIP_Real)MAX(liftrhs, 1))) )
   {
      SCIP_ROW* row;
      char name[SCIP_MAXSTRLEN];
      int j;

      /* creates LP row */
      assert( cons == NULL || sepa == NULL );
      if ( cons != NULL )
      {
         (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "%s_mcsup%" SCIP_LONGINT_FORMAT "", SCIPconsGetName(cons), SCIPconshdlrGetNCutsFound(SCIPconsGetHdlr(cons)));
         SCIP_CALL( SCIPcreateEmptyRowConshdlr(scip, &row, SCIPconsGetHdlr(cons), name, -SCIPinfinity(scip), (SCIP_Real)liftrhs,
               cons != NULL ? SCIPconsIsLocal(cons) : FALSE, FALSE,
               cons != NULL ? SCIPconsIsRemovable(cons) : TRUE) );
      }
      else if ( sepa != NULL )
      {
         (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "%s_mcsup%" SCIP_LONGINT_FORMAT "", SCIPsepaGetName(sepa), SCIPsepaGetNCutsFound(sepa));
         SCIP_CALL( SCIPcreateEmptyRowSepa(scip, &row, sepa, name, -SCIPinfinity(scip), (SCIP_Real)liftrhs, FALSE, FALSE, TRUE) );
      }
      else
      {
         (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "nn_mcsup_%d", *ncuts);
         SCIP_CALL( SCIPcreateEmptyRowUnspec(scip, &row, name, -SCIPinfinity(scip), (SCIP_Real)liftrhs, FALSE, FALSE, TRUE) );
      }

      /* adds all variables in the knapsack constraint with calculated lifting coefficient to the cut */
      SCIP_CALL( SCIPcacheRowExtensions(scip, row) );
      assert(nmincovervars + nnonmincovervars == nvars - ntightened);
      for( j = 0; j < nmincovervars; j++ )
      {
         SCIP_CALL( SCIPaddVarToRow(scip, row, vars[mincovervars[j]], 1.0) );
      }
      for( j = 0; j < nnonmincovervars; j++ )
      {
         assert(SCIPisFeasGE(scip, realliftcoefs[nonmincovervars[j]], 0.0));
         if( SCIPisFeasGT(scip, realliftcoefs[nonmincovervars[j]], 0.0) )
         {
            SCIP_CALL( SCIPaddVarToRow(scip, row, vars[nonmincovervars[j]], realliftcoefs[nonmincovervars[j]]) );
         }
      }
      SCIP_CALL( SCIPflushRowExtensions(scip, row) );

      /* checks, if cut is violated enough */
      if( SCIPisCutEfficacious(scip, sol, row) )
      {
         if( cons != NULL )
         {
            SCIP_CALL( SCIPresetConsAge(scip, cons) );
         }
         SCIP_CALL( SCIPaddRow(scip, row, FALSE, cutoff) );
         (*ncuts)++;
      }
      SCIP_CALL( SCIPreleaseRow(scip, &row) );
   }

   /* frees temporary memory */
   SCIPfreeBufferArray(scip, &realliftcoefs);

   return SCIP_OKAY;
}

/** converts given cover C to a minimal cover by removing variables in the reverse order in which the variables were chosen
 *  to be in C, i.e. in the order of non-increasing (1 - x*_j)/a_j, if the transformed separation problem was used to find
 *  C and in the order of non-increasing (1 - x*_j), if the modified transformed separation problem was used to find C; 
 *  note that all variables with x*_j = 1 will be removed last
 */
static
SCIP_RETCODE makeCoverMinimal(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_Longint*         weights,            /**< weights of variables in knapsack constraint */
   SCIP_Longint          capacity,           /**< capacity of knapsack */
   SCIP_Real*            solvals,            /**< solution values of all problem variables */
   int*                  covervars,          /**< pointer to store cover variables */
   int*                  noncovervars,       /**< pointer to store noncover variables */
   int*                  ncovervars,         /**< pointer to store number of cover variables */
   int*                  nnoncovervars,      /**< pointer to store number of noncover variables */
   SCIP_Longint*         coverweight,        /**< pointer to store weight of cover */
   SCIP_Bool             modtransused        /**< TRUE if mod trans sepa prob was used to find cover */    
   )
{
   SORTKEYPAIR** sortkeypairs;
   SORTKEYPAIR** sortkeypairssorted;
   SCIP_Longint minweight;
   int nsortkeypairs;
   int minweightidx;
   int j;
   int k;

   assert(scip != NULL);
   assert(covervars != NULL);
   assert(noncovervars != NULL);
   assert(ncovervars != NULL);
   assert(*ncovervars > 0);
   assert(nnoncovervars != NULL);
   assert(*nnoncovervars >= 0);
   assert(coverweight != NULL);
   assert(*coverweight > 0);
   assert(*coverweight > capacity);

   /* allocates temporary memory; we need two arrays for the keypairs in order to be able to free them in the correct
    * order */
   nsortkeypairs = *ncovervars;
   SCIP_CALL( SCIPallocBufferArray(scip, &sortkeypairs, nsortkeypairs) );
   SCIP_CALL( SCIPallocBufferArray(scip, &sortkeypairssorted, nsortkeypairs) );

   /* sorts C in the reverse order in which the variables were chosen to be in the cover, i.e. 
    *   such that (1 - x*_1)/a_1 >= ... >= (1 - x*_|C|)/a_|C|,  if          trans separation problem was used to find C 
    *   such that (1 - x*_1)     >= ... >= (1 - x*_|C|),        if modified trans separation problem was used to find C 
    * note that all variables with x*_j = 1 are in the end of the sorted C, so they will be removed last from C
    */
   assert(*ncovervars == nsortkeypairs);
   if( modtransused )
   {
      for( j = 0; j < *ncovervars; j++ )
      {
         SCIP_CALL( SCIPallocBuffer(scip, &(sortkeypairs[j])) ); /*lint !e866 */
         sortkeypairssorted[j] = sortkeypairs[j];

         sortkeypairs[j]->key1 = solvals[covervars[j]]; 
         sortkeypairs[j]->key2 = (SCIP_Real) weights[covervars[j]]; 
      }
   }
   else
   {
      for( j = 0; j < *ncovervars; j++ )
      {
         SCIP_CALL( SCIPallocBuffer(scip, &(sortkeypairs[j])) ); /*lint !e866 */
         sortkeypairssorted[j] = sortkeypairs[j];

         sortkeypairs[j]->key1 = (solvals[covervars[j]] - 1.0) / ((SCIP_Real) weights[covervars[j]]);
         sortkeypairs[j]->key2 = (SCIP_Real) (-weights[covervars[j]]);
      }
   }
   SCIPsortPtrInt((void**)sortkeypairssorted, covervars, compSortkeypairs, *ncovervars);

   /* gets j' with a_j' = min{ a_j : j in C } */
   minweightidx = 0;
   minweight = weights[covervars[minweightidx]];
   for( j = 1; j < *ncovervars; j++ )
   {
      if( weights[covervars[j]] <= minweight )
      {
         minweightidx = j;
         minweight = weights[covervars[minweightidx]];
      }
   }
   assert(minweightidx >= 0 && minweightidx < *ncovervars);
   assert(minweight > 0 && minweight <= *coverweight);

   j = 0;
   /* removes variables from C until the remaining variables form a minimal cover */
   while( j < *ncovervars && ((*coverweight) - minweight > capacity) )
   {
      assert(minweightidx >= j);
      assert(checkMinweightidx(weights, capacity, covervars, *ncovervars, *coverweight, minweightidx, j));

      /* if sum_{i in C} a_i - a_j <= a_0, j cannot be removed from C */
      if( (*coverweight) - weights[covervars[j]] <= capacity )
      {
         ++j;
         continue;
      }

      /* adds j to N\C */
      noncovervars[*nnoncovervars] = covervars[j];
      (*nnoncovervars)++;

      /* removes j from C */
      (*coverweight) -= weights[covervars[j]];
      for( k = j; k < (*ncovervars) - 1; k++ )
         covervars[k] = covervars[k+1];
      (*ncovervars)--;

      /* updates j' with a_j' = min{ a_j : j in C } */
      if( j == minweightidx )
      {
         minweightidx = 0;
         minweight = weights[covervars[minweightidx]];
         for( k = 1; k < *ncovervars; k++ )
         {
            if( weights[covervars[k]] <= minweight )
            {
               minweightidx = k;
               minweight = weights[covervars[minweightidx]];
            }
         }
         assert(minweight > 0 && minweight <= *coverweight);
         assert(minweightidx >= 0 && minweightidx < *ncovervars);
      }
      else
      {
         assert(minweightidx > j);
         minweightidx--;
      }
      /* j needs to stay the same */
   }
   assert((*coverweight) > capacity);
   assert((*coverweight) - minweight <= capacity);

   /* frees temporary memory */
   for( j = nsortkeypairs-1; j >= 0; j-- )
      SCIPfreeBuffer(scip, &(sortkeypairs[j])); /*lint !e866 */
   SCIPfreeBufferArray(scip, &sortkeypairssorted);
   SCIPfreeBufferArray(scip, &sortkeypairs);

   return SCIP_OKAY;
}

/** converts given initial cover C_init to a feasible set by removing variables in the reverse order in which
 *  they were chosen to be in C_init:
 *   non-increasing (1 - x*_j)/a_j,   if          transformed separation problem was used to find C_init
 *   non-increasing (1 - x*_j),       if modified transformed separation problem was used to find C_init.
 *  separates lifted extended weight inequalities using sequential up- and down-lifting for this feasible set
 *  and all subsequent feasible sets.
 */
static
SCIP_RETCODE getFeasibleSet(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONS*            cons,               /**< constraint that originates the knapsack problem */
   SCIP_SEPA*            sepa,               /**< originating separator of the knapsack problem, or NULL */
   SCIP_VAR**            vars,               /**< variables in knapsack constraint */
   int                   nvars,              /**< number of variables in knapsack constraint */
   int                   ntightened,         /**< number of variables with tightened upper bound */
   SCIP_Longint*         weights,            /**< weights of variables in knapsack constraint */
   SCIP_Longint          capacity,           /**< capacity of knapsack */
   SCIP_Real*            solvals,            /**< solution values of all problem variables */
   int*                  covervars,          /**< pointer to store cover variables */
   int*                  noncovervars,       /**< pointer to store noncover variables */
   int*                  ncovervars,         /**< pointer to store number of cover variables */
   int*                  nnoncovervars,      /**< pointer to store number of noncover variables */
   SCIP_Longint*         coverweight,        /**< pointer to store weight of cover */
   SCIP_Bool             modtransused,       /**< TRUE if mod trans sepa prob was used to find cover */
   SCIP_SOL*             sol,                /**< primal SCIP solution to separate, NULL for current LP solution */
   SCIP_Bool*            cutoff,             /**< whether a cutoff has been detected */
   int*                  ncuts               /**< pointer to add up the number of found cuts */
   )
{
   SCIP_Real* sortkeys;
   int j;
   int k;

   assert(scip != NULL);
   assert(covervars != NULL);
   assert(noncovervars != NULL);
   assert(ncovervars != NULL);
   assert(*ncovervars > 0);
   assert(nnoncovervars != NULL);
   assert(*nnoncovervars >= 0);
   assert(coverweight != NULL);
   assert(*coverweight > 0);
   assert(*coverweight > capacity);
   assert(*ncovervars + *nnoncovervars == nvars - ntightened);
   assert(cutoff != NULL);

   *cutoff = FALSE;

   /* allocates temporary memory */
   SCIP_CALL( SCIPallocBufferArray(scip, &sortkeys, *ncovervars) );

   /* sorts C in the reverse order in which the variables were chosen to be in the cover, i.e.
    *   such that (1 - x*_1)/a_1 >= ... >= (1 - x*_|C|)/a_|C|,  if          trans separation problem was used to find C
    *   such that (1 - x*_1)     >= ... >= (1 - x*_|C|),        if modified trans separation problem was used to find C
    * note that all variables with x*_j = 1 are in the end of the sorted C, so they will be removed last from C
    */
   if( modtransused )
   {
      for( j = 0; j < *ncovervars; j++ )
      {
         sortkeys[j] = solvals[covervars[j]];
         assert(SCIPisFeasGE(scip, sortkeys[j], 0.0));
      }
   }
   else
   {
      for( j = 0; j < *ncovervars; j++ )
      {
         sortkeys[j] = (solvals[covervars[j]] - 1.0) / ((SCIP_Real) weights[covervars[j]]);
         assert(SCIPisFeasLE(scip, sortkeys[j], 0.0));
      }
   }
   SCIPsortRealInt(sortkeys, covervars, *ncovervars);

   /* removes variables from C_init and separates lifted extended weight inequalities using sequential up- and down-lifting;
    * in addition to an extended weight inequality this gives cardinality inequalities */
   while( *ncovervars >= 2 )
   {
      /* adds first element of C_init to N\C_init */
      noncovervars[*nnoncovervars] = covervars[0];
      (*nnoncovervars)++;

      /* removes first element from C_init */
      (*coverweight) -= weights[covervars[0]];
      for( k = 0; k < (*ncovervars) - 1; k++ )
         covervars[k] = covervars[k+1];
      (*ncovervars)--;

      assert(*ncovervars + *nnoncovervars == nvars - ntightened);
      if( (*coverweight) <= capacity )
      {
         SCIP_CALL( separateSequLiftedExtendedWeightInequality(scip, cons, sepa, vars, nvars, ntightened, weights, capacity, solvals,
               covervars, noncovervars, *ncovervars, *nnoncovervars, sol, cutoff, ncuts) );
      }

      /* stop if cover is too large */
      if ( *ncovervars >= MAXCOVERSIZEITERLEWI )
         break;
   }

   /* frees temporary memory */
   SCIPfreeBufferArray(scip, &sortkeys);

   return SCIP_OKAY;
}

/** separates different classes of valid inequalities for the 0-1 knapsack problem */
SCIP_RETCODE SCIPseparateKnapsackCuts(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONS*            cons,               /**< originating constraint of the knapsack problem, or NULL */
   SCIP_SEPA*            sepa,               /**< originating separator of the knapsack problem, or NULL */
   SCIP_VAR**            vars,               /**< variables in knapsack constraint */
   int                   nvars,              /**< number of variables in knapsack constraint */
   SCIP_Longint*         weights,            /**< weights of variables in knapsack constraint */
   SCIP_Longint          capacity,           /**< capacity of knapsack */
   SCIP_SOL*             sol,                /**< primal SCIP solution to separate, NULL for current LP solution */
   SCIP_Bool             usegubs,            /**< should GUB information be used for separation? */
   SCIP_Bool*            cutoff,             /**< pointer to store whether a cutoff has been detected */
   int*                  ncuts               /**< pointer to add up the number of found cuts */
   )
{
   SCIP_Real* solvals;
   int* covervars;
   int* noncovervars;
   SCIP_Bool coverfound;
   SCIP_Bool fractional;
   SCIP_Bool modtransused;
   SCIP_Longint coverweight;
   int ncovervars;
   int nnoncovervars;
   int ntightened;

   assert(scip != NULL);
   assert(capacity >= 0);
   assert(cutoff != NULL);
   assert(ncuts != NULL);

   *cutoff = FALSE;

   if( nvars == 0 )
      return SCIP_OKAY;

   assert(vars != NULL);
   assert(nvars > 0);
   assert(weights != NULL);

   /* increase age of constraint (age is reset to zero, if a cut was found) */
   if( cons != NULL )
   {
      SCIP_CALL( SCIPincConsAge(scip, cons) );
   }

   /* allocates temporary memory */
   SCIP_CALL( SCIPallocBufferArray(scip, &solvals, nvars) );
   SCIP_CALL( SCIPallocBufferArray(scip, &covervars, nvars) );
   SCIP_CALL( SCIPallocBufferArray(scip, &noncovervars, nvars) );

   /* gets solution values of all problem variables */
   SCIP_CALL( SCIPgetSolVals(scip, sol, nvars, vars, solvals) );

#ifdef SCIP_DEBUG
   {
      int i;

      SCIPdebugMsg(scip, "separate cuts for knapsack constraint originated by cons <%s>:\n",
         cons == NULL ? "-" : SCIPconsGetName(cons));
      for( i = 0; i < nvars; ++i )
      {
         SCIPdebugMsgPrint(scip, "%+" SCIP_LONGINT_FORMAT "<%s>(%g)", weights[i], SCIPvarGetName(vars[i]), solvals[i]);
      }
      SCIPdebugMsgPrint(scip, " <= %" SCIP_LONGINT_FORMAT "\n", capacity);
   }
#endif

   /* LMCI1 (lifted minimal cover inequalities using sequential up- and down-lifting) using GUB information
    */
   if( usegubs )
   {
      SCIP_GUBSET* gubset;

      SCIPdebugMsg(scip, "separate LMCI1-GUB cuts:\n");

      /* initializes partion of knapsack variables into nonoverlapping GUB constraints */
      SCIP_CALL( GUBsetCreate(scip, &gubset, nvars, weights, capacity) );

      /* constructs sophisticated partition of knapsack variables into nonoverlapping GUBs */
      SCIP_CALL( GUBsetGetCliquePartition(scip, gubset, vars, solvals) );
      assert(gubset->ngubconss <= nvars);

      /* gets a most violated initial cover C_init ( sum_{j in C_init} a_j > a_0 ) by using the
       * MODIFIED transformed separation problem and taking into account the following fixing:
       *   j in C_init,           if j in N_1 = {j in N : x*_j = 1} and
       *   j in N\C_init,         if j in N_0 = {j in N : x*_j = 0},
       * if one exists
       */
      modtransused = TRUE;
      SCIP_CALL( getCover(scip, vars, nvars, weights, capacity, solvals, covervars, noncovervars, &ncovervars,
            &nnoncovervars, &coverweight, &coverfound, modtransused, &ntightened, &fractional) );

      assert(!coverfound || !fractional || ncovervars + nnoncovervars == nvars - ntightened);

      /* if x* is not fractional we stop the separation routine */
      if( !fractional )
      {
         SCIPdebugMsg(scip, "   LMCI1-GUB terminated by no variable with fractional LP value.\n");

         /* frees memory for GUB set data structure */
         GUBsetFree(scip, &gubset);

         goto TERMINATE;
      }

      /* if no cover was found we stop the separation routine for lifted minimal cover inequality */
      if( coverfound )
      {
         /* converts initial cover C_init to a minimal cover C by removing variables in the reverse order in which the
          * variables were chosen to be in C_init; note that variables with x*_j = 1 will be removed last
          */
         SCIP_CALL( makeCoverMinimal(scip, weights, capacity, solvals, covervars, noncovervars, &ncovervars,
               &nnoncovervars, &coverweight, modtransused) );

         /* only separate with GUB information if we have at least one nontrivial GUB (with more than one variable) */
         if( gubset->ngubconss < nvars )
         {
            /* separates lifted minimal cover inequalities using sequential up- and down-lifting and GUB information */
            SCIP_CALL( separateSequLiftedMinimalCoverInequality(scip, cons, sepa, vars, nvars, ntightened, weights, capacity,
                  solvals, covervars, noncovervars, ncovervars, nnoncovervars, sol, gubset, cutoff, ncuts) );
         }
         else
         {
            /* separates lifted minimal cover inequalities using sequential up- and down-lifting, but do not use trivial
             * GUB information
             */
            SCIP_CALL( separateSequLiftedMinimalCoverInequality(scip, cons, sepa, vars, nvars, ntightened, weights, capacity,
                  solvals, covervars, noncovervars, ncovervars, nnoncovervars, sol, NULL, cutoff, ncuts) );
         }
      }

      /* frees memory for GUB set data structure */
      GUBsetFree(scip, &gubset);
   }
   else
   {
      /* LMCI1 (lifted minimal cover inequalities using sequential up- and down-lifting)
       * (and LMCI2 (lifted minimal cover inequalities using superadditive up-lifting))
       */

      /* gets a most violated initial cover C_init ( sum_{j in C_init} a_j > a_0 ) by using the
       * MODIFIED transformed separation problem and taking into account the following fixing:
       *   j in C_init,           if j in N_1 = {j in N : x*_j = 1} and
       *   j in N\C_init,         if j in N_0 = {j in N : x*_j = 0},
       * if one exists
       */
      SCIPdebugMsg(scip, "separate LMCI1 cuts:\n");
      modtransused = TRUE;
      SCIP_CALL( getCover(scip, vars, nvars, weights, capacity, solvals, covervars, noncovervars, &ncovervars,
            &nnoncovervars, &coverweight, &coverfound, modtransused, &ntightened, &fractional) );
      assert(!coverfound || !fractional || ncovervars + nnoncovervars == nvars - ntightened);

      /* if x* is not fractional we stop the separation routine */
      if( !fractional )
         goto TERMINATE;

      /* if no cover was found we stop the separation routine for lifted minimal cover inequality */
      if( coverfound )
      {
         /* converts initial cover C_init to a minimal cover C by removing variables in the reverse order in which the
          * variables were chosen to be in C_init; note that variables with x*_j = 1 will be removed last
          */
         SCIP_CALL( makeCoverMinimal(scip, weights, capacity, solvals, covervars, noncovervars, &ncovervars,
               &nnoncovervars, &coverweight, modtransused) );

         /* separates lifted minimal cover inequalities using sequential up- and down-lifting */
         SCIP_CALL( separateSequLiftedMinimalCoverInequality(scip, cons, sepa, vars, nvars, ntightened, weights, capacity,
               solvals, covervars, noncovervars, ncovervars, nnoncovervars, sol, NULL, cutoff, ncuts) );

         if( USESUPADDLIFT ) /*lint !e506 !e774*/
         {
            SCIPdebugMsg(scip, "separate LMCI2 cuts:\n");
            /* separates lifted minimal cover inequalities using superadditive up-lifting */
            SCIP_CALL( separateSupLiftedMinimalCoverInequality(scip, cons, sepa, vars, nvars, ntightened, weights, capacity,
                  solvals, covervars, noncovervars, ncovervars, nnoncovervars, coverweight, sol, cutoff, ncuts) );
         }
      }
   }

   /* LEWI (lifted extended weight inequalities using sequential up- and down-lifting) */
   if ( ! (*cutoff) )
   {
      /* gets a most violated initial cover C_init ( sum_{j in C_init} a_j > a_0 ) by using the
       * transformed separation problem and taking into account the following fixing:
       *   j in C_init,           if j in N_1 = {j in N : x*_j = 1} and
       *   j in N\C_init,         if j in N_0 = {j in N : x*_j = 0},
       * if one exists
       */
      SCIPdebugMsg(scip, "separate LEWI cuts:\n");
      modtransused = FALSE;
      SCIP_CALL( getCover(scip, vars, nvars, weights, capacity, solvals, covervars, noncovervars, &ncovervars,
            &nnoncovervars, &coverweight, &coverfound, modtransused, &ntightened, &fractional) );
      assert(fractional);
      assert(!coverfound || ncovervars + nnoncovervars == nvars - ntightened);

      /* if no cover was found we stop the separation routine */
      if( coverfound )
      {
         /* converts initial cover C_init to a feasible set by removing variables in the reverse order in which
          * they were chosen to be in C_init and separates lifted extended weight inequalities using sequential
          * up- and down-lifting for this feasible set and all subsequent feasible sets.
          */
         SCIP_CALL( getFeasibleSet(scip, cons, sepa, vars, nvars, ntightened, weights, capacity, solvals, covervars, noncovervars,
               &ncovervars, &nnoncovervars, &coverweight, modtransused, sol, cutoff, ncuts) );
      }
   }

 TERMINATE:
   /* frees temporary memory */
   SCIPfreeBufferArray(scip, &noncovervars);
   SCIPfreeBufferArray(scip, &covervars);
   SCIPfreeBufferArray(scip, &solvals);

   return SCIP_OKAY;
}

/* relaxes given general linear constraint into a knapsack constraint and separates lifted knapsack cover inequalities */
SCIP_RETCODE SCIPseparateRelaxedKnapsack(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONS*            cons,               /**< originating constraint of the knapsack problem, or NULL */
   SCIP_SEPA*            sepa,               /**< originating separator of the knapsack problem, or NULL */
   int                   nknapvars,          /**< number of variables in the continuous knapsack constraint */
   SCIP_VAR**            knapvars,           /**< variables in the continuous knapsack constraint */
   SCIP_Real*            knapvals,           /**< coefficients of the variables in the continuous knapsack constraint */
   SCIP_Real             valscale,           /**< -1.0 if lhs of row is used as rhs of c. k. constraint, +1.0 otherwise */
   SCIP_Real             rhs,                /**< right hand side of the continuous knapsack constraint */
   SCIP_SOL*             sol,                /**< primal CIP solution, NULL for current LP solution */
   SCIP_Bool*            cutoff,             /**< pointer to store whether a cutoff was found */
   int*                  ncuts               /**< pointer to add up the number of found cuts */
   )
{
   SCIP_VAR** binvars;
   SCIP_VAR** consvars;
   SCIP_Real* binvals;
   SCIP_Longint* consvals;
   SCIP_Longint minact;
   SCIP_Longint maxact;
   SCIP_Real intscalar;
   SCIP_Bool success;
   int nbinvars;
   int nconsvars;
   int i;

   int* tmpindices;
   int tmp;
   SCIP_CONSHDLR* conshdlr;
   SCIP_CONSHDLRDATA* conshdlrdata;
   SCIP_Bool noknapsackconshdlr;
   SCIP_Bool usegubs;

   assert(nknapvars > 0);
   assert(knapvars != NULL);
   assert(cutoff != NULL);

   tmpindices = NULL;

   SCIPdebugMsg(scip, "separate linear constraint <%s> relaxed to knapsack\n", cons != NULL ? SCIPconsGetName(cons) : "-");
   SCIPdebug( if( cons != NULL ) { SCIPdebugPrintCons(scip, cons, NULL); } );

   binvars = SCIPgetVars(scip);

   /* all variables which are of integral type can be potentially of binary type; this can be checked via the method SCIPvarIsBinary(var) */
   nbinvars = SCIPgetNVars(scip) - SCIPgetNContVars(scip);

   *cutoff = FALSE;

   if( nbinvars == 0 )
      return SCIP_OKAY;

   /* set up data structures */
   SCIP_CALL( SCIPallocBufferArray(scip, &consvars, nbinvars) );
   SCIP_CALL( SCIPallocBufferArray(scip, &consvals, nbinvars) );

   /* get conshdlrdata to use cleared memory */
   conshdlr = SCIPfindConshdlr(scip, CONSHDLR_NAME);
   if( conshdlr == NULL )
   {
      noknapsackconshdlr = TRUE;
      usegubs = DEFAULT_USEGUBS;

      SCIP_CALL( SCIPallocBufferArray(scip, &binvals, nbinvars) );
      BMSclearMemoryArray(binvals, nbinvars);
   }
   else
   {
      noknapsackconshdlr = FALSE;
      conshdlrdata = SCIPconshdlrGetData(conshdlr);
      assert(conshdlrdata != NULL);
      usegubs = conshdlrdata->usegubs;

      SCIP_CALL( SCIPallocBufferArray(scip, &tmpindices, nknapvars) );

      /* increase array size to avoid an endless loop in the next block; this might happen if continuous variables
       * change their types to SCIP_VARTYPE_BINARY during presolving
       */
      if( conshdlrdata->reals1size == 0 )
      {
         SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &conshdlrdata->reals1, conshdlrdata->reals1size, 1) );
         conshdlrdata->reals1size = 1;
         conshdlrdata->reals1[0] = 0.0;
      }

      assert(conshdlrdata->reals1size > 0);

      /* next if condition should normally not be true, because it means that presolving has created more binary
       * variables than binary + integer variables existed at the constraint initialization method, but for example if you would
       * transform all integers into their binary representation then it maybe happens
       */
      if( conshdlrdata->reals1size < nbinvars )
      {
         int oldsize = conshdlrdata->reals1size;

         conshdlrdata->reals1size = nbinvars;
         SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &conshdlrdata->reals1, oldsize, conshdlrdata->reals1size) );
         BMSclearMemoryArray(&(conshdlrdata->reals1[oldsize]), conshdlrdata->reals1size - oldsize); /*lint !e866 */
      }
      binvals = conshdlrdata->reals1;

      /* check for cleared array, all entries have to be zero */
#ifndef NDEBUG
      for( tmp = nbinvars - 1; tmp >= 0; --tmp )
      {
         assert(binvals[tmp] == 0);
      }
#endif
   }

   tmp = 0;

   /* relax continuous knapsack constraint:
    * 1. make all variables binary:
    *    if x_j is continuous or integer variable substitute:
    *      - a_j < 0: x_j = lb  or  x_j = b*z + d with variable lower bound b*z + d with binary variable z
    *      - a_j > 0: x_j = ub  or  x_j = b*z + d with variable upper bound b*z + d with binary variable z
    * 2. convert coefficients of all variables to positive integers:
    *      - scale all coefficients a_j to a~_j integral
    *      - substitute  x~_j = 1 - x_j if a~_j < 0
    */

   /* replace integer and continuous variables with binary variables */
   for( i = 0; i < nknapvars; i++ )
   {
      SCIP_VAR* var;

      var = knapvars[i];

      if( SCIPvarIsBinary(var) && SCIPvarIsActive(var) )
      {
         SCIP_Real solval;
         assert(0 <= SCIPvarGetProbindex(var) && SCIPvarGetProbindex(var) < nbinvars);

         solval = SCIPgetSolVal(scip, sol, var);

         /* knapsack relaxation assumes solution values between 0.0 and 1.0 for binary variables */
         if( SCIPisFeasLT(scip, solval, 0.0 )
               || SCIPisFeasGT(scip, solval, 1.0) )
         {
            SCIPdebugMsg(scip, "Solution value %.15g <%s> outside domain [0.0, 1.0]\n",
                  solval, SCIPvarGetName(var));
            goto TERMINATE;
         }

         binvals[SCIPvarGetProbindex(var)] += valscale * knapvals[i];
         if( !noknapsackconshdlr )
         {
            assert(tmpindices != NULL);

            tmpindices[tmp] = SCIPvarGetProbindex(var);
            ++tmp;
         }
         SCIPdebugMsg(scip, " -> binary variable %+.15g<%s>(%.15g)\n", valscale * knapvals[i], SCIPvarGetName(var), SCIPgetSolVal(scip, sol, var));
      }
      else if( valscale * knapvals[i] > 0.0 )
      {
         SCIP_VAR** zvlb;
         SCIP_Real* bvlb;
         SCIP_Real* dvlb;
         SCIP_Real bestlbsol;
         int bestlbtype;
         int nvlb;
         int j;

         /* a_j > 0: substitution with lb or vlb */
         nvlb = SCIPvarGetNVlbs(var);
         zvlb = SCIPvarGetVlbVars(var);
         bvlb = SCIPvarGetVlbCoefs(var);
         dvlb = SCIPvarGetVlbConstants(var);

         /* search for lb or vlb with maximal bound value */
         bestlbsol = SCIPvarGetLbGlobal(var);
         bestlbtype = -1;
         for( j = 0; j < nvlb; j++ )
         {
            /* use only numerical stable vlb with binary variable z */
            if( SCIPvarIsBinary(zvlb[j]) && SCIPvarIsActive(zvlb[j]) && REALABS(bvlb[j]) <= MAXABSVBCOEF )
            {
               SCIP_Real vlbsol;

               if( (bvlb[j] >= 0.0 && SCIPisGT(scip, bvlb[j] * SCIPvarGetLbLocal(zvlb[j]) + dvlb[j], SCIPvarGetUbLocal(var))) ||
                   (bvlb[j] <= 0.0 && SCIPisGT(scip, bvlb[j] * SCIPvarGetUbLocal(zvlb[j]) + dvlb[j], SCIPvarGetUbLocal(var))) )
               {
                  *cutoff = TRUE;
                  SCIPdebugMsg(scip, "variable bound <%s>[%g,%g] >= %g<%s>[%g,%g] + %g implies local cutoff\n",
                     SCIPvarGetName(var), SCIPvarGetLbLocal(var), SCIPvarGetUbLocal(var),
                     bvlb[j], SCIPvarGetName(zvlb[j]), SCIPvarGetLbLocal(zvlb[j]), SCIPvarGetUbLocal(zvlb[j]), dvlb[j]);
                  goto TERMINATE;
               }

               assert(0 <= SCIPvarGetProbindex(zvlb[j]) && SCIPvarGetProbindex(zvlb[j]) < nbinvars);
               vlbsol = bvlb[j] * SCIPgetSolVal(scip, sol, zvlb[j]) + dvlb[j];
               if( SCIPisGE(scip, vlbsol, bestlbsol) )
               {
                  bestlbsol = vlbsol;
                  bestlbtype = j;
               }
            }
         }

         /* if no lb or vlb with binary variable was found, we have to abort */
         if( SCIPisInfinity(scip, -bestlbsol) )
            goto TERMINATE;

         if( bestlbtype == -1 )
         {
            rhs -= valscale * knapvals[i] * bestlbsol;
            SCIPdebugMsg(scip, " -> non-binary variable %+.15g<%s>(%.15g) replaced with lower bound %.15g (rhs=%.15g)\n",
               valscale * knapvals[i], SCIPvarGetName(var), SCIPgetSolVal(scip, sol, var), SCIPvarGetLbGlobal(var), rhs);
         }
         else
         {
            assert(0 <= SCIPvarGetProbindex(zvlb[bestlbtype]) && SCIPvarGetProbindex(zvlb[bestlbtype]) < nbinvars);
            rhs -= valscale * knapvals[i] * dvlb[bestlbtype];
            binvals[SCIPvarGetProbindex(zvlb[bestlbtype])] += valscale * knapvals[i] * bvlb[bestlbtype];

            if( SCIPisInfinity(scip, REALABS(binvals[SCIPvarGetProbindex(zvlb[bestlbtype])])) )
               goto TERMINATE;

            if( !noknapsackconshdlr )
            {
               assert(tmpindices != NULL);

               tmpindices[tmp] = SCIPvarGetProbindex(zvlb[bestlbtype]);
               ++tmp;
            }
            SCIPdebugMsg(scip, " -> non-binary variable %+.15g<%s>(%.15g) replaced with variable lower bound %+.15g<%s>(%.15g) %+.15g (rhs=%.15g)\n",
               valscale * knapvals[i], SCIPvarGetName(var), SCIPgetSolVal(scip, sol, var),
               bvlb[bestlbtype], SCIPvarGetName(zvlb[bestlbtype]),
               SCIPgetSolVal(scip, sol, zvlb[bestlbtype]), dvlb[bestlbtype], rhs);
         }