sepa_rlt.c

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*                  This file is part of the program and library             */
/*         SCIP --- Solving Constraint Integer Programs                      */
/*                                                                           */
/*  Copyright (c) 2002-2024 Zuse Institute Berlin (ZIB)                      */
/*                                                                           */
/*  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   sepa_rlt.c
 * @ingroup DEFPLUGINS_SEPA
 * @brief  separator for cuts generated by Reformulation-Linearization-Technique (RLT)
 * @author Fabian Wegscheider
 * @author Ksenia Bestuzheva
 *
 * @todo implement the possibility to add extra auxiliary variables for RLT (like in DOI 10.1080/10556788.2014.916287)
 * @todo add RLT cuts for the product of equality constraints
 * @todo implement dynamic addition of RLT cuts during branching (see DOI 10.1007/s10898-012-9874-7)
 * @todo use SCIPvarIsBinary instead of SCIPvarGetType() == SCIP_VARTYPE_BINARY ?
 * @todo parameter maxusedvars seems arbitrary (too large for small problems; too small for large problems); something more adaptive we can do? (e.g., all variables with priority >= x% of highest prio)
 */

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

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

#include "scip/sepa_rlt.h"
#include "scip/cons_nonlinear.h"
#include "scip/pub_lp.h"
#include "scip/expr_pow.h"
#include "scip/nlhdlr_bilinear.h"
#include "scip/cutsel_hybrid.h"


#define SEPA_NAME              "rlt"
#define SEPA_DESC              "reformulation-linearization-technique separator"
#define SEPA_PRIORITY                10 /**< priority for separation */
#define SEPA_FREQ                     0 /**< frequency for separating cuts; zero means to separate only in the root node */
#define SEPA_MAXBOUNDDIST           1.0 /**< maximal relative distance from the current node's dual bound to primal bound
                                         *   compared to best node's dual bound for applying separation.*/
#define SEPA_USESSUBSCIP          FALSE /**< does the separator use a secondary SCIP instance? */
#define SEPA_DELAY                FALSE /**< should separation method be delayed, if other separators found cuts? */

#define DEFAULT_MAXUNKNOWNTERMS       0 /**< maximum number of unknown bilinear terms a row can have to be used */
#define DEFAULT_MAXUSEDVARS         100 /**< maximum number of variables that will be used to compute rlt cuts */
#define DEFAULT_MAXNCUTS             -1 /**< maximum number of cuts that will be added per round */
#define DEFAULT_MAXROUNDS             1 /**< maximum number of separation rounds per node (-1: unlimited) */
#define DEFAULT_MAXROUNDSROOT        10 /**< maximum number of separation rounds in the root node (-1: unlimited) */
#define DEFAULT_ONLYEQROWS        FALSE /**< whether only equality rows should be used for rlt cuts */
#define DEFAULT_ONLYCONTROWS      FALSE /**< whether only continuous rows should be used for rlt cuts */
#define DEFAULT_ONLYORIGINAL       TRUE /**< whether only original variables and rows should be used for rlt cuts */
#define DEFAULT_USEINSUBSCIP      FALSE /**< whether the separator should also be used in sub-scips */
#define DEFAULT_USEPROJECTION     FALSE /**< whether the separator should first check projected rows */
#define DEFAULT_DETECTHIDDEN      FALSE /**< whether implicit products should be detected and separated by McCormick */
#define DEFAULT_HIDDENRLT         FALSE /**< whether RLT cuts should be added for hidden products */
#define DEFAULT_ADDTOPOOL          TRUE /**< whether globally valid RLT cuts are added to the global cut pool */

#define DEFAULT_GOODSCORE           1.0 /**< threshold for score of cut relative to best score to be considered good,
                                         *   so that less strict filtering is applied */
#define DEFAULT_BADSCORE            0.5 /**< threshold for score of cut relative to best score to be discarded */
#define DEFAULT_OBJPARALWEIGHT      0.0 /**< weight of objective parallelism in cut score calculation */
#define DEFAULT_EFFICACYWEIGHT      1.0 /**< weight of efficacy in cut score calculation */
#define DEFAULT_DIRCUTOFFDISTWEIGHT 0.0 /**< weight of directed cutoff distance in cut score calculation */
#define DEFAULT_GOODMAXPARALL       0.1 /**< maximum parallelism for good cuts */
#define DEFAULT_MAXPARALL           0.1 /**< maximum parallelism for non-good cuts */

#define MAXVARBOUND                1e+5 /**< maximum allowed variable bound for computing an RLT-cut */

/*
 * Data structures
 */

/** data object for pairs and triples of variables */
struct HashData
{
   SCIP_VAR*             vars[3];            /**< variables in the pair or triple, used for hash comparison */
   int                   nvars;              /**< number of variables */
   int                   nrows;              /**< number of rows */
   int                   firstrow;           /**< beginning of the corresponding row linked list */
};
typedef struct HashData HASHDATA;

/** data structure representing an array of variables together with number of elements and size;
 *  used for storing variables that are in some sense adjacent to a given variable
 */
struct AdjacentVarData
{
   SCIP_VAR**            adjacentvars;       /**< adjacent vars */
   int                   nadjacentvars;      /**< number of vars in adjacentvars */
   int                   sadjacentvars;      /**< size of adjacentvars */
};
typedef struct AdjacentVarData ADJACENTVARDATA;

/** separator data */
struct SCIP_SepaData
{
   SCIP_CONSHDLR*        conshdlr;           /**< nonlinear constraint handler */
   SCIP_Bool             iscreated;          /**< indicates whether the sepadata has been initialized yet */
   SCIP_Bool             isinitialround;     /**< indicates that this is the first round and original rows are used */

   /* bilinear variables */
   SCIP_VAR**            varssorted;         /**< variables that occur in bilinear terms sorted by priority */
   SCIP_HASHMAP*         bilinvardatamap;    /**< maps each bilinear var to ADJACENTVARDATA containing vars appearing
                                                  together with it in bilinear products */
   int*                  varpriorities;      /**< priorities of variables */
   int                   nbilinvars;         /**< total number of variables occurring in bilinear terms */
   int                   sbilinvars;         /**< size of arrays for variables occurring in bilinear terms */

   /* information about bilinear terms */
   int*                  eqauxexpr;          /**< position of the auxexpr that is equal to the product (-1 if none) */
   int                   nbilinterms;        /**< total number of bilinear terms */

   /* parameters */
   int                   maxunknownterms;    /**< maximum number of unknown bilinear terms a row can have to be used (-1: unlimited) */
   int                   maxusedvars;        /**< maximum number of variables that will be used to compute rlt cuts (-1: unlimited) */
   int                   maxncuts;           /**< maximum number of cuts that will be added per round (-1: unlimited) */
   int                   maxrounds;          /**< maximum number of separation rounds per node (-1: unlimited) */
   int                   maxroundsroot;      /**< maximum number of separation rounds in the root node (-1: unlimited) */
   SCIP_Bool             onlyeqrows;         /**< whether only equality rows should be used for rlt cuts */
   SCIP_Bool             onlycontrows;       /**< whether only continuous rows should be used for rlt cuts */
   SCIP_Bool             onlyoriginal;       /**< whether only original rows and variables should be used for rlt cuts */
   SCIP_Bool             useinsubscip;       /**< whether the separator should also be used in sub-scips */
   SCIP_Bool             useprojection;      /**< whether the separator should first check projected rows */
   SCIP_Bool             detecthidden;       /**< whether implicit products should be detected and separated by McCormick */
   SCIP_Bool             hiddenrlt;          /**< whether RLT cuts should be added for hidden products */
   SCIP_Bool             addtopool;          /**< whether globally valid RLT cuts are added to the global cut pool */

   /* cut selection parameters */
   SCIP_Real             goodscore;          /**< threshold for score of cut relative to best score to be considered good,
                                              *   so that less strict filtering is applied */
   SCIP_Real             badscore;           /**< threshold for score of cut relative to best score to be discarded */
   SCIP_Real             objparalweight;     /**< weight of objective parallelism in cut score calculation */
   SCIP_Real             efficacyweight;     /**< weight of efficacy in cut score calculation */
   SCIP_Real             dircutoffdistweight;/**< weight of directed cutoff distance in cut score calculation */
   SCIP_Real             goodmaxparall;      /**< maximum parallelism for good cuts */
   SCIP_Real             maxparall;          /**< maximum parallelism for non-good cuts */
};

/* a simplified representation of an LP row */
struct RLT_SimpleRow
{
   const char*           name;               /**< name of the row */
   SCIP_Real*            coefs;              /**< coefficients */
   SCIP_VAR**            vars;               /**< variables */
   SCIP_Real             rhs;                /**< right hand side */
   SCIP_Real             lhs;                /**< left hand side */
   SCIP_Real             cst;                /**< constant */
   int                   nnonz;              /**< number of nonzeroes */
   int                   size;               /**< size of the coefs and vars arrays */
};
typedef struct RLT_SimpleRow RLT_SIMPLEROW;

/*
 * Local methods
 */

/** returns TRUE iff both keys are equal
 *
 * two variable pairs/triples are equal if the variables are equal
 */
static
SCIP_DECL_HASHKEYEQ(hashdataKeyEqConss)
{  /*lint --e{715}*/
   HASHDATA* hashdata1;
   HASHDATA* hashdata2;
   int v;

   hashdata1 = (HASHDATA*)key1;
   hashdata2 = (HASHDATA*)key2;

   /* check data structure */
   assert(hashdata1->nvars == hashdata2->nvars);
   assert(hashdata1->firstrow != -1 || hashdata2->firstrow != -1);

   for( v = hashdata1->nvars-1; v >= 0; --v )
   {
      /* tests if variables are equal */
      if( hashdata1->vars[v] != hashdata2->vars[v] )
         return FALSE;

      assert(SCIPvarCompare(hashdata1->vars[v], hashdata2->vars[v]) == 0);
   }

   /* if two hashdata objects have the same variables, then either one of them doesn't have a row list yet
    * (firstrow == -1) or they both point to the same row list
    */
   assert(hashdata1->firstrow == -1 || hashdata2->firstrow == -1 || hashdata1->firstrow == hashdata2->firstrow);

   return TRUE;
}

/** returns the hash value of the key */
static
SCIP_DECL_HASHKEYVAL(hashdataKeyValConss)
{  /*lint --e{715}*/
   HASHDATA* hashdata;
   int minidx;
   int mididx;
   int maxidx;
   int idx[3];

   hashdata = (HASHDATA*)key;
   assert(hashdata != NULL);
   assert(hashdata->nvars == 3 || hashdata->nvars == 2);

   idx[0] = SCIPvarGetIndex(hashdata->vars[0]);
   idx[1] = SCIPvarGetIndex(hashdata->vars[1]);
   idx[2] = SCIPvarGetIndex(hashdata->vars[hashdata->nvars - 1]);

   minidx = MIN3(idx[0], idx[1], idx[2]);
   maxidx = MAX3(idx[0], idx[1], idx[2]);
   if( idx[0] == maxidx )
      mididx = MAX(idx[1], idx[2]);
   else
      mididx = MAX(idx[0], MIN(idx[1], idx[2]));

   /* vars should already be sorted by index */
   assert(minidx <= mididx && mididx <= maxidx);

   return SCIPhashFour(hashdata->nvars, minidx, mididx, maxidx);
}

/** store a pair of adjacent variables */
static
SCIP_RETCODE addAdjacentVars(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_HASHMAP*         adjvarmap,          /**< hashmap mapping variables to their ADJACENTVARDATAs */
   SCIP_VAR**            vars                /**< variable pair to be stored */
   )
{
   int v1;
   int v2;
   int i;
   ADJACENTVARDATA* adjacentvardata;

   assert(adjvarmap != NULL);

   /* repeat for each variable of the new pair */
   for( v1 = 0; v1 < 2; ++v1 )
   {
      v2 = 1 - v1;

      /* look for data of the first variable */
      adjacentvardata = (ADJACENTVARDATA*) SCIPhashmapGetImage(adjvarmap, (void*)(size_t) SCIPvarGetIndex(vars[v1]));

      /* if the first variable has not been added to adjvarmap yet, add it here */
      if( adjacentvardata == NULL )
      {
         SCIP_CALL( SCIPallocClearBlockMemory(scip, &adjacentvardata) );
         SCIP_CALL( SCIPhashmapInsert(adjvarmap, (void*)(size_t) SCIPvarGetIndex(vars[v1]), adjacentvardata) );
      }

      assert(adjacentvardata != NULL);

      /* look for second variable in adjacentvars of the first variable */
      if( adjacentvardata->adjacentvars == NULL )
      {
         /* we don't know how many adjacent vars there will be - take a guess */
         SCIP_CALL( SCIPallocBlockMemoryArray(scip, &adjacentvardata->adjacentvars, 4) );
         adjacentvardata->adjacentvars[0] = vars[v2];
         ++adjacentvardata->nadjacentvars;
         adjacentvardata->sadjacentvars = 4;
      }
      else
      {
         SCIP_Bool found;
         int pos2;

         found = SCIPsortedvecFindPtr((void**) adjacentvardata->adjacentvars, SCIPvarComp, vars[v2],
               adjacentvardata->nadjacentvars, &pos2);

         /* add second var to adjacentvardata->adjacentvars, if not already added */
         if( !found )
         {
            /* ensure size of adjacentvardata->adjacentvars */
            SCIP_CALL( SCIPensureBlockMemoryArray(scip, &adjacentvardata->adjacentvars, &adjacentvardata->sadjacentvars,
                  adjacentvardata->nadjacentvars + 1) );

            /* insert second var at the correct position */
            for( i = adjacentvardata->nadjacentvars; i > pos2; --i )
            {
               adjacentvardata->adjacentvars[i] = adjacentvardata->adjacentvars[i-1];
            }
            adjacentvardata->adjacentvars[pos2] = vars[v2];
            ++adjacentvardata->nadjacentvars;
         }
      }

      /* if this is a self-adjacent var, only need to add the connection once */
      if( vars[v1] == vars[v2] )
         break;
   }

   return SCIP_OKAY;
}

/** returns the array of adjacent variables for a given variable */
static
SCIP_VAR** getAdjacentVars(
   SCIP_HASHMAP*         adjvarmap,          /**< hashmap mapping variables to their ADJACENTVARDATAs */
   SCIP_VAR*             var,                /**< variable */
   int*                  nadjacentvars       /**< buffer to store the number of variables in the returned array */
   )
{
   ADJACENTVARDATA* adjacentvardata;

   assert(adjvarmap != NULL);

   *nadjacentvars = 0;
   adjacentvardata = (ADJACENTVARDATA*) SCIPhashmapGetImage(adjvarmap, (void*)(size_t) SCIPvarGetIndex(var));

   if( adjacentvardata == NULL )
      return NULL;

   *nadjacentvars = adjacentvardata->nadjacentvars;

   return adjacentvardata->adjacentvars;
}

/** frees all ADJACENTVARDATAs stored in a hashmap */
static
void clearVarAdjacency(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_HASHMAP*         adjvarmap           /**< hashmap mapping variables to their ADJACENTVARDATAs */
   )
{
   int i;
   SCIP_HASHMAPENTRY* entry;
   ADJACENTVARDATA* adjacentvardata;

   assert(adjvarmap != NULL);

   for( i = 0; i < SCIPhashmapGetNEntries(adjvarmap); ++i )
   {
      entry = SCIPhashmapGetEntry(adjvarmap, i);

      if( entry == NULL )
         continue;

      adjacentvardata = (ADJACENTVARDATA*) SCIPhashmapEntryGetImage(entry);

      /* if adjacentvardata has been added to the hashmap, it can't be empty */
      assert(adjacentvardata->adjacentvars != NULL);

      SCIPfreeBlockMemoryArray(scip, &adjacentvardata->adjacentvars, adjacentvardata->sadjacentvars);
      SCIPfreeBlockMemory(scip, &adjacentvardata);
   }
}

/** free separator data */
static
SCIP_RETCODE freeSepaData(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_SEPADATA*        sepadata            /**< separation data */
   )
{  /*lint --e{715}*/
   int i;

   assert(sepadata->iscreated);

   if( sepadata->nbilinvars != 0 )
   {
      /* release bilinvars that were captured for rlt and free all related arrays */

      /* if there are bilinear vars, some of them must also participate in the same product */
      assert(sepadata->bilinvardatamap != NULL);

      clearVarAdjacency(scip, sepadata->bilinvardatamap);

      for( i = 0; i < sepadata->nbilinvars; ++i )
      {
         assert(sepadata->varssorted[i] != NULL);
         SCIP_CALL( SCIPreleaseVar(scip, &(sepadata->varssorted[i])) );
      }

      SCIPhashmapFree(&sepadata->bilinvardatamap);
      SCIPfreeBlockMemoryArray(scip, &sepadata->varssorted, sepadata->sbilinvars);
      SCIPfreeBlockMemoryArray(scip, &sepadata->varpriorities, sepadata->sbilinvars);
      sepadata->nbilinvars = 0;
      sepadata->sbilinvars = 0;
   }

   /* free the remaining array */
   if( sepadata->nbilinterms > 0 )
   {
      SCIPfreeBlockMemoryArray(scip, &sepadata->eqauxexpr, sepadata->nbilinterms);
   }

   sepadata->iscreated = FALSE;

   return SCIP_OKAY;
}

/** creates and returns rows of original linear constraints */
static
SCIP_RETCODE getOriginalRows(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_ROW***           rows,               /**< buffer to store the rows */
   int*                  nrows               /**< buffer to store the number of linear rows */
   )
{
   SCIP_CONS** conss;
   int nconss;
   int i;

   assert(rows != NULL);
   assert(nrows != NULL);

   conss = SCIPgetConss(scip);
   nconss = SCIPgetNConss(scip);
   *nrows = 0;

   SCIP_CALL( SCIPallocBufferArray(scip, rows, nconss) );

   for( i = 0; i < nconss; ++i )
   {
      SCIP_ROW *row;

      row = SCIPconsGetRow(scip, conss[i]);

      if( row != NULL )
      {
         (*rows)[*nrows] = row;
         ++*nrows;
      }
   }

   return SCIP_OKAY;
}

