
#include <assert.h>
#include <string.h>

#include <scip/cons_knapsack.h>
#include "handle_localcuts.h"
#include "sepa_localcuts.h"
#include "oracle/lmo_api.h"

#ifdef USE_FW_SEPA
#include "scalar-fw/fwknapcut.h"
#endif


#define SEPA_NAME              "localcuts"
#define SEPA_DESC              "local cut separator"
#define SEPA_PRIORITY           -100000
#define SEPA_FREQ                    -1
#define SEPA_MAXBOUNDDIST           1.0
#define SEPA_USESUBSCIP            TRUE /**< to avoid being called in subscips */
#define SEPA_DELAY                 TRUE /**< should separation method be delayed, if other separators found cuts? */


/* default parameters */
#define DEFAULT_THRESHOLD                    0
#define DEFAULT_THRESHOLDROOT                0
#define DEFAULT_MAXROUNDS                    5
#define DEFAULT_MAXROUNDSROOT               40
#define DEFAULT_NICECUTS                  TRUE
#define DEFAULT_TILTING                   TRUE
#define DEFAULT_LIFTING                  FALSE
#define DEFAULT_UPLIFTING                FALSE
#define DEFAULT_USECMIR                  FALSE
#define DEFAULT_ONLYCHECKCONSS            TRUE

/* other parameters */
#define MAX_TABLE_SIZE                 1000000 /**< maximal table size for lfiting */
#define BOUNDSWITCH                       0.51 /**< threshold for bound switching - see cuts.c */
#define POSTPROCESS                      FALSE /**< apply postprocessing to the cut - see cuts.c */
#define USEVBDS                          FALSE /**< use variable bounds - see cuts.c */
#define ALLOWLOCAL                       FALSE /**< allow to generate local cuts - see cuts. */
#define MINFRAC                          0.05  /**< minimal fractionality of floor(rhs) - see cuts.c */
#define MAXFRAC                          0.999 /**< maximal fractionality of floor(rhs) - see cuts.c */

/** separator data */
struct SCIP_SepaData
{
   int                   threshold;          /**< threshold for local cut separation */
   int                   thresholdroot;      /**< threshold for local cut separation (root) */
   int                   maxrounds;          /**< max. number of local cut rounds */
   int                   maxroundsroot;      /**< max. number of local cut rounds in the root */
   SCIP_Bool             nicecuts;           /**< Try to make a nice cut by finding a rational approximation? */
   SCIP_Bool             tilting;            /**< Try tilting to strengthen cuts? */
   SCIP_Bool             lifting;            /**< Use lifting? */
   SCIP_Bool             uplifting;          /**< Perfrom uplifting (if lifting is true)? */
   SCIP_Bool             usecmir;            /**< Additionally apply CMIR functions to generated cuts? */
   SCIP_Bool             onlycheckconss;     /**< Only run separation on constraints that are checked? */

   SCIP_CONSHDLR*        knapsackconshdlr;   /**< Knapsack constraint handler */
   int                   nrounds;            /**< counter for local cut rounds per node */
   SCIP_CUTPOOL*         pool;               /**< pool of generated local cuts (globally valid) */
   SCIP_Longint          nsepacuts;          /**< total number of separated local cuts */
   int                   ngencmir;           /**< number of generated CMIR cuts */
   SCIP_Longint          currentnode;        /**< id of current node */
   int                   ntotalgen;          /**< total number of generated cuts */

   SCIP_Longint          cummulativereddim;  /**< helper to compute average submatrix dimension */
   SCIP_Longint          cummulativeredcapacity; /**< helper to compute average submatrix capacity */
   SCIP_Longint          maxreddim;          /**< helper to compute maximum submatrix dimension */
   SCIP_Longint          maxredcapacity;     /**< helper to compute maximum submatrix capacity */
   SCIP_Longint          nreddimcalls;       /**< helper to compute averages */

   int                   ntotalsepacalls;    /**< total number of separation calls that actually perform work */
   SCIP_Longint          ntotaloraclecalls;  /**< total number of calls to the oracle */
   SCIP_Longint          ntotalfwiters;      /**< total number of FW iterations */
   int                   nexit_max_iter;     /**< number of times FW terminated reaching the maximal iteration bound */
   int                   nexit_early;        /**< number of times FW terminated early */
   int                   nexit_zero_gradient;/**< number of times FW terminated with a zero gradient */
   int                   nexit_term_check;   /**< number of times FW terminated by the termination check */
   int                   nexit_primal_gap;  /**< number of times FW terminated reaching the primal gap bound */
#ifdef USE_FW_SEPA
   LC_CLOCK*             fwclock;            /**< clock for measuring FW time */
   LC_CLOCK*             oracleclock;        /**< clock for measuring oracle time */
#endif
   int                   npoolcuts;          /**< total number of cuts in the pool */
   SCIP_Real             pooltime;           /**< total separation time for pool */
   SCIP_Longint          npoolcalls;         /**< number of calls to pool separation */
   SCIP_Longint          npoolcutsfound;     /**< number of cuts found by pool separation */
   SCIP_Longint          npoolcutsadded;     /**< number of cuts added by pool separation */
   SCIP_Longint          nmaxpoolcuts;       /**< maximal number of cuts in the pool */
};


/** data for local cut separator */
struct LC_Data
{
   int                   nvars;
   SCIP_Longint*         weights;
   SCIP_Longint          capacity;
   SCIP_VAR**            vars;
};



/*
 * Local methods
 */



#if ! defined(NDEBUG) || ! defined(USE_FW_SEPA)

/** Optimization oracle maximizing knapsack
 *
 * LC_ORACLE(x) SCIP_RETCODE x (SCIP* scip, unsigned int dim, const LC_DATA* data, \
 *     const SCIP_Real* obj, SCIP_Real cutoff, unsigned int maxnsols, unsigned int* nsols, SCIP_Real** sols, SCIP_Real* objvals, \
 *     LC_STATUS* status)
 */ /*lint --e{715}*/
static
LC_ORACLE(OracleKnapsack)
{
   int i;

   assert( scip != NULL );
   assert( data != NULL );
   assert( obj != NULL );
   assert( status != NULL );

   assert( maxnsols >= 1 );
   assert( nsols != NULL );
   assert( objvals != NULL );

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

   *status = LC_STATUS_UNKNOWN;
   *nsols = 0;

#if 1
   {
      double* vertex;
      SCIP_CALL( SCIPallocBufferArray(scip, &vertex, dim) );
      *objvals = callKnapsackOracle(scip, dim, data->weights, data->capacity, obj, vertex);

      if ( *objvals > cutoff )
      {
         *nsols = 1;
         for (i = 0; i < dim; ++i)
         {
            assert( SCIPisFeasEQ(scip, vertex[i], 0.0) || SCIPisFeasEQ(scip, vertex[i], 1.0) );
            sols[0][i] = vertex[i];
         }
         *status = LC_STATUS_OPTIMAL;
      }
      else
         *status = LC_STATUS_CUTOFF;

      SCIPfreeBufferArray(scip, &vertex);
   }
#else
   {
      SCIP_Bool success;
      SCIP_Real solval;
      int* solitems;
      int* nonsolitems;
      int* items;
      int nsolitems;
      int nnonsolitems;

      SCIP_CALL( SCIPallocBufferArray(scip, &items, data->nvars) );
      for (i = 0; i < data->nvars; ++i)
         items[i] = i;
      SCIP_CALL( SCIPallocBufferArray(scip, &solitems, data->nvars) );
      SCIP_CALL( SCIPallocBufferArray(scip, &nonsolitems, data->nvars) );

      SCIP_CALL( SCIPsolveKnapsackExactly(scip, data->nvars, data->weights, (double*) obj, data->capacity, items, solitems, nonsolitems, &nsolitems, &nnonsolitems, &solval, &success) );

      if ( success )
      {
         if ( solval > cutoff )
         {
            *nsols = 1;
            for (i = 0; i < data->nvars; ++i)
               sols[0][i] = 0.0;
            for (i = 0; i < nsolitems; ++i)
               sols[0][solitems[i]] = 1.0;
            *objvals = solval;

            *status = LC_STATUS_OPTIMAL;
         }
         else
            *status = LC_STATUS_CUTOFF;
      }

      SCIPfreeBufferArray(scip, &solitems);
      SCIPfreeBufferArray(scip, &nonsolitems);
      SCIPfreeBufferArray(scip, &items);
   }
#endif

   return SCIP_OKAY;
}
#endif

