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