Scippy

SCIP

Solving Constraint Integer Programs

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