/** create and add local cut */
static
SCIP_RETCODE addLocalCut(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_SEPA*            sepa,               /**< separator */
   SCIP_SEPADATA*        sepadata,           /**< separator data */
   SCIP_SOL*             sol,                /**< solution to be separated */
   SCIP_VAR**            knapvars,           /**< variables in knapsack constraint */
   int                   nknapvars,          /**< number of variables in knapsack constraint */
   SCIP_Real*            cutcoefs,           /**< array containing coefficients of cut */
   SCIP_Real             cutrhs,             /**< rhs of cut */
   SCIP_Bool             cutislocal,         /**< whether cut is local */
   SCIP_Real             cutact,             /**< activity of cut */
   SCIP_Real*            cmircoefs,          /**< temporary array for creating CMIR cuts (or NULL) */
   int*                  cmirinds,           /**< temporary array for creating CMIR cuts (or NULL) */
   int*                  nlc,                /**< number of generated local cuts */
   SCIP_Bool*            infeasible,         /**< pointer to store whether we are infeasible */
   SCIP_Bool*            success             /**< whether cut has actually been lifted */
   )
{
   char name[SCIP_MAXSTRLEN];
   SCIP_Bool generateliftcut = TRUE;
   SCIP_ROW* row;
   int i;

   assert( infeasible != NULL );
   assert( success != NULL );

   *infeasible = FALSE;
   *success = FALSE;

   /* possibly apply CMIR routine */
   if ( sepadata->usecmir )
   {
      SCIP_AGGRROW* aggrrow;
      SCIP_VAR** scipvars;
      SCIP_Real cmirefficacy;
      SCIP_Real cmirrhs;
      int cmirnnz = 0;

      SCIP_CALL( SCIPaggrRowCreate(scip, &aggrrow) );

      scipvars = SCIPgetVars(scip);

      /* set up row data */
      cmirrhs = cutrhs;
      for (i = 0; i < nknapvars; ++i)
      {
         if ( ! SCIPisZero(scip, cutcoefs[i]) )
         {
            SCIP_VAR* var;

            var = knapvars[i];
            if ( SCIPvarIsNegated(var) )
            {
               var = SCIPvarGetNegationVar(var);
               cmircoefs[cmirnnz] = - cutcoefs[i];
               cmirrhs -= cutcoefs[i];
            }
            else
               cmircoefs[cmirnnz] = cutcoefs[i];
            cmirinds[cmirnnz] = SCIPvarGetProbindex(var);
            assert( cmirinds[cmirnnz] >= 0 );
            ++cmirnnz;
         }
      }
      assert( cmirnnz > 0 );

      /* add produced row as custom row */
      SCIP_CALL( SCIPaggrRowAddCustomCons(scip, aggrrow, cmirinds, cmircoefs, cmirnnz, cmirrhs, 1.0, 0, cutislocal) );

      /* try to generate CMIR inequality */
      cmirefficacy = (cutact - cmirrhs) / SCIPaggrRowCalcEfficacyNorm(scip, aggrrow);
      SCIP_CALL( SCIPcutGenerationHeuristicCMIR(scip, sol, POSTPROCESS, BOUNDSWITCH, USEVBDS, ALLOWLOCAL, INT_MAX, NULL, NULL,
            MINFRAC, MAXFRAC, aggrrow, cmircoefs, &cmirrhs, cmirinds, &cmirnnz, &cmirefficacy, NULL, &cutislocal, success) );

      if ( success )
      {
         /* the lifted cuts are globally valid, i.e., "local = FALSE" */
         (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "lcliftcmir%d", sepadata->ntotalgen);
         SCIP_CALL( SCIPcreateEmptyRowSepa(scip, &row, sepa, name, -SCIPinfinity(scip), cmirrhs, cutislocal, FALSE, TRUE) );
         SCIP_CALL( SCIPcacheRowExtensions(scip, row) );
         for (i = 0; i < cmirnnz; ++i)
         {
            SCIP_CALL( SCIPaddVarToRow(scip, row, scipvars[cmirinds[i]], cmircoefs[i]) );
         }
         SCIP_CALL( SCIPflushRowExtensions(scip, row) );

         /* only add row if it is efficacious with respect to norm */
         if ( SCIPisCutEfficacious(scip, sol, row) )
         {
            SCIPdebugMsg(scip, "Found efficious lifted CMIR cut with efficacy %f (feasibility = %f, norm = %f, nnz = %d, nknapvars = %d)\n",
               SCIPgetCutEfficacy(scip, sol, row), SCIPgetRowLPFeasibility(scip, row), SCIProwGetNorm(row), SCIProwGetNNonz(row), nknapvars);

            SCIP_CALL( SCIPaddRow(scip, row, FALSE, infeasible) );
#ifdef SCIP_MORE_DEBUG
            SCIPdebug( SCIPprintRow(scip, row, NULL) );
#endif

            ++sepadata->ntotalgen;
            ++sepadata->ngencmir;
            ++(*nlc);

            /* add cut to pool if required (pool != NULL) */
            if ( sepadata->pool != NULL && ! cutislocal )
            {
               SCIP_CALL( SCIPaddRowCutpool(scip, sepadata->pool, row) );
            }

            /* do not generate original cut, if we were successful */
            generateliftcut = FALSE;
         }
         SCIP_CALL( SCIPreleaseRow(scip, &row) );
      }

      SCIPaggrRowFree(scip, &aggrrow);
   }

   if ( generateliftcut && !(*infeasible) )
   {
      /* the lifted cuts are globally valid, i.e., "local = FALSE" */
      (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "lclift%d", sepadata->ntotalgen);
      SCIP_CALL( SCIPcreateEmptyRowSepa(scip, &row, sepa, name, -SCIPinfinity(scip), cutrhs, cutislocal, FALSE, TRUE) );
      SCIP_CALL( SCIPaddVarsToRow(scip, row, nknapvars, knapvars, cutcoefs) );

      /* only add row if it is efficacious with respect to norm */
      if ( SCIPisCutEfficacious(scip, sol, row) )
      {
         SCIPdebugMsg(scip, "Found efficious cut with efficacy %f (feasibility = %f, norm = %f, nnz = %d, nknapvars = %d)\n",
            SCIPgetCutEfficacy(scip, sol, row), SCIPgetRowLPFeasibility(scip, row), SCIProwGetNorm(row), SCIProwGetNNonz(row), nknapvars);

         SCIP_CALL( SCIPaddRow(scip, row, FALSE, infeasible) );
#ifdef SCIP_MORE_DEBUG
         SCIPdebug( SCIPprintRow(scip, row, NULL) );
#endif

         ++sepadata->ntotalgen;
         ++(*nlc);

         /* add cut to pool if required (pool != NULL) */
         if ( sepadata->pool != NULL && ! cutislocal )
         {
            SCIP_CALL( SCIPaddRowCutpool(scip, sepadata->pool, row) );
         }
      }
      SCIP_CALL( SCIPreleaseRow(scip, &row) );
   }

   return SCIP_OKAY;
}