/** fills an array of rows suitable for RLT cut generation */
static
SCIP_RETCODE storeSuitableRows(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_SEPA*            sepa,               /**< separator */
   SCIP_SEPADATA*        sepadata,           /**< separator data */
   SCIP_ROW**            prob_rows,          /**< problem rows */
   SCIP_ROW**            rows,               /**< an array to be filled with suitable rows */
   int*                  nrows,              /**< buffer to store the number of suitable rows */
   SCIP_HASHMAP*         row_to_pos,         /**< hashmap linking row indices to positions in rows */
   SCIP_Bool             allowlocal          /**< are local rows allowed? */
   )
{
   int new_nrows;
   int r;
   int j;
   SCIP_Bool iseqrow;
   SCIP_COL** cols;
   SCIP_Bool iscontrow;

   new_nrows = 0;

   for( r = 0; r < *nrows; ++r )
   {
      iseqrow = SCIPisEQ(scip, SCIProwGetLhs(prob_rows[r]), SCIProwGetRhs(prob_rows[r]));

      /* if equality rows are requested, only those can be used */
      if( sepadata->onlyeqrows && !iseqrow )
         continue;

      /* if global cuts are requested, only globally valid rows can be used */
      if( !allowlocal && SCIProwIsLocal(prob_rows[r]) )
         continue;

      /* if continuous rows are requested, only those can be used */
      if( sepadata->onlycontrows )
      {
         cols = SCIProwGetCols(prob_rows[r]);
         iscontrow = TRUE;

         /* check row for integral variables */
         for( j = 0; j < SCIProwGetNNonz(prob_rows[r]); ++j )
         {
            if( SCIPcolIsIntegral(cols[j]) )
            {
               iscontrow = FALSE;
               break;
            }
         }

         if( !iscontrow )
            continue;
      }

      /* don't try to use rows that have been generated by the RLT separator */
      if( SCIProwGetOriginSepa(prob_rows[r]) == sepa )
         continue;

      /* if we are here, the row has passed all checks and should be added to rows */
      rows[new_nrows] = prob_rows[r];
      SCIP_CALL( SCIPhashmapSetImageInt(row_to_pos, (void*)(size_t)SCIProwGetIndex(prob_rows[r]), new_nrows) ); /*lint !e571 */
      ++new_nrows;
   }

   *nrows = new_nrows;

   return SCIP_OKAY;
}

/** make sure that the arrays in sepadata are large enough to store information on n variables */
static
SCIP_RETCODE ensureVarsSize(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_SEPADATA*        sepadata,           /**< separator data */
   int                   n                   /**< number of variables that we need to store */
   )
{
   int newsize;

   /* check whether array is large enough */
   if( n <= sepadata->sbilinvars )
      return SCIP_OKAY;

   /* compute new size */
   newsize = SCIPcalcMemGrowSize(scip, n);
   assert(n <= newsize);

   /* realloc arrays */
   SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &sepadata->varssorted, sepadata->sbilinvars, newsize) );
   SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &sepadata->varpriorities, sepadata->sbilinvars, newsize) );

   sepadata->sbilinvars = newsize;

   return SCIP_OKAY;
}

/** saves variables x and y to separator data and stores information about their connection
 *
 *  variables must be captured separately
 */
static
SCIP_RETCODE addProductVars(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_SEPADATA*        sepadata,           /**< separator data */
   SCIP_VAR*             x,                  /**< x variable */
   SCIP_VAR*             y,                  /**< y variable */
   SCIP_HASHMAP*         varmap,             /**< hashmap linking var index to position */
   int                   nlocks              /**< number of locks */
   )
{
   int xpos;
   int ypos;
   int xidx;
   int yidx;
   SCIP_VAR* vars[2];

   if( sepadata->bilinvardatamap == NULL )
   {
      int varmapsize;
      int nvars;

      /* the number of variables participating in bilinear products cannot exceed twice the number of bilinear terms;
       * however, if we detect hidden products, the number of terms is yet unknown, so use the number of variables
       */
      nvars = SCIPgetNVars(scip);
      varmapsize = sepadata->detecthidden ? nvars : MIN(nvars, sepadata->nbilinterms * 2);

      SCIP_CALL( SCIPhashmapCreate(&sepadata->bilinvardatamap, SCIPblkmem(scip), varmapsize) );
   }

   xidx = SCIPvarGetIndex(x);
   yidx = SCIPvarGetIndex(y);

   xpos = SCIPhashmapGetImageInt(varmap, (void*)(size_t) xidx); /*lint !e571 */

   if( xpos == INT_MAX )
   {
      /* add x to sepadata and initialise its priority */
      SCIP_CALL( SCIPhashmapInsertInt(varmap, (void*)(size_t) xidx, sepadata->nbilinvars) ); /*lint !e571*/
      SCIP_CALL( ensureVarsSize(scip, sepadata, sepadata->nbilinvars + 1) );
      sepadata->varssorted[sepadata->nbilinvars] = x;
      sepadata->varpriorities[sepadata->nbilinvars] = 0;
      xpos = sepadata->nbilinvars;
      ++sepadata->nbilinvars;
   }

   assert(xpos >= 0 && xpos < sepadata->nbilinvars );
   assert(xpos == SCIPhashmapGetImageInt(varmap, (void*)(size_t) xidx)); /*lint !e571 */

   /* add locks to priority of x */
   sepadata->varpriorities[xpos] += nlocks;

   if( xidx != yidx )
   {
      ypos = SCIPhashmapGetImageInt(varmap, (void*)(size_t) yidx); /*lint !e571 */

      if( ypos == INT_MAX )
      {
         /* add y to sepadata and initialise its priority */
         SCIP_CALL( SCIPhashmapInsertInt(varmap, (void*)(size_t) yidx, sepadata->nbilinvars) ); /*lint !e571*/
         SCIP_CALL( ensureVarsSize(scip, sepadata, sepadata->nbilinvars + 1) );
         sepadata->varssorted[sepadata->nbilinvars] = y;
         sepadata->varpriorities[sepadata->nbilinvars] = 0;
         ypos = sepadata->nbilinvars;
         ++sepadata->nbilinvars;
      }

      assert(ypos >= 0 && ypos < sepadata->nbilinvars);
      assert(ypos == SCIPhashmapGetImageInt(varmap, (void*)(size_t) yidx)); /*lint !e571 */

      /* add locks to priority of y */
      sepadata->varpriorities[ypos] += nlocks;
   }

   /* remember the connection between x and y */
   vars[0] = x;
   vars[1] = y;
   SCIP_CALL( addAdjacentVars(scip, sepadata->bilinvardatamap, vars) );

   return SCIP_OKAY;
}

/** extract a bilinear product from two linear relations, if possible
 *
 * First, the two given rows are brought to the form:
 * \f[
 *   a_1x + b_1w + c_1y \leq/\geq d_1,\\
 *   a_2x + b_2w + c_2y \leq/\geq d_2,
 * \f]
 * where \f$ a_1a_2 \leq 0 \f$ and the first implied relation is enabled when \f$ x = 1 \f$
 * and the second when \f$ x = 0 \f$, and \f$ b_1, b_2 > 0 \f$, the product relation can be written as:
 * \f[
 *   \frac{b_1b_2w + (b_2(a_1 - d_1) + b_1d_2)x + b_1c_2y - b_1d_2}{b_1c_2 - c_1b_2} \leq/\geq xy.
 * \f]
 * The inequality sign in the product relation is similar to that in the given linear relations if
 * \f$ b_1c_2 - c_1b_2 > 0 \f$ and opposite if \f$ b_1c_2 - c_1b_2 > 0 \f$.
 *
 * To obtain this formula, the given relations are first multiplied by scaling factors \f$ \alpha \f$
 * and \f$ \beta \f$, which is necessary in order for the solution to always exist, and written as
 * implications:
 * \f{align}{
 *   x = 1 & ~\Rightarrow~ \alpha b_1w + \alpha c_1y \leq/\geq \alpha(d_1 - a_1), \\
 *   x = 0 & ~\Rightarrow~ \beta b_2w + \beta c_2y \leq/\geq \beta d_2.
 * \f}
 * Then a linear system is solved which ensures that the coefficients of the two implications of the product
 * relation are equal to the corresponding coefficients in the linear relations.
 * If the product relation is written as:
 * \f[
 *   Ax + Bw + Cy + D \leq/\geq xy,
 * \f]
 * then the system is
 * \f[
 *   B = \alpha b_1, ~C - 1 = \alpha c_1, ~D+A = \alpha(a_1-d_1),\\
 *   B = \beta b_2, ~C = \beta c_2, ~D = -\beta d_2.
 * \f]
 */
static
SCIP_RETCODE extractProducts(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_SEPADATA*        sepadata,           /**< separator data */
   SCIP_VAR**            vars_xwy,           /**< 3 variables involved in the inequalities in the order x,w,y */
   SCIP_Real*            coefs1,             /**< coefficients of the first inequality (always implied, i.e. has x) */
   SCIP_Real*            coefs2,             /**< coefficients of the second inequality (can be unconditional) */
   SCIP_Real             d1,                 /**< side of the first inequality */
   SCIP_Real             d2,                 /**< side of the second inequality */
   SCIP_SIDETYPE         sidetype1,          /**< side type (lhs or rls) in the first inequality */
   SCIP_SIDETYPE         sidetype2,          /**< side type (lhs or rhs) in the second inequality */
   SCIP_HASHMAP*         varmap,             /**< variable map */
   SCIP_Bool             f                   /**< the first relation is an implication x == f */
   )
{
   SCIP_Real mult;

   /* coefficients and constant of the auxexpr */
   SCIP_Real A; /* coefficient of x */
   SCIP_Real B; /* coefficient of w */
   SCIP_Real C; /* coefficient of y */
   SCIP_Real D; /* constant */

   /* variables */
   SCIP_VAR* w;
   SCIP_VAR* x;
   SCIP_VAR* y;

   /* does auxexpr overestimate the product? */
   SCIP_Bool overestimate;

   /* coefficients in given relations: a for x, b for w, c for y; 1 and 2 for 1st and 2nd relation, respectively */
   SCIP_Real a1 = coefs1[0];
   SCIP_Real b1 = coefs1[1];
   SCIP_Real c1 = coefs1[2];
   SCIP_Real a2 = coefs2[0];
   SCIP_Real b2 = coefs2[1];
   SCIP_Real c2 = coefs2[2];

   x = vars_xwy[0];
   w = vars_xwy[1];
   y = vars_xwy[2];

   /* check given linear relations and decide if to continue */

   assert(SCIPvarGetType(x) == SCIP_VARTYPE_BINARY);  /* x must be binary */
   assert(a1 != 0.0); /* the first relation is always conditional */
   assert(b1 != 0.0 || b2 != 0.0); /* at least one w coefficient must be nonzero */

   SCIPdebugMsg(scip, "Extracting product from two implied relations:\n");
   SCIPdebugMsg(scip, "Relation 1: <%s> == %u => %g<%s> + %g<%s> %s %g\n", SCIPvarGetName(x), f, b1,
      SCIPvarGetName(w), c1, SCIPvarGetName(y), sidetype1 == SCIP_SIDETYPE_LEFT ? ">=" : "<=",
      f ? d1 - a1 : d1);
   SCIPdebugMsg(scip, "Relation 2: <%s> == %d => %g<%s> + %g<%s> %s %g\n", SCIPvarGetName(x), !f, b2,
      SCIPvarGetName(w), c2, SCIPvarGetName(y), sidetype2 == SCIP_SIDETYPE_LEFT ? ">=" : "<=",
      f ? d2 : d2 - a2);

   /* cannot use a global bound on x to detect a product */
   if( (b1 == 0.0 && c1 == 0.0) || (b2 == 0.0 && c2 == 0.0) )
      return SCIP_OKAY;

   /* cannot use a global bound on y to detect a non-redundant product relation */
   if( a2 == 0.0 && b2 == 0.0 ) /* only check the 2nd relation because the 1st at least has x */
   {
      SCIPdebugMsg(scip, "Ignoring a global bound on y\n");
      return SCIP_OKAY;
   }

   SCIPdebugMsg(scip, "binary var = <%s>, product of its coefs: %g\n", SCIPvarGetName(x), a1*a2);

   /* rewrite the linear relations in a standard form:
    * a1x + b1w + c1y <=/>= d1,
    * a2x + b2w + c2y <=/>= d2,
    * where b1 > 0, b2 > 0 and first implied relation is activated when x == 1
    */

   /* if needed, multiply the rows by -1 so that coefs of w are positive */
   if( b1 < 0 )
   {
      a1 *= -1.0;
      b1 *= -1.0;
      c1 *= -1.0;
      d1 *= -1.0;
      sidetype1 = sidetype1 == SCIP_SIDETYPE_LEFT ? SCIP_SIDETYPE_RIGHT : SCIP_SIDETYPE_LEFT;
   }
   if( b2 < 0 )
   {
      a2 *= -1.0;
      b2 *= -1.0;
      c2 *= -1.0;
      d2 *= -1.0;
      sidetype2 = sidetype2 == SCIP_SIDETYPE_LEFT ? SCIP_SIDETYPE_RIGHT : SCIP_SIDETYPE_LEFT;
   }

   /* the linear relations imply a product only if the inequality signs are similar */
   if( sidetype1 != sidetype2 )
      return SCIP_OKAY;

   /* when b1c2 = b2c1, the linear relations do not imply a product relation */
   if( SCIPisRelEQ(scip, b2*c1, c2*b1) )
   {
      SCIPdebugMsg(scip, "Ignoring a pair of linear relations because b1c2 = b2c1\n");
      return SCIP_OKAY;
   }

   if( !f )
   {
      /* swap the linear relations so that the relation implied by x == TRUE goes first */
      SCIPswapReals(&a1, &a2);
      SCIPswapReals(&b1, &b2);
      SCIPswapReals(&c1, &c2);
      SCIPswapReals(&d1, &d2);
   }

   /* all conditions satisfied, we can extract the product and write it as:
    * (1/(b1c2 - c1b2))*(b1b2w + (b2(a1 - d1) + b1d2)x + b1c2y - b1d2) >=/<= xy,
    * where the inequality sign in the product relation is similar to that in the given linear relations
    * if b1c2 - c1b2 > 0 and opposite if b1c2 - c1b2 > 0
    */

   /* compute the multiplier */
   mult = 1/(b1*c2 - c1*b2);

   /* determine the inequality sign; only check sidetype1 because sidetype2 is equal to it */
   overestimate = (sidetype1 == SCIP_SIDETYPE_LEFT && mult > 0.0) || (sidetype1 == SCIP_SIDETYPE_RIGHT && mult < 0.0);

   SCIPdebugMsg(scip, "found suitable implied rels (w,x,y): %g<%s> + %g<%s> + %g<%s> <= %g\n", a1,
      SCIPvarGetName(x), b1, SCIPvarGetName(w), c1, SCIPvarGetName(y), d1);
   SCIPdebugMsg(scip, "  and %g<%s> + %g<%s> + %g<%s> <= %g\n", a2, SCIPvarGetName(x),
      b2, SCIPvarGetName(w), c2, SCIPvarGetName(y), d2);

   /* compute the coefficients for x, w and y and the constant in auxexpr */
   A = (b2*a1 - d1*b2 + d2*b1)*mult;
   B = b1*b2*mult;
   C = b1*c2*mult;
   D = -b1*d2*mult;

   SCIPdebugMsg(scip, "product: <%s><%s> %s %g<%s> + %g<%s> + %g<%s> + %g\n", SCIPvarGetName(x), SCIPvarGetName(y),
      overestimate ? "<=" : ">=", A, SCIPvarGetName(x), B, SCIPvarGetName(w), C, SCIPvarGetName(y), D);

   SCIP_CALL( addProductVars(scip, sepadata, x, y, varmap, 1) );
   SCIP_CALL( SCIPinsertBilinearTermImplicitNonlinear(scip, sepadata->conshdlr, x, y, w, A, C, B, D, overestimate) );

   return SCIP_OKAY;
}

/** convert an implied bound: `binvar` = `binval`  &rArr;  `implvar` &le;/&ge; `implbnd` into a big-M constraint */
static
void implBndToBigM(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_VAR**            vars_xwy,           /**< variables in order x,w,y */
   int                   binvarpos,          /**< position of binvar in vars_xwy */
   int                   implvarpos,         /**< position of implvar in vars_xwy */
   SCIP_BOUNDTYPE        bndtype,            /**< type of implied bound */
   SCIP_Bool             binval,             /**< value of binvar which implies the bound */
   SCIP_Real             implbnd,            /**< value of the implied bound */
   SCIP_Real*            coefs,              /**< coefficients of the big-M constraint */
   SCIP_Real*            side                /**< side of the big-M constraint */
   )
{
   SCIP_VAR* implvar;
   SCIP_Real globbnd;

   assert(vars_xwy != NULL);
   assert(coefs != NULL);
   assert(side != NULL);
   assert(binvarpos != implvarpos);

   implvar = vars_xwy[implvarpos];
   globbnd = bndtype == SCIP_BOUNDTYPE_LOWER ? SCIPvarGetLbGlobal(implvar) : SCIPvarGetUbGlobal(implvar);

   /* Depending on the bound type and binval, there are four possibilities:
    * binvar == 1  =>  implvar >= implbnd   <=>   (implvar^l - implbnd)binvar + implvar >= implvar^l;
    * binvar == 0  =>  implvar >= implbnd   <=>   (implbnd - implvar^l)binvar + implvar >= implbnd;
    * binvar == 1  =>  implvar <= implbnd   <=>   (implvar^u - implbnd)binvar + implvar <= implvar^u;
    * binvar == 0  =>  implvar <= implbnd   <=>   (implbnd - implvar^u)binvar + implvar <= implbnd.
    */

   coefs[0] = 0.0;
   coefs[1] = 0.0;
   coefs[2] = 0.0;
   coefs[binvarpos] = binval ? globbnd - implbnd : implbnd - globbnd;
   coefs[implvarpos] = 1.0;
   *side = binval ? globbnd : implbnd;

   SCIPdebugMsg(scip, "Got an implied relation with binpos = %d, implpos = %d, implbnd = %g, "
                   "bnd type = %s, binval = %u, glbbnd = %g\n", binvarpos, implvarpos, implbnd,
                   bndtype == SCIP_BOUNDTYPE_LOWER ? "lower" : "upper", binval, globbnd);
   SCIPdebugMsg(scip, "Constructed big-M: %g*bvar + implvar %s %g\n", coefs[binvarpos],
                   bndtype == SCIP_BOUNDTYPE_LOWER ? ">=" : "<=", *side);
}

