1    	/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2    	/*                                                                           */
3    	/*                  This file is part of the program and library             */
4    	/*         SCIP --- Solving Constraint Integer Programs                      */
5    	/*                                                                           */
6    	/*  Copyright (c) 2002-2023 Zuse Institute Berlin (ZIB)                      */
7    	/*                                                                           */
8    	/*  Licensed under the Apache License, Version 2.0 (the "License");          */
9    	/*  you may not use this file except in compliance with the License.         */
10   	/*  You may obtain a copy of the License at                                  */
11   	/*                                                                           */
12   	/*      http://www.apache.org/licenses/LICENSE-2.0                           */
13   	/*                                                                           */
14   	/*  Unless required by applicable law or agreed to in writing, software      */
15   	/*  distributed under the License is distributed on an "AS IS" BASIS,        */
16   	/*  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. */
17   	/*  See the License for the specific language governing permissions and      */
18   	/*  limitations under the License.                                           */
19   	/*                                                                           */
20   	/*  You should have received a copy of the Apache-2.0 license                */
21   	/*  along with SCIP; see the file LICENSE. If not visit scipopt.org.         */
22   	/*                                                                           */
23   	/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
24   	
25   	/**@file   nlhdlr_quadratic.c
26   	 * @ingroup DEFPLUGINS_NLHDLR
27   	 * @brief  nonlinear handler to handle quadratic expressions
28   	 * @author Felipe Serrano
29   	 * @author Antonia Chmiela
30   	 *
31   	 * Some definitions:
32   	 * - a `BILINEXPRTERM` is a product of two expressions
33   	 * - a `QUADEXPRTERM` stores an expression `expr` that is known to appear in a nonlinear, quadratic term, that is
34   	 *   `expr^2` or `expr*other_expr`. It stores its `sqrcoef` (that can be 0), its linear coef and all the bilinear expression
35   	 *   terms in which expr appears.
36   	 */
37   	
38   	/*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
39   	
40   	/* #define DEBUG_INTERSECTIONCUT */
41   	/* #define DEBUG_MONOIDAL */
42   	/* #define INTERCUT_MOREDEBUG */
43   	/* #define INTERCUTS_VERBOSE */
44   	
45   	#ifdef INTERCUTS_VERBOSE
46   	#define INTER_LOG
47   	#endif
48   	
49   	#ifdef INTER_LOG
50   	#define INTERLOG(x) if( SCIPgetSubscipDepth(scip) == 0 && SCIPgetVerbLevel(scip) >= SCIP_VERBLEVEL_NORMAL ) { x }
51   	#else
52   	#define INTERLOG(x)
53   	#endif
54   	
55   	#include "scip/cons_nonlinear.h"
56   	#include "scip/pub_nlhdlr.h"
57   	#include "scip/nlhdlr_quadratic.h"
58   	#include "scip/expr_pow.h"
59   	#include "scip/expr_sum.h"
60   	#include "scip/expr_var.h"
61   	#include "scip/expr_product.h"
62   	#include "scip/pub_misc_rowprep.h"
63   	
64   	/* fundamental nonlinear handler properties */
65   	#define NLHDLR_NAME                    "quadratic"
66   	#define NLHDLR_DESC                    "handler for quadratic expressions"
67   	#define NLHDLR_DETECTPRIORITY          1
68   	#define NLHDLR_ENFOPRIORITY            100
69   	
70   	/* properties of the quadratic nlhdlr statistics table */
71   	#define TABLE_NAME_QUADRATIC           "nlhdlr_quadratic"
72   	#define TABLE_DESC_QUADRATIC           "quadratic nlhdlr statistics table"
73   	#define TABLE_POSITION_QUADRATIC       14700                  /**< the position of the statistics table */
74   	#define TABLE_EARLIEST_STAGE_QUADRATIC SCIP_STAGE_TRANSFORMED /**< output of the statistics table is only printed from this stage onwards */
75   	
76   	/* some default values */
77   	#define INTERCUTS_MINVIOL              1e-4
78   	#define DEFAULT_USEINTERCUTS           FALSE
79   	#define DEFAULT_USESTRENGTH            FALSE
80   	#define DEFAULT_USEMONOIDAL            TRUE
81   	#define DEFAULT_USEMINREP              TRUE
82   	#define DEFAULT_USEBOUNDS              FALSE
83   	#define BINSEARCH_MAXITERS             120
84   	#define DEFAULT_NCUTSROOT              20
85   	#define DEFAULT_NCUTS                  2
86   	
87   	/*
88   	 * Data structures
89   	 */
90   	
91   	/** nonlinear handler expression data */
92   	struct SCIP_NlhdlrExprData
93   	{
94   	   SCIP_EXPR*            qexpr;              /**< quadratic expression (stored here again for convenient access) */
95   	
96   	   SCIP_EXPRCURV         curvature;          /**< curvature of the quadratic representation of the expression */
97   	
98   	   SCIP_INTERVAL         linactivity;        /**< activity of linear part */
99   	
100  	   /* activities of quadratic parts as defined in nlhdlrIntevalQuadratic */
101  	   SCIP_Real             minquadfiniteact;   /**< minimum activity of quadratic part where only terms with finite min
102  	                                                  activity contribute */
103  	   SCIP_Real             maxquadfiniteact;   /**< maximum activity of quadratic part where only terms with finite max
104  	                                                  activity contribute */
105  	   int                   nneginfinityquadact;/**< number of quadratic terms contributing -infinity to activity */
106  	   int                   nposinfinityquadact;/**< number of quadratic terms contributing +infinity to activity */
107  	   SCIP_INTERVAL*        quadactivities;     /**< activity of each quadratic term as defined in nlhdlrIntevalQuadratic */
108  	   SCIP_INTERVAL         quadactivity;       /**< activity of quadratic part (sum of quadactivities) */
109  	   SCIP_Longint          activitiestag;      /**< value of activities tag when activities were computed */
110  	
111  	   SCIP_CONS*            cons;               /**< if expr is the root of constraint cons, store cons; otherwise NULL */
112  	   SCIP_Bool             separating;         /**< whether we are using the nlhdlr also for separation */
113  	   SCIP_Bool             origvars;           /**< whether the quad expr in qexpr is in original (non-aux) variables */
114  	
115  	   int                   ncutsadded;         /**< number of intersection cuts added for this quadratic */
116  	};
117  	
118  	/** nonlinear handler data */
119  	struct SCIP_NlhdlrData
120  	{
121  	   int                   ncutsgenerated;     /**< total number of cuts that where generated by separateQuadratic */
122  	   int                   ncutsadded;         /**< total number of cuts that where generated by separateQuadratic and actually added */
123  	   SCIP_Longint          lastnodenumber;     /**< number of last node for which cuts were (allowed to be) generated */
124  	   int                   lastncuts;          /**< number of cuts already generated */
125  	
126  	   /* parameter */
127  	   SCIP_Bool             useintersectioncuts; /**< whether to use intersection cuts for quadratic constraints or not */
128  	   SCIP_Bool             usestrengthening;   /**< whether the strengthening should be used */
129  	   SCIP_Bool             usemonoidal;        /**< whether monoidal strengthening should be used */
130  	   SCIP_Bool             useminrep;          /**< whether the minimal representation of the S-free set should be used (instead of the gauge) */
131  	   SCIP_Bool             useboundsasrays;    /**< use bounds of variables in quadratic as rays for intersection cuts */
132  	   int                   ncutslimit;         /**< limit for number of cuts generated consecutively */
133  	   int                   ncutslimitroot;     /**< limit for number of cuts generated at root node */
134  	   int                   maxrank;            /**< maximal rank a slackvar can have */
135  	   SCIP_Real             mincutviolation;    /**< minimal cut violation the generated cuts must fulfill to be added to the LP */
136  	   SCIP_Real             minviolation;       /**< minimal violation the constraint must fulfill such that a cut can be generated */
137  	   int                   atwhichnodes;       /**< determines at which nodes cut is used (if it's -1, it's used only at the root node,
138  	                                                  if it's n >= 0, it's used at every multiple of n) */
139  	   int                   nstrengthlimit;     /**< limit for number of rays we do the strengthening for */
140  	   SCIP_Bool             sparsifycuts;       /**< should we try to sparisfy the intersection cuts? */
141  	   SCIP_Bool             ignorebadrayrestriction; /**< should cut be generated even with bad numerics when restricting to ray? */
142  	   SCIP_Bool             ignorehighre;       /**< should cut be added even when range / efficacy is large? */
143  	   SCIP_Bool             trackmore;          /**< for monoidal strengthening, should we track more statistics (more expensive) */
144  	
145  	   /* statistics */
146  	   int                   ncouldimprovedcoef; /**< number of times a coefficient could improve but didn't because of numerics */
147  	   int                   nbadrayrestriction; /**< number of times a cut was aborted because of numerics when restricting to ray */
148  	   int                   nbadnonbasic;       /**< number of times a cut was aborted because the nonbasic row was not nonbasic enough */
149  	   int                   nhighre;            /**< number of times a cut was not added because range / efficacy was too large */
150  	   int                   nphinonneg;         /**< number of times a cut was aborted because phi is nonnegative at 0 */
151  	   int                   nstrengthenings;    /**< number of successful strengthenings */
152  	   int                   nboundcuts;         /**< number of successful bound cuts */
153  	   int                   nmonoidal;          /**< number of successful monoidal strengthenings */
154  	   SCIP_Real             ncalls;             /**< number of calls to separation */
155  	   SCIP_Real             densitysum;         /**< sum of density of cuts */
156  	   SCIP_Real             cutcoefsum;         /**< sum of average cutcoefs of a cut */
157  	   SCIP_Real             monoidalimprovementsum; /**< sum of average improvement of a cut when using monoidal strengthening */
158  	   SCIP_Real             efficacysum;        /**< sum of efficacy of cuts */
159  	   SCIP_Real             currentavecutcoef;  /**< average cutcoef of current cut */
160  	   SCIP_Real             currentavemonoidalimprovement;/**< average improvement of current cut when using monoidal strengthening */
161  	};
162  	
163  	/* structure to store rays. note that for a given ray, the entries in raysidx are sorted. */
164  	struct Rays
165  	{
166  	   SCIP_Real*            rays;               /**< coefficients of rays */
167  	   int*                  raysidx;            /**< to which var the coef belongs; vars are [quadvars, linvars, auxvar] */
168  	   int*                  raysbegin;          /**< positions of rays: coefs of i-th ray [raybegin[i], raybegin[i+1]) */
169  	   int*                  lpposray;           /**< lp pos of var associated with the ray;
170  	                                               >= 0 -> lppos of var; < 0 -> var is row -lpposray -1 */
171  	
172  	   int                   rayssize;           /**< size of rays and rays idx */
173  	   int                   nrays;              /**< size of raysbegin is nrays + 1; size of lpposray */
174  	};
175  	typedef struct Rays RAYS;
176  	
177  	
178  	/*
179  	 * Callback methods of the table
180  	 */
181  	
182  	/** output method of statistics table to output file stream 'file' */
183  	static
184  	SCIP_DECL_TABLEOUTPUT(tableOutputQuadratic)
185  	{ /*lint --e{715}*/
186  	   SCIP_NLHDLR* nlhdlr;
187  	   SCIP_NLHDLRDATA* nlhdlrdata;
188  	   SCIP_CONSHDLR* conshdlr;
189  	
190  	   conshdlr = SCIPfindConshdlr(scip, "nonlinear");
191  	   assert(conshdlr != NULL);
192  	   nlhdlr = SCIPfindNlhdlrNonlinear(conshdlr, NLHDLR_NAME);
193  	   assert(nlhdlr != NULL);
194  	   nlhdlrdata = SCIPnlhdlrGetData(nlhdlr);
195  	   assert(nlhdlrdata);
196  	
197  	   /* print statistics */
198  	   SCIPinfoMessage(scip, file, "Quadratic Nlhdlr   : %10s %10s %10s %10s %10s %10s %10s %10s %10s %10s %20s %10s %10s %10s \n", "GenCuts", "AddCuts", "CouldImpr", "NLargeRE",
199  	         "AbrtBadRay", "AbrtPosPhi", "AbrtNonBas", "NStrength", "NMonoidal", "AveCutcoef", "AveMonoidalImprov", "AveDensity", "AveEfficacy", "AveBCutsFrac");
200  	   SCIPinfoMessage(scip, file, "  %-17s:", "Quadratic Nlhdlr");
201  	   SCIPinfoMessage(scip, file, " %10d", nlhdlrdata->ncutsgenerated);
202  	   SCIPinfoMessage(scip, file, " %10d", nlhdlrdata->ncutsadded);
203  	   SCIPinfoMessage(scip, file, " %10d", nlhdlrdata->ncouldimprovedcoef);
204  	   SCIPinfoMessage(scip, file, " %10d", nlhdlrdata->nhighre);
205  	   SCIPinfoMessage(scip, file, " %10d", nlhdlrdata->nbadrayrestriction);
206  	   SCIPinfoMessage(scip, file, " %10d", nlhdlrdata->nphinonneg);
207  	   SCIPinfoMessage(scip, file, " %10d", nlhdlrdata->nbadnonbasic);
208  	   SCIPinfoMessage(scip, file, " %10d", nlhdlrdata->nstrengthenings);
209  	   SCIPinfoMessage(scip, file, " %10d", nlhdlrdata->nmonoidal);
210  	   SCIPinfoMessage(scip, file, " %10g", nlhdlrdata->ncutsgenerated > 0 ? nlhdlrdata->cutcoefsum / nlhdlrdata->ncutsgenerated : 0.0);
211  	   SCIPinfoMessage(scip, file, " %20g", (nlhdlrdata->nmonoidal > 0 && nlhdlrdata->trackmore) ? nlhdlrdata->monoidalimprovementsum / nlhdlrdata->nmonoidal : -1.0);
212  	   SCIPinfoMessage(scip, file, " %10g", nlhdlrdata->ncutsgenerated > 0 ? nlhdlrdata->densitysum / nlhdlrdata->ncutsgenerated : 0.0);
213  	   SCIPinfoMessage(scip, file, " %10g", nlhdlrdata->ncutsgenerated > 0 ? nlhdlrdata->efficacysum / nlhdlrdata->ncutsgenerated : 0.0);
214  	   SCIPinfoMessage(scip, file, " %10g", nlhdlrdata->ncalls > 0 ? nlhdlrdata->nboundcuts / nlhdlrdata->ncalls : 0.0);
215  	   SCIPinfoMessage(scip, file, "\n");
216  	
217  	   return SCIP_OKAY;
218  	}
219  	
220  	
221  	/*
222  	 * static methods
223  	 */
224  	
225  	/** adds cutcoef * (col - col*) to rowprep */
226  	static
227  	SCIP_RETCODE addColToCut(
228  	   SCIP*                 scip,               /**< SCIP data structure */
229  	   SCIP_ROWPREP*         rowprep,            /**< rowprep to store intersection cut */
230  	   SCIP_SOL*             sol,                /**< solution to separate */
231  	   SCIP_Real             cutcoef,            /**< cut coefficient */
232  	   SCIP_COL*             col                 /**< column to add to rowprep */
233  	   )
234  	{
235  	   assert(col != NULL);
236  	
237  	#ifdef DEBUG_INTERCUTS_NUMERICS
238  	   SCIPinfoMessage(scip, NULL, "adding col %s to cut. %g <= col <= %g\n", SCIPvarGetName(SCIPcolGetVar(col)),
239  	      SCIPvarGetLbLocal(SCIPcolGetVar(col)), SCIPvarGetUbLocal(SCIPcolGetVar(col)));
240  	   SCIPinfoMessage(scip, NULL, "col is active at %s. Value %.15f\n", SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_LOWER ? "lower bound" :
241  	      "upper bound" , SCIPcolGetPrimsol(col));
242  	#endif
243  	
244  	   SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, SCIPcolGetVar(col), cutcoef) );
245  	   SCIProwprepAddConstant(rowprep, -cutcoef * SCIPgetSolVal(scip, sol, SCIPcolGetVar(col)) );
246  	
247  	   return SCIP_OKAY;
248  	}
249  	
250  	/** adds cutcoef * (slack - slack*) to rowprep
251  	 *
252  	  * row is lhs &le; <coefs, vars> + constant &le; rhs, thus slack is defined by
253  	  * slack + <coefs.vars> + constant = side
254  	  *
255  	  * If row (slack) is at upper, it means that <coefs,vars*> + constant = rhs, and so
256  	  * slack* = side - rhs --> slack - slack* = rhs - <coefs, vars> - constant.
257  	  *
258  	  * If row (slack) is at lower, then <coefs,vars*> + constant = lhs, and so
259  	  * slack* = side - lhs --> slack - slack* = lhs - <coefs, vars> - constant.
260  	  */
261  	static
262  	SCIP_RETCODE addRowToCut(
263  	   SCIP*                 scip,               /**< SCIP data structure */
264  	   SCIP_ROWPREP*         rowprep,            /**< rowprep to store intersection cut */
265  	   SCIP_Real             cutcoef,            /**< cut coefficient */
266  	   SCIP_ROW*             row,                /**< row, whose slack we are ading to rowprep */
267  	   SCIP_Bool*            success             /**< if the row is nonbasic enough */
268  	   )
269  	{
270  	   int i;
271  	   SCIP_COL** rowcols;
272  	   SCIP_Real* rowcoefs;
273  	   int nnonz;
274  	
275  	   assert(row != NULL);
276  	
277  	   rowcols = SCIProwGetCols(row);
278  	   rowcoefs = SCIProwGetVals(row);
279  	   nnonz = SCIProwGetNLPNonz(row);
280  	
281  	#ifdef DEBUG_INTERCUTS_NUMERICS
282  	   SCIPinfoMessage(scip, NULL, "adding slack var row_%d to cut. %g <= row <= %g\n", SCIProwGetLPPos(row), SCIProwGetLhs(row), SCIProwGetRhs(row));
283  	   SCIPinfoMessage(scip, NULL, "row is active at %s = %.15f Activity %.15f\n", SCIProwGetBasisStatus(row) == SCIP_BASESTAT_LOWER ? "lhs" :
284  	   "rhs" , SCIProwGetBasisStatus(row) == SCIP_BASESTAT_LOWER ? SCIProwGetLhs(row) : SCIProwGetRhs(row),
285  	   SCIPgetRowActivity(scip, row));
286  	#endif
287  	
288  	   if( SCIProwGetBasisStatus(row) == SCIP_BASESTAT_LOWER )
289  	   {
290  	      assert(!SCIPisInfinity(scip, -SCIProwGetLhs(row)));
291  	      if( ! SCIPisFeasEQ(scip, SCIProwGetLhs(row), SCIPgetRowActivity(scip, row)) )
292  	      {
293  	         *success = FALSE;
294  	         return SCIP_OKAY;
295  	      }
296  	
297  	      SCIProwprepAddConstant(rowprep, SCIProwGetLhs(row) * cutcoef);
298  	   }
299  	   else
300  	   {
301  	      assert(!SCIPisInfinity(scip, SCIProwGetRhs(row)));
302  	      if( ! SCIPisFeasEQ(scip, SCIProwGetRhs(row), SCIPgetRowActivity(scip, row)) )
303  	      {
304  	         *success = FALSE;
305  	         return SCIP_OKAY;
306  	      }
307  	
308  	      SCIProwprepAddConstant(rowprep, SCIProwGetRhs(row) * cutcoef);
309  	   }
310  	
311  	   for( i = 0; i < nnonz; i++ )
312  	   {
313  	      SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, SCIPcolGetVar(rowcols[i]), -rowcoefs[i] * cutcoef) );
314  	   }
315  	
316  	   SCIProwprepAddConstant(rowprep, -SCIProwGetConstant(row) * cutcoef);
317  	
318  	   return SCIP_OKAY;
319  	}
320  	
321  	/** constructs map between lp position of a basic variable and its row in the tableau */
322  	static
323  	SCIP_RETCODE constructBasicVars2TableauRowMap(
324  	   SCIP*                 scip,               /**< SCIP data structure */
325  	   int*                  map                 /**< buffer to store the map */
326  	   )
327  	{
328  	   int* basisind;
329  	   int nrows;
330  	   int i;
331  	
332  	   nrows = SCIPgetNLPRows(scip);
333  	   SCIP_CALL( SCIPallocBufferArray(scip, &basisind, nrows) );
334  	
335  	   SCIP_CALL( SCIPgetLPBasisInd(scip, basisind) );
336  	   for( i = 0; i < nrows; ++i )
337  	   {
338  	      if( basisind[i] >= 0 )
339  	         map[basisind[i]] = i;
340  	   }
341  	
342  	   SCIPfreeBufferArray(scip, &basisind);
343  	
344  	   return SCIP_OKAY;
345  	}
346  	
347  	/** counts the number of basic variables in the quadratic expr */
348  	static
349  	int countBasicVars(
350  	   SCIP_NLHDLREXPRDATA*  nlhdlrexprdata,     /**< nlhdlr expression data */
351  	   SCIP_VAR*             auxvar,             /**< aux var of expr or NULL if not needed (e.g. separating real cons) */
352  	   SCIP_Bool*            nozerostat          /**< whether there is no variable with basis status zero */
353  	   )
354  	{
355  	   SCIP_EXPR* qexpr;
356  	   SCIP_EXPR** linexprs;
357  	   SCIP_COL* col;
358  	   int i;
359  	   int nbasicvars = 0;
360  	   int nquadexprs;
361  	   int nlinexprs;
362  	
363  	   *nozerostat = TRUE;
364  	
365  	   qexpr = nlhdlrexprdata->qexpr;
366  	   SCIPexprGetQuadraticData(qexpr, NULL, &nlinexprs, &linexprs, NULL, &nquadexprs, NULL, NULL, NULL);
367  	
368  	   /* loop over quadratic vars */
369  	   for( i = 0; i < nquadexprs; ++i )
370  	   {
371  	      SCIP_EXPR* expr;
372  	
373  	      SCIPexprGetQuadraticQuadTerm(qexpr, i, &expr, NULL, NULL, NULL, NULL, NULL);
374  	
375  	      col = SCIPvarGetCol(SCIPgetExprAuxVarNonlinear(expr));
376  	      if( SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_BASIC )
377  	         nbasicvars += 1;
378  	      else if( SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_ZERO )
379  	      {
380  	         *nozerostat = FALSE;
381  	         return 0;
382  	      }
383  	   }
384  	
385  	   /* loop over linear vars */
386  	   for( i = 0; i < nlinexprs; ++i )
387  	   {
388  	      col = SCIPvarGetCol(SCIPgetExprAuxVarNonlinear(linexprs[i]));
389  	      if( SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_BASIC )
390  	         nbasicvars += 1;
391  	      else if( SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_ZERO )
392  	      {
393  	         *nozerostat = FALSE;
394  	         return 0;
395  	      }
396  	   }
397  	
398  	   /* finally consider the aux var (if it exists) */
399  	   if( auxvar != NULL )
400  	   {
401  	      col = SCIPvarGetCol(auxvar);
402  	      if( SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_BASIC )
403  	         nbasicvars += 1;
404  	      else if( SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_ZERO )
405  	      {
406  	         *nozerostat = FALSE;
407  	         return 0;
408  	      }
409  	   }
410  	
411  	   return nbasicvars;
412  	}
413  	
414  	/** stores the row of the tableau where `col` is basic
415  	 *
416  	 *  In general, we will have
417  	 *
418  	 *      basicvar1 = tableaurow var1
419  	 *      basicvar2 = tableaurow var2
420  	 *      ...
421  	 *      basicvarn = tableaurow varn
422  	 *
423  	 *  However, we want to store the the tableau row by columns. Thus, we need to know which of the basic vars `col` is.
424  	 *
425  	 *  Note we only store the entries of the nonbasic variables
426  	 */
427  	static
428  	SCIP_RETCODE storeDenseTableauRow(
429  	   SCIP*                 scip,               /**< SCIP data structure */
430  	   SCIP_COL*             col,                /**< basic columns to store its tableau row */
431  	   int*                  basicvarpos2tableaurow,/**< map between basic var and its tableau row */
432  	   int                   nbasiccol,          /**< which basic var this is */
433  	   int                   raylength,          /**< the length of a ray (the total number of basic vars) */
434  	   SCIP_Real*            binvrow,            /**< buffer to store row of Binv */
435  	   SCIP_Real*            binvarow,           /**< buffer to store row of Binv A */
436  	   SCIP_Real*            tableaurows         /**< pointer to store the tableau rows */
437  	   )
438  	{
439  	   int nrows;
440  	   int ncols;
441  	   int lppos;
442  	   int i;
443  	   SCIP_COL** cols;
444  	   SCIP_ROW** rows;
445  	   int nray;
446  	
447  	   assert(nbasiccol < raylength);
448  	   assert(col != NULL);
449  	   assert(binvrow != NULL);
450  	   assert(binvarow != NULL);
451  	   assert(tableaurows != NULL);
452  	   assert(basicvarpos2tableaurow != NULL);
453  	   assert(SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_BASIC);
454  	
455  	   SCIP_CALL( SCIPgetLPRowsData(scip, &rows, &nrows) );
456  	   SCIP_CALL( SCIPgetLPColsData(scip, &cols, &ncols) );
457  	
458  	   lppos = SCIPcolGetLPPos(col);
459  	
460  	   assert(basicvarpos2tableaurow[lppos] >= 0);
461  	
462  	   SCIP_CALL( SCIPgetLPBInvRow(scip, basicvarpos2tableaurow[lppos], binvrow, NULL, NULL) );
463  	   SCIP_CALL( SCIPgetLPBInvARow(scip, basicvarpos2tableaurow[lppos], binvrow, binvarow, NULL, NULL) );
464  	
465  	   nray = 0;
466  	   for( i = 0; i < ncols; ++i )
467  	      if( SCIPcolGetBasisStatus(cols[i]) != SCIP_BASESTAT_BASIC )
468  	      {
469  	         tableaurows[nbasiccol + nray * raylength] = binvarow[i];
470  	         nray++;
471  	      }
472  	   for( ; i < ncols + nrows; ++i )
473  	      if( SCIProwGetBasisStatus(rows[i - ncols]) != SCIP_BASESTAT_BASIC )
474  	      {
475  	         tableaurows[nbasiccol + nray * raylength] = binvrow[i - ncols];
476  	         nray++;
477  	      }
478  	
479  	   return SCIP_OKAY;
480  	}
481  	
482  	/** stores the rows of the tableau corresponding to the basic variables in the quadratic expression
483  	 *
484  	 * Also return a map storing to which var the entry of a ray corresponds, i.e., if the tableau is
485  	 *
486  	 *     basicvar_1 = ray1_1 nonbasicvar_1 + ...
487  	 *     basicvar_2 = ray1_2 nonbasicvar_1 + ...
488  	 *     ...
489  	 *     basicvar_n = ray1_n nonbasicvar_1 + ...
490  	 *
491  	 * The map maps k to the position of basicvar_k in the variables of the constraint assuming the variables are sorted as
492  	 * [quadratic vars, linear vars, auxvar].
493  	 */
494  	static
495  	SCIP_RETCODE storeDenseTableauRowsByColumns(
496  	   SCIP*                 scip,               /**< SCIP data structure */
497  	   SCIP_NLHDLREXPRDATA*  nlhdlrexprdata,     /**< nlhdlr expression data */
498  	   int                   raylength,          /**< length of a ray of the tableau */
499  	   SCIP_VAR*             auxvar,             /**< aux var of expr or NULL if not needed (e.g. separating real cons) */
500  	   SCIP_Real*            tableaurows,        /**< buffer to store the tableau rows */
501  	   int*                  rayentry2conspos    /**< buffer to store the map */
502  	   )
503  	{
504  	   SCIP_EXPR* qexpr;
505  	   SCIP_EXPR** linexprs;
506  	   SCIP_Real* binvarow;
507  	   SCIP_Real* binvrow;
508  	   SCIP_COL* col;
509  	   int* basicvarpos2tableaurow; /* map between basic var and its tableau row */
510  	   int nrayentries;
511  	   int nquadexprs;
512  	   int nlinexprs;
513  	   int nrows;
514  	   int ncols;
515  	   int i;
516  	
517  	   nrows = SCIPgetNLPRows(scip);
518  	   ncols = SCIPgetNLPCols(scip);
519  	
520  	   SCIP_CALL( SCIPallocBufferArray(scip, &basicvarpos2tableaurow, ncols) );
521  	   SCIP_CALL( SCIPallocBufferArray(scip, &binvrow, nrows) );
522  	   SCIP_CALL( SCIPallocBufferArray(scip, &binvarow, ncols) );
523  	
524  	   for( i = 0; i < ncols; ++i )
525  	      basicvarpos2tableaurow[i] = -1;
526  	   SCIP_CALL( constructBasicVars2TableauRowMap(scip, basicvarpos2tableaurow) );
527  	
528  	   qexpr = nlhdlrexprdata->qexpr;
529  	   SCIPexprGetQuadraticData(qexpr, NULL, &nlinexprs, &linexprs, NULL, &nquadexprs, NULL, NULL, NULL);
530  	
531  	   /* entries of quadratic basic vars */
532  	   nrayentries = 0;
533  	   for( i = 0; i < nquadexprs; ++i )
534  	   {
535  	      SCIP_EXPR* expr;
536  	      SCIPexprGetQuadraticQuadTerm(qexpr, i, &expr, NULL, NULL, NULL, NULL, NULL);
537  	
538  	      col = SCIPvarGetCol(SCIPgetExprAuxVarNonlinear(expr));
539  	      if( SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_BASIC )
540  	      {
541  	         SCIP_CALL( storeDenseTableauRow(scip, col, basicvarpos2tableaurow, nrayentries, raylength, binvrow, binvarow,
542  	                  tableaurows) );
543  	
544  	         rayentry2conspos[nrayentries] = i;
545  	         nrayentries++;
546  	      }
547  	   }
548  	   /* entries of linear vars */
549  	   for( i = 0; i < nlinexprs; ++i )
550  	   {
551  	      col = SCIPvarGetCol(SCIPgetExprAuxVarNonlinear(linexprs[i]));
552  	      if( SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_BASIC )
553  	      {
554  	         SCIP_CALL( storeDenseTableauRow(scip, col, basicvarpos2tableaurow, nrayentries, raylength, binvrow, binvarow,
555  	                  tableaurows) );
556  	
557  	         rayentry2conspos[nrayentries] = nquadexprs + i;
558  	         nrayentries++;
559  	      }
560  	   }
561  	   /* entry of aux var (if it exists) */
562  	   if( auxvar != NULL )
563  	   {
564  	      col = SCIPvarGetCol(auxvar);
565  	      if( SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_BASIC )
566  	      {
567  	         SCIP_CALL( storeDenseTableauRow(scip, col, basicvarpos2tableaurow, nrayentries, raylength, binvrow, binvarow,
568  	                  tableaurows) );
569  	
570  	         rayentry2conspos[nrayentries] = nquadexprs + nlinexprs;
571  	         nrayentries++;
572  	      }
573  	   }
574  	   assert(nrayentries == raylength);
575  	
576  	#ifdef  DEBUG_INTERSECTIONCUT
577  	   for( i = 0; i < ncols; ++i )
578  	   {
579  	      SCIPinfoMessage(scip, NULL, "%d column of tableau is: ",i);
580  	      for( int j = 0; j < raylength; ++j )
581  	         SCIPinfoMessage(scip, NULL, "%g ", tableaurows[i * raylength + j]);
582  	      SCIPinfoMessage(scip, NULL, "\n");
583  	   }
584  	#endif
585  	
586  	   SCIPfreeBufferArray(scip, &binvarow);
587  	   SCIPfreeBufferArray(scip, &binvrow);
588  	   SCIPfreeBufferArray(scip, &basicvarpos2tableaurow);
589  	
590  	   return SCIP_OKAY;
591  	}
592  	
593  	/** initializes rays data structure */
594  	static
595  	SCIP_RETCODE createRays(
596  	   SCIP*                 scip,               /**< SCIP data structure */
597  	   RAYS**                rays                /**< rays data structure */
598  	   )
599  	{
600  	   SCIP_CALL( SCIPallocBuffer(scip, rays) );
601  	   BMSclearMemory(*rays);
602  	
603  	   SCIP_CALL( SCIPallocBufferArray(scip, &(*rays)->rays, SCIPgetNLPCols(scip)) );
604  	   SCIP_CALL( SCIPallocBufferArray(scip, &(*rays)->raysidx, SCIPgetNLPCols(scip)) );
605  	
606  	   /* overestimate raysbegin and lpposray */
607  	   SCIP_CALL( SCIPallocBufferArray(scip, &(*rays)->raysbegin, SCIPgetNLPCols(scip) + SCIPgetNLPRows(scip) + 1) );
608  	   SCIP_CALL( SCIPallocBufferArray(scip, &(*rays)->lpposray, SCIPgetNLPCols(scip) + SCIPgetNLPRows(scip)) );
609  	   (*rays)->raysbegin[0] = 0;
610  	
611  	   (*rays)->rayssize = SCIPgetNLPCols(scip);
612  	
613  	   return SCIP_OKAY;
614  	}
615  	
616  	/** initializes rays data structure for bound rays */
617  	static
618  	SCIP_RETCODE createBoundRays(
619  	   SCIP*                 scip,               /**< SCIP data structure */
620  	   RAYS**                rays,               /**< rays data structure */
621  	   int                   size                /**< number of rays to allocate */
622  	   )
623  	{
624  	   SCIP_CALL( SCIPallocBuffer(scip, rays) );
625  	   BMSclearMemory(*rays);
626  	
627  	   SCIP_CALL( SCIPallocBufferArray(scip, &(*rays)->rays, size) );
628  	   SCIP_CALL( SCIPallocBufferArray(scip, &(*rays)->raysidx, size) );
629  	
630  	   /* overestimate raysbegin and lpposray */
631  	   SCIP_CALL( SCIPallocBufferArray(scip, &(*rays)->raysbegin, size + 1) );
632  	   SCIP_CALL( SCIPallocBufferArray(scip, &(*rays)->lpposray, size) );
633  	   (*rays)->raysbegin[0] = 0;
634  	
635  	   (*rays)->rayssize = size;
636  	
637  	   return SCIP_OKAY;
638  	}
639  	
640  	/** frees rays data structure */
641  	static
642  	void freeRays(
643  	   SCIP*                 scip,               /**< SCIP data structure */
644  	   RAYS**                rays                /**< rays data structure */
645  	   )
646  	{
647  	   if( *rays == NULL )
648  	      return;
649  	
650  	   SCIPfreeBufferArray(scip, &(*rays)->lpposray);
651  	   SCIPfreeBufferArray(scip, &(*rays)->raysbegin);
652  	   SCIPfreeBufferArray(scip, &(*rays)->raysidx);
653  	   SCIPfreeBufferArray(scip, &(*rays)->rays);
654  	
655  	   SCIPfreeBuffer(scip, rays);
656  	}
657  	
658  	/** inserts entry to rays data structure; checks and resizes if more space is needed */
659  	static
660  	SCIP_RETCODE insertRayEntry(
661  	   SCIP*                 scip,               /**< SCIP data structure */
662  	   RAYS*                 rays,               /**< rays data structure */
663  	   SCIP_Real             coef,               /**< coefficient to insert */
664  	   int                   coefidx,            /**< index of coefficient (conspos of var it corresponds to) */
665  	   int                   coefpos             /**< where to insert coefficient */
666  	   )
667  	{
668  	   /* check for size */
669  	   if( rays->rayssize <= coefpos + 1 )
670  	   {
671  	      int newsize;
672  	
673  	      newsize = SCIPcalcMemGrowSize(scip, coefpos + 1);
674  	      SCIP_CALL( SCIPreallocBufferArray(scip, &(rays->rays), newsize) );
675  	      SCIP_CALL( SCIPreallocBufferArray(scip, &(rays->raysidx), newsize) );
676  	      rays->rayssize = newsize;
677  	   }
678  	
679  	   /* insert entry */
680  	   rays->rays[coefpos] = coef;
681  	   rays->raysidx[coefpos] = coefidx;
682  	
683  	   return SCIP_OKAY;
684  	}
685  	
686  	/** constructs map between the lppos of a variables and its position in the constraint assuming the constraint variables
687  	 * are sorted as [quad vars, lin vars, aux var (if it exists)]
688  	 *
689  	 * If a variable doesn't appear in the constraint, then its position is -1.
690  	 */
691  	static
692  	void constructLPPos2ConsPosMap(
693  	   SCIP_NLHDLREXPRDATA*  nlhdlrexprdata,     /**< nlhdlr expression data */
694  	   SCIP_VAR*             auxvar,             /**< aux var of the expr */
695  	   int*                  map                 /**< buffer to store the mapping */
696  	   )
697  	{
698  	   SCIP_EXPR* qexpr;
699  	   SCIP_EXPR** linexprs;
700  	   int nquadexprs;
701  	   int nlinexprs;
702  	   int lppos;
703  	   int i;
704  	
705  	   qexpr = nlhdlrexprdata->qexpr;
706  	   SCIPexprGetQuadraticData(qexpr, NULL, &nlinexprs, &linexprs, NULL, &nquadexprs, NULL, NULL, NULL);
707  	
708  	   /* set pos of quadratic vars */
709  	   for( i = 0; i < nquadexprs; ++i )
710  	   {
711  	      SCIP_EXPR* expr;
712  	      SCIPexprGetQuadraticQuadTerm(qexpr, i, &expr, NULL, NULL, NULL, NULL, NULL);
713  	
714  	      lppos = SCIPcolGetLPPos(SCIPvarGetCol(SCIPgetExprAuxVarNonlinear(expr)));
715  	      map[lppos] = i;
716  	   }
717  	   /* set pos of lin vars */
718  	   for( i = 0; i < nlinexprs; ++i )
719  	   {
720  	      lppos = SCIPcolGetLPPos(SCIPvarGetCol(SCIPgetExprAuxVarNonlinear(linexprs[i])));
721  	      map[lppos] = nquadexprs + i;
722  	   }
723  	   /* set pos of aux var (if it exists) */
724  	   if( auxvar != NULL )
725  	   {
726  	      lppos = SCIPcolGetLPPos(SCIPvarGetCol(auxvar));
727  	      map[lppos] = nquadexprs + nlinexprs;
728  	   }
729  	
730  	   return;
731  	}
732  	
733  	/** inserts entries of factor * nray-th column of densetableaucols into rays data structure */
734  	static
735  	SCIP_RETCODE insertRayEntries(
736  	   SCIP*                 scip,               /**< SCIP data structure */
737  	   RAYS*                 rays,               /**< rays data structure */
738  	   SCIP_Real*            densetableaucols,   /**< column of the tableau in dense format */
739  	   int*                  rayentry2conspos,   /**< map between rayentry and conspos of associated var */
740  	   int                   raylength,          /**< length of a tableau column */
741  	   int                   nray,               /**< which tableau column to insert */
742  	   int                   conspos,            /**< conspos of ray's nonbasic var in the cons; -1 if not in the cons */
743  	   SCIP_Real             factor,             /**< factor to multiply each tableau col */
744  	   int*                  nnonz,              /**< position to start adding the ray in rays and buffer to store nnonz */
745  	   SCIP_Bool*            success             /**< we can't separate if there is a nonzero ray with basis status ZERO */
746  	   )
747  	{
748  	   int i;
749  	
750  	   *success = TRUE;
751  	
752  	   for( i = 0; i < raylength; ++i )
753  	   {
754  	      SCIP_Real coef;
755  	
756  	      /* we have a nonzero ray with base stat zero -> can't generate cut */
757  	      if( factor == 0.0 && ! SCIPisZero(scip, densetableaucols[nray * raylength + i]) )
758  	      {
759  	         *success = FALSE;
760  	         return SCIP_OKAY;
761  	      }
762  	
763  	      coef = factor * densetableaucols[nray * raylength + i];
764  	
765  	      /* this might be a source of numerical issues
766  	       * TODO: in case of problems, an idea would be to scale the ray entries; compute the cut coef and scale it back down
767  	       * another idea would be to check against a smaller epsilion.
768  	       * The problem is that if the cut coefficient is 1/t where lpsol + t*ray intersects the S-free set.
769  	       * Now if t is super big, then a super small coefficient would have had an impact...
770  	       */
771  	      if( SCIPisZero(scip, coef) )
772  	         continue;
773  	
774  	      /* check if nonbasic var entry should come before this one */
775  	      if( conspos > -1 && conspos < rayentry2conspos[i] )
776  	      {
777  	         /* add nonbasic entry */
778  	         assert(factor != 0.0);
779  	
780  	#ifdef  DEBUG_INTERSECTIONCUT
781  	         SCIPinfoMessage(scip, NULL, "ray belongs to nonbasic var pos %d in cons\n", conspos);
782  	#endif
783  	
784  	         SCIP_CALL( insertRayEntry(scip, rays, -factor, conspos, *nnonz) );
785  	         (*nnonz)++;
786  	
787  	         /* we are done with nonbasic entry */
788  	         conspos = -1;
789  	      }
790  	
791  	      SCIP_CALL( insertRayEntry(scip, rays, coef, rayentry2conspos[i], *nnonz) );
792  	      (*nnonz)++;
793  	   }
794  	
795  	   /* if nonbasic entry was not added and should still be added, then it should go at the end */
796  	   if( conspos > -1 )
797  	   {
798  	      /* add nonbasic entry */
799  	      assert(factor != 0.0);
800  	
801  	#ifdef  DEBUG_INTERSECTIONCUT
802  	      SCIPinfoMessage(scip, NULL, "ray belongs to nonbasic var pos %d in cons\n", conspos);
803  	#endif
804  	
805  	      SCIP_CALL( insertRayEntry(scip, rays, -factor, conspos, *nnonz) );
806  	      (*nnonz)++;
807  	   }
808  	
809  	   /* finished ray entry; store its end */
810  	   rays->raysbegin[rays->nrays + 1] = *nnonz;
811  	
812  	#ifdef  DEBUG_INTERSECTIONCUT
813  	   SCIPinfoMessage(scip, NULL, "entries of ray %d are between [%d, %d):\n", rays->nrays, rays->raysbegin[rays->nrays], *nnonz);
814  	   for( i = rays->raysbegin[rays->nrays]; i < *nnonz; ++i )
815  	      SCIPinfoMessage(scip, NULL, "(%d, %g), ", rays->raysidx[i], rays->rays[i]);
816  	   SCIPinfoMessage(scip, NULL, "\n");
817  	#endif
818  	
819  	   return SCIP_OKAY;
820  	}
821  	
822  	/** stores rays in sparse form
823  	 *
824  	 * The first rays correspond to the nonbasic variables
825  	 * and the last rays to the nonbasic slack variables.
826  	 *
827  	 * More details: The LP tableau is of the form
828  	 *
829  	 *     basicvar_1 = ray1_1 nonbasicvar_1 + ... + raym_1 nonbasicvar_m
830  	 *     basicvar_2 = ray1_2 nonbasicvar_1 + ... + raym_2 nonbasicvar_m
831  	 *     ...
832  	 *     basicvar_n = ray1_n nonbasicvar_1 + ... + raym_n nonbasicvar_m
833  	 *     nonbasicvar_1 = 1.0 nonbasicvar_1 + ... +    0.0 nonbasicvar_m
834  	 *     ...
835  	 *     nonbasicvar_m = 0.0 nonbasicvar_1 + ... +    1.0 nonbasicvar_m
836  	 *
837  	 *  so rayk = (rayk_1, ... rayk_n, e_k)
838  	 *  We store the entries of the rays associated to the variables present in the quadratic expr.
839  	 *  We do not store zero rays.
840  	 *
841  	 *  Also, we store the rays as if every nonbasic variable was at lower (so that all rays moves to infinity)
842  	 *  Since the tableau is:
843  	 *
844  	 *      basicvar + Binv L (nonbasic_lower - lb) + Binv U (nonbasic_upper - ub) = basicvar_sol
845  	 *
846  	 *  then:
847  	 *
848  	 *      basicvar = basicvar_sol - Binv L (nonbasic_lower - lb) + Binv U (ub - nonbasic_upper)
849  	 *
850  	 *  and so the entries of the rays associated with the basic variables are:
851  	 *  rays_basicvars = [-BinvL, BinvU].
852  	 *
853  	 *  So we flip the sign of the rays associated to nonbasic vars at lower.
854  	 *  In constrast, the nonbasic part of the ray has a 1.0 for nonbasic at lower and a -1.0 for nonbasic at upper, i.e.
855  	 *  nonbasic_lower = lb + 1.0(nonbasic_lower - lb) and
856  	 *  nonbasic_upper = ub - 1.0(ub - nonbasic_upper)
857  	 */
858  	static
859  	SCIP_RETCODE createAndStoreSparseRays(
860  	   SCIP*                 scip,               /**< SCIP data structure */
861  	   SCIP_NLHDLREXPRDATA*  nlhdlrexprdata,     /**< nlhdlr expression data */
862  	   SCIP_VAR*             auxvar,             /**< aux var of expr or NULL if not needed (e.g. separating real cons) */
863  	   RAYS**                raysptr,            /**< buffer to store rays datastructure */
864  	   SCIP_Bool*            success             /**< we can't separate if there is a var with basis status ZERO */
865  	   )
866  	{
867  	   SCIP_Real* densetableaucols;
868  	   SCIP_COL** cols;
869  	   SCIP_ROW** rows;
870  	   RAYS* rays;
871  	   int* rayentry2conspos;
872  	   int* lppos2conspos;
873  	   int nnonbasic;
874  	   int nrows;
875  	   int ncols;
876  	   int nnonz;
877  	   int raylength;
878  	   int i;
879  	
880  	   nrows = SCIPgetNLPRows(scip);
881  	   ncols = SCIPgetNLPCols(scip);
882  	
883  	   *success = TRUE;
884  	
885  	   raylength = countBasicVars(nlhdlrexprdata, auxvar, success);
886  	   if( ! *success )
887  	   {
888  	      SCIPdebugMsg(scip, "failed to store sparse rays: there is a var with base status zero\n");
889  	      return SCIP_OKAY;
890  	   }
891  	
892  	   SCIP_CALL( SCIPallocBufferArray(scip, &densetableaucols, raylength * (ncols + nrows)) );
893  	   SCIP_CALL( SCIPallocBufferArray(scip, &rayentry2conspos, raylength) );
894  	
895  	   /* construct dense tableau and map between ray entries and position of corresponding var in quad cons */
896  	   SCIP_CALL( storeDenseTableauRowsByColumns(scip, nlhdlrexprdata, raylength, auxvar,
897  	            densetableaucols, rayentry2conspos) );
898  	
899  	   /* build rays sparsely now */
900  	   SCIP_CALL( SCIPallocBufferArray(scip, &lppos2conspos, ncols) );
901  	   for( i = 0; i < ncols; ++i )
902  	      lppos2conspos[i] = -1;
903  	
904  	   constructLPPos2ConsPosMap(nlhdlrexprdata, auxvar, lppos2conspos);
905  	
906  	   /* store sparse rays */
907  	   SCIP_CALL( createRays(scip, raysptr) );
908  	   rays = *raysptr;
909  	
910  	   nnonz = 0;
911  	   nnonbasic = 0;
912  	
913  	   /* go through nonbasic variables */
914  	   cols = SCIPgetLPCols(scip);
915  	   for( i = 0; i < ncols; ++i )
916  	   {
917  	      int oldnnonz = nnonz;
918  	      SCIP_COL* col;
919  	      SCIP_Real factor;
920  	
921  	      col = cols[i];
922  	
923  	      /* set factor to store basic entries of ray as = [-BinvL, BinvU] */
924  	      if( SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_LOWER )
925  	         factor = -1.0;
926  	      else if( SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_UPPER )
927  	         factor = 1.0;
928  	      else if( SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_ZERO )
929  	         factor = 0.0;
930  	      else
931  	         continue;
932  	
933  	      SCIP_CALL( insertRayEntries(scip, rays, densetableaucols, rayentry2conspos, raylength, nnonbasic,
934  	               lppos2conspos[SCIPcolGetLPPos(col)], factor, &nnonz, success) );
935  	
936  	#ifdef  DEBUG_INTERSECTIONCUT
937  	      SCIPinfoMessage(scip, NULL, "looked at ray of var %s with basestat %d, it has %d nonzeros\n-----------------\n",
938  	            SCIPvarGetName(SCIPcolGetVar(col)), SCIPcolGetBasisStatus(col), nnonz - oldnnonz);
939  	#endif
940  	      if( ! (*success) )
941  	      {
942  	#ifdef  DEBUG_INTERSECTIONCUT
943  	         SCIPdebugMsg(scip, "nonzero ray associated with variable <%s> has base status zero -> abort storing rays\n",
944  	               SCIPvarGetName(SCIPcolGetVar(col)));
945  	#endif
946  	         goto CLEANUP;
947  	      }
948  	
949  	      /* if ray is non zero remember who it belongs to */
950  	      assert(oldnnonz <= nnonz);
951  	      if( oldnnonz < nnonz )
952  	      {
953  	         rays->lpposray[rays->nrays] = SCIPcolGetLPPos(col);
954  	         (rays->nrays)++;
955  	      }
956  	      nnonbasic++;
957  	   }
958  	
959  	   /* go through nonbasic slack variables */
960  	   rows = SCIPgetLPRows(scip);
961  	   for( i = 0; i < nrows; ++i )
962  	   {
963  	      int oldnnonz = nnonz;
964  	      SCIP_ROW* row;
965  	      SCIP_Real factor;
966  	
967  	      row = rows[i];
968  	
969  	      /* set factor to store basic entries of ray as = [-BinvL, BinvU]; basic status of rows are flipped! See lpi.h! */
970  	      assert(SCIProwGetBasisStatus(row) != SCIP_BASESTAT_ZERO);
971  	      if( SCIProwGetBasisStatus(row) == SCIP_BASESTAT_LOWER )
972  	         factor = 1.0;
973  	      else if( SCIProwGetBasisStatus(row) == SCIP_BASESTAT_UPPER )
974  	         factor =-1.0;
975  	      else
976  	         continue;
977  	
978  	      SCIP_CALL( insertRayEntries(scip, rays, densetableaucols, rayentry2conspos, raylength, nnonbasic, -1, factor,
979  	               &nnonz, success) );
980  	      assert(*success);
981  	
982  	      /* if ray is non zero remember who it belongs to */
983  	#ifdef  DEBUG_INTERSECTIONCUT
984  	      SCIPinfoMessage(scip, NULL, "looked at ray of row %d, it has %d nonzeros\n-----------------\n", i, nnonz - oldnnonz);
985  	#endif
986  	      assert(oldnnonz <= nnonz);
987  	      if( oldnnonz < nnonz )
988  	      {
989  	         rays->lpposray[rays->nrays] = -SCIProwGetLPPos(row) - 1;
990  	         (rays->nrays)++;
991  	      }
992  	      nnonbasic++;
993  	   }
994  	
995  	CLEANUP:
996  	   SCIPfreeBufferArray(scip, &lppos2conspos);
997  	   SCIPfreeBufferArray(scip, &rayentry2conspos);
998  	   SCIPfreeBufferArray(scip, &densetableaucols);
999  	
1000 	   if( ! *success )
1001 	   {
1002 	      freeRays(scip, &rays);
1003 	   }
1004 	
1005 	   return SCIP_OKAY;
1006 	}
1007 	
1008 	/* TODO: which function this comment belongs to? */
1009 	/* this function determines how the maximal S-free set is going to look like
1010 	 *
1011 	 * There are 4 possibilities: after writing the quadratic constraint
1012 	 * \f$q(z) \leq 0\f$
1013 	 * as
1014 	 * \f$\Vert x(z)\Vert^2 - \Vert y\Vert^2 + w(z) + kappa \leq 0\f$,
1015 	 * the cases are determined according to the following:
1016 	 * - Case 1: w = 0 and kappa = 0
1017 	 * - Case 2: w = 0 and kappa > 0
1018 	 * - Case 3: w = 0 and kappa < 0
1019 	 * - Case 4: w != 0
1020 	 */
1021 	
1022 	/** compute quantities for intersection cuts
1023 	 *
1024 	 * Assume the quadratic is stored as
1025 	 * \f[ q(z) = z_q^T Q z_q + b_q^T z_q + b_l z_l + c - z_a \f]
1026 	 * where:
1027 	 *  - \f$z_q\f$ are the quadratic vars
1028 	 *  - \f$z_l\f$ are the linear vars
1029 	 *  - \f$z_a\f$ is the aux var if it exists
1030 	 *
1031 	 * We can rewrite it as
1032 	 * \f[ \Vert x(z)\Vert^2 - \Vert y\Vert^2 + w(z) + \kappa \leq 0. \f]
1033 	 * To do this transformation and later to compute the actual cut we need to compute and store some quantities.
1034 	 * Let
1035 	 *    - \f$I_0\f$, \f$I_+\f$, and \f$I_-\f$ be the index set of zero, positive, and negative eigenvalues, respectively
1036 	 *    - \f$v_i\f$ be the i-th eigenvector of \f$Q\f$
1037 	 *    - \f$zlp\f$ be the lp value of the variables \f$z\f$
1038 	 *
1039 	 * The quantities we need are:
1040 	 *    - \f$vb_i = v_i^T b\f$ for \f$i \in I_+ \cup I_-\f$
1041 	 *    - \f$vzlp_i = v_i^T zlp_q\f$ for \f$i \in I_+ \cup I_-\f$
1042 	 *    - \f$\kappa = c - 1/4 \sum_{i \in I_+ \cup I_-} (v_i^T b_q)^2 / eigval_i\f$
1043 	 *    - \f$w(z) = (\sum_{i \in I_0} v_i^T b_q v_i^T) z_q + b_l^T z_l - z_a\f$
1044 	 *    - \f$w(zlp)\f$
1045 	 *
1046 	 * @return \f$\kappa\f$ and the vector \f$\sum_{i \in I_0} v_i^T b_q v_i^T\f$
1047 	 *
1048 	 * @note if the constraint is q(z) &le; rhs, then the constant when writing the constraint as quad &le; 0 is c - rhs.
1049 	 * @note if the quadratic constraint we are separating is q(z) &ge; lhs, then we multiply by -1.
1050 	 * In practice, what changes is
1051 	 *    - the sign of the eigenvalues
1052 	 *    - the sign of \f$b_q\f$ and \f$b_l\f$
1053 	 *    - the sign of the coefficient of the auxvar (if it exists)
1054 	 *    - the constant of the quadratic written as quad &le; 0 is lhs - c
1055 	 * @note The eigenvectors _do not_ change sign!
1056 	 */
1057 	static
1058 	SCIP_RETCODE intercutsComputeCommonQuantities(
1059 	   SCIP*                 scip,               /**< SCIP data structure */
1060 	   SCIP_NLHDLREXPRDATA*  nlhdlrexprdata,     /**< nlhdlr expression data */
1061 	   SCIP_VAR*             auxvar,             /**< aux var of expr or NULL if not needed (e.g. separating real cons) */
1062 	   SCIP_Real             sidefactor,         /**< 1.0 if the violated constraint is q &le; rhs, -1.0 otherwise */
1063 	   SCIP_SOL*             sol,                /**< solution to separate */
1064 	   SCIP_Real*            vb,                 /**< buffer to store \f$v_i^T b\f$ for \f$i \in I_+ \cup I_-\f$ */
1065 	   SCIP_Real*            vzlp,               /**< buffer to store \f$v_i^T zlp_q\f$ for \f$i \in I_+ \cup I_-\f$ */
1066 	   SCIP_Real*            wcoefs,             /**< buffer to store the coefs of quad vars of w */
1067 	   SCIP_Real*            wzlp,               /**< pointer to store the value of w at zlp */
1068 	   SCIP_Real*            kappa               /**< pointer to store the value of kappa */
1069 	   )
1070 	{
1071 	   SCIP_EXPR* qexpr;
1072 	   SCIP_EXPR** linexprs;
1073 	   SCIP_Real* eigenvectors;
1074 	   SCIP_Real* eigenvalues;
1075 	   SCIP_Real* lincoefs;
1076 	   SCIP_Real constant; /* constant of the quadratic when written as <= 0 */
1077 	   int nquadexprs;
1078 	   int nlinexprs;
1079 	   int i;
1080 	   int j;
1081 	
1082 	   assert(sidefactor == 1.0 || sidefactor == -1.0);
1083 	   assert(nlhdlrexprdata->cons != NULL || auxvar != NULL );
1084 	
1085 	   qexpr = nlhdlrexprdata->qexpr;
1086 	   SCIPexprGetQuadraticData(qexpr, &constant, &nlinexprs, &linexprs, &lincoefs, &nquadexprs, NULL, &eigenvalues,
1087 	         &eigenvectors);
1088 	
1089 	   assert( eigenvalues != NULL );
1090 	
1091 	   /* first get constant of quadratic when written as quad <= 0 */
1092 	   if( nlhdlrexprdata->cons != NULL )
1093 	      constant = (sidefactor == 1.0) ? constant - SCIPgetRhsNonlinear(nlhdlrexprdata->cons) :
1094 	         SCIPgetLhsNonlinear(nlhdlrexprdata->cons) - constant;
1095 	   else
1096 	      constant = (sidefactor * constant);
1097 	
1098 	   *kappa = 0.0;
1099 	   *wzlp = 0.0;
1100 	   BMSclearMemoryArray(wcoefs, nquadexprs);
1101 	
1102 	   for( i = 0; i < nquadexprs; ++i )
1103 	   {
1104 	      SCIP_Real vdotb;
1105 	      SCIP_Real vdotzlp;
1106 	      int offset;
1107 	
1108 	      offset = i * nquadexprs;
1109 	
1110 	      /* compute v_i^T b and v_i^T zlp */
1111 	      vdotb = 0;
1112 	      vdotzlp = 0;
1113 	      for( j = 0; j < nquadexprs; ++j )
1114 	      {
1115 	         SCIP_EXPR* expr;
1116 	         SCIP_Real lincoef;
1117 	
1118 	         SCIPexprGetQuadraticQuadTerm(qexpr, j, &expr, &lincoef, NULL, NULL, NULL, NULL);
1119 	
1120 	         vdotb += (sidefactor * lincoef) * eigenvectors[offset + j];
1121 	
1122 	#ifdef INTERCUT_MOREDEBUG
1123 	         printf("vdotb: offset %d, eigenvector %d = %g, lincoef quad %g\n", offset, j,
1124 	               eigenvectors[offset + j], lincoef);
1125 	#endif
1126 	         vdotzlp += SCIPgetSolVal(scip, sol, SCIPgetExprAuxVarNonlinear(expr)) * eigenvectors[offset + j];
1127 	      }
1128 	      vb[i] = vdotb;
1129 	      vzlp[i] = vdotzlp;
1130 	
1131 	      if( ! SCIPisZero(scip, eigenvalues[i]) )
1132 	      {
1133 	         /* nonzero eigenvalue: compute kappa */
1134 	         *kappa += SQR(vdotb) / (sidefactor * eigenvalues[i]);
1135 	      }
1136 	      else
1137 	      {
1138 	         /* compute coefficients of w and compute w at zlp */
1139 	         for( j = 0; j < nquadexprs; ++j )
1140 	            wcoefs[j] += vdotb * eigenvectors[offset + j];
1141 	
1142 	         *wzlp += vdotb * vdotzlp;
1143 	      }
1144 	   }
1145 	
1146 	   /* finish kappa computation */
1147 	   *kappa *= -0.25;
1148 	   *kappa += constant;
1149 	
1150 	   if( SCIPisZero(scip, *kappa) )
1151 	      *kappa = 0.0;
1152 	
1153 	   /* finish w(zlp) computation: linear part (including auxvar, if applicable) */
1154 	   for( i = 0; i < nlinexprs; ++i )
1155 	   {
1156 	      *wzlp += (sidefactor * lincoefs[i]) * SCIPgetSolVal(scip, sol, SCIPgetExprAuxVarNonlinear(linexprs[i]));
1157 	   }
1158 	   if( auxvar != NULL )
1159 	   {
1160 	      *wzlp += (sidefactor * -1.0) * SCIPgetSolVal(scip, sol, auxvar);
1161 	   }
1162 	
1163 	#ifdef  DEBUG_INTERSECTIONCUT
1164 	   SCIPinfoMessage(scip, NULL, "Computed common quantities needed for intercuts:\n");
1165 	   SCIPinfoMessage(scip, NULL, "   kappa = %g\n   quad part w = ", *kappa);
1166 	   for( i = 0; i < nquadexprs; ++i )
1167 	   {
1168 	      SCIPinfoMessage(scip, NULL, "%g ", wcoefs[i]);
1169 	   }
1170 	   SCIPinfoMessage(scip, NULL, "\n");
1171 	#endif
1172 	   return SCIP_OKAY;
1173 	}
1174 	
1175 	/** computes eigenvec^T ray */
1176 	static
1177 	SCIP_Real computeEigenvecDotRay(
1178 	   SCIP_Real*            eigenvec,           /**< eigenvector */
1179 	   int                   nquadvars,          /**< number of quadratic vars (length of eigenvec) */
1180 	   SCIP_Real*            raycoefs,           /**< coefficients of ray */
1181 	   int*                  rayidx,             /**< index of consvar the ray coef is associated to */
1182 	   int                   raynnonz            /**< length of raycoefs and rayidx */
1183 	   )
1184 	{
1185 	   SCIP_Real retval;
1186 	   int i;
1187 	
1188 	   retval = 0.0;
1189 	   for( i = 0; i < raynnonz; ++i )
1190 	   {
1191 	      /* rays are sorted; the first entries correspond to the quad vars, so we stop after first nonquad var entry */
1192 	      if( rayidx[i] >= nquadvars )
1193 	         break;
1194 	
1195 	      retval += eigenvec[rayidx[i]] * raycoefs[i];
1196 	   }
1197 	
1198 	   return retval;
1199 	}
1200 	
1201 	/** computes linear part of evaluation of w(ray): \f$b_l^T ray_l - ray_a\f$
1202 	 *
1203 	 * @note we can know whether the auxiliary variable appears by the entries of the ray
1204 	 */
1205 	static
1206 	SCIP_Real computeWRayLinear(
1207 	   SCIP_NLHDLREXPRDATA*  nlhdlrexprdata,     /**< nlhdlr expression data */
1208 	   SCIP_Real             sidefactor,         /**< 1.0 if the violated constraint is q &le; rhs, -1.0 otherwise */
1209 	   SCIP_Real*            raycoefs,           /**< coefficients of ray */
1210 	   int*                  rayidx,             /**< ray coef[i] affects var at pos rayidx[i] in consvar */
1211 	   int                   raynnonz            /**< length of raycoefs and rayidx */
1212 	   )
1213 	{
1214 	   SCIP_EXPR* qexpr;
1215 	   SCIP_Real* lincoefs;
1216 	   SCIP_Real retval;
1217 	   int nquadexprs;
1218 	   int nlinexprs;
1219 	
1220 	   int i;
1221 	   int start;
1222 	
1223 	#ifdef INTERCUT_MOREDEBUG
1224 	   printf("Computing w(ray) \n");
1225 	#endif
1226 	   retval = 0.0;
1227 	   start = raynnonz - 1;
1228 	
1229 	   qexpr = nlhdlrexprdata->qexpr;
1230 	   SCIPexprGetQuadraticData(qexpr, NULL, &nlinexprs, NULL, &lincoefs, &nquadexprs, NULL, NULL, NULL);
1231 	
1232 	   /* process ray entry associated to the auxvar if it applies */
1233 	   if( rayidx[raynnonz - 1] == nquadexprs + nlinexprs )
1234 	   {
1235 	#ifdef INTERCUT_MOREDEBUG
1236 	      printf("wray auxvar term %g \n", (sidefactor * -1.0) * raycoefs[raynnonz - 1]);
1237 	#endif
1238 	      retval += (sidefactor * -1.0) * raycoefs[raynnonz - 1];
1239 	      start--;
1240 	   }
1241 	
1242 	   /* process the rest of the entries */
1243 	   for( i = start; i >= 0; --i )
1244 	   {
1245 	      /* rays are sorted; last entries correspond to the lin vars, so we stop after first quad var entry */
1246 	      if( rayidx[i] < nquadexprs )
1247 	         break;
1248 	
1249 	#ifdef INTERCUT_MOREDEBUG
1250 	      printf("wray var in pos %d term %g:: lincoef %g raycoef %g\n", rayidx[i], (sidefactor *
1251 	            lincoefs[rayidx[i] - nquadexprs]) * raycoefs[i], lincoefs[rayidx[i] - nquadexprs] ,raycoefs[i]);
1252 	#endif
1253 	      retval += (sidefactor * lincoefs[rayidx[i] - nquadexprs]) * raycoefs[i] ;
1254 	   }
1255 	
1256 	   return retval;
1257 	}
1258 	
1259 	/** computes the dot product of v_i and the current ray as well as of v_i and the apex where v_i
1260 	 * is the i-th eigenvalue
1261 	 */
1262 	static
1263 	void computeVApexAndVRay(
1264 	   SCIP_NLHDLREXPRDATA*  nlhdlrexprdata,     /**< nlhdlr expression data */
1265 	   SCIP_Real*            apex,               /**< array containing the apex of the S-free set in the original space */
1266 	   SCIP_Real*            raycoefs,           /**< coefficients of ray */
1267 	   int*                  rayidx,             /**< index of consvar the ray coef is associated to */
1268 	   int                   raynnonz,           /**< length of raycoefs and rayidx */
1269 	   SCIP_Real*            vapex,              /**< array to store \f$v_i^T apex\f$ */
1270 	   SCIP_Real*            vray                /**< array to store \f$v_i^T ray\f$ */
1271 	   )
1272 	{
1273 	   SCIP_EXPR* qexpr;
1274 	   int nquadexprs;
1275 	   SCIP_Real* eigenvectors;
1276 	   SCIP_Real* eigenvalues;
1277 	   int i;
1278 	
1279 	   qexpr = nlhdlrexprdata->qexpr;
1280 	   SCIPexprGetQuadraticData(qexpr, NULL, NULL, NULL, NULL, &nquadexprs, NULL, &eigenvalues, &eigenvectors);
1281 	
1282 	   for( i = 0; i < nquadexprs; ++i )
1283 	   {
1284 	      SCIP_Real vdotapex;
1285 	      SCIP_Real vdotray;
1286 	      int offset;
1287 	      int j;
1288 	      int k;
1289 	
1290 	      offset = i * nquadexprs;
1291 	
1292 	      /* compute v_i^T apex and v_i^T ray */
1293 	      vdotapex = 0.0;
1294 	      vdotray = 0.0;
1295 	      k = 0;
1296 	      for( j = 0; j < nquadexprs; ++j )
1297 	      {
1298 	         SCIP_Real rayentry;
1299 	
1300 	         /* get entry of ray -> check if current var index corresponds to a non-zero entry in ray */
1301 	         if( k < raynnonz && j == rayidx[k] )
1302 	         {
1303 	            rayentry = raycoefs[k];
1304 	            ++k;
1305 	         }
1306 	         else
1307 	            rayentry = 0.0;
1308 	
1309 	         vdotray += rayentry * eigenvectors[offset + j];
1310 	         vdotapex += apex[j] * eigenvectors[offset + j];
1311 	      }
1312 	
1313 	      vray[i] = vdotray;
1314 	      vapex[i] = vdotapex;
1315 	   }
1316 	}
1317 	
1318 	/** calculate coefficients of restriction of the function to given ray.
1319 	 *
1320 	 * The restriction of the function representing the maximal S-free set to zlp + t * ray has the form
1321 	 * SQRT(A t^2 + B t + C) - (D t + E) for cases 1, 2, and 3.
1322 	 * For case 4 it is a piecewise defined function and each piece is of the aforementioned form.
1323 	 *
1324 	 * This function computes the coefficients A, B, C, D, E for the given ray.
1325 	 * In case 4, it computes the coefficients for both pieces, in addition to coefficients needed to evaluate the condition
1326 	 * in the piecewise definition of the function.
1327 	 *
1328 	 * The parameter iscase4 tells the function if it is case 4 or not.
1329 	 */
1330 	static
1331 	SCIP_RETCODE computeRestrictionToLine(
1332 	   SCIP*                 scip,               /**< SCIP data structure */
1333 	   SCIP_NLHDLREXPRDATA*  nlhdlrexprdata,     /**< nlhdlr expression data */
1334 	   SCIP_Real             sidefactor,         /**< 1.0 if the violated constraint is q &le; rhs, -1.0 otherwise */
1335 	   SCIP_Real*            raycoefs,           /**< coefficients of ray */
1336 	   int*                  rayidx,             /**< index of consvar the ray coef is associated to */
1337 	   int                   raynnonz,           /**< length of raycoefs and rayidx */
1338 	   SCIP_Real*            vb,                 /**< array containing \f$v_i^T b\f$ for \f$i \in I_+ \cup I_-\f$ */
1339 	   SCIP_Real*            vzlp,               /**< array containing \f$v_i^T zlp_q\f$ for \f$i \in I_+ \cup I_-\f$ */
1340 	   SCIP_Real             kappa,              /**< value of kappa */
1341 	   SCIP_Real*            apex,               /**< apex of C */
1342 	   SCIP_Real*            coefs2,             /**< buffer to store A, B, C, D, and E of case 2 */
1343 	   SCIP_Bool*            success             /**< did we successfully compute the coefficients? */
1344 	   )
1345 	{
1346 	   SCIP_EXPR* qexpr;
1347 	   int nquadexprs;
1348 	   SCIP_Real* eigenvectors;
1349 	   SCIP_Real* eigenvalues;
1350 	   SCIP_Real* a;
1351 	   SCIP_Real* b;
1352 	   SCIP_Real* c;
1353 	   SCIP_Real* d;
1354 	   SCIP_Real* e;
1355 	   SCIP_Real* vray;
1356 	   SCIP_Real* vapex;
1357 	   SCIP_Real norm;
1358 	   int i;
1359 	
1360 	   *success = TRUE;
1361 	
1362 	   qexpr = nlhdlrexprdata->qexpr;
1363 	   SCIPexprGetQuadraticData(qexpr, NULL, NULL, NULL, NULL, &nquadexprs, NULL, &eigenvalues, &eigenvectors);
1364 	
1365 	#ifdef  DEBUG_INTERSECTIONCUT
1366 	   SCIPinfoMessage(scip, NULL, "\n############################################\n");
1367 	   SCIPinfoMessage(scip, NULL, "Restricting to line:\n");
1368 	#endif
1369 	
1370 	   assert(coefs2 != NULL);
1371 	
1372 	   /* set all coefficients to zero */
1373 	   BMSclearMemoryArray(coefs2, 5);
1374 	
1375 	   /* compute v_i^T apex in vapex[i] and v_i^T ray in vray[i] */
1376 	   SCIP_CALL( SCIPallocBufferArray(scip, &vapex, nquadexprs) );
1377 	   SCIP_CALL( SCIPallocBufferArray(scip, &vray, nquadexprs) );
1378 	   computeVApexAndVRay(nlhdlrexprdata, apex, raycoefs, rayidx, raynnonz, vapex, vray);
1379 	
1380 	   a = coefs2;
1381 	   b = coefs2 + 1;
1382 	   c = coefs2 + 2;
1383 	   d = coefs2 + 3;
1384 	   e = coefs2 + 4;
1385 	
1386 	   norm = 0.0;
1387 	   for( i = 0; i < nquadexprs; ++i )
1388 	   {
1389 	      SCIP_Real dot;
1390 	      SCIP_Real eigval;
1391 	
1392 	      eigval = sidefactor * eigenvalues[i];
1393 	
1394 	      if( SCIPisZero(scip, eigval) )
1395 	         continue;
1396 	
1397 	      dot = vzlp[i] + vb[i] / (2.0 * eigval);
1398 	
1399 	      if( eigval > 0.0 )
1400 	      {
1401 	         *d += eigval * dot * (vzlp[i] - vapex[i]);
1402 	         *e += eigval * dot * (vray[i] + vapex[i] + vb[i] / (2.0 * eigval));
1403 	         norm += eigval * SQR(dot);
1404 	      }
1405 	      else
1406 	      {
1407 	         *a -= eigval * SQR(vzlp[i] - vapex[i]);
1408 	         *b -= eigval * (vzlp[i] - vapex[i]) * (vray[i] + vapex[i] + vb[i] / (2.0 * eigval));
1409 	         *c -= eigval * SQR(vray[i] + vapex[i] + vb[i] / (2.0 * eigval));
1410 	      }
1411 	   }
1412 	
1413 	   norm = SQRT(norm / kappa + 1.0);
1414 	   *a /= kappa;
1415 	   *b /= kappa;
1416 	   *c /= kappa;
1417 	   *d /= (kappa * norm);
1418 	   *e /= (kappa * norm);
1419 	   *e += 1.0 / norm;
1420 	
1421 	   /* In theory, the function at 0 must be negative. Because of bad numerics this might not always hold, so we abort
1422 	    * the generation of the cut in this case.
1423 	    */
1424 	   if( SQRT( *c ) - *e >= 0 )
1425 	   {
1426 	      /* check if it's really a numerical problem */
1427 	      assert(SQRT( *c ) > 10e+15 || *e > 10e+15 || SQRT( *c ) - *e < 10e+9);
1428 	
1429 	      INTERLOG(printf("Bad numerics: phi(0) >= 0\n"); )
1430 	      *success = FALSE;
1431 	      goto TERMINATE;
1432 	   }
1433 	
1434 	#ifdef  DEBUG_INTERSECTIONCUT
1435 	   SCIPinfoMessage(scip, NULL, "Restriction yields case 2: a,b,c,d,e %g %g %g %g %g\n", coefs1234a[0], coefs1234a[1], coefs1234a[2],
1436 	         coefs1234a[3], coefs1234a[4]);
1437 	#endif
1438 	
1439 	   /* some sanity check applicable to all cases */
1440 	   assert(*a >= 0); /* the function inside the root is convex */
1441 	   assert(*c >= 0); /* radicand at zero */
1442 	
1443 	TERMINATE:
1444 	   /* free memory */
1445 	   SCIPfreeBufferArray(scip, &vray);
1446 	   SCIPfreeBufferArray(scip, &vapex);
1447 	
1448 	   return SCIP_OKAY;
1449 	}
1450 	
1451 	/** calculate coefficients of restriction of the function to given ray.
1452 	 *
1453 	 * The restriction of the function representing the maximal S-free set to zlp + t * ray has the form
1454 	 * SQRT(A t^2 + B t + C) - (D t + E) for cases 1, 2, and 3.
1455 	 * For case 4 it is a piecewise defined function and each piece is of the aforementioned form.
1456 	 *
1457 	 * This function computes the coefficients A, B, C, D, E for the given ray.
1458 	 * In case 4, it computes the coefficients for both pieces, in addition to coefficients needed to evaluate the condition
1459 	 * in the piecewise definition of the function.
1460 	 *
1461 	 * The parameter iscase4 tells the function if it is case 4 or not.
1462 	 */
1463 	static
1464 	SCIP_RETCODE computeRestrictionToRay(
1465 	   SCIP*                 scip,               /**< SCIP data structure */
1466 	   SCIP_NLHDLREXPRDATA*  nlhdlrexprdata,     /**< nlhdlr expression data */
1467 	   SCIP_Real             sidefactor,         /**< 1.0 if the violated constraint is q &le; rhs, -1.0 otherwise */
1468 	   SCIP_Bool             iscase4,            /**< whether we are in case 4 */
1469 	   SCIP_Real*            raycoefs,           /**< coefficients of ray */
1470 	   int*                  rayidx,             /**< index of consvar the ray coef is associated to */
1471 	   int                   raynnonz,           /**< length of raycoefs and rayidx */
1472 	   SCIP_Real*            vb,                 /**< array containing \f$v_i^T b\f$ for \f$i \in I_+ \cup I_-\f$ */
1473 	   SCIP_Real*            vzlp,               /**< array containing \f$v_i^T zlp_q\f$ for \f$i \in I_+ \cup I_-\f$ */
1474 	   SCIP_Real*            wcoefs,             /**< coefficients of w for the qud vars or NULL if w is 0 */
1475 	   SCIP_Real             wzlp,               /**< value of w at zlp */
1476 	   SCIP_Real             kappa,              /**< value of kappa */
1477 	   SCIP_Real*            coefs1234a,         /**< buffer to store A, B, C, D, and E of cases 1, 2, 3, or 4a */
1478 	   SCIP_Real*            coefs4b,            /**< buffer to store A, B, C, D, and E of case 4b (or NULL if not needed) */
1479 	   SCIP_Real*            coefscondition,     /**< buffer to store data to evaluate condition to decide case 4a or 4b */
1480 	   SCIP_Bool*            success             /**< did we successfully compute the coefficients? */
1481 	   )
1482 	{
1483 	   SCIP_EXPR* qexpr;
1484 	   int nquadexprs;
1485 	   SCIP_Real* eigenvectors;
1486 	   SCIP_Real* eigenvalues;
1487 	   SCIP_Real* a;
1488 	   SCIP_Real* b;
1489 	   SCIP_Real* c;
1490 	   SCIP_Real* d;
1491 	   SCIP_Real* e;
1492 	   SCIP_Real wray;
1493 	   int i;
1494 	
1495 	   *success = TRUE;
1496 	
1497 	   qexpr = nlhdlrexprdata->qexpr;
1498 	   SCIPexprGetQuadraticData(qexpr, NULL, NULL, NULL, NULL, &nquadexprs, NULL, &eigenvalues, &eigenvectors);
1499 	
1500 	#ifdef  DEBUG_INTERSECTIONCUT
1501 	   SCIPinfoMessage(scip, NULL, "\n############################################\n");
1502 	   SCIPinfoMessage(scip, NULL, "Restricting to ray:\n");
1503 	   for( i = 0; i < raynnonz; ++i )
1504 	   {
1505 	      SCIPinfoMessage(scip, NULL, "(%d, %g), ", rayidx[i], raycoefs[i]);
1506 	   }
1507 	   SCIPinfoMessage(scip, NULL, "\n");
1508 	#endif
1509 	
1510 	   assert(coefs1234a != NULL);
1511 	
1512 	   /* set all coefficients to zero */
1513 	   BMSclearMemoryArray(coefs1234a, 5);
1514 	   if( iscase4 )
1515 	   {
1516 	      assert(coefs4b != NULL);
1517 	      assert(coefscondition != NULL);
1518 	      assert(wcoefs != NULL);
1519 	
1520 	      BMSclearMemoryArray(coefs4b, 5);
1521 	      BMSclearMemoryArray(coefscondition, 3);
1522 	   }
1523 	
1524 	   a = coefs1234a;
1525 	   b = coefs1234a + 1;
1526 	   c = coefs1234a + 2;
1527 	   d = coefs1234a + 3;
1528 	   e = coefs1234a + 4;
1529 	   wray = 0;
1530 	
1531 	   for( i = 0; i < nquadexprs; ++i )
1532 	   {
1533 	      SCIP_Real dot = 0.0;
1534 	      SCIP_Real vdotray;
1535 	
1536 	      if( SCIPisZero(scip, eigenvalues[i]) )
1537 	      {
1538 	         if( wcoefs == NULL )
1539 	            continue;
1540 	      }
1541 	      else
1542 	      {
1543 	         dot = vzlp[i] + vb[i] / (2.0 * (sidefactor * eigenvalues[i]));
1544 	      }
1545 	      vdotray = computeEigenvecDotRay(&eigenvectors[i * nquadexprs], nquadexprs, raycoefs, rayidx, raynnonz);
1546 	
1547 	      if( SCIPisZero(scip, eigenvalues[i]) )
1548 	      {
1549 	         /* zero eigenvalue (and wcoefs not null) -> case 4: compute w(r) */
1550 	         assert(wcoefs != NULL);
1551 	         wray += vb[i] * vdotray;
1552 	#ifdef INTERCUT_MOREDEBUG
1553 	         printf(" wray += %g, vb[i] %g and vdotray %g\n", vb[i] * vdotray,vb[i],vdotray);
1554 	#endif
1555 	      }
1556 	      else if( sidefactor * eigenvalues[i] > 0 )
1557 	      {
1558 	         /* positive eigenvalue: compute common part of D and E */
1559 	         *d += (sidefactor * eigenvalues[i]) * dot * vdotray;
1560 	         *e += (sidefactor * eigenvalues[i]) * SQR( dot );
1561 	
1562 	#ifdef INTERCUT_MOREDEBUG
1563 	         printf("Positive eigenvalue: computing D: v^T ray %g, v^T( zlp + b/theta ) %g and theta %g \n", vdotray, dot, (sidefactor * eigenvalues[i]));
1564 	#endif
1565 	      }
1566 	      else
1567 	      {
1568 	         /* negative eigenvalue: compute common part of A, B, and C */
1569 	         *a -= (sidefactor * eigenvalues[i]) * SQR( vdotray );
1570 	         *b -= 2.0 * (sidefactor * eigenvalues[i]) *  dot * vdotray;
1571 	         *c -= (sidefactor * eigenvalues[i]) * SQR( dot );
1572 	
1573 	#ifdef INTERCUT_MOREDEBUG
1574 	         printf("Negative eigenvalue: computing A: v^T ray %g, and theta %g \n", vdotray, (sidefactor * eigenvalues[i]));
1575 	#endif
1576 	      }
1577 	   }
1578 	
1579 	   if( ! iscase4 )
1580 	   {
1581 	      /* We are in one of the first 3 cases */
1582 	      *e += MAX(kappa, 0.0);
1583 	      *c -= MIN(kappa, 0.0);
1584 	
1585 	      /* finish computation of D and E */
1586 	      assert(*e > 0);
1587 	      *e = SQRT( *e );
1588 	      *d /= *e;
1589 	
1590 	      /* some sanity checks only applicable to these cases (more at the end) */
1591 	      assert(*c >= 0);
1592 	
1593 	      /* In theory, the function at 0 must be negative. Because of bad numerics this might not always hold, so we abort
1594 	       * the generation of the cut in this case.
1595 	       */
1596 	      if( SQRT( *c ) - *e >= 0 )
1597 	      {
1598 	         /* check if it's really a numerical problem */
1599 	         assert(SQRT( *c ) > 10e+15 || *e > 10e+15 || SQRT( *c ) - *e < 10e+9);
1600 	
1601 	         INTERLOG(printf("Bad numerics: phi(0) >= 0\n"); )
1602 	         *success = FALSE;
1603 	         return SCIP_OKAY;
1604 	      }
1605 	   }
1606 	   else
1607 	   {
1608 	      SCIP_Real norm;
1609 	      SCIP_Real xextra;
1610 	      SCIP_Real yextra;
1611 	
1612 	      norm = SQRT( 1 + SQR( kappa ) );
1613 	      xextra = wzlp + kappa + norm;
1614 	      yextra = wzlp + kappa - norm;
1615 	
1616 	      /* finish computing w(ray), the linear part is missing */
1617 	      wray += computeWRayLinear(nlhdlrexprdata, sidefactor, raycoefs, rayidx, raynnonz);
1618 	
1619 	      /*
1620 	       * coefficients of case 4b
1621 	       */
1622 	      /* at this point E is \|x(zlp)\|^2, so we can finish A, B, and C */
1623 	      coefs4b[0] = (*a) * (*e);
1624 	      coefs4b[1] = (*b) * (*e);
1625 	      coefs4b[2] = (*c) * (*e);
1626 	
1627 	      /* finish D and E */
1628 	      coefs4b[3] = *d;
1629 	      coefs4b[4] = (*e) + xextra / 2.0;
1630 	
1631 	      /* when \|x(zlp)\|^2 is too large, we can divide everything by \|x(zlp)\| */
1632 	      if( *e > 100 )
1633 	      {
1634 	         coefs4b[0] = (*a);
1635 	         coefs4b[1] = (*b);
1636 	         coefs4b[2] = (*c);
1637 	
1638 	         /* finish D and E */
1639 	         coefs4b[3] = *d / SQRT( *e );
1640 	         coefs4b[4] = SQRT( *e ) + (xextra / (2.0 * SQRT( *e )));
1641 	      }
1642 	
1643 	      /*
1644 	       * coefficients of case 4a
1645 	       */
1646 	      /* finish A, B, and C */
1647 	      *a += SQR( wray ) / (4.0 * norm);
1648 	      *b += 2.0 * yextra * (wray) / (4.0 * norm);
1649 	      *c += SQR( yextra ) / (4.0 * norm);
1650 	
1651 	      /* finish D and E */
1652 	      *e +=  SQR( xextra ) / (4.0 * norm);
1653 	      *e = SQRT( *e );
1654 	
1655 	      *d += xextra * (wray) / (4.0 * norm);
1656 	      *d /= *e;
1657 	
1658 	      /*
1659 	       * coefficients of condition: stores -numerator of x_{r+1}/ norm xhat, w(ray), and numerator of y_{s+1} at zlp
1660 	       *
1661 	       */
1662 	      /* at this point E is \| \hat{x} (zlp)\| */
1663 	      coefscondition[0] = - xextra / (*e);
1664 	      coefscondition[1] = wray;
1665 	      coefscondition[2] = yextra;
1666 	   }
1667 	
1668 	#ifdef  DEBUG_INTERSECTIONCUT
1669 	   if( ! iscase4 )
1670 	   {
1671 	      SCIPinfoMessage(scip, NULL, "Restriction yields case 1,2 or 3: a,b,c,d,e %g %g %g %g %g\n", coefs1234a[0], coefs1234a[1], coefs1234a[2],
1672 	            coefs1234a[3], coefs1234a[4]);
1673 	   }
1674 	   else
1675 	   {
1676 	      SCIPinfoMessage(scip, NULL, "Restriction yields\n   Case 4a: a,b,c,d,e %g %g %g %g %g\n", coefs1234a[0],
1677 	            coefs1234a[1], coefs1234a[2], coefs1234a[3], coefs1234a[4]);
1678 	      SCIPinfoMessage(scip, NULL, "   Case 4b: a,b,c,d,e %g %g %g %g %g\n", coefs4b[0], coefs4b[1], coefs4b[2],
1679 	            coefs4b[3], coefs4b[4]);
1680 	      SCIPinfoMessage(scip, NULL, "   Condition: xextra/e, wray, yextra %g %g %g g\n", coefscondition[0],
1681 	            coefscondition[1], coefscondition[2]);
1682 	   }
1683 	#endif
1684 	
1685 	   /* some sanity check applicable to all cases */
1686 	   assert(*a >= 0); /* the function inside the root is convex */
1687 	   assert(*c >= 0); /* radicand at zero */
1688 	
1689 	   if( iscase4 )
1690 	   {
1691 	      assert(coefs4b[0] >= 0); /* the function inside the root is convex */
1692 	      assert(coefs4b[2] >= 0); /* radicand at zero */
1693 	   }
1694 	   /*assert(4.0 * (*a) * (*c) >= SQR( *b ) ); *//* the function is defined everywhere, so minimum of radicand must be nonnegative */
1695 	
1696 	   return SCIP_OKAY;
1697 	}
1698 	
1699 	/** returns phi(zlp + t * ray) = SQRT(A t^2 + B t + C) - (D t + E) */
1700 	static
1701 	SCIP_Real evalPhiAtRay(
1702 	   SCIP_Real             t,                  /**< argument of phi restricted to ray */
1703 	   SCIP_Real             a,                  /**< value of A */
1704 	   SCIP_Real             b,                  /**< value of B */
1705 	   SCIP_Real             c,                  /**< value of C */
1706 	   SCIP_Real             d,                  /**< value of D */
1707 	   SCIP_Real             e                   /**< value of E */
1708 	   )
1709 	{
1710 	#ifdef INTERCUTS_DBLDBL
1711 	   SCIP_Real QUAD(lin);
1712 	   SCIP_Real QUAD(disc);
1713 	   SCIP_Real QUAD(tmp);
1714 	   SCIP_Real QUAD(root);
1715 	
1716 	   /* d * t + e */
1717 	   SCIPquadprecProdDD(lin, d, t);
1718 	   SCIPquadprecSumQD(lin, lin, e);
1719 	
1720 	   /* a * t * t */
1721 	   SCIPquadprecSquareD(disc, t);
1722 	   SCIPquadprecProdQD(disc, disc, a);
1723 	
1724 	   /* b * t */
1725 	   SCIPquadprecProdDD(tmp, b, t);
1726 	
1727 	   /* a * t * t + b * t */
1728 	   SCIPquadprecSumQQ(disc, disc, tmp);
1729 	
1730 	   /* a * t * t + b * t + c */
1731 	   SCIPquadprecSumQD(disc, disc, c);
1732 	
1733 	   /* sqrt(above): can't take sqrt of 0! */
1734 	   if( QUAD_TO_DBL(disc) == 0 )
1735 	   {
1736 	      QUAD_ASSIGN(root, 0.0);
1737 	   }
1738 	   else
1739 	   {
1740 	      SCIPquadprecSqrtQ(root, disc);
1741 	   }
1742 	
1743 	   /* final result */
1744 	   QUAD_SCALE(lin, -1.0);
1745 	   SCIPquadprecSumQQ(tmp, root, lin);
1746 	
1747 	   assert(t != 1e20 || QUAD_TO_DBL(tmp) <= 0);
1748 	
1749 	   return  QUAD_TO_DBL(tmp);
1750 	#else
1751 	   return SQRT( a * t * t + b * t + c ) - ( d * t + e );
1752 	#endif
1753 	}
1754 	
1755 	/** checks whether case 4a applies
1756 	 *
1757 	 * The condition for being in case 4a is
1758 	 * \f[ -\lambda_{r+1} \Vert \hat y(zlp + tsol\, ray)\Vert + \hat y_{s+1}(zlp + tsol\, ray) \leq 0\f]
1759 	 *
1760 	 * This reduces to
1761 	 * \f[ -num(\hat x_{r+1}(zlp)) \sqrt{A t^2 + B t + C} / E  + w(ray) \cdot t + num(\hat y_{s+1}(zlp)) \leq 0\f]
1762 	 * where num is the numerator.
1763 	 */
1764 	static
1765 	SCIP_Real isCase4a(
1766 	   SCIP_Real             tsol,               /**< t in the above formula */
1767 	   SCIP_Real*            coefs4a,            /**< coefficients A, B, C, D, and E of case 4a */
1768 	   SCIP_Real*            coefscondition      /**< extra coefficients needed for the evaluation of the condition:
1769 	                                              *   \f$num(\hat x_{r+1}(zlp)) / E\f$; \f$w(ray)\f$; \f$num(\hat y_{s+1}(zlp))\f$ */
1770 	   )
1771 	{
1772 	   return (coefscondition[0] * SQRT( coefs4a[0] * SQR( tsol ) + coefs4a[1] * tsol + coefs4a[2] ) + coefscondition[1] *
1773 	         tsol + coefscondition[2]) <= 0.0;
1774 	}
1775 	
1776 	/** helper function of computeRoot: we want phi to be &le; 0 */
1777 	static
1778 	void doBinarySearch(
1779 	   SCIP*                 scip,               /**< SCIP data structure */
1780 	   SCIP_Real             a,                  /**< value of A */
1781 	   SCIP_Real             b,                  /**< value of B */
1782 	   SCIP_Real             c,                  /**< value of C */
1783 	   SCIP_Real             d,                  /**< value of D */
1784 	   SCIP_Real             e,                  /**< value of E */
1785 	   SCIP_Real*            sol                 /**< buffer to store solution; also gives initial point */
1786 	   )
1787 	{
1788 	   SCIP_Real lb = 0.0;
1789 	   SCIP_Real ub = *sol;
1790 	   SCIP_Real curr;
1791 	   int i;
1792 	
1793 	   for( i = 0; i < BINSEARCH_MAXITERS; ++i )
1794 	   {
1795 	      SCIP_Real phival;
1796 	
1797 	      curr = (lb + ub) / 2.0;
1798 	      phival = evalPhiAtRay(curr, a, b, c, d, e);
1799 	#ifdef INTERCUT_MOREDEBUG
1800 	      printf("%d: lb,ub %.10f, %.10f. curr = %g -> phi at curr %g -> phi at lb %g \n", i, lb, ub, curr, phival, evalPhiAtRay(lb, a, b, c, d, e));
1801 	#endif
1802 	
1803 	      if( phival <= 0.0 )
1804 	      {
1805 	         lb = curr;
1806 	         if( SCIPisFeasZero(scip, phival) || SCIPisFeasEQ(scip, ub, lb) )
1807 	            break;
1808 	      }
1809 	      else
1810 	         ub = curr;
1811 	   }
1812 	
1813 	   *sol = lb;
1814 	}
1815 	
1816 	/** finds smallest positive root phi by finding the smallest positive root of
1817 	 * (A - D^2) t^2 + (B - 2 D*E) t + (C - E^2) = 0
1818 	 *
1819 	 * However, we are conservative and want a solution such that phi is negative, but close to 0.
1820 	 * Thus, we correct the result with a binary search.
1821 	 */
1822 	static
1823 	SCIP_Real computeRoot(
1824 	   SCIP*                 scip,               /**< SCIP data structure */
1825 	   SCIP_Real*            coefs               /**< value of A */
1826 	   )
1827 	{
1828 	   SCIP_Real sol;
1829 	   SCIP_INTERVAL bounds;
1830 	   SCIP_INTERVAL result;
1831 	   SCIP_Real a = coefs[0];
1832 	   SCIP_Real b = coefs[1];
1833 	   SCIP_Real c = coefs[2];
1834 	   SCIP_Real d = coefs[3];
1835 	   SCIP_Real e = coefs[4];
1836 	
1837 	   /* there is an intersection point if and only if SQRT(A) > D: here we are beliving in math, this might cause
1838 	    * numerical issues
1839 	    */
1840 	   if( SQRT( a ) <= d )
1841 	   {
1842 	      sol = SCIPinfinity(scip);
1843 	      return sol;
1844 	   }
1845 	
1846 	   SCIPintervalSetBounds(&bounds, 0.0, SCIPinfinity(scip));
1847 	
1848 	   /* SCIPintervalSolveUnivariateQuadExpressionPositiveAllScalar finds all x such that a x^2 + b x >= c and x in bounds.
1849 	    * it is known that if tsol is the root we are looking for, then gamma(zlp + t * ray) <= 0 between 0 and tsol, thus
1850 	    * tsol is the smallest t such that (A - D^2) t^2 + (B - 2 D*E) t + (C - E^2) >= 0
1851 	    */
1852 	   SCIPintervalSolveUnivariateQuadExpressionPositiveAllScalar(SCIP_INTERVAL_INFINITY, &result, a - d * d, b - 2.0 * d *
1853 	         e, -(c - e * e), bounds);
1854 	
1855 	   /* it can still be empty because of our infinity, I guess... */
1856 	   sol = SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, result) ? SCIPinfinity(scip) : SCIPintervalGetInf(result);
1857 	
1858 	#ifdef INTERCUT_MOREDEBUG
1859 	   {
1860 	      SCIP_Real binsol;
1861 	      binsol = SCIPinfinity(scip);
1862 	      doBinarySearch(scip, a, b, c, d, e, &binsol);
1863 	      printf("got root %g: with binsearch get %g\n", sol, binsol);
1864 	   }
1865 	#endif
1866 	
1867 	   /* check that solution is acceptable, ideally it should be <= 0, however when it is positive, we trigger a binary
1868 	    * search to make it negative. This binary search might return a solution point that is not at accurately 0 as the
1869 	    * one obtained from the function above. Thus, it might fail to satisfy the condition of case 4b in some cases, e.g.,
1870 	    * ex8_3_1, bchoco05, etc
1871 	    */
1872 	   if( evalPhiAtRay(sol, a, b, c, d, e) <= 1e-10 )
1873 	   {
1874 	#ifdef INTERCUT_MOREDEBUG
1875 	      printf("interval solution returned %g -> phival = %g, believe it\n", sol, evalPhiAtRay(sol, a, b, c, d, e));
1876 	      printf("don't do bin search\n");
1877 	#endif
1878 	      return sol;
1879 	   }
1880 	   else
1881 	   {
1882 	      /* perform a binary search to make it negative: this might correct a wrong infinity (e.g. crudeoil_lee1_05) */
1883 	#ifdef INTERCUT_MOREDEBUG
1884 	      printf("do bin search because phival is %g\n", evalPhiAtRay(sol, a, b, c, d, e));
1885 	#endif
1886 	      doBinarySearch(scip, a, b, c, d, e, &sol);
1887 	   }
1888 	
1889 	   return sol;
1890 	}
1891 	
1892 	/** The maximal S-free set is \f$\gamma(z) \leq 0\f$; we find the intersection point of the ray `ray` starting from zlp with the
1893 	 * boundary of the S-free set.
1894 	 * That is, we find \f$t \geq 0\f$ such that \f$\gamma(zlp + t \cdot \text{ray}) = 0\f$.
1895 	 *
1896 	 * In cases 1,2, and 3, gamma is of the form
1897 	 *    \f[ \gamma(zlp + t \cdot \text{ray}) = \sqrt{A t^2 + B t + C} - (D t + E) \f]
1898 	 *
1899 	 * In the case 4 gamma is of the form
1900 	 *    \f[ \gamma(zlp + t \cdot \text{ray}) =
1901 	 *      \begin{cases}
1902 	 *        \sqrt{A t^2 + B t + C} - (D t + E), & \text{if some condition holds}, \\
1903 	 *        \sqrt{A' t^2 + B' t + C'} - (D' t + E'), & \text{otherwise.}
1904 	 *      \end{cases}
1905 	 *    \f]
1906 	 *
1907 	 * It can be shown (given the special properties of \f$\gamma\f$) that the smallest positive root of each function of the form
1908 	 * \f$\sqrt{a t^2 + b t + c} - (d t + e)\f$
1909 	 * is the same as the smallest positive root of the quadratic equation:
1910 	 * \f{align}{
1911 	 *     & \sqrt{a t^2 + b t + c} - (d t + e)) (\sqrt{a t^2 + b t + c} + (d t + e)) = 0 \\  \Leftrightarrow
1912 	 *     & (a - d^2) t^2 + (b - 2 d\,e) t + (c - e^2) = 0
1913 	 * \f}
1914 	 *
1915 	 * So, in cases 1, 2, and 3, this function just returns the solution of the above equation.
1916 	 * In case 4, it first solves the equation assuming we are in the first piece.
1917 	 * If there is no solution, then the second piece can't have a solution (first piece &ge; second piece for all t)
1918 	 * Then we check if the solution satisfies the condition.
1919 	 * If it doesn't then we solve the equation for the second piece.
1920 	 * If it has a solution, then it _has_ to be the solution.
1921 	 */
1922 	static
1923 	SCIP_Real computeIntersectionPoint(
1924 	   SCIP*                 scip,               /**< SCIP data structure */
1925 	   SCIP_NLHDLRDATA*      nlhdlrdata,         /**< nlhdlr data */
1926 	   SCIP_Bool             iscase4,            /**< whether we are in case 4 or not */
1927 	   SCIP_Real*            coefs1234a,         /**< values of A, B, C, D, and E of cases 1, 2, 3, or 4a */
1928 	   SCIP_Real*            coefs4b,            /**< values of A, B, C, D, and E of case 4b */
1929 	   SCIP_Real*            coefscondition      /**< values needed to evaluate condition of case 4 */
1930 	   )
1931 	{
1932 	   SCIP_Real sol1234a;
1933 	   SCIP_Real sol4b;
1934 	
1935 	   assert(coefs1234a != NULL);
1936 	
1937 	#ifdef  DEBUG_INTERSECTIONCUT
1938 	   SCIPinfoMessage(scip, NULL, "Computing intersection point for case 4? %d\n", iscase4);
1939 	#endif
1940 	   if( ! iscase4 )
1941 	      return computeRoot(scip, coefs1234a);
1942 	
1943 	   assert(coefs4b != NULL);
1944 	   assert(coefscondition != NULL);
1945 	
1946 	   /* compute solution of first piece */
1947 	#ifdef  DEBUG_INTERSECTIONCUT
1948 	   SCIPinfoMessage(scip, NULL, "Compute root in 4a\n");
1949 	#endif
1950 	   sol1234a = computeRoot(scip, coefs1234a);
1951 	
1952 	   /* if there is no solution --> second piece doesn't have solution */
1953 	   if( SCIPisInfinity(scip, sol1234a) )
1954 	   {
1955 	      /* this assert fails on multiplants_mtg5 the problem is that sqrt(A) <= D in 4a but not in 4b,
1956 	       * now, this is impossible since the phi4a >= phi4b, so actually sqrt(A) is 10e-15 away from
1957 	       * D in 4b
1958 	       */
1959 	      /* assert(SCIPisInfinity(scip, computeRoot(scip, coefs4b))); */
1960 	      return sol1234a;
1961 	   }
1962 	
1963 	   /* if solution of 4a is in 4a, then return */
1964 	   if( isCase4a(sol1234a, coefs1234a, coefscondition) )
1965 	      return sol1234a;
1966 	
1967 	#ifdef  DEBUG_INTERSECTIONCUT
1968 	   SCIPinfoMessage(scip, NULL, "Root not in 4a -> Compute root in 4b\n");
1969 	#endif
1970 	
1971 	   /* not on 4a --> then the intersection point is whatever 4b says: as phi4a >= phi4b, the solution of phi4b should
1972 	    * always be larger (but shouldn't be equal at this point given that isCase4a failed, and the condition function
1973 	    * evaluates to 0 when phi4a == phi4b) than the solution of phi4a; However, because of numerics (or limits in the
1974 	    * binary search) we can find a slightly smaller solution; thus, we just keep the larger one
1975 	    */
1976 	   sol4b = computeRoot(scip, coefs4b);
1977 	
1978 	   /* this assert fails in many instances, e.g. water, because sol4b < sol1234a  */
1979 	   /* assert(SCIPisInfinity(scip, sol4b) || !isCase4a(sol4b, coefs1234a, coefscondition)); */
1980 	   /* count number of times we could have improved the coefficient by 10% */
1981 	   if( sol4b < sol1234a && evalPhiAtRay(1.1 * sol1234a, coefs4b[0], coefs4b[1], coefs4b[2], coefs4b[3], coefs4b[4]) <=
1982 	         0.0 )
1983 	      nlhdlrdata->ncouldimprovedcoef++;
1984 	
1985 	   return MAX(sol1234a, sol4b);
1986 	}
1987 	
1988 	/** checks if numerics of the coefficients are not too bad */
1989 	static
1990 	SCIP_Bool areCoefsNumericsGood(
1991 	   SCIP*                 scip,               /**< SCIP data structure */
1992 	   SCIP_NLHDLRDATA*      nlhdlrdata,         /**< nlhdlr data */
1993 	   SCIP_Real*            coefs1234a,         /**< coefficients for case 1-3 and 4a */
1994 	   SCIP_Real*            coefs4b,            /**< coefficients for case 4b */
1995 	   SCIP_Bool             iscase4             /**< whether we are in case 4 */
1996 	   )
1997 	{
1998 	   SCIP_Real max;
1999 	   SCIP_Real min;
2000 	   int j;
2001 	
2002 	   /* check at phi at 0 is negative (note; this could be checked before restricting to the ray) also, if this
2003 	    * succeeds for one ray, it should suceed for every ray
2004 	    */
2005 	   if( SQRT( coefs1234a[2] ) - coefs1234a[4] >= 0.0 )
2006 	   {
2007 	      INTERLOG(printf("Bad numerics: phi(0) >= 0\n"); )
2008 	      nlhdlrdata->nphinonneg++;
2009 	      return FALSE;
2010 	   }
2011 	
2012 	   /* maybe we want to avoid a large dynamism between A, B and C */
2013 	   if( nlhdlrdata->ignorebadrayrestriction )
2014 	   {
2015 	      max = 0.0;
2016 	      min = SCIPinfinity(scip);
2017 	      for( j = 0; j < 3; ++j )
2018 	      {
2019 	         SCIP_Real absval;
2020 	
2021 	         absval = REALABS(coefs1234a[j]);
2022 	         if( max < absval )
2023 	            max = absval;
2024 	         if( absval != 0.0 && absval < min )
2025 	            min = absval;
2026 	      }
2027 	
2028 	      if( SCIPisHugeValue(scip, max / min) )
2029 	      {
2030 	         INTERLOG(printf("Bad numerics 1 2 3 or 4a: max(A,B,C)/min(A,B,C) is too large (%g)\n", max / min); )
2031 	         nlhdlrdata->nbadrayrestriction++;
2032 	         return FALSE;
2033 	      }
2034 	
2035 	      if( iscase4 )
2036 	      {
2037 	         max = 0.0;
2038 	         min = SCIPinfinity(scip);
2039 	         for( j = 0; j < 3; ++j )
2040 	         {
2041 	            SCIP_Real absval;
2042 	
2043 	            absval = ABS(coefs4b[j]);
2044 	            if( max < absval )
2045 	               max = absval;
2046 	            if( absval != 0.0 && absval < min )
2047 	               min = absval;
2048 	         }
2049 	
2050 	         if( SCIPisHugeValue(scip, max / min) )
2051 	         {
2052 	            INTERLOG(printf("Bad numeric 4b: max(A,B,C)/min(A,B,C) is too large (%g)\n", max / min); )
2053 	            nlhdlrdata->nbadrayrestriction++;
2054 	            return FALSE;
2055 	         }
2056 	      }
2057 	   }
2058 	
2059 	   return TRUE;
2060 	}
2061 	
2062 	
2063 	/** computes the coefficients a, b, c defining the quadratic function defining set S restricted to the line
2064 	 *  theta * apex.
2065 	 *
2066 	 *  The solution to the monoidal strengthening problem is then given by the smallest root of the function
2067 	 *  a * theta^2 + b * theta + c
2068 	 */
2069 	static
2070 	SCIP_RETCODE computeMonoidalQuadCoefs(
2071 	   SCIP*                 scip,               /**< SCIP data structure */
2072 	   SCIP_NLHDLREXPRDATA*  nlhdlrexprdata,     /**< nlhdlr expression data */
2073 	   SCIP_Real*            vb,                 /**< array containing \f$v_i^T b\f$ for \f$i \in I_+ \cup I_-\f$ */
2074 	   SCIP_Real*            vzlp,               /**< array containing \f$v_i^T zlp_q\f$ for \f$i \in I_+ \cup I_-\f$ */
2075 	   SCIP_Real*            vapex,              /**< array containing \f$v_i^T apex\f$ */
2076 	   SCIP_Real*            vray,               /**< array containing \f$v_i^T ray\f$ */
2077 	   SCIP_Real             kappa,              /**< value of kappa */
2078 	   SCIP_Real             sidefactor,         /**< 1.0 if the violated constraint is q &le; rhs, -1.0 otherwise */
2079 	   SCIP_Real*            a,                  /**< pointer to store quadratic coefficient */
2080 	   SCIP_Real*            b,                  /**< pointer to store linear coefficient */
2081 	   SCIP_Real*            c                   /**< pointer to store constant */
2082 	   )
2083 	{
2084 	   SCIP_EXPR* qexpr;
2085 	   int nquadexprs;
2086 	   SCIP_Real* eigenvectors;
2087 	   SCIP_Real* eigenvalues;
2088 	   int i;
2089 	
2090 	   qexpr = nlhdlrexprdata->qexpr;
2091 	   SCIPexprGetQuadraticData(qexpr, NULL, NULL, NULL, NULL, &nquadexprs, NULL, &eigenvalues, &eigenvectors);
2092 	
2093 	   *a = 0.0;
2094 	   *b = 0.0;
2095 	   *c = 0.0;
2096 	   for( i = 0; i < nquadexprs; ++i )
2097 	   {
2098 	      SCIP_Real dot;
2099 	
2100 	      if( SCIPisZero(scip, sidefactor * eigenvalues[i]) )
2101 	         continue;
2102 	
2103 	      dot = vray[i] + vapex[i] + vb[i] / (2.0 * sidefactor * eigenvalues[i]);
2104 	
2105 	      *a += sidefactor * eigenvalues[i] * SQR(vzlp[i] - vapex[i]);
2106 	      *b += sidefactor * eigenvalues[i] * (vzlp[i] - vapex[i]) * dot;
2107 	      *c += sidefactor * eigenvalues[i] * SQR(dot);
2108 	   }
2109 	
2110 	   *b *= 2.0;
2111 	   *c += kappa;
2112 	
2113 	   return SCIP_OKAY;
2114 	}
2115 	
2116 	/** check if ray was in strip by checking if the point in the monoid corresponding to the cutcoef we just found
2117 	 * is "on the wrong side" of the hyperplane -(a - lambda^Ta lambda)^T x
2118 	 */
2119 	static
2120 	SCIP_Bool isRayInStrip(
2121 	   SCIP_NLHDLREXPRDATA*  nlhdlrexprdata,     /**< nlhdlr expression data */
2122 	   SCIP_Real*            vb,                 /**< array containing \f$v_i^T b\f$ for \f$i \in I_+ \cup I_-\f$ */
2123 	   SCIP_Real*            vzlp,               /**< array containing \f$v_i^T zlp_q\f$ for \f$i \in I_+ \cup I_-\f$ */
2124 	   SCIP_Real*            vapex,              /**< array containing \f$v_i^T apex\f$ */
2125 	   SCIP_Real*            vray,               /**< array containing \f$v_i^T ray\f$ */
2126 	   SCIP_Real             kappa,              /**< value of kappa */
2127 	   SCIP_Real             sidefactor,         /**< 1.0 if the violated constraint is q &le; rhs, -1.0 otherwise */
2128 	   SCIP_Real             cutcoef             /**< optimal solution of the monoidal quadratic */
2129 	   )
2130 	{
2131 	   SCIP_EXPR* qexpr;
2132 	   int nquadexprs;
2133 	   SCIP_Real* eigenvectors;
2134 	   SCIP_Real* eigenvalues;
2135 	   SCIP_Real num;
2136 	   SCIP_Real denom;
2137 	   int i;
2138 	
2139 	   qexpr = nlhdlrexprdata->qexpr;
2140 	   SCIPexprGetQuadraticData(qexpr, NULL, NULL, NULL, NULL, &nquadexprs, NULL, &eigenvalues, &eigenvectors);
2141 	
2142 	   num = 0.0;
2143 	   denom = 0.0;
2144 	   for( i = 0; i < nquadexprs; ++i )
2145 	   {
2146 	      SCIP_Real dot;
2147 	
2148 	      if( sidefactor * eigenvalues[i] <= 0.0 )
2149 	         continue;
2150 	
2151 	      dot = vzlp[i] + vb[i] / (2.0 * sidefactor * eigenvalues[i]);
2152 	
2153 	      num += sidefactor * eigenvalues[i] * dot * (cutcoef * (vzlp[i] - vapex[i]) + vray[i] + vapex[i]);
2154 	      denom += sidefactor * eigenvalues[i] * SQR(dot);
2155 	   }
2156 	
2157 	   num /= kappa;
2158 	   num += 1.0;
2159 	   denom /= kappa;
2160 	   denom += 1.0;
2161 	
2162 	   assert(denom > 0);
2163 	
2164 	   return num / denom < 1;
2165 	}
2166 	
2167 	/** computes the smallest root of the quadratic function a*x^2 + b*x + c with a > 0
2168 	 *  and b^2 - 4ac >= 0. We use binary search between -inf and minimum at -b/2a.
2169 	 */
2170 	static
2171 	SCIP_Real findMonoidalQuadRoot(
2172 	   SCIP*                 scip,               /**< SCIP data structure */
2173 	   SCIP_Real             a,                  /**< quadratic coefficient */
2174 	   SCIP_Real             b,                  /**< linear coefficient */
2175 	   SCIP_Real             c                   /**< constant */
2176 	   )
2177 	{
2178 	   SCIP_Real sol;
2179 	   SCIP_INTERVAL bounds;
2180 	   SCIP_INTERVAL result;
2181 	
2182 	   assert(SQR(b) - 4 * a * c >= 0.0);
2183 	
2184 	   if( SCIPisZero(scip, a) )
2185 	   {
2186 	      assert(b != 0.0);
2187 	      sol = - c / b;
2188 	   }
2189 	   else
2190 	   {
2191 	      SCIPintervalSetBounds(&bounds, - b / (2.0 * a), SCIPinfinity(scip));
2192 	
2193 	      /* find all positive x such that a x^2 + b x >= -c and x in bounds.*/
2194 	      SCIPintervalSolveUnivariateQuadExpressionPositiveAllScalar(SCIP_INTERVAL_INFINITY, &result, a, b, -c, bounds);
2195 	      sol = SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, result) ? SCIPinfinity(scip) : SCIPintervalGetInf(result);
2196 	
2197 	      /* if we didn't find any positive solutions, negate quadratic and find negative solutions */
2198 	      if( SCIPisInfinity(scip, sol) )
2199 	      {
2200 	         SCIPintervalSetBounds(&bounds, b / (2.0 * a), SCIPinfinity(scip));
2201 	
2202 	         /* find all positive x such that a x^2 - b x >= -c and x in bounds.*/
2203 	         SCIPintervalSolveUnivariateQuadExpressionPositiveAllScalar(SCIP_INTERVAL_INFINITY, &result, a, -b, -c, bounds);
2204 	         sol = SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, result) ? SCIPinfinity(scip) : -SCIPintervalGetInf(result);
2205 	      }
2206 	   }
2207 	
2208 	   /* check if that solution is close enough or if we need to improve it more with binary search */
2209 	   if( a * SQR(sol) + sol * b + c > 1e-10 )
2210 	   {
2211 	      SCIP_Real val;
2212 	      SCIP_Real lb;
2213 	      SCIP_Real ub;
2214 	      SCIP_Real lastposval;
2215 	      SCIP_Real lastpossol;
2216 	      int niter;
2217 	
2218 	      lb = - b / (2.0 * a);
2219 	      ub = sol;
2220 	      lastposval = SCIPinfinity(scip);
2221 	      lastpossol = SCIPinfinity(scip);
2222 	      val = SCIPinfinity(scip);
2223 	      niter = 0;
2224 	      while( niter < BINSEARCH_MAXITERS && ABS(val) > 1e-10 )
2225 	      {
2226 	         sol = (ub + lb) / 2.0;
2227 	         val = a * SQR(sol) + b * sol + c;
2228 	
2229 	         if( val < 0.0 )
2230 	            lb = sol;
2231 	         else
2232 	            ub = sol;
2233 	
2234 	         /* if we are close enough, return with (feasible) solution */
2235 	         if( val > 0.0 && val < 1e-6 )
2236 	            break;
2237 	
2238 	         if( val > 0.0 && lastposval > val )
2239 	         {
2240 	            lastposval = val;
2241 	            lastpossol = sol;
2242 	         }
2243 	
2244 	         ++niter;
2245 	      }
2246 	      if( val < 0.0 && ! SCIPisZero(scip, val) )
2247 	      {
2248 	         sol = lastpossol;
2249 	         val = lastposval;
2250 	      }
2251 	
2252 	      assert( val > 0.0 || SCIPisZero(scip, val) );
2253 	   }
2254 	
2255 	   return sol;
2256 	}
2257 	
2258 	/** computes the apex of the S-free set (if it exists) */
2259 	static
2260 	void computeApex(
2261 	   SCIP_NLHDLREXPRDATA*  nlhdlrexprdata,     /**< nlhdlr expression data */
2262 	   SCIP_Real*            vb,                 /**< array containing \f$v_i^T b\f$ for \f$i \in I_+ \cup I_-\f$ */
2263 	   SCIP_Real*            vzlp,               /**< array containing \f$v_i^T zlp_q\f$ for \f$i \in I_+ \cup I_-\f$ */
2264 	   SCIP_Real             kappa,              /**< value of kappa */
2265 	   SCIP_Real             sidefactor,         /**< 1.0 if the violated constraint is q &le; rhs, -1.0 otherwise */
2266 	   SCIP_Real*            apex,               /**< array to store apex */
2267 	   SCIP_Bool*            success             /**< TRUE if computation of apex was successful */
2268 	   )
2269 	{
2270 	   SCIP_EXPR* qexpr;
2271 	   int nquadexprs;
2272 	   SCIP_Real* eigenvectors;
2273 	   SCIP_Real* eigenvalues;
2274 	   int i;
2275 	
2276 	   qexpr = nlhdlrexprdata->qexpr;
2277 	   SCIPexprGetQuadraticData(qexpr, NULL, NULL, NULL, NULL, &nquadexprs, NULL, &eigenvalues, &eigenvectors);
2278 	
2279 	   *success = TRUE;
2280 	
2281 	   for( i = 0; i < nquadexprs; ++i )
2282 	   {
2283 	      SCIP_Real entry;
2284 	      SCIP_Real denom;
2285 	      SCIP_Real num;
2286 	      int j;
2287 	
2288 	      entry = 0.0;
2289 	      num = 0.0;
2290 	      denom = 0.0;
2291 	      for( j = 0; j < nquadexprs; ++j )
2292 	      {
2293 	         if( sidefactor * eigenvalues[j] == 0.0 )
2294 	            continue;
2295 	
2296 	         entry -= eigenvectors[j * nquadexprs + i] * vb[j] / (2.0 * sidefactor * eigenvalues[j]);
2297 	
2298 	         if( sidefactor * eigenvalues[j] > 0.0 )
2299 	         {
2300 	            SCIP_Real dot;
2301 	
2302 	            dot = vzlp[j] + vb[j] / (2.0 * sidefactor * eigenvalues[j]);
2303 	
2304 	            num += eigenvectors[j * nquadexprs + i] * dot;
2305 	            denom += sidefactor * eigenvalues[j] * SQR(dot);
2306 	         }
2307 	      }
2308 	
2309 	      /* if denom = 0, the S-free set is just the strip, so we don't have an apex -> monoidal not possible */
2310 	      if( denom == 0.0 )
2311 	      {
2312 	         *success = FALSE;
2313 	         return;
2314 	      }
2315 	      assert(denom > 0.0);
2316 	
2317 	      num *= -kappa;
2318 	      entry += num / denom;
2319 	
2320 	      apex[i] = entry;
2321 	   }
2322 	}
2323 	
2324 	/** for a given ray, computes the cut coefficient using monoidal strengthening (if possible) */
2325 	static
2326 	SCIP_RETCODE computeMonoidalStrengthCoef(
2327 	   SCIP*                 scip,               /**< SCIP data structure */
2328 	   SCIP_NLHDLREXPRDATA*  nlhdlrexprdata,     /**< nlhdlr expression data */
2329 	   int                   lppos,              /**< lp pos of current ray */
2330 	   SCIP_Real*            raycoefs,           /**< coefficients of ray */
2331 	   int*                  rayidx,             /**< index of consvar the ray coef is associated to */
2332 	   int                   raynnonz,           /**< length of raycoefs and rayidx */
2333 	   SCIP_Real*            vb,                 /**< array containing \f$v_i^T b\f$ for \f$i \in I_+ \cup I_-\f$ */
2334 	   SCIP_Real*            vzlp,               /**< array containing \f$v_i^T zlp_q\f$ for \f$i \in I_+ \cup I_-\f$ */
2335 	   SCIP_Real*            wcoefs,             /**< coefficients of w for the qud vars or NULL if w is 0 */
2336 	   SCIP_Real             kappa,              /**< value of kappa */
2337 	   SCIP_Real*            apex,               /**< array containing the apex of the S-free set in the original space */
2338 	   SCIP_Real             sidefactor,         /**< 1.0 if the violated constraint is q &le; rhs, -1.0 otherwise */
2339 	   SCIP_Real*            cutcoef,            /**< pointer to store cut coef */
2340 	   SCIP_Bool*            success             /**< TRUE if monoidal strengthening could be applied */
2341 	   )
2342 	{
2343 	   SCIP_COL** cols;
2344 	   SCIP_ROW** rows;
2345 	
2346 	   *success = FALSE;
2347 	
2348 	   /* check if we are in the correct case (case 2) */
2349 	   assert(wcoefs == NULL && kappa > 0.0);
2350 	
2351 	   cols = SCIPgetLPCols(scip);
2352 	   rows = SCIPgetLPRows(scip);
2353 	
2354 	   /* if var corresponding to current ray is integer, we can do monoidal */
2355 	   if( (lppos >= 0 && SCIPvarGetType(SCIPcolGetVar(cols[lppos])) != SCIP_VARTYPE_CONTINUOUS) ||
2356 	       (lppos < 0 && SCIProwIsIntegral(rows[- lppos - 1])) )
2357 	   {
2358 	      SCIP_Real* vapex;
2359 	      SCIP_Real* vray;
2360 	      SCIP_Real a;
2361 	      SCIP_Real b;
2362 	      SCIP_Real c;
2363 	      int nquadexprs;
2364 	
2365 	      SCIPexprGetQuadraticData(nlhdlrexprdata->qexpr, NULL, NULL, NULL, NULL, &nquadexprs, NULL, NULL, NULL);
2366 	
2367 	      /* compute v_i^T apex in vapex[i] and v_i^T ray in vray[i] */
2368 	      SCIP_CALL( SCIPallocBufferArray(scip, &vapex, nquadexprs) );
2369 	      SCIP_CALL( SCIPallocBufferArray(scip, &vray, nquadexprs) );
2370 	      computeVApexAndVRay(nlhdlrexprdata, apex, raycoefs, rayidx, raynnonz, vapex, vray);
2371 	
2372 	      /* compute coefficients of the quadratic monoidal problem function */
2373 	      SCIP_CALL( computeMonoidalQuadCoefs(scip, nlhdlrexprdata, vb, vzlp, vapex, vray, kappa,
2374 	         sidefactor, &a, &b, &c) );
2375 	
2376 	      /* check if ray is in strip */
2377 	      if( SQR(b) - (4 * a * c) >= 0.0 )
2378 	      {
2379 	         /* find smallest root of quadratic function a * x^2 + b * x + c -> this is the cut coef */
2380 	         *cutcoef = findMonoidalQuadRoot(scip, a, b, c);
2381 	
2382 	         /* if the cutcoef is negative, we have to set it to 0 */
2383 	         *cutcoef = MAX(*cutcoef, 0.0);
2384 	
2385 	         /* check if ray is in strip. If not, monoidal is possible and cutcoef is the strengthened cut coef */
2386 	         if( ! isRayInStrip(nlhdlrexprdata, vb, vzlp, vapex, vray, kappa, sidefactor, *cutcoef) )
2387 	         {
2388 	            *success = TRUE;
2389 	         }
2390 	      }
2391 	
2392 	      SCIPfreeBufferArray(scip, &vray);
2393 	      SCIPfreeBufferArray(scip, &vapex);
2394 	   }
2395 	
2396 	   return SCIP_OKAY;
2397 	}
2398 	
2399 	/** sparsify intersection cut by replacing non-basic variables with their bounds if their coefficient allows it */
2400 	static
2401 	void sparsifyIntercut(
2402 	   SCIP*                 scip,               /**< SCIP data structure */
2403 	   SCIP_ROWPREP*         rowprep             /**< rowprep for the generated cut */
2404 	   )
2405 	{
2406 	   int i;
2407 	   int nvars;
2408 	
2409 	   /* get number of variables in rowprep */
2410 	   nvars = SCIProwprepGetNVars(rowprep);
2411 	
2412 	   /* go though all the variables in rowprep */
2413 	   for( i = 0; i < nvars; ++i )
2414 	   {
2415 	      SCIP_VAR* var;
2416 	      SCIP_Real coef;
2417 	      SCIP_Real lb;
2418 	      SCIP_Real ub;
2419 	      SCIP_Real solval;
2420 	
2421 	      /* get variable and its coefficient */
2422 	      var = SCIProwprepGetVars(rowprep)[i];
2423 	      coef = SCIProwprepGetCoefs(rowprep)[i];
2424 	
2425 	      /* get bounds of variable */
2426 	      lb = SCIPvarGetLbLocal(var);
2427 	      ub = SCIPvarGetUbLocal(var);
2428 	
2429 	      /* get LP solution value of variable */
2430 	      solval = SCIPgetSolVal(scip, NULL, var);
2431 	
2432 	      /* if the variable is at its lower or upper bound and the coefficient has the correct sign, we can
2433 	       * set the cutcoef to 0
2434 	       */
2435 	      if( SCIPisZero(scip, ub - solval) && coef > 0.0 )
2436 	      {
2437 	         SCIProwprepAddSide(rowprep, -coef * ub);
2438 	         SCIProwprepSetCoef(rowprep, i, 0.0);
2439 	      }
2440 	      else if( SCIPisZero(scip, solval - lb) && coef < 0.0 )
2441 	      {
2442 	         SCIProwprepAddSide(rowprep, -coef * lb);
2443 	         SCIProwprepSetCoef(rowprep, i, 0.0);
2444 	      }
2445 	   }
2446 	
2447 	   return;
2448 	}
2449 	
2450 	/** computes intersection cut cuts off sol (because solution sol violates the quadratic constraint cons)
2451 	 * and stores it in rowprep. Here, we don't use any strengthening.
2452 	 */
2453 	static
2454 	SCIP_RETCODE computeIntercut(
2455 	   SCIP*                 scip,               /**< SCIP data structure */
2456 	   SCIP_NLHDLRDATA*      nlhdlrdata,         /**< nlhdlr data */
2457 	   SCIP_NLHDLREXPRDATA*  nlhdlrexprdata,     /**< nlhdlr expression data */
2458 	   RAYS*                 rays,               /**< rays */
2459 	   SCIP_Real             sidefactor,         /**< 1.0 if the violated constraint is q &le; rhs, -1.0 otherwise */
2460 	   SCIP_Bool             iscase4,            /**< whether we are in case 4 */
2461 	   SCIP_Real*            vb,                 /**< array containing \f$v_i^T b\f$ for \f$i \in I_+ \cup I_-\f$ */
2462 	   SCIP_Real*            vzlp,               /**< array containing \f$v_i^T zlp_q\f$ for \f$i \in I_+ \cup I_-\f$ */
2463 	   SCIP_Real*            wcoefs,             /**< coefficients of w for the qud vars or NULL if w is 0 */
2464 	   SCIP_Real             wzlp,               /**< value of w at zlp */
2465 	   SCIP_Real             kappa,              /**< value of kappa */
2466 	   SCIP_ROWPREP*         rowprep,            /**< rowprep for the generated cut */
2467 	   SCIP_Real*            interpoints,        /**< array to store intersection points for all rays or NULL if nothing
2468 	                                              *   needs to be stored */
2469 	   SCIP_SOL*             sol,                /**< solution we want to separate */
2470 	   SCIP_Bool*            success             /**< if a cut candidate could be computed */
2471 	   )
2472 	{
2473 	   SCIP_COL** cols;
2474 	   SCIP_ROW** rows;
2475 	   SCIP_Real* apex;
2476 	   SCIP_Real avecutcoefsum;
2477 	   SCIP_Real avemonoidalimprovsum;
2478 	   int monoidalcounter;
2479 	   int counter;
2480 	   int i;
2481 	   SCIP_Bool usemonoidal;
2482 	   SCIP_Bool monoidalwasused;
2483 	   SCIP_Bool case2;
2484 	
2485 	   cols = SCIPgetLPCols(scip);
2486 	   rows = SCIPgetLPRows(scip);
2487 	
2488 	   case2 = wcoefs == NULL && kappa > 0.0;
2489 	   monoidalwasused = FALSE;
2490 	
2491 	   counter = 0;
2492 	   monoidalcounter = 0;
2493 	   avecutcoefsum = 0.0;
2494 	   avemonoidalimprovsum = 0.0;
2495 	
2496 	   /* if we use monoidal and we are in the right case for it, compute the apex of the S-free set */
2497 	   if( case2 )
2498 	   {
2499 	      int nquadexprs;
2500 	
2501 	      SCIPexprGetQuadraticData(nlhdlrexprdata->qexpr, NULL, NULL, NULL, NULL, &nquadexprs, NULL, NULL, NULL);
2502 	
2503 	      /* allocate memory for apex */
2504 	      SCIP_CALL( SCIPallocBufferArray(scip, &apex, nquadexprs) );
2505 	
2506 	      computeApex(nlhdlrexprdata, vb, vzlp, kappa, sidefactor, apex, success);
2507 	
2508 	      /* if computation of apex was not successful, don't apply monoidal strengthening */
2509 	      if( ! *success )
2510 	         case2 = FALSE;
2511 	   }
2512 	   else
2513 	   {
2514 	      apex = NULL;
2515 	   }
2516 	
2517 	   usemonoidal = nlhdlrdata->usemonoidal && case2;
2518 	
2519 	   /* for every ray: compute cut coefficient and add var associated to ray into cut */
2520 	   for( i = 0; i < rays->nrays; ++i )
2521 	   {
2522 	      SCIP_Real interpoint;
2523 	      SCIP_Real cutcoef;
2524 	      int lppos;
2525 	      SCIP_Real coefs1234a[5];
2526 	      SCIP_Real coefs4b[5];
2527 	      SCIP_Real coefscondition[3];
2528 	      SCIP_Bool monoidalsuccess;
2529 	
2530 	      monoidalsuccess = FALSE;
2531 	      cutcoef = SCIPinfinity(scip);
2532 	
2533 	      /* if we (can) use monoidal strengthening, compute the cut coefficient with that */
2534 	      if( usemonoidal )
2535 	      {
2536 	         SCIP_CALL( computeMonoidalStrengthCoef(scip, nlhdlrexprdata, rays->lpposray[i], &rays->rays[rays->raysbegin[i]],
2537 	               &rays->raysidx[rays->raysbegin[i]], rays->raysbegin[i + 1] - rays->raysbegin[i], vb, vzlp, wcoefs, kappa,
2538 	               apex, sidefactor, &cutcoef, &monoidalsuccess) );
2539 	      }
2540 	
2541 	      /* track whether monoidal was successful at least once if it is used */
2542 	      if( usemonoidal && ! monoidalwasused && monoidalsuccess )
2543 	         monoidalwasused = TRUE;
2544 	
2545 	      /* if we don't use monoidal or if monoidal couldn't be applied, use gauge to compute coef  */
2546 	      if( ! usemonoidal || ! monoidalsuccess || nlhdlrdata->trackmore )
2547 	      {
2548 	         /* restrict phi to ray */
2549 	         SCIP_CALL( computeRestrictionToRay(scip, nlhdlrexprdata, sidefactor, iscase4,
2550 	                  &rays->rays[rays->raysbegin[i]], &rays->raysidx[rays->raysbegin[i]], rays->raysbegin[i + 1] -
2551 	                  rays->raysbegin[i], vb, vzlp, wcoefs, wzlp, kappa, coefs1234a, coefs4b, coefscondition, success) );
2552 	
2553 	         if( ! *success )
2554 	            goto TERMINATE;
2555 	
2556 	         /* if restriction to ray is numerically nasty -> abort cut separation */
2557 	         *success = areCoefsNumericsGood(scip, nlhdlrdata, coefs1234a, coefs4b, iscase4);
2558 	
2559 	         if( ! *success )
2560 	            goto TERMINATE;
2561 	
2562 	         /* compute intersection point */
2563 	         interpoint = computeIntersectionPoint(scip, nlhdlrdata, iscase4, coefs1234a, coefs4b, coefscondition);
2564 	
2565 	#ifdef  DEBUG_INTERSECTIONCUT
2566 	         SCIPinfoMessage(scip, NULL, "interpoint for ray %d is %g\n", i, interpoint);
2567 	#endif
2568 	
2569 	         /* store intersection point */
2570 	         if( interpoints != NULL )
2571 	            interpoints[i] = interpoint;
2572 	
2573 	         /* if we are only computing the "normal" cut coefficient to track statistics (and we have been successful
2574 	          * with monoidal strengthening), don't overwrite cutcoef variable, but just track statistics */
2575 	         if( nlhdlrdata->trackmore && monoidalsuccess )
2576 	         {
2577 	            SCIP_Real normalcutcoef;
2578 	
2579 	            normalcutcoef = SCIPisInfinity(scip, interpoint) ? 0.0 : 1.0 / interpoint;
2580 	
2581 	            if( cutcoef >= 0 )
2582 	               avemonoidalimprovsum += cutcoef / normalcutcoef;
2583 	            ++monoidalcounter;
2584 	         }
2585 	         else
2586 	         {
2587 	            /* compute cut coef */
2588 	            cutcoef = SCIPisInfinity(scip, interpoint) ? 0.0 : 1.0 / interpoint;
2589 	         }
2590 	
2591 	         if( cutcoef == 0.0 && case2 && nlhdlrdata->useminrep )
2592 	         {
2593 	            /* restrict phi to the line defined by ray + apex + t*(f - apex) */
2594 	            SCIP_CALL( computeRestrictionToLine(scip, nlhdlrexprdata, sidefactor,
2595 	               &rays->rays[rays->raysbegin[i]], &rays->raysidx[rays->raysbegin[i]], rays->raysbegin[i + 1] -
2596 	               rays->raysbegin[i], vb, vzlp, kappa, apex, coefs1234a, success) );
2597 	
2598 	            if( ! *success )
2599 	               goto TERMINATE;
2600 	
2601 	            /* if restriction to ray is numerically nasty -> abort cut separation */
2602 	            *success = areCoefsNumericsGood(scip, nlhdlrdata, coefs1234a, NULL, FALSE);
2603 	
2604 	            if( ! *success )
2605 	               goto TERMINATE;
2606 	
2607 	            coefs1234a[1] *= -1.0;
2608 	            coefs1234a[3] *= -1.0;
2609 	
2610 	            /* compute intersection point */
2611 	            cutcoef = - computeRoot(scip, coefs1234a);
2612 	
2613 	            assert(cutcoef <= 0.0);
2614 	         }
2615 	      }
2616 	
2617 	      /* keep track of average cut coefficient */
2618 	      ++counter;
2619 	      avecutcoefsum += cutcoef;
2620 	      assert( ! SCIPisInfinity(scip, cutcoef) );
2621 	
2622 	      /* add var to cut: if variable is nonbasic at upper we have to flip sign of cutcoef */
2623 	      lppos = rays->lpposray[i];
2624 	      if( lppos < 0 )
2625 	      {
2626 	         lppos = -lppos - 1;
2627 	
2628 	         assert(SCIProwGetBasisStatus(rows[lppos]) == SCIP_BASESTAT_LOWER || SCIProwGetBasisStatus(rows[lppos]) ==
2629 	               SCIP_BASESTAT_UPPER);
2630 	
2631 	         /* flip cutcoef when necessary. Note: rows have flipped base status! */
2632 	         cutcoef = SCIProwGetBasisStatus(rows[lppos]) == SCIP_BASESTAT_UPPER ? cutcoef : -cutcoef;
2633 	
2634 	         SCIP_CALL( addRowToCut(scip, rowprep, cutcoef, rows[lppos], success) );
2635 	
2636 	         if( ! *success )
2637 	         {
2638 	            INTERLOG(printf("Bad numeric: now not nonbasic enough\n");)
2639 	            nlhdlrdata->nbadnonbasic++;
2640 	            goto TERMINATE;
2641 	         }
2642 	      }
2643 	      else
2644 	      {
2645 	         if( ! nlhdlrdata->useboundsasrays )
2646 	         {
2647 	            assert(SCIPcolGetBasisStatus(cols[lppos]) == SCIP_BASESTAT_UPPER || SCIPcolGetBasisStatus(cols[lppos]) ==
2648 	                  SCIP_BASESTAT_LOWER);
2649 	
2650 	            /* flip cutcoef when necessary */
2651 	            cutcoef = SCIPcolGetBasisStatus(cols[lppos]) == SCIP_BASESTAT_UPPER ? -cutcoef : cutcoef;
2652 	
2653 	            SCIP_CALL( addColToCut(scip, rowprep, sol, cutcoef, cols[lppos]) );
2654 	         }
2655 	         else
2656 	         {
2657 	            /* flip cutcoef when necessary */
2658 	            cutcoef = rays->rays[i] == -1 ? -cutcoef : cutcoef;
2659 	
2660 	            SCIP_CALL( addColToCut(scip, rowprep, sol, cutcoef, cols[lppos]) );
2661 	         }
2662 	      }
2663 	   }
2664 	
2665 	   /* track statistics */
2666 	   if( counter > 0 )
2667 	      nlhdlrdata->currentavecutcoef = avecutcoefsum / counter;
2668 	   if( monoidalwasused )
2669 	      nlhdlrdata->nmonoidal += 1;
2670 	   if( monoidalcounter > 0 )
2671 	      nlhdlrdata->currentavemonoidalimprovement = avemonoidalimprovsum / monoidalcounter;
2672 	
2673 	TERMINATE:
2674 	   SCIPfreeBufferArrayNull(scip, &apex);
2675 	
2676 	   return SCIP_OKAY;
2677 	}
2678 	
2679 	/** combine ray 1 and 2 to obtain new ray coef1 * ray1 - coef2 * ray2 in sparse format */
2680 	static
2681 	void combineRays(
2682 	   SCIP_Real*            raycoefs1,          /**< coefficients of ray 1 */
2683 	   int*                  rayidx1,            /**< index of consvar of the ray 1 coef is associated to */
2684 	   int                   raynnonz1,          /**< length of raycoefs1 and rayidx1 */
2685 	   SCIP_Real*            raycoefs2,          /**< coefficients of ray 2 */
2686 	   int*                  rayidx2,            /**< index of consvar of the ray 2 coef is associated to */
2687 	   int                   raynnonz2,          /**< length of raycoefs2 and rayidx2 */
2688 	   SCIP_Real*            newraycoefs,        /**< coefficients of combined ray */
2689 	   int*                  newrayidx,          /**< index of consvar of the combined ray coef is associated to */
2690 	   int*                  newraynnonz,        /**< pointer to length of newraycoefs and newrayidx */
2691 	   SCIP_Real             coef1,              /**< coef of ray 1 */
2692 	   SCIP_Real             coef2               /**< coef of ray 2 */
2693 	   )
2694 	{
2695 	   int idx1;
2696 	   int idx2;
2697 	
2698 	   idx1 = 0;
2699 	   idx2 = 0;
2700 	   *newraynnonz = 0;
2701 	
2702 	   while( idx1 < raynnonz1 || idx2 < raynnonz2 )
2703 	   {
2704 	      /* if the pointers look at different variables (or one already arrieved at the end), only one pointer can move
2705 	       * on
2706 	       */
2707 	      if( idx1 >= raynnonz1 || (idx2 < raynnonz2 && rayidx1[idx1] > rayidx2[idx2]) )
2708 	      {
2709 	         /*printf("case 1 \n"); */
2710 	         newraycoefs[*newraynnonz] = - coef2 * raycoefs2[idx2];
2711 	         newrayidx[*newraynnonz] = rayidx2[idx2];
2712 	         ++(*newraynnonz);
2713 	         ++idx2;
2714 	      }
2715 	      else if( idx2 >= raynnonz2 || rayidx1[idx1] < rayidx2[idx2] )
2716 	      {
2717 	         /*printf("case 2 \n"); */
2718 	         newraycoefs[*newraynnonz] = coef1 * raycoefs1[idx1];
2719 	         newrayidx[*newraynnonz] = rayidx1[idx1];
2720 	         ++(*newraynnonz);
2721 	         ++idx1;
2722 	      }
2723 	      /* if both pointers look at the same variable, just compute the difference and move both pointers */
2724 	      else if( rayidx1[idx1] == rayidx2[idx2] )
2725 	      {
2726 	         /*printf("case 3 \n"); */
2727 	         newraycoefs[*newraynnonz] = coef1 * raycoefs1[idx1] - coef2 * raycoefs2[idx2];
2728 	         newrayidx[*newraynnonz] = rayidx1[idx1];
2729 	         ++(*newraynnonz);
2730 	         ++idx1;
2731 	         ++idx2;
2732 	      }
2733 	   }
2734 	}
2735 	
2736 	/** checks if two rays are linearly dependent */
2737 	static
2738 	SCIP_Bool raysAreDependent(
2739 	   SCIP*                 scip,               /**< SCIP data structure */
2740 	   SCIP_Real*            raycoefs1,          /**< coefficients of ray 1 */
2741 	   int*                  rayidx1,            /**< index of consvar of the ray 1 coef is associated to */
2742 	   int                   raynnonz1,          /**< length of raycoefs1 and rayidx1 */
2743 	   SCIP_Real*            raycoefs2,          /**< coefficients of ray 2 */
2744 	   int*                  rayidx2,            /**< index of consvar of the ray 2 coef is associated to */
2745 	   int                   raynnonz2,          /**< length of raycoefs2 and rayidx2 */
2746 	   SCIP_Real*            coef                /**< pointer to store coef (s.t. r1 = coef * r2) in case rays are
2747 	                                              *   dependent */
2748 	   )
2749 	{
2750 	   int i;
2751 	
2752 	   /* cannot be dependent if they have different number of non-zero entries */
2753 	   if( raynnonz1 != raynnonz2 )
2754 	      return FALSE;
2755 	
2756 	   *coef = 0.0;
2757 	
2758 	   for( i = 0; i < raynnonz1; ++i )
2759 	   {
2760 	      /* cannot be dependent if different variables have non-zero entries */
2761 	      if( rayidx1[i] != rayidx2[i] ||
2762 	         (SCIPisZero(scip, raycoefs1[i]) && !SCIPisZero(scip, raycoefs2[i])) ||
2763 	         (!SCIPisZero(scip, raycoefs1[i]) && SCIPisZero(scip, raycoefs2[i])) )
2764 	      {
2765 	         return FALSE;
2766 	      }
2767 	
2768 	      if( *coef != 0.0 )
2769 	      {
2770 	         /* cannot be dependent if the coefs aren't equal for all entries */
2771 	         if( ! SCIPisEQ(scip, *coef, raycoefs1[i] / raycoefs2[i]) )
2772 	         {
2773 	            return FALSE;
2774 	         }
2775 	      }
2776 	      else
2777 	         *coef = raycoefs1[i] / raycoefs2[i];
2778 	   }
2779 	
2780 	   return TRUE;
2781 	}
2782 	
2783 	/** checks if the ray alpha * ray_i + (1 - alpha) * ray_j is in the recession cone of the S-free set. To do so,
2784 	  * we check if phi restricted to the ray has a positive root.
2785 	  */
2786 	static
2787 	SCIP_RETCODE rayInRecessionCone(
2788 	   SCIP*                 scip,               /**< SCIP data structure */
2789 	   SCIP_NLHDLRDATA*      nlhdlrdata,         /**< nlhdlr data */
2790 	   SCIP_NLHDLREXPRDATA*  nlhdlrexprdata,     /**< nlhdlr expression data */
2791 	   RAYS*                 rays,               /**< rays */
2792 	   int                   j,                  /**< index of current ray in recession cone */
2793 	   int                   i,                  /**< index of current ray not in recession cone */
2794 	   SCIP_Real             sidefactor,         /**< 1.0 if the violated constraint is q &le; rhs, -1.0 otherwise */
2795 	   SCIP_Bool             iscase4,            /**< whether we are in case 4 */
2796 	   SCIP_Real*            vb,                 /**< array containing \f$v_i^T b\f$ for \f$i \in I_+ \cup I_-\f$ */
2797 	   SCIP_Real*            vzlp,               /**< array containing \f$v_i^T zlp_q\f$ for \f$i \in I_+ \cup I_-\f$ */
2798 	   SCIP_Real*            wcoefs,             /**< coefficients of w for the quad vars or NULL if w is 0 */
2799 	   SCIP_Real             wzlp,               /**< value of w at zlp */
2800 	   SCIP_Real             kappa,              /**< value of kappa */
2801 	   SCIP_Real             alpha,              /**< coef for combining the two rays */
2802 	   SCIP_Bool*            inreccone,          /**< pointer to store whether the ray is in the recession cone or not */
2803 	   SCIP_Bool*            success             /**< Did numerical troubles occur? */
2804 	   )
2805 	{
2806 	   SCIP_Real coefs1234a[5];
2807 	   SCIP_Real coefs4b[5];
2808 	   SCIP_Real coefscondition[3];
2809 	   SCIP_Real interpoint;
2810 	   SCIP_Real* newraycoefs;
2811 	   int* newrayidx;
2812 	   int newraynnonz;
2813 	
2814 	  *inreccone = FALSE;
2815 	
2816 	   /* allocate memory for new ray */
2817 	   newraynnonz = (rays->raysbegin[i + 1] - rays->raysbegin[i]) + (rays->raysbegin[j + 1] - rays->raysbegin[j]);
2818 	   SCIP_CALL( SCIPallocBufferArray(scip, &newraycoefs, newraynnonz) );
2819 	   SCIP_CALL( SCIPallocBufferArray(scip, &newrayidx, newraynnonz) );
2820 	
2821 	   /* build the ray alpha * ray_i + (1 - alpha) * ray_j */
2822 	   combineRays(&rays->rays[rays->raysbegin[i]], &rays->raysidx[rays->raysbegin[i]], rays->raysbegin[i + 1] -
2823 	           rays->raysbegin[i], &rays->rays[rays->raysbegin[j]], &rays->raysidx[rays->raysbegin[j]],
2824 	           rays->raysbegin[j + 1] - rays->raysbegin[j], newraycoefs, newrayidx, &newraynnonz, alpha,
2825 	           -1 + alpha);
2826 	
2827 	   /* restrict phi to the "new" ray */
2828 	   SCIP_CALL( computeRestrictionToRay(scip, nlhdlrexprdata, sidefactor, iscase4, newraycoefs, newrayidx,
2829 	           newraynnonz, vb, vzlp, wcoefs, wzlp, kappa, coefs1234a, coefs4b, coefscondition, success) );
2830 	
2831 	   if( ! *success )
2832 	      goto CLEANUP;
2833 	
2834 	   /* check if restriction to "new" ray is numerically nasty. If so, treat the corresponding rho as if phi is
2835 	    * positive
2836 	    */
2837 	
2838 	   /* compute intersection point */
2839 	   interpoint = computeIntersectionPoint(scip, nlhdlrdata, iscase4, coefs1234a, coefs4b, coefscondition);
2840 	
2841 	   /* no root exists */
2842 	   if( SCIPisInfinity(scip, interpoint) )
2843 	      *inreccone = TRUE;
2844 	
2845 	CLEANUP:
2846 	   SCIPfreeBufferArray(scip, &newrayidx);
2847 	   SCIPfreeBufferArray(scip, &newraycoefs);
2848 	
2849 	   return SCIP_OKAY;
2850 	}
2851 	
2852 	/** finds the smallest negative steplength for the current ray r_idx such that the combination
2853 	 * of r_idx with all rays not in the recession cone is in the recession cone
2854 	 */
2855 	static
2856 	SCIP_RETCODE findRho(
2857 	   SCIP*                 scip,               /**< SCIP data structure */
2858 	   SCIP_NLHDLRDATA*      nlhdlrdata,         /**< nlhdlr data */
2859 	   SCIP_NLHDLREXPRDATA*  nlhdlrexprdata,     /**< nlhdlr expression data */
2860 	   RAYS*                 rays,               /**< rays */
2861 	   int                   idx,                /**< index of current ray we want to find rho for */
2862 	   SCIP_Real             sidefactor,         /**< 1.0 if the violated constraint is q &le; rhs, -1.0 otherwise */
2863 	   SCIP_Bool             iscase4,            /**< whether we are in case 4 */
2864 	   SCIP_Real*            vb,                 /**< array containing \f$v_i^T b\f$ for \f$i \in I_+ \cup I_-\f$ */
2865 	   SCIP_Real*            vzlp,               /**< array containing \f$v_i^T zlp_q\f$ for \f$i \in I_+ \cup I_-\f$ */
2866 	   SCIP_Real*            wcoefs,             /**< coefficients of w for the quad vars or NULL if w is 0 */
2867 	   SCIP_Real             wzlp,               /**< value of w at zlp */
2868 	   SCIP_Real             kappa,              /**< value of kappa */
2869 	   SCIP_Real*            interpoints,        /**< array to store intersection points for all rays or NULL if nothing
2870 	                                              *   needs to be stored */
2871 	   SCIP_Real*            rho,                /**< pointer to store the optimal rho */
2872 	   SCIP_Bool*            success             /**< could we successfully find the right rho? */
2873 	   )
2874 	{
2875 	   int i;
2876 	
2877 	   /* go through all rays not in the recession cone and compute the largest negative steplength possible. The
2878 	    * smallest of them is then the steplength rho we use for the current ray */
2879 	   *rho = 0.0;
2880 	   for( i = 0; i < rays->nrays; ++i )
2881 	   {
2882 	      SCIP_Real currentrho;
2883 	      SCIP_Real coef;
2884 	
2885 	      if( SCIPisInfinity(scip, interpoints[i]) )
2886 	         continue;
2887 	
2888 	      /* if we cannot strengthen enough, we don't strengthen at all */
2889 	      if( SCIPisInfinity(scip, -*rho) )
2890 	      {
2891 	         *rho = -SCIPinfinity(scip);
2892 	         return SCIP_OKAY;
2893 	      }
2894 	
2895 	      /* if the rays are linearly independent, we don't need to search for rho */
2896 	      if( raysAreDependent(scip, &rays->rays[rays->raysbegin[i]], &rays->raysidx[rays->raysbegin[i]],
2897 	          rays->raysbegin[i + 1] - rays->raysbegin[i], &rays->rays[rays->raysbegin[idx]],
2898 	          &rays->raysidx[rays->raysbegin[idx]], rays->raysbegin[idx + 1] - rays->raysbegin[idx], &coef) )
2899 	      {
2900 	         currentrho = coef * interpoints[i];
2901 	      }
2902 	      else
2903 	      {
2904 	         /*  since the two rays are linearly independent, we need to find the biggest alpha such that
2905 	          *  alpha * ray_i + (1 - alpha) * ray_idx in the recession cone is. For every ray i, we compute
2906 	          *  such a alpha but take the smallest one of them. We use "maxalpha" to keep track of this.
2907 	          *  Since we know that we can only use alpha < maxalpha, we don't need to do the whole binary search
2908 	          *  for every ray i. We only need to search the intervall [0, maxalpha]. Thereby, we start by checking
2909 	          *  if alpha = maxalpha is already feasable */
2910 	
2911 	         SCIP_Bool inreccone;
2912 	         SCIP_Real alpha;
2913 	         SCIP_Real lb;
2914 	         SCIP_Real ub;
2915 	         int j;
2916 	
2917 	         lb = 0.0;
2918 	         ub = 1.0;
2919 	         for( j = 0; j < BINSEARCH_MAXITERS; ++j )
2920 	         {
2921 	            alpha = (lb + ub) / 2.0;
2922 	
2923 	            if( SCIPisZero(scip, alpha) )
2924 	            {
2925 	               alpha = 0.0;
2926 	               break;
2927 	            }
2928 	
2929 	            SCIP_CALL( rayInRecessionCone(scip, nlhdlrdata, nlhdlrexprdata, rays, idx, i, sidefactor, iscase4, vb,
2930 	                  vzlp, wcoefs, wzlp, kappa, alpha, &inreccone, success) );
2931 	
2932 	            if( ! *success )
2933 	               return SCIP_OKAY;
2934 	
2935 	            /* no root exists */
2936 	            if( inreccone )
2937 	            {
2938 	               lb = alpha;
2939 	               if( SCIPisEQ(scip, ub, lb) )
2940 	                  break;
2941 	            }
2942 	            else
2943 	               ub = alpha;
2944 	         }
2945 	
2946 	         /* now we found the best convex combination which we use to derive the corresponding coef. If alpha = 0, we
2947 	          * cannot move the ray in the recession cone, i.e. strengthening is not possible */
2948 	         if( SCIPisZero(scip, alpha) )
2949 	         {
2950 	            *rho = -SCIPinfinity(scip);
2951 	            return SCIP_OKAY;
2952 	         }
2953 	         else
2954 	            currentrho = (alpha - 1) * interpoints[i] / alpha;  /*lint !e795*/
2955 	      }
2956 	
2957 	      if( currentrho < *rho )
2958 	         *rho = currentrho;
2959 	
2960 	      if( *rho < -10e+06 )
2961 	         *rho = -SCIPinfinity(scip);
2962 	
2963 	      /* if rho is too small, don't add it */
2964 	      if( SCIPisZero(scip, *rho) )
2965 	         *success = FALSE;
2966 	   }
2967 	
2968 	   return SCIP_OKAY;
2969 	}
2970 	
2971 	/** computes intersection cut using negative edge extension to strengthen rays that do not intersect
2972 	 * (i.e., rays in the recession cone)
2973 	 */
2974 	static
2975 	SCIP_RETCODE computeStrengthenedIntercut(
2976 	   SCIP*                 scip,               /**< SCIP data structure */
2977 	   SCIP_NLHDLRDATA*      nlhdlrdata,         /**< nlhdlr data */
2978 	   SCIP_NLHDLREXPRDATA*  nlhdlrexprdata,     /**< nlhdlr expression data */
2979 	   RAYS*                 rays,               /**< rays */
2980 	   SCIP_Real             sidefactor,         /**< 1.0 if the violated constraint is q &le; rhs, -1.0 otherwise */
2981 	   SCIP_Bool             iscase4,            /**< whether we are in case 4 */
2982 	   SCIP_Real*            vb,                 /**< array containing \f$v_i^T b\f$ for \f$i \in I_+ \cup I_-\f$ */
2983 	   SCIP_Real*            vzlp,               /**< array containing \f$v_i^T zlp_q\f$ for \f$i \in I_+ \cup I_-\f$ */
2984 	   SCIP_Real*            wcoefs,             /**< coefficients of w for the qud vars or NULL if w is 0 */
2985 	   SCIP_Real             wzlp,               /**< value of w at zlp */
2986 	   SCIP_Real             kappa,              /**< value of kappa */
2987 	   SCIP_ROWPREP*         rowprep,            /**< rowprep for the generated cut */
2988 	   SCIP_SOL*             sol,                /**< solution to separate */
2989 	   SCIP_Bool*            success,            /**< if a cut candidate could be computed */
2990 	   SCIP_Bool*            strengthsuccess     /**< if strengthening was successfully applied */
2991 	   )
2992 	{
2993 	   SCIP_COL** cols;
2994 	   SCIP_ROW** rows;
2995 	   SCIP_Real* interpoints;
2996 	   SCIP_Real avecutcoef;
2997 	   int counter;
2998 	   int i;
2999 	
3000 	   *success = TRUE;
3001 	   *strengthsuccess = FALSE;
3002 	
3003 	   cols = SCIPgetLPCols(scip);
3004 	   rows = SCIPgetLPRows(scip);
3005 	
3006 	   /* allocate memory for intersection points */
3007 	   SCIP_CALL( SCIPallocBufferArray(scip, &interpoints, rays->nrays) );
3008 	
3009 	   /* compute all intersection points and store them in interpoints; build not-stregthened intersection cut */
3010 	   SCIP_CALL( computeIntercut(scip, nlhdlrdata, nlhdlrexprdata, rays, sidefactor, iscase4, vb, vzlp, wcoefs, wzlp, kappa,
3011 	            rowprep, interpoints, sol, success) );
3012 	
3013 	   if( ! *success )
3014 	      goto CLEANUP;
3015 	
3016 	   /* keep track of the number of attempted strengthenings and average cutcoef */
3017 	   counter = 0;
3018 	   avecutcoef = 0.0;
3019 	
3020 	   /* go through all intersection points that are equal to infinity -> these correspond to the rays which are in the
3021 	    * recession cone of C, i.e. the rays for which we (possibly) can compute a negative steplength */
3022 	   for( i = 0; i < rays->nrays; ++i )
3023 	   {
3024 	      SCIP_Real rho;
3025 	      SCIP_Real cutcoef;
3026 	      int lppos;
3027 	
3028 	      if( !SCIPisInfinity(scip, interpoints[i]) )
3029 	         continue;
3030 	
3031 	      /* if we reached the limit of strengthenings, we stop */
3032 	      if( counter >= nlhdlrdata->nstrengthlimit )
3033 	         break;
3034 	
3035 	      /* compute the smallest rho */
3036 	      SCIP_CALL( findRho(scip, nlhdlrdata, nlhdlrexprdata, rays, i, sidefactor, iscase4, vb, vzlp, wcoefs, wzlp, kappa,
3037 	               interpoints, &rho, success));
3038 	
3039 	      /* compute cut coef */
3040 	      if( ! *success  || SCIPisInfinity(scip, -rho) )
3041 	         cutcoef = 0.0;
3042 	      else
3043 	         cutcoef = 1.0 / rho;
3044 	
3045 	      /* track average cut coef */
3046 	      counter += 1;
3047 	      avecutcoef += cutcoef;
3048 	
3049 	      if( ! SCIPisZero(scip, cutcoef) )
3050 	         *strengthsuccess = TRUE;
3051 	
3052 	      /* add var to cut: if variable is nonbasic at upper we have to flip sign of cutcoef */
3053 	      lppos = rays->lpposray[i];
3054 	      if( lppos < 0 )
3055 	      {
3056 	         lppos = -lppos - 1;
3057 	
3058 	         assert(SCIProwGetBasisStatus(rows[lppos]) == SCIP_BASESTAT_LOWER || SCIProwGetBasisStatus(rows[lppos]) ==
3059 	               SCIP_BASESTAT_UPPER);
3060 	
3061 	         SCIP_CALL( addRowToCut(scip, rowprep, SCIProwGetBasisStatus(rows[lppos]) == SCIP_BASESTAT_UPPER ? cutcoef :
3062 	                  -cutcoef, rows[lppos], success) ); /* rows have flipper base status! */
3063 	
3064 	         if( ! *success )
3065 	         {
3066 	            INTERLOG(printf("Bad numeric: row not nonbasic enough\n");)
3067 	            nlhdlrdata->nbadnonbasic++;
3068 	            return SCIP_OKAY;
3069 	         }
3070 	      }
3071 	      else
3072 	      {
3073 	         assert(SCIPcolGetBasisStatus(cols[lppos]) == SCIP_BASESTAT_UPPER || SCIPcolGetBasisStatus(cols[lppos]) ==
3074 	               SCIP_BASESTAT_LOWER);
3075 	         SCIP_CALL( addColToCut(scip, rowprep, sol, SCIPcolGetBasisStatus(cols[lppos]) == SCIP_BASESTAT_UPPER ? -cutcoef :
3076 	                  cutcoef, cols[lppos]) );
3077 	      }
3078 	   }
3079 	
3080 	   if( counter > 0 )
3081 	      nlhdlrdata->cutcoefsum += avecutcoef / counter;
3082 	
3083 	CLEANUP:
3084 	   SCIPfreeBufferArray(scip, &interpoints);
3085 	
3086 	   return SCIP_OKAY;
3087 	}
3088 	
3089 	/** sets variable in solution "vertex" to its nearest bound */
3090 	static
3091 	SCIP_RETCODE setVarToNearestBound(
3092 	   SCIP*                 scip,               /**< SCIP data structure */
3093 	   SCIP_SOL*             sol,                /**< solution to separate */
3094 	   SCIP_SOL*             vertex,             /**< new solution to separate */
3095 	   SCIP_VAR*             var,                /**< var we want to find nearest bound to */
3096 	   SCIP_Real*            factor,             /**< is vertex for current var at lower or upper? */
3097 	   SCIP_Bool*            success             /**< TRUE if no variable is bounded */
3098 	   )
3099 	{
3100 	   SCIP_Real solval;
3101 	   SCIP_Real bound;
3102 	
3103 	   solval = SCIPgetSolVal(scip, sol, var);
3104 	   *success = TRUE;
3105 	
3106 	   /* find nearest bound */
3107 	   if( SCIPisInfinity(scip, SCIPvarGetLbLocal(var)) && SCIPisInfinity(scip, SCIPvarGetUbLocal(var)) )
3108 	   {
3109 	      *success = FALSE;
3110 	      return SCIP_OKAY;
3111 	   }
3112 	   else if( solval - SCIPvarGetLbLocal(var) < SCIPvarGetUbLocal(var) - solval )
3113 	   {
3114 	      bound = SCIPvarGetLbLocal(var);
3115 	      *factor = 1.0;
3116 	   }
3117 	   else
3118 	   {
3119 	      bound = SCIPvarGetUbLocal(var);
3120 	      *factor = -1.0;
3121 	   }
3122 	
3123 	   /* set val to bound in solution */
3124 	   SCIP_CALL( SCIPsetSolVal(scip, vertex, var, bound) );
3125 	   return SCIP_OKAY;
3126 	}
3127 	
3128 	/** This function finds vertex (w.r.t. bounds of variables appearing in the quadratic) that is closest to the current
3129 	 * solution we want to separate.
3130 	 *
3131 	 * Furthermore, we store the rays corresponding to the unit vectors, i.e.,
3132 	 *    - if \f$x_i\f$ is at its lower bound in vertex --> \f$r_i =  e_i\f$
3133 	 *    - if \f$x_i\f$ is at its upper bound in vertex --> \f$r_i = -e_i\f$
3134 	 */
3135 	static
3136 	SCIP_RETCODE findVertexAndGetRays(
3137 	   SCIP*                 scip,               /**< SCIP data structure */
3138 	   SCIP_NLHDLREXPRDATA*  nlhdlrexprdata,     /**< nlhdlr expression data */
3139 	   SCIP_SOL*             sol,                /**< solution to separate */
3140 	   SCIP_SOL*             vertex,             /**< new 'vertex' (w.r.t. bounds) solution to separate */
3141 	   SCIP_VAR*             auxvar,             /**< aux var of expr or NULL if not needed (e.g. separating real cons) */
3142 	   RAYS**                raysptr,            /**< pointer to new bound rays */
3143 	   SCIP_Bool*            success             /**< TRUE if no variable is unbounded */
3144 	   )
3145 	{
3146 	   SCIP_EXPR* qexpr;
3147 	   SCIP_EXPR** linexprs;
3148 	   RAYS* rays;
3149 	   int nquadexprs;
3150 	   int nlinexprs;
3151 	   int raylength;
3152 	   int i;
3153 	
3154 	   *success = TRUE;
3155 	
3156 	   qexpr = nlhdlrexprdata->qexpr;
3157 	   SCIPexprGetQuadraticData(qexpr, NULL, &nlinexprs, &linexprs, NULL, &nquadexprs, NULL, NULL, NULL);
3158 	
3159 	   raylength = (auxvar == NULL) ? nquadexprs + nlinexprs : nquadexprs + nlinexprs + 1;
3160 	
3161 	   /* create rays */
3162 	   SCIP_CALL( createBoundRays(scip, raysptr, raylength) );
3163 	   rays = *raysptr;
3164 	
3165 	   rays->rayssize = raylength;
3166 	   rays->nrays = raylength;
3167 	
3168 	   /* go through quadratic variables */
3169 	   for( i = 0; i < nquadexprs; ++i )
3170 	   {
3171 	      SCIP_EXPR* expr;
3172 	      SCIPexprGetQuadraticQuadTerm(qexpr, i, &expr, NULL, NULL, NULL, NULL, NULL);
3173 	
3174 	      rays->raysbegin[i] = i;
3175 	      rays->raysidx[i] = i;
3176 	      rays->lpposray[i] = SCIPcolGetLPPos(SCIPvarGetCol(SCIPgetExprAuxVarNonlinear(expr)));
3177 	
3178 	      SCIP_CALL( setVarToNearestBound(scip, sol, vertex, SCIPgetExprAuxVarNonlinear(expr),
3179 	         &rays->rays[i], success) );
3180 	
3181 	      if( ! *success )
3182 	         return SCIP_OKAY;
3183 	   }
3184 	
3185 	   /* go through linear variables */
3186 	   for( i = 0; i < nlinexprs; ++i )
3187 	   {
3188 	      rays->raysbegin[i + nquadexprs] = i + nquadexprs;
3189 	      rays->raysidx[i + nquadexprs] = i + nquadexprs;
3190 	      rays->lpposray[i + nquadexprs] = SCIPcolGetLPPos(SCIPvarGetCol(SCIPgetExprAuxVarNonlinear(linexprs[i])));
3191 	
3192 	      SCIP_CALL( setVarToNearestBound(scip, sol, vertex, SCIPgetExprAuxVarNonlinear(linexprs[i]),
3193 	         &rays->rays[i + nquadexprs], success) );
3194 	
3195 	      if( ! *success )
3196 	         return SCIP_OKAY;
3197 	   }
3198 	
3199 	   /* consider auxvar if it exists */
3200 	   if( auxvar != NULL )
3201 	   {
3202 	      rays->raysbegin[nquadexprs + nlinexprs] = nquadexprs + nlinexprs;
3203 	      rays->raysidx[nquadexprs + nlinexprs] = nquadexprs + nlinexprs;
3204 	      rays->lpposray[nquadexprs + nlinexprs] = SCIPcolGetLPPos(SCIPvarGetCol(auxvar));
3205 	
3206 	      SCIP_CALL( setVarToNearestBound(scip, sol, vertex, auxvar, &rays->rays[nquadexprs + nlinexprs], success) );
3207 	
3208 	      if( ! *success )
3209 	         return SCIP_OKAY;
3210 	   }
3211 	
3212 	   rays->raysbegin[raylength] = raylength;
3213 	
3214 	   return SCIP_OKAY;
3215 	}
3216 	
3217 	/** checks if the quadratic constraint is violated by sol */
3218 	static
3219 	SCIP_Bool isQuadConsViolated(
3220 	   SCIP*                 scip,               /**< SCIP data structure */
3221 	   SCIP_NLHDLREXPRDATA*  nlhdlrexprdata,     /**< nlhdlr expression data */
3222 	   SCIP_VAR*             auxvar,             /**< aux var of expr or NULL if not needed (e.g. separating real cons) */
3223 	   SCIP_SOL*             sol,                /**< solution to check feasibility for */
3224 	   SCIP_Real             sidefactor          /**< 1.0 if the violated constraint is q &le; rhs, -1.0 otherwise */
3225 	   )
3226 	{
3227 	   SCIP_EXPR* qexpr;
3228 	   SCIP_EXPR** linexprs;
3229 	   SCIP_Real* lincoefs;
3230 	   SCIP_Real constant;
3231 	   SCIP_Real val;
3232 	   int nquadexprs;
3233 	   int nlinexprs;
3234 	   int nbilinexprs;
3235 	   int i;
3236 	
3237 	   qexpr = nlhdlrexprdata->qexpr;
3238 	   SCIPexprGetQuadraticData(qexpr, &constant, &nlinexprs, &linexprs, &lincoefs, &nquadexprs,
3239 	      &nbilinexprs, NULL, NULL);
3240 	
3241 	   val = 0.0;
3242 	
3243 	   /* go through quadratic terms */
3244 	   for( i = 0; i < nquadexprs; i++ )
3245 	   {
3246 	      SCIP_EXPR* expr;
3247 	      SCIP_Real quadlincoef;
3248 	      SCIP_Real sqrcoef;
3249 	      SCIP_Real solval;
3250 	
3251 	      SCIPexprGetQuadraticQuadTerm(qexpr, i, &expr, &quadlincoef, &sqrcoef, NULL, NULL, NULL);
3252 	
3253 	      solval = SCIPgetSolVal(scip, sol, SCIPgetExprAuxVarNonlinear(expr));
3254 	
3255 	      /* add square term */
3256 	      val += sqrcoef * SQR(solval);
3257 	
3258 	      /* add linear term */
3259 	      val += quadlincoef * solval;
3260 	   }
3261 	
3262 	   /* go through bilinear terms */
3263 	   for( i = 0; i < nbilinexprs; i++ )
3264 	   {
3265 	      SCIP_EXPR* expr1;
3266 	      SCIP_EXPR* expr2;
3267 	      SCIP_Real bilincoef;
3268 	
3269 	      SCIPexprGetQuadraticBilinTerm(qexpr, i, &expr1, &expr2, &bilincoef, NULL, NULL);
3270 	
3271 	      val += bilincoef * SCIPgetSolVal(scip, sol, SCIPgetExprAuxVarNonlinear(expr1))
3272 	         * SCIPgetSolVal(scip, sol, SCIPgetExprAuxVarNonlinear(expr2));
3273 	   }
3274 	
3275 	   /* go through linear terms */
3276 	   for( i = 0; i < nlinexprs; i++ )
3277 	   {
3278 	      val += lincoefs[i] * SCIPgetSolVal(scip, sol, SCIPgetExprAuxVarNonlinear(linexprs[i]));
3279 	   }
3280 	
3281 	   /* add auxvar if exists and get constant */
3282 	   if( auxvar != NULL )
3283 	   {
3284 	      val -= SCIPgetSolVal(scip, sol, auxvar);
3285 	
3286 	      constant = (sidefactor == 1.0) ? constant - SCIPgetRhsNonlinear(nlhdlrexprdata->cons) :
3287 	         SCIPgetLhsNonlinear(nlhdlrexprdata->cons) - constant;
3288 	   }
3289 	   else
3290 	      constant = (sidefactor * constant);
3291 	
3292 	   val = (sidefactor * val);
3293 	
3294 	   /* now constraint is q(z) <= const */
3295 	   if( val <= constant )
3296 	      return FALSE;
3297 	   else
3298 	      return TRUE;
3299 	}
3300 	
3301 	/** generates intersection cut that cuts off sol (which violates the quadratic constraint cons) */
3302 	static
3303 	SCIP_RETCODE generateIntercut(
3304 	   SCIP*                 scip,               /**< SCIP data structure */
3305 	   SCIP_EXPR*            expr,               /**< expr */
3306 	   SCIP_NLHDLRDATA*      nlhdlrdata,         /**< nlhdlr data */
3307 	   SCIP_NLHDLREXPRDATA*  nlhdlrexprdata,     /**< nlhdlr expression data */
3308 	   SCIP_CONS*            cons,               /**< violated constraint that contains expr */
3309 	   SCIP_SOL*             sol,                /**< solution to separate */
3310 	   SCIP_ROWPREP*         rowprep,            /**< rowprep for the generated cut */
3311 	   SCIP_Bool             overestimate,       /**< TRUE if viol cons is q(z) &ge; lhs; FALSE if q(z) &le; rhs */
3312 	   SCIP_Bool*            success             /**< whether separation was successfull or not */
3313 	   )
3314 	{
3315 	   SCIP_EXPR* qexpr;
3316 	   RAYS* rays;
3317 	   SCIP_VAR* auxvar;
3318 	   SCIP_Real sidefactor;
3319 	   SCIP_Real* vb;      /* eigenvectors * b */
3320 	   SCIP_Real* vzlp;    /* eigenvectors * lpsol */
3321 	   SCIP_Real* wcoefs;  /* coefficients affecting quadterms in w */
3322 	   SCIP_Real wzlp;     /* w(lpsol) */
3323 	   SCIP_Real kappa;
3324 	   SCIP_Bool iscase4;
3325 	   SCIP_SOL* vertex;
3326 	   SCIP_SOL* soltoseparate;
3327 	   int nquadexprs;
3328 	   int nlinexprs;
3329 	   int i;
3330 	
3331 	   /* count number of calls */
3332 	   (nlhdlrdata->ncalls++);
3333 	
3334 	   qexpr = nlhdlrexprdata->qexpr;
3335 	   SCIPexprGetQuadraticData(qexpr, NULL, &nlinexprs, NULL, NULL, &nquadexprs, NULL, NULL, NULL);
3336 	
3337 	#ifdef  DEBUG_INTERSECTIONCUT
3338 	   SCIPinfoMessage(scip, NULL, "Generating intersection cut for quadratic expr %p aka", (void*)expr);
3339 	   SCIP_CALL( SCIPprintExpr(scip, expr, NULL) );
3340 	   SCIPinfoMessage(scip, NULL, "\n");
3341 	#endif
3342 	
3343 	   *success = TRUE;
3344 	   iscase4 = TRUE;
3345 	
3346 	   /* in nonbasic space cut is >= 1 */
3347 	   assert(SCIProwprepGetSide(rowprep) == 0.0);
3348 	   SCIProwprepAddSide(rowprep, 1.0);
3349 	   SCIProwprepSetSidetype(rowprep, SCIP_SIDETYPE_LEFT);
3350 	   assert(SCIProwprepGetSide(rowprep) == 1.0);
3351 	
3352 	   auxvar = (nlhdlrexprdata->cons != cons) ? SCIPgetExprAuxVarNonlinear(expr) : NULL;
3353 	   sidefactor = overestimate ? -1.0 : 1.0;
3354 	
3355 	   rays = NULL;
3356 	
3357 	   /* check if we use tableau or bounds as rays */
3358 	   if( ! nlhdlrdata->useboundsasrays )
3359 	   {
3360 	      SCIP_CALL( createAndStoreSparseRays(scip, nlhdlrexprdata, auxvar, &rays, success) );
3361 	
3362 	      if( ! *success )
3363 	      {
3364 	         INTERLOG(printf("Failed to get rays: there is a var with base status ZERO!\n"); )
3365 	         return SCIP_OKAY;
3366 	      }
3367 	
3368 	      soltoseparate = sol;
3369 	   }
3370 	   else
3371 	   {
3372 	      SCIP_Bool violated;
3373 	
3374 	      if( auxvar != NULL )
3375 	      {
3376 	         *success = FALSE;
3377 	         return SCIP_OKAY;
3378 	      }
3379 	
3380 	      /* create new solution */
3381 	      SCIP_CALL( SCIPcreateSol(scip, &vertex, NULL) );
3382 	
3383 	      /* find nearest vertex of the box to separate and compute rays */
3384 	      SCIP_CALL( findVertexAndGetRays(scip, nlhdlrexprdata, sol, vertex, auxvar, &rays, success) );
3385 	
3386 	      if( ! *success )
3387 	      {
3388 	         INTERLOG(printf("Failed to use bounds as rays: variable is unbounded!\n"); )
3389 	         freeRays(scip, &rays);
3390 	         SCIP_CALL( SCIPfreeSol(scip, &vertex) );
3391 	         return SCIP_OKAY;
3392 	      }
3393 	
3394 	      /* check if vertex is violated */
3395 	      violated = isQuadConsViolated(scip, nlhdlrexprdata, auxvar, vertex, sidefactor);
3396 	
3397 	      if( ! violated )
3398 	      {
3399 	         INTERLOG(printf("Failed to use bounds as rays: nearest vertex is not violated!\n"); )
3400 	         freeRays(scip, &rays);
3401 	         SCIP_CALL( SCIPfreeSol(scip, &vertex) );
3402 	         *success = FALSE;
3403 	         return SCIP_OKAY;
3404 	      }
3405 	
3406 	      soltoseparate = vertex;
3407 	   }
3408 	
3409 	   SCIP_CALL( SCIPallocBufferArray(scip, &vb, nquadexprs) );
3410 	   SCIP_CALL( SCIPallocBufferArray(scip, &vzlp, nquadexprs) );
3411 	   SCIP_CALL( SCIPallocBufferArray(scip, &wcoefs, nquadexprs) );
3412 	
3413 	   SCIP_CALL( intercutsComputeCommonQuantities(scip, nlhdlrexprdata, auxvar, sidefactor, soltoseparate,
3414 	         vb, vzlp, wcoefs, &wzlp, &kappa) );
3415 	
3416 	   /* check if we are in case 4 */
3417 	   if( nlinexprs == 0 && auxvar == NULL )
3418 	   {
3419 	      for( i = 0; i < nquadexprs; ++i )
3420 	         if( wcoefs[i] != 0.0 )
3421 	            break;
3422 	
3423 	      if( i == nquadexprs )
3424 	      {
3425 	         /* from now on wcoefs is going to be NULL --> case 1, 2 or 3 */
3426 	         SCIPfreeBufferArray(scip, &wcoefs);
3427 	         iscase4 = FALSE;
3428 	      }
3429 	   }
3430 	
3431 	   /* compute (strengthened) intersection cut */
3432 	   if( nlhdlrdata->usestrengthening )
3433 	   {
3434 	      SCIP_Bool strengthsuccess;
3435 	
3436 	      SCIP_CALL( computeStrengthenedIntercut(scip, nlhdlrdata, nlhdlrexprdata, rays, sidefactor, iscase4, vb, vzlp, wcoefs,
3437 	            wzlp, kappa, rowprep, soltoseparate, success, &strengthsuccess) );
3438 	
3439 	      if( *success && strengthsuccess )
3440 	         nlhdlrdata->nstrengthenings++;
3441 	   }
3442 	   else
3443 	   {
3444 	      SCIP_CALL( computeIntercut(scip, nlhdlrdata, nlhdlrexprdata, rays, sidefactor, iscase4, vb, vzlp, wcoefs, wzlp, kappa,
3445 	            rowprep, NULL, soltoseparate, success) );
3446 	   }
3447 	
3448 	   SCIPfreeBufferArrayNull(scip, &wcoefs);
3449 	   SCIPfreeBufferArray(scip, &vzlp);
3450 	   SCIPfreeBufferArray(scip, &vb);
3451 	   freeRays(scip, &rays);
3452 	
3453 	   if( nlhdlrdata->useboundsasrays )
3454 	   {
3455 	      SCIP_CALL( SCIPfreeSol(scip, &vertex) );
3456 	   }
3457 	
3458 	   return SCIP_OKAY;
3459 	}
3460 	
3461 	/** returns whether a quadratic form is "propagable"
3462 	 *
3463 	 * It is propagable, if a variable (aka child expr) appears at least twice, which is the case if at least two of the following hold:
3464 	 * - it appears as a linear term (coef*expr)
3465 	 * - it appears as a square term (coef*expr^2)
3466 	 * - it appears in a bilinear term
3467 	 * - it appears in another bilinear term
3468 	 */
3469 	static
3470 	SCIP_Bool isPropagable(
3471 	   SCIP_EXPR*            qexpr               /**< quadratic representation data */
3472 	   )
3473 	{
3474 	   int nquadexprs;
3475 	   int i;
3476 	
3477 	   SCIPexprGetQuadraticData(qexpr, NULL, NULL, NULL, NULL, &nquadexprs, NULL, NULL, NULL);
3478 	
3479 	   for( i = 0; i < nquadexprs; ++i )
3480 	   {
3481 	      SCIP_Real lincoef;
3482 	      SCIP_Real sqrcoef;
3483 	      int nadjbilin;
3484 	
3485 	      SCIPexprGetQuadraticQuadTerm(qexpr, i, NULL, &lincoef, &sqrcoef, &nadjbilin, NULL, NULL);
3486 	
3487 	      if( (lincoef != 0.0) + (sqrcoef != 0.0) + nadjbilin >= 2 )  /*lint !e514*/ /* actually MIN(2, nadjbilin), but we check >= 2 */
3488 	         return TRUE;
3489 	   }
3490 	
3491 	   return FALSE;
3492 	}
3493 	
3494 	/** returns whether a quadratic term is "propagable"
3495 	 *
3496 	 * A term is propagable, if its variable (aka child expr) appears at least twice, which is the case if at least two of the following hold:
3497 	 * - it appears as a linear term (coef*expr)
3498 	 * - it appears as a square term (coef*expr^2)
3499 	 * - it appears in a bilinear term
3500 	 * - it appears in another bilinear term
3501 	 */
3502 	static
3503 	SCIP_Bool isPropagableTerm(
3504 	   SCIP_EXPR*            qexpr,              /**< quadratic representation data */
3505 	   int                   idx                 /**< index of quadratic term to consider */
3506 	   )
3507 	{
3508 	   SCIP_Real lincoef;
3509 	   SCIP_Real sqrcoef;
3510 	   int nadjbilin;
3511 	
3512 	   SCIPexprGetQuadraticQuadTerm(qexpr, idx, NULL, &lincoef, &sqrcoef, &nadjbilin, NULL,  NULL);
3513 	
3514 	   return (lincoef != 0.0) + (sqrcoef != 0.0) + nadjbilin >= 2;  /*lint !e514*/ /* actually MIN(2, nadjbilin), but we check >= 2 */
3515 	}
3516 	
3517 	/** solves a quadratic equation \f$ a\, \text{expr}^2 + b\, \text{expr} \in \text{rhs} \f$ (with \f$b\f$ an interval)
3518 	 * and reduces bounds on `expr` or deduces infeasibility if possible
3519 	 */
3520 	static
3521 	SCIP_RETCODE propagateBoundsQuadExpr(
3522 	   SCIP*                 scip,               /**< SCIP data structure */
3523 	   SCIP_EXPR*            expr,               /**< expression for which to solve */
3524 	   SCIP_Real             sqrcoef,            /**< square coefficient */
3525 	   SCIP_INTERVAL         b,                  /**< interval acting as linear coefficient */
3526 	   SCIP_INTERVAL         rhs,                /**< interval acting as rhs */
3527 	   SCIP_Bool*            infeasible,         /**< buffer to store if propagation produced infeasibility */
3528 	   int*                  nreductions         /**< buffer to store the number of interval reductions */
3529 	   )
3530 	{
3531 	   SCIP_INTERVAL a;
3532 	   SCIP_INTERVAL exprbounds;
3533 	   SCIP_INTERVAL newrange;
3534 	
3535 	   assert(scip != NULL);
3536 	   assert(expr != NULL);
3537 	   assert(infeasible != NULL);
3538 	   assert(nreductions != NULL);
3539 	
3540 	#ifdef DEBUG_PROP
3541 	   SCIPinfoMessage(scip, NULL, "Propagating <expr> by solving a <expr>^2 + b <expr> in rhs, where <expr> is: ");
3542 	   SCIP_CALL( SCIPprintExpr(scip, expr, NULL) );
3543 	   SCIPinfoMessage(scip, NULL, "\n");
3544 	   SCIPinfoMessage(scip, NULL, "expr in [%g, %g], a = %g, b = [%g, %g] and rhs = [%g, %g]\n",
3545 	         SCIPintervalGetInf(SCIPgetExprBoundsNonlinear(scip, expr)),
3546 	         SCIPintervalGetSup(SCIPgetExprBoundsNonlinear(scip, expr)), sqrcoef, b.inf, b.sup,
3547 	         rhs.inf, rhs.sup);
3548 	#endif
3549 	
3550 	   exprbounds = SCIPgetExprBoundsNonlinear(scip, expr);
3551 	   if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, exprbounds) )
3552 	   {
3553 	      *infeasible = TRUE;
3554 	      return SCIP_OKAY;
3555 	   }
3556 	
3557 	   /* compute solution of a*x^2 + b*x in rhs */
3558 	   SCIPintervalSet(&a, sqrcoef);
3559 	   SCIPintervalSolveUnivariateQuadExpression(SCIP_INTERVAL_INFINITY, &newrange, a, b, rhs, exprbounds);
3560 	
3561 	#ifdef DEBUG_PROP
3562 	   SCIPinfoMessage(scip, NULL, "Solution [%g, %g]\n", newrange.inf, newrange.sup);
3563 	#endif
3564 	
3565 	   SCIP_CALL( SCIPtightenExprIntervalNonlinear(scip, expr, newrange, infeasible, nreductions) );
3566 	
3567 	   return SCIP_OKAY;
3568 	}
3569 	
3570 	/** solves a linear equation \f$ b\, \text{expr} \in \text{rhs} \f$ (with \f$b\f$ a scalar) and reduces bounds on `expr` or deduces infeasibility if possible */
3571 	static
3572 	SCIP_RETCODE propagateBoundsLinExpr(
3573 	   SCIP*                 scip,               /**< SCIP data structure */
3574 	   SCIP_EXPR*            expr,               /**< expression for which to solve */
3575 	   SCIP_Real             b,                  /**< linear coefficient */
3576 	   SCIP_INTERVAL         rhs,                /**< interval acting as rhs */
3577 	   SCIP_Bool*            infeasible,         /**< buffer to store if propagation produced infeasibility */
3578 	   int*                  nreductions         /**< buffer to store the number of interval reductions */
3579 	   )
3580 	{
3581 	   SCIP_INTERVAL newrange;
3582 	
3583 	   assert(scip != NULL);
3584 	   assert(expr != NULL);
3585 	   assert(infeasible != NULL);
3586 	   assert(nreductions != NULL);
3587 	
3588 	#ifdef DEBUG_PROP
3589 	   SCIPinfoMessage(scip, NULL, "Propagating <expr> by solving %g <expr> in [%g, %g], where <expr> is: ", b, rhs.inf, rhs.sup);
3590 	   SCIP_CALL( SCIPprintExpr(scip, expr, NULL) );
3591 	   SCIPinfoMessage(scip, NULL, "\n");
3592 	#endif
3593 	
3594 	   /* compute solution of b*x in rhs */
3595 	   SCIPintervalDivScalar(SCIP_INTERVAL_INFINITY, &newrange, rhs, b);
3596 	
3597 	#ifdef DEBUG_PROP
3598 	   SCIPinfoMessage(scip, NULL, "Solution [%g, %g]\n", newrange.inf, newrange.sup);
3599 	#endif
3600 	
3601 	   SCIP_CALL( SCIPtightenExprIntervalNonlinear(scip, expr, newrange, infeasible, nreductions) );
3602 	
3603 	   return SCIP_OKAY;
3604 	}
3605 	
3606 	/** returns max of a/x - c*x for x in {x1, x2} with x1, x2 > 0 */
3607 	static
3608 	SCIP_Real computeMaxBoundaryForBilinearProp(
3609 	   SCIP_Real             a,                  /**< coefficient a */
3610 	   SCIP_Real             c,                  /**< coefficient c */
3611 	   SCIP_Real             x1,                 /**< coefficient x1 > 0 */
3612 	   SCIP_Real             x2                  /**< coefficient x2 > 0 */
3613 	   )
3614 	{
3615 	   SCIP_Real cneg;
3616 	   SCIP_Real cand1;
3617 	   SCIP_Real cand2;
3618 	   SCIP_ROUNDMODE roundmode;
3619 	
3620 	   assert(x1 > 0.0);
3621 	   assert(x2 > 0.0);
3622 	
3623 	   cneg = SCIPintervalNegateReal(c);
3624 	
3625 	   roundmode = SCIPintervalGetRoundingMode();
3626 	   SCIPintervalSetRoundingModeUpwards();
3627 	   cand1 = a/x1 + cneg*x1;
3628 	   cand2 = a/x2 + cneg*x2;
3629 	   SCIPintervalSetRoundingMode(roundmode);
3630 	
3631 	   return MAX(cand1, cand2);
3632 	}
3633 	
3634 	/** returns max of a/x - c*x for x in dom; it assumes that dom is contained in (0, +inf) */
3635 	static
3636 	SCIP_Real computeMaxForBilinearProp(
3637 	   SCIP_Real             a,                  /**< coefficient a */
3638 	   SCIP_Real             c,                  /**< coefficient c */
3639 	   SCIP_INTERVAL         dom                 /**< domain of x */
3640 	   )
3641 	{
3642 	   SCIP_ROUNDMODE roundmode;
3643 	   SCIP_INTERVAL argmax;
3644 	   SCIP_Real negunresmax;
3645 	   SCIP_Real boundarymax;
3646 	   assert(dom.inf > 0);
3647 	
3648 	   /* if a >= 0, then the function is convex which means the maximum is at one of the boundaries
3649 	    *
3650 	    * if c = 0, then the function is monotone which means the maximum is also at one of the boundaries
3651 	    *
3652 	    * if a < 0, then the function is concave. The function then has a maximum if and only if there is a point with derivative 0,
3653 	    * that is, iff -a/x^2 - c = 0 has a solution; i.e. if -a/c >= 0, i.e. (using a<0 and c != 0), c > 0.
3654 	    * Otherwise (that is, c<0), the maximum is at one of the boundaries.
3655 	    */
3656 	   if( a >= 0.0 || c <= 0.0 )
3657 	      return computeMaxBoundaryForBilinearProp(a, c, dom.inf, dom.sup);
3658 	
3659 	   /* now, the (unrestricted) maximum is at sqrt(-a/c).
3660 	    * if the argmax is not in the interior of dom then the solution is at a boundary, too
3661 	    * we check this by computing an interval that contains sqrt(-a/c) first
3662 	    */
3663 	   SCIPintervalSet(&argmax, -a);
3664 	   SCIPintervalDivScalar(SCIP_INTERVAL_INFINITY, &argmax, argmax, c);
3665 	   SCIPintervalSquareRoot(SCIP_INTERVAL_INFINITY, &argmax, argmax);
3666 	
3667 	   /* if the interval containing sqrt(-a/c) does not intersect with the interior of dom, then
3668 	    * the (restricted) maximum is at a boundary (we could even say at which boundary, but that doesn't save much)
3669 	    */
3670 	   if( argmax.sup <= dom.inf || argmax.inf >= dom.sup )
3671 	      return computeMaxBoundaryForBilinearProp(a, c, dom.inf, dom.sup);
3672 	
3673 	   /* the maximum at sqrt(-a/c) is -2*sqrt(-a*c), so we compute an upper bound for that by computing a lower bound for 2*sqrt(-a*c) */
3674 	   roundmode = SCIPintervalGetRoundingMode();
3675 	   SCIPintervalSetRoundingModeDownwards();
3676 	   negunresmax = 2.0*SCIPnextafter(sqrt(SCIPintervalNegateReal(a)*c), 0.0);
3677 	   SCIPintervalSetRoundingMode(roundmode);
3678 	
3679 	   /* if the interval containing sqrt(-a/c) is contained in dom, then we can return -negunresmax */
3680 	   if( argmax.inf >= dom.inf && argmax.sup <= dom.sup )
3681 	      return -negunresmax;
3682 	
3683 	   /* now what is left is the case where we cannot say for sure whether sqrt(-a/c) is contained in dom or not
3684 	    * so we are conservative and return the max of both cases, i.e.,
3685 	    * the max of the upper bounds on -2*sqrt(-a*c), a/dom.inf-c*dom.inf, a/dom.sup-c*dom.sup.
3686 	    */
3687 	   boundarymax = computeMaxBoundaryForBilinearProp(a, c, dom.inf, dom.sup);
3688 	   return MAX(boundarymax, -negunresmax);
3689 	}
3690 	
3691 	/** computes the range of rhs/x - coef * x for x in exprdom; this is used for the propagation of bilinear terms
3692 	 *
3693 	 * If 0 is in the exprdom, we set range to \f$\mathbb{R}\f$ (even though this is not quite correct, it is correct for the
3694 	 * intended use of the function).
3695 	 * TODO: maybe check before calling it whether 0 is in the domain and then just avoid calling it
3696 	 *
3697 	 * If rhs is [A,B] and x > 0, then we want the min of A/x - coef*x and max of B/x - coef*x for x in [exprdom].
3698 	 * If rhs is [A,B] and x < 0, then we want the min of B/x - coef*x and max of A/x - coef*x for x in [exprdom].
3699 	 * However, this is the same as min of -B/x + coef*x and max of -A/x + coef*x for x in -[exprdom].
3700 	 * Thus, we can always reduce to x > 0 by multiplying [exprdom], rhs, and coef by -1.
3701 	 */
3702 	static
3703 	void computeRangeForBilinearProp(
3704 	   SCIP_INTERVAL         exprdom,            /**< expression for which to solve */
3705 	   SCIP_Real             coef,               /**< expression for which to solve */
3706 	   SCIP_INTERVAL         rhs,                /**< rhs used for computation */
3707 	   SCIP_INTERVAL*        range               /**< storage for the resulting range */
3708 	   )
3709 	{
3710 	   SCIP_Real max;
3711 	   SCIP_Real min;
3712 	
3713 	   if( exprdom.inf <= 0.0 && 0.0 <= exprdom.sup )
3714 	   {
3715 	      SCIPintervalSetEntire(SCIP_INTERVAL_INFINITY, range);
3716 	      return;
3717 	   }
3718 	
3719 	   /* reduce to positive case */
3720 	   if( exprdom.sup < 0 )
3721 	   {
3722 	      SCIPintervalMulScalar(SCIP_INTERVAL_INFINITY, &exprdom, exprdom, -1.0);
3723 	      SCIPintervalMulScalar(SCIP_INTERVAL_INFINITY, &rhs, rhs, -1.0);
3724 	      coef *= -1.0;
3725 	   }
3726 	   assert(exprdom.inf > 0.0);
3727 	
3728 	   /* compute maximum and minimum */
3729 	   max = computeMaxForBilinearProp(rhs.sup, coef, exprdom);
3730 	   min = -computeMaxForBilinearProp(-rhs.inf, -coef, exprdom);
3731 	
3732 	   /* set interval */
3733 	   SCIPintervalSetBounds(range, min, max);
3734 	}
3735 	
3736 	/** reverse propagates coef_i expr_i + constant in rhs */
3737 	static
3738 	SCIP_RETCODE reversePropagateLinearExpr(
3739 	   SCIP*                 scip,               /**< SCIP data structure */
3740 	   SCIP_EXPR**           linexprs,           /**< linear expressions */
3741 	   int                   nlinexprs,          /**< number of linear expressions */
3742 	   SCIP_Real*            lincoefs,           /**< coefficients of linear expressions */
3743 	   SCIP_Real             constant,           /**< constant */
3744 	   SCIP_INTERVAL         rhs,                /**< rhs */
3745 	   SCIP_Bool*            infeasible,         /**< buffer to store whether an exps' bounds were propagated to an empty interval */
3746 	   int*                  nreductions         /**< buffer to store the number of interval reductions of all exprs */
3747 	   )
3748 	{
3749 	   SCIP_INTERVAL* oldboundslin;
3750 	   SCIP_INTERVAL* newboundslin;
3751 	   int i;
3752 	
3753 	   if( nlinexprs == 0 )
3754 	      return SCIP_OKAY;
3755 	
3756 	   SCIP_CALL( SCIPallocBufferArray(scip, &oldboundslin, nlinexprs) );
3757 	   SCIP_CALL( SCIPallocBufferArray(scip, &newboundslin, nlinexprs) );
3758 	
3759 	   for( i = 0; i < nlinexprs; ++i )
3760 	      oldboundslin[i] = SCIPexprGetActivity(linexprs[i]);  /* TODO use SCIPgetExprBoundsNonlinear(scip, linexprs[i]) ? */
3761 	
3762 	   *nreductions = SCIPintervalPropagateWeightedSum(SCIP_INTERVAL_INFINITY, nlinexprs,
3763 	            oldboundslin, lincoefs, constant, rhs, newboundslin, infeasible);
3764 	
3765 	   if( *nreductions > 0 && !*infeasible )
3766 	   {
3767 	      /* SCIP is more conservative with what constitutes a reduction than interval arithmetic so we follow SCIP */
3768 	      *nreductions = 0;
3769 	      for( i = 0; i < nlinexprs && ! (*infeasible); ++i )
3770 	      {
3771 	         SCIP_CALL( SCIPtightenExprIntervalNonlinear(scip, linexprs[i], newboundslin[i], infeasible, nreductions) );
3772 	      }
3773 	   }
3774 	
3775 	   SCIPfreeBufferArray(scip, &newboundslin);
3776 	   SCIPfreeBufferArray(scip, &oldboundslin);
3777 	
3778 	   return SCIP_OKAY;
3779 	}
3780 	
3781 	
3782 	/*
3783 	 * Callback methods of nonlinear handler
3784 	 */
3785 	
3786 	/** callback to free expression specific data */
3787 	static
3788 	SCIP_DECL_NLHDLRFREEEXPRDATA(nlhdlrFreeexprdataQuadratic)
3789 	{  /*lint --e{715}*/
3790 	   assert(nlhdlrexprdata != NULL);
3791 	   assert(*nlhdlrexprdata != NULL);
3792 	
3793 	   if( (*nlhdlrexprdata)->quadactivities != NULL )
3794 	   {
3795 	      int nquadexprs;
3796 	      SCIPexprGetQuadraticData((*nlhdlrexprdata)->qexpr, NULL, NULL, NULL, NULL, &nquadexprs, NULL, NULL, NULL);
3797 	      SCIPfreeBlockMemoryArray(scip, &(*nlhdlrexprdata)->quadactivities, nquadexprs);
3798 	   }
3799 	
3800 	   SCIPfreeBlockMemory(scip, nlhdlrexprdata);
3801 	
3802 	   return SCIP_OKAY;
3803 	}
3804 	
3805 	/** callback to detect structure in expression tree
3806 	 *
3807 	 * A term is quadratic if
3808 	 * - it is a product expression of two expressions, or
3809 	 * - it is power expression of an expression with exponent 2.0.
3810 	 *
3811 	 * We define a _propagable_ quadratic expression as a quadratic expression whose termwise propagation does not yield the
3812 	 * best propagation. In other words, is a quadratic expression that suffers from the dependency problem.
3813 	 *
3814 	 * Specifically, a propagable quadratic expression is a sum expression such that there is at least one expr that appears
3815 	 * at least twice (because of simplification, this means it appears in a quadratic terms and somewhere else).
3816 	 * For example: \f$x^2 + y^2\f$ is not a propagable quadratic expression; \f$x^2 + x\f$ is a propagable quadratic expression;
3817 	 * \f$x^2 + x y\f$ is also a propagable quadratic expression
3818 	 *
3819 	 * Furthermore, we distinguish between propagable and non-propagable terms. A term is propagable if any of the expressions
3820 	 * involved in it appear somewhere else. For example, \f$xy + z^2 + z\f$ is a propagable quadratic, the term \f$xy\f$ is
3821 	 * non-propagable, and \f$z^2\f$ is propagable. For propagation, non-propagable terms are handled as if they were linear
3822 	 * terms, that is, we do not use the activity of \f$x\f$ and \f$y\f$ to compute the activity of \f$xy\f$ but rather we use directly
3823 	 * the activity of \f$xy\f$. Similarly, we do not backward propagate to \f$x\f$ and \f$y\f$ (the product expr handler will do this),
3824 	 * but we backward propagate to \f$x*y\f$. More technically, we register \f$xy\f$ for its activity usage, rather than\f$x\f$ and \f$y\f$.
3825 	 *
3826 	 * For propagation, we store the quadratic in our data structure in the following way: We count how often a variable
3827 	 * appears. Then, a bilinear product expr_i * expr_j is stored as expr_i * expr_j if # expr_i appears > # expr_j
3828 	 * appears. When # expr_i appears = # expr_j appears, it then it will be stored as expr_i * expr_j if and only if
3829 	 * expr_i < expr_j, where '<' is the expression order (see \ref EXPR_ORDER "Ordering Rules" in \ref scip_expr.h).
3830 	 * Heuristically, this should be useful for propagation. The intuition is that by factoring out the variable that
3831 	 * appears most often we should be able to take care of the dependency problem better.
3832 	 *
3833 	 * Simple convex quadratics like \f$x^2 + y^2\f$ are ignored since the default nlhdlr will take care of them.
3834 	 *
3835 	 * @note The expression needs to be simplified (in particular, it is assumed to be sorted).
3836 	 * @note Common subexpressions are also assumed to have been identified, the hashing will fail otherwise!
3837 	 *
3838 	 * Sorted implies that:
3839 	 *  - expr < expr^2: bases are the same, but exponent 1 < 2
3840 	 *  - expr < expr * other_expr: u*v < w holds if and only if v < w (OR8), but here w = u < v, since expr comes before
3841 	 *  other_expr in the product
3842 	 *  - expr < other_expr * expr: u*v < w holds if and only if v < w (OR8), but here v = w
3843 	 *
3844 	 *  Thus, if we see somebody twice, it is a propagable quadratic.
3845 	 *
3846 	 * It also implies that
3847 	 *  - expr^2 < expr * other_expr
3848 	 *  - other_expr * expr < expr^2
3849 	 *
3850 	 * It also implies that x^-2 < x^-1, but since, so far, we do not interpret x^-2 as (x^-1)^2, it is not a problem.
3851 	 */
3852 	static
3853 	SCIP_DECL_NLHDLRDETECT(nlhdlrDetectQuadratic)
3854 	{  /*lint --e{715,774}*/
3855 	   SCIP_NLHDLREXPRDATA* nlexprdata;
3856 	   SCIP_NLHDLRDATA* nlhdlrdata;
3857 	   SCIP_Real* eigenvalues;
3858 	   SCIP_Bool isquadratic;
3859 	   SCIP_Bool propagable;
3860 	
3861 	   assert(scip != NULL);
3862 	   assert(nlhdlr != NULL);
3863 	   assert(expr != NULL);
3864 	   assert(enforcing != NULL);
3865 	   assert(participating != NULL);
3866 	   assert(nlhdlrexprdata != NULL);
3867 	
3868 	   nlhdlrdata = SCIPnlhdlrGetData(nlhdlr);
3869 	   assert(nlhdlrdata != NULL);
3870 	
3871 	   /* don't check if all enforcement methods are already ensured */
3872 	   if( (*enforcing & SCIP_NLHDLR_METHOD_ALL) == SCIP_NLHDLR_METHOD_ALL )
3873 	      return SCIP_OKAY;
3874 	
3875 	   /* if it is not a sum of at least two terms, it is not interesting */
3876 	   /* TODO: constraints of the form l<= x*y <= r ? */
3877 	   if( ! SCIPisExprSum(scip, expr) || SCIPexprGetNChildren(expr) < 2 )
3878 	      return SCIP_OKAY;
3879 	
3880 	   /* If we are in a subSCIP we don't want to separate intersection cuts */
3881 	   if( SCIPgetSubscipDepth(scip) > 0 )
3882 	      nlhdlrdata->useintersectioncuts = FALSE;
3883 	
3884 	#ifdef SCIP_DEBUG
3885 	   SCIPinfoMessage(scip, NULL, "Nlhdlr quadratic detecting expr %p aka ", (void*)expr);
3886 	   SCIP_CALL( SCIPprintExpr(scip, expr, NULL) );
3887 	   SCIPinfoMessage(scip, NULL, "\n");
3888 	   SCIPinfoMessage(scip, NULL, "Have to enforce %d\n", *enforcing);
3889 	#endif
3890 	
3891 	   /* check whether expression is quadratic (a sum with at least one square or bilinear term) */
3892 	   SCIP_CALL( SCIPcheckExprQuadratic(scip, expr, &isquadratic) );
3893 	
3894 	   /* not quadratic -> nothing for us */
3895 	   if( !isquadratic )
3896 	   {
3897 	      SCIPdebugMsg(scip, "expr %p is not quadratic -> abort detect\n", (void*)expr);
3898 	      return SCIP_OKAY;
3899 	   }
3900 	
3901 	   propagable = isPropagable(expr);
3902 	
3903 	   /* if we are not propagable and are in presolving, return */
3904 	   if( !propagable && SCIPgetStage(scip) == SCIP_STAGE_PRESOLVING )
3905 	   {
3906 	      SCIPdebugMsg(scip, "expr %p is not propagable and in presolving -> abort detect\n", (void*)expr);
3907 	      return SCIP_OKAY;
3908 	   }
3909 	
3910 	   /* if we do not use intersection cuts and are not propagable, then we do not want to handle it at all;
3911 	    * if not propagable, then we need to check the curvature to decide if we want to generate intersection cuts
3912 	    */
3913 	   if( !propagable && !nlhdlrdata->useintersectioncuts )
3914 	   {
3915 	      SCIPdebugMsg(scip, "expr %p is not propagable -> abort detect\n", (void*)expr);
3916 	      return SCIP_OKAY;
3917 	   }
3918 	
3919 	   /* store quadratic in nlhdlrexprdata */
3920 	   SCIP_CALL( SCIPallocClearBlockMemory(scip, nlhdlrexprdata) );
3921 	   nlexprdata = *nlhdlrexprdata;
3922 	   nlexprdata->qexpr = expr;
3923 	   nlexprdata->cons = cons;
3924 	
3925 	#ifdef DEBUG_DETECT
3926 	   SCIPinfoMessage(scip, NULL, "Nlhdlr quadratic detected:\n");
3927 	   SCIP_CALL( SCIPprintExprQuadratic(scip, conshdlr, qexpr) );
3928 	#endif
3929 	
3930 	   /* every propagable quadratic expression will be handled since we can propagate */
3931 	   if( propagable )
3932 	   {
3933 	      SCIP_EXPR** linexprs;
3934 	      int nlinexprs;
3935 	      int nquadexprs;
3936 	      int nbilin;
3937 	      int i;
3938 	
3939 	      *participating |= SCIP_NLHDLR_METHOD_ACTIVITY;
3940 	      *enforcing |= SCIP_NLHDLR_METHOD_ACTIVITY;
3941 	
3942 	      SCIPexprGetQuadraticData(expr, NULL, &nlinexprs, &linexprs, NULL, &nquadexprs, &nbilin, NULL, NULL);
3943 	      SCIP_CALL( SCIPallocBlockMemoryArray(scip, &nlexprdata->quadactivities, nquadexprs) );
3944 	
3945 	      /* notify children of quadratic that we will need their activity for propagation */
3946 	      for( i = 0; i < nlinexprs; ++i )
3947 	         SCIP_CALL( SCIPregisterExprUsageNonlinear(scip, linexprs[i], FALSE, TRUE, FALSE, FALSE) );
3948 	
3949 	      for( i = 0; i < nquadexprs; ++i )
3950 	      {
3951 	         SCIP_EXPR* argexpr;
3952 	         if( isPropagableTerm(expr, i) )
3953 	         {
3954 	            SCIPexprGetQuadraticQuadTerm(expr, i, &argexpr, NULL, NULL, &nbilin, NULL, NULL);
3955 	            SCIP_CALL( SCIPregisterExprUsageNonlinear(scip, argexpr, FALSE, TRUE, FALSE, FALSE) );
3956 	
3957 	#ifdef DEBUG_DETECT
3958 	            SCIPinfoMessage(scip, NULL, "quadterm %d propagable, using %p, unbounded=%d\n", i, (void*)argexpr, nbilin >
3959 	                  0 && SCIPintervalIsEntire(SCIP_INTERVAL_INFINITY, SCIPexprGetActivity(scip, argexpr)));
3960 	#endif
3961 	         }
3962 	         else
3963 	         {
3964 	            /* non-propagable quadratic is either a single square term or a single bilinear term
3965 	             * we should make use nlhdlrs in pow or product for this term, so we register usage of the square or product
3966 	             * expr instead of argexpr
3967 	             */
3968 	            SCIP_EXPR* sqrexpr;
3969 	            int* adjbilin;
3970 	
3971 	            SCIPexprGetQuadraticQuadTerm(expr, i, &argexpr, NULL, NULL, &nbilin, &adjbilin, &sqrexpr);
3972 	
3973 	            if( sqrexpr != NULL )
3974 	            {
3975 	               SCIP_CALL( SCIPregisterExprUsageNonlinear(scip, sqrexpr, FALSE, TRUE, FALSE, FALSE) );
3976 	               assert(nbilin == 0);
3977 	
3978 	#ifdef DEBUG_DETECT
3979 	               SCIPinfoMessage(scip, NULL, "quadterm %d non-propagable square, using %p\n", i, (void*)sqrexpr);
3980 	#endif
3981 	            }
3982 	            else
3983 	            {
3984 	               /* we have expr1 * other_expr or other_expr * expr1; know that expr1 is non propagable, but to decide if
3985 	                * we want the bounds of expr1 or of the product expr1 * other_expr (or other_expr * expr1), we have to
3986 	                * decide whether other_expr is also non propagable; due to the way we sort bilinear terms (by
3987 	                * frequency), we can deduce that other_expr doesn't appear anywhere else (i.e. is non propagable) if the
3988 	                * product is of the form expr1 * other_expr; however, if we see other_expr * expr1 we need to find
3989 	                * other_expr and check whether it is propagable
3990 	                */
3991 	               SCIP_EXPR* expr1;
3992 	               SCIP_EXPR* prodexpr;
3993 	
3994 	               assert(nbilin == 1);
3995 	               SCIPexprGetQuadraticBilinTerm(expr, adjbilin[0], &expr1, NULL, NULL, NULL, &prodexpr);
3996 	
3997 	               if( expr1 == argexpr )
3998 	               {
3999 	                  SCIP_CALL( SCIPregisterExprUsageNonlinear(scip, prodexpr, FALSE, TRUE, FALSE, FALSE) );
4000 	
4001 	#ifdef DEBUG_DETECT
4002 	                  SCIPinfoMessage(scip, NULL, "quadterm %d non-propagable product, using %p\n", i, (void*)prodexpr);
4003 	#endif
4004 	               }
4005 	               else
4006 	               {
4007 	                  int j;
4008 	                  /* check if other_expr is propagable in which case we need the bounds of expr1; otherwise we just need
4009 	                   * the bounds of the product and this will be (or was) registered when the loop takes us to the
4010 	                   * quadexpr other_expr.
4011 	                   * TODO this should be done faster, maybe store pos1 in bilinexprterm or store quadexprterm's in bilinexprterm
4012 	                   */
4013 	                  for( j = 0; j < nquadexprs; ++j )
4014 	                  {
4015 	                     SCIP_EXPR* exprj;
4016 	                     SCIPexprGetQuadraticQuadTerm(expr, j, &exprj, NULL, NULL, NULL, NULL, NULL);
4017 	                     if( expr1 == exprj )
4018 	                     {
4019 	                        if( isPropagableTerm(expr, j) )
4020 	                        {
4021 	                           SCIP_CALL( SCIPregisterExprUsageNonlinear(scip, argexpr, FALSE, TRUE, FALSE, FALSE) );
4022 	#ifdef DEBUG_DETECT
4023 	                           SCIPinfoMessage(scip, NULL, "quadterm %d non-propagable alien product, using %p\n", i, (void*)argexpr);
4024 	#endif
4025 	                        }
4026 	                        break;
4027 	                     }
4028 	                  }
4029 	               }
4030 	            }
4031 	         }
4032 	      }
4033 	   }
4034 	
4035 	   /* check if we are going to separate or not */
4036 	   nlexprdata->curvature = SCIP_EXPRCURV_UNKNOWN;
4037 	
4038 	   /* for now, we do not care about separation if it is not required */
4039 	   if( (*enforcing & SCIP_NLHDLR_METHOD_SEPABOTH) == SCIP_NLHDLR_METHOD_SEPABOTH )
4040 	   {
4041 	      /* if nobody can do anything, remove data */
4042 	      if( *participating == SCIP_NLHDLR_METHOD_NONE ) /*lint !e845*/
4043 	      {
4044 	         SCIP_CALL( nlhdlrFreeexprdataQuadratic(scip, nlhdlr, expr, nlhdlrexprdata) );
4045 	      }
4046 	      else
4047 	      {
4048 	         SCIPdebugMsg(scip, "expr %p is quadratic and propagable -> propagate\n", (void*)expr);
4049 	      }
4050 	      return SCIP_OKAY;
4051 	   }
4052 	
4053 	   assert(SCIPgetStage(scip) >= SCIP_STAGE_INITSOLVE);  /* separation should only be required in (init)solving stage */
4054 	
4055 	   /* check if we can do something more: check curvature of quadratic function stored in nlexprdata
4056 	    * this is currently only used to decide whether we want to separate, so it can be skipped if in presolve
4057 	    */
4058 	   SCIPdebugMsg(scip, "checking curvature of expr %p\n", (void*)expr);
4059 	   SCIP_CALL( SCIPcomputeExprQuadraticCurvature(scip, expr, &nlexprdata->curvature, NULL, nlhdlrdata->useintersectioncuts) );
4060 	
4061 	   /* get eigenvalues to be able to check whether they were computed */
4062 	   SCIPexprGetQuadraticData(expr, NULL, NULL, NULL, NULL, NULL, NULL, &eigenvalues, NULL);
4063 	
4064 	   /* if we use intersection cuts then we can handle any non-convex quadratic */
4065 	   if( nlhdlrdata->useintersectioncuts && eigenvalues != NULL && (*enforcing & SCIP_NLHDLR_METHOD_SEPABELOW) ==
4066 	         FALSE && nlexprdata->curvature != SCIP_EXPRCURV_CONVEX )
4067 	   {
4068 	      *participating |= SCIP_NLHDLR_METHOD_SEPABELOW;
4069 	   }
4070 	
4071 	   if( nlhdlrdata->useintersectioncuts && eigenvalues != NULL && (*enforcing & SCIP_NLHDLR_METHOD_SEPAABOVE) == FALSE &&
4072 	         nlexprdata->curvature != SCIP_EXPRCURV_CONCAVE )
4073 	   {
4074 	      *participating |= SCIP_NLHDLR_METHOD_SEPAABOVE;
4075 	   }
4076 	
4077 	   /* if nobody can do anything, remove data */
4078 	   if( *participating == SCIP_NLHDLR_METHOD_NONE ) /*lint !e845*/
4079 	   {
4080 	      SCIP_CALL( nlhdlrFreeexprdataQuadratic(scip, nlhdlr, expr, nlhdlrexprdata) );
4081 	      return SCIP_OKAY;
4082 	   }
4083 	
4084 	   /* we only need auxiliary variables if we are going to separate */
4085 	   if( *participating & SCIP_NLHDLR_METHOD_SEPABOTH )
4086 	   {
4087 	      SCIP_EXPR** linexprs;
4088 	      int nquadexprs;
4089 	      int nlinexprs;
4090 	      int i;
4091 	
4092 	      SCIPexprGetQuadraticData(expr, NULL, &nlinexprs, &linexprs, NULL, &nquadexprs, NULL, NULL, NULL);
4093 	
4094 	      for( i = 0; i < nlinexprs; ++i ) /* expressions appearing linearly */
4095 	      {
4096 	         SCIP_CALL( SCIPregisterExprUsageNonlinear(scip, linexprs[i], TRUE, FALSE, FALSE, FALSE) );
4097 	      }
4098 	      for( i = 0; i < nquadexprs; ++i ) /* expressions appearing quadratically */
4099 	      {
4100 	         SCIP_EXPR* quadexpr;
4101 	         SCIPexprGetQuadraticQuadTerm(expr, i, &quadexpr, NULL, NULL, NULL, NULL, NULL);
4102 	         SCIP_CALL( SCIPregisterExprUsageNonlinear(scip, quadexpr, TRUE, FALSE, FALSE, FALSE) );
4103 	      }
4104 	
4105 	      SCIPdebugMsg(scip, "expr %p is quadratic and propagable -> propagate and separate\n", (void*)expr);
4106 	
4107 	      nlexprdata->separating = TRUE;
4108 	   }
4109 	   else
4110 	   {
4111 	      SCIPdebugMsg(scip, "expr %p is quadratic and propagable -> propagate only\n", (void*)expr);
4112 	   }
4113 	
4114 	   if( SCIPexprAreQuadraticExprsVariables(expr) )
4115 	   {
4116 	      SCIPexprSetCurvature(expr, nlexprdata->curvature);
4117 	      SCIPdebugMsg(scip, "expr is %s in the original variables\n", nlexprdata->curvature == SCIP_EXPRCURV_CONCAVE ? "concave" : "convex");
4118 	      nlexprdata->origvars = TRUE;
4119 	   }
4120 	
4121 	   return SCIP_OKAY;
4122 	}
4123 	
4124 	/** nonlinear handler auxiliary evaluation callback */
4125 	static
4126 	SCIP_DECL_NLHDLREVALAUX(nlhdlrEvalauxQuadratic)
4127 	{  /*lint --e{715}*/
4128 	   int i;
4129 	   int nlinexprs;
4130 	   int nquadexprs;
4131 	   int nbilinexprs;
4132 	   SCIP_Real constant;
4133 	   SCIP_Real* lincoefs;
4134 	   SCIP_EXPR** linexprs;
4135 	
4136 	   assert(scip != NULL);
4137 	   assert(expr != NULL);
4138 	   assert(auxvalue != NULL);
4139 	   assert(nlhdlrexprdata->separating);
4140 	   assert(nlhdlrexprdata->qexpr == expr);
4141 	
4142 	   /* if the quadratic is in the original variable we can just evaluate the expression */
4143 	   if( nlhdlrexprdata->origvars )
4144 	   {
4145 	      *auxvalue = SCIPexprGetEvalValue(expr);
4146 	      return SCIP_OKAY;
4147 	   }
4148 	
4149 	   /* TODO there was a
4150 	     *auxvalue = SCIPevalExprQuadratic(scip, nlhdlrexprdata->qexpr, sol);
4151 	     here; any reason why not using this anymore?
4152 	   */
4153 	
4154 	   SCIPexprGetQuadraticData(expr, &constant, &nlinexprs, &linexprs, &lincoefs, &nquadexprs, &nbilinexprs, NULL, NULL);
4155 	
4156 	   *auxvalue = constant;
4157 	
4158 	   for( i = 0; i < nlinexprs; ++i ) /* linear exprs */
4159 	      *auxvalue += lincoefs[i] * SCIPgetSolVal(scip, sol, SCIPgetExprAuxVarNonlinear(linexprs[i]));
4160 	
4161 	   for( i = 0; i < nquadexprs; ++i ) /* quadratic terms */
4162 	   {
4163 	      SCIP_Real solval;
4164 	      SCIP_Real lincoef;
4165 	      SCIP_Real sqrcoef;
4166 	      SCIP_EXPR* qexpr;
4167 	
4168 	      SCIPexprGetQuadraticQuadTerm(expr, i, &qexpr, &lincoef, &sqrcoef, NULL, NULL, NULL);
4169 	
4170 	      solval = SCIPgetSolVal(scip, sol, SCIPgetExprAuxVarNonlinear(qexpr));
4171 	      *auxvalue += (lincoef + sqrcoef * solval) * solval;
4172 	   }
4173 	
4174 	   for( i = 0; i < nbilinexprs; ++i ) /* bilinear terms */
4175 	   {
4176 	      SCIP_EXPR* expr1;
4177 	      SCIP_EXPR* expr2;
4178 	      SCIP_Real coef;
4179 	
4180 	      SCIPexprGetQuadraticBilinTerm(expr, i, &expr1, &expr2, &coef, NULL, NULL);
4181 	
4182 	      *auxvalue += coef * SCIPgetSolVal(scip, sol, SCIPgetExprAuxVarNonlinear(expr1)) * SCIPgetSolVal(scip, sol,
4183 	            SCIPgetExprAuxVarNonlinear(expr2));
4184 	   }
4185 	
4186 	   return SCIP_OKAY;
4187 	}
4188 	
4189 	/** nonlinear handler enforcement callback */
4190 	static
4191 	SCIP_DECL_NLHDLRENFO(nlhdlrEnfoQuadratic)
4192 	{  /*lint --e{715}*/
4193 	   SCIP_NLHDLRDATA* nlhdlrdata;
4194 	   SCIP_ROWPREP* rowprep;
4195 	   SCIP_Bool success = FALSE;
4196 	   SCIP_NODE* node;
4197 	   int depth;
4198 	   SCIP_Longint nodenumber;
4199 	   SCIP_Real* eigenvalues;
4200 	   SCIP_Real violation;
4201 	
4202 	   assert(nlhdlrexprdata != NULL);
4203 	   assert(nlhdlrexprdata->qexpr == expr);
4204 	
4205 	   INTERLOG(printf("Starting interesection cuts!\n");)
4206 	
4207 	   nlhdlrdata = SCIPnlhdlrGetData(nlhdlr);
4208 	   assert(nlhdlrdata != NULL);
4209 	
4210 	   assert(result != NULL);
4211 	   *result = SCIP_DIDNOTRUN;
4212 	
4213 	   /* estimate should take care of convex and concave quadratics */
4214 	   if( nlhdlrexprdata->curvature == SCIP_EXPRCURV_CONCAVE || nlhdlrexprdata->curvature == SCIP_EXPRCURV_CONVEX )
4215 	   {
4216 	      INTERLOG(printf("Convex or concave, no need of interesection cuts!\n");)
4217 	      return SCIP_OKAY;
4218 	   }
4219 	
4220 	   /* nothing to do if we can't use intersection cuts */
4221 	   if( ! nlhdlrdata->useintersectioncuts )
4222 	   {
4223 	      INTERLOG(printf("We don't use intersection cuts!\n");)
4224 	      return SCIP_OKAY;
4225 	   }
4226 	
4227 	   /* right now can use interesction cuts only if a basic LP solution is at hand; TODO: in principle we can do something
4228 	    * even if it is not optimal
4229 	    */
4230 	   if( sol != NULL || SCIPgetLPSolstat(scip) != SCIP_LPSOLSTAT_OPTIMAL || !SCIPisLPSolBasic(scip) )
4231 	   {
4232 	      INTERLOG(printf("LP solutoin not good!\n");)
4233 	      return SCIP_OKAY;
4234 	   }
4235 	
4236 	   /* only separate at selected nodes */
4237 	   node = SCIPgetCurrentNode(scip);
4238 	   depth = SCIPnodeGetDepth(node);
4239 	   if( (nlhdlrdata->atwhichnodes == -1 && depth != 0) || (nlhdlrdata->atwhichnodes != -1 && depth % nlhdlrdata->atwhichnodes != 0) )
4240 	   {
4241 	      INTERLOG(printf("Don't separate at this node\n");)
4242 	      return SCIP_OKAY;
4243 	   }
4244 	
4245 	   /* do not add more than ncutslimitroot cuts in root node and ncutslimit cuts in the non-root nodes */
4246 	   nodenumber = SCIPnodeGetNumber(node);
4247 	   if( nlhdlrdata->lastnodenumber != nodenumber )
4248 	   {
4249 	      nlhdlrdata->lastnodenumber = nodenumber;
4250 	      nlhdlrdata->lastncuts = nlhdlrdata->ncutsadded;
4251 	   }
4252 	   /*else if( (depth > 0 && nlhdlrdata->ncutsadded - nlhdlrdata->lastncuts >= nlhdlrdata->ncutslimit) || (depth == 0 &&
4253 	            nlhdlrdata->ncutsadded - nlhdlrdata->lastncuts >= nlhdlrdata->ncutslimitroot)) */
4254 	   /* allow the addition of a certain number of cuts per quadratic */
4255 	   if( (depth > 0 && nlhdlrexprdata->ncutsadded >= nlhdlrdata->ncutslimit) || (depth == 0 &&
4256 	      nlhdlrexprdata->ncutsadded >= nlhdlrdata->ncutslimitroot) )
4257 	   {
4258 	      INTERLOG(printf("Too many cuts added already\n");)
4259 	      return SCIP_OKAY;
4260 	   }
4261 	
4262 	   /* can't separate if we do not have an eigendecomposition */
4263 	   SCIPexprGetQuadraticData(expr, NULL, NULL, NULL, NULL, NULL, NULL, &eigenvalues, NULL);
4264 	   if( eigenvalues == NULL )
4265 	   {
4266 	      INTERLOG(printf("No known eigenvalues!\n");)
4267 	      return SCIP_OKAY;
4268 	   }
4269 	
4270 	   /* if constraint is not sufficiently violated -> do nothing */
4271 	   if( cons != nlhdlrexprdata->cons )
4272 	   {
4273 	      /* constraint is w.r.t auxvar */
4274 	      violation = auxvalue - SCIPgetSolVal(scip, NULL, SCIPgetExprAuxVarNonlinear(expr));
4275 	      violation = ABS( violation );
4276 	   }
4277 	   else
4278 	      /* quadratic is a constraint */
4279 	      violation = MAX( SCIPgetLhsNonlinear(nlhdlrexprdata->cons) - auxvalue, auxvalue -
4280 	            SCIPgetRhsNonlinear(nlhdlrexprdata->cons)); /*lint !e666*/
4281 	
4282 	   if( violation < nlhdlrdata->minviolation )
4283 	   {
4284 	      INTERLOG(printf("Violation %g is just too small\n", violation); )
4285 	      return SCIP_OKAY;
4286 	   }
4287 	
4288 	   /* we can't build an intersection cut when the expr is the root of some constraint and also a subexpression of
4289 	    * another constraint because we initialize data differently */
4290 	   if( nlhdlrexprdata->cons != NULL && cons != nlhdlrexprdata->cons )
4291 	   {
4292 	      INTERLOG(printf("WARNING!! expr is root of one constraint and subexpr of another!\n"); )
4293 	      return SCIP_OKAY;
4294 	   }
4295 	
4296 	   /* if we are the root of a constraint and we are feasible w.r.t our auxiliary variables, that is, auxvalue is
4297 	    * actually feasible for the sides of the constraint, then do not separate
4298 	    */
4299 	   if( cons == nlhdlrexprdata->cons && ((overestimate && (SCIPgetLhsNonlinear(cons)) - auxvalue < SCIPfeastol(scip)) ||
4300 	            (! overestimate && (auxvalue - SCIPgetRhsNonlinear(cons) < SCIPfeastol(scip)))) )
4301 	   {
4302 	      INTERLOG(printf("We are actually feasible for the sides of the constraint\n"); )
4303 	      return SCIP_OKAY;
4304 	   }
4305 	
4306 	#ifdef  DEBUG_INTERSECTIONCUT
4307 	   SCIPinfoMessage(scip, NULL, "Build intersection cut for \n");
4308 	   if( cons == nlhdlrexprdata->cons )
4309 	   {
4310 	      SCIP_CALL( SCIPprintCons(scip, cons, NULL) );
4311 	      SCIPinfoMessage(scip, NULL, "\n");
4312 	   }
4313 	   else
4314 	   {
4315 	      SCIP_CALL( SCIPprintExpr(scip, expr, NULL) );
4316 	      SCIPinfoMessage(scip, NULL, " == %s\n", SCIPvarGetName(SCIPgetExprAuxVarNonlinear(expr)));
4317 	   }
4318 	   SCIPinfoMessage(scip, NULL, "We need to %sestimate\n", overestimate ? "over" : "under" );
4319 	   SCIPinfoMessage(scip, NULL, "LP sol: \n");
4320 	   SCIP_CALL( SCIPprintTransSol(scip, NULL, NULL, FALSE) );
4321 	#endif
4322 	   *result = SCIP_DIDNOTFIND;
4323 	
4324 	   /* cut (in the nonbasic space) is of the form alpha^T x >= 1 */
4325 	   SCIP_CALL( SCIPcreateRowprep(scip, &rowprep, SCIP_SIDETYPE_LEFT, TRUE) );
4326 	   INTERLOG(printf("Generating inter cut\n"); )
4327 	
4328 	   SCIP_CALL( generateIntercut(scip, expr, nlhdlrdata, nlhdlrexprdata, cons, sol, rowprep, overestimate, &success) );
4329 	   INTERLOG(if( !success) printf("Generation failed\n"); )
4330 	
4331 	   /* we generated something, let us see if it survives the clean up */
4332 	   if( success )
4333 	   {
4334 	      assert(sol == NULL);
4335 	      nlhdlrdata->ncutsgenerated += 1;
4336 	      nlhdlrexprdata->ncutsadded += 1;
4337 	
4338 	      /* merge coefficients that belong to same variable */
4339 	      SCIPmergeRowprepTerms(scip, rowprep);
4340 	
4341 	      /* sparsify cut */
4342 	      if( nlhdlrdata->sparsifycuts )
4343 	         sparsifyIntercut(scip, rowprep);
4344 	
4345 	      SCIP_CALL( SCIPcleanupRowprep(scip, rowprep, sol, nlhdlrdata->mincutviolation, &violation, &success) );
4346 	      INTERLOG(if( !success) printf("Clean up failed\n"); )
4347 	   }
4348 	
4349 	   /* if cut looks good (numerics ok and cutting off solution), then turn into row and add to sepastore */
4350 	   if( success )
4351 	   {
4352 	      SCIP_ROW* row;
4353 	      SCIP_Bool infeasible;
4354 	
4355 	      /* count number of bound cuts */
4356 	      if( nlhdlrdata->useboundsasrays )
4357 	         nlhdlrdata->nboundcuts += 1;
4358 	
4359 	      (void) SCIPsnprintf(SCIProwprepGetName(rowprep), SCIP_MAXSTRLEN, "%s_intersection_quadratic%p_lp%" SCIP_LONGINT_FORMAT,
4360 	         overestimate ? "over" : "under",
4361 	         (void*)expr,
4362 	         SCIPgetNLPs(scip));
4363 	
4364 	      SCIP_CALL( SCIPgetRowprepRowCons(scip, &row, rowprep, cons) );
4365 	
4366 	      /* printf("## New cut\n");
4367 	      printf(" -> found maxquad-free cut <%s>: act=%f, lhs=%f, norm=%f, eff=%f, min=%f, max=%f (range=%f)\n\n",
4368 	            SCIProwGetName(row), SCIPgetRowLPActivity(scip, row), SCIProwGetLhs(row), SCIProwGetNorm(row),
4369 	            SCIPgetCutEfficacy(scip, NULL, row),
4370 	            SCIPgetRowMinCoef(scip, row), SCIPgetRowMaxCoef(scip, row),
4371 	            SCIPgetRowMaxCoef(scip, row)/SCIPgetRowMinCoef(scip, row)); */
4372 	
4373 	      /* intersection cuts can be numerically nasty; we do some extra numerical checks here */
4374 	      /*printf("SCIP DEPTH %d got a cut with violation %g, efficacy %g and r/e %g\n", SCIPgetSubscipDepth(scip),
4375 	       * violation, SCIPgetCutEfficacy(scip, NULL, row), SCIPgetRowMaxCoef(scip, row) / SCIPgetRowMinCoef(scip, row) /
4376 	       * SCIPgetCutEfficacy(scip, NULL, row));
4377 	       */
4378 	      assert(SCIPgetCutEfficacy(scip, NULL, row) > 0.0);
4379 	      if( ! nlhdlrdata->ignorehighre || SCIPgetRowMaxCoef(scip, row) / SCIPgetRowMinCoef(scip, row) / SCIPgetCutEfficacy(scip, NULL, row) < 1e9 )
4380 	      {
4381 	#ifdef SCIP_DEBUG
4382 	         SCIPdebugMsg(scip, "adding cut ");
4383 	         SCIP_CALL( SCIPprintRow(scip, row, NULL) );
4384 	#endif
4385 	
4386 	         SCIP_CALL( SCIPaddRow(scip, row, FALSE, &infeasible) );
4387 	
4388 	         if( infeasible )
4389 	         {
4390 	            *result = SCIP_CUTOFF;
4391 	         }
4392 	         else
4393 	         {
4394 	            *result = SCIP_SEPARATED;
4395 	            nlhdlrdata->cutcoefsum += nlhdlrdata->currentavecutcoef;
4396 	            nlhdlrdata->monoidalimprovementsum += nlhdlrdata->currentavemonoidalimprovement;
4397 	            nlhdlrdata->ncutsadded += 1;
4398 	            nlhdlrdata->densitysum += (SCIP_Real) SCIProwprepGetNVars(rowprep) / (SCIP_Real) SCIPgetNVars(scip);
4399 	            nlhdlrdata->efficacysum += SCIPgetCutEfficacy(scip, NULL, row);
4400 	         }
4401 	      }
4402 	      else
4403 	      {
4404 	         nlhdlrdata->nhighre++;
4405 	      }
4406 	      SCIP_CALL( SCIPreleaseRow(scip, &row) );
4407 	   }
4408 	
4409 	   SCIPfreeRowprep(scip, &rowprep);
4410 	
4411 	   return SCIP_OKAY;
4412 	}
4413 	
4414 	/** nonlinear handler forward propagation callback
4415 	 *
4416 	 * This method should solve the problem
4417 	 * <pre>
4418 	 *    max/min quad expression over box constraints
4419 	 * </pre>
4420 	 * However, this problem is difficult so we are satisfied with a proxy.
4421 	 * Interval arithmetic suffices when no variable appears twice, however this is seldom the case, so we try
4422 	 * to take care of the dependency problem to some extent:
4423 	 * Let \f$P_l = \{i : \text{expr}_l \text{expr}_i \,\text{is a bilinear expr}\}\f$.
4424 	 * 1. partition the quadratic expression as sum of quadratic functions \f$\sum_l q_l\f$
4425 	 *    where \f$q_l = a_l \text{expr}_l^2 + c_l \text{expr}_l + \sum_{i \in P_l} b_{il} \text{expr}_i \text{expr}_l\f$
4426 	 * 2. build interval quadratic functions, i.e., \f$a x^2 + b x\f$ where \f$b\f$ is an interval, i.e.,
4427 	 *    \f$a_l \text{expr}_l^2 + [\sum_{i \in P_l} b_{il} \text{expr}_i + c_l] \text{expr}_l\f$
4428 	 * 3. compute \f$\min/\max \{ a x^2 + b x : x \in [x] \}\f$ for each interval quadratic, i.e.,
4429 	 *    \f$\min/\max a_l \text{expr}_l^2 + \text{expr}_l [\sum_{i \in P_l} b_{il} \text{expr}_i + c_l] : \text{expr}_l \in [\text{expr}_l]\f$
4430 	 *
4431 	 * Notes:
4432 	 * 1. The \f$l\f$-th quadratic expr (expressions that appear quadratically) is associated with \f$q_l\f$.
4433 	 * 2. `nlhdlrdata->quadactivities[l]` is the activity of \f$q_l\f$ as computed in the description above.
4434 	 * 3. The \f$q_l\f$ of a quadratic term might be empty, in which case `nlhdlrdata->quadactivities[l]` is [0,0].\n
4435 	 *    For example, consider \f$x^2 + xy\f$. There are two quadratic expressions, \f$x\f$ and \f$y\f$.
4436 	 *    The \f$q\f$ associated to \f$x\f$ is \f$x^2 + xy\f$, while the \f$q\f$ associated to \f$y\f$ is empty.
4437 	 *    Thus, `nlhdlrdata->quadactivities[1]` is [0,0] in this case.
4438 	 *    The logic is to avoid considering the term \f$xy\f$ twice.
4439 	 *
4440 	 * @note The order matters! If \f$\text{expr}_i\, \text{expr}_l\f$ is a term in the quadratic, then \f$i\f$ is *not* in \f$P_l\f$
4441 	 */
4442 	static
4443 	SCIP_DECL_NLHDLRINTEVAL(nlhdlrIntevalQuadratic)
4444 	{ /*lint --e{715}*/
4445 	   SCIP_EXPR** linexprs;
4446 	   SCIP_Real* lincoefs;
4447 	   SCIP_Real constant;
4448 	   int nquadexprs;
4449 	   int nlinexprs;
4450 	
4451 	   assert(scip != NULL);
4452 	   assert(expr != NULL);
4453 	
4454 	   assert(nlhdlrexprdata != NULL);
4455 	   assert(nlhdlrexprdata->quadactivities != NULL);
4456 	   assert(nlhdlrexprdata->qexpr == expr);
4457 	
4458 	   SCIPdebugMsg(scip, "Interval evaluation of quadratic expr\n");
4459 	
4460 	   SCIPexprGetQuadraticData(expr, &constant, &nlinexprs, &linexprs, &lincoefs, &nquadexprs, NULL, NULL, NULL);
4461 	
4462 	   /*
4463 	    * compute activity of linear part, if some linear term has changed
4464 	    */
4465 	   {
4466 	      int i;
4467 	
4468 	      SCIPdebugMsg(scip, "Computing activity of linear part\n");
4469 	
4470 	      SCIPintervalSet(&nlhdlrexprdata->linactivity, constant);
4471 	      for( i = 0; i < nlinexprs; ++i )
4472 	      {
4473 	         SCIP_INTERVAL linterminterval;
4474 	
4475 	         linterminterval = SCIPexprGetActivity(linexprs[i]);
4476 	         if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, linterminterval) )
4477 	         {
4478 	            SCIPdebugMsg(scip, "Activity of linear part is empty due to child %d\n", i);
4479 	            SCIPintervalSetEmpty(interval);
4480 	            return SCIP_OKAY;
4481 	         }
4482 	         SCIPintervalMulScalar(SCIP_INTERVAL_INFINITY, &linterminterval, linterminterval, lincoefs[i]);
4483 	         SCIPintervalAdd(SCIP_INTERVAL_INFINITY, &nlhdlrexprdata->linactivity, nlhdlrexprdata->linactivity, linterminterval);
4484 	      }
4485 	
4486 	      SCIPdebugMsg(scip, "Activity of linear part is [%g, %g]\n", nlhdlrexprdata->linactivity.inf,
4487 	            nlhdlrexprdata->linactivity.sup);
4488 	   }
4489 	
4490 	   /*
4491 	    * compute activity of quadratic part
4492 	    */
4493 	   {
4494 	      int i;
4495 	
4496 	      SCIPdebugMsg(scip, "Computing activity of quadratic part\n");
4497 	
4498 	      nlhdlrexprdata->nneginfinityquadact = 0;
4499 	      nlhdlrexprdata->nposinfinityquadact = 0;
4500 	      nlhdlrexprdata->minquadfiniteact = 0.0;
4501 	      nlhdlrexprdata->maxquadfiniteact = 0.0;
4502 	      SCIPintervalSet(&nlhdlrexprdata->quadactivity, 0.0);
4503 	
4504 	      for( i = 0; i < nquadexprs; ++i )
4505 	      {
4506 	         SCIP_Real quadlb;
4507 	         SCIP_Real quadub;
4508 	         SCIP_EXPR* qexpr;
4509 	         SCIP_Real lincoef;
4510 	         SCIP_Real sqrcoef;
4511 	         int nadjbilin;
4512 	         int* adjbilin;
4513 	         SCIP_EXPR* sqrexpr;
4514 	
4515 	         SCIPexprGetQuadraticQuadTerm(expr, i, &qexpr, &lincoef, &sqrcoef, &nadjbilin, &adjbilin, &sqrexpr);
4516 	
4517 	         if( !isPropagableTerm(expr, i) )
4518 	         {
4519 	            /* term is not propagable, i.e., the exprs involved in term only appear once; thus use the activity of the
4520 	             * quadratic term directly and not the activity of the exprs involed in the term. See also documentation of
4521 	             * DETECT
4522 	             */
4523 	            SCIP_INTERVAL tmp;
4524 	
4525 	            assert(lincoef == 0.0);
4526 	
4527 	            if( sqrcoef != 0.0 )
4528 	            {
4529 	               assert(sqrexpr != NULL);
4530 	               assert(nadjbilin == 0);
4531 	
4532 	               tmp = SCIPexprGetActivity(sqrexpr);
4533 	               if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, tmp) )
4534 	               {
4535 	                  SCIPintervalSetEmpty(interval);
4536 	                  return SCIP_OKAY;
4537 	               }
4538 	
4539 	               SCIPintervalMulScalar(SCIP_INTERVAL_INFINITY, &tmp, tmp, sqrcoef);
4540 	               quadlb = tmp.inf;
4541 	               quadub = tmp.sup;
4542 	
4543 	#ifdef DEBUG_PROP
4544 	               SCIPinfoMessage(scip, NULL, "Computing activity for quadratic term %g <expr>, where <expr> is: ", sqrcoef);
4545 	               SCIP_CALL( SCIPprintExpr(scip, sqrexpr, NULL) );
4546 	#endif
4547 	            }
4548 	            else
4549 	            {
4550 	               SCIP_EXPR* expr1;
4551 	               SCIP_EXPR* prodexpr;
4552 	               SCIP_Real prodcoef;
4553 	
4554 	               assert(nadjbilin == 1);
4555 	               SCIPexprGetQuadraticBilinTerm(expr, adjbilin[0], &expr1, NULL, &prodcoef, NULL, &prodexpr);
4556 	
4557 	               if( expr1 == qexpr )
4558 	               {
4559 	                  /* the quadratic expression expr1 appears only as expr1 * expr2, so its 'q' is expr1 * expr2 */
4560 	                  tmp = SCIPexprGetActivity(prodexpr);
4561 	                  if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, tmp) )
4562 	                  {
4563 	                     SCIPintervalSetEmpty(interval);
4564 	                     return SCIP_OKAY;
4565 	                  }
4566 	
4567 	                  SCIPintervalMulScalar(SCIP_INTERVAL_INFINITY, &tmp, tmp, prodcoef);
4568 	                  quadlb = tmp.inf;
4569 	                  quadub = tmp.sup;
4570 	
4571 	#ifdef DEBUG_PROP
4572 	                  SCIPinfoMessage(scip, NULL, "Computing activity for quadratic term %g <expr>, where <expr> is: ", prodcoef);
4573 	                  SCIP_CALL( SCIPprintExpr(scip, prodexpr, NULL) );
4574 	#endif
4575 	               }
4576 	               else
4577 	               {
4578 	                  /* the quadratic expression expr1 appears as expr2 * expr1, thus its 'q' is empty, see also the Notes
4579 	                   * in the documentation of the function
4580 	                   */
4581 	                  SCIPintervalSet(&nlhdlrexprdata->quadactivities[i], 0.0);
4582 	                  continue;
4583 	               }
4584 	            }
4585 	         }
4586 	         else
4587 	         {
4588 	            int j;
4589 	            SCIP_INTERVAL b;
4590 	
4591 	            SCIPexprGetQuadraticQuadTerm(expr, i, &qexpr, &lincoef, &sqrcoef, &nadjbilin, &adjbilin, NULL);
4592 	
4593 	            if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, SCIPexprGetActivity(qexpr)) )
4594 	            {
4595 	               SCIPintervalSetEmpty(interval);
4596 	               return SCIP_OKAY;
4597 	            }
4598 	
4599 	            /* b = [c_l] */
4600 	            SCIPintervalSet(&b, lincoef);
4601 	#ifdef DEBUG_PROP
4602 	            SCIPinfoMessage(scip, NULL, "b := %g\n", lincoef);
4603 	#endif
4604 	            for( j = 0; j < nadjbilin; ++j )
4605 	            {
4606 	               SCIP_INTERVAL bterm;
4607 	               SCIP_EXPR* expr1;
4608 	               SCIP_EXPR* expr2;
4609 	               SCIP_Real bilincoef;
4610 	
4611 	               SCIPexprGetQuadraticBilinTerm(expr, adjbilin[j], &expr1, &expr2, &bilincoef, NULL, NULL);
4612 	
4613 	               if( expr1 != qexpr )
4614 	                  continue;
4615 	
4616 	               bterm = SCIPexprGetActivity(expr2);
4617 	               if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, bterm) )
4618 	               {
4619 	                  SCIPintervalSetEmpty(interval);
4620 	                  return SCIP_OKAY;
4621 	               }
4622 	
4623 	               /* b += [b_jl * expr_j] for j \in P_l */
4624 	               SCIPintervalMulScalar(SCIP_INTERVAL_INFINITY, &bterm, bterm, bilincoef);
4625 	               SCIPintervalAdd(SCIP_INTERVAL_INFINITY, &b, b, bterm);
4626 	
4627 	#ifdef DEBUG_PROP
4628 	               SCIPinfoMessage(scip, NULL, "b += %g * [expr2], where <expr2> is: ", bilincoef);
4629 	               SCIP_CALL( SCIPprintExpr(scip, expr2, NULL) );
4630 	               SCIPinfoMessage(scip, NULL, " [%g,%g]\n", SCIPexprGetActivity(expr2).inf, SCIPexprGetActivity(expr2).sup);
4631 	#endif
4632 	            }
4633 	
4634 	            /* TODO: under which assumptions do we know that we just need to compute min or max? its probably the locks that give some information here */
4635 	            quadub = SCIPintervalQuadUpperBound(SCIP_INTERVAL_INFINITY, sqrcoef, b,
4636 	               SCIPexprGetActivity(qexpr));
4637 	
4638 	            /* TODO: implement SCIPintervalQuadLowerBound */
4639 	            {
4640 	               SCIP_INTERVAL minusb;
4641 	               SCIPintervalSetBounds(&minusb, -SCIPintervalGetSup(b), -SCIPintervalGetInf(b));
4642 	
4643 	               quadlb = -SCIPintervalQuadUpperBound(SCIP_INTERVAL_INFINITY, -sqrcoef, minusb,
4644 	                  SCIPexprGetActivity(qexpr));
4645 	            }
4646 	
4647 	#ifdef DEBUG_PROP
4648 	            SCIPinfoMessage(scip, NULL, "Computing activity for quadratic term %g <expr>^2 + [%g,%g] <expr>, where <expr> is: ", sqrcoef, b.inf, b.sup);
4649 	            SCIP_CALL( SCIPprintExpr(scip, qexpr, NULL) );
4650 	#endif
4651 	         }
4652 	#ifdef DEBUG_PROP
4653 	         SCIPinfoMessage(scip, NULL, " -> [%g, %g]\n", quadlb, quadub);
4654 	#endif
4655 	
4656 	         SCIPintervalSetBounds(&nlhdlrexprdata->quadactivities[i], quadlb, quadub);
4657 	         SCIPintervalAdd(SCIP_INTERVAL_INFINITY, &nlhdlrexprdata->quadactivity, nlhdlrexprdata->quadactivity, nlhdlrexprdata->quadactivities[i]);
4658 	
4659 	         /* get number of +/-infinity contributions and compute finite activity */
4660 	         if( quadlb <= -SCIP_INTERVAL_INFINITY )
4661 	            nlhdlrexprdata->nneginfinityquadact++;
4662 	         else
4663 	         {
4664 	            SCIP_ROUNDMODE roundmode;
4665 	
4666 	            roundmode = SCIPintervalGetRoundingMode();
4667 	            SCIPintervalSetRoundingModeDownwards();
4668 	
4669 	            nlhdlrexprdata->minquadfiniteact += quadlb;
4670 	
4671 	            SCIPintervalSetRoundingMode(roundmode);
4672 	         }
4673 	         if( quadub >= SCIP_INTERVAL_INFINITY )
4674 	            nlhdlrexprdata->nposinfinityquadact++;
4675 	         else
4676 	         {
4677 	            SCIP_ROUNDMODE roundmode;
4678 	
4679 	            roundmode = SCIPintervalGetRoundingMode();
4680 	            SCIPintervalSetRoundingModeUpwards();
4681 	
4682 	            nlhdlrexprdata->maxquadfiniteact += quadub;
4683 	
4684 	            SCIPintervalSetRoundingMode(roundmode);
4685 	         }
4686 	      }
4687 	
4688 	      SCIPdebugMsg(scip, "Activity of quadratic part is [%g, %g]\n", nlhdlrexprdata->quadactivity.inf, nlhdlrexprdata->quadactivity.sup);
4689 	   }
4690 	
4691 	   /* interval evaluation is linear activity + quadactivity */
4692 	   SCIPintervalAdd(SCIP_INTERVAL_INFINITY, interval, nlhdlrexprdata->linactivity,  nlhdlrexprdata->quadactivity);
4693 	
4694 	   nlhdlrexprdata->activitiestag = SCIPgetCurBoundsTagNonlinear(SCIPfindConshdlr(scip, "nonlinear"));
4695 	
4696 	   return SCIP_OKAY;
4697 	}
4698 	
4699 	/** nonlinear handler reverse propagation callback
4700 	 *
4701 	 * @note the implemented technique is a proxy for solving the problem min/max{ x_i : quad expr in [quad expr] }
4702 	 * and as such can be improved.
4703 	 */
4704 	static
4705 	SCIP_DECL_NLHDLRREVERSEPROP(nlhdlrReversepropQuadratic)
4706 	{ /*lint --e{715}*/
4707 	   SCIP_EXPR** linexprs;
4708 	   SCIP_EXPR** bilinexprs; /* TODO: should this be stored in the nlhdlr expr data? */
4709 	   SCIP_Real* bilincoefs;
4710 	   SCIP_Real* lincoefs;
4711 	   SCIP_Real constant;
4712 	   int nquadexprs;
4713 	   int nlinexprs;
4714 	
4715 	   SCIP_INTERVAL rhs;
4716 	   SCIP_INTERVAL quadactivity;
4717 	   int i;
4718 	
4719 	   SCIPdebugMsg(scip, "Reverse propagation of quadratic expr given bounds = [%g,%g]\n", bounds.inf, bounds.sup);
4720 	
4721 	   assert(scip != NULL);
4722 	   assert(expr != NULL);
4723 	   assert(infeasible != NULL);
4724 	   assert(nreductions != NULL);
4725 	   assert(nlhdlrexprdata != NULL);
4726 	   assert(nlhdlrexprdata->quadactivities != NULL);
4727 	   assert(nlhdlrexprdata->qexpr == expr);
4728 	
4729 	   *nreductions = 0;
4730 	
4731 	   /* not possible to conclude finite bounds if the interval of the expression is [-inf,inf] */
4732 	   if( SCIPintervalIsEntire(SCIP_INTERVAL_INFINITY, bounds) )
4733 	   {
4734 	      SCIPdebugMsg(scip, "expr's range is R -> cannot reverse propagate\n");
4735 	      return SCIP_OKAY;
4736 	   }
4737 	
4738 	   /* ensure that partial activities as stored in nlhdlrexprdata are uptodate
4739 	    * if the activity stored in expr is more recent than the partial activities stored in this nlhdlrexprdata,
4740 	    * then we should reevaluate the partial activities
4741 	    */
4742 	   if( SCIPexprGetActivityTag(expr) > nlhdlrexprdata->activitiestag )
4743 	   {
4744 	      SCIP_CALL( nlhdlrIntevalQuadratic(scip, nlhdlr, expr, nlhdlrexprdata, &quadactivity, NULL, NULL) );
4745 	   }
4746 	
4747 	   SCIPexprGetQuadraticData(expr, &constant, &nlinexprs, &linexprs, &lincoefs, &nquadexprs, NULL, NULL, NULL);
4748 	
4749 	   /* propagate linear part in rhs = expr's interval - quadratic activity; first, reconstruct the quadratic activity */
4750 	   SCIPintervalSetBounds(&quadactivity,
4751 	         nlhdlrexprdata->nneginfinityquadact > 0 ? -SCIP_INTERVAL_INFINITY : nlhdlrexprdata->minquadfiniteact,
4752 	         nlhdlrexprdata->nposinfinityquadact > 0 ?  SCIP_INTERVAL_INFINITY : nlhdlrexprdata->maxquadfiniteact);
4753 	
4754 	   SCIPintervalSub(SCIP_INTERVAL_INFINITY, &rhs, bounds, quadactivity);
4755 	
4756 	   SCIP_CALL( reversePropagateLinearExpr(scip, linexprs, nlinexprs, lincoefs, constant, rhs, infeasible, nreductions) );
4757 	
4758 	   /* stop if we find infeasibility */
4759 	   if( *infeasible )
4760 	      return SCIP_OKAY;
4761 	
4762 	   /* propagate quadratic part in expr's interval - linear activity, where linear activity was computed in INTEVAL.
4763 	    * The idea is basically to write interval quadratics for each expr and then solve for expr.
4764 	    *
4765 	    * One way of achieving this is:
4766 	    * - for each expression expr_i, write the quadratic expression as a_i expr^2_i + expr_i ( \sum_{j \in J_i} b_ij
4767 	    *   expr_j + c_i ) + quadratic expression in expr_k for k \neq i
4768 	    * - compute the interval b = [\sum_{j \in J_i} b_ij expr_j + c_i], where J_i are all the indices j such that the
4769 	    *   bilinear expression expr_i expr_j appears
4770 	    * - use some technique (like the one in nlhdlrIntevalQuadratic), to evaluate the activity of rest_i = [quadratic
4771 	    *   expression in expr_k for k \neq i].
4772 	    * - solve a_i expr_i^2 + b expr_i \in rhs_i := [expr activity] - rest_i
4773 	    *
4774 	    * However, this might be expensive, especially computing rest_i. Hence, we implement a simpler version.
4775 	    * - we use the same partition as in nlhdlrIntevalQuadratic for the bilinear terms. This way, b = [\sum_{j \in P_i}
4776 	    *   b_ij expr_j + c_i], where P_i is the set of indices j such that expr_i * expr_j appears in that order
4777 	    * - we evaluate the activity of rest_i as sum_{k \neq i} [\min q_k, \max q_k] where q_k = a_k expr_k^2 + [\sum_{j
4778 	    *   \in P_k} b_jk expr_j + c_k] expr_k. The intervals [\min q_k, \max q_k] were already computed in
4779 	    *   nlhdlrIntevalQuadratic, so we just reuse them.
4780 	    *
4781 	    * A downside of the above is that we might not deduce any bounds for variables that appear less often. For example,
4782 	    * consider x^2 + x * y + x * z + y * z + z. This quadratic gets partitioned as (x^2 + x*y + x*z) + (z*y + z). The
4783 	    * first parenthesis is interpreted as a function of x, while the second one as a function of z.
4784 	    * To also get bounds on y, after reverse propagating x in x^2 + x*y + x*z \in rhs, we rewrite this as y + z \in rhs/x -
4785 	    * x and propagate the y + z).
4786 	    * In general, after reverse propagating expr_i, we consider
4787 	    *   \sum_{j \in J_i} b_ij expr_j in ([expr activity] - quadratic expression in expr_k for k \neq i - c_i) / expr_i - a_i expr_i,
4788 	    * compute an interval for the right hand side (see computeRangeForBilinearProp) and use that to propagate the
4789 	    * linear sum on the left hand side.
4790 	    *
4791 	    * Note: this last step generalizes a technique that appeared in the classic cons_quadratic.
4792 	    * The idea of that technique was to borrow a bilinear term expr_k expr_l when propagating expr_l and the quadratic
4793 	    * function for expr_k was simple enough.
4794 	    * Since in P_l we only consider the indices of expressions that appear multiplying expr_l as _second_ factor, we
4795 	    * would lose the bilinear terms expr_k * expr_l, which contributes to the dependency problem.
4796 	    * The problem is that the contribution of b_kl * expr_k * expr_l to rest_i is not just [b_kl * expr_k * expr_l], but
4797 	    * rather quadactivities[k] (= max/min of a_k expr_k^2 + expr_k * [c_k + sum_i \in P_k b_ki expr_i]).
4798 	    * Thus, we _cannot_ just substract [b_kl * expr_k * expr_l] from rest_i.
4799 	    * But, if expr_k only appears as expr_k * expr_l, then  quadactivities[k] = [b_kl * expr_k * expr_l]. So this
4800 	    * case was handled in old cons_quadratic.
4801 	    *
4802 	    *
4803 	    * TODO: handle simple cases
4804 	    * TODO: identify early when there is nothing to be gain
4805 	    */
4806 	   SCIPintervalSub(SCIP_INTERVAL_INFINITY, &rhs, bounds, nlhdlrexprdata->linactivity);
4807 	   SCIP_CALL( SCIPallocBufferArray(scip, &bilinexprs, nquadexprs) );
4808 	   SCIP_CALL( SCIPallocBufferArray(scip, &bilincoefs, nquadexprs) );
4809 	
4810 	   for( i = 0; i < nquadexprs; ++i )
4811 	   {
4812 	      SCIP_INTERVAL rhs_i;
4813 	      SCIP_INTERVAL rest_i;
4814 	      SCIP_EXPR* qexpr;
4815 	      SCIP_Real lincoef;
4816 	      SCIP_Real sqrcoef;
4817 	      int nadjbilin;
4818 	      int* adjbilin;
4819 	      SCIP_EXPR* sqrexpr;
4820 	
4821 	      SCIPexprGetQuadraticQuadTerm(expr, i, &qexpr, &lincoef, &sqrcoef, &nadjbilin, &adjbilin, &sqrexpr);
4822 	
4823 	      /* rhs_i = rhs - rest_i.
4824 	       * to compute rest_i = [\sum_{k \neq i} q_k] we just have to substract
4825 	       * the activity of q_i from quadactivity; however, care must be taken about infinities;
4826 	       * if [q_i].sup = +infinity and there is = 1 contributing +infinity -> rest_i.sup = maxquadfiniteact
4827 	       * if [q_i].sup = +infinity and there is > 1 contributing +infinity -> rest_i.sup = +infinity
4828 	       * if [q_i].sup = finite and there is > 0 contributing +infinity -> rest_i.sup = +infinity
4829 	       * if [q_i].sup = finite and there is = 0 contributing +infinity -> rest_i.sup = maxquadfiniteact - [q_i].sup
4830 	       *
4831 	       * the same holds when replacing sup with inf, + with - and max(quadfiniteact) with min(...)
4832 	       */
4833 	      /* compute rest_i.sup */
4834 	      if( SCIPintervalGetSup(nlhdlrexprdata->quadactivities[i]) < SCIP_INTERVAL_INFINITY &&
4835 	         nlhdlrexprdata->nposinfinityquadact == 0 )
4836 	      {
4837 	         SCIP_ROUNDMODE roundmode;
4838 	
4839 	         roundmode = SCIPintervalGetRoundingMode();
4840 	         SCIPintervalSetRoundingModeUpwards();
4841 	         rest_i.sup = nlhdlrexprdata->maxquadfiniteact - SCIPintervalGetSup(nlhdlrexprdata->quadactivities[i]);
4842 	
4843 	         SCIPintervalSetRoundingMode(roundmode);
4844 	      }
4845 	      else if( SCIPintervalGetSup(nlhdlrexprdata->quadactivities[i]) >= SCIP_INTERVAL_INFINITY &&
4846 	         nlhdlrexprdata->nposinfinityquadact == 1 )
4847 	         rest_i.sup = nlhdlrexprdata->maxquadfiniteact;
4848 	      else
4849 	         rest_i.sup = SCIP_INTERVAL_INFINITY;
4850 	
4851 	      /* compute rest_i.inf */
4852 	      if( SCIPintervalGetInf(nlhdlrexprdata->quadactivities[i]) > -SCIP_INTERVAL_INFINITY &&
4853 	         nlhdlrexprdata->nneginfinityquadact == 0 )
4854 	      {
4855 	         SCIP_ROUNDMODE roundmode;
4856 	
4857 	         roundmode = SCIPintervalGetRoundingMode();
4858 	         SCIPintervalSetRoundingModeDownwards();
4859 	         rest_i.inf = nlhdlrexprdata->minquadfiniteact - SCIPintervalGetInf(nlhdlrexprdata->quadactivities[i]);
4860 	
4861 	         SCIPintervalSetRoundingMode(roundmode);
4862 	      }
4863 	      else if( SCIPintervalGetInf(nlhdlrexprdata->quadactivities[i]) <= -SCIP_INTERVAL_INFINITY &&
4864 	         nlhdlrexprdata->nneginfinityquadact == 1 )
4865 	         rest_i.inf = nlhdlrexprdata->minquadfiniteact;
4866 	      else
4867 	         rest_i.inf = -SCIP_INTERVAL_INFINITY;
4868 	
4869 	#ifdef SCIP_DISABLED_CODE  /* I (SV) added the following in cons_quadratic to fix/workaround some bug. Maybe we'll need this here, too? */
4870 	      /* FIXME in theory, rest_i should not be empty here
4871 	       * what we tried to do here is to remove the contribution of the i'th bilinear term (=bilinterm) to [minquadactivity,maxquadactivity] from rhs
4872 	       * however, quadactivity is computed differently (as x*(a1*y1+...+an*yn)) than q_i (a*ak*yk) and since interval arithmetics do overestimation,
4873 	       * it can happen that q_i is actually slightly larger than quadactivity, which results in rest_i being (slightly) empty
4874 	       * a proper fix could be to compute the quadactivity also as x*a1*y1+...+x*an*yn if sqrcoef=0, but due to taking
4875 	       * also infinite bounds into account, this complicates the code even further
4876 	       * instead, I'll just work around this by turning an empty rest_i into a small non-empty one
4877 	       */
4878 	      if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, rest_i) )
4879 	      {
4880 	         assert(SCIPisSumRelEQ(scip, rest_i.inf, rest_i.sup));
4881 	         SCIPswapReals(&rest_i.inf, &rest_i.sup);
4882 	      }
4883 	#endif
4884 	      assert(!SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, rest_i));
4885 	
4886 	      /* compute rhs_i */
4887 	      SCIPintervalSub(SCIP_INTERVAL_INFINITY, &rhs_i, rhs, rest_i);
4888 	
4889 	      if( SCIPintervalIsEntire(SCIP_INTERVAL_INFINITY, rhs_i) )
4890 	         continue;
4891 	
4892 	      /* try to propagate */
4893 	      if( !isPropagableTerm(expr, i) )
4894 	      {
4895 	         assert(lincoef == 0.0);
4896 	
4897 	         if( sqrcoef != 0.0 )
4898 	         {
4899 	            assert(sqrexpr != NULL);
4900 	            assert(nadjbilin == 0);
4901 	
4902 	            /* solve sqrcoef sqrexpr in rhs_i */
4903 	            SCIP_CALL( propagateBoundsLinExpr(scip, sqrexpr, sqrcoef, rhs_i, infeasible, nreductions) );
4904 	         }
4905 	         else
4906 	         {
4907 	            /* qexpr only appears in a term of the form qexpr * other_expr (or other_expr * qexpr); we only care about
4908 	             * getting bounds for the product, thus we will compute these bounds when qexpr appears as qexpr *
4909 	             * other_expr; note that if it appears as other_expr * qexpr, then when we process other_expr bounds for the
4910 	             * product will be computed
4911 	             * TODO: we can actually avoid computing rhs_i in the case that qexpr is not propagable and it appears as
4912 	             * other_expr * qexpr
4913 	             */
4914 	            SCIP_EXPR* expr1;
4915 	            SCIP_EXPR* prodexpr;
4916 	            SCIP_Real prodcoef;
4917 	
4918 	            assert(nadjbilin == 1);
4919 	            SCIPexprGetQuadraticBilinTerm(expr, adjbilin[0], &expr1, NULL, &prodcoef, NULL, &prodexpr);
4920 	
4921 	            if( expr1 == qexpr )
4922 	            {
4923 	               /* solve prodcoef prodexpr in rhs_i */
4924 	               SCIP_CALL( propagateBoundsLinExpr(scip, prodexpr, prodcoef, rhs_i, infeasible, nreductions) );
4925 	            }
4926 	         }
4927 	      }
4928 	      else
4929 	      {
4930 	         SCIP_INTERVAL b;
4931 	         SCIP_EXPR* expr1 = NULL;
4932 	         SCIP_EXPR* expr2 = NULL;
4933 	         SCIP_Real bilincoef = 0.0;
4934 	         int nbilin = 0;
4935 	         int pos2 = 0;
4936 	         int j;
4937 	
4938 	         /* set b to [c_l] */
4939 	         SCIPintervalSet(&b, lincoef);
4940 	
4941 	         /* add [\sum_{j \in P_l} b_lj expr_j + c_l] into b */
4942 	         for( j = 0; j < nadjbilin; ++j )
4943 	         {
4944 	            SCIP_INTERVAL bterm;
4945 	            SCIP_INTERVAL expr2bounds;
4946 	
4947 	            SCIPexprGetQuadraticBilinTerm(expr, adjbilin[j], &expr1, &expr2, &bilincoef, &pos2, NULL);
4948 	
4949 	            if( expr1 != qexpr )
4950 	               continue;
4951 	
4952 	            expr2bounds = SCIPgetExprBoundsNonlinear(scip, expr2);
4953 	            if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, expr2bounds) )
4954 	            {
4955 	               *infeasible = TRUE;
4956 	               break;
4957 	            }
4958 	
4959 	            /* b += [b_lj * expr_j] for j \in P_l */
4960 	            SCIPintervalMulScalar(SCIP_INTERVAL_INFINITY, &bterm, expr2bounds, bilincoef);
4961 	            SCIPintervalAdd(SCIP_INTERVAL_INFINITY, &b, b, bterm);
4962 	
4963 	            /* remember b_lj and expr_j to propagate them too */
4964 	            bilinexprs[nbilin] = expr2;
4965 	            bilincoefs[nbilin] = bilincoef;
4966 	            nbilin++;
4967 	         }
4968 	
4969 	         if( !*infeasible )
4970 	         {
4971 	            /* solve a_i expr_i^2 + b expr_i in rhs_i */
4972 	            SCIP_CALL( propagateBoundsQuadExpr(scip, qexpr, sqrcoef, b, rhs_i, infeasible, nreductions) );
4973 	         }
4974 	
4975 	         if( nbilin > 0 && !*infeasible )
4976 	         {
4977 	            /* if 0 is not in [expr_i], then propagate bilincoefs^T bilinexpr in rhs_i/expr_i - a_i expr_i - c_i */
4978 	            SCIP_INTERVAL bilinrhs;
4979 	            SCIP_INTERVAL qexprbounds;
4980 	
4981 	            qexprbounds = SCIPgetExprBoundsNonlinear(scip, qexpr);
4982 	            if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, qexprbounds) )
4983 	            {
4984 	               *infeasible = TRUE;
4985 	            }
4986 	            else
4987 	            {
4988 	               /* compute bilinrhs := [rhs_i/expr_i - a_i expr_i] */
4989 	               computeRangeForBilinearProp(qexprbounds, sqrcoef, rhs_i, &bilinrhs);
4990 	
4991 	               if( !SCIPintervalIsEntire(SCIP_INTERVAL_INFINITY, bilinrhs) )
4992 	               {
4993 	                  int nreds;
4994 	
4995 	                  /* propagate \sum_{j \in P_i} b_ij expr_j + c_i in bilinrhs */
4996 	                  SCIP_CALL( reversePropagateLinearExpr(scip, bilinexprs, nbilin, bilincoefs, lincoef, bilinrhs,
4997 	                           infeasible, &nreds) );
4998 	
4999 	                  /* TODO FIXME: we are overestimating the number of reductions: an expr might be tightened many times! */
5000 	                  *nreductions += nreds;
5001 	               }
5002 	            }
5003 	         }
5004 	      }
5005 	
5006 	      /* stop if we find infeasibility */
5007 	      if( *infeasible )
5008 	         break;
5009 	   }
5010 	
5011 	   SCIPfreeBufferArray(scip, &bilincoefs);
5012 	   SCIPfreeBufferArray(scip, &bilinexprs);
5013 	
5014 	   return SCIP_OKAY;
5015 	}
5016 	
5017 	/** callback to free data of handler */
5018 	static
5019 	SCIP_DECL_NLHDLRFREEHDLRDATA(nlhdlrFreehdlrdataQuadratic)
5020 	{ /*lint --e{715}*/
5021 	   assert(nlhdlrdata != NULL);
5022 	
5023 	   SCIPfreeBlockMemory(scip, nlhdlrdata);
5024 	
5025 	   return SCIP_OKAY;
5026 	}
5027 	
5028 	/** nonlinear handler copy callback */
5029 	static
5030 	SCIP_DECL_NLHDLRCOPYHDLR(nlhdlrCopyhdlrQuadratic)
5031 	{  /*lint --e{715}*/
5032 	   assert(targetscip != NULL);
5033 	   assert(sourcenlhdlr != NULL);
5034 	   assert(strcmp(SCIPnlhdlrGetName(sourcenlhdlr), NLHDLR_NAME) == 0);
5035 	
5036 	   SCIP_CALL( SCIPincludeNlhdlrQuadratic(targetscip) );
5037 	
5038 	   return SCIP_OKAY;
5039 	}
5040 	
5041 	/** includes quadratic nonlinear handler in nonlinear constraint handler */
5042 	SCIP_RETCODE SCIPincludeNlhdlrQuadratic(
5043 	   SCIP*                 scip                /**< SCIP data structure */
5044 	   )
5045 	{
5046 	   SCIP_NLHDLRDATA* nlhdlrdata;
5047 	   SCIP_NLHDLR* nlhdlr;
5048 	
5049 	   assert(scip != NULL);
5050 	
5051 	   /* create nonlinear handler specific data */
5052 	   SCIP_CALL( SCIPallocBlockMemory(scip, &nlhdlrdata) );
5053 	   BMSclearMemory(nlhdlrdata);
5054 	
5055 	   SCIP_CALL( SCIPincludeNlhdlrNonlinear(scip, &nlhdlr, NLHDLR_NAME, NLHDLR_DESC, NLHDLR_DETECTPRIORITY,
5056 	      NLHDLR_ENFOPRIORITY, nlhdlrDetectQuadratic, nlhdlrEvalauxQuadratic, nlhdlrdata) );
5057 	
5058 	   SCIPnlhdlrSetCopyHdlr(nlhdlr, nlhdlrCopyhdlrQuadratic);
5059 	   SCIPnlhdlrSetFreeHdlrData(nlhdlr, nlhdlrFreehdlrdataQuadratic);
5060 	   SCIPnlhdlrSetFreeExprData(nlhdlr, nlhdlrFreeexprdataQuadratic);
5061 	   SCIPnlhdlrSetSepa(nlhdlr, NULL, nlhdlrEnfoQuadratic, NULL, NULL);
5062 	   SCIPnlhdlrSetProp(nlhdlr, nlhdlrIntevalQuadratic, nlhdlrReversepropQuadratic);
5063 	
5064 	   /* parameters */
5065 	   SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" NLHDLR_NAME "/useintersectioncuts",
5066 	         "whether to use intersection cuts for quadratic constraints to separate",
5067 	         &nlhdlrdata->useintersectioncuts, FALSE, DEFAULT_USEINTERCUTS, NULL, NULL) );
5068 	
5069 	   SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" NLHDLR_NAME "/usestrengthening",
5070 	         "whether the strengthening should be used",
5071 	         &nlhdlrdata->usestrengthening, FALSE, DEFAULT_USESTRENGTH, NULL, NULL) );
5072 	
5073 	   SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" NLHDLR_NAME "/usemonoidal",
5074 	         "whether monoidal strengthening should be used",
5075 	         &nlhdlrdata->usemonoidal, FALSE, DEFAULT_USEMONOIDAL, NULL, NULL) );
5076 	
5077 	   SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" NLHDLR_NAME "/useminrep",
5078 	         "whether the minimal representation of the S-free set should be used (instead of the gauge)",
5079 	         &nlhdlrdata->useminrep, FALSE, DEFAULT_USEMINREP, NULL, NULL) );
5080 	
5081 	   SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" NLHDLR_NAME "/useboundsasrays",
5082 	         "use bounds of variables in quadratic as rays for intersection cuts",
5083 	         &nlhdlrdata->useboundsasrays, FALSE, DEFAULT_USEBOUNDS, NULL, NULL) );
5084 	
5085 	   SCIP_CALL( SCIPaddIntParam(scip, "nlhdlr/" NLHDLR_NAME "/ncutslimit",
5086 	         "limit for number of cuts generated consecutively",
5087 	         &nlhdlrdata->ncutslimit, FALSE, DEFAULT_NCUTS, 0, INT_MAX, NULL, NULL) );
5088 	
5089 	   SCIP_CALL( SCIPaddIntParam(scip, "nlhdlr/" NLHDLR_NAME "/ncutslimitroot",
5090 	         "limit for number of cuts generated at root node",
5091 	         &nlhdlrdata->ncutslimitroot, FALSE, DEFAULT_NCUTSROOT, 0, INT_MAX, NULL, NULL) );
5092 	
5093 	   SCIP_CALL( SCIPaddIntParam(scip, "nlhdlr/" NLHDLR_NAME "/maxrank",
5094 	         "maximal rank a slackvar can have",
5095 	         &nlhdlrdata->maxrank, FALSE, INT_MAX, 0, INT_MAX, NULL, NULL) );
5096 	
5097 	   SCIP_CALL( SCIPaddRealParam(scip, "nlhdlr/" NLHDLR_NAME "/mincutviolation",
5098 	         "minimal cut violation the generated cuts must fulfill to be added to the LP",
5099 	         &nlhdlrdata->mincutviolation, FALSE, 1e-4, 0.0, SCIPinfinity(scip), NULL, NULL) );
5100 	
5101 	   SCIP_CALL( SCIPaddRealParam(scip, "nlhdlr/" NLHDLR_NAME "/minviolation",
5102 	         "minimal violation the constraint must fulfill such that a cut is generated",
5103 	         &nlhdlrdata->minviolation, FALSE, INTERCUTS_MINVIOL, 0.0, SCIPinfinity(scip), NULL, NULL) );
5104 	
5105 	   SCIP_CALL( SCIPaddIntParam(scip, "nlhdlr/" NLHDLR_NAME "/atwhichnodes",
5106 	         "determines at which nodes cut is used (if it's -1, it's used only at the root node, if it's n >= 0, it's used at every multiple of n",
5107 	         &nlhdlrdata->atwhichnodes, FALSE, 1, -1, INT_MAX, NULL, NULL) );
5108 	
5109 	   SCIP_CALL( SCIPaddIntParam(scip, "nlhdlr/" NLHDLR_NAME "/nstrengthlimit",
5110 	         "limit for number of rays we do the strengthening for",
5111 	         &nlhdlrdata->nstrengthlimit, FALSE, INT_MAX, 0, INT_MAX, NULL, NULL) );
5112 	
5113 	   SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" NLHDLR_NAME "/sparsifycuts",
5114 	         "should we try to sparisfy the intersection cut?",
5115 	         &nlhdlrdata->sparsifycuts, FALSE, FALSE, NULL, NULL) );
5116 	
5117 	   SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" NLHDLR_NAME "/ignorebadrayrestriction",
5118 	         "should cut be generated even with bad numerics when restricting to ray?",
5119 	         &nlhdlrdata->ignorebadrayrestriction, FALSE, TRUE, NULL, NULL) );
5120 	
5121 	   SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" NLHDLR_NAME "/ignorenhighre",
5122 	         "should cut be added even when range / efficacy is large?",
5123 	         &nlhdlrdata->ignorehighre, FALSE, TRUE, NULL, NULL) );
5124 	
5125 	   SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" NLHDLR_NAME "/trackmore",
5126 	         "for monoidal strengthening, should we track more statistics (more expensive)?",
5127 	         &nlhdlrdata->trackmore, FALSE, FALSE, NULL, NULL) );
5128 	
5129 	   /* statistic table */
5130 	   assert(SCIPfindTable(scip, TABLE_NAME_QUADRATIC) == NULL);
5131 	   SCIP_CALL( SCIPincludeTable(scip, TABLE_NAME_QUADRATIC, TABLE_DESC_QUADRATIC, TRUE,
5132 	         NULL, NULL, NULL, NULL, NULL, NULL, tableOutputQuadratic,
5133 	         NULL, TABLE_POSITION_QUADRATIC, TABLE_EARLIEST_STAGE_QUADRATIC) );
5134 	   return SCIP_OKAY;
5135 	}
5136