/** Lifting function
 *
 * Lifts given inequality \f$ \sum_{j \in S} \alpha_j x_j \leq \alpha_0 \f$ valid for
 * \f[
 *     S^0 = \{ x \in \{0,1\}^{S} : \sum_{j \in S} a_j x_j \leq a_0 - \sum_{j \in F_1} a_j \}
 * \f]
 * to a valid inequality
 * \f[
 *     \sum_{j \in S} \alpha_j x_j + \sum_{j \in F_1} \beta_j x_j + \sum_{j \in F_0} \beta_j x_j
 *     \leq \alpha_0 + \sum_{j \in F_1} \beta_j
 * \f]
 * for \f$ K = \{ x \in \{0,1\}^{N} : \sum_{j \in N} a_j x_j \leq a_0 \} \f$.  Uses sequential down-lifting for the
 * variable in \f$F_1\f$, and sequential up-lifting for the variables in \f$F_0\f$.
 *
 * This function is adapted from cons_knapsack.
 */
static
SCIP_RETCODE sequentialUpAndDownLifting(
   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,           /**< original capacity of knapsack */
   SCIP_Real*            solvals,            /**< solution values of all problem variables */
   int*                  varsS,              /**< variables in S (already in the cut) */
   SCIP_Real*            coefsS,             /**< cut coefficients for variables in S */
   int*                  varsF1,             /**< variables in F_1 (down-lifted) */
   int*                  varsF0,             /**< variables in F_0 (up-lifted) */
   int                   nvarsS,             /**< number of variables in S (already in the cut) */
   int                   nvarsF1,            /**< number of variables in F_1 */
   int                   nvarsF0,            /**< number of variables in F_2 */
   SCIP_Real             alpha0,             /**< right hand side of given valid inequality */
   SCIP_Real*            liftcoefs,          /**< array to store lifting coefficient of vars in knapsack constraint */
   SCIP_Real*            cutact,             /**< pointer to store activity of lifted valid inequality */
   SCIP_Real*            liftrhs,            /**< pointer to store right hand side of the lifted valid inequality */
   SCIP_Bool*            success             /**< whether cut has actually been lifted */
   )
{
   SCIP_Real* table;
   SCIP_Longint fixedonesweight;
   SCIP_Longint tablesize;
   int j;

   assert( scip != NULL );
   assert( vars != NULL );
   assert( nvars >= 0 );
   assert( weights != NULL );
   assert( capacity >= 0 );
   assert( solvals != NULL );
   assert( varsS != NULL );
   assert( coefsS != NULL );
   assert( varsF1 != NULL );
   assert( 0 <= nvarsS && nvarsS <= nvars );
   assert( 0 <= nvarsF1 && nvarsF1 <= nvars );
   assert( 0 <= nvarsF0 && nvarsF0 <= nvars );
   assert( alpha0 >= 0.0 );
   assert( liftcoefs != NULL );
   assert( cutact != NULL );
   assert( liftrhs != NULL );
   assert( success != NULL );

   *success = TRUE;

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

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

   if ( capacity > MAX_TABLE_SIZE )
   {
      *success = FALSE;
      return SCIP_OKAY;
   }

   /* allocates memory for DP table */
   tablesize = capacity + 1;
   SCIP_CALL( SCIPallocBufferArray(scip, &table, tablesize) );
   BMSclearMemoryArray(table, tablesize);

   /* fill in table using items in S */
   for (j = 0; j < nvarsS; ++j)
   {
      SCIP_Longint weight;
      SCIP_Real coef;
      SCIP_Longint w;

      assert( 0 <= varsS[j] && varsS[j] < nvars );
      weight = weights[varsS[j]];
      assert( weight > 0);

      coef = coefsS[j];

      for (w = tablesize - weight - 1; w >= 0; --w)
      {
         if ( table[w] + coef > table[w + weight] )
            table[w + weight] = table[w] + coef;
      }

      assert( liftcoefs[varsS[j]] == 0.0 );
      liftcoefs[varsS[j]] = coef;
      *cutact += coef * solvals[varsS[j]];
   }

   /* gets sum of weights of variables fixed to 1, i.e., sum of weights of variables in F_1 */
   fixedonesweight = 0;
   for (j = 0; j < nvarsF1; ++j)
   {
      assert( 0 <= varsF1[j] && varsF1[j] < nvars );
      fixedonesweight += weights[varsF1[j]];
   }
   assert( fixedonesweight >= 0 );

   /* sequentially down-lifts all variables in F_1: */
   for (j = 0; j < nvarsF1; ++j)
   {
      SCIP_Longint weight;
      SCIP_Real liftcoef;
      SCIP_Longint w;
      int liftvar;

      liftvar = varsF1[j];
      assert( 0 <= liftvar && liftvar < nvars );

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

      assert( liftcoefs[liftvar] == 0.0 );
      assert( capacity - fixedonesweight + weight >= 0 );
      liftcoef = table[capacity - fixedonesweight + weight] - (*liftrhs);
      liftcoefs[liftvar] = liftcoef;
      assert( liftcoef >= 0 );

      /* update rhs */
      *liftrhs += liftcoef;
      fixedonesweight -= weight;

      /* update table */
      if ( nvarsF0 > 0 || j < nvarsF1 - 1 )
      {
         for (w = tablesize - weight - 1; w >= 0; --w)
         {
            if ( table[w] + liftcoef > table[w + weight] )
               table[w + weight] = table[w] + liftcoef;
         }
      }
   }
   assert( *liftrhs >= alpha0 );
   assert( fixedonesweight == 0 );

   /* sequentially up-lifts all variables in F_0: */
   for (j = 0; j < nvarsF0; ++j)
   {
      SCIP_Longint weight;
      SCIP_Real liftcoef;
      SCIP_Longint w;
      int liftvar;

      assert( varsF0 != NULL );
      liftvar = varsF0[j];
      assert( 0 <= liftvar && liftvar < nvars);
      weight = weights[liftvar];

      /* variables are in varsF0 if their value is 0 or if they do not fit into the reduced knapsack */
      assert( 0 < weight && weight <= capacity );

      liftcoef = (*liftrhs) - table[capacity - weight];
      liftcoefs[liftvar] = liftcoef;
      assert( 0 <= liftcoef && liftcoef <= *liftrhs );

      /* update table */
      if ( j < nvarsF0 - 1 )
      {
         for (w = tablesize - weight - 1; w >= 0; --w)
         {
            if ( table[w] + liftcoef > table[w + weight] )
               table[w + weight] = table[w] + liftcoef;
         }
      }
   }

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

   return SCIP_OKAY;
}


/** separate local cuts with lifting
 *
 *  We create a reduced knapsack by removing variables of each knapsack that have integral values (0/1) in the
 *  LP-relaxation. In general, cuts on the reduced knapsack are only valid if we do not reduce the capacity, i.e., if no
 *  solution value of a variable is 1. In this function, we always down-lift these variables, so that the resulting cut
 *  is globally valid.
 */
static
SCIP_RETCODE separateLocalCutsKnapsackLifting(
   SCIP*                 scip,               /**< scip pointer */
   SCIP_SEPA*            sepa,               /**< separator */
   SCIP_SEPADATA*        sepadata,           /**< data of separator */
   SCIP_SOL*             sol,                /**< SCIP solution to be separated (NULL if current LP-solution) */
   int*                  nlc,                /**< number of generated local cuts */
   SCIP_Bool*            infeasible          /**< pointer to store whether infeasibility was detected */
   )
{
   SCIP_CONS** conss;
   SCIP_Real* xrelint = NULL;
   SCIP_Real* x = NULL;
   SCIP_Longint* weights = NULL;
   SCIP_Real* solvals;
   SCIP_Real* cutcoefs;
   SCIP_VAR** vars;
   SCIP_Real* liftcoefs;
   SCIP_Real* cmircoefs = NULL;
   int* cmirinds = NULL;
   int* redvars = NULL;   /* variable indices in reduced knapsack (i.e., fractional variables in knapsack) */
   int* varsF1 = NULL;    /* variable indices in knapsack that are fixed to 1 */
   int* varsF0 = NULL;    /* variable indices in knapsack that are fixed to 0 */
   SCIP_Real liftrhs;
   int nconss;
   int nvars;
   int c;

   assert( scip != NULL );
   assert( sepa != NULL );
   assert( sepadata != NULL );
   assert( nlc != NULL );
   assert( infeasible != NULL );

   *nlc = 0;
   *infeasible = FALSE;

   nvars = SCIPgetNVars(scip);
   nconss = SCIPconshdlrGetNConss(sepadata->knapsackconshdlr);
   if ( nconss == 0 || nvars == 0 )
      return SCIP_OKAY;
   assert( nconss > 0 );
   assert( nvars > 0 );

   SCIP_CALL( SCIPallocBufferArray(scip, &liftcoefs, nvars) );
   SCIP_CALL( SCIPallocBufferArray(scip, &cutcoefs, nvars) );
   SCIP_CALL( SCIPallocBufferArray(scip, &vars, nvars) );
   SCIP_CALL( SCIPallocBufferArray(scip, &redvars, nvars) );
   SCIP_CALL( SCIPallocBufferArray(scip, &varsF1, nvars) );
   SCIP_CALL( SCIPallocBufferArray(scip, &varsF0, nvars) );
   SCIP_CALL( SCIPallocBufferArray(scip, &weights, nvars) );
   SCIP_CALL( SCIPallocBufferArray(scip, &x, nvars) );
   SCIP_CALL( SCIPallocBufferArray(scip, &xrelint, nvars) );
   SCIP_CALL( SCIPallocBufferArray(scip, &solvals, nvars) );

   if ( sepadata->usecmir )
   {
      SCIP_CALL( SCIPallocBufferArray(scip, &cmircoefs, nvars) );
      SCIP_CALL( SCIPallocBufferArray(scip, &cmirinds, nvars) );
   }

   /* loop through knapsack constraints */
   conss = SCIPconshdlrGetConss(sepadata->knapsackconshdlr);
   for (c = 0; c < nconss && ! SCIPisStopped(scip); ++c)
   {
      SCIP_CONS* cons;
      SCIP_VAR** knapvars;
      SCIP_Longint* knapweights;
      SCIP_Longint redcapacity;
      SCIP_Longint capacity;
      SCIP_Longint weightsum = 0;
      SCIP_Bool somefractional = FALSE;
      SCIP_Real cutact;
      SCIP_Bool success;
      SCIP_Real cutrhs = SCIP_INVALID;
      int infeasibleitems = 0;
      int nknapvars;
      int nvarsF1 = 0;
      int nvarsF0 = 0;
      int dim = 0;
      int i;

      cons = conss[c];
      assert( cons != NULL );

      if ( sepadata->onlycheckconss && ! SCIPconsIsChecked(cons) )
         continue;

      nknapvars = SCIPgetNVarsKnapsack(scip, cons);
      knapvars = SCIPgetVarsKnapsack(scip, cons);
      knapweights = SCIPgetWeightsKnapsack(scip, cons);
      capacity = SCIPgetCapacityKnapsack(scip, cons);

      /* first compute reduced capacity */
      redcapacity = capacity;
      for (i = 0; i < nknapvars; ++i)
      {
         assert( SCIPvarIsBinary(knapvars[i]) );
         assert( knapweights[i] >= 0 );

         solvals[i] = SCIPgetSolVal(scip, sol, knapvars[i]);
         assert( SCIPisFeasGE(scip, solvals[i], 0.0) && SCIPisFeasLE(scip, solvals[i], 1.0) );

         /* if value is 1, reduce the capacity */
         if ( SCIPisFeasEQ(scip, solvals[i], 1.0) )
            redcapacity -= knapweights[i];
      }

      /* if reduced problem is empty, skip separation */
      if ( redcapacity <= 0 )
         continue;

      /* now collect reduced knapsack */
      for (i = 0; i < nknapvars; ++i)
      {
         if ( SCIPisFeasEQ(scip, solvals[i], 1.0) )
            varsF1[nvarsF1++] = i;
         else if ( SCIPisFeasEQ(scip, solvals[i], 0.0) )
            varsF0[nvarsF0++] = i;
         else
         {
            /* if weights are larger than the reduced capacity, continue to collect varsF0/varsF1 and create a special cut below */
            if ( knapweights[i] > redcapacity )
               ++infeasibleitems;
            else
            {
               somefractional = TRUE;
               weightsum += knapweights[i];

               x[dim] = solvals[i];
               vars[dim] = knapvars[i];
               weights[dim] = knapweights[i];
               redvars[dim++] = i;
            }
         }
      }
      assert( dim + nvarsF0 + nvarsF1 + infeasibleitems == nknapvars );

      /* skip knapsacks with all variables having integral values */
      if ( ! somefractional )
         continue;

      /* skip trivially feasible knapsacks */
      if ( weightsum <= redcapacity )
         continue;

      if ( infeasibleitems > 0 )
      {
         /* collect variables that are infeasible to create the cut \sum_{i infeasible} x_i \leq 0, see Vasilyev et al. [2016] */
         dim = 0;
         for (i = 0; i < nknapvars; ++i)
         {
            if ( knapweights[i] > redcapacity )
            {
               if ( ! SCIPisFeasEQ(scip, solvals[i], 1.0) && ! SCIPisFeasEQ(scip, solvals[i], 0.0) )
               {
                  cutcoefs[dim] = 1.0;
                  redvars[dim++] = i;
               }
            }
         }
         cutrhs = 0.0;
         success = TRUE;
         SCIPdebugMsg(scip, "Found cut that fixes some variables to 0 for knapsack constraint <%s> (dim = %d) ...\n", SCIPconsGetName(cons), dim);
         assert( dim == infeasibleitems );
      }
      else if ( dim == 2 )
      {
         /* if there are only two variables and the knapsack is not trivial, there is only one cut x1 + x2 \leq 1. */
         for (i = 0; i < dim; ++i)
            cutcoefs[i] = 1.0;
         cutrhs = 1.0;
         success = TRUE;
         SCIPdebugMsg(scip, "Found reduced knapsack with two variables for knapsack constraint <%s> (dim = %d) ...\n", SCIPconsGetName(cons), dim);
      }
      else
      {
         /* Skip separation in trivial cases:
          * 1) there are no variables left (dim = 0).
          * 2) dim = 1, then the problem is either trivially valid or we have infeasible items (this case is treated above)
          */
         if ( dim <= 1 )
            continue;
         assert( dim > 1 );

         SCIPdebugMsg(scip, "Separating local cuts for knapsack constraint <%s> (dim = %d) ...\n", SCIPconsGetName(cons), dim);

         /* set up interior point (reduced problem should be full-dimensional now) */
         for (i = 0; i < dim; ++i)
            xrelint[i] = 1.0 / ((SCIP_Real)dim * redcapacity + 1.0);

#ifdef USE_FW_SEPA
         {
            FWTuple* fwtpl;
            uint64_t max_fw_iters = 1000;

            /* increase number of iterations in root node */
            if ( SCIPgetDepth(scip) == 0 )
               max_fw_iters = 10000;

            /* SCIP_CALL( SCIPgetIntParam(scip, "separating/" SEPA_NAME "/maxfwiters", &max_fw_iters) ); */

            fwtpl = scalar_afw_knap_gen_cut(scip, dim, weights, redcapacity, x, xrelint, sepadata->fwclock, sepadata->oracleclock, 1e-6, max_fw_iters);

            // record statistics
            sepadata->ntotaloraclecalls += fwtpl->noraclecalls;
            sepadata->ntotalfwiters += fwtpl->niters;
            sepadata->cummulativeredcapacity += redcapacity;
            sepadata->cummulativereddim += dim;
            sepadata->nreddimcalls++;

            if ( redcapacity > sepadata->maxredcapacity )
               sepadata->maxredcapacity = redcapacity;

            if ( dim > sepadata->maxreddim )
               sepadata->maxreddim = dim;

            switch ( fwtpl->exitcode )
            {
            case exit_max_iters: ++sepadata->nexit_max_iter; break;
            case exit_early: ++sepadata->nexit_early; break;
            case exit_zero_gradient: ++sepadata->nexit_zero_gradient; break;
            case exit_term_check: ++sepadata->nexit_term_check; break;
            case exit_primal_gap: ++sepadata->nexit_primal_gap; break;
            case exit_not_run:
            default: abort();
            }

            if ( fwtpl->success )
            {
               success = TRUE;
               for (i = 0; i < dim; ++i)
                  cutcoefs[i] = fwtpl->a[i];
               cutrhs = fwtpl->b;
            }
            else
               success = FALSE;

            free(fwtpl->a);
            free(fwtpl);
         }
#else
         {
            SCIP_Real usedtime;
            LC_DATA data;

            data.nvars = dim;
            data.vars = vars;
            data.weights = weights;
            data.capacity = redcapacity;

            /* call separation */
            SCIP_CALL( LC_generateCut(scip, data.nvars, x, xrelint, data.vars, OracleKnapsack, &data, 0, NULL, NULL,
                  0, NULL, NULL, NULL, NULL, NULL, FALSE, sepadata->tilting, FALSE, sepadata->nicecuts,
                  cutcoefs, &cutrhs, &success, NULL, &usedtime) );
         }
#endif
      }

      /* call lifting */
      if ( success )
      {
         assert( cutrhs != SCIP_INVALID );

         if ( sepadata->uplifting )
         {
            SCIPdebugMsg(scip, "Up- and downlifting cut ...\n");
            SCIP_CALL( sequentialUpAndDownLifting(scip, knapvars, nknapvars, knapweights, capacity, solvals, redvars, cutcoefs, varsF1,
                  varsF0, dim, nvarsF1, nvarsF0, cutrhs, liftcoefs, &cutact, &liftrhs, &success) );
         }
         else
         {
            SCIPdebugMsg(scip, "Downlift cut ...\n");
            SCIP_CALL( sequentialUpAndDownLifting(scip, knapvars, nknapvars, knapweights, capacity, solvals, redvars, cutcoefs, varsF1,
                  NULL, dim, nvarsF1, 0, cutrhs, liftcoefs, &cutact, &liftrhs, &success) );
         }

         /* create cut if successfully lifted */
         if ( success )
         {
#ifndef NDEBUG
            int noraclesols = 0;
            SCIP_Real* oraclesol;
            SCIP_Real oracleobj;
            LC_STATUS status;
            LC_DATA data;

            data.nvars = nknapvars;
            data.vars = knapvars;
            data.weights = knapweights;
            data.capacity = capacity;
            SCIP_CALL( SCIPallocBufferArray(scip, &oraclesol, data.nvars) );

            /* check validity of generated cut */
            SCIP_CALL( OracleKnapsack(scip, data.nvars, &data, liftcoefs, -SCIPinfinity(scip), 1, &noraclesols, &oraclesol, &oracleobj, &status) );
            SCIPdebugMsg(scip, "Lifted rhs = %g, optimized rhs = %g\n", liftrhs, oracleobj);
            assert( SCIPisEQ(scip, oracleobj, liftrhs) );
            SCIPfreeBufferArray(scip, &oraclesol);
#endif

            /* add (CMIR strengthened) local cut if efficacious (cut is not local) */
            SCIP_CALL( addLocalCut(scip, sepa, sepadata, sol, knapvars, nknapvars, liftcoefs, liftrhs, FALSE, cutact, cmircoefs, cmirinds, nlc, infeasible, &success) );

            if ( *infeasible )
               break;
         }
      }
   }

   SCIPfreeBufferArrayNull(scip, &cmirinds);
   SCIPfreeBufferArrayNull(scip, &cmircoefs);

   SCIPfreeBufferArrayNull(scip, &solvals);
   SCIPfreeBufferArrayNull(scip, &xrelint);
   SCIPfreeBufferArrayNull(scip, &x);
   SCIPfreeBufferArrayNull(scip, &weights);
   SCIPfreeBufferArrayNull(scip, &varsF1);
   SCIPfreeBufferArrayNull(scip, &varsF0);
   SCIPfreeBufferArrayNull(scip, &redvars);
   SCIPfreeBufferArrayNull(scip, &vars);
   SCIPfreeBufferArrayNull(scip, &cutcoefs);
   SCIPfreeBufferArrayNull(scip, &liftcoefs);

   return SCIP_OKAY;
}