/** extract products from a relation given by coefs1, vars, side1 and sidetype1 and
 *  implied bounds of the form `binvar` = `!f` &rArr; `implvar` &ge;/&le; `implbnd`
 */
static
SCIP_RETCODE detectProductsImplbnd(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_SEPADATA*        sepadata,           /**< separator data */
   SCIP_Real*            coefs1,             /**< coefficients of the first linear relation */
   SCIP_VAR**            vars_xwy,           /**< variables in the order x, w, y */
   SCIP_Real             side1,              /**< side of the first relation */
   SCIP_SIDETYPE         sidetype1,          /**< is the left or right hand side given for the first relation? */
   int                   binvarpos,          /**< position of the indicator variable in the vars_xwy array */
   int                   implvarpos,         /**< position of the variable that is bounded */
   SCIP_HASHMAP*         varmap,             /**< variable map */
   SCIP_Bool             f                   /**< the value of x that activates the first relation */
   )
{
   SCIP_Real coefs2[3] = { 0., 0., 0. };
   SCIP_Real impllb;
   SCIP_Real implub;
   SCIP_VAR* binvar;
   SCIP_VAR* implvar;
   SCIP_Real side2;
   int i;
   SCIP_Bool binvals[2] = {!f, f};

   assert(binvarpos != implvarpos);
   assert(implvarpos != 0); /* implied variable must be continuous, therefore it can't be x */

   binvar = vars_xwy[binvarpos];
   implvar = vars_xwy[implvarpos];

   assert(SCIPvarGetType(binvar) == SCIP_VARTYPE_BINARY);
   assert(SCIPvarGetType(implvar) != SCIP_VARTYPE_BINARY);

   /* loop over binvals; if binvar is x (case binvarpos == 0), then we want to use only implications from
    * binvar == !f (which is the option complementing the first relation, which is implied from f); if
    * binvar is not x, this doesn't matter since the implbnd doesn't depend on x, therefore try both !f and f
    */
   for( i = 0; i < (binvarpos == 0 ? 1 : 2); ++i )
   {
      /* get implications binvar == binval  =>  implvar <=/>= implbnd */
      SCIPvarGetImplicVarBounds(binvar, binvals[i], implvar, &impllb, &implub);

      if( impllb != SCIP_INVALID ) /*lint !e777*/
      {
         /* write the implied bound as a big-M constraint */
         implBndToBigM(scip, vars_xwy, binvarpos, implvarpos, SCIP_BOUNDTYPE_LOWER, binvals[i], impllb, coefs2, &side2);

         SCIP_CALL( extractProducts(scip, sepadata, vars_xwy, coefs1, coefs2, side1, side2, sidetype1,
               SCIP_SIDETYPE_LEFT, varmap, f) );
      }

      if( implub != SCIP_INVALID ) /*lint !e777*/
      {
         /* write the implied bound as a big-M constraint */
         implBndToBigM(scip, vars_xwy, binvarpos, implvarpos, SCIP_BOUNDTYPE_UPPER, binvals[i], implub, coefs2, &side2);

         SCIP_CALL( extractProducts(scip, sepadata, vars_xwy, coefs1, coefs2, side1, side2, sidetype1,
               SCIP_SIDETYPE_RIGHT, varmap, f) );
      }
   }

   return SCIP_OKAY;
}

/** extract products from a relation given by `coefs1`, `vars_xwy`, `side1` and `sidetype1` and
 *  cliques containing `vars_xwy[varpos1]` and `vars_xwy[varpos2]`
 */
static
SCIP_RETCODE detectProductsClique(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_SEPADATA*        sepadata,           /**< separator data */
   SCIP_Real*            coefs1,             /**< coefficients of the first linear relation */
   SCIP_VAR**            vars_xwy,           /**< variables of the first relation in the order x, w, y */
   SCIP_Real             side1,              /**< side of the first relation */
   SCIP_SIDETYPE         sidetype1,          /**< is the left or right hand side given for the first relation? */
   int                   varpos1,            /**< position of the first variable in the vars_xwy array */
   int                   varpos2,            /**< position of the second variable in the vars_xwy array */
   SCIP_HASHMAP*         varmap,             /**< variable map */
   SCIP_Bool             f                   /**< the value of x that activates the first relation */
   )
{
   SCIP_Real coefs2[3] = { 0., 0., 0. };
   SCIP_VAR* var1;
   SCIP_VAR* var2;
   SCIP_Real side2;
   int i;
   int imax;
   SCIP_Bool binvals[2] = {!f, f};

   var1 = vars_xwy[varpos1];
   var2 = vars_xwy[varpos2];

   /* this decides whether we do one or two iterations of the loop for binvals: if var1
    * or var2 is x, we only want cliques with x = !f (which is the option complementing
    * the first relation, which is implied from f); otherwise this doesn't matter since
    * the clique doesn't depend on x, therefore try both !f and f
    */
   imax = (varpos1 == 0 || varpos2 == 0) ? 1 : 2;

   assert(SCIPvarGetType(var1) == SCIP_VARTYPE_BINARY);
   assert(SCIPvarGetType(var2) == SCIP_VARTYPE_BINARY);

   for( i = 0; i < imax; ++i )
   {
      /* if var1=TRUE and var2=TRUE are in a clique (binvals[i] == TRUE), the relation var1 + var2 <= 1 is implied
       * if var1=FALSE and var2=TRUE are in a clique (binvals[i] == FALSE), the relation (1 - var1) + var2 <= 1 is implied
       */
      if( SCIPvarsHaveCommonClique(var1, binvals[i], var2, TRUE, TRUE) )
      {
         SCIPdebugMsg(scip, "vars %s<%s> and <%s> are in a clique\n", binvals[i] ? "" : "!", SCIPvarGetName(var1), SCIPvarGetName(var2));
         coefs2[varpos1] = binvals[i] ? 1.0 : -1.0;
         coefs2[varpos2] = 1.0;
         side2 = binvals[i] ? 1.0 : 0.0;

         SCIP_CALL( extractProducts(scip, sepadata, vars_xwy, coefs1, coefs2, side1, side2, sidetype1,
               SCIP_SIDETYPE_RIGHT, varmap, f) );
      }

      /* if var1=TRUE and var2=FALSE are in the same clique, the relation var1 + (1-var2) <= 1 is implied
       * if var1=FALSE and var2=FALSE are in the same clique, the relation (1-var1) + (1-var2) <= 1 is implied
       */
      if( SCIPvarsHaveCommonClique(var1, binvals[i], var2, FALSE, TRUE) )
      {
         SCIPdebugMsg(scip, "vars %s<%s> and !<%s> are in a clique\n", binvals[i] ? "" : "!", SCIPvarGetName(var1), SCIPvarGetName(var2));
         coefs2[varpos1] = binvals[i] ? 1.0 : -1.0;
         coefs2[varpos2] = -1.0;
         side2 = binvals[i] ? 0.0 : -1.0;

         SCIP_CALL( extractProducts(scip, sepadata, vars_xwy, coefs1, coefs2, side1, side2, sidetype1,
               SCIP_SIDETYPE_RIGHT, varmap, f) );
      }
   }

   return SCIP_OKAY;
}


/** extract products from a relation given by `coefs1`, `vars`, `side1` and `sidetype1` and unconditional relations
 * (inequalities with 2 nonzeros) containing `vars[varpos1]` and `vars[varpos2]`
 */
static
SCIP_RETCODE detectProductsUnconditional(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_SEPADATA*        sepadata,           /**< separator data */
   SCIP_ROW**            rows,               /**< problem rows */
   int*                  row_list,           /**< linked list of rows corresponding to 2 or 3 var sets */
   SCIP_HASHTABLE*       hashtable,          /**< hashtable storing unconditional relations */
   SCIP_Real*            coefs1,             /**< coefficients of the first linear relation */
   SCIP_VAR**            vars_xwy,           /**< variables of the first relation in the order x, w, y */
   SCIP_Real             side1,              /**< side of the first relation */
   SCIP_SIDETYPE         sidetype1,          /**< is the left or right hand side given for the first relation? */
   int                   varpos1,            /**< position of the first unconditional variable in the vars_xwy array */
   int                   varpos2,            /**< position of the second unconditional variable in the vars_xwy array */
   SCIP_HASHMAP*         varmap,             /**< variable map */
   SCIP_Bool             f                   /**< the value of x that activates the first relation */
   )
{
   HASHDATA hashdata;
   HASHDATA* foundhashdata;
   SCIP_ROW* row2;
   int r2;
   int pos1;
   int pos2;
   SCIP_Real coefs2[3] = { 0., 0., 0. };
   SCIP_VAR* var1;
   SCIP_VAR* var2;

   /* always unconditional, therefore x must not be one of the two variables */
   assert(varpos1 != 0);
   assert(varpos2 != 0);

   var1 = vars_xwy[varpos1];
   var2 = vars_xwy[varpos2];

   hashdata.nvars = 2;
   hashdata.firstrow = -1;
   if( SCIPvarGetIndex(var1) < SCIPvarGetIndex(var2) )
   {
      pos1 = 0;
      pos2 = 1;
   }
   else
   {
      pos1 = 1;
      pos2 = 0;
   }

   hashdata.vars[pos1] = var1;
   hashdata.vars[pos2] = var2;

   foundhashdata = (HASHDATA*)SCIPhashtableRetrieve(hashtable, &hashdata);

   if( foundhashdata != NULL )
   {
      /* if the var pair exists, use all corresponding rows */
      r2 = foundhashdata->firstrow;

      while( r2 != -1 )
      {
         row2 = rows[r2];
         assert(SCIProwGetNNonz(row2) == 2);
         assert(var1 == SCIPcolGetVar(SCIProwGetCols(row2)[pos1]));
         assert(var2 == SCIPcolGetVar(SCIProwGetCols(row2)[pos2]));

         coefs2[varpos1] = SCIProwGetVals(row2)[pos1];
         coefs2[varpos2] = SCIProwGetVals(row2)[pos2];

         SCIPdebugMsg(scip, "Unconditional:\n");
         if( !SCIPisInfinity(scip, -SCIProwGetLhs(row2)) )
         {
            SCIP_CALL( extractProducts(scip, sepadata, vars_xwy, coefs1, coefs2, side1,
                  SCIProwGetLhs(row2) - SCIProwGetConstant(row2), sidetype1, SCIP_SIDETYPE_LEFT, varmap, f) );
         }
         if( !SCIPisInfinity(scip,  SCIProwGetRhs(row2)) )
         {
            SCIP_CALL( extractProducts(scip, sepadata, vars_xwy, coefs1, coefs2, side1,
                  SCIProwGetRhs(row2) - SCIProwGetConstant(row2), sidetype1, SCIP_SIDETYPE_RIGHT, varmap, f) );
         }

         r2 = row_list[r2];
      }
   }

   return SCIP_OKAY;
}

/** finds and stores implied relations (x = f &rArr; ay + bw &le; c, f can be 0 or 1) and 2-variable relations
 *
 *  Fills the following:
 *
 * - An array of variables that participate in two variable relations; for each such variable, ADJACENTVARDATA
 *   containing an array of variables that participate in two variable relations together with it; and a hashmap
 *   mapping variables to ADJACENTVARDATAs.
 *
 * - Hashtables storing hashdata objects with the two or three variables and the position of the first row in the
 *   `prob_rows` array, which in combination with the linked list (described below) will allow access to all rows that
 *   depend only on the corresponding variables.
 *
 * - Linked lists of row indices. Each list corresponds to a pair or triple of variables and contains positions of rows
 *   which depend only on those variables. All lists are stored in `row_list`, an array of length `nrows`, which is
 *   possible because each row belongs to at most one list. The array indices of `row_list` represent the positions of
 *   rows in `prob_rows`, and a value in the `row_list` array represents the next index in the list (-1 if there is no next
 *   list element). The first index of each list is stored in one of the hashdata objects as firstrow.
 */
static
SCIP_RETCODE fillRelationTables(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_ROW**            prob_rows,          /**< linear rows of the problem */
   int                   nrows,              /**< number of rows */
   SCIP_HASHTABLE*       hashtable2,         /**< hashtable to store 2-variable relations */
   SCIP_HASHTABLE*       hashtable3,         /**< hashtable to store implied relations */
   SCIP_HASHMAP*         vars_in_2rels,      /**< connections between variables that appear in 2-variable relations */
   int*                  row_list            /**< linked lists of row positions for each 2 or 3 variable set */
   )
{
   int r;
   SCIP_COL** cols;
   HASHDATA searchhashdata;
   HASHDATA* elementhashdata;

   assert(prob_rows != NULL);
   assert(nrows > 0);
   assert(hashtable2 != NULL);
   assert(hashtable3 != NULL);
   assert(vars_in_2rels != NULL);
   assert(row_list != NULL);

   for( r = 0; r < nrows; ++r )
   {
      assert(prob_rows[r] != NULL);

      cols = SCIProwGetCols(prob_rows[r]);
      assert(cols != NULL);

      /* initialise with the "end of list" value */
      row_list[r] = -1;

      /* look for unconditional relations with 2 variables */
      if( SCIProwGetNNonz(prob_rows[r]) == 2 )
      {
         /* if at least one of the variables is binary, this is either an implied bound
          * or a clique; these are covered separately */
         if( SCIPvarGetType(SCIPcolGetVar(cols[0])) == SCIP_VARTYPE_BINARY ||
             SCIPvarGetType(SCIPcolGetVar(cols[1])) == SCIP_VARTYPE_BINARY )
         {
            SCIPdebugMsg(scip, "ignoring relation <%s> because a var is binary\n", SCIProwGetName(prob_rows[r]));
            continue;
         }

         /* fill in searchhashdata so that to search for the two variables in hashtable2 */
         searchhashdata.nvars = 2;
         searchhashdata.firstrow = -1;
         searchhashdata.vars[0] = SCIPcolGetVar(cols[0]);
         searchhashdata.vars[1] = SCIPcolGetVar(cols[1]);

         /* get the element corresponding to the two variables */
         elementhashdata = (HASHDATA*)SCIPhashtableRetrieve(hashtable2, &searchhashdata);

         if( elementhashdata != NULL )
         {
            /* if element exists, update it by adding the row */
            row_list[r] = elementhashdata->firstrow;
            elementhashdata->firstrow = r;
            ++elementhashdata->nrows;
         }
         else
         {
            /* create an element for the combination of two variables */
            SCIP_CALL( SCIPallocBuffer(scip, &elementhashdata) );

            elementhashdata->nvars = 2;
            elementhashdata->nrows = 1;
            elementhashdata->vars[0] = searchhashdata.vars[0];
            elementhashdata->vars[1] = searchhashdata.vars[1];
            elementhashdata->firstrow = r;

            SCIP_CALL( SCIPhashtableInsert(hashtable2, (void*)elementhashdata) );

            /* hashdata.vars are two variables participating together in a two variable relation, therefore update
             * these variables' adjacency data
             */
            SCIP_CALL( addAdjacentVars(scip, vars_in_2rels, searchhashdata.vars) );
         }
      }

      /* look for implied relations (three variables, at least one binary variable) */
      if( SCIProwGetNNonz(prob_rows[r]) == 3 )
      {
         /* an implied relation contains at least one binary variable */
         if( SCIPvarGetType(SCIPcolGetVar(cols[0])) != SCIP_VARTYPE_BINARY &&
             SCIPvarGetType(SCIPcolGetVar(cols[1])) != SCIP_VARTYPE_BINARY &&
             SCIPvarGetType(SCIPcolGetVar(cols[2])) != SCIP_VARTYPE_BINARY )
            continue;

         /* fill in hashdata so that to search for the three variables in hashtable3 */
         searchhashdata.nvars = 3;
         searchhashdata.firstrow = -1;
         searchhashdata.vars[0] = SCIPcolGetVar(cols[0]);
         searchhashdata.vars[1] = SCIPcolGetVar(cols[1]);
         searchhashdata.vars[2] = SCIPcolGetVar(cols[2]);

         /* get the element corresponding to the three variables */
         elementhashdata = (HASHDATA*)SCIPhashtableRetrieve(hashtable3, &searchhashdata);

         if( elementhashdata != NULL )
         {
            /* if element exists, update it by adding the row */
            row_list[r] = elementhashdata->firstrow;
            elementhashdata->firstrow = r;
            ++elementhashdata->nrows;
         }
         else
         {
            /* create an element for the combination of three variables */
            SCIP_CALL( SCIPallocBuffer(scip, &elementhashdata) );

            elementhashdata->nvars = 3;
            elementhashdata->nrows = 1;
            elementhashdata->vars[0] = searchhashdata.vars[0];
            elementhashdata->vars[1] = searchhashdata.vars[1];
            elementhashdata->vars[2] = searchhashdata.vars[2];
            elementhashdata->firstrow = r;

            SCIP_CALL( SCIPhashtableInsert(hashtable3, (void*)elementhashdata) );
         }
      }
   }

   return SCIP_OKAY;
}

