/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
	/*                                                                           */
	/*                  This file is part of the program and library             */
	/*         SCIP --- Solving Constraint Integer Programs                      */
	/*                                                                           */
	/*    Copyright (C) 2002-2022 Konrad-Zuse-Zentrum                            */
	/*                            fuer Informationstechnik Berlin                */
	/*                                                                           */
	/*  SCIP is distributed under the terms of the ZIB Academic License.         */
	/*                                                                           */
	/*  You should have received a copy of the ZIB Academic License              */
	/*  along with SCIP; see the file COPYING. If not visit scip.zib.de.         */
	/*                                                                           */
	/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
	
	/**@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;
	   SCIP_Bool found;
	   int pos2;
	   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 */

	      if( adjacentvardata->adjacentvars == NULL )
	      {
	         found = FALSE;
	         pos2 = 0;

	      }
	      else
	      {
	         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;
	
	      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;
	
	      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);
	      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_%d", 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_%d",
	                             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;
	}