/** separate local cuts without lifting
 *
 *  In this function we reduce each knapsack by removing variables that are *fixed*. Cuts for the reduced knapsack are
 *  then locally valid. We do not perform lifting here, but mark the rows as locally valid only.
 */
static
SCIP_RETCODE separateLocalCutsKnapsack(
   SCIP*                 scip,               /**< scip pointer */
   SCIP_SEPA*            sepa,               /**< separator */
   SCIP_SEPADATA*        sepadata,           /**< data of separator */
   SCIP_SOL*             sol,                /**< SCIP solution to be separated (NULL if current LP-solution) */
   int*                  nlc,                /**< number of generated local cuts */
   SCIP_Bool*            infeasible          /**< pointer to store whether infeasibility was detected */
   )
{
   SCIP_CONS** conss;
   SCIP_Real* xrelint = NULL;
   SCIP_Real* x = NULL;
   SCIP_VAR** vars = NULL;
   SCIP_Longint* weights = NULL;
   SCIP_Real* cmircoefs = NULL;
   int* cmirinds = NULL;
   SCIP_Real* cutcoefs;
   SCIP_Real cutrhs = SCIP_INVALID;
   int nvars;
   int nconss;
   int c;

   assert( scip != NULL );
   assert( sepa != NULL );
   assert( sepadata != NULL );
   assert( nlc != NULL );
   assert( infeasible != NULL );

   *nlc = 0;
   *infeasible = FALSE;

   nconss = SCIPconshdlrGetNConss(sepadata->knapsackconshdlr);
   conss = SCIPconshdlrGetConss(sepadata->knapsackconshdlr);
   nvars = SCIPgetNVars(scip);
   if ( nconss == 0 || nvars == 0 )
      return SCIP_OKAY;
   assert( nconss > 0 );
   assert( nvars > 0 );

   SCIP_CALL( SCIPallocBufferArray(scip, &vars, nvars) );
   SCIP_CALL( SCIPallocBufferArray(scip, &weights, nvars) );
   SCIP_CALL( SCIPallocBufferArray(scip, &x, nvars) );
   SCIP_CALL( SCIPallocBufferArray(scip, &xrelint, nvars) );
   SCIP_CALL( SCIPallocBufferArray(scip, &cutcoefs, nvars) );

   if ( sepadata->usecmir )
   {
      SCIP_CALL( SCIPallocBufferArray(scip, &cmircoefs, nvars) );
      SCIP_CALL( SCIPallocBufferArray(scip, &cmirinds, nvars) );
   }

   /* loop through knapsack constraints */
   for (c = 0; c < nconss && ! SCIPisStopped(scip); ++c)
   {
      SCIP_CONS* cons;
      SCIP_VAR** knapvars;
      SCIP_Longint* knapweights;
      SCIP_Longint redcapacity;
      SCIP_Longint weightsum = 0;
      SCIP_Bool somefractional = FALSE;
      SCIP_Bool cutislocal = FALSE;
      SCIP_Bool success;
      int nknapvars;
      int dim = 0;
      int i;

      cons = conss[c];
      assert( cons != NULL );

      if ( sepadata->onlycheckconss && ! SCIPconsIsChecked(cons) )
         continue;

      nknapvars = SCIPgetNVarsKnapsack(scip, cons);
      knapvars = SCIPgetVarsKnapsack(scip, cons);
      knapweights = SCIPgetWeightsKnapsack(scip, cons);

      /* compute reduced capacity */
      redcapacity = SCIPgetCapacityKnapsack(scip, cons);
      for (i = 0; i < nknapvars; ++i)
      {
         assert( SCIPvarIsBinary(knapvars[i]) );
         assert( SCIPisIntegral(scip, SCIPvarGetLbLocal(knapvars[i])) );
         assert( SCIPisIntegral(scip, SCIPvarGetUbLocal(knapvars[i])) );

         /* if the variable is fixed to 1, reduce the capacity */
         if ( SCIPvarGetLbLocal(knapvars[i]) > 0.5 )
         {
            redcapacity -= knapweights[i];
            cutislocal = TRUE;  /* cut is local if there are variables fixed to 1 */
         }
      }

      /* if reduced problem is empty or trival, skip separation */
      if ( redcapacity <= 0 )
         continue;

      /* collect reduced knapsack */
      for (i = 0; i < nknapvars; ++i)
      {
         /* if variable is not fixed */
         if ( SCIPvarGetLbLocal(knapvars[i]) < 0.5 && SCIPvarGetUbLocal(knapvars[i]) > 0.5 )
         {
            if ( knapweights[i] <= redcapacity )
            {
               x[dim] = SCIPgetSolVal(scip, sol, knapvars[i]);
               if ( ! SCIPisFeasZero(scip, x[dim]) )
               {
                  /* keep variables not having value 0 - note that they do not have to be lifted to be valid */
                  vars[dim] = knapvars[i];
                  weights[dim] = knapweights[i];
                  weightsum += knapweights[i];
                  if ( ! SCIPisFeasIntegral(scip, x[dim]) )
                     somefractional = TRUE;
                  ++dim;
               }
            }
         }
      }

      /* skip all integral knapsacks */
      if ( ! somefractional )
         continue;

      /* skip trivially feasible knapsacks */
      if ( weightsum <= redcapacity )
         continue;

      /* Skip separation in trivial cases:
       * 1) there are no variables left (dim = 0).
       * 2) dim = 1, then the problem is either trivially valid or we have infeasible items (this case is treated above)
       */
      if ( dim <= 1 )
         continue;
      assert( dim > 0 );

      SCIPdebugMsg(scip, "Separating local cuts for knapsack constraint <%s> (dim = %d) ...\n", SCIPconsGetName(cons), dim);

      /* set up relative interior point */
      for (i = 0; i < dim; ++i)
         xrelint[i] = 1.0 / ((SCIP_Real)dim * redcapacity + 1.0);

#ifdef USE_FW_SEPA
      {
         FWTuple* fwtpl;
         uint64_t max_fw_iters = 1000;

         /* increase number of iterations in root node */
         if ( SCIPgetDepth(scip) == 0 )
            max_fw_iters = 10000;

         /* SCIP_CALL( SCIPgetIntParam(scip, "separating/" SEPA_NAME "/maxfwiters", &max_fw_iters) ); */

         fwtpl = scalar_afw_knap_gen_cut(scip, dim, weights, redcapacity, x, xrelint, sepadata->oracleclock, sepadata->fwclock, 1e-6, max_fw_iters);

         // record statistics
         sepadata->ntotaloraclecalls += fwtpl->noraclecalls;
         sepadata->ntotalfwiters += fwtpl->niters;
         sepadata->cummulativeredcapacity += redcapacity;
         sepadata->cummulativereddim += dim;
         sepadata->nreddimcalls++;

         if ( redcapacity > sepadata->maxredcapacity )
            sepadata->maxredcapacity = redcapacity;

         if ( dim > sepadata->maxreddim )
            sepadata->maxreddim = dim;

         switch ( fwtpl->exitcode )
         {
         case exit_max_iters: ++sepadata->nexit_max_iter; break;
         case exit_early: ++sepadata->nexit_early; break;
         case exit_zero_gradient: ++sepadata->nexit_zero_gradient; break;
         case exit_term_check: ++sepadata->nexit_term_check; break;
         case exit_primal_gap: ++sepadata->nexit_primal_gap; break;
         case exit_not_run:
         default: abort();
         }

         if ( fwtpl->success )
         {
            success = TRUE;
            for (i = 0; i < dim; ++i)
               cutcoefs[i] = fwtpl->a[i];
            cutrhs = fwtpl->b;
         }
         else
            success = FALSE;

         free(fwtpl->a);
         free(fwtpl);
      }
#else
      {
         LC_DATA data;
         SCIP_Real usedtime;

         data.nvars = dim;
         data.vars = vars;
         data.weights = weights;
         data.capacity = redcapacity;

         /* call separation */
         SCIP_CALL( LC_generateCut(scip, data.nvars, x, xrelint, data.vars, OracleKnapsack, &data, 0, NULL, NULL,
               0, NULL, NULL, NULL, NULL, NULL, FALSE, sepadata->tilting, FALSE, sepadata->nicecuts,
               cutcoefs, &cutrhs, &success, NULL, &usedtime) );
      }
#endif

      if ( success )
      {
         assert( cutrhs != SCIP_INVALID );

         /* add (CMIR strengthened) local cut if efficacious (use cutrhs as activity, because we have not computed it) */
         SCIP_CALL( addLocalCut(scip, sepa, sepadata, sol, vars, dim, cutcoefs, cutrhs, cutislocal, cutrhs, cmircoefs, cmirinds, nlc, infeasible, &success) );

         if ( *infeasible )
            break;
      }
   }
   SCIPfreeBufferArrayNull(scip, &cmirinds);
   SCIPfreeBufferArrayNull(scip, &cmircoefs);

   SCIPfreeBufferArrayNull(scip, &cutcoefs);
   SCIPfreeBufferArrayNull(scip, &xrelint);
   SCIPfreeBufferArrayNull(scip, &x);
   SCIPfreeBufferArrayNull(scip, &weights);
   SCIPfreeBufferArrayNull(scip, &vars);

   return SCIP_OKAY;
}