/** detect bilinear products encoded in linear constraints */
static
SCIP_RETCODE detectHiddenProducts(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_SEPADATA*        sepadata,           /**< separation data */
   SCIP_HASHMAP*         varmap              /**< variable map */
   )
{
   int r1; /* first relation index */
   int r2; /* second relation index */
   int i; /* outer loop counter */
   int permwy; /* index for permuting w and y */
   int nrows;
   SCIP_ROW** prob_rows;
   SCIP_HASHTABLE* hashtable3;
   SCIP_HASHTABLE* hashtable2;
   HASHDATA* foundhashdata;
   SCIP_VAR* vars_xwy[3];
   SCIP_Real coefs1[3];
   SCIP_Real coefs2[3];
   SCIP_ROW* row1;
   SCIP_ROW* row2;
   int xpos;
   int ypos;
   int wpos;
   int f; /* value of the binary variable */
   SCIP_VAR** relatedvars;
   int nrelatedvars;
   SCIP_Bool xfixing;
   SCIP_SIDETYPE sidetype1;
   SCIP_SIDETYPE sidetype2;
   SCIP_Real side1;
   SCIP_Real side2;
   int* row_list;
   SCIP_HASHMAP* vars_in_2rels;
   int nvars;

   /* get the (original) rows */
   SCIP_CALL( getOriginalRows(scip, &prob_rows, &nrows) );

   if( nrows == 0 )
   {
      SCIPfreeBufferArray(scip, &prob_rows);
      return SCIP_OKAY;
   }

   /* create tables of implied and unconditional relations */
   SCIP_CALL( SCIPhashtableCreate(&hashtable3, SCIPblkmem(scip), nrows, SCIPhashGetKeyStandard,
      hashdataKeyEqConss, hashdataKeyValConss, NULL) );
   SCIP_CALL( SCIPhashtableCreate(&hashtable2, SCIPblkmem(scip), nrows, SCIPhashGetKeyStandard,
      hashdataKeyEqConss, hashdataKeyValConss, NULL) );
   SCIP_CALL( SCIPallocBufferArray(scip, &row_list, nrows) );

   /* allocate the adjacency data map for variables that appear in 2-var relations */
   nvars = SCIPgetNVars(scip);
   SCIP_CALL( SCIPhashmapCreate(&vars_in_2rels, SCIPblkmem(scip), MIN(nvars, nrows * 2)) );

   /* fill the data structures that will be used for product detection: hashtables and linked lists allowing to access
    * two and three variable relations by the variables; and the hashmap for accessing variables participating in two
    * variable relations with each given variable */
   SCIP_CALL( fillRelationTables(scip, prob_rows, nrows, hashtable2, hashtable3, vars_in_2rels, row_list) );

   /* start actually looking for products */
   /* go through all sets of three variables */
   for( i = 0; i < SCIPhashtableGetNEntries(hashtable3); ++i )
   {
      foundhashdata = (HASHDATA*)SCIPhashtableGetEntry(hashtable3, i);
      if( foundhashdata == NULL )
         continue;

      SCIPdebugMsg(scip, "(<%s>, <%s>, <%s>): ", SCIPvarGetName(foundhashdata->vars[0]),
         SCIPvarGetName(foundhashdata->vars[1]), SCIPvarGetName(foundhashdata->vars[2]));

      /* An implied relation has the form: x == f  =>  l(w,y) <=/>= side (f is 0 or 1, l is a linear function). Given
       * a linear relation with three variables, any binary var can be x: we try them all here because this can
       * produce different products.
       */
      for( xpos = 0; xpos < 3; ++xpos )
      {
         /* in vars_xwy, the order of variables is always as in the name: x, w, y */
         vars_xwy[0] = foundhashdata->vars[xpos];

         /* x must be binary */
         if( SCIPvarGetType(vars_xwy[0]) != SCIP_VARTYPE_BINARY )
            continue;

         /* the first row might be an implication from f == 0 or f == 1: try both */
         for( f = 0; f <= 1; ++f )
         {
            xfixing = f == 1;

            /* go through implied relations for the corresponding three variables */
            for( r1 = foundhashdata->firstrow; r1 != -1; r1 = row_list[r1] )
            {
               /* get the implied relation */
               row1 = prob_rows[r1];

               assert(SCIProwGetNNonz(row1) == 3);
               /* the order of variables in all rows should be the same, and similar to the order in hashdata->vars,
                * therefore the x variable from vars_xwy should be similar to the column variable at xpos
                */
               assert(vars_xwy[0] == SCIPcolGetVar(SCIProwGetCols(row1)[xpos]));

               coefs1[0] = SCIProwGetVals(row1)[xpos];

               /* use the side for which the inequality becomes tighter when x == xfixing than when x == !xfixing */
               if( (!xfixing && coefs1[0] > 0.0) || (xfixing && coefs1[0] < 0.0) )
               {
                  sidetype1 = SCIP_SIDETYPE_LEFT;
                  side1 = SCIProwGetLhs(row1);
               }
               else
               {
                  sidetype1 = SCIP_SIDETYPE_RIGHT;
                  side1 = SCIProwGetRhs(row1);
               }

               if( SCIPisInfinity(scip, REALABS(side1)) )
                  continue;

               side1 -= SCIProwGetConstant(row1);

               /* permute w and y */
               for( permwy = 1; permwy <= 2; ++permwy )
               {
                  wpos = (xpos + permwy) % 3;
                  ypos = (xpos - permwy + 3) % 3;
                  vars_xwy[1] = foundhashdata->vars[wpos];
                  vars_xwy[2] = foundhashdata->vars[ypos];

                  assert(vars_xwy[1] == SCIPcolGetVar(SCIProwGetCols(row1)[wpos]));
                  assert(vars_xwy[2] == SCIPcolGetVar(SCIProwGetCols(row1)[ypos]));

                  coefs1[1] = SCIProwGetVals(row1)[wpos];
                  coefs1[2] = SCIProwGetVals(row1)[ypos];

                  /* look for the second relation: it should be tighter when x == !xfixing than when x == xfixing
                   * and can be either another implied relation or one of several types of two and one variable
                   * relations
                   */

                  /* go through the remaining rows (implied relations) for these three variables */
                  for( r2 = row_list[r1]; r2 != -1; r2 = row_list[r2] )
                  {
                     /* get the second implied relation */
                     row2 = prob_rows[r2];

                     assert(SCIProwGetNNonz(row2) == 3);
                     assert(vars_xwy[0] == SCIPcolGetVar(SCIProwGetCols(row2)[xpos]));
                     assert(vars_xwy[1] == SCIPcolGetVar(SCIProwGetCols(row2)[wpos]));
                     assert(vars_xwy[2] == SCIPcolGetVar(SCIProwGetCols(row2)[ypos]));

                     coefs2[0] = SCIProwGetVals(row2)[xpos];
                     coefs2[1] = SCIProwGetVals(row2)[wpos];
                     coefs2[2] = SCIProwGetVals(row2)[ypos];

                     /* use the side for which the inequality becomes tighter when x == !xfixing than when x == xfixing */
                     if( (!xfixing && coefs2[0] > 0.0) || (xfixing && coefs2[0] < 0.0) )
                     {
                        sidetype2 = SCIP_SIDETYPE_RIGHT;
                        side2 = SCIProwGetRhs(row2);
                     }
                     else
                     {
                        sidetype2 = SCIP_SIDETYPE_LEFT;
                        side2 = SCIProwGetLhs(row2);
                     }

                     if( SCIPisInfinity(scip, REALABS(side2)) )
                        continue;

                     side2 -= SCIProwGetConstant(row2);

                     SCIPdebugMsg(scip, "Two implied relations:\n");
                     SCIP_CALL( extractProducts(scip, sepadata, vars_xwy, coefs1, coefs2, side1, side2, sidetype1,
                        sidetype2, varmap, xfixing) );
                  }

                  /* use global bounds on w */
                  coefs2[0] = 0.0;
                  coefs2[1] = 1.0;
                  coefs2[2] = 0.0;
                  SCIPdebugMsg(scip, "w global bounds:\n");
                  if( !SCIPisInfinity(scip, -SCIPvarGetLbGlobal(vars_xwy[1])) )
                  {
                     SCIP_CALL( extractProducts(scip, sepadata, vars_xwy, coefs1, coefs2, side1,
                        SCIPvarGetLbGlobal(vars_xwy[1]), sidetype1, SCIP_SIDETYPE_LEFT, varmap, xfixing) );
                  }

                  if( !SCIPisInfinity(scip, SCIPvarGetUbGlobal(vars_xwy[1])) )
                  {
                     SCIP_CALL( extractProducts(scip, sepadata, vars_xwy, coefs1, coefs2, side1,
                        SCIPvarGetUbGlobal(vars_xwy[1]), sidetype1, SCIP_SIDETYPE_RIGHT, varmap, xfixing) );
                  }

                  /* use implied bounds and cliques with w */
                  if( SCIPvarGetType(vars_xwy[1]) != SCIP_VARTYPE_BINARY )
                  {
                     /* w is non-binary - look for implied bounds x == !f => w >=/<= bound */
                     SCIPdebugMsg(scip, "Implied relation + implied bounds on w:\n");
                     SCIP_CALL( detectProductsImplbnd(scip, sepadata, coefs1, vars_xwy, side1, sidetype1, 0, 1,
                           varmap, xfixing) );
                  }
                  else
                  {
                     /* w is binary - look for cliques containing x and w */
                     SCIPdebugMsg(scip, "Implied relation + cliques with x and w:\n");
                     SCIP_CALL( detectProductsClique(scip, sepadata, coefs1, vars_xwy, side1, sidetype1, 0, 1,
                           varmap, xfixing) );
                  }

                  /* use unconditional relations (i.e. relations of w and y) */

                  /* implied bound w == 0/1 => y >=/<= bound */
                  if( SCIPvarGetType(vars_xwy[1]) == SCIP_VARTYPE_BINARY && SCIPvarGetType(vars_xwy[2]) != SCIP_VARTYPE_BINARY )
                  {
                     SCIPdebugMsg(scip, "Implied relation + implied bounds with w and y:\n");
                     SCIP_CALL( detectProductsImplbnd(scip, sepadata, coefs1, vars_xwy, side1, sidetype1, 1, 2, varmap, xfixing) );
                  }

                  /* implied bound y == 0/1 => w >=/<= bound */
                  if( SCIPvarGetType(vars_xwy[2]) == SCIP_VARTYPE_BINARY && SCIPvarGetType(vars_xwy[1]) != SCIP_VARTYPE_BINARY )
                  {
                     SCIPdebugMsg(scip, "Implied relation + implied bounds with y and w:\n");
                     SCIP_CALL( detectProductsImplbnd(scip, sepadata, coefs1, vars_xwy, side1, sidetype1, 2, 1, varmap, xfixing) );
                  }

                  /* cliques containing w and y */
                  if( SCIPvarGetType(vars_xwy[1]) == SCIP_VARTYPE_BINARY && SCIPvarGetType(vars_xwy[2]) == SCIP_VARTYPE_BINARY )
                  {
                     SCIPdebugMsg(scip, "Implied relation + cliques with w and y:\n");
                     SCIP_CALL( detectProductsClique(scip, sepadata, coefs1, vars_xwy, side1, sidetype1, 1, 2, varmap, xfixing) );
                  }

                  /* inequalities containing w and y */
                  if( SCIPvarGetType(vars_xwy[1]) != SCIP_VARTYPE_BINARY && SCIPvarGetType(vars_xwy[2]) != SCIP_VARTYPE_BINARY )
                  {
                     SCIPdebugMsg(scip, "Implied relation + unconditional with w and y:\n");
                     SCIP_CALL( detectProductsUnconditional(scip, sepadata, prob_rows, row_list, hashtable2, coefs1,
                        vars_xwy, side1, sidetype1, 1, 2, varmap, xfixing) );
                  }
               }
            }
         }
      }
      SCIPfreeBuffer(scip, &foundhashdata);
   }

   /* also loop through implied bounds to look for products */
   for( i = 0; i < SCIPgetNBinVars(scip); ++i )
   {
      /* first choose the x variable: it can be any binary variable in the problem */
      vars_xwy[0] = SCIPgetVars(scip)[i];

      assert(SCIPvarGetType(vars_xwy[0]) == SCIP_VARTYPE_BINARY);

      /* consider both possible values of x */
      for( f = 0; f <= 1; ++f )
      {
         xfixing = f == 1;

         /* go through implications of x */
         for( r1 = 0; r1 < SCIPvarGetNImpls(vars_xwy[0], xfixing); ++r1 )
         {
            /* w is the implication var */
            vars_xwy[1] = SCIPvarGetImplVars(vars_xwy[0], xfixing)[r1];
            assert(SCIPvarGetType(vars_xwy[1]) != SCIP_VARTYPE_BINARY);

            /* write the implication as a big-M constraint */
            implBndToBigM(scip, vars_xwy, 0, 1, SCIPvarGetImplTypes(vars_xwy[0], xfixing)[r1], xfixing,
                  SCIPvarGetImplBounds(vars_xwy[0], xfixing)[r1], coefs1, &side1);
            sidetype1 = SCIPvarGetImplTypes(vars_xwy[0], xfixing)[r1] == SCIP_BOUNDTYPE_LOWER ?
                     SCIP_SIDETYPE_LEFT : SCIP_SIDETYPE_RIGHT;

            /* if the global bound is equal to the implied bound, there is nothing to do */
            if( SCIPisZero(scip, coefs1[0]) )
               continue;

            SCIPdebugMsg(scip, "Implication %s == %u  =>  %s %s %g\n", SCIPvarGetName(vars_xwy[0]), xfixing,
                  SCIPvarGetName(vars_xwy[1]), sidetype1 == SCIP_SIDETYPE_LEFT ? ">=" : "<=",
                  SCIPvarGetImplBounds(vars_xwy[0], xfixing)[r1]);
            SCIPdebugMsg(scip, "Written as big-M: %g%s + %s %s %g\n", coefs1[0], SCIPvarGetName(vars_xwy[0]),
                  SCIPvarGetName(vars_xwy[1]), sidetype1 == SCIP_SIDETYPE_LEFT ? ">=" : "<=", side1);

            /* the second relation is in w and y (y could be anything, but must be in relation with w) */

            /* x does not participate in the second relation, so we immediately set its coefficient to 0.0 */
            coefs2[0] = 0.0;

            SCIPdebugMsg(scip, "Implic of x = <%s> + implied lb on w = <%s>:\n", SCIPvarGetName(vars_xwy[0]), SCIPvarGetName(vars_xwy[1]));

            /* use implied lower bounds on w: w >= b*y + d */
            for( r2 = 0; r2 < SCIPvarGetNVlbs(vars_xwy[1]); ++r2 )
            {
               vars_xwy[2] = SCIPvarGetVlbVars(vars_xwy[1])[r2];
               if( vars_xwy[2] == vars_xwy[0] )
                  continue;

               coefs2[1] = 1.0;
               coefs2[2] = -SCIPvarGetVlbCoefs(vars_xwy[1])[r2];

               SCIP_CALL( extractProducts(scip, sepadata, vars_xwy, coefs1, coefs2, side1,
                  SCIPvarGetVlbConstants(vars_xwy[1])[r2], sidetype1, SCIP_SIDETYPE_LEFT, varmap, xfixing) );
            }

            SCIPdebugMsg(scip, "Implic of x = <%s> + implied ub on w = <%s>:\n", SCIPvarGetName(vars_xwy[0]), SCIPvarGetName(vars_xwy[1]));

            /* use implied upper bounds on w: w <= b*y + d */
            for( r2 = 0; r2 < SCIPvarGetNVubs(vars_xwy[1]); ++r2 )
            {
               vars_xwy[2] = SCIPvarGetVubVars(vars_xwy[1])[r2];
               if( vars_xwy[2] == vars_xwy[0] )
                  continue;

               coefs2[1] = 1.0;
               coefs2[2] = -SCIPvarGetVubCoefs(vars_xwy[1])[r2];

               SCIP_CALL( extractProducts(scip, sepadata, vars_xwy, coefs1, coefs2, side1,
                  SCIPvarGetVubConstants(vars_xwy[1])[r2], sidetype1, SCIP_SIDETYPE_RIGHT, varmap, xfixing) );
            }

            /* use unconditional relations containing w */
            relatedvars = getAdjacentVars(vars_in_2rels, vars_xwy[1], &nrelatedvars);
            if( relatedvars == NULL )
               continue;

            for( r2 = 0; r2 < nrelatedvars; ++r2 )
            {
               vars_xwy[2] = relatedvars[r2];
               SCIPdebugMsg(scip, "Implied bound + unconditional with w and y:\n");
               SCIP_CALL( detectProductsUnconditional(scip, sepadata, prob_rows, row_list, hashtable2, coefs1,
                  vars_xwy, side1, sidetype1, 1, 2, varmap, xfixing) );
            }
         }
      }
   }

   /* free memory */
   clearVarAdjacency(scip, vars_in_2rels);
   SCIPhashmapFree(&vars_in_2rels);

   SCIPdebugMsg(scip, "Unconditional relations table:\n");
   for( i = 0; i < SCIPhashtableGetNEntries(hashtable2); ++i )
   {
      foundhashdata = (HASHDATA*)SCIPhashtableGetEntry(hashtable2, i);
      if( foundhashdata == NULL )
         continue;

      SCIPdebugMsg(scip, "(%s, %s): ", SCIPvarGetName(foundhashdata->vars[0]),
                   SCIPvarGetName(foundhashdata->vars[1]));

      SCIPfreeBuffer(scip, &foundhashdata);
   }

   SCIPfreeBufferArray(scip, &row_list);

   SCIPhashtableFree(&hashtable2);
   SCIPhashtableFree(&hashtable3);

   SCIPfreeBufferArray(scip, &prob_rows);

   return SCIP_OKAY;
}

/** helper method to create separation data */
static
SCIP_RETCODE createSepaData(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_SEPADATA*        sepadata            /**< separation data */
   )
{
   SCIP_HASHMAP* varmap;
   int i;
   SCIP_CONSNONLINEAR_BILINTERM* bilinterms;
   int varmapsize;
   int nvars;

   assert(sepadata != NULL);

   /* initialize some fields of sepadata */
   sepadata->varssorted = NULL;
   sepadata->varpriorities = NULL;
   sepadata->bilinvardatamap = NULL;
   sepadata->eqauxexpr = NULL;
   sepadata->nbilinvars = 0;
   sepadata->sbilinvars = 0;

   /* get total number of bilinear terms */
   sepadata->nbilinterms = SCIPgetNBilinTermsNonlinear(sepadata->conshdlr);

   /* skip if there are no bilinear terms and implicit product detection is off */
   if( sepadata->nbilinterms == 0 && !sepadata->detecthidden )
      return SCIP_OKAY;

   /* the number of variables participating in bilinear products cannot exceed twice the number of bilinear terms;
    * however, if we detect hidden products, the number of terms is yet unknown, so use the number of variables
    */
   nvars = SCIPgetNVars(scip);
   varmapsize = sepadata->detecthidden ? nvars : MIN(nvars, sepadata->nbilinterms * 2);

   /* create variable map */
   SCIP_CALL( SCIPhashmapCreate(&varmap, SCIPblkmem(scip), varmapsize) );

   /* get all bilinear terms from the nonlinear constraint handler */
   bilinterms = SCIPgetBilinTermsNonlinear(sepadata->conshdlr);

   /* store the information of all variables that appear bilinearly */
   for( i = 0; i < sepadata->nbilinterms; ++i )
   {
      assert(bilinterms[i].x != NULL);
      assert(bilinterms[i].y != NULL);
      assert(bilinterms[i].nlockspos + bilinterms[i].nlocksneg > 0);

      /* skip bilinear term if it does not have an auxiliary variable */
      if( bilinterms[i].aux.var == NULL )
         continue;

      /* if only original variables should be used, skip products that contain at least one auxiliary variable */
      if( sepadata->onlyoriginal && (SCIPvarIsRelaxationOnly(bilinterms[i].x) ||
          SCIPvarIsRelaxationOnly(bilinterms[i].y)) )
         continue;

      /* coverity[var_deref_model] */
      SCIP_CALL( addProductVars(scip, sepadata, bilinterms[i].x, bilinterms[i].y, varmap,
            bilinterms[i].nlockspos + bilinterms[i].nlocksneg) );
   }

   if( sepadata->detecthidden )
   {
      int oldnterms = sepadata->nbilinterms;

      /* coverity[var_deref_model] */
      SCIP_CALL( detectHiddenProducts(scip, sepadata, varmap) );

      /* update nbilinterms and bilinterms, as detectHiddenProducts might have found new terms */
      sepadata->nbilinterms = SCIPgetNBilinTermsNonlinear(sepadata->conshdlr);
      bilinterms = SCIPgetBilinTermsNonlinear(sepadata->conshdlr);

      if( sepadata->nbilinterms > oldnterms )
      {
         SCIPstatisticMessage(" Number of hidden products: %d\n", sepadata->nbilinterms - oldnterms);
      }
   }

   SCIPhashmapFree(&varmap);

   if( sepadata->nbilinterms == 0 )
   {
      return SCIP_OKAY;
   }

   /* mark positions of aux.exprs that must be equal to the product */
   SCIP_CALL( SCIPallocBlockMemoryArray(scip, &sepadata->eqauxexpr, sepadata->nbilinterms) );

   for( i = 0; i < sepadata->nbilinterms; ++i )
   {
      int j;

      sepadata->eqauxexpr[i] = -1;
      for( j = 0; j < bilinterms[i].nauxexprs; ++j )
      {
         assert(bilinterms[i].aux.exprs[j] != NULL);

         if( bilinterms[i].aux.exprs[j]->underestimate && bilinterms[i].aux.exprs[j]->overestimate )
         {
            sepadata->eqauxexpr[i] = j;
            break;
         }
      }
   }

   /* find maxnumber of variables that occur most often and sort them by number of occurrences
    * (same as normal sort, except that entries at positions maxusedvars..nbilinvars may be unsorted at end)
    */
   SCIPselectDownIntPtr(sepadata->varpriorities, (void**) sepadata->varssorted, MIN(sepadata->maxusedvars,sepadata->nbilinvars-1),
         sepadata->nbilinvars);

   /* capture all variables */
   for( i = 0; i < sepadata->nbilinvars; ++i )
   {
      assert(sepadata->varssorted[i] != NULL);
      SCIP_CALL( SCIPcaptureVar(scip, sepadata->varssorted[i]) );
   }

   /* mark that separation data has been created */
   sepadata->iscreated = TRUE;
   sepadata->isinitialround = TRUE;

   if( SCIPgetNBilinTermsNonlinear(sepadata->conshdlr) > 0 )
      SCIPstatisticMessage(" Found bilinear terms\n");
   else
      SCIPstatisticMessage(" No bilinear terms\n");

   return SCIP_OKAY;
}

/** get the positions of the most violated auxiliary under- and overestimators for each product
 *
 * -1 means no relation with given product is violated
 */
static
void getBestEstimators(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_SEPADATA*        sepadata,           /**< separator data */
   SCIP_SOL*             sol,                /**< solution at which to evaluate the expressions */
   int*                  bestunderestimators,/**< array of indices of best underestimators for each term */
   int*                  bestoverestimators  /**< array of indices of best overestimators for each term */
   )
{
   SCIP_Real prodval;
   SCIP_Real auxval;
   SCIP_Real prodviol;
   SCIP_Real viol_below;
   SCIP_Real viol_above;
   int i;
   int j;
   SCIP_CONSNONLINEAR_BILINTERM* terms;

   assert(bestunderestimators != NULL);
   assert(bestoverestimators != NULL);

   terms = SCIPgetBilinTermsNonlinear(sepadata->conshdlr);

   for( j = 0; j < SCIPgetNBilinTermsNonlinear(sepadata->conshdlr); ++j )
   {
      viol_below = 0.0;
      viol_above = 0.0;

      /* evaluate the product expression */
      prodval = SCIPgetSolVal(scip, sol, terms[j].x) * SCIPgetSolVal(scip, sol, terms[j].y);

      bestunderestimators[j] = -1;
      bestoverestimators[j] = -1;

      /* if there are any auxexprs, look there */
      for( i = 0; i < terms[j].nauxexprs; ++i )
      {
         auxval = SCIPevalBilinAuxExprNonlinear(scip, terms[j].x, terms[j].y, terms[j].aux.exprs[i], sol);
         prodviol = auxval - prodval;

         if( terms[j].aux.exprs[i]->underestimate && SCIPisFeasGT(scip, auxval, prodval) && prodviol > viol_below )
         {
            viol_below = prodviol;
            bestunderestimators[j] = i;
         }
         if( terms[j].aux.exprs[i]->overestimate && SCIPisFeasGT(scip, prodval, auxval) && -prodviol > viol_above )
         {
            viol_above = -prodviol;
            bestoverestimators[j] = i;
         }
      }

      /* if the term has a plain auxvar, it will be treated differently - do nothing here */
   }
}

/** tests if a row contains too many unknown bilinear terms w.r.t. the parameters */
static
SCIP_RETCODE isAcceptableRow(
   SCIP_SEPADATA*        sepadata,           /**< separation data */
   SCIP_ROW*             row,                /**< the row to be tested */
   SCIP_VAR*             var,                /**< the variable that is to be multiplied with row */
   int*                  currentnunknown,    /**< buffer to store number of unknown terms in current row if acceptable */
   SCIP_Bool*            acceptable          /**< buffer to store the result */
   )
{
   int i;
   int idx;
   SCIP_CONSNONLINEAR_BILINTERM* terms;

   assert(row != NULL);
   assert(var != NULL);

   *currentnunknown = 0;
   terms = SCIPgetBilinTermsNonlinear(sepadata->conshdlr);

   for( i = 0; (i < SCIProwGetNNonz(row)) && (sepadata->maxunknownterms < 0 || *currentnunknown <= sepadata->maxunknownterms); ++i )
   {
      idx = SCIPgetBilinTermIdxNonlinear(sepadata->conshdlr, var, SCIPcolGetVar(SCIProwGetCols(row)[i]));

      /* if the product hasn't been found, no auxiliary expressions for it are known */
      if( idx < 0 )
      {
         ++(*currentnunknown);
         continue;
      }

      /* known terms are only those that have an aux.var or equality estimators */
      if( sepadata->eqauxexpr[idx] == -1 && !(terms[idx].nauxexprs == 0 && terms[idx].aux.var != NULL) )
      {
         ++(*currentnunknown);
      }
   }

   *acceptable = sepadata->maxunknownterms < 0 || *currentnunknown <= sepadata->maxunknownterms;

   return SCIP_OKAY;
}

/** adds coefficients and constant of an auxiliary expression
 *
 *  the variables the pointers are pointing to must already be initialized
 */
static
void addAuxexprCoefs(
   SCIP_VAR*             var1,               /**< first product variable */
   SCIP_VAR*             var2,               /**< second product variable */
   SCIP_CONSNONLINEAR_AUXEXPR* auxexpr,      /**< auxiliary expression to be added */
   SCIP_Real             coef,               /**< coefficient of the auxiliary expression */
   SCIP_Real*            coefaux,            /**< pointer to add the coefficient of the auxiliary variable */
   SCIP_Real*            coef1,              /**< pointer to add the coefficient of the first variable */
   SCIP_Real*            coef2,              /**< pointer to add the coefficient of the second variable */
   SCIP_Real*            cst                 /**< pointer to add the constant */
   )
{
   assert(auxexpr != NULL);
   assert(auxexpr->auxvar != NULL);
   assert(coefaux != NULL);
   assert(coef1 != NULL);
   assert(coef2 != NULL);
   assert(cst != NULL);

   *coefaux += auxexpr->coefs[0] * coef;

   /* in auxexpr, x goes before y and has the smaller index,
    * so compare vars to figure out which one is x and which is y
    */
   if( SCIPvarCompare(var1, var2) < 1 )
   {
      *coef1 += auxexpr->coefs[1] * coef;
      *coef2 += auxexpr->coefs[2] * coef;
   }
   else
   {
      *coef1 += auxexpr->coefs[2] * coef;
      *coef2 += auxexpr->coefs[1] * coef;
   }
   *cst += coef * auxexpr->cst;
}

/** add a linear term `coef`*`colvar` multiplied by a bound factor (var - lb(var)) or (ub(var) - var)
 *
 *  adds the linear term with `colvar` to `cut` and updates `coefvar` and `cst`
 */
static
SCIP_RETCODE addRltTerm(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_SEPADATA*        sepadata,           /**< separator data */
   SCIP_SOL*             sol,                /**< the point to be separated (can be NULL) */
   int*                  bestunderest,       /**< positions of most violated underestimators for each product term */
   int*                  bestoverest,        /**< positions of most violated overestimators for each product term */
   SCIP_ROW*             cut,                /**< cut to which the term is to be added */
   SCIP_VAR*             var,                /**< multiplier variable */
   SCIP_VAR*             colvar,             /**< row variable to be multiplied */
   SCIP_Real             coef,               /**< coefficient of the bilinear term */
   SCIP_Bool             uselb,              /**< whether we multiply with (var - lb) or (ub - var) */
   SCIP_Bool             uselhs,             /**< whether to create a cut for the lhs or rhs */
   SCIP_Bool             local,              /**< whether local or global cuts should be computed */
   SCIP_Bool             computeEqCut,       /**< whether conditions are fulfilled to compute equality cuts */
   SCIP_Real*            coefvar,            /**< coefficient of var */
   SCIP_Real*            cst,                /**< buffer to store the constant part of the cut */
   SCIP_Bool*            success             /**< buffer to store whether cut was updated successfully */
   )
{
   SCIP_Real lbvar;
   SCIP_Real ubvar;
   SCIP_Real refpointvar;
   SCIP_Real signfactor;
   SCIP_Real boundfactor;
   SCIP_Real coefauxvar;
   SCIP_Real coefcolvar;
   SCIP_Real coefterm;
   int auxpos;
   int idx;
   SCIP_CONSNONLINEAR_BILINTERM* terms;
   SCIP_VAR* auxvar;

   terms = SCIPgetBilinTermsNonlinear(sepadata->conshdlr);

   if( computeEqCut )
   {
      lbvar = 0.0;
      ubvar = 0.0;
   }
   else
   {
      lbvar = local ? SCIPvarGetLbLocal(var) : SCIPvarGetLbGlobal(var);
      ubvar = local ? SCIPvarGetUbLocal(var) : SCIPvarGetUbGlobal(var);
   }

   refpointvar = MAX(lbvar, MIN(ubvar, SCIPgetSolVal(scip, sol, var))); /*lint !e666*/

   signfactor = (uselb ? 1.0 : -1.0);
   boundfactor = (uselb ? -lbvar : ubvar);

   coefterm = coef * signfactor; /* coefficient of the bilinear term */
   coefcolvar = coef * boundfactor; /* coefficient of the linear term */
   coefauxvar = 0.0; /* coefficient of the auxiliary variable corresponding to the bilinear term */
   auxvar = NULL;

   assert(!SCIPisInfinity(scip, REALABS(coefterm)));

   /* first, add the linearisation of the bilinear term */

   idx = SCIPgetBilinTermIdxNonlinear(sepadata->conshdlr, var, colvar);
   auxpos = -1;

   /* for an implicit term, get the position of the best estimator */
   if( idx >= 0 && terms[idx].nauxexprs > 0 )
   {
      if( computeEqCut )
      {
         /* use an equality auxiliary expression (which should exist for computeEqCut to be TRUE) */
         assert(sepadata->eqauxexpr[idx] >= 0);
         auxpos = sepadata->eqauxexpr[idx];
      }
      else if( (uselhs && coefterm > 0.0) || (!uselhs && coefterm < 0.0) )
      {
         /* use an overestimator */
         auxpos = bestoverest[idx];
      }
      else
      {
         /* use an underestimator */
         auxpos = bestunderest[idx];
      }
   }

   /* if the term is implicit and a suitable auxiliary expression for var*colvar exists, add the coefficients
    * of the auxiliary expression for coefterm*var*colvar to coefauxvar, coefcolvar, coefvar and cst
    */
   if( auxpos >= 0 )
   {
      SCIPdebugMsg(scip, "auxiliary expression for <%s> and <%s> found, will be added to cut:\n",
                          SCIPvarGetName(colvar), SCIPvarGetName(var));
      addAuxexprCoefs(var, colvar, terms[idx].aux.exprs[auxpos], coefterm, &coefauxvar, coefvar, &coefcolvar, cst);
      auxvar = terms[idx].aux.exprs[auxpos]->auxvar;
   }
   /* for an existing term, use the auxvar if there is one */
   else if( idx >= 0 && terms[idx].nauxexprs == 0 && terms[idx].aux.var != NULL )
   {
      SCIPdebugMsg(scip, "auxvar for <%s> and <%s> found, will be added to cut:\n",
                   SCIPvarGetName(colvar), SCIPvarGetName(var));
      coefauxvar += coefterm;
      auxvar = terms[idx].aux.var;
   }

   /* otherwise, use clique information or the McCormick estimator in place of the bilinear term */
   else if( colvar != var )
   {
      SCIP_Bool found_clique = FALSE;
      SCIP_Real lbcolvar = local ? SCIPvarGetLbLocal(colvar) : SCIPvarGetLbGlobal(colvar);
      SCIP_Real ubcolvar = local ? SCIPvarGetUbLocal(colvar) : SCIPvarGetUbGlobal(colvar);
      SCIP_Real refpointcolvar = MAX(lbcolvar, MIN(ubcolvar, SCIPgetSolVal(scip, sol, colvar))); /*lint !e666*/

      assert(!computeEqCut);

      if( REALABS(lbcolvar) > MAXVARBOUND || REALABS(ubcolvar) > MAXVARBOUND )
      {
         *success = FALSE;
         return SCIP_OKAY;
      }

      SCIPdebugMsg(scip, "auxvar for <%s> and <%s> not found, will linearize the product\n", SCIPvarGetName(colvar), SCIPvarGetName(var));

      /* if both variables are binary, check if they are contained together in some clique */
      if( SCIPvarGetType(var) == SCIP_VARTYPE_BINARY && SCIPvarGetType(colvar) == SCIP_VARTYPE_BINARY )
      {
         int c;
         SCIP_CLIQUE** varcliques;

         varcliques = SCIPvarGetCliques(var, TRUE);

         /* look through cliques containing var */
         for( c = 0; c < SCIPvarGetNCliques(var, TRUE); ++c )
         {
            if( SCIPcliqueHasVar(varcliques[c], colvar, TRUE) ) /* var + colvar <= 1 => var*colvar = 0 */
            {
               /* product is zero, add nothing */
               found_clique = TRUE;
               break;
            }

            if( SCIPcliqueHasVar(varcliques[c], colvar, FALSE) ) /* var + (1-colvar) <= 1 => var*colvar = var */
            {
               *coefvar += coefterm;
               found_clique = TRUE;
               break;
            }
         }

         if( !found_clique )
         {
            varcliques = SCIPvarGetCliques(var, FALSE);

            /* look through cliques containing complement of var */
            for( c = 0; c < SCIPvarGetNCliques(var, FALSE); ++c )
            {
               if( SCIPcliqueHasVar(varcliques[c], colvar, TRUE) ) /* (1-var) + colvar <= 1 => var*colvar = colvar */
               {
                  coefcolvar += coefterm;
                  found_clique = TRUE;
                  break;
               }

               if( SCIPcliqueHasVar(varcliques[c], colvar, FALSE) ) /* (1-var) + (1-colvar) <= 1 => var*colvar = var + colvar - 1 */
               {
                  *coefvar += coefterm;
                  coefcolvar += coefterm;
                  *cst -= coefterm;
                  found_clique = TRUE;
                  break;
               }
            }
         }
      }

      if( !found_clique )
      {
         SCIPdebugMsg(scip, "clique for <%s> and <%s> not found or at least one of them is not binary, will use McCormick\n", SCIPvarGetName(colvar), SCIPvarGetName(var));
         SCIPaddBilinMcCormick(scip, coefterm, lbvar, ubvar, refpointvar, lbcolvar,
            ubcolvar, refpointcolvar, uselhs, coefvar, &coefcolvar, cst, success);
         if( !*success )
            return SCIP_OKAY;
      }
   }

   /* or, if it's a quadratic term, use a secant for overestimation and a gradient for underestimation */
   else
   {
      SCIPdebugMsg(scip, "auxvar for <%s>^2 not found, will use gradient and secant estimators\n", SCIPvarGetName(colvar));

      assert(!computeEqCut);

      /* for a binary var, var^2 = var */
      if( SCIPvarGetType(var) == SCIP_VARTYPE_BINARY )
      {
         *coefvar += coefterm;
      }
      else
      {
         /* depending on over-/underestimation and the sign of the column variable, compute secant or tangent */
         if( (uselhs && coefterm > 0.0) || (!uselhs && coefterm < 0.0) )
            SCIPaddSquareSecant(scip, coefterm, lbvar, ubvar, coefvar, cst, success);
         else
            SCIPaddSquareLinearization(scip, coefterm, refpointvar, SCIPvarIsIntegral(var), coefvar, cst, success);

         if( !*success )
            return SCIP_OKAY;
      }
   }

   /* add the auxiliary variable if its coefficient is nonzero */
   if( !SCIPisZero(scip, coefauxvar) )
   {
      assert(auxvar != NULL);
      /* coverity[var_deref_model] */
      SCIP_CALL( SCIPaddVarToRow(scip, cut, auxvar, coefauxvar) );
   }

   /* we are done with the product linearisation, now add the term which comes from multiplying
    * coef*colvar by the constant part of the bound factor
    */

   if( colvar != var )
   {
      assert(!SCIPisInfinity(scip, REALABS(coefcolvar)));
      SCIP_CALL( SCIPaddVarToRow(scip, cut, colvar, coefcolvar) );
   }
   else
      *coefvar += coefcolvar;

   return SCIP_OKAY;
}