/*
 * Callback methods of separator
 */

/** destructor of separator to free user data (called when SCIP is exiting) */
static
SCIP_DECL_SEPAFREE(sepaFreeLC)
{
   SCIP_SEPADATA* sepadata;

   assert( scip != NULL );
   assert( sepa != NULL );
   assert( strcmp(SCIPsepaGetName(sepa), SEPA_NAME) == 0 );

   sepadata = SCIPsepaGetData(sepa);
   assert( sepadata != 0 );
   assert( sepadata->pool == NULL );

   SCIPinfoMessage(scip, NULL, "\n");
   if ( sepadata->npoolcuts > 0 )
   {
      SCIPinfoMessage(scip, NULL, "Separator <%s> local cut pool statistics:\n", SCIPsepaGetName(sepa));
      SCIPinfoMessage(scip, NULL, "Number of stored cuts: %7d\n", sepadata->npoolcuts);
      SCIPinfoMessage(scip, NULL, "Pool separating time:  %7.2f\n", sepadata->pooltime);
      SCIPinfoMessage(scip, NULL, "# of pool sep. calls:  %7" SCIP_LONGINT_FORMAT "\n", sepadata->npoolcalls);
      SCIPinfoMessage(scip, NULL, "# of sep. pool cuts:   %7" SCIP_LONGINT_FORMAT "\n", sepadata->npoolcutsfound);
      SCIPinfoMessage(scip, NULL, "# of added pool cuts:  %7" SCIP_LONGINT_FORMAT "\n", sepadata->npoolcutsadded);
      SCIPinfoMessage(scip, NULL, "maximal # of cuts:     %7" SCIP_LONGINT_FORMAT "\n", sepadata->nmaxpoolcuts);
      SCIPinfoMessage(scip, NULL, "\n");
   }

   if ( sepadata->nsepacuts > 0 )
   {
      SCIPinfoMessage(scip, NULL, "Local cuts separation calls:\t %7d\n", sepadata->ntotalsepacalls);
      SCIPinfoMessage(scip, NULL, "Local cuts separated:\t\t %7" SCIP_LONGINT_FORMAT "\n", sepadata->nsepacuts);
      SCIPinfoMessage(scip, NULL, "Local cuts Applied:\t\t %7" SCIP_LONGINT_FORMAT "\n", SCIPsepaGetNCutsApplied(sepa));
      if ( sepadata->usecmir > 0 )
         SCIPinfoMessage(scip, NULL, "Generated CMIR cuts:\t\t %7d\n", sepadata->ngencmir);
      SCIPinfoMessage(scip, NULL, "\n");
   }

#ifdef USE_FW_SEPA
   SCIPinfoMessage(scip, NULL, "Total time for oracle calls:\t %6.2f\n", LC_clockGetTime(sepadata->oracleclock));
   SCIPinfoMessage(scip, NULL, "Total time for FW calls:\t %6.2f\n", LC_clockGetTime(sepadata->fwclock));
   SCIPinfoMessage(scip, NULL, "Total number of oracle calls:\t %6" SCIP_LONGINT_FORMAT "\n", sepadata->ntotaloraclecalls);
   SCIPinfoMessage(scip, NULL, "Total number of FW iterations:\t %6" SCIP_LONGINT_FORMAT "\n", sepadata->ntotalfwiters);
   SCIPinfoMessage(scip, NULL, "\n");
   SCIPinfoMessage(scip, NULL, "Total number of FW runs reaching maximal number iterations:\t %6d\n", sepadata->nexit_max_iter);
   SCIPinfoMessage(scip, NULL, "Total number of FW runs with early termination:\t\t\t %6d\n", sepadata->nexit_early);
   SCIPinfoMessage(scip, NULL, "Total number of FW runs with zero gradient termination:\t\t %6d\n", sepadata->nexit_zero_gradient);
   SCIPinfoMessage(scip, NULL, "Total number of FW runs with termination check:\t\t\t %6d\n", sepadata->nexit_term_check);
   SCIPinfoMessage(scip, NULL, "Total number of FW runs reaching primal gap:\t\t\t %6d\n", sepadata->nexit_primal_gap);
   SCIPinfoMessage(scip, NULL, "\n");

   LC_clockFree(&sepadata->oracleclock);
   LC_clockFree(&sepadata->fwclock);
   if ( sepadata->ntotaloraclecalls > 0 )
   {
      SCIPinfoMessage(scip, NULL, "average reduced problem size:\t %6.2f\n", (double) sepadata->cummulativereddim / sepadata->nreddimcalls);
      SCIPinfoMessage(scip, NULL, "average reduced capacity:\t %6.2f\n", (double) sepadata->cummulativeredcapacity / sepadata->nreddimcalls);
      SCIPinfoMessage(scip, NULL, "Max reduced problem size:\t %6" SCIP_LONGINT_FORMAT "\n", sepadata->maxreddim);
      SCIPinfoMessage(scip, NULL, "Max reduced problem capacity:\t %6" SCIP_LONGINT_FORMAT "\n", sepadata->maxredcapacity);
   }
#endif

   SCIPdebugMsg(scip, "Freeing data of local separator <%s>\n", SCIPsepaGetName(sepa));

   SCIPfreeMemory(scip, &sepadata);

   return SCIP_OKAY;
}