/** creates the RLT cut formed by multiplying a given row with (x - lb) or (ub - x)
 *
 * In detail:
 * - The row is multiplied either with (x - lb(x)) or with (ub(x) - x), depending on parameter `uselb`, or by x if
 *   this is an equality cut
 * - The (inequality) cut is computed either for lhs or rhs, depending on parameter `uselhs`.
 * - Terms for which no auxiliary variable and no clique relation exists are replaced by either McCormick, secants,
 *   or gradient linearization cuts.
 */
static
SCIP_RETCODE computeRltCut(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_SEPA*            sepa,               /**< separator */
   SCIP_SEPADATA*        sepadata,           /**< separation data */
   SCIP_ROW**            cut,                /**< buffer to store the cut */
   SCIP_ROW*             row,                /**< the row that is used for the rlt cut (NULL if using projected row) */
   RLT_SIMPLEROW*        projrow,            /**< projected row that is used for the rlt cut (NULL if using row) */
   SCIP_SOL*             sol,                /**< the point to be separated (can be NULL) */
   int*                  bestunderest,       /**< positions of most violated underestimators for each product term */
   int*                  bestoverest,        /**< positions of most violated overestimators for each product term */
   SCIP_VAR*             var,                /**< the variable that is used for the rlt cuts */
   SCIP_Bool*            success,            /**< buffer to store whether cut was created successfully */
   SCIP_Bool             uselb,              /**< whether we multiply with (var - lb) or (ub - var) */
   SCIP_Bool             uselhs,             /**< whether to create a cut for the lhs or rhs */
   SCIP_Bool             local,              /**< whether local or global cuts should be computed */
   SCIP_Bool             computeEqCut,       /**< whether conditions are fulfilled to compute equality cuts */
   SCIP_Bool             useprojrow          /**< whether to use projected row instead of normal row */
   )
{ /*lint --e{413}*/
   SCIP_Real signfactor;
   SCIP_Real boundfactor;
   SCIP_Real lbvar;
   SCIP_Real ubvar;
   SCIP_Real coefvar;
   SCIP_Real consside;
   SCIP_Real finalside;
   SCIP_Real cstterm;
   SCIP_Real lhs;
   SCIP_Real rhs;
   SCIP_Real rowcst;
   int i;
   const char* rowname;
   char cutname[SCIP_MAXSTRLEN];

   assert(sepadata != NULL);
   assert(cut != NULL);
   assert(useprojrow || row != NULL);
   assert(!useprojrow || projrow != NULL);
   assert(var != NULL);
   assert(success != NULL);

   lhs = useprojrow ? projrow->lhs : SCIProwGetLhs(row);
   rhs = useprojrow ? projrow->rhs : SCIProwGetRhs(row);
   rowname = useprojrow ? projrow->name : SCIProwGetName(row);
   rowcst = useprojrow ? projrow ->cst : SCIProwGetConstant(row);

   assert(!computeEqCut || SCIPisEQ(scip, lhs, rhs));

   *cut = NULL;

   /* get data for given variable */
   if( computeEqCut )
   {
      lbvar = 0.0;
      ubvar = 0.0;
   }
   else
   {
      lbvar = local ? SCIPvarGetLbLocal(var) : SCIPvarGetLbGlobal(var);
      ubvar = local ? SCIPvarGetUbLocal(var) : SCIPvarGetUbGlobal(var);
   }

   /* get row side */
   consside = uselhs ? lhs : rhs;

   /* if the bounds are too large or the respective side is infinity, skip this cut */
   if( (uselb && REALABS(lbvar) > MAXVARBOUND) || (!uselb && REALABS(ubvar) > MAXVARBOUND)
         || SCIPisInfinity(scip, REALABS(consside)) )
   {
      SCIPdebugMsg(scip, "cut generation for %srow <%s>, %s, and variable <%s> with its %s %g not possible\n",
         useprojrow ? "projected " : "", rowname, uselhs ? "lhs" : "rhs", SCIPvarGetName(var),
         uselb ? "lower bound" : "upper bound", uselb ? lbvar : ubvar);

      if( REALABS(lbvar) > MAXVARBOUND )
         SCIPdebugMsg(scip, " because of lower bound\n");
      if( REALABS(ubvar) > MAXVARBOUND )
         SCIPdebugMsg(scip, " because of upper bound\n");
      if( SCIPisInfinity(scip, REALABS(consside)) )
         SCIPdebugMsg(scip, " because of side %g\n", consside);

      *success = FALSE;
      return SCIP_OKAY;
   }

   /* initialize some factors needed for computation */
   coefvar = 0.0;
   cstterm = 0.0;
   signfactor = (uselb ? 1.0 : -1.0);
   boundfactor = (uselb ? -lbvar : ubvar);
   *success = TRUE;

   /* create an empty row which we then fill with variables step by step */
   (void) SCIPsnprintf(cutname, SCIP_MAXSTRLEN, "rlt_%scut_%s_%s_%s_%s_%" SCIP_LONGINT_FORMAT, useprojrow ? "proj" : "", rowname,
                       uselhs ? "lhs" : "rhs", SCIPvarGetName(var), uselb ? "lb" : "ub", SCIPgetNLPs(scip));
   SCIP_CALL( SCIPcreateEmptyRowSepa(scip, cut, sepa, cutname, -SCIPinfinity(scip), SCIPinfinity(scip),
         SCIPgetDepth(scip) > 0 && local, FALSE, FALSE) );  /* TODO SCIPgetDepth() should be replaced by depth that is passed on to the SEPAEXEC calls (?) */

   SCIP_CALL( SCIPcacheRowExtensions(scip, *cut) );

   /* iterate over all variables in the row and add the corresponding terms coef*colvar*(bound factor) to the cuts */
   for( i = 0; i < (useprojrow ? projrow->nnonz : SCIProwGetNNonz(row)); ++i )
   {
      SCIP_VAR* colvar;

      colvar = useprojrow ? projrow->vars[i] : SCIPcolGetVar(SCIProwGetCols(row)[i]);
      SCIP_CALL( addRltTerm(scip, sepadata, sol, bestunderest, bestoverest, *cut, var, colvar,
            useprojrow ? projrow->coefs[i] : SCIProwGetVals(row)[i], uselb, uselhs, local, computeEqCut,
            &coefvar, &cstterm, success) );
   }

   if( REALABS(cstterm) > MAXVARBOUND )
   {
      *success = FALSE;
      return SCIP_OKAY;
   }

   /* multiply (x-lb) or (ub -x) with the lhs and rhs of the row */
   coefvar += signfactor * (rowcst - consside);
   finalside = boundfactor * (consside - rowcst) - cstterm;

   assert(!SCIPisInfinity(scip, REALABS(coefvar)));
   assert(!SCIPisInfinity(scip, REALABS(finalside)));

   /* set the coefficient of var and update the side */
   SCIP_CALL( SCIPaddVarToRow(scip, *cut, var, coefvar) );
   SCIP_CALL( SCIPflushRowExtensions(scip, *cut) );
   if( uselhs || computeEqCut )
   {
      SCIP_CALL( SCIPchgRowLhs(scip, *cut, finalside) );
   }
   if( !uselhs || computeEqCut )
   {
      SCIP_CALL( SCIPchgRowRhs(scip, *cut, finalside) );
   }

   SCIPdebugMsg(scip, "%scut was generated successfully:\n", useprojrow ? "projected " : "");
#ifdef SCIP_DEBUG
   SCIP_CALL( SCIPprintRow(scip, *cut, NULL) );
#endif

   return SCIP_OKAY;
}

/** store a row projected by fixing all variables that are at bound at sol; the result is a simplified row */
static
SCIP_RETCODE createProjRow(
   SCIP*                 scip,               /**< SCIP data structure */
   RLT_SIMPLEROW*        simplerow,          /**< pointer to the simplified row */
   SCIP_ROW*             row,                /**< row to be projected */
   SCIP_SOL*             sol,                /**< the point to be separated (can be NULL) */
   SCIP_Bool             local               /**< whether local bounds should be checked */
   )
{
   int i;
   SCIP_VAR* var;
   SCIP_Real val;
   SCIP_Real vlb;
   SCIP_Real vub;

   assert(simplerow != NULL);

   SCIP_CALL( SCIPduplicateBlockMemoryArray(scip,  &(simplerow->name), SCIProwGetName(row),
         strlen(SCIProwGetName(row))+1) ); /*lint !e666*/
   simplerow->nnonz = 0;
   simplerow->size = 0;
   simplerow->vars = NULL;
   simplerow->coefs = NULL;
   simplerow->lhs = SCIProwGetLhs(row);
   simplerow->rhs = SCIProwGetRhs(row);
   simplerow->cst = SCIProwGetConstant(row);

   for( i = 0; i < SCIProwGetNNonz(row); ++i )
   {
      var = SCIPcolGetVar(SCIProwGetCols(row)[i]);
      val = SCIPgetSolVal(scip, sol, var);
      vlb = local ? SCIPvarGetLbLocal(var) : SCIPvarGetLbGlobal(var);
      vub = local ? SCIPvarGetUbLocal(var) : SCIPvarGetUbGlobal(var);
      if( SCIPisFeasEQ(scip, vlb, val) || SCIPisFeasEQ(scip, vub, val) )
      {
         /* if we are projecting and the var is at bound, add var as a constant to simplerow */
         if( !SCIPisInfinity(scip, -simplerow->lhs) )
            simplerow->lhs -= SCIProwGetVals(row)[i]*val;
         if( !SCIPisInfinity(scip, simplerow->rhs) )
            simplerow->rhs -= SCIProwGetVals(row)[i]*val;
      }
      else
      {
         if( simplerow->nnonz + 1 > simplerow->size )
         {
            int newsize;

            newsize = SCIPcalcMemGrowSize(scip, simplerow->nnonz + 1);
            SCIP_CALL( SCIPreallocBufferArray(scip, &simplerow->coefs, newsize) );
            SCIP_CALL( SCIPreallocBufferArray(scip, &simplerow->vars, newsize) );
            simplerow->size = newsize;
         }

         /* add the term to simplerow */
         simplerow->vars[simplerow->nnonz] = var;
         simplerow->coefs[simplerow->nnonz] = SCIProwGetVals(row)[i];
         ++(simplerow->nnonz);
      }
   }

   return SCIP_OKAY;
}

/** free the projected row */
static
void freeProjRow(
   SCIP*                 scip,               /**< SCIP data structure */
   RLT_SIMPLEROW*        simplerow           /**< simplified row to be freed */
   )
{
   assert(simplerow != NULL);

   if( simplerow->size > 0 )
   {
      assert(simplerow->vars != NULL);
      assert(simplerow->coefs != NULL);

      SCIPfreeBufferArray(scip, &simplerow->vars);
      SCIPfreeBufferArray(scip, &simplerow->coefs);
   }
   SCIPfreeBlockMemoryArray(scip, &simplerow->name, strlen(simplerow->name)+1);
}

/** creates the projected problem
 *
 *  All variables that are at their bounds at the current solution are added
 *  to left and/or right hand sides as constant values.
 */
static
SCIP_RETCODE createProjRows(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_ROW**            rows,               /**< problem rows */
   int                   nrows,              /**< number of rows */
   SCIP_SOL*             sol,                /**< the point to be separated (can be NULL) */
   RLT_SIMPLEROW**       projrows,           /**< the projected rows to be filled */
   SCIP_Bool             local,              /**< are local cuts allowed? */
   SCIP_Bool*            allcst              /**< buffer to store whether all projected rows have only constants */
   )
{
   int i;

   assert(scip != NULL);
   assert(rows != NULL);
   assert(projrows != NULL);
   assert(allcst != NULL);

   *allcst = TRUE;
   SCIP_CALL( SCIPallocBufferArray(scip, projrows, nrows) );

   for( i = 0; i < nrows; ++i )
   {
      /* get a simplified and projected row */
      SCIP_CALL( createProjRow(scip, &(*projrows)[i], rows[i], sol, local) );
      if( (*projrows)[i].nnonz > 0 )
         *allcst = FALSE;
   }

   return SCIP_OKAY;
}

#ifdef SCIP_DEBUG
/* prints the projected LP */
static
void printProjRows(
   SCIP*                 scip,               /**< SCIP data structure */
   RLT_SIMPLEROW*        projrows,           /**< the projected rows */
   int                   nrows,              /**< number of projected rows */
   FILE*                 file                /**< output file (or NULL for standard output) */
   )
{
   int i;
   int j;

   assert(projrows != NULL);

   for( i = 0; i < nrows; ++i )
   {
      SCIPinfoMessage(scip, file, "\nproj_row[%d]: ", i);
      if( !SCIPisInfinity(scip, -projrows[i].lhs) )
         SCIPinfoMessage(scip, file, "%.15g <= ", projrows[i].lhs);
      for( j = 0; j < projrows[i].nnonz; ++j )
      {
         if( j == 0 )
         {
            if( projrows[i].coefs[j] < 0 )
               SCIPinfoMessage(scip, file, "-");
         }
         else
         {
            if( projrows[i].coefs[j] < 0 )
               SCIPinfoMessage(scip, file, " - ");
            else
               SCIPinfoMessage(scip, file, " + ");
         }

         if( projrows[i].coefs[j] != 1.0 )
            SCIPinfoMessage(scip, file, "%.15g*", REALABS(projrows[i].coefs[j]));
         SCIPinfoMessage(scip, file, "<%s>", SCIPvarGetName(projrows[i].vars[j]));
      }
      if( projrows[i].cst > 0 )
         SCIPinfoMessage(scip, file, " + %.15g", projrows[i].cst);
      else if( projrows[i].cst < 0 )
         SCIPinfoMessage(scip, file, " - %.15g", REALABS(projrows[i].cst));

      if( !SCIPisInfinity(scip, projrows[i].rhs) )
         SCIPinfoMessage(scip, file, " <= %.15g", projrows[i].rhs);
   }
   SCIPinfoMessage(scip, file, "\n");
}
#endif

/** frees the projected rows */
static
void freeProjRows(
   SCIP*                 scip,               /**< SCIP data structure */
   RLT_SIMPLEROW**       projrows,           /**< the projected LP */
   int                   nrows               /**< number of rows in projrows */
   )
{
   int i;

   for( i = 0; i < nrows; ++i )
      freeProjRow(scip, &(*projrows)[i]);

   SCIPfreeBufferArray(scip, projrows);
}

/** mark a row for rlt cut selection
 *
 * depending on the sign of the coefficient and violation, set or update mark which cut is required:
 * - 1 - cuts for axy < aw case,
 * - 2 - cuts for axy > aw case,
 * - 3 - cuts for both cases
 */
static
void addRowMark(
   int                   ridx,               /**< row index */
   SCIP_Real             a,                  /**< coefficient of x in the row */
   SCIP_Bool             violatedbelow,      /**< whether the relation auxexpr <= xy is violated */
   SCIP_Bool             violatedabove,      /**< whether the relation xy <= auxexpr is violated */
   int*                  row_idcs,           /**< sparse array with indices of marked rows */
   unsigned int*         row_marks,          /**< sparse array to store the marks */
   int*                  nmarked             /**< number of marked rows */
   )
{
   unsigned int newmark;
   int pos;
   SCIP_Bool exists;

   assert(a != 0.0);

   if( (a > 0.0 && violatedbelow) || (a < 0.0 && violatedabove) )
      newmark = 1; /* axy < aw case */
   else
      newmark = 2; /* axy > aw case */

   /* find row idx in row_idcs */
   exists = SCIPsortedvecFindInt(row_idcs, ridx, *nmarked, &pos);

   if( exists )
   {
      /* we found the row index: update the mark at pos */
      row_marks[pos] |= newmark;
   }
   else /* the given row index does not yet exist in row_idcs */
   {
      int i;

      /* insert row index at the correct position */
      for( i = *nmarked; i > pos; --i )
      {
         row_idcs[i] = row_idcs[i-1];
         row_marks[i] = row_marks[i-1];
      }
      row_idcs[pos] = ridx;
      row_marks[pos] = newmark;
      (*nmarked)++;
   }
}

/** mark all rows that should be multiplied by xj */
static
SCIP_RETCODE markRowsXj(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_SEPADATA*        sepadata,           /**< separator data */
   SCIP_CONSHDLR*        conshdlr,           /**< nonlinear constraint handler */
   SCIP_SOL*             sol,                /**< point to be separated (can be NULL) */
   int                   j,                  /**< index of the multiplier variable in sepadata */
   SCIP_Bool             local,              /**< are local cuts allowed? */
   SCIP_HASHMAP*         row_to_pos,         /**< hashmap linking row indices to positions in array */
   int*                  bestunderest,       /**< positions of most violated underestimators for each product term */
   int*                  bestoverest,        /**< positions of most violated overestimators for each product term */
   unsigned int*         row_marks,          /**< sparse array storing the row marks */
   int*                  row_idcs,           /**< sparse array storing the marked row positions */
   int*                  nmarked             /**< number of marked rows */
   )
{
   int i;
   int idx;
   int ncolrows;
   int r;
   int ridx;
   SCIP_VAR* xi;
   SCIP_VAR* xj;
   SCIP_Real vlb;
   SCIP_Real vub;
   SCIP_Real vali;
   SCIP_Real valj;
   SCIP_Real a;
   SCIP_COL* coli;
   SCIP_Real* colvals;
   SCIP_ROW** colrows;
   SCIP_CONSNONLINEAR_BILINTERM* terms;
   SCIP_Bool violatedbelow;
   SCIP_Bool violatedabove;
   SCIP_VAR** bilinadjvars;
   int nbilinadjvars;

   *nmarked = 0;

   xj = sepadata->varssorted[j];
   assert(xj != NULL);

   valj = SCIPgetSolVal(scip, sol, xj);
   vlb = local ? SCIPvarGetLbLocal(xj) : SCIPvarGetLbGlobal(xj);
   vub = local ? SCIPvarGetUbLocal(xj) : SCIPvarGetUbGlobal(xj);

   if( sepadata->useprojection && (SCIPisFeasEQ(scip, vlb, valj) || SCIPisFeasEQ(scip, vub, valj)) )
   {
      /* we don't want to multiply by variables that are at bound */
      SCIPdebugMsg(scip, "Rejected multiplier <%s> in [%g,%g] because it is at bound (current value %g)\n", SCIPvarGetName(xj), vlb, vub, valj);
      return SCIP_OKAY;
   }

   terms = SCIPgetBilinTermsNonlinear(conshdlr);
   bilinadjvars = getAdjacentVars(sepadata->bilinvardatamap, xj, &nbilinadjvars);
   assert(bilinadjvars != NULL);

   /* for each var which appears in a bilinear product together with xj, mark rows */
   for( i = 0; i < nbilinadjvars; ++i )
   {
      xi = bilinadjvars[i];

      if( SCIPvarGetStatus(xi) != SCIP_VARSTATUS_COLUMN )
         continue;

      vali = SCIPgetSolVal(scip, sol, xi);
      vlb = local ? SCIPvarGetLbLocal(xi) : SCIPvarGetLbGlobal(xi);
      vub = local ? SCIPvarGetUbLocal(xi) : SCIPvarGetUbGlobal(xi);

      /* if we use projection, we aren't interested in products with variables that are at bound */
      if( sepadata->useprojection && (SCIPisFeasEQ(scip, vlb, vali) || SCIPisFeasEQ(scip, vub, vali)) )
         continue;

      /* get the index of the bilinear product */
      idx = SCIPgetBilinTermIdxNonlinear(conshdlr, xj, xi);
      assert(idx >= 0 && idx < SCIPgetNBilinTermsNonlinear(conshdlr));

      /* skip implicit products if we don't want to add RLT cuts for them */
      if( !sepadata->hiddenrlt && !terms[idx].existing )
         continue;

      /* use the most violated under- and overestimators for this product;
       * if equality cuts are computed, we might end up using a different auxiliary expression;
       * so this is an optimistic (i.e. taking the largest possible violation) estimation
       */
      if( bestunderest == NULL || bestunderest[idx] == -1 )
      { /* no violated implicit underestimation relations -> either use auxvar or set violatedbelow to FALSE */
         if( terms[idx].nauxexprs == 0 && terms[idx].aux.var != NULL )
         {
            assert(terms[idx].existing);
            violatedbelow = SCIPisFeasPositive(scip, SCIPgetSolVal(scip, sol, terms[idx].aux.var) - valj * vali);
         }
         else
         {
            assert(bestunderest != NULL);
            violatedbelow = FALSE;
         }
      }
      else
      {
         assert(bestunderest[idx] >= 0 && bestunderest[idx] < terms[idx].nauxexprs);

         /* if we are here, the relation with the best underestimator must be violated */
         assert(SCIPisFeasPositive(scip, SCIPevalBilinAuxExprNonlinear(scip, terms[idx].x, terms[idx].y,
               terms[idx].aux.exprs[bestunderest[idx]], sol) - valj * vali));
         violatedbelow = TRUE;
      }

      if( bestoverest == NULL || bestoverest[idx] == -1 )
      { /* no violated implicit overestimation relations -> either use auxvar or set violatedabove to FALSE */
         if( terms[idx].nauxexprs == 0 && terms[idx].aux.var != NULL )
         {
            assert(terms[idx].existing);
            violatedabove = SCIPisFeasPositive(scip, valj * vali - SCIPgetSolVal(scip, sol, terms[idx].aux.var));
         }
         else
         {
            assert(bestoverest != NULL);
            violatedabove = FALSE;
         }
      }
      else
      {
         assert(bestoverest[idx] >= 0 && bestoverest[idx] < terms[idx].nauxexprs);

         /* if we are here, the relation with the best overestimator must be violated */
         assert(SCIPisFeasPositive(scip, valj * vali - SCIPevalBilinAuxExprNonlinear(scip, terms[idx].x, terms[idx].y,
               terms[idx].aux.exprs[bestoverest[idx]], sol)));
         violatedabove = TRUE;
      }

      /* only violated products contribute to row marks */
      if( !violatedbelow && !violatedabove )
      {
         SCIPdebugMsg(scip, "the product for vars <%s> and <%s> is not violated\n", SCIPvarGetName(xj), SCIPvarGetName(xi));
         continue;
      }

      /* get the column of xi */
      coli = SCIPvarGetCol(xi);
      colvals = SCIPcolGetVals(coli);
      ncolrows = SCIPcolGetNNonz(coli);
      colrows = SCIPcolGetRows(coli);

      SCIPdebugMsg(scip, "marking rows for xj = <%s>, xi = <%s>\n", SCIPvarGetName(xj), SCIPvarGetName(xi));

      /* mark the rows */
      for( r = 0; r < ncolrows; ++r )
      {
         ridx = SCIProwGetIndex(colrows[r]);

         if( !SCIPhashmapExists(row_to_pos, (void*)(size_t)ridx) )
            continue; /* if row index is not in row_to_pos, it means that storeSuitableRows decided to ignore this row */

         a = colvals[r];
         if( a == 0.0 )
            continue;

         SCIPdebugMsg(scip, "Marking row %d\n", ridx);
         addRowMark(ridx, a, violatedbelow, violatedabove, row_idcs, row_marks, nmarked);
      }
   }

   return SCIP_OKAY;
}

/** adds McCormick inequalities for implicit products */
static
SCIP_RETCODE separateMcCormickImplicit(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_SEPA*            sepa,               /**< separator */
   SCIP_SEPADATA*        sepadata,           /**< separator data */
   SCIP_SOL*             sol,                /**< the point to be separated (can be NULL) */
   int*                  bestunderestimators,/**< indices of auxiliary underestimators with largest violation in sol */
   int*                  bestoverestimators, /**< indices of auxiliary overestimators with largest violation in sol */
   SCIP_RESULT*          result              /**< pointer to store the result */
   )
{
   int i;
   int j;
   SCIP_CONSNONLINEAR_BILINTERM* terms;
   SCIP_ROW* cut;
   char name[SCIP_MAXSTRLEN];
   SCIP_Bool underestimate;
   SCIP_Real xcoef;
   SCIP_Real ycoef;
   SCIP_Real auxcoef;
   SCIP_Real constant;
   SCIP_Bool success;
   SCIP_CONSNONLINEAR_AUXEXPR* auxexpr;
   SCIP_Bool cutoff;
   SCIP_Real refpointx;
   SCIP_Real refpointy;
   SCIP_INTERVAL bndx;
   SCIP_INTERVAL bndy;
#ifndef NDEBUG
   SCIP_Real productval;
   SCIP_Real auxval;
#endif

   assert(sepadata->nbilinterms == SCIPgetNBilinTermsNonlinear(sepadata->conshdlr));
   assert(bestunderestimators != NULL && bestoverestimators != NULL);

   cutoff = FALSE;
   terms = SCIPgetBilinTermsNonlinear(sepadata->conshdlr);

   for( i = 0; i < sepadata->nbilinterms; ++i )
   {
      if( terms[i].existing )
         continue;

      assert(terms[i].nauxexprs > 0);

      bndx.inf = SCIPvarGetLbLocal(terms[i].x);
      bndx.sup = SCIPvarGetUbLocal(terms[i].x);
      bndy.inf = SCIPvarGetLbLocal(terms[i].y);
      bndy.sup = SCIPvarGetUbLocal(terms[i].y);
      refpointx = SCIPgetSolVal(scip, sol, terms[i].x);
      refpointy = SCIPgetSolVal(scip, sol, terms[i].y);

      /* adjust the reference points */
      refpointx = MIN(MAX(refpointx, bndx.inf), bndx.sup); /*lint !e666*/
      refpointy = MIN(MAX(refpointy, bndy.inf), bndy.sup); /*lint !e666*/

      /* one iteration for underestimation and one for overestimation */
      for( j = 0; j < 2; ++j )
      {
         /* if underestimate, separate xy <= auxexpr; if !underestimate, separate xy >= auxexpr;
          * the cuts will be:
          * if underestimate: McCormick_under(xy) - auxexpr <= 0,
          * if !underestimate: McCormick_over(xy) - auxexpr >= 0
          */
         underestimate = j == 0;
         if( underestimate && bestoverestimators[i] != -1 )
            auxexpr = terms[i].aux.exprs[bestoverestimators[i]];
         else if( !underestimate && bestunderestimators[i] != -1 )
            auxexpr = terms[i].aux.exprs[bestunderestimators[i]];
         else
            continue;
         assert(!underestimate || auxexpr->overestimate);
         assert(underestimate || auxexpr->underestimate);

#ifndef NDEBUG
         /* make sure that the term is violated */
         productval = SCIPgetSolVal(scip, sol, terms[i].x) * SCIPgetSolVal(scip, sol, terms[i].y);
         auxval = SCIPevalBilinAuxExprNonlinear(scip, terms[i].x, terms[i].y, auxexpr, sol);

         /* if underestimate, then xy <= aux must be violated; otherwise aux <= xy must be violated */
         assert((underestimate && SCIPisFeasLT(scip, auxval, productval)) ||
               (!underestimate && SCIPisFeasLT(scip, productval, auxval)));
#endif

         /* create an empty row */
         (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "mccormick_%sestimate_implicit_%s*%s_%" SCIP_LONGINT_FORMAT,
                             underestimate ? "under" : "over", SCIPvarGetName(terms[i].x), SCIPvarGetName(terms[i].y),
                             SCIPgetNLPs(scip));

         SCIP_CALL(SCIPcreateEmptyRowSepa(scip, &cut, sepa, name, -SCIPinfinity(scip), SCIPinfinity(scip), TRUE,
               FALSE, FALSE) );

         xcoef = 0.0;
         ycoef = 0.0;
         auxcoef = 0.0;
         constant = 0.0;
         success = TRUE;

         /* subtract auxexpr from the cut */
         addAuxexprCoefs(terms[i].x, terms[i].y, auxexpr, -1.0, &auxcoef, &xcoef, &ycoef, &constant);

         /* add McCormick terms: ask for an underestimator if relation is xy <= auxexpr, and vice versa */
         SCIPaddBilinMcCormick(scip, 1.0, bndx.inf, bndx.sup, refpointx, bndy.inf, bndy.sup, refpointy, !underestimate,
               &xcoef, &ycoef, &constant, &success);

         if( REALABS(constant) > MAXVARBOUND )
            success = FALSE;

         if( success )
         {
            assert(!SCIPisInfinity(scip, REALABS(xcoef)));
            assert(!SCIPisInfinity(scip, REALABS(ycoef)));
            assert(!SCIPisInfinity(scip, REALABS(constant)));

            SCIP_CALL( SCIPaddVarToRow(scip, cut, terms[i].x, xcoef) );
            SCIP_CALL( SCIPaddVarToRow(scip, cut, terms[i].y, ycoef) );
            SCIP_CALL( SCIPaddVarToRow(scip, cut, auxexpr->auxvar, auxcoef) );

            /* set side */
            if( underestimate )
               SCIP_CALL( SCIPchgRowRhs(scip, cut, -constant) );
            else
               SCIP_CALL( SCIPchgRowLhs(scip, cut, -constant) );

            /* if the cut is violated, add it to SCIP */
            if( SCIPisFeasNegative(scip, SCIPgetRowFeasibility(scip, cut)) )
            {
               SCIP_CALL( SCIPaddRow(scip, cut, FALSE, &cutoff) );
               *result = SCIP_SEPARATED;
            }
            else
            {
               SCIPdebugMsg(scip, "\nMcCormick cut for hidden product <%s>*<%s> was created successfully, but is not violated",
                            SCIPvarGetName(terms[i].x), SCIPvarGetName(terms[i].y));
            }
         }

         /* release the cut */
         if( cut != NULL )
         {
            SCIP_CALL( SCIPreleaseRow(scip, &cut) );
         }

         if( cutoff )
         {
            *result = SCIP_CUTOFF;
            SCIPdebugMsg(scip, "exit separator because we found a cutoff -> skip\n");
            return SCIP_OKAY;
         }
      }
   }

   return SCIP_OKAY;
}

/** builds and adds the RLT cuts */
static
SCIP_RETCODE separateRltCuts(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_SEPA*            sepa,               /**< separator */
   SCIP_SEPADATA*        sepadata,           /**< separator data */
   SCIP_CONSHDLR*        conshdlr,           /**< nonlinear constraint handler */
   SCIP_SOL*             sol,                /**< the point to be separated (can be NULL) */
   SCIP_HASHMAP*         row_to_pos,         /**< hashmap linking row indices to positions in array */
   RLT_SIMPLEROW*        projrows,           /**< projected rows */
   SCIP_ROW**            rows,               /**< problem rows */
   int                   nrows,              /**< number of problem rows */
   SCIP_Bool             allowlocal,         /**< are local cuts allowed? */
   int*                  bestunderestimators,/**< indices of auxiliary underestimators with largest violation in sol */
   int*                  bestoverestimators, /**< indices of auxiliary overestimators with largest violation in sol */
   SCIP_RESULT*          result              /**< buffer to store whether separation was successful */
   )
{
   int j;
   int r;
   int k;
   int nmarked;
   int cutssize;
   int ncuts;
   SCIP_VAR* xj;
   unsigned int* row_marks;
   int* row_idcs;
   SCIP_ROW* cut;
   SCIP_ROW** cuts;
   SCIP_Bool uselb[4] = {TRUE, TRUE, FALSE, FALSE};
   SCIP_Bool uselhs[4] = {TRUE, FALSE, TRUE, FALSE};
   SCIP_Bool success;
   SCIP_Bool infeasible;
   SCIP_Bool accepted;
   SCIP_Bool buildeqcut;
   SCIP_Bool iseqrow;

   assert(!sepadata->useprojection || projrows != NULL);
   assert(!sepadata->detecthidden || (bestunderestimators != NULL && bestoverestimators != NULL));

   ncuts = 0;
   cutssize = 0;
   cuts = NULL;
   *result = SCIP_DIDNOTFIND;

   SCIP_CALL( SCIPallocCleanBufferArray(scip, &row_marks, nrows) );
   SCIP_CALL( SCIPallocBufferArray(scip, &row_idcs, nrows) );

   /* loop through all variables that appear in bilinear products */
   for( j = 0; j < sepadata->nbilinvars && (sepadata->maxusedvars < 0 || j < sepadata->maxusedvars); ++j )
   {
      xj = sepadata->varssorted[j];

      /* mark all rows for multiplier xj */
      SCIP_CALL( markRowsXj(scip, sepadata, conshdlr, sol, j, allowlocal, row_to_pos, bestunderestimators,
         bestoverestimators, row_marks, row_idcs, &nmarked) );

      assert(nmarked <= nrows);

      /* generate the projected cut and if it is violated, generate the actual cut */
      for( r = 0; r < nmarked; ++r )
      {
         int pos;
         int currentnunknown;
         SCIP_ROW* row;

         assert(row_marks[r] != 0);
         assert(SCIPhashmapExists(row_to_pos, (void*)(size_t) row_idcs[r])); /*lint !e571 */

         pos = SCIPhashmapGetImageInt(row_to_pos, (void*)(size_t) row_idcs[r]); /*lint !e571 */
         row = rows[pos];
         assert(SCIProwGetIndex(row) == row_idcs[r]);

         /* check whether this row and var fulfill the conditions */
         SCIP_CALL( isAcceptableRow(sepadata, row, xj, &currentnunknown, &accepted) );
         if( !accepted )
         {
            SCIPdebugMsg(scip, "rejected row <%s> for variable <%s> (introduces too many new products)\n", SCIProwGetName(row), SCIPvarGetName(xj));
            row_marks[r] = 0;
            continue;
         }

         SCIPdebugMsg(scip, "accepted row <%s> for variable <%s>\n", SCIProwGetName(rows[r]), SCIPvarGetName(xj));
#ifdef SCIP_DEBUG
         SCIP_CALL( SCIPprintRow(scip, rows[r], NULL) );
#endif
         iseqrow = SCIPisEQ(scip, SCIProwGetLhs(row), SCIProwGetRhs(row));

         /* if all terms are known and it is an equality row, compute equality cut, that is, multiply row with (x-lb) and/or (ub-x) (but see also @todo at top)
          * otherwise, multiply row w.r.t. lhs and/or rhs with (x-lb) and/or (ub-x) and estimate product terms that have no aux.var or aux.expr
          */
         buildeqcut = (currentnunknown == 0 && iseqrow);

         /* go over all suitable combinations of sides and bounds and compute the respective cuts */
         for( k = 0; k < 4; ++k )
         {
            /* if equality cuts are possible, lhs and rhs cuts are equal so skip rhs */
            if( buildeqcut )
            {
               if( k != 1 )
                  continue;
            }
            /* otherwise which cuts are generated depends on the marks */
            else
            {
               if( row_marks[r] == 1 && uselb[k] == uselhs[k] )
                  continue;

               if( row_marks[r] == 2 && uselb[k] != uselhs[k] )
                  continue;
            }

            success = TRUE;
            cut = NULL;

            SCIPdebugMsg(scip, "row <%s>, uselb = %u, uselhs = %u\n", SCIProwGetName(row), uselb[k], uselhs[k]);

            if( sepadata->useprojection )
            {
               /* if no variables are left in the projected row, the RLT cut will not be violated */
               if( projrows[pos].nnonz == 0 )
                  continue;

               /* compute the rlt cut for a projected row first */
               SCIP_CALL( computeRltCut(scip, sepa, sepadata, &cut, NULL, &(projrows[pos]), sol, bestunderestimators,
                     bestoverestimators, xj, &success, uselb[k], uselhs[k], allowlocal, buildeqcut, TRUE) );

               /* if the projected cut is not violated, set success to FALSE */
               if( cut != NULL )
               {
                  SCIPdebugMsg(scip, "proj cut viol = %g\n", -SCIPgetRowFeasibility(scip, cut));
               }
               if( cut != NULL && !SCIPisFeasPositive(scip, -SCIPgetRowFeasibility(scip, cut)) )
               {
                  SCIPdebugMsg(scip, "projected cut is not violated, feasibility = %g\n", SCIPgetRowFeasibility(scip, cut));
                  success = FALSE;
               }

               /* release the projected cut */
               if( cut != NULL )
                  SCIP_CALL( SCIPreleaseRow(scip, &cut) );
            }

            /* if we don't use projection or if the projected cut was generated successfully and is violated,
             * generate the actual cut */
            if( success )
            {
               SCIP_CALL( computeRltCut(scip, sepa, sepadata, &cut, row, NULL, sol, bestunderestimators,
                     bestoverestimators, xj, &success, uselb[k], uselhs[k], allowlocal, buildeqcut, FALSE) );
            }

            if( success )
            {
               success = SCIPisFeasNegative(scip, SCIPgetRowFeasibility(scip, cut)) || (sepadata->addtopool &&
                       !SCIProwIsLocal(cut));
            }

            /* if the cut was created successfully and is violated or (if addtopool == TRUE) globally valid,
             * it is added to the cuts array */
            if( success )
            {
               if( ncuts + 1 > cutssize )
               {
                  int newsize;

                  newsize = SCIPcalcMemGrowSize(scip, ncuts + 1);
                  SCIP_CALL( SCIPreallocBufferArray(scip, &cuts, newsize) );
                  cutssize = newsize;
               }
               cuts[ncuts] = cut;
               (ncuts)++;
            }
            else
            {
               SCIPdebugMsg(scip, "the generation of the cut failed or cut not violated and not added to cutpool\n");
               /* release the cut */
               if( cut != NULL )
               {
                  SCIP_CALL( SCIPreleaseRow(scip, &cut) );
               }
            }
         }

         /* clear row_marks[r] since it will be used for the next multiplier */
         row_marks[r] = 0;
      }
   }

   /* if cuts were found, we apply an additional filtering procedure, which is similar to sepastore */
   if( ncuts > 0  )
   {
      int nselectedcuts;
      int i;

      assert(cuts != NULL);

      SCIP_CALL( SCIPselectCutsHybrid(scip, cuts, NULL, NULL, sepadata->goodscore, sepadata->badscore, sepadata->goodmaxparall,
            sepadata->maxparall, sepadata->dircutoffdistweight, sepadata->efficacyweight, sepadata->objparalweight,
            0.0, ncuts, 0, sepadata->maxncuts == -1 ? ncuts : sepadata->maxncuts, &nselectedcuts) );

      for( i = 0; i < ncuts; ++i )
      {
         assert(cuts[i] != NULL);

         if( i < nselectedcuts )
         {
            /* if selected, add global cuts to the pool and local cuts to the sepastore */
            if( SCIProwIsLocal(cuts[i]) || !sepadata->addtopool )
            {
               SCIP_CALL( SCIPaddRow(scip, cuts[i], FALSE, &infeasible) );

               if( infeasible )
               {
                  SCIPdebugMsg(scip, "CUTOFF! The cut <%s> revealed infeasibility\n", SCIProwGetName(cuts[i]));
                  *result = SCIP_CUTOFF;
               }
               else
               {
                  SCIPdebugMsg(scip, "SEPARATED: added cut to scip\n");
                  *result = SCIP_SEPARATED;
               }
            }
            else
            {
               SCIP_CALL( SCIPaddPoolCut(scip, cuts[i]) );
            }
         }

         /* release current cut */
         SCIP_CALL( SCIPreleaseRow(scip, &cuts[i]) );
      }
   }

   SCIPdebugMsg(scip, "exit separator because cut calculation is finished\n");

   SCIPfreeBufferArrayNull(scip, &cuts);
   SCIPfreeBufferArray(scip, &row_idcs);
   SCIPfreeCleanBufferArray(scip, &row_marks);

   return SCIP_OKAY;
}