/** initialization method of separator (called after problem was transformed) */
static
SCIP_DECL_SEPAINIT(sepaInitLC)
{
#ifdef USE_FW_SEPA
   SCIP_SEPADATA* sepadata;
   int clocktype;

   assert( scip != NULL );
   assert( sepa != NULL );
   assert( strcmp(SCIPsepaGetName(sepa), SEPA_NAME) == 0 );

   sepadata = SCIPsepaGetData(sepa);
   assert( sepadata != NULL );
   assert( sepadata->oracleclock == NULL );
   assert( sepadata->fwclock == NULL );

   sepadata->oracleclock = LC_clockCreate();
   sepadata->fwclock = LC_clockCreate();

   SCIP_CALL( SCIPgetIntParam(scip, "timing/clocktype", &clocktype) );

   /* wall clock time is the default */
   if ( clocktype == SCIP_CLOCKTYPE_CPU )
   {
      LC_clockSetType(sepadata->oracleclock, LC_CLOCKTYPE_CPU);
      LC_clockSetType(sepadata->fwclock, LC_CLOCKTYPE_CPU);
   }
#endif

   return SCIP_OKAY;
}


/** solving process initialization method of separator (called when branch and bound process is about to begin) */
static
SCIP_DECL_SEPAINITSOL(sepaInitsolLC)
{
   SCIP_SEPADATA* sepadata;
   int cutagelimit;

   assert( scip != NULL );
   assert( sepa != NULL );
   assert( strcmp(SCIPsepaGetName(sepa), SEPA_NAME) == 0 );

   sepadata = SCIPsepaGetData(sepa);
   assert( sepadata != 0 );
   assert( sepadata->pool == NULL );

   SCIPdebugMsg(scip, "creating cutpool for local cut separator <%s>\n", SCIPsepaGetName(sepa));

   SCIP_CALL( SCIPgetIntParam(scip, "separating/cutagelimit", &cutagelimit) );

   /* create cutpool (-1 = no age limit) */
   SCIP_CALL( SCIPcreateCutpool(scip, &(sepadata->pool), cutagelimit) );

   /* find knapsack constraint handler */
   sepadata->knapsackconshdlr = SCIPfindConshdlr(scip, "knapsack");
   if ( sepadata->knapsackconshdlr == NULL )
   {
      SCIPwarningMessage(scip, "Could not find knapsack constraint handler. Local cuts cannot be separated.\n");
   }

   return SCIP_OKAY;
}


/** solving process deinitialization method of separator (called before branch and bound process data is freed) */
static
SCIP_DECL_SEPAEXITSOL(sepaExitsolLC)
{
   SCIP_SEPADATA* sepadata;
   SCIP_CUTPOOL* pool;

   assert( scip != NULL );
   assert( sepa != NULL );
   assert( strcmp(SCIPsepaGetName(sepa), SEPA_NAME) == 0 );

   sepadata = SCIPsepaGetData(sepa);
   assert( sepadata != 0 );
   assert( sepadata->pool != NULL );

   pool = sepadata->pool;
   sepadata->npoolcuts += SCIPcutpoolGetNCuts(pool);
   sepadata->pooltime += SCIPcutpoolGetTime(pool);
   sepadata->npoolcalls += SCIPcutpoolGetNCalls(pool);
   sepadata->npoolcutsfound += SCIPcutpoolGetNCutsFound(pool);
   sepadata->npoolcutsadded += SCIPcutpoolGetNCutsAdded(pool);
   sepadata->nmaxpoolcuts += SCIPcutpoolGetMaxNCuts(pool);

   SCIPdebugMsg(scip, "deleting pool of separator <%s>\n", SCIPsepaGetName(sepa));
   SCIP_CALL( SCIPfreeCutpool(scip, &(sepadata->pool)) );

   return SCIP_OKAY;
}


/** LP solution separation method of separator */
static
SCIP_DECL_SEPAEXECLP(sepaExeclpLC)
{  /*lint --e{715}*/
   SCIP_SEPADATA* sepadata;
   SCIP_CUTPOOL* pool;
   SCIP_Bool infeasible = FALSE;
   int ncutsfound;
   int ncuts = 0;
   int ngen = 0;

   assert( scip != NULL );
   assert( sepa != NULL );
   assert( strcmp(SCIPsepaGetName(sepa), SEPA_NAME) == 0 );
   assert( result != NULL );

   sepadata = SCIPsepaGetData(sepa);
   assert( sepadata != 0 );

   *result = SCIP_DIDNOTRUN;

   if ( sepadata->knapsackconshdlr == NULL )
      return SCIP_OKAY;

   /* double check to avoid running in sub-SCIPs */
   if ( SCIPgetSubscipDepth(scip) > 0 )
      return SCIP_OKAY;

   SCIPdebugMsg(scip, "separating local cuts ...\n");

   /* if current node has changed reset counter */
   if ( SCIPnodeGetNumber(SCIPgetCurrentNode(scip)) != sepadata->currentnode )
   {
      sepadata->currentnode = SCIPnodeGetNumber(SCIPgetCurrentNode(scip));
      sepadata->nrounds = 0;
   }

   /* get depth of tree and number of cuts found in this separation round */
   ncutsfound = SCIPgetNCutsFoundRound(scip);
   SCIPdebugMsg(scip, "Cuts previously found: %d\n", ncutsfound);

   /* separate only if we have not reached the maximal number of rounds in this node */
   if ( (depth == 0 && (sepadata->maxroundsroot == -1 || sepadata->nrounds < sepadata->maxroundsroot)) ||
        (depth > 0  && (sepadata->maxrounds == -1 || sepadata->nrounds < sepadata->maxrounds)) )
   {
      /* only separate if fewer cuts than the given thresholds have been found in this separation round */
      if ( (depth == 0 && ncutsfound <= sepadata->thresholdroot ) ||
  	   (depth > 0 && ncutsfound <= sepadata->threshold ) )
      {
         SCIP_Longint noldcutspool;
         SCIP_Longint nfoundpool;

         *result = SCIP_DIDNOTFIND;

	 /* separate local cuts from the pool */
         pool = sepadata->pool;
         noldcutspool = SCIPcutpoolGetNCutsFound(pool);
         SCIP_CALL( SCIPseparateCutpool(scip, pool, result) );
         nfoundpool = SCIPcutpoolGetNCutsFound(pool) - noldcutspool;
         SCIPdebugMsg(scip, "separated %"SCIP_LONGINT_FORMAT" local cuts from the pool [%d]\n", nfoundpool, SCIPcutpoolGetNCuts(pool));
         ngen += (int) nfoundpool;
         sepadata->nsepacuts += nfoundpool;
         ncutsfound += nfoundpool;

	 /* separate local cuts - only if the number of cuts already found is not above the given thresholds */
	 if ( (depth == 0 && ncutsfound <= sepadata->thresholdroot) ||
	      (depth > 0 && ncutsfound <= sepadata->threshold) )
	 {
            ++sepadata->nrounds;
            ++sepadata->ntotalsepacalls;

            if ( sepadata->lifting )
            {
               SCIP_CALL( separateLocalCutsKnapsackLifting(scip, sepa, sepadata, NULL, &ncuts, &infeasible) );
            }
            else
            {
               SCIP_CALL( separateLocalCutsKnapsack(scip, sepa, sepadata, NULL, &ncuts, &infeasible) );
            }

	    SCIPdebugMsg(scip, "separated local cuts: %d  (round: %d).\n", ncuts, sepadata->nrounds);
	    sepadata->nsepacuts += ncuts;
	    ngen += ncuts;
	 }
      }
      else
         *result = SCIP_DELAYED;
   }
   if ( infeasible )
      *result = SCIP_CUTOFF;
   else if ( ngen > 0 )
      *result = SCIP_SEPARATED;

   return SCIP_OKAY;
}


/** arbitrary primal solution separation method of separator */
static
SCIP_DECL_SEPAEXECSOL(sepaExecsolLC)
{  /*lint --e{715}*/
   assert( result != NULL );

   *result = SCIP_DIDNOTRUN;

   return SCIP_OKAY;
}


/** creates the open position Local Cut separator and includes it in SCIP */
SCIP_RETCODE SCIPincludeSepaLocalcuts(
   SCIP*                 scip           /**< SCIP data structure */
   )
{
   SCIP_SEPADATA* sepadata;
   SCIP_SEPA* sepa;

   /* create separator data */
   SCIP_CALL( SCIPallocMemory(scip, &sepadata) );
   sepadata->pool = NULL;
   sepadata->nsepacuts = 0;
   sepadata->currentnode = -1;
   sepadata->knapsackconshdlr = NULL;
   sepadata->nrounds = 0;
   sepadata->ntotalgen = 0;
   sepadata->nexit_max_iter = 0;
   sepadata->nexit_early = 0;
   sepadata->nexit_zero_gradient = 0;
   sepadata->nexit_term_check = 0;
   sepadata->nexit_primal_gap = 0;

   sepadata->cummulativereddim = 0LL;
   sepadata->cummulativeredcapacity = 0LL;
   sepadata->maxreddim = 0LL;
   sepadata->maxredcapacity = 0LL;
   sepadata->nreddimcalls = 0LL;
   sepadata->ntotalsepacalls = 0;
   sepadata->ntotalfwiters = 0LL;
   sepadata->ngencmir = 0;
   sepadata->ntotaloraclecalls = 0LL;
#ifdef USE_FW_SEPA
   sepadata->oracleclock = NULL;
   sepadata->fwclock = NULL;
#endif
   sepadata->npoolcuts = 0;
   sepadata->pooltime = 0.0;
   sepadata->npoolcalls = 0;
   sepadata->npoolcutsfound = 0;
   sepadata->npoolcutsadded = 0;
   sepadata->nmaxpoolcuts = 0;

   /* include separator */
   SCIP_CALL( SCIPincludeSepaBasic(scip, &sepa, SEPA_NAME, SEPA_DESC, SEPA_PRIORITY, SEPA_FREQ, SEPA_MAXBOUNDDIST, SEPA_USESUBSCIP, SEPA_DELAY,
	 sepaExeclpLC, sepaExecsolLC, sepadata) );

   SCIP_CALL( SCIPsetSepaFree(scip, sepa, sepaFreeLC) );
   SCIP_CALL( SCIPsetSepaInit(scip, sepa, sepaInitLC) );
   SCIP_CALL( SCIPsetSepaInitsol(scip, sepa, sepaInitsolLC) );
   SCIP_CALL( SCIPsetSepaExitsol(scip, sepa, sepaExitsolLC) );

   /* add separator parameters */
   SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/threshold", "threshold when local cuts are separated",
	 &sepadata->threshold, TRUE, DEFAULT_THRESHOLD, 0, 100000, NULL, NULL) );

   SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/thresholdroot", "threshold when local cuts are separated (root)",
	 &sepadata->thresholdroot, TRUE, DEFAULT_THRESHOLDROOT, 0, 100000, NULL, NULL) );

   SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/maxrounds", "max. number of local cuts rounds per node",
	 &sepadata->maxrounds, TRUE, DEFAULT_MAXROUNDS, -1, 100000, NULL, NULL) );

   SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/maxroundsroot", "max. number of local cuts rounds in the root",
	 &sepadata->maxroundsroot, TRUE, DEFAULT_MAXROUNDSROOT, -1, 100000, NULL, NULL) );

   SCIP_CALL( SCIPaddBoolParam(scip, "separating/" SEPA_NAME "/nicecuts", "Try to make a nice cut by finding a rational approximation?",
         &sepadata->nicecuts, TRUE, DEFAULT_NICECUTS, NULL, NULL) );

   SCIP_CALL( SCIPaddBoolParam(scip, "separating/" SEPA_NAME "/tilting", "Try tilting to strengthen cuts?",
         &sepadata->tilting, TRUE, DEFAULT_TILTING, NULL, NULL) );

   SCIP_CALL( SCIPaddBoolParam(scip, "separating/" SEPA_NAME "/lifting", "Use lifting?",
         &sepadata->lifting, TRUE, DEFAULT_LIFTING, NULL, NULL) );

   SCIP_CALL( SCIPaddBoolParam(scip, "separating/" SEPA_NAME "/uplifting", "Perfrom uplifting (if lifting is true)?",
         &sepadata->uplifting, TRUE, DEFAULT_UPLIFTING, NULL, NULL) );

   SCIP_CALL( SCIPaddBoolParam(scip, "separating/" SEPA_NAME "/usecmir", "Additionally apply CMIR functions to generated cuts?",
         &sepadata->usecmir, TRUE, DEFAULT_USECMIR, NULL, NULL) );

   SCIP_CALL( SCIPaddBoolParam(scip, "separating/" SEPA_NAME "/onlycheckconss", "Only run separation on constraints that are checked?",
         &sepadata->onlycheckconss, TRUE, DEFAULT_ONLYCHECKCONSS, NULL, NULL) );

   return SCIP_OKAY;
}