/*
 * Callback methods of separator
 */

/** copy method for separator plugins (called when SCIP copies plugins) */
static
SCIP_DECL_SEPACOPY(sepaCopyRlt)
{  /*lint --e{715}*/
   assert(scip != NULL);
   assert(sepa != NULL);
   assert(strcmp(SCIPsepaGetName(sepa), SEPA_NAME) == 0);

   /* call inclusion method of separator */
   SCIP_CALL( SCIPincludeSepaRlt(scip) );

   return SCIP_OKAY;
}

/** destructor of separator to free user data (called when SCIP is exiting) */
static
SCIP_DECL_SEPAFREE(sepaFreeRlt)
{  /*lint --e{715}*/
   SCIP_SEPADATA* sepadata;

   assert(strcmp(SCIPsepaGetName(sepa), SEPA_NAME) == 0);

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

   /* free separator data */
   SCIPfreeBlockMemory(scip, &sepadata);

   SCIPsepaSetData(sepa, NULL);

   return SCIP_OKAY;
}

/** solving process deinitialization method of separator (called before branch and bound process data is freed) */
static
SCIP_DECL_SEPAEXITSOL(sepaExitsolRlt)
{  /*lint --e{715}*/
   SCIP_SEPADATA* sepadata;

   assert(strcmp(SCIPsepaGetName(sepa), SEPA_NAME) == 0);

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

   if( sepadata->iscreated )
   {
      SCIP_CALL( freeSepaData(scip, sepadata) );
   }

   return SCIP_OKAY;
}

/** LP solution separation method of separator */
static
SCIP_DECL_SEPAEXECLP(sepaExeclpRlt)
{  /*lint --e{715}*/
   SCIP_ROW** prob_rows;
   SCIP_ROW** rows;
   SCIP_SEPADATA* sepadata;
   int ncalls;
   int nrows;
   SCIP_HASHMAP* row_to_pos;
   RLT_SIMPLEROW* projrows;

   assert(strcmp(SCIPsepaGetName(sepa), SEPA_NAME) == 0);

   sepadata = SCIPsepaGetData(sepa);
   *result = SCIP_DIDNOTRUN;

   if( sepadata->maxncuts == 0 )
   {
      SCIPdebugMsg(scip, "exit separator because maxncuts is set to 0\n");
      return SCIP_OKAY;
   }

   /* don't run in a sub-SCIP or in probing */
   if( SCIPgetSubscipDepth(scip) > 0 && !sepadata->useinsubscip )
   {
      SCIPdebugMsg(scip, "exit separator because in sub-SCIP\n");
      return SCIP_OKAY;
   }

   /* don't run in probing */
   if( SCIPinProbing(scip) )
   {
      SCIPdebugMsg(scip, "exit separator because in probing\n");
      return SCIP_OKAY;
   }

   /* only call separator a given number of times at each node */
   ncalls = SCIPsepaGetNCallsAtNode(sepa);
   if( (depth == 0 && sepadata->maxroundsroot >= 0 && ncalls >= sepadata->maxroundsroot)
        || (depth > 0 && sepadata->maxrounds >= 0 && ncalls >= sepadata->maxrounds) )
   {
      SCIPdebugMsg(scip, "exit separator because round limit for this node is reached\n");
      return SCIP_OKAY;
   }

   /* if this is called for the first time, create the sepadata and start the initial separation round */
   if( !sepadata->iscreated )
   {
      *result = SCIP_DIDNOTFIND;
      SCIP_CALL( createSepaData(scip, sepadata) );
   }
   assert(sepadata->iscreated || (sepadata->nbilinvars == 0 && sepadata->nbilinterms == 0));
   assert(sepadata->nbilinterms == SCIPgetNBilinTermsNonlinear(sepadata->conshdlr));

   /* no bilinear terms available -> skip */
   if( sepadata->nbilinvars == 0 )
   {
      SCIPdebugMsg(scip, "exit separator because there are no known bilinear terms\n");
      return SCIP_OKAY;
   }

   /* only call separator, if we are not close to terminating */
   if( SCIPisStopped(scip) )
   {
      SCIPdebugMsg(scip, "exit separator because we are too close to terminating\n");
      return SCIP_OKAY;
   }

   /* only call separator, if an optimal LP solution is at hand */
   if( SCIPgetLPSolstat(scip) != SCIP_LPSOLSTAT_OPTIMAL )
   {
      SCIPdebugMsg(scip, "exit separator because there is no LP solution at hand\n");
      return SCIP_OKAY;
   }

   /* get the rows, depending on settings */
   if( sepadata->isinitialround || sepadata->onlyoriginal )
   {
      SCIP_CALL( getOriginalRows(scip, &prob_rows, &nrows) );
   }
   else
   {
      SCIP_CALL( SCIPgetLPRowsData(scip, &prob_rows, &nrows) );
   }

   /* save the suitable rows */
   SCIP_CALL( SCIPallocBufferArray(scip, &rows, nrows) );
   SCIP_CALL( SCIPhashmapCreate(&row_to_pos, SCIPblkmem(scip), nrows) );

   SCIP_CALL( storeSuitableRows(scip, sepa, sepadata, prob_rows, rows, &nrows, row_to_pos, allowlocal) );

   if( nrows == 0 ) /* no suitable rows found, free memory and exit */
   {
      SCIPhashmapFree(&row_to_pos);
      SCIPfreeBufferArray(scip, &rows);
      if( sepadata->isinitialround || sepadata->onlyoriginal )
      {
         SCIPfreeBufferArray(scip, &prob_rows);
         sepadata->isinitialround = FALSE;
      }
      return SCIP_OKAY;
   }

   /* create the projected problem */
   if( sepadata->useprojection )
   {
      SCIP_Bool allcst;

      SCIP_CALL( createProjRows(scip, rows, nrows, NULL, &projrows, allowlocal, &allcst) );

      /* if all projected rows have only constants left, quit */
      if( allcst )
         goto TERMINATE;

#ifdef SCIP_DEBUG
      printProjRows(scip, projrows, nrows, NULL);
#endif
   }
   else
   {
      projrows = NULL;
   }

   /* separate the cuts */
   if( sepadata->detecthidden )
   {
      int* bestunderestimators;
      int* bestoverestimators;

      /* if we detect implicit products, a term might have more than one estimator in each direction;
       * save the indices of the most violated estimators
       */
      SCIP_CALL( SCIPallocBufferArray(scip, &bestunderestimators, sepadata->nbilinterms) );
      SCIP_CALL( SCIPallocBufferArray(scip, &bestoverestimators, sepadata->nbilinterms) );
      getBestEstimators(scip, sepadata, NULL, bestunderestimators, bestoverestimators);

      /* also separate McCormick cuts for implicit products */
      SCIP_CALL( separateMcCormickImplicit(scip, sepa, sepadata, NULL, bestunderestimators, bestoverestimators,
            result) );

      if( *result != SCIP_CUTOFF )
      {
         SCIP_CALL( separateRltCuts(scip, sepa, sepadata, sepadata->conshdlr, NULL, row_to_pos, projrows, rows, nrows,
               allowlocal, bestunderestimators, bestoverestimators, result) );
      }

      SCIPfreeBufferArray(scip, &bestoverestimators);
      SCIPfreeBufferArray(scip, &bestunderestimators);
   }
   else
   {
      SCIP_CALL( separateRltCuts(scip, sepa, sepadata, sepadata->conshdlr, NULL, row_to_pos, projrows, rows, nrows,
            allowlocal, NULL, NULL, result) );
   }

 TERMINATE:
   /* free the projected problem */
   if( sepadata->useprojection )
   {
      freeProjRows(scip, &projrows, nrows);
   }

   SCIPhashmapFree(&row_to_pos);
   SCIPfreeBufferArray(scip, &rows);

   if( sepadata->isinitialround || sepadata->onlyoriginal )
   {
      SCIPfreeBufferArray(scip, &prob_rows);
      sepadata->isinitialround = FALSE;
   }

   return SCIP_OKAY;
}

/*
 * separator specific interface methods
 */

/** creates the RLT separator and includes it in SCIP */
SCIP_RETCODE SCIPincludeSepaRlt(
   SCIP*                 scip                /**< SCIP data structure */
   )
{
   SCIP_SEPADATA* sepadata;
   SCIP_SEPA* sepa;

   /* create RLT separator data */
   SCIP_CALL( SCIPallocClearBlockMemory(scip, &sepadata) );
   sepadata->conshdlr = SCIPfindConshdlr(scip, "nonlinear");
   assert(sepadata->conshdlr != NULL);

   /* include separator */
   SCIP_CALL( SCIPincludeSepaBasic(scip, &sepa, SEPA_NAME, SEPA_DESC, SEPA_PRIORITY, SEPA_FREQ, SEPA_MAXBOUNDDIST,
         SEPA_USESSUBSCIP, SEPA_DELAY, sepaExeclpRlt, NULL, sepadata) );

   /* set non fundamental callbacks via setter functions */
   SCIP_CALL( SCIPsetSepaCopy(scip, sepa, sepaCopyRlt) );
   SCIP_CALL( SCIPsetSepaFree(scip, sepa, sepaFreeRlt) );
   SCIP_CALL( SCIPsetSepaExitsol(scip, sepa, sepaExitsolRlt) );

   /* add RLT separator parameters */
   SCIP_CALL( SCIPaddIntParam(scip,
         "separating/" SEPA_NAME "/maxncuts",
         "maximal number of rlt-cuts that are added per round (-1: unlimited)",
         &sepadata->maxncuts, FALSE, DEFAULT_MAXNCUTS, -1, INT_MAX, NULL, NULL) );

   SCIP_CALL( SCIPaddIntParam(scip,
         "separating/" SEPA_NAME "/maxunknownterms",
         "maximal number of unknown bilinear terms a row is still used with (-1: unlimited)",
         &sepadata->maxunknownterms, FALSE, DEFAULT_MAXUNKNOWNTERMS, -1, INT_MAX, NULL, NULL) );

   SCIP_CALL( SCIPaddIntParam(scip,
         "separating/" SEPA_NAME "/maxusedvars",
         "maximal number of variables used to compute rlt cuts (-1: unlimited)",
         &sepadata->maxusedvars, FALSE, DEFAULT_MAXUSEDVARS, -1, INT_MAX, NULL, NULL) );

   SCIP_CALL( SCIPaddIntParam(scip,
      "separating/" SEPA_NAME "/maxrounds",
      "maximal number of separation rounds per node (-1: unlimited)",
      &sepadata->maxrounds, FALSE, DEFAULT_MAXROUNDS, -1, INT_MAX, NULL, NULL) );

   SCIP_CALL( SCIPaddIntParam(scip,
      "separating/" SEPA_NAME "/maxroundsroot",
      "maximal number of separation rounds in the root node (-1: unlimited)",
      &sepadata->maxroundsroot, FALSE, DEFAULT_MAXROUNDSROOT, -1, INT_MAX, NULL, NULL) );

   SCIP_CALL( SCIPaddBoolParam(scip,
      "separating/" SEPA_NAME "/onlyeqrows",
      "if set to true, only equality rows are used for rlt cuts",
      &sepadata->onlyeqrows, FALSE, DEFAULT_ONLYEQROWS, NULL, NULL) );

   SCIP_CALL( SCIPaddBoolParam(scip,
      "separating/" SEPA_NAME "/onlycontrows",
      "if set to true, only continuous rows are used for rlt cuts",
      &sepadata->onlycontrows, FALSE, DEFAULT_ONLYCONTROWS, NULL, NULL) );

   SCIP_CALL( SCIPaddBoolParam(scip,
      "separating/" SEPA_NAME "/onlyoriginal",
      "if set to true, only original rows and variables are used",
      &sepadata->onlyoriginal, FALSE, DEFAULT_ONLYORIGINAL, NULL, NULL) );

   SCIP_CALL( SCIPaddBoolParam(scip,
      "separating/" SEPA_NAME "/useinsubscip",
      "if set to true, rlt is also used in sub-scips",
      &sepadata->useinsubscip, FALSE, DEFAULT_USEINSUBSCIP, NULL, NULL) );

   SCIP_CALL( SCIPaddBoolParam(scip,
                               "separating/" SEPA_NAME "/useprojection",
      "if set to true, projected rows are checked first",
      &sepadata->useprojection, FALSE, DEFAULT_USEPROJECTION, NULL, NULL) );

   SCIP_CALL( SCIPaddBoolParam(scip,
                               "separating/" SEPA_NAME "/detecthidden",
      "if set to true, hidden products are detected and separated by McCormick cuts",
      &sepadata->detecthidden, FALSE, DEFAULT_DETECTHIDDEN, NULL, NULL) );

   SCIP_CALL( SCIPaddBoolParam(scip,
                               "separating/" SEPA_NAME "/hiddenrlt",
      "whether RLT cuts (TRUE) or only McCormick inequalities (FALSE) should be added for hidden products",
      &sepadata->hiddenrlt, FALSE, DEFAULT_HIDDENRLT, NULL, NULL) );

   SCIP_CALL( SCIPaddBoolParam(scip,
                               "separating/" SEPA_NAME "/addtopool",
      "if set to true, globally valid RLT cuts are added to the global cut pool",
      &sepadata->addtopool, FALSE, DEFAULT_ADDTOPOOL, NULL, NULL) );

   SCIP_CALL( SCIPaddRealParam(scip,
                               "separating/" SEPA_NAME "/goodscore",
      "threshold for score of cut relative to best score to be considered good, so that less strict filtering is applied",
      &sepadata->goodscore, TRUE, DEFAULT_GOODSCORE, 0.0, 1.0, NULL, NULL) );

   SCIP_CALL( SCIPaddRealParam(scip,
                               "separating/" SEPA_NAME "/badscore",
      "threshold for score of cut relative to best score to be discarded",
      &sepadata->badscore, TRUE, DEFAULT_BADSCORE, 0.0, 1.0, NULL, NULL) );

   SCIP_CALL( SCIPaddRealParam(scip,
                               "separating/" SEPA_NAME "/objparalweight",
      "weight of objective parallelism in cut score calculation",
      &sepadata->objparalweight, TRUE, DEFAULT_OBJPARALWEIGHT, 0.0, 1.0, NULL, NULL) );

   SCIP_CALL( SCIPaddRealParam(scip,
                               "separating/" SEPA_NAME "/efficacyweight",
      "weight of efficacy in cut score calculation",
      &sepadata->efficacyweight, TRUE, DEFAULT_EFFICACYWEIGHT, 0.0, 1.0, NULL, NULL) );

   SCIP_CALL( SCIPaddRealParam(scip,
                               "separating/" SEPA_NAME "/dircutoffdistweight",
      "weight of directed cutoff distance in cut score calculation",
      &sepadata->dircutoffdistweight, TRUE, DEFAULT_DIRCUTOFFDISTWEIGHT, 0.0, 1.0, NULL, NULL) );

   SCIP_CALL( SCIPaddRealParam(scip,
                               "separating/" SEPA_NAME "/goodmaxparall",
      "maximum parallelism for good cuts",
      &sepadata->goodmaxparall, TRUE, DEFAULT_GOODMAXPARALL, 0.0, 1.0, NULL, NULL) );

   SCIP_CALL( SCIPaddRealParam(scip,
                               "separating/" SEPA_NAME "/maxparall",
      "maximum parallelism for non-good cuts",
      &sepadata->maxparall, TRUE, DEFAULT_MAXPARALL, 0.0, 1.0, NULL, NULL) );

   return SCIP_OKAY;
}