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  assert(coefs4b != NULL);
2040  for( j = 0; j < 3; ++j )
2041  {
2042  SCIP_Real absval;
2043 
2044  absval = ABS(coefs4b[j]);
2045  if( max < absval )
2046  max = absval;
2047  if( absval != 0.0 && absval < min )
2048  min = absval;
2049  }
2050 
2051  if( SCIPisHugeValue(scip, max / min) )
2052  {
2053  INTERLOG(printf("Bad numeric 4b: max(A,B,C)/min(A,B,C) is too large (%g)\n", max / min); )
2054  nlhdlrdata->nbadrayrestriction++;
2055  return FALSE;
2056  }
2057  }
2058  }
2059 
2060  return TRUE;
2061 }
2062 
2063 
2064 /** computes the coefficients a, b, c defining the quadratic function defining set S restricted to the line
2065  * theta * apex.
2066  *
2067  * The solution to the monoidal strengthening problem is then given by the smallest root of the function
2068  * a * theta^2 + b * theta + c
2069  */
2070 static
2072  SCIP* scip, /**< SCIP data structure */
2073  SCIP_NLHDLREXPRDATA* nlhdlrexprdata, /**< nlhdlr expression data */
2074  SCIP_Real* vb, /**< array containing \f$v_i^T b\f$ for \f$i \in I_+ \cup I_-\f$ */
2075  SCIP_Real* vzlp, /**< array containing \f$v_i^T zlp_q\f$ for \f$i \in I_+ \cup I_-\f$ */
2076  SCIP_Real* vapex, /**< array containing \f$v_i^T apex\f$ */
2077  SCIP_Real* vray, /**< array containing \f$v_i^T ray\f$ */
2078  SCIP_Real kappa, /**< value of kappa */
2079  SCIP_Real sidefactor, /**< 1.0 if the violated constraint is q &le; rhs, -1.0 otherwise */
2080  SCIP_Real* a, /**< pointer to store quadratic coefficient */
2081  SCIP_Real* b, /**< pointer to store linear coefficient */
2082  SCIP_Real* c /**< pointer to store constant */
2083  )
2084 {
2085  SCIP_EXPR* qexpr;
2086  int nquadexprs;
2087  SCIP_Real* eigenvectors;
2088  SCIP_Real* eigenvalues;
2089  int i;
2090 
2091  qexpr = nlhdlrexprdata->qexpr;
2092  SCIPexprGetQuadraticData(qexpr, NULL, NULL, NULL, NULL, &nquadexprs, NULL, &eigenvalues, &eigenvectors);
2093 
2094  *a = 0.0;
2095  *b = 0.0;
2096  *c = 0.0;
2097  for( i = 0; i < nquadexprs; ++i )
2098  {
2099  SCIP_Real dot;
2100 
2101  if( SCIPisZero(scip, sidefactor * eigenvalues[i]) )
2102  continue;
2103 
2104  dot = vray[i] + vapex[i] + vb[i] / (2.0 * sidefactor * eigenvalues[i]);
2105 
2106  *a += sidefactor * eigenvalues[i] * SQR(vzlp[i] - vapex[i]);
2107  *b += sidefactor * eigenvalues[i] * (vzlp[i] - vapex[i]) * dot;
2108  *c += sidefactor * eigenvalues[i] * SQR(dot);
2109  }
2110 
2111  *b *= 2.0;
2112  *c += kappa;
2113 
2114  return SCIP_OKAY;
2115 }
2116 
2117 /** check if ray was in strip by checking if the point in the monoid corresponding to the cutcoef we just found
2118  * is "on the wrong side" of the hyperplane -(a - lambda^Ta lambda)^T x
2119  */
2120 static
2122  SCIP_NLHDLREXPRDATA* nlhdlrexprdata, /**< nlhdlr expression data */
2123  SCIP_Real* vb, /**< array containing \f$v_i^T b\f$ for \f$i \in I_+ \cup I_-\f$ */
2124  SCIP_Real* vzlp, /**< array containing \f$v_i^T zlp_q\f$ for \f$i \in I_+ \cup I_-\f$ */
2125  SCIP_Real* vapex, /**< array containing \f$v_i^T apex\f$ */
2126  SCIP_Real* vray, /**< array containing \f$v_i^T ray\f$ */
2127  SCIP_Real kappa, /**< value of kappa */
2128  SCIP_Real sidefactor, /**< 1.0 if the violated constraint is q &le; rhs, -1.0 otherwise */
2129  SCIP_Real cutcoef /**< optimal solution of the monoidal quadratic */
2130  )
2131 {
2132  SCIP_EXPR* qexpr;
2133  int nquadexprs;
2134  SCIP_Real* eigenvectors;
2135  SCIP_Real* eigenvalues;
2136  SCIP_Real num;
2137  SCIP_Real denom;
2138  int i;
2139 
2140  qexpr = nlhdlrexprdata->qexpr;
2141  SCIPexprGetQuadraticData(qexpr, NULL, NULL, NULL, NULL, &nquadexprs, NULL, &eigenvalues, &eigenvectors);
2142 
2143  num = 0.0;
2144  denom = 0.0;
2145  for( i = 0; i < nquadexprs; ++i )
2146  {
2147  SCIP_Real dot;
2148 
2149  if( sidefactor * eigenvalues[i] <= 0.0 )
2150  continue;
2151 
2152  dot = vzlp[i] + vb[i] / (2.0 * sidefactor * eigenvalues[i]);
2153 
2154  num += sidefactor * eigenvalues[i] * dot * (cutcoef * (vzlp[i] - vapex[i]) + vray[i] + vapex[i]);
2155  denom += sidefactor * eigenvalues[i] * SQR(dot);
2156  }
2157 
2158  num /= kappa;
2159  num += 1.0;
2160  denom /= kappa;
2161  denom += 1.0;
2162 
2163  assert(denom > 0);
2164 
2165  return num / denom < 1;
2166 }
2167 
2168 /** computes the smallest root of the quadratic function a*x^2 + b*x + c with a > 0
2169  * and b^2 - 4ac >= 0. We use binary search between -inf and minimum at -b/2a.
2170  */
2171 static
2173  SCIP* scip, /**< SCIP data structure */
2174  SCIP_Real a, /**< quadratic coefficient */
2175  SCIP_Real b, /**< linear coefficient */
2176  SCIP_Real c /**< constant */
2177  )
2178 {
2179  SCIP_Real sol;
2180  SCIP_INTERVAL bounds;
2181  SCIP_INTERVAL result;
2182 
2183  assert(SQR(b) - 4 * a * c >= 0.0);
2184 
2185  if( SCIPisZero(scip, a) )
2186  {
2187  assert(b != 0.0);
2188  sol = - c / b;
2189  }
2190  else
2191  {
2192  SCIPintervalSetBounds(&bounds, - b / (2.0 * a), SCIPinfinity(scip));
2193 
2194  /* find all positive x such that a x^2 + b x >= -c and x in bounds.*/
2197 
2198  /* if we didn't find any positive solutions, negate quadratic and find negative solutions */
2199  if( SCIPisInfinity(scip, sol) )
2200  {
2201  SCIPintervalSetBounds(&bounds, b / (2.0 * a), SCIPinfinity(scip));
2202 
2203  /* find all positive x such that a x^2 - b x >= -c and x in bounds.*/
2206  }
2207  }
2208 
2209  /* check if that solution is close enough or if we need to improve it more with binary search */
2210  if( a * SQR(sol) + sol * b + c > 1e-10 )
2211  {
2212  SCIP_Real val;
2213  SCIP_Real lb;
2214  SCIP_Real ub;
2215  SCIP_Real lastposval;
2216  SCIP_Real lastpossol;
2217  int niter;
2218 
2219  lb = - b / (2.0 * a);
2220  ub = sol;
2221  lastposval = SCIPinfinity(scip);
2222  lastpossol = SCIPinfinity(scip);
2223  val = SCIPinfinity(scip);
2224  niter = 0;
2225  while( niter < BINSEARCH_MAXITERS && ABS(val) > 1e-10 )
2226  {
2227  sol = (ub + lb) / 2.0;
2228  val = a * SQR(sol) + b * sol + c;
2229 
2230  if( val < 0.0 )
2231  lb = sol;
2232  else
2233  ub = sol;
2234 
2235  /* if we are close enough, return with (feasible) solution */
2236  if( val > 0.0 && val < 1e-6 )
2237  break;
2238 
2239  if( val > 0.0 && lastposval > val )
2240  {
2241  lastposval = val;
2242  lastpossol = sol;
2243  }
2244 
2245  ++niter;
2246  }
2247  if( val < 0.0 && ! SCIPisZero(scip, val) )
2248  {
2249  sol = lastpossol;
2250  val = lastposval;
2251  }
2252 
2253  assert( val > 0.0 || SCIPisZero(scip, val) );
2254  }
2255 
2256  return sol;
2257 }
2258 
2259 /** computes the apex of the S-free set (if it exists) */
2260 static
2262  SCIP_NLHDLREXPRDATA* nlhdlrexprdata, /**< nlhdlr expression data */
2263  SCIP_Real* vb, /**< array containing \f$v_i^T b\f$ for \f$i \in I_+ \cup I_-\f$ */
2264  SCIP_Real* vzlp, /**< array containing \f$v_i^T zlp_q\f$ for \f$i \in I_+ \cup I_-\f$ */
2265  SCIP_Real kappa, /**< value of kappa */
2266  SCIP_Real sidefactor, /**< 1.0 if the violated constraint is q &le; rhs, -1.0 otherwise */
2267  SCIP_Real* apex, /**< array to store apex */
2268  SCIP_Bool* success /**< TRUE if computation of apex was successful */
2269  )
2270 {
2271  SCIP_EXPR* qexpr;
2272  int nquadexprs;
2273  SCIP_Real* eigenvectors;
2274  SCIP_Real* eigenvalues;
2275  int i;
2276 
2277  qexpr = nlhdlrexprdata->qexpr;
2278  SCIPexprGetQuadraticData(qexpr, NULL, NULL, NULL, NULL, &nquadexprs, NULL, &eigenvalues, &eigenvectors);
2279 
2280  *success = TRUE;
2281 
2282  for( i = 0; i < nquadexprs; ++i )
2283  {
2284  SCIP_Real entry;
2285  SCIP_Real denom;
2286  SCIP_Real num;
2287  int j;
2288 
2289  entry = 0.0;
2290  num = 0.0;
2291  denom = 0.0;
2292  for( j = 0; j < nquadexprs; ++j )
2293  {
2294  if( sidefactor * eigenvalues[j] == 0.0 )
2295  continue;
2296 
2297  entry -= eigenvectors[j * nquadexprs + i] * vb[j] / (2.0 * sidefactor * eigenvalues[j]);
2298 
2299  if( sidefactor * eigenvalues[j] > 0.0 )
2300  {
2301  SCIP_Real dot;
2302 
2303  dot = vzlp[j] + vb[j] / (2.0 * sidefactor * eigenvalues[j]);
2304 
2305  num += eigenvectors[j * nquadexprs + i] * dot;
2306  denom += sidefactor * eigenvalues[j] * SQR(dot);
2307  }
2308  }
2309 
2310  /* if denom = 0, the S-free set is just the strip, so we don't have an apex -> monoidal not possible */
2311  if( denom == 0.0 )
2312  {
2313  *success = FALSE;
2314  return;
2315  }
2316  assert(denom > 0.0);
2317 
2318  num *= -kappa;
2319  entry += num / denom;
2320 
2321  apex[i] = entry;
2322  }
2323 }
2324 
2325 /** for a given ray, computes the cut coefficient using monoidal strengthening (if possible) */
2326 static
2328  SCIP* scip, /**< SCIP data structure */
2329  SCIP_NLHDLREXPRDATA* nlhdlrexprdata, /**< nlhdlr expression data */
2330  int lppos, /**< lp pos of current ray */
2331  SCIP_Real* raycoefs, /**< coefficients of ray */
2332  int* rayidx, /**< index of consvar the ray coef is associated to */
2333  int raynnonz, /**< length of raycoefs and rayidx */
2334  SCIP_Real* vb, /**< array containing \f$v_i^T b\f$ for \f$i \in I_+ \cup I_-\f$ */
2335  SCIP_Real* vzlp, /**< array containing \f$v_i^T zlp_q\f$ for \f$i \in I_+ \cup I_-\f$ */
2336  SCIP_Real* wcoefs, /**< coefficients of w for the qud vars or NULL if w is 0 */
2337  SCIP_Real kappa, /**< value of kappa */
2338  SCIP_Real* apex, /**< array containing the apex of the S-free set in the original space */
2339  SCIP_Real sidefactor, /**< 1.0 if the violated constraint is q &le; rhs, -1.0 otherwise */
2340  SCIP_Real* cutcoef, /**< pointer to store cut coef */
2341  SCIP_Bool* success /**< TRUE if monoidal strengthening could be applied */
2342  )
2343 {
2344  SCIP_COL** cols;
2345  SCIP_ROW** rows;
2346 
2347  *success = FALSE;
2348 
2349  /* check if we are in the correct case (case 2) */
2350  assert(wcoefs == NULL && kappa > 0.0);
2351 
2352  cols = SCIPgetLPCols(scip);
2353  rows = SCIPgetLPRows(scip);
2354 
2355  /* if var corresponding to current ray is integer, we can do monoidal */
2356  if( (lppos >= 0 && SCIPvarGetType(SCIPcolGetVar(cols[lppos])) != SCIP_VARTYPE_CONTINUOUS) ||
2357  (lppos < 0 && SCIProwIsIntegral(rows[- lppos - 1])) )
2358  {
2359  SCIP_Real* vapex;
2360  SCIP_Real* vray;
2361  SCIP_Real a;
2362  SCIP_Real b;
2363  SCIP_Real c;
2364  int nquadexprs;
2365 
2366  SCIPexprGetQuadraticData(nlhdlrexprdata->qexpr, NULL, NULL, NULL, NULL, &nquadexprs, NULL, NULL, NULL);
2367 
2368  /* compute v_i^T apex in vapex[i] and v_i^T ray in vray[i] */
2369  SCIP_CALL( SCIPallocBufferArray(scip, &vapex, nquadexprs) );
2370  SCIP_CALL( SCIPallocBufferArray(scip, &vray, nquadexprs) );
2371  computeVApexAndVRay(nlhdlrexprdata, apex, raycoefs, rayidx, raynnonz, vapex, vray);
2372 
2373  /* compute coefficients of the quadratic monoidal problem function */
2374  SCIP_CALL( computeMonoidalQuadCoefs(scip, nlhdlrexprdata, vb, vzlp, vapex, vray, kappa,
2375  sidefactor, &a, &b, &c) );
2376 
2377  /* check if ray is in strip */
2378  if( SQR(b) - (4 * a * c) >= 0.0 )
2379  {
2380  /* find smallest root of quadratic function a * x^2 + b * x + c -> this is the cut coef */
2381  *cutcoef = findMonoidalQuadRoot(scip, a, b, c);
2382 
2383  /* if the cutcoef is negative, we have to set it to 0 */
2384  *cutcoef = MAX(*cutcoef, 0.0);
2385 
2386  /* check if ray is in strip. If not, monoidal is possible and cutcoef is the strengthened cut coef */
2387  if( ! isRayInStrip(nlhdlrexprdata, vb, vzlp, vapex, vray, kappa, sidefactor, *cutcoef) )
2388  {
2389  *success = TRUE;
2390  }
2391  }
2392 
2393  SCIPfreeBufferArray(scip, &vray);
2394  SCIPfreeBufferArray(scip, &vapex);
2395  }
2396 
2397  return SCIP_OKAY;
2398 }
2399 
2400 /** sparsify intersection cut by replacing non-basic variables with their bounds if their coefficient allows it */
2401 static
2403  SCIP* scip, /**< SCIP data structure */
2404  SCIP_ROWPREP* rowprep /**< rowprep for the generated cut */
2405  )
2406 {
2407  int i;
2408  int nvars;
2409 
2410  /* get number of variables in rowprep */
2411  nvars = SCIProwprepGetNVars(rowprep);
2412 
2413  /* go though all the variables in rowprep */
2414  for( i = 0; i < nvars; ++i )
2415  {
2416  SCIP_VAR* var;
2417  SCIP_Real coef;
2418  SCIP_Real lb;
2419  SCIP_Real ub;
2420  SCIP_Real solval;
2421 
2422  /* get variable and its coefficient */
2423  var = SCIProwprepGetVars(rowprep)[i];
2424  coef = SCIProwprepGetCoefs(rowprep)[i];
2425 
2426  /* get bounds of variable */
2427  lb = SCIPvarGetLbLocal(var);
2428  ub = SCIPvarGetUbLocal(var);
2429 
2430  /* get LP solution value of variable */
2431  solval = SCIPgetSolVal(scip, NULL, var);
2432 
2433  /* if the variable is at its lower or upper bound and the coefficient has the correct sign, we can
2434  * set the cutcoef to 0
2435  */
2436  if( SCIPisZero(scip, ub - solval) && coef > 0.0 )
2437  {
2438  SCIProwprepAddSide(rowprep, -coef * ub);
2439  SCIProwprepSetCoef(rowprep, i, 0.0);
2440  }
2441  else if( SCIPisZero(scip, solval - lb) && coef < 0.0 )
2442  {
2443  SCIProwprepAddSide(rowprep, -coef * lb);
2444  SCIProwprepSetCoef(rowprep, i, 0.0);
2445  }
2446  }
2447 
2448  return;
2449 }
2450 
2451 /** computes intersection cut cuts off sol (because solution sol violates the quadratic constraint cons)
2452  * and stores it in rowprep. Here, we don't use any strengthening.
2453  */
2454 static
2456  SCIP* scip, /**< SCIP data structure */
2457  SCIP_NLHDLRDATA* nlhdlrdata, /**< nlhdlr data */
2458  SCIP_NLHDLREXPRDATA* nlhdlrexprdata, /**< nlhdlr expression data */
2459  RAYS* rays, /**< rays */
2460  SCIP_Real sidefactor, /**< 1.0 if the violated constraint is q &le; rhs, -1.0 otherwise */
2461  SCIP_Bool iscase4, /**< whether we are in case 4 */
2462  SCIP_Real* vb, /**< array containing \f$v_i^T b\f$ for \f$i \in I_+ \cup I_-\f$ */
2463  SCIP_Real* vzlp, /**< array containing \f$v_i^T zlp_q\f$ for \f$i \in I_+ \cup I_-\f$ */
2464  SCIP_Real* wcoefs, /**< coefficients of w for the qud vars or NULL if w is 0 */
2465  SCIP_Real wzlp, /**< value of w at zlp */
2466  SCIP_Real kappa, /**< value of kappa */
2467  SCIP_ROWPREP* rowprep, /**< rowprep for the generated cut */
2468  SCIP_Real* interpoints, /**< array to store intersection points for all rays or NULL if nothing
2469  * needs to be stored */
2470  SCIP_SOL* sol, /**< solution we want to separate */
2471  SCIP_Bool* success /**< if a cut candidate could be computed */
2472  )
2473 {
2474  SCIP_COL** cols;
2475  SCIP_ROW** rows;
2476  SCIP_Real* apex;
2477  SCIP_Real avecutcoefsum;
2478  SCIP_Real avemonoidalimprovsum;
2479  int monoidalcounter;
2480  int counter;
2481  int i;
2482  SCIP_Bool usemonoidal;
2483  SCIP_Bool monoidalwasused;
2484  SCIP_Bool case2;
2485 
2486  cols = SCIPgetLPCols(scip);
2487  rows = SCIPgetLPRows(scip);
2488 
2489  case2 = wcoefs == NULL && kappa > 0.0;
2490  monoidalwasused = FALSE;
2491 
2492  counter = 0;
2493  monoidalcounter = 0;
2494  avecutcoefsum = 0.0;
2495  avemonoidalimprovsum = 0.0;
2496 
2497  /* if we use monoidal and we are in the right case for it, compute the apex of the S-free set */
2498  if( case2 )
2499  {
2500  int nquadexprs;
2501 
2502  SCIPexprGetQuadraticData(nlhdlrexprdata->qexpr, NULL, NULL, NULL, NULL, &nquadexprs, NULL, NULL, NULL);
2503 
2504  /* allocate memory for apex */
2505  SCIP_CALL( SCIPallocBufferArray(scip, &apex, nquadexprs) );
2506 
2507  computeApex(nlhdlrexprdata, vb, vzlp, kappa, sidefactor, apex, success);
2508 
2509  /* if computation of apex was not successful, don't apply monoidal strengthening */
2510  if( ! *success )
2511  case2 = FALSE;
2512  }
2513  else
2514  {
2515  apex = NULL;
2516  }
2517 
2518  usemonoidal = nlhdlrdata->usemonoidal && case2;
2519 
2520  /* for every ray: compute cut coefficient and add var associated to ray into cut */
2521  for( i = 0; i < rays->nrays; ++i )
2522  {
2523  SCIP_Real interpoint;
2524  SCIP_Real cutcoef;
2525  int lppos;
2526  SCIP_Real coefs1234a[5];
2527  SCIP_Real coefs4b[5];
2528  SCIP_Real coefscondition[3];
2529  SCIP_Bool monoidalsuccess;
2530 
2531  monoidalsuccess = FALSE;
2532  cutcoef = SCIPinfinity(scip);
2533 
2534  /* if we (can) use monoidal strengthening, compute the cut coefficient with that */
2535  if( usemonoidal )
2536  {
2537  SCIP_CALL( computeMonoidalStrengthCoef(scip, nlhdlrexprdata, rays->lpposray[i], &rays->rays[rays->raysbegin[i]],
2538  &rays->raysidx[rays->raysbegin[i]], rays->raysbegin[i + 1] - rays->raysbegin[i], vb, vzlp, wcoefs, kappa,
2539  apex, sidefactor, &cutcoef, &monoidalsuccess) );
2540  }
2541 
2542  /* track whether monoidal was successful at least once if it is used */
2543  if( usemonoidal && ! monoidalwasused && monoidalsuccess )
2544  monoidalwasused = TRUE;
2545 
2546  /* if we don't use monoidal or if monoidal couldn't be applied, use gauge to compute coef */
2547  if( ! usemonoidal || ! monoidalsuccess || nlhdlrdata->trackmore )
2548  {
2549  /* restrict phi to ray */
2550  SCIP_CALL( computeRestrictionToRay(scip, nlhdlrexprdata, sidefactor, iscase4,
2551  &rays->rays[rays->raysbegin[i]], &rays->raysidx[rays->raysbegin[i]], rays->raysbegin[i + 1] -
2552  rays->raysbegin[i], vb, vzlp, wcoefs, wzlp, kappa, coefs1234a, coefs4b, coefscondition, success) );
2553 
2554  if( ! *success )
2555  goto TERMINATE;
2556 
2557  /* if restriction to ray is numerically nasty -> abort cut separation */
2558  *success = areCoefsNumericsGood(scip, nlhdlrdata, coefs1234a, coefs4b, iscase4);
2559 
2560  if( ! *success )
2561  goto TERMINATE;
2562 
2563  /* compute intersection point */
2564  interpoint = computeIntersectionPoint(scip, nlhdlrdata, iscase4, coefs1234a, coefs4b, coefscondition);
2565 
2566 #ifdef DEBUG_INTERSECTIONCUT
2567  SCIPinfoMessage(scip, NULL, "interpoint for ray %d is %g\n", i, interpoint);
2568 #endif
2569 
2570  /* store intersection point */
2571  if( interpoints != NULL )
2572  interpoints[i] = interpoint;
2573 
2574  /* if we are only computing the "normal" cut coefficient to track statistics (and we have been successful
2575  * with monoidal strengthening), don't overwrite cutcoef variable, but just track statistics */
2576  if( nlhdlrdata->trackmore && monoidalsuccess )
2577  {
2578  SCIP_Real normalcutcoef;
2579 
2580  normalcutcoef = SCIPisInfinity(scip, interpoint) ? 0.0 : 1.0 / interpoint;
2581 
2582  if( cutcoef >= 0 )
2583  avemonoidalimprovsum += cutcoef / normalcutcoef;
2584  ++monoidalcounter;
2585  }
2586  else
2587  {
2588  /* compute cut coef */
2589  cutcoef = SCIPisInfinity(scip, interpoint) ? 0.0 : 1.0 / interpoint;
2590  }
2591 
2592  if( cutcoef == 0.0 && case2 && nlhdlrdata->useminrep )
2593  {
2594  /* restrict phi to the line defined by ray + apex + t*(f - apex) */
2595  SCIP_CALL( computeRestrictionToLine(scip, nlhdlrexprdata, sidefactor,
2596  &rays->rays[rays->raysbegin[i]], &rays->raysidx[rays->raysbegin[i]], rays->raysbegin[i + 1] -
2597  rays->raysbegin[i], vb, vzlp, kappa, apex, coefs1234a, success) );
2598 
2599  if( ! *success )
2600  goto TERMINATE;
2601 
2602  /* if restriction to ray is numerically nasty -> abort cut separation */
2603  *success = areCoefsNumericsGood(scip, nlhdlrdata, coefs1234a, NULL, FALSE);
2604 
2605  if( ! *success )
2606  goto TERMINATE;
2607 
2608  coefs1234a[1] *= -1.0;
2609  coefs1234a[3] *= -1.0;
2610 
2611  /* compute intersection point */
2612  cutcoef = - computeRoot(scip, coefs1234a);
2613 
2614  assert(cutcoef <= 0.0);
2615  }
2616  }
2617 
2618  /* keep track of average cut coefficient */
2619  ++counter;
2620  avecutcoefsum += cutcoef;
2621  assert( ! SCIPisInfinity(scip, cutcoef) );
2622 
2623  /* add var to cut: if variable is nonbasic at upper we have to flip sign of cutcoef */
2624  lppos = rays->lpposray[i];
2625  if( lppos < 0 )
2626  {
2627  lppos = -lppos - 1;
2628 
2629  assert(SCIProwGetBasisStatus(rows[lppos]) == SCIP_BASESTAT_LOWER || SCIProwGetBasisStatus(rows[lppos]) ==
2631 
2632  /* flip cutcoef when necessary. Note: rows have flipped base status! */
2633  cutcoef = SCIProwGetBasisStatus(rows[lppos]) == SCIP_BASESTAT_UPPER ? cutcoef : -cutcoef;
2634 
2635  SCIP_CALL( addRowToCut(scip, rowprep, cutcoef, rows[lppos], success) );
2636 
2637  if( ! *success )
2638  {
2639  INTERLOG(printf("Bad numeric: now not nonbasic enough\n");)
2640  nlhdlrdata->nbadnonbasic++;
2641  goto TERMINATE;
2642  }
2643  }
2644  else
2645  {
2646  if( ! nlhdlrdata->useboundsasrays )
2647  {
2648  assert(SCIPcolGetBasisStatus(cols[lppos]) == SCIP_BASESTAT_UPPER || SCIPcolGetBasisStatus(cols[lppos]) ==
2650 
2651  /* flip cutcoef when necessary */
2652  cutcoef = SCIPcolGetBasisStatus(cols[lppos]) == SCIP_BASESTAT_UPPER ? -cutcoef : cutcoef;
2653 
2654  SCIP_CALL( addColToCut(scip, rowprep, sol, cutcoef, cols[lppos]) );
2655  }
2656  else
2657  {
2658  /* flip cutcoef when necessary */
2659  cutcoef = rays->rays[i] == -1 ? -cutcoef : cutcoef;
2660 
2661  SCIP_CALL( addColToCut(scip, rowprep, sol, cutcoef, cols[lppos]) );
2662  }
2663  }
2664  }
2665 
2666  /* track statistics */
2667  if( counter > 0 )
2668  nlhdlrdata->currentavecutcoef = avecutcoefsum / counter;
2669  if( monoidalwasused )
2670  nlhdlrdata->nmonoidal += 1;
2671  if( monoidalcounter > 0 )
2672  nlhdlrdata->currentavemonoidalimprovement = avemonoidalimprovsum / monoidalcounter;
2673 
2674 TERMINATE:
2675  SCIPfreeBufferArrayNull(scip, &apex);
2676 
2677  return SCIP_OKAY;
2678 }
2679 
2680 /** combine ray 1 and 2 to obtain new ray coef1 * ray1 - coef2 * ray2 in sparse format */
2681 static
2683  SCIP_Real* raycoefs1, /**< coefficients of ray 1 */
2684  int* rayidx1, /**< index of consvar of the ray 1 coef is associated to */
2685  int raynnonz1, /**< length of raycoefs1 and rayidx1 */
2686  SCIP_Real* raycoefs2, /**< coefficients of ray 2 */
2687  int* rayidx2, /**< index of consvar of the ray 2 coef is associated to */
2688  int raynnonz2, /**< length of raycoefs2 and rayidx2 */
2689  SCIP_Real* newraycoefs, /**< coefficients of combined ray */
2690  int* newrayidx, /**< index of consvar of the combined ray coef is associated to */
2691  int* newraynnonz, /**< pointer to length of newraycoefs and newrayidx */
2692  SCIP_Real coef1, /**< coef of ray 1 */
2693  SCIP_Real coef2 /**< coef of ray 2 */
2694  )
2695 {
2696  int idx1;
2697  int idx2;
2698 
2699  idx1 = 0;
2700  idx2 = 0;
2701  *newraynnonz = 0;
2702 
2703  while( idx1 < raynnonz1 || idx2 < raynnonz2 )
2704  {
2705  /* if the pointers look at different variables (or one already arrieved at the end), only one pointer can move
2706  * on
2707  */
2708  if( idx1 >= raynnonz1 || (idx2 < raynnonz2 && rayidx1[idx1] > rayidx2[idx2]) )
2709  {
2710  /*printf("case 1 \n"); */
2711  newraycoefs[*newraynnonz] = - coef2 * raycoefs2[idx2];
2712  newrayidx[*newraynnonz] = rayidx2[idx2];
2713  ++(*newraynnonz);
2714  ++idx2;
2715  }
2716  else if( idx2 >= raynnonz2 || rayidx1[idx1] < rayidx2[idx2] )
2717  {
2718  /*printf("case 2 \n"); */
2719  newraycoefs[*newraynnonz] = coef1 * raycoefs1[idx1];
2720  newrayidx[*newraynnonz] = rayidx1[idx1];
2721  ++(*newraynnonz);
2722  ++idx1;
2723  }
2724  /* if both pointers look at the same variable, just compute the difference and move both pointers */
2725  else if( rayidx1[idx1] == rayidx2[idx2] )
2726  {
2727  /*printf("case 3 \n"); */
2728  newraycoefs[*newraynnonz] = coef1 * raycoefs1[idx1] - coef2 * raycoefs2[idx2];
2729  newrayidx[*newraynnonz] = rayidx1[idx1];
2730  ++(*newraynnonz);
2731  ++idx1;
2732  ++idx2;
2733  }
2734  }
2735 }
2736 
2737 /** checks if two rays are linearly dependent */
2738 static
2740  SCIP* scip, /**< SCIP data structure */
2741  SCIP_Real* raycoefs1, /**< coefficients of ray 1 */
2742  int* rayidx1, /**< index of consvar of the ray 1 coef is associated to */
2743  int raynnonz1, /**< length of raycoefs1 and rayidx1 */
2744  SCIP_Real* raycoefs2, /**< coefficients of ray 2 */
2745  int* rayidx2, /**< index of consvar of the ray 2 coef is associated to */
2746  int raynnonz2, /**< length of raycoefs2 and rayidx2 */
2747  SCIP_Real* coef /**< pointer to store coef (s.t. r1 = coef * r2) in case rays are
2748  * dependent */
2749  )
2750 {
2751  int i;
2752 
2753  /* cannot be dependent if they have different number of non-zero entries */
2754  if( raynnonz1 != raynnonz2 )
2755  return FALSE;
2756 
2757  *coef = 0.0;
2758 
2759  for( i = 0; i < raynnonz1; ++i )
2760  {
2761  /* cannot be dependent if different variables have non-zero entries */
2762  if( rayidx1[i] != rayidx2[i] ||
2763  (SCIPisZero(scip, raycoefs1[i]) && !SCIPisZero(scip, raycoefs2[i])) ||
2764  (!SCIPisZero(scip, raycoefs1[i]) && SCIPisZero(scip, raycoefs2[i])) )
2765  {
2766  return FALSE;
2767  }
2768 
2769  if( *coef != 0.0 )
2770  {
2771  /* cannot be dependent if the coefs aren't equal for all entries */
2772  if( ! SCIPisEQ(scip, *coef, raycoefs1[i] / raycoefs2[i]) )
2773  {
2774  return FALSE;
2775  }
2776  }
2777  else
2778  *coef = raycoefs1[i] / raycoefs2[i];
2779  }
2780 
2781  return TRUE;
2782 }
2783 
2784 /** checks if the ray alpha * ray_i + (1 - alpha) * ray_j is in the recession cone of the S-free set. To do so,
2785  * we check if phi restricted to the ray has a positive root.
2786  */
2787 static
2789  SCIP* scip, /**< SCIP data structure */
2790  SCIP_NLHDLRDATA* nlhdlrdata, /**< nlhdlr data */
2791  SCIP_NLHDLREXPRDATA* nlhdlrexprdata, /**< nlhdlr expression data */
2792  RAYS* rays, /**< rays */
2793  int j, /**< index of current ray in recession cone */
2794  int i, /**< index of current ray not in recession cone */
2795  SCIP_Real sidefactor, /**< 1.0 if the violated constraint is q &le; rhs, -1.0 otherwise */
2796  SCIP_Bool iscase4, /**< whether we are in case 4 */
2797  SCIP_Real* vb, /**< array containing \f$v_i^T b\f$ for \f$i \in I_+ \cup I_-\f$ */
2798  SCIP_Real* vzlp, /**< array containing \f$v_i^T zlp_q\f$ for \f$i \in I_+ \cup I_-\f$ */
2799  SCIP_Real* wcoefs, /**< coefficients of w for the quad vars or NULL if w is 0 */
2800  SCIP_Real wzlp, /**< value of w at zlp */
2801  SCIP_Real kappa, /**< value of kappa */
2802  SCIP_Real alpha, /**< coef for combining the two rays */
2803  SCIP_Bool* inreccone, /**< pointer to store whether the ray is in the recession cone or not */
2804  SCIP_Bool* success /**< Did numerical troubles occur? */
2805  )
2806 {
2807  SCIP_Real coefs1234a[5];
2808  SCIP_Real coefs4b[5];
2809  SCIP_Real coefscondition[3];
2810  SCIP_Real interpoint;
2811  SCIP_Real* newraycoefs;
2812  int* newrayidx;
2813  int newraynnonz;
2814 
2815  *inreccone = FALSE;
2816 
2817  /* allocate memory for new ray */
2818  newraynnonz = (rays->raysbegin[i + 1] - rays->raysbegin[i]) + (rays->raysbegin[j + 1] - rays->raysbegin[j]);
2819  SCIP_CALL( SCIPallocBufferArray(scip, &newraycoefs, newraynnonz) );
2820  SCIP_CALL( SCIPallocBufferArray(scip, &newrayidx, newraynnonz) );
2821 
2822  /* build the ray alpha * ray_i + (1 - alpha) * ray_j */
2823  combineRays(&rays->rays[rays->raysbegin[i]], &rays->raysidx[rays->raysbegin[i]], rays->raysbegin[i + 1] -
2824  rays->raysbegin[i], &rays->rays[rays->raysbegin[j]], &rays->raysidx[rays->raysbegin[j]],
2825  rays->raysbegin[j + 1] - rays->raysbegin[j], newraycoefs, newrayidx, &newraynnonz, alpha,
2826  -1 + alpha);
2827 
2828  /* restrict phi to the "new" ray */
2829  SCIP_CALL( computeRestrictionToRay(scip, nlhdlrexprdata, sidefactor, iscase4, newraycoefs, newrayidx,
2830  newraynnonz, vb, vzlp, wcoefs, wzlp, kappa, coefs1234a, coefs4b, coefscondition, success) );
2831 
2832  if( ! *success )
2833  goto CLEANUP;
2834 
2835  /* check if restriction to "new" ray is numerically nasty. If so, treat the corresponding rho as if phi is
2836  * positive
2837  */
2838 
2839  /* compute intersection point */
2840  interpoint = computeIntersectionPoint(scip, nlhdlrdata, iscase4, coefs1234a, coefs4b, coefscondition);
2841 
2842  /* no root exists */
2843  if( SCIPisInfinity(scip, interpoint) )
2844  *inreccone = TRUE;
2845 
2846 CLEANUP:
2847  SCIPfreeBufferArray(scip, &newrayidx);
2848  SCIPfreeBufferArray(scip, &newraycoefs);
2849 
2850  return SCIP_OKAY;
2851 }
2852 
2853 /** finds the smallest negative steplength for the current ray r_idx such that the combination
2854  * of r_idx with all rays not in the recession cone is in the recession cone
2855  */
2856 static
2858  SCIP* scip, /**< SCIP data structure */
2859  SCIP_NLHDLRDATA* nlhdlrdata, /**< nlhdlr data */
2860  SCIP_NLHDLREXPRDATA* nlhdlrexprdata, /**< nlhdlr expression data */
2861  RAYS* rays, /**< rays */
2862  int idx, /**< index of current ray we want to find rho for */
2863  SCIP_Real sidefactor, /**< 1.0 if the violated constraint is q &le; rhs, -1.0 otherwise */
2864  SCIP_Bool iscase4, /**< whether we are in case 4 */
2865  SCIP_Real* vb, /**< array containing \f$v_i^T b\f$ for \f$i \in I_+ \cup I_-\f$ */
2866  SCIP_Real* vzlp, /**< array containing \f$v_i^T zlp_q\f$ for \f$i \in I_+ \cup I_-\f$ */
2867  SCIP_Real* wcoefs, /**< coefficients of w for the quad vars or NULL if w is 0 */
2868  SCIP_Real wzlp, /**< value of w at zlp */
2869  SCIP_Real kappa, /**< value of kappa */
2870  SCIP_Real* interpoints, /**< array to store intersection points for all rays or NULL if nothing
2871  * needs to be stored */
2872  SCIP_Real* rho, /**< pointer to store the optimal rho */
2873  SCIP_Bool* success /**< could we successfully find the right rho? */
2874  )
2875 {
2876  int i;
2877 
2878  /* go through all rays not in the recession cone and compute the largest negative steplength possible. The
2879  * smallest of them is then the steplength rho we use for the current ray */
2880  *rho = 0.0;
2881  for( i = 0; i < rays->nrays; ++i )
2882  {
2883  SCIP_Real currentrho;
2884  SCIP_Real coef;
2885 
2886  if( SCIPisInfinity(scip, interpoints[i]) )
2887  continue;
2888 
2889  /* if we cannot strengthen enough, we don't strengthen at all */
2890  if( SCIPisInfinity(scip, -*rho) )
2891  {
2892  *rho = -SCIPinfinity(scip);
2893  return SCIP_OKAY;
2894  }
2895 
2896  /* if the rays are linearly independent, we don't need to search for rho */
2897  if( raysAreDependent(scip, &rays->rays[rays->raysbegin[i]], &rays->raysidx[rays->raysbegin[i]],
2898  rays->raysbegin[i + 1] - rays->raysbegin[i], &rays->rays[rays->raysbegin[idx]],
2899  &rays->raysidx[rays->raysbegin[idx]], rays->raysbegin[idx + 1] - rays->raysbegin[idx], &coef) )
2900  {
2901  currentrho = coef * interpoints[i];
2902  }
2903  else
2904  {
2905  /* since the two rays are linearly independent, we need to find the biggest alpha such that
2906  * alpha * ray_i + (1 - alpha) * ray_idx in the recession cone is. For every ray i, we compute
2907  * such a alpha but take the smallest one of them. We use "maxalpha" to keep track of this.
2908  * Since we know that we can only use alpha < maxalpha, we don't need to do the whole binary search
2909  * for every ray i. We only need to search the intervall [0, maxalpha]. Thereby, we start by checking
2910  * if alpha = maxalpha is already feasable */
2911 
2912  SCIP_Bool inreccone;
2913  SCIP_Real alpha;
2914  SCIP_Real lb;
2915  SCIP_Real ub;
2916  int j;
2917 
2918  lb = 0.0;
2919  ub = 1.0;
2920  for( j = 0; j < BINSEARCH_MAXITERS; ++j )
2921  {
2922  alpha = (lb + ub) / 2.0;
2923 
2924  if( SCIPisZero(scip, alpha) )
2925  {
2926  alpha = 0.0;
2927  break;
2928  }
2929 
2930  SCIP_CALL( rayInRecessionCone(scip, nlhdlrdata, nlhdlrexprdata, rays, idx, i, sidefactor, iscase4, vb,
2931  vzlp, wcoefs, wzlp, kappa, alpha, &inreccone, success) );
2932 
2933  if( ! *success )
2934  return SCIP_OKAY;
2935 
2936  /* no root exists */
2937  if( inreccone )
2938  {
2939  lb = alpha;
2940  if( SCIPisEQ(scip, ub, lb) )
2941  break;
2942  }
2943  else
2944  ub = alpha;
2945  }
2946 
2947  /* now we found the best convex combination which we use to derive the corresponding coef. If alpha = 0, we
2948  * cannot move the ray in the recession cone, i.e. strengthening is not possible */
2949  if( SCIPisZero(scip, alpha) )
2950  {
2951  *rho = -SCIPinfinity(scip);
2952  return SCIP_OKAY;
2953  }
2954  else
2955  currentrho = (alpha - 1) * interpoints[i] / alpha; /*lint !e795*/
2956  }
2957 
2958  if( currentrho < *rho )
2959  *rho = currentrho;
2960 
2961  if( *rho < -10e+06 )
2962  *rho = -SCIPinfinity(scip);
2963 
2964  /* if rho is too small, don't add it */
2965  if( SCIPisZero(scip, *rho) )
2966  *success = FALSE;
2967  }
2968 
2969  return SCIP_OKAY;
2970 }
2971 
2972 /** computes intersection cut using negative edge extension to strengthen rays that do not intersect
2973  * (i.e., rays in the recession cone)
2974  */
2975 static
2977  SCIP* scip, /**< SCIP data structure */
2978  SCIP_NLHDLRDATA* nlhdlrdata, /**< nlhdlr data */
2979  SCIP_NLHDLREXPRDATA* nlhdlrexprdata, /**< nlhdlr expression data */
2980  RAYS* rays, /**< rays */
2981  SCIP_Real sidefactor, /**< 1.0 if the violated constraint is q &le; rhs, -1.0 otherwise */
2982  SCIP_Bool iscase4, /**< whether we are in case 4 */
2983  SCIP_Real* vb, /**< array containing \f$v_i^T b\f$ for \f$i \in I_+ \cup I_-\f$ */
2984  SCIP_Real* vzlp, /**< array containing \f$v_i^T zlp_q\f$ for \f$i \in I_+ \cup I_-\f$ */
2985  SCIP_Real* wcoefs, /**< coefficients of w for the qud vars or NULL if w is 0 */
2986  SCIP_Real wzlp, /**< value of w at zlp */
2987  SCIP_Real kappa, /**< value of kappa */
2988  SCIP_ROWPREP* rowprep, /**< rowprep for the generated cut */
2989  SCIP_SOL* sol, /**< solution to separate */
2990  SCIP_Bool* success, /**< if a cut candidate could be computed */
2991  SCIP_Bool* strengthsuccess /**< if strengthening was successfully applied */
2992  )
2993 {
2994  SCIP_COL** cols;
2995  SCIP_ROW** rows;
2996  SCIP_Real* interpoints;
2997  SCIP_Real avecutcoef;
2998  int counter;
2999  int i;
3000 
3001  *success = TRUE;
3002  *strengthsuccess = FALSE;
3003 
3004  cols = SCIPgetLPCols(scip);
3005  rows = SCIPgetLPRows(scip);
3006 
3007  /* allocate memory for intersection points */
3008  SCIP_CALL( SCIPallocBufferArray(scip, &interpoints, rays->nrays) );
3009 
3010  /* compute all intersection points and store them in interpoints; build not-stregthened intersection cut */
3011  SCIP_CALL( computeIntercut(scip, nlhdlrdata, nlhdlrexprdata, rays, sidefactor, iscase4, vb, vzlp, wcoefs, wzlp, kappa,
3012  rowprep, interpoints, sol, success) );
3013 
3014  if( ! *success )
3015  goto CLEANUP;
3016 
3017  /* keep track of the number of attempted strengthenings and average cutcoef */
3018  counter = 0;
3019  avecutcoef = 0.0;
3020 
3021  /* go through all intersection points that are equal to infinity -> these correspond to the rays which are in the
3022  * recession cone of C, i.e. the rays for which we (possibly) can compute a negative steplength */
3023  for( i = 0; i < rays->nrays; ++i )
3024  {
3025  SCIP_Real rho;
3026  SCIP_Real cutcoef;
3027  int lppos;
3028 
3029  if( !SCIPisInfinity(scip, interpoints[i]) )
3030  continue;
3031 
3032  /* if we reached the limit of strengthenings, we stop */
3033  if( counter >= nlhdlrdata->nstrengthlimit )
3034  break;
3035 
3036  /* compute the smallest rho */
3037  SCIP_CALL( findRho(scip, nlhdlrdata, nlhdlrexprdata, rays, i, sidefactor, iscase4, vb, vzlp, wcoefs, wzlp, kappa,
3038  interpoints, &rho, success));
3039 
3040  /* compute cut coef */
3041  if( ! *success || SCIPisInfinity(scip, -rho) )
3042  cutcoef = 0.0;
3043  else
3044  cutcoef = 1.0 / rho;
3045 
3046  /* track average cut coef */
3047  counter += 1;
3048  avecutcoef += cutcoef;
3049 
3050  if( ! SCIPisZero(scip, cutcoef) )
3051  *strengthsuccess = TRUE;
3052 
3053  /* add var to cut: if variable is nonbasic at upper we have to flip sign of cutcoef */
3054  lppos = rays->lpposray[i];
3055  if( lppos < 0 )
3056  {
3057  lppos = -lppos - 1;
3058 
3059  assert(SCIProwGetBasisStatus(rows[lppos]) == SCIP_BASESTAT_LOWER || SCIProwGetBasisStatus(rows[lppos]) ==
3061 
3062  SCIP_CALL( addRowToCut(scip, rowprep, SCIProwGetBasisStatus(rows[lppos]) == SCIP_BASESTAT_UPPER ? cutcoef :
3063  -cutcoef, rows[lppos], success) ); /* rows have flipper base status! */
3064 
3065  if( ! *success )
3066  {
3067  INTERLOG(printf("Bad numeric: row not nonbasic enough\n");)
3068  nlhdlrdata->nbadnonbasic++;
3069  return SCIP_OKAY;
3070  }
3071  }
3072  else
3073  {
3074  assert(SCIPcolGetBasisStatus(cols[lppos]) == SCIP_BASESTAT_UPPER || SCIPcolGetBasisStatus(cols[lppos]) ==
3076  SCIP_CALL( addColToCut(scip, rowprep, sol, SCIPcolGetBasisStatus(cols[lppos]) == SCIP_BASESTAT_UPPER ? -cutcoef :
3077  cutcoef, cols[lppos]) );
3078  }
3079  }
3080 
3081  if( counter > 0 )
3082  nlhdlrdata->cutcoefsum += avecutcoef / counter;
3083 
3084 CLEANUP:
3085  SCIPfreeBufferArray(scip, &interpoints);
3086 
3087  return SCIP_OKAY;
3088 }
3089 
3090 /** sets variable in solution "vertex" to its nearest bound */
3091 static
3093  SCIP* scip, /**< SCIP data structure */
3094  SCIP_SOL* sol, /**< solution to separate */
3095  SCIP_SOL* vertex, /**< new solution to separate */
3096  SCIP_VAR* var, /**< var we want to find nearest bound to */
3097  SCIP_Real* factor, /**< is vertex for current var at lower or upper? */
3098  SCIP_Bool* success /**< TRUE if no variable is bounded */
3099  )
3100 {
3101  SCIP_Real solval;
3102  SCIP_Real bound;
3103 
3104  solval = SCIPgetSolVal(scip, sol, var);
3105  *success = TRUE;
3106 
3107  /* find nearest bound */
3108  if( SCIPisInfinity(scip, SCIPvarGetLbLocal(var)) && SCIPisInfinity(scip, SCIPvarGetUbLocal(var)) )
3109  {
3110  *success = FALSE;
3111  return SCIP_OKAY;
3112  }
3113  else if( solval - SCIPvarGetLbLocal(var) < SCIPvarGetUbLocal(var) - solval )
3114  {
3115  bound = SCIPvarGetLbLocal(var);
3116  *factor = 1.0;
3117  }
3118  else
3119  {
3120  bound = SCIPvarGetUbLocal(var);
3121  *factor = -1.0;
3122  }
3123 
3124  /* set val to bound in solution */
3125  SCIP_CALL( SCIPsetSolVal(scip, vertex, var, bound) );
3126  return SCIP_OKAY;
3127 }
3128 
3129 /** This function finds vertex (w.r.t. bounds of variables appearing in the quadratic) that is closest to the current
3130  * solution we want to separate.
3131  *
3132  * Furthermore, we store the rays corresponding to the unit vectors, i.e.,
3133  * - if \f$x_i\f$ is at its lower bound in vertex --> \f$r_i = e_i\f$
3134  * - if \f$x_i\f$ is at its upper bound in vertex --> \f$r_i = -e_i\f$
3135  */
3136 static
3138  SCIP* scip, /**< SCIP data structure */
3139  SCIP_NLHDLREXPRDATA* nlhdlrexprdata, /**< nlhdlr expression data */
3140  SCIP_SOL* sol, /**< solution to separate */
3141  SCIP_SOL* vertex, /**< new 'vertex' (w.r.t. bounds) solution to separate */
3142  SCIP_VAR* auxvar, /**< aux var of expr or NULL if not needed (e.g. separating real cons) */
3143  RAYS** raysptr, /**< pointer to new bound rays */
3144  SCIP_Bool* success /**< TRUE if no variable is unbounded */
3145  )
3146 {
3147  SCIP_EXPR* qexpr;
3148  SCIP_EXPR** linexprs;
3149  RAYS* rays;
3150  int nquadexprs;
3151  int nlinexprs;
3152  int raylength;
3153  int i;
3154 
3155  *success = TRUE;
3156 
3157  qexpr = nlhdlrexprdata->qexpr;
3158  SCIPexprGetQuadraticData(qexpr, NULL, &nlinexprs, &linexprs, NULL, &nquadexprs, NULL, NULL, NULL);
3159 
3160  raylength = (auxvar == NULL) ? nquadexprs + nlinexprs : nquadexprs + nlinexprs + 1;
3161 
3162  /* create rays */
3163  SCIP_CALL( createBoundRays(scip, raysptr, raylength) );
3164  rays = *raysptr;
3165 
3166  rays->rayssize = raylength;
3167  rays->nrays = raylength;
3168 
3169  /* go through quadratic variables */
3170  for( i = 0; i < nquadexprs; ++i )
3171  {
3172  SCIP_EXPR* expr;
3173  SCIPexprGetQuadraticQuadTerm(qexpr, i, &expr, NULL, NULL, NULL, NULL, NULL);
3174 
3175  rays->raysbegin[i] = i;
3176  rays->raysidx[i] = i;
3178 
3179  SCIP_CALL( setVarToNearestBound(scip, sol, vertex, SCIPgetExprAuxVarNonlinear(expr),
3180  &rays->rays[i], success) );
3181 
3182  if( ! *success )
3183  return SCIP_OKAY;
3184  }
3185 
3186  /* go through linear variables */
3187  for( i = 0; i < nlinexprs; ++i )
3188  {
3189  rays->raysbegin[i + nquadexprs] = i + nquadexprs;
3190  rays->raysidx[i + nquadexprs] = i + nquadexprs;
3191  rays->lpposray[i + nquadexprs] = SCIPcolGetLPPos(SCIPvarGetCol(SCIPgetExprAuxVarNonlinear(linexprs[i])));
3192 
3193  SCIP_CALL( setVarToNearestBound(scip, sol, vertex, SCIPgetExprAuxVarNonlinear(linexprs[i]),
3194  &rays->rays[i + nquadexprs], success) );
3195 
3196  if( ! *success )
3197  return SCIP_OKAY;
3198  }
3199 
3200  /* consider auxvar if it exists */
3201  if( auxvar != NULL )
3202  {
3203  rays->raysbegin[nquadexprs + nlinexprs] = nquadexprs + nlinexprs;
3204  rays->raysidx[nquadexprs + nlinexprs] = nquadexprs + nlinexprs;
3205  rays->lpposray[nquadexprs + nlinexprs] = SCIPcolGetLPPos(SCIPvarGetCol(auxvar));
3206 
3207  SCIP_CALL( setVarToNearestBound(scip, sol, vertex, auxvar, &rays->rays[nquadexprs + nlinexprs], success) );
3208 
3209  if( ! *success )
3210  return SCIP_OKAY;
3211  }
3212 
3213  rays->raysbegin[raylength] = raylength;
3214 
3215  return SCIP_OKAY;
3216 }
3217 
3218 /** checks if the quadratic constraint is violated by sol */
3219 static
3221  SCIP* scip, /**< SCIP data structure */
3222  SCIP_NLHDLREXPRDATA* nlhdlrexprdata, /**< nlhdlr expression data */
3223  SCIP_VAR* auxvar, /**< aux var of expr or NULL if not needed (e.g. separating real cons) */
3224  SCIP_SOL* sol, /**< solution to check feasibility for */
3225  SCIP_Real sidefactor /**< 1.0 if the violated constraint is q &le; rhs, -1.0 otherwise */
3226  )
3227 {
3228  SCIP_EXPR* qexpr;
3229  SCIP_EXPR** linexprs;
3230  SCIP_Real* lincoefs;
3231  SCIP_Real constant;
3232  SCIP_Real val;
3233  int nquadexprs;
3234  int nlinexprs;
3235  int nbilinexprs;
3236  int i;
3237 
3238  qexpr = nlhdlrexprdata->qexpr;
3239  SCIPexprGetQuadraticData(qexpr, &constant, &nlinexprs, &linexprs, &lincoefs, &nquadexprs,
3240  &nbilinexprs, NULL, NULL);
3241 
3242  val = 0.0;
3243 
3244  /* go through quadratic terms */
3245  for( i = 0; i < nquadexprs; i++ )
3246  {
3247  SCIP_EXPR* expr;
3248  SCIP_Real quadlincoef;
3249  SCIP_Real sqrcoef;
3250  SCIP_Real solval;
3251 
3252  SCIPexprGetQuadraticQuadTerm(qexpr, i, &expr, &quadlincoef, &sqrcoef, NULL, NULL, NULL);
3253 
3254  solval = SCIPgetSolVal(scip, sol, SCIPgetExprAuxVarNonlinear(expr));
3255 
3256  /* add square term */
3257  val += sqrcoef * SQR(solval);
3258 
3259  /* add linear term */
3260  val += quadlincoef * solval;
3261  }
3262 
3263  /* go through bilinear terms */
3264  for( i = 0; i < nbilinexprs; i++ )
3265  {
3266  SCIP_EXPR* expr1;
3267  SCIP_EXPR* expr2;
3268  SCIP_Real bilincoef;
3269 
3270  SCIPexprGetQuadraticBilinTerm(qexpr, i, &expr1, &expr2, &bilincoef, NULL, NULL);
3271 
3272  val += bilincoef * SCIPgetSolVal(scip, sol, SCIPgetExprAuxVarNonlinear(expr1))
3273  * SCIPgetSolVal(scip, sol, SCIPgetExprAuxVarNonlinear(expr2));
3274  }
3275 
3276  /* go through linear terms */
3277  for( i = 0; i < nlinexprs; i++ )
3278  {
3279  val += lincoefs[i] * SCIPgetSolVal(scip, sol, SCIPgetExprAuxVarNonlinear(linexprs[i]));
3280  }
3281 
3282  /* add auxvar if exists and get constant */
3283  if( auxvar != NULL )
3284  {
3285  val -= SCIPgetSolVal(scip, sol, auxvar);
3286 
3287  constant = (sidefactor == 1.0) ? constant - SCIPgetRhsNonlinear(nlhdlrexprdata->cons) :
3288  SCIPgetLhsNonlinear(nlhdlrexprdata->cons) - constant;
3289  }
3290  else
3291  constant = (sidefactor * constant);
3292 
3293  val = (sidefactor * val);
3294 
3295  /* now constraint is q(z) <= const */
3296  if( val <= constant )
3297  return FALSE;
3298  else
3299  return TRUE;
3300 }
3301 
3302 /** generates intersection cut that cuts off sol (which violates the quadratic constraint cons) */
3303 static
3305  SCIP* scip, /**< SCIP data structure */
3306  SCIP_EXPR* expr, /**< expr */
3307  SCIP_NLHDLRDATA* nlhdlrdata, /**< nlhdlr data */
3308  SCIP_NLHDLREXPRDATA* nlhdlrexprdata, /**< nlhdlr expression data */
3309  SCIP_CONS* cons, /**< violated constraint that contains expr */
3310  SCIP_SOL* sol, /**< solution to separate */
3311  SCIP_ROWPREP* rowprep, /**< rowprep for the generated cut */
3312  SCIP_Bool overestimate, /**< TRUE if viol cons is q(z) &ge; lhs; FALSE if q(z) &le; rhs */
3313  SCIP_Bool* success /**< whether separation was successfull or not */
3314  )
3315 {
3316  SCIP_EXPR* qexpr;
3317  RAYS* rays;
3318  SCIP_VAR* auxvar;
3319  SCIP_Real sidefactor;
3320  SCIP_Real* vb; /* eigenvectors * b */
3321  SCIP_Real* vzlp; /* eigenvectors * lpsol */
3322  SCIP_Real* wcoefs; /* coefficients affecting quadterms in w */
3323  SCIP_Real wzlp; /* w(lpsol) */
3324  SCIP_Real kappa;
3325  SCIP_Bool iscase4;
3326  SCIP_SOL* vertex;
3327  SCIP_SOL* soltoseparate;
3328  int nquadexprs;
3329  int nlinexprs;
3330  int i;
3331 
3332  /* count number of calls */
3333  (nlhdlrdata->ncalls++);
3334 
3335  qexpr = nlhdlrexprdata->qexpr;
3336  SCIPexprGetQuadraticData(qexpr, NULL, &nlinexprs, NULL, NULL, &nquadexprs, NULL, NULL, NULL);
3337 
3338 #ifdef DEBUG_INTERSECTIONCUT
3339  SCIPinfoMessage(scip, NULL, "Generating intersection cut for quadratic expr %p aka", (void*)expr);
3340  SCIP_CALL( SCIPprintExpr(scip, expr, NULL) );
3341  SCIPinfoMessage(scip, NULL, "\n");
3342 #endif
3343 
3344  *success = TRUE;
3345  iscase4 = TRUE;
3346 
3347  /* in nonbasic space cut is >= 1 */
3348  assert(SCIProwprepGetSide(rowprep) == 0.0);
3349  SCIProwprepAddSide(rowprep, 1.0);
3351  assert(SCIProwprepGetSide(rowprep) == 1.0);
3352 
3353  auxvar = (nlhdlrexprdata->cons != cons) ? SCIPgetExprAuxVarNonlinear(expr) : NULL;
3354  sidefactor = overestimate ? -1.0 : 1.0;
3355 
3356  rays = NULL;
3357 
3358  /* check if we use tableau or bounds as rays */
3359  if( ! nlhdlrdata->useboundsasrays )
3360  {
3361  SCIP_CALL( createAndStoreSparseRays(scip, nlhdlrexprdata, auxvar, &rays, success) );
3362 
3363  if( ! *success )
3364  {
3365  INTERLOG(printf("Failed to get rays: there is a var with base status ZERO!\n"); )
3366  return SCIP_OKAY;
3367  }
3368 
3369  soltoseparate = sol;
3370  }
3371  else
3372  {
3373  SCIP_Bool violated;
3374 
3375  if( auxvar != NULL )
3376  {
3377  *success = FALSE;
3378  return SCIP_OKAY;
3379  }
3380 
3381  /* create new solution */
3382  SCIP_CALL( SCIPcreateSol(scip, &vertex, NULL) );
3383 
3384  /* find nearest vertex of the box to separate and compute rays */
3385  SCIP_CALL( findVertexAndGetRays(scip, nlhdlrexprdata, sol, vertex, auxvar, &rays, success) );
3386 
3387  if( ! *success )
3388  {
3389  INTERLOG(printf("Failed to use bounds as rays: variable is unbounded!\n"); )
3390  freeRays(scip, &rays);
3391  SCIP_CALL( SCIPfreeSol(scip, &vertex) );
3392  return SCIP_OKAY;
3393  }
3394 
3395  /* check if vertex is violated */
3396  violated = isQuadConsViolated(scip, nlhdlrexprdata, auxvar, vertex, sidefactor);
3397 
3398  if( ! violated )
3399  {
3400  INTERLOG(printf("Failed to use bounds as rays: nearest vertex is not violated!\n"); )
3401  freeRays(scip, &rays);
3402  SCIP_CALL( SCIPfreeSol(scip, &vertex) );
3403  *success = FALSE;
3404  return SCIP_OKAY;
3405  }
3406 
3407  soltoseparate = vertex;
3408  }
3409 
3410  SCIP_CALL( SCIPallocBufferArray(scip, &vb, nquadexprs) );
3411  SCIP_CALL( SCIPallocBufferArray(scip, &vzlp, nquadexprs) );
3412  SCIP_CALL( SCIPallocBufferArray(scip, &wcoefs, nquadexprs) );
3413 
3414  SCIP_CALL( intercutsComputeCommonQuantities(scip, nlhdlrexprdata, auxvar, sidefactor, soltoseparate,
3415  vb, vzlp, wcoefs, &wzlp, &kappa) );
3416 
3417  /* check if we are in case 4 */
3418  if( nlinexprs == 0 && auxvar == NULL )
3419  {
3420  for( i = 0; i < nquadexprs; ++i )
3421  if( wcoefs[i] != 0.0 )
3422  break;
3423 
3424  if( i == nquadexprs )
3425  {
3426  /* from now on wcoefs is going to be NULL --> case 1, 2 or 3 */
3427  SCIPfreeBufferArray(scip, &wcoefs);
3428  iscase4 = FALSE;
3429  }
3430  }
3431 
3432  /* compute (strengthened) intersection cut */
3433  if( nlhdlrdata->usestrengthening )
3434  {
3435  SCIP_Bool strengthsuccess;
3436 
3437  SCIP_CALL( computeStrengthenedIntercut(scip, nlhdlrdata, nlhdlrexprdata, rays, sidefactor, iscase4, vb, vzlp, wcoefs,
3438  wzlp, kappa, rowprep, soltoseparate, success, &strengthsuccess) );
3439 
3440  if( *success && strengthsuccess )
3441  nlhdlrdata->nstrengthenings++;
3442  }
3443  else
3444  {
3445  SCIP_CALL( computeIntercut(scip, nlhdlrdata, nlhdlrexprdata, rays, sidefactor, iscase4, vb, vzlp, wcoefs, wzlp, kappa,
3446  rowprep, NULL, soltoseparate, success) );
3447  }
3448 
3449  SCIPfreeBufferArrayNull(scip, &wcoefs);
3450  SCIPfreeBufferArray(scip, &vzlp);
3451  SCIPfreeBufferArray(scip, &vb);
3452  freeRays(scip, &rays);
3453 
3454  if( nlhdlrdata->useboundsasrays )
3455  {
3456  SCIP_CALL( SCIPfreeSol(scip, &vertex) );
3457  }
3458 
3459  return SCIP_OKAY;
3460 }
3461 
3462 /** returns whether a quadratic form is "propagable"
3463  *
3464  * 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:
3465  * - it appears as a linear term (coef*expr)
3466  * - it appears as a square term (coef*expr^2)
3467  * - it appears in a bilinear term
3468  * - it appears in another bilinear term
3469  */
3470 static
3472  SCIP_EXPR* qexpr /**< quadratic representation data */
3473  )
3474 {
3475  int nquadexprs;
3476  int i;
3477 
3478  SCIPexprGetQuadraticData(qexpr, NULL, NULL, NULL, NULL, &nquadexprs, NULL, NULL, NULL);
3479 
3480  for( i = 0; i < nquadexprs; ++i )
3481  {
3482  SCIP_Real lincoef;
3483  SCIP_Real sqrcoef;
3484  int nadjbilin;
3485 
3486  SCIPexprGetQuadraticQuadTerm(qexpr, i, NULL, &lincoef, &sqrcoef, &nadjbilin, NULL, NULL);
3487 
3488  if( (lincoef != 0.0) + (sqrcoef != 0.0) + nadjbilin >= 2 ) /*lint !e514*/ /* actually MIN(2, nadjbilin), but we check >= 2 */
3489  return TRUE;
3490  }
3491 
3492  return FALSE;
3493 }
3494 
3495 /** returns whether a quadratic term is "propagable"
3496  *
3497  * 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:
3498  * - it appears as a linear term (coef*expr)
3499  * - it appears as a square term (coef*expr^2)
3500  * - it appears in a bilinear term
3501  * - it appears in another bilinear term
3502  */
3503 static
3505  SCIP_EXPR* qexpr, /**< quadratic representation data */
3506  int idx /**< index of quadratic term to consider */
3507  )
3508 {
3509  SCIP_Real lincoef;
3510  SCIP_Real sqrcoef;
3511  int nadjbilin;
3512 
3513  SCIPexprGetQuadraticQuadTerm(qexpr, idx, NULL, &lincoef, &sqrcoef, &nadjbilin, NULL, NULL);
3514 
3515  return (lincoef != 0.0) + (sqrcoef != 0.0) + nadjbilin >= 2; /*lint !e514*/ /* actually MIN(2, nadjbilin), but we check >= 2 */
3516 }
3517 
3518 /** solves a quadratic equation \f$ a\, \text{expr}^2 + b\, \text{expr} \in \text{rhs} \f$ (with \f$b\f$ an interval)
3519  * and reduces bounds on `expr` or deduces infeasibility if possible
3520  */
3521 static
3523  SCIP* scip, /**< SCIP data structure */
3524  SCIP_EXPR* expr, /**< expression for which to solve */
3525  SCIP_Real sqrcoef, /**< square coefficient */
3526  SCIP_INTERVAL b, /**< interval acting as linear coefficient */
3527  SCIP_INTERVAL rhs, /**< interval acting as rhs */
3528  SCIP_Bool* infeasible, /**< buffer to store if propagation produced infeasibility */
3529  int* nreductions /**< buffer to store the number of interval reductions */
3530  )
3531 {
3532  SCIP_INTERVAL a;
3533  SCIP_INTERVAL exprbounds;
3534  SCIP_INTERVAL newrange;
3535 
3536  assert(scip != NULL);
3537  assert(expr != NULL);
3538  assert(infeasible != NULL);
3539  assert(nreductions != NULL);
3540 
3541 #ifdef DEBUG_PROP
3542  SCIPinfoMessage(scip, NULL, "Propagating <expr> by solving a <expr>^2 + b <expr> in rhs, where <expr> is: ");
3543  SCIP_CALL( SCIPprintExpr(scip, expr, NULL) );
3544  SCIPinfoMessage(scip, NULL, "\n");
3545  SCIPinfoMessage(scip, NULL, "expr in [%g, %g], a = %g, b = [%g, %g] and rhs = [%g, %g]\n",
3547  SCIPintervalGetSup(SCIPgetExprBoundsNonlinear(scip, expr)), sqrcoef, b.inf, b.sup,
3548  rhs.inf, rhs.sup);
3549 #endif
3550 
3551  exprbounds = SCIPgetExprBoundsNonlinear(scip, expr);
3552  if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, exprbounds) )
3553  {
3554  *infeasible = TRUE;
3555  return SCIP_OKAY;
3556  }
3557 
3558  /* compute solution of a*x^2 + b*x in rhs */
3559  SCIPintervalSet(&a, sqrcoef);
3560  SCIPintervalSolveUnivariateQuadExpression(SCIP_INTERVAL_INFINITY, &newrange, a, b, rhs, exprbounds);
3561 
3562 #ifdef DEBUG_PROP
3563  SCIPinfoMessage(scip, NULL, "Solution [%g, %g]\n", newrange.inf, newrange.sup);
3564 #endif
3565 
3566  SCIP_CALL( SCIPtightenExprIntervalNonlinear(scip, expr, newrange, infeasible, nreductions) );
3567 
3568  return SCIP_OKAY;
3569 }
3570 
3571 /** 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 */
3572 static
3574  SCIP* scip, /**< SCIP data structure */
3575  SCIP_EXPR* expr, /**< expression for which to solve */
3576  SCIP_Real b, /**< linear coefficient */
3577  SCIP_INTERVAL rhs, /**< interval acting as rhs */
3578  SCIP_Bool* infeasible, /**< buffer to store if propagation produced infeasibility */
3579  int* nreductions /**< buffer to store the number of interval reductions */
3580  )
3581 {
3582  SCIP_INTERVAL newrange;
3583 
3584  assert(scip != NULL);
3585  assert(expr != NULL);
3586  assert(infeasible != NULL);
3587  assert(nreductions != NULL);
3588 
3589 #ifdef DEBUG_PROP
3590  SCIPinfoMessage(scip, NULL, "Propagating <expr> by solving %g <expr> in [%g, %g], where <expr> is: ", b, rhs.inf, rhs.sup);
3591  SCIP_CALL( SCIPprintExpr(scip, expr, NULL) );
3592  SCIPinfoMessage(scip, NULL, "\n");
3593 #endif
3594 
3595  /* compute solution of b*x in rhs */
3596  SCIPintervalDivScalar(SCIP_INTERVAL_INFINITY, &newrange, rhs, b);
3597 
3598 #ifdef DEBUG_PROP
3599  SCIPinfoMessage(scip, NULL, "Solution [%g, %g]\n", newrange.inf, newrange.sup);
3600 #endif
3601 
3602  SCIP_CALL( SCIPtightenExprIntervalNonlinear(scip, expr, newrange, infeasible, nreductions) );
3603 
3604  return SCIP_OKAY;
3605 }
3606 
3607 /** returns max of a/x - c*x for x in {x1, x2} with x1, x2 > 0 */
3608 static
3610  SCIP_Real a, /**< coefficient a */
3611  SCIP_Real c, /**< coefficient c */
3612  SCIP_Real x1, /**< coefficient x1 > 0 */
3613  SCIP_Real x2 /**< coefficient x2 > 0 */
3614  )
3615 {
3616  SCIP_Real cneg;
3617  SCIP_Real cand1;
3618  SCIP_Real cand2;
3619  SCIP_ROUNDMODE roundmode;
3620 
3621  assert(x1 > 0.0);
3622  assert(x2 > 0.0);
3623 
3624  cneg = SCIPintervalNegateReal(c);
3625 
3626  roundmode = SCIPintervalGetRoundingMode();
3628  cand1 = a/x1 + cneg*x1;
3629  cand2 = a/x2 + cneg*x2;
3630  SCIPintervalSetRoundingMode(roundmode);
3631 
3632  return MAX(cand1, cand2);
3633 }
3634 
3635 /** returns max of a/x - c*x for x in dom; it assumes that dom is contained in (0, +inf) */
3636 static
3638  SCIP_Real a, /**< coefficient a */
3639  SCIP_Real c, /**< coefficient c */
3640  SCIP_INTERVAL dom /**< domain of x */
3641  )
3642 {
3643  SCIP_ROUNDMODE roundmode;
3644  SCIP_INTERVAL argmax;
3645  SCIP_Real negunresmax;
3646  SCIP_Real boundarymax;
3647  assert(dom.inf > 0);
3648 
3649  /* if a >= 0, then the function is convex which means the maximum is at one of the boundaries
3650  *
3651  * if c = 0, then the function is monotone which means the maximum is also at one of the boundaries
3652  *
3653  * 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,
3654  * 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.
3655  * Otherwise (that is, c<0), the maximum is at one of the boundaries.
3656  */
3657  if( a >= 0.0 || c <= 0.0 )
3658  return computeMaxBoundaryForBilinearProp(a, c, dom.inf, dom.sup);
3659 
3660  /* now, the (unrestricted) maximum is at sqrt(-a/c).
3661  * if the argmax is not in the interior of dom then the solution is at a boundary, too
3662  * we check this by computing an interval that contains sqrt(-a/c) first
3663  */
3664  SCIPintervalSet(&argmax, -a);
3665  SCIPintervalDivScalar(SCIP_INTERVAL_INFINITY, &argmax, argmax, c);
3667 
3668  /* if the interval containing sqrt(-a/c) does not intersect with the interior of dom, then
3669  * the (restricted) maximum is at a boundary (we could even say at which boundary, but that doesn't save much)
3670  */
3671  if( argmax.sup <= dom.inf || argmax.inf >= dom.sup )
3672  return computeMaxBoundaryForBilinearProp(a, c, dom.inf, dom.sup);
3673 
3674  /* 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) */
3675  roundmode = SCIPintervalGetRoundingMode();
3677  negunresmax = 2.0*SCIPnextafter(sqrt(SCIPintervalNegateReal(a)*c), 0.0);
3678  SCIPintervalSetRoundingMode(roundmode);
3679 
3680  /* if the interval containing sqrt(-a/c) is contained in dom, then we can return -negunresmax */
3681  if( argmax.inf >= dom.inf && argmax.sup <= dom.sup )
3682  return -negunresmax;
3683 
3684  /* now what is left is the case where we cannot say for sure whether sqrt(-a/c) is contained in dom or not
3685  * so we are conservative and return the max of both cases, i.e.,
3686  * the max of the upper bounds on -2*sqrt(-a*c), a/dom.inf-c*dom.inf, a/dom.sup-c*dom.sup.
3687  */
3688  boundarymax = computeMaxBoundaryForBilinearProp(a, c, dom.inf, dom.sup);
3689  return MAX(boundarymax, -negunresmax);
3690 }
3691 
3692 /** computes the range of rhs/x - coef * x for x in exprdom; this is used for the propagation of bilinear terms
3693  *
3694  * 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
3695  * intended use of the function).
3696  * TODO: maybe check before calling it whether 0 is in the domain and then just avoid calling it
3697  *
3698  * 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].
3699  * 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].
3700  * However, this is the same as min of -B/x + coef*x and max of -A/x + coef*x for x in -[exprdom].
3701  * Thus, we can always reduce to x > 0 by multiplying [exprdom], rhs, and coef by -1.
3702  */
3703 static
3705  SCIP_INTERVAL exprdom, /**< expression for which to solve */
3706  SCIP_Real coef, /**< expression for which to solve */
3707  SCIP_INTERVAL rhs, /**< rhs used for computation */
3708  SCIP_INTERVAL* range /**< storage for the resulting range */
3709  )
3710 {
3711  SCIP_Real max;
3712  SCIP_Real min;
3713 
3714  if( exprdom.inf <= 0.0 && 0.0 <= exprdom.sup )
3715  {
3717  return;
3718  }
3719 
3720  /* reduce to positive case */
3721  if( exprdom.sup < 0 )
3722  {
3723  SCIPintervalMulScalar(SCIP_INTERVAL_INFINITY, &exprdom, exprdom, -1.0);
3725  coef *= -1.0;
3726  }
3727  assert(exprdom.inf > 0.0);
3728 
3729  /* compute maximum and minimum */
3730  max = computeMaxForBilinearProp(rhs.sup, coef, exprdom);
3731  min = -computeMaxForBilinearProp(-rhs.inf, -coef, exprdom);
3732 
3733  /* set interval */
3734  SCIPintervalSetBounds(range, min, max);
3735 }
3736 
3737 /** reverse propagates coef_i expr_i + constant in rhs */
3738 static
3740  SCIP* scip, /**< SCIP data structure */
3741  SCIP_EXPR** linexprs, /**< linear expressions */
3742  int nlinexprs, /**< number of linear expressions */
3743  SCIP_Real* lincoefs, /**< coefficients of linear expressions */
3744  SCIP_Real constant, /**< constant */
3745  SCIP_INTERVAL rhs, /**< rhs */
3746  SCIP_Bool* infeasible, /**< buffer to store whether an exps' bounds were propagated to an empty interval */
3747  int* nreductions /**< buffer to store the number of interval reductions of all exprs */
3748  )
3749 {
3750  SCIP_INTERVAL* oldboundslin;
3751  SCIP_INTERVAL* newboundslin;
3752  int i;
3753 
3754  if( nlinexprs == 0 )
3755  return SCIP_OKAY;
3756 
3757  SCIP_CALL( SCIPallocBufferArray(scip, &oldboundslin, nlinexprs) );
3758  SCIP_CALL( SCIPallocBufferArray(scip, &newboundslin, nlinexprs) );
3759 
3760  for( i = 0; i < nlinexprs; ++i )
3761  oldboundslin[i] = SCIPexprGetActivity(linexprs[i]); /* TODO use SCIPgetExprBoundsNonlinear(scip, linexprs[i]) ? */
3762 
3763  *nreductions = SCIPintervalPropagateWeightedSum(SCIP_INTERVAL_INFINITY, nlinexprs,
3764  oldboundslin, lincoefs, constant, rhs, newboundslin, infeasible);
3765 
3766  if( *nreductions > 0 && !*infeasible )
3767  {
3768  /* SCIP is more conservative with what constitutes a reduction than interval arithmetic so we follow SCIP */
3769  *nreductions = 0;
3770  for( i = 0; i < nlinexprs && ! (*infeasible); ++i )
3771  {
3772  SCIP_CALL( SCIPtightenExprIntervalNonlinear(scip, linexprs[i], newboundslin[i], infeasible, nreductions) );
3773  }
3774  }
3775 
3776  SCIPfreeBufferArray(scip, &newboundslin);
3777  SCIPfreeBufferArray(scip, &oldboundslin);
3778 
3779  return SCIP_OKAY;
3780 }
3781 
3782 
3783 /*
3784  * Callback methods of nonlinear handler
3785  */
3786 
3787 /** callback to free expression specific data */
3788 static
3789 SCIP_DECL_NLHDLRFREEEXPRDATA(nlhdlrFreeexprdataQuadratic)
3790 { /*lint --e{715}*/
3791  assert(nlhdlrexprdata != NULL);
3792  assert(*nlhdlrexprdata != NULL);
3793 
3794  if( (*nlhdlrexprdata)->quadactivities != NULL )
3795  {
3796  int nquadexprs;
3797  SCIPexprGetQuadraticData((*nlhdlrexprdata)->qexpr, NULL, NULL, NULL, NULL, &nquadexprs, NULL, NULL, NULL);
3798  SCIPfreeBlockMemoryArray(scip, &(*nlhdlrexprdata)->quadactivities, nquadexprs);
3799  }
3800 
3801  SCIPfreeBlockMemory(scip, nlhdlrexprdata);
3802 
3803  return SCIP_OKAY;
3804 }
3805 
3806 /** callback to detect structure in expression tree
3807  *
3808  * A term is quadratic if
3809  * - it is a product expression of two expressions, or
3810  * - it is power expression of an expression with exponent 2.0.
3811  *
3812  * We define a _propagable_ quadratic expression as a quadratic expression whose termwise propagation does not yield the
3813  * best propagation. In other words, is a quadratic expression that suffers from the dependency problem.
3814  *
3815  * Specifically, a propagable quadratic expression is a sum expression such that there is at least one expr that appears
3816  * at least twice (because of simplification, this means it appears in a quadratic terms and somewhere else).
3817  * For example: \f$x^2 + y^2\f$ is not a propagable quadratic expression; \f$x^2 + x\f$ is a propagable quadratic expression;
3818  * \f$x^2 + x y\f$ is also a propagable quadratic expression
3819  *
3820  * Furthermore, we distinguish between propagable and non-propagable terms. A term is propagable if any of the expressions
3821  * involved in it appear somewhere else. For example, \f$xy + z^2 + z\f$ is a propagable quadratic, the term \f$xy\f$ is
3822  * non-propagable, and \f$z^2\f$ is propagable. For propagation, non-propagable terms are handled as if they were linear
3823  * 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
3824  * 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),
3825  * 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$.
3826  *
3827  * For propagation, we store the quadratic in our data structure in the following way: We count how often a variable
3828  * appears. Then, a bilinear product expr_i * expr_j is stored as expr_i * expr_j if # expr_i appears > # expr_j
3829  * appears. When # expr_i appears = # expr_j appears, it then it will be stored as expr_i * expr_j if and only if
3830  * expr_i < expr_j, where '<' is the expression order (see \ref EXPR_ORDER "Ordering Rules" in \ref scip_expr.h).
3831  * Heuristically, this should be useful for propagation. The intuition is that by factoring out the variable that
3832  * appears most often we should be able to take care of the dependency problem better.
3833  *
3834  * Simple convex quadratics like \f$x^2 + y^2\f$ are ignored since the default nlhdlr will take care of them.
3835  *
3836  * @note The expression needs to be simplified (in particular, it is assumed to be sorted).
3837  * @note Common subexpressions are also assumed to have been identified, the hashing will fail otherwise!
3838  *
3839  * Sorted implies that:
3840  * - expr < expr^2: bases are the same, but exponent 1 < 2
3841  * - expr < expr * other_expr: u*v < w holds if and only if v < w (OR8), but here w = u < v, since expr comes before
3842  * other_expr in the product
3843  * - expr < other_expr * expr: u*v < w holds if and only if v < w (OR8), but here v = w
3844  *
3845  * Thus, if we see somebody twice, it is a propagable quadratic.
3846  *
3847  * It also implies that
3848  * - expr^2 < expr * other_expr
3849  * - other_expr * expr < expr^2
3850  *
3851  * 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.
3852  */
3853 static
3854 SCIP_DECL_NLHDLRDETECT(nlhdlrDetectQuadratic)
3855 { /*lint --e{715,774}*/
3856  SCIP_NLHDLREXPRDATA* nlexprdata;
3857  SCIP_NLHDLRDATA* nlhdlrdata;
3858  SCIP_Real* eigenvalues;
3859  SCIP_Bool isquadratic;
3860  SCIP_Bool propagable;
3861 
3862  assert(scip != NULL);
3863  assert(nlhdlr != NULL);
3864  assert(expr != NULL);
3865  assert(enforcing != NULL);
3866  assert(participating != NULL);
3867  assert(nlhdlrexprdata != NULL);
3868 
3869  nlhdlrdata = SCIPnlhdlrGetData(nlhdlr);
3870  assert(nlhdlrdata != NULL);
3871 
3872  /* don't check if all enforcement methods are already ensured */
3873  if( (*enforcing & SCIP_NLHDLR_METHOD_ALL) == SCIP_NLHDLR_METHOD_ALL )
3874  return SCIP_OKAY;
3875 
3876  /* if it is not a sum of at least two terms, it is not interesting */
3877  /* TODO: constraints of the form l<= x*y <= r ? */
3878  if( ! SCIPisExprSum(scip, expr) || SCIPexprGetNChildren(expr) < 2 )
3879  return SCIP_OKAY;
3880 
3881  /* If we are in a subSCIP we don't want to separate intersection cuts */
3882  if( SCIPgetSubscipDepth(scip) > 0 )
3883  nlhdlrdata->useintersectioncuts = FALSE;
3884 
3885 #ifdef SCIP_DEBUG
3886  SCIPinfoMessage(scip, NULL, "Nlhdlr quadratic detecting expr %p aka ", (void*)expr);
3887  SCIP_CALL( SCIPprintExpr(scip, expr, NULL) );
3888  SCIPinfoMessage(scip, NULL, "\n");
3889  SCIPinfoMessage(scip, NULL, "Have to enforce %d\n", *enforcing);
3890 #endif
3891 
3892  /* check whether expression is quadratic (a sum with at least one square or bilinear term) */
3893  SCIP_CALL( SCIPcheckExprQuadratic(scip, expr, &isquadratic) );
3894 
3895  /* not quadratic -> nothing for us */
3896  if( !isquadratic )
3897  {
3898  SCIPdebugMsg(scip, "expr %p is not quadratic -> abort detect\n", (void*)expr);
3899  return SCIP_OKAY;
3900  }
3901 
3902  propagable = isPropagable(expr);
3903 
3904  /* if we are not propagable and are in presolving, return */
3905  if( !propagable && SCIPgetStage(scip) == SCIP_STAGE_PRESOLVING )
3906  {
3907  SCIPdebugMsg(scip, "expr %p is not propagable and in presolving -> abort detect\n", (void*)expr);
3908  return SCIP_OKAY;
3909  }
3910 
3911  /* if we do not use intersection cuts and are not propagable, then we do not want to handle it at all;
3912  * if not propagable, then we need to check the curvature to decide if we want to generate intersection cuts
3913  */
3914  if( !propagable && !nlhdlrdata->useintersectioncuts )
3915  {
3916  SCIPdebugMsg(scip, "expr %p is not propagable -> abort detect\n", (void*)expr);
3917  return SCIP_OKAY;
3918  }
3919 
3920  /* store quadratic in nlhdlrexprdata */
3921  SCIP_CALL( SCIPallocClearBlockMemory(scip, nlhdlrexprdata) );
3922  nlexprdata = *nlhdlrexprdata;
3923  nlexprdata->qexpr = expr;
3924  nlexprdata->cons = cons;
3925 
3926 #ifdef DEBUG_DETECT
3927  SCIPinfoMessage(scip, NULL, "Nlhdlr quadratic detected:\n");
3928  SCIP_CALL( SCIPprintExprQuadratic(scip, conshdlr, qexpr) );
3929 #endif
3930 
3931  /* every propagable quadratic expression will be handled since we can propagate */
3932  if( propagable )
3933  {
3934  SCIP_EXPR** linexprs;
3935  int nlinexprs;
3936  int nquadexprs;
3937  int nbilin;
3938  int i;
3939 
3940  *participating |= SCIP_NLHDLR_METHOD_ACTIVITY;
3941  *enforcing |= SCIP_NLHDLR_METHOD_ACTIVITY;
3942 
3943  SCIPexprGetQuadraticData(expr, NULL, &nlinexprs, &linexprs, NULL, &nquadexprs, &nbilin, NULL, NULL);
3944  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &nlexprdata->quadactivities, nquadexprs) );
3945 
3946  /* notify children of quadratic that we will need their activity for propagation */
3947  for( i = 0; i < nlinexprs; ++i )
3949 
3950  for( i = 0; i < nquadexprs; ++i )
3951  {
3952  SCIP_EXPR* argexpr;
3953  if( isPropagableTerm(expr, i) )
3954  {
3955  SCIPexprGetQuadraticQuadTerm(expr, i, &argexpr, NULL, NULL, &nbilin, NULL, NULL);
3957 
3958 #ifdef DEBUG_DETECT
3959  SCIPinfoMessage(scip, NULL, "quadterm %d propagable, using %p, unbounded=%d\n", i, (void*)argexpr, nbilin >
3961 #endif
3962  }
3963  else
3964  {
3965  /* non-propagable quadratic is either a single square term or a single bilinear term
3966  * we should make use nlhdlrs in pow or product for this term, so we register usage of the square or product
3967  * expr instead of argexpr
3968  */
3969  SCIP_EXPR* sqrexpr;
3970  int* adjbilin;
3971 
3972  SCIPexprGetQuadraticQuadTerm(expr, i, &argexpr, NULL, NULL, &nbilin, &adjbilin, &sqrexpr);
3973 
3974  if( sqrexpr != NULL )
3975  {
3977  assert(nbilin == 0);
3978 
3979 #ifdef DEBUG_DETECT
3980  SCIPinfoMessage(scip, NULL, "quadterm %d non-propagable square, using %p\n", i, (void*)sqrexpr);
3981 #endif
3982  }
3983  else
3984  {
3985  /* we have expr1 * other_expr or other_expr * expr1; know that expr1 is non propagable, but to decide if
3986  * we want the bounds of expr1 or of the product expr1 * other_expr (or other_expr * expr1), we have to
3987  * decide whether other_expr is also non propagable; due to the way we sort bilinear terms (by
3988  * frequency), we can deduce that other_expr doesn't appear anywhere else (i.e. is non propagable) if the
3989  * product is of the form expr1 * other_expr; however, if we see other_expr * expr1 we need to find
3990  * other_expr and check whether it is propagable
3991  */
3992  SCIP_EXPR* expr1;
3993  SCIP_EXPR* prodexpr;
3994 
3995  assert(nbilin == 1);
3996  SCIPexprGetQuadraticBilinTerm(expr, adjbilin[0], &expr1, NULL, NULL, NULL, &prodexpr);
3997 
3998  if( expr1 == argexpr )
3999  {
4001 
4002 #ifdef DEBUG_DETECT
4003  SCIPinfoMessage(scip, NULL, "quadterm %d non-propagable product, using %p\n", i, (void*)prodexpr);
4004 #endif
4005  }
4006  else
4007  {
4008  int j;
4009  /* check if other_expr is propagable in which case we need the bounds of expr1; otherwise we just need
4010  * the bounds of the product and this will be (or was) registered when the loop takes us to the
4011  * quadexpr other_expr.
4012  * TODO this should be done faster, maybe store pos1 in bilinexprterm or store quadexprterm's in bilinexprterm
4013  */
4014  for( j = 0; j < nquadexprs; ++j )
4015  {
4016  SCIP_EXPR* exprj;
4017  SCIPexprGetQuadraticQuadTerm(expr, j, &exprj, NULL, NULL, NULL, NULL, NULL);
4018  if( expr1 == exprj )
4019  {
4020  if( isPropagableTerm(expr, j) )
4021  {
4023 #ifdef DEBUG_DETECT
4024  SCIPinfoMessage(scip, NULL, "quadterm %d non-propagable alien product, using %p\n", i, (void*)argexpr);
4025 #endif
4026  }
4027  break;
4028  }
4029  }
4030  }
4031  }
4032  }
4033  }
4034  }
4035 
4036  /* check if we are going to separate or not */
4037  nlexprdata->curvature = SCIP_EXPRCURV_UNKNOWN;
4038 
4039  /* for now, we do not care about separation if it is not required */
4040  if( (*enforcing & SCIP_NLHDLR_METHOD_SEPABOTH) == SCIP_NLHDLR_METHOD_SEPABOTH )
4041  {
4042  /* if nobody can do anything, remove data */
4043  if( *participating == SCIP_NLHDLR_METHOD_NONE ) /*lint !e845*/
4044  {
4045  SCIP_CALL( nlhdlrFreeexprdataQuadratic(scip, nlhdlr, expr, nlhdlrexprdata) );
4046  }
4047  else
4048  {
4049  SCIPdebugMsg(scip, "expr %p is quadratic and propagable -> propagate\n", (void*)expr);
4050  }
4051  return SCIP_OKAY;
4052  }
4053 
4054  assert(SCIPgetStage(scip) >= SCIP_STAGE_INITSOLVE); /* separation should only be required in (init)solving stage */
4055 
4056  /* check if we can do something more: check curvature of quadratic function stored in nlexprdata
4057  * this is currently only used to decide whether we want to separate, so it can be skipped if in presolve
4058  */
4059  SCIPdebugMsg(scip, "checking curvature of expr %p\n", (void*)expr);
4060  SCIP_CALL( SCIPcomputeExprQuadraticCurvature(scip, expr, &nlexprdata->curvature, NULL, nlhdlrdata->useintersectioncuts) );
4061 
4062  /* get eigenvalues to be able to check whether they were computed */
4063  SCIPexprGetQuadraticData(expr, NULL, NULL, NULL, NULL, NULL, NULL, &eigenvalues, NULL);
4064 
4065  /* if we use intersection cuts then we can handle any non-convex quadratic */
4066  if( nlhdlrdata->useintersectioncuts && eigenvalues != NULL && (*enforcing & SCIP_NLHDLR_METHOD_SEPABELOW) ==
4067  FALSE && nlexprdata->curvature != SCIP_EXPRCURV_CONVEX )
4068  {
4069  *participating |= SCIP_NLHDLR_METHOD_SEPABELOW;
4070  }
4071 
4072  if( nlhdlrdata->useintersectioncuts && eigenvalues != NULL && (*enforcing & SCIP_NLHDLR_METHOD_SEPAABOVE) == FALSE &&
4073  nlexprdata->curvature != SCIP_EXPRCURV_CONCAVE )
4074  {
4075  *participating |= SCIP_NLHDLR_METHOD_SEPAABOVE;
4076  }
4077 
4078  /* if nobody can do anything, remove data */
4079  if( *participating == SCIP_NLHDLR_METHOD_NONE ) /*lint !e845*/
4080  {
4081  SCIP_CALL( nlhdlrFreeexprdataQuadratic(scip, nlhdlr, expr, nlhdlrexprdata) );
4082  return SCIP_OKAY;
4083  }
4084 
4085  /* we only need auxiliary variables if we are going to separate */
4086  if( *participating & SCIP_NLHDLR_METHOD_SEPABOTH )
4087  {
4088  SCIP_EXPR** linexprs;
4089  int nquadexprs;
4090  int nlinexprs;
4091  int i;
4092 
4093  SCIPexprGetQuadraticData(expr, NULL, &nlinexprs, &linexprs, NULL, &nquadexprs, NULL, NULL, NULL);
4094 
4095  for( i = 0; i < nlinexprs; ++i ) /* expressions appearing linearly */
4096  {
4098  }
4099  for( i = 0; i < nquadexprs; ++i ) /* expressions appearing quadratically */
4100  {
4101  SCIP_EXPR* quadexpr;
4102  SCIPexprGetQuadraticQuadTerm(expr, i, &quadexpr, NULL, NULL, NULL, NULL, NULL);
4104  }
4105 
4106  SCIPdebugMsg(scip, "expr %p is quadratic and propagable -> propagate and separate\n", (void*)expr);
4107 
4108  nlexprdata->separating = TRUE;
4109  }
4110  else
4111  {
4112  SCIPdebugMsg(scip, "expr %p is quadratic and propagable -> propagate only\n", (void*)expr);
4113  }
4114 
4116  {
4117  SCIPexprSetCurvature(expr, nlexprdata->curvature);
4118  SCIPdebugMsg(scip, "expr is %s in the original variables\n", nlexprdata->curvature == SCIP_EXPRCURV_CONCAVE ? "concave" : "convex");
4119  nlexprdata->origvars = TRUE;
4120  }
4121 
4122  return SCIP_OKAY;
4123 }
4124 
4125 /** nonlinear handler auxiliary evaluation callback */
4126 static
4127 SCIP_DECL_NLHDLREVALAUX(nlhdlrEvalauxQuadratic)
4128 { /*lint --e{715}*/
4129  int i;
4130  int nlinexprs;
4131  int nquadexprs;
4132  int nbilinexprs;
4133  SCIP_Real constant;
4134  SCIP_Real* lincoefs;
4135  SCIP_EXPR** linexprs;
4136 
4137  assert(scip != NULL);
4138  assert(expr != NULL);
4139  assert(auxvalue != NULL);
4140  assert(nlhdlrexprdata->separating);
4141  assert(nlhdlrexprdata->qexpr == expr);
4142 
4143  /* if the quadratic is in the original variable we can just evaluate the expression */
4144  if( nlhdlrexprdata->origvars )
4145  {
4146  *auxvalue = SCIPexprGetEvalValue(expr);
4147  return SCIP_OKAY;
4148  }
4149 
4150  /* TODO there was a
4151  *auxvalue = SCIPevalExprQuadratic(scip, nlhdlrexprdata->qexpr, sol);
4152  here; any reason why not using this anymore?
4153  */
4154 
4155  SCIPexprGetQuadraticData(expr, &constant, &nlinexprs, &linexprs, &lincoefs, &nquadexprs, &nbilinexprs, NULL, NULL);
4156 
4157  *auxvalue = constant;
4158 
4159  for( i = 0; i < nlinexprs; ++i ) /* linear exprs */
4160  *auxvalue += lincoefs[i] * SCIPgetSolVal(scip, sol, SCIPgetExprAuxVarNonlinear(linexprs[i]));
4161 
4162  for( i = 0; i < nquadexprs; ++i ) /* quadratic terms */
4163  {
4164  SCIP_Real solval;
4165  SCIP_Real lincoef;
4166  SCIP_Real sqrcoef;
4167  SCIP_EXPR* qexpr;
4168 
4169  SCIPexprGetQuadraticQuadTerm(expr, i, &qexpr, &lincoef, &sqrcoef, NULL, NULL, NULL);
4170 
4171  solval = SCIPgetSolVal(scip, sol, SCIPgetExprAuxVarNonlinear(qexpr));
4172  *auxvalue += (lincoef + sqrcoef * solval) * solval;
4173  }
4174 
4175  for( i = 0; i < nbilinexprs; ++i ) /* bilinear terms */
4176  {
4177  SCIP_EXPR* expr1;
4178  SCIP_EXPR* expr2;
4179  SCIP_Real coef;
4180 
4181  SCIPexprGetQuadraticBilinTerm(expr, i, &expr1, &expr2, &coef, NULL, NULL);
4182 
4183  *auxvalue += coef * SCIPgetSolVal(scip, sol, SCIPgetExprAuxVarNonlinear(expr1)) * SCIPgetSolVal(scip, sol,
4185  }
4186 
4187  return SCIP_OKAY;
4188 }
4189 
4190 /** nonlinear handler enforcement callback */
4191 static
4192 SCIP_DECL_NLHDLRENFO(nlhdlrEnfoQuadratic)
4193 { /*lint --e{715}*/
4194  SCIP_NLHDLRDATA* nlhdlrdata;
4195  SCIP_ROWPREP* rowprep;
4196  SCIP_Bool success = FALSE;
4197  SCIP_NODE* node;
4198  int depth;
4199  SCIP_Longint nodenumber;
4200  SCIP_Real* eigenvalues;
4201  SCIP_Real violation;
4202 
4203  assert(nlhdlrexprdata != NULL);
4204  assert(nlhdlrexprdata->qexpr == expr);
4205 
4206  INTERLOG(printf("Starting interesection cuts!\n");)
4207 
4208  nlhdlrdata = SCIPnlhdlrGetData(nlhdlr);
4209  assert(nlhdlrdata != NULL);
4210 
4211  assert(result != NULL);
4212  *result = SCIP_DIDNOTRUN;
4213 
4214  if( branchcandonly )
4215  return SCIP_OKAY;
4216 
4217  /* estimate should take care of convex quadratics */
4218  if( ( overestimate && nlhdlrexprdata->curvature == SCIP_EXPRCURV_CONCAVE) ||
4219  (!overestimate && nlhdlrexprdata->curvature == SCIP_EXPRCURV_CONVEX) )
4220  {
4221  INTERLOG(printf("Convex or concave, no need of interesection cuts!\n");)
4222  return SCIP_OKAY;
4223  }
4224 
4225  /* nothing to do if we can't use intersection cuts */
4226  if( ! nlhdlrdata->useintersectioncuts )
4227  {
4228  INTERLOG(printf("We don't use intersection cuts!\n");)
4229  return SCIP_OKAY;
4230  }
4231 
4232  /* right now can use interesction cuts only if a basic LP solution is at hand; TODO: in principle we can do something
4233  * even if it is not optimal
4234  */
4236  {
4237  INTERLOG(printf("LP solutoin not good!\n");)
4238  return SCIP_OKAY;
4239  }
4240 
4241  /* only separate at selected nodes */
4242  node = SCIPgetCurrentNode(scip);
4243  depth = SCIPnodeGetDepth(node);
4244  if( (nlhdlrdata->atwhichnodes == -1 && depth != 0) || (nlhdlrdata->atwhichnodes != -1 && depth % nlhdlrdata->atwhichnodes != 0) )
4245  {
4246  INTERLOG(printf("Don't separate at this node\n");)
4247  return SCIP_OKAY;
4248  }
4249 
4250  /* do not add more than ncutslimitroot cuts in root node and ncutslimit cuts in the non-root nodes */
4251  nodenumber = SCIPnodeGetNumber(node);
4252  if( nlhdlrdata->lastnodenumber != nodenumber )
4253  {
4254  nlhdlrdata->lastnodenumber = nodenumber;
4255  nlhdlrdata->lastncuts = nlhdlrdata->ncutsadded;
4256  }
4257  /*else if( (depth > 0 && nlhdlrdata->ncutsadded - nlhdlrdata->lastncuts >= nlhdlrdata->ncutslimit) || (depth == 0 &&
4258  nlhdlrdata->ncutsadded - nlhdlrdata->lastncuts >= nlhdlrdata->ncutslimitroot)) */
4259  /* allow the addition of a certain number of cuts per quadratic */
4260  if( (depth > 0 && nlhdlrexprdata->ncutsadded >= nlhdlrdata->ncutslimit) || (depth == 0 &&
4261  nlhdlrexprdata->ncutsadded >= nlhdlrdata->ncutslimitroot) )
4262  {
4263  INTERLOG(printf("Too many cuts added already\n");)
4264  return SCIP_OKAY;
4265  }
4266 
4267  /* can't separate if we do not have an eigendecomposition */
4268  SCIPexprGetQuadraticData(expr, NULL, NULL, NULL, NULL, NULL, NULL, &eigenvalues, NULL);
4269  if( eigenvalues == NULL )
4270  {
4271  INTERLOG(printf("No known eigenvalues!\n");)
4272  return SCIP_OKAY;
4273  }
4274 
4275  /* if constraint is not sufficiently violated -> do nothing */
4276  if( cons != nlhdlrexprdata->cons )
4277  {
4278  /* constraint is w.r.t auxvar */
4279  violation = auxvalue - SCIPgetSolVal(scip, NULL, SCIPgetExprAuxVarNonlinear(expr));
4280  violation = ABS( violation );
4281  }
4282  else
4283  /* quadratic is a constraint */
4284  violation = MAX( SCIPgetLhsNonlinear(nlhdlrexprdata->cons) - auxvalue, auxvalue -
4285  SCIPgetRhsNonlinear(nlhdlrexprdata->cons)); /*lint !e666*/
4286 
4287  if( violation < nlhdlrdata->minviolation )
4288  {
4289  INTERLOG(printf("Violation %g is just too small\n", violation); )
4290  return SCIP_OKAY;
4291  }
4292 
4293  /* we can't build an intersection cut when the expr is the root of some constraint and also a subexpression of
4294  * another constraint because we initialize data differently */
4295  if( nlhdlrexprdata->cons != NULL && cons != nlhdlrexprdata->cons )
4296  {
4297  INTERLOG(printf("WARNING!! expr is root of one constraint and subexpr of another!\n"); )
4298  return SCIP_OKAY;
4299  }
4300 
4301  /* if we are the root of a constraint and we are feasible w.r.t our auxiliary variables, that is, auxvalue is
4302  * actually feasible for the sides of the constraint, then do not separate
4303  */
4304  if( cons == nlhdlrexprdata->cons && ((overestimate && (SCIPgetLhsNonlinear(cons)) - auxvalue < SCIPfeastol(scip)) ||
4305  (! overestimate && (auxvalue - SCIPgetRhsNonlinear(cons) < SCIPfeastol(scip)))) )
4306  {
4307  INTERLOG(printf("We are actually feasible for the sides of the constraint\n"); )
4308  return SCIP_OKAY;
4309  }
4310 
4311 #ifdef DEBUG_INTERSECTIONCUT
4312  SCIPinfoMessage(scip, NULL, "Build intersection cut for \n");
4313  if( cons == nlhdlrexprdata->cons )
4314  {
4315  SCIP_CALL( SCIPprintCons(scip, cons, NULL) );
4316  SCIPinfoMessage(scip, NULL, "\n");
4317  }
4318  else
4319  {
4320  SCIP_CALL( SCIPprintExpr(scip, expr, NULL) );
4322  }
4323  SCIPinfoMessage(scip, NULL, "We need to %sestimate\n", overestimate ? "over" : "under" );
4324  SCIPinfoMessage(scip, NULL, "LP sol: \n");
4326 #endif
4327  *result = SCIP_DIDNOTFIND;
4328 
4329  /* cut (in the nonbasic space) is of the form alpha^T x >= 1 */
4331  INTERLOG(printf("Generating inter cut\n"); )
4332 
4333  SCIP_CALL( generateIntercut(scip, expr, nlhdlrdata, nlhdlrexprdata, cons, sol, rowprep, overestimate, &success) );
4334  INTERLOG(if( !success) printf("Generation failed\n"); )
4335 
4336  /* we generated something, let us see if it survives the clean up */
4337  if( success )
4338  {
4339  assert(sol == NULL);
4340  nlhdlrdata->ncutsgenerated += 1;
4341  nlhdlrexprdata->ncutsadded += 1;
4342 
4343  /* merge coefficients that belong to same variable */
4344  SCIPmergeRowprepTerms(scip, rowprep);
4345 
4346  /* sparsify cut */
4347  if( nlhdlrdata->sparsifycuts )
4348  sparsifyIntercut(scip, rowprep);
4349 
4350  SCIP_CALL( SCIPcleanupRowprep(scip, rowprep, sol, nlhdlrdata->mincutviolation, &violation, &success) );
4351  INTERLOG(if( !success) printf("Clean up failed\n"); )
4352  }
4353 
4354  /* if cut looks good (numerics ok and cutting off solution), then turn into row and add to sepastore */
4355  if( success )
4356  {
4357  SCIP_ROW* row;
4358  SCIP_Bool infeasible;
4359 
4360  /* count number of bound cuts */
4361  if( nlhdlrdata->useboundsasrays )
4362  nlhdlrdata->nboundcuts += 1;
4363 
4364  (void) SCIPsnprintf(SCIProwprepGetName(rowprep), SCIP_MAXSTRLEN, "%s_intersection_quadratic%p_lp%" SCIP_LONGINT_FORMAT,
4365  overestimate ? "over" : "under",
4366  (void*)expr,
4367  SCIPgetNLPs(scip));
4368 
4369  SCIP_CALL( SCIPgetRowprepRowCons(scip, &row, rowprep, cons) );
4370 
4371  /* printf("## New cut\n");
4372  printf(" -> found maxquad-free cut <%s>: act=%f, lhs=%f, norm=%f, eff=%f, min=%f, max=%f (range=%f)\n\n",
4373  SCIProwGetName(row), SCIPgetRowLPActivity(scip, row), SCIProwGetLhs(row), SCIProwGetNorm(row),
4374  SCIPgetCutEfficacy(scip, NULL, row),
4375  SCIPgetRowMinCoef(scip, row), SCIPgetRowMaxCoef(scip, row),
4376  SCIPgetRowMaxCoef(scip, row)/SCIPgetRowMinCoef(scip, row)); */
4377 
4378  /* intersection cuts can be numerically nasty; we do some extra numerical checks here */
4379  /*printf("SCIP DEPTH %d got a cut with violation %g, efficacy %g and r/e %g\n", SCIPgetSubscipDepth(scip),
4380  * violation, SCIPgetCutEfficacy(scip, NULL, row), SCIPgetRowMaxCoef(scip, row) / SCIPgetRowMinCoef(scip, row) /
4381  * SCIPgetCutEfficacy(scip, NULL, row));
4382  */
4383  assert(SCIPgetCutEfficacy(scip, NULL, row) > 0.0);
4384  if( ! nlhdlrdata->ignorehighre || SCIPgetRowMaxCoef(scip, row) / SCIPgetRowMinCoef(scip, row) / SCIPgetCutEfficacy(scip, NULL, row) < 1e9 )
4385  {
4386 #ifdef SCIP_DEBUG
4387  SCIPdebugMsg(scip, "adding cut ");
4388  SCIP_CALL( SCIPprintRow(scip, row, NULL) );
4389 #endif
4390 
4391  SCIP_CALL( SCIPaddRow(scip, row, FALSE, &infeasible) );
4392 
4393  if( infeasible )
4394  {
4395  *result = SCIP_CUTOFF;
4396  }
4397  else
4398  {
4399  *result = SCIP_SEPARATED;
4400  nlhdlrdata->cutcoefsum += nlhdlrdata->currentavecutcoef;
4401  nlhdlrdata->monoidalimprovementsum += nlhdlrdata->currentavemonoidalimprovement;
4402  nlhdlrdata->ncutsadded += 1;
4403  nlhdlrdata->densitysum += (SCIP_Real) SCIProwprepGetNVars(rowprep) / (SCIP_Real) SCIPgetNVars(scip);
4404  nlhdlrdata->efficacysum += SCIPgetCutEfficacy(scip, NULL, row);
4405  }
4406  }
4407  else
4408  {
4409  nlhdlrdata->nhighre++;
4410  }
4411  SCIP_CALL( SCIPreleaseRow(scip, &row) );
4412  }
4413 
4414  SCIPfreeRowprep(scip, &rowprep);
4415 
4416  return SCIP_OKAY;
4417 }
4418 
4419 /** nonlinear handler forward propagation callback
4420  *
4421  * This method should solve the problem
4422  * <pre>
4423  * max/min quad expression over box constraints
4424  * </pre>
4425  * However, this problem is difficult so we are satisfied with a proxy.
4426  * Interval arithmetic suffices when no variable appears twice, however this is seldom the case, so we try
4427  * to take care of the dependency problem to some extent:
4428  * Let \f$P_l = \{i : \text{expr}_l \text{expr}_i \,\text{is a bilinear expr}\}\f$.
4429  * 1. partition the quadratic expression as sum of quadratic functions \f$\sum_l q_l\f$
4430  * 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$
4431  * 2. build interval quadratic functions, i.e., \f$a x^2 + b x\f$ where \f$b\f$ is an interval, i.e.,
4432  * \f$a_l \text{expr}_l^2 + [\sum_{i \in P_l} b_{il} \text{expr}_i + c_l] \text{expr}_l\f$
4433  * 3. compute \f$\min/\max \{ a x^2 + b x : x \in [x] \}\f$ for each interval quadratic, i.e.,
4434  * \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$
4435  *
4436  * Notes:
4437  * 1. The \f$l\f$-th quadratic expr (expressions that appear quadratically) is associated with \f$q_l\f$.
4438  * 2. `nlhdlrdata->quadactivities[l]` is the activity of \f$q_l\f$ as computed in the description above.
4439  * 3. The \f$q_l\f$ of a quadratic term might be empty, in which case `nlhdlrdata->quadactivities[l]` is [0,0].\n
4440  * For example, consider \f$x^2 + xy\f$. There are two quadratic expressions, \f$x\f$ and \f$y\f$.
4441  * 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.
4442  * Thus, `nlhdlrdata->quadactivities[1]` is [0,0] in this case.
4443  * The logic is to avoid considering the term \f$xy\f$ twice.
4444  *
4445  * @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$
4446  */
4447 static
4448 SCIP_DECL_NLHDLRINTEVAL(nlhdlrIntevalQuadratic)
4449 { /*lint --e{715}*/
4450  SCIP_EXPR** linexprs;
4451  SCIP_Real* lincoefs;
4452  SCIP_Real constant;
4453  int nquadexprs;
4454  int nlinexprs;
4455 
4456  assert(scip != NULL);
4457  assert(expr != NULL);
4458 
4459  assert(nlhdlrexprdata != NULL);
4460  assert(nlhdlrexprdata->quadactivities != NULL);
4461  assert(nlhdlrexprdata->qexpr == expr);
4462 
4463  SCIPdebugMsg(scip, "Interval evaluation of quadratic expr\n");
4464 
4465  SCIPexprGetQuadraticData(expr, &constant, &nlinexprs, &linexprs, &lincoefs, &nquadexprs, NULL, NULL, NULL);
4466 
4467  /*
4468  * compute activity of linear part, if some linear term has changed
4469  */
4470  {
4471  int i;
4472 
4473  SCIPdebugMsg(scip, "Computing activity of linear part\n");
4474 
4475  SCIPintervalSet(&nlhdlrexprdata->linactivity, constant);
4476  for( i = 0; i < nlinexprs; ++i )
4477  {
4478  SCIP_INTERVAL linterminterval;
4479 
4480  linterminterval = SCIPexprGetActivity(linexprs[i]);
4481  if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, linterminterval) )
4482  {
4483  SCIPdebugMsg(scip, "Activity of linear part is empty due to child %d\n", i);
4484  SCIPintervalSetEmpty(interval);
4485  return SCIP_OKAY;
4486  }
4487  SCIPintervalMulScalar(SCIP_INTERVAL_INFINITY, &linterminterval, linterminterval, lincoefs[i]);
4488  SCIPintervalAdd(SCIP_INTERVAL_INFINITY, &nlhdlrexprdata->linactivity, nlhdlrexprdata->linactivity, linterminterval);
4489  }
4490 
4491  SCIPdebugMsg(scip, "Activity of linear part is [%g, %g]\n", nlhdlrexprdata->linactivity.inf,
4492  nlhdlrexprdata->linactivity.sup);
4493  }
4494 
4495  /*
4496  * compute activity of quadratic part
4497  */
4498  {
4499  int i;
4500 
4501  SCIPdebugMsg(scip, "Computing activity of quadratic part\n");
4502 
4503  nlhdlrexprdata->nneginfinityquadact = 0;
4504  nlhdlrexprdata->nposinfinityquadact = 0;
4505  nlhdlrexprdata->minquadfiniteact = 0.0;
4506  nlhdlrexprdata->maxquadfiniteact = 0.0;
4507  SCIPintervalSet(&nlhdlrexprdata->quadactivity, 0.0);
4508 
4509  for( i = 0; i < nquadexprs; ++i )
4510  {
4511  SCIP_Real quadlb;
4512  SCIP_Real quadub;
4513  SCIP_EXPR* qexpr;
4514  SCIP_Real lincoef;
4515  SCIP_Real sqrcoef;
4516  int nadjbilin;
4517  int* adjbilin;
4518  SCIP_EXPR* sqrexpr;
4519 
4520  SCIPexprGetQuadraticQuadTerm(expr, i, &qexpr, &lincoef, &sqrcoef, &nadjbilin, &adjbilin, &sqrexpr);
4521 
4522  if( !isPropagableTerm(expr, i) )
4523  {
4524  /* term is not propagable, i.e., the exprs involved in term only appear once; thus use the activity of the
4525  * quadratic term directly and not the activity of the exprs involed in the term. See also documentation of
4526  * DETECT
4527  */
4528  SCIP_INTERVAL tmp;
4529 
4530  assert(lincoef == 0.0);
4531 
4532  if( sqrcoef != 0.0 )
4533  {
4534  assert(sqrexpr != NULL);
4535  assert(nadjbilin == 0);
4536 
4537  tmp = SCIPexprGetActivity(sqrexpr);
4539  {
4540  SCIPintervalSetEmpty(interval);
4541  return SCIP_OKAY;
4542  }
4543 
4544  SCIPintervalMulScalar(SCIP_INTERVAL_INFINITY, &tmp, tmp, sqrcoef);
4545  quadlb = tmp.inf;
4546  quadub = tmp.sup;
4547 
4548 #ifdef DEBUG_PROP
4549  SCIPinfoMessage(scip, NULL, "Computing activity for quadratic term %g <expr>, where <expr> is: ", sqrcoef);
4550  SCIP_CALL( SCIPprintExpr(scip, sqrexpr, NULL) );
4551 #endif
4552  }
4553  else
4554  {
4555  SCIP_EXPR* expr1;
4556  SCIP_EXPR* prodexpr;
4557  SCIP_Real prodcoef;
4558 
4559  assert(nadjbilin == 1);
4560  SCIPexprGetQuadraticBilinTerm(expr, adjbilin[0], &expr1, NULL, &prodcoef, NULL, &prodexpr);
4561 
4562  if( expr1 == qexpr )
4563  {
4564  /* the quadratic expression expr1 appears only as expr1 * expr2, so its 'q' is expr1 * expr2 */
4565  tmp = SCIPexprGetActivity(prodexpr);
4567  {
4568  SCIPintervalSetEmpty(interval);
4569  return SCIP_OKAY;
4570  }
4571 
4572  SCIPintervalMulScalar(SCIP_INTERVAL_INFINITY, &tmp, tmp, prodcoef);
4573  quadlb = tmp.inf;
4574  quadub = tmp.sup;
4575 
4576 #ifdef DEBUG_PROP
4577  SCIPinfoMessage(scip, NULL, "Computing activity for quadratic term %g <expr>, where <expr> is: ", prodcoef);
4578  SCIP_CALL( SCIPprintExpr(scip, prodexpr, NULL) );
4579 #endif
4580  }
4581  else
4582  {
4583  /* the quadratic expression expr1 appears as expr2 * expr1, thus its 'q' is empty, see also the Notes
4584  * in the documentation of the function
4585  */
4586  SCIPintervalSet(&nlhdlrexprdata->quadactivities[i], 0.0);
4587  continue;
4588  }
4589  }
4590  }
4591  else
4592  {
4593  int j;
4594  SCIP_INTERVAL b;
4595 
4596  SCIPexprGetQuadraticQuadTerm(expr, i, &qexpr, &lincoef, &sqrcoef, &nadjbilin, &adjbilin, NULL);
4597 
4599  {
4600  SCIPintervalSetEmpty(interval);
4601  return SCIP_OKAY;
4602  }
4603 
4604  /* b = [c_l] */
4605  SCIPintervalSet(&b, lincoef);
4606 #ifdef DEBUG_PROP
4607  SCIPinfoMessage(scip, NULL, "b := %g\n", lincoef);
4608 #endif
4609  for( j = 0; j < nadjbilin; ++j )
4610  {
4611  SCIP_INTERVAL bterm;
4612  SCIP_EXPR* expr1;
4613  SCIP_EXPR* expr2;
4614  SCIP_Real bilincoef;
4615 
4616  SCIPexprGetQuadraticBilinTerm(expr, adjbilin[j], &expr1, &expr2, &bilincoef, NULL, NULL);
4617 
4618  if( expr1 != qexpr )
4619  continue;
4620 
4621  bterm = SCIPexprGetActivity(expr2);
4623  {
4624  SCIPintervalSetEmpty(interval);
4625  return SCIP_OKAY;
4626  }
4627 
4628  /* b += [b_jl * expr_j] for j \in P_l */
4629  SCIPintervalMulScalar(SCIP_INTERVAL_INFINITY, &bterm, bterm, bilincoef);
4630  SCIPintervalAdd(SCIP_INTERVAL_INFINITY, &b, b, bterm);
4631 
4632 #ifdef DEBUG_PROP
4633  SCIPinfoMessage(scip, NULL, "b += %g * [expr2], where <expr2> is: ", bilincoef);
4634  SCIP_CALL( SCIPprintExpr(scip, expr2, NULL) );
4635  SCIPinfoMessage(scip, NULL, " [%g,%g]\n", SCIPexprGetActivity(expr2).inf, SCIPexprGetActivity(expr2).sup);
4636 #endif
4637  }
4638 
4639  /* 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 */
4641  SCIPexprGetActivity(qexpr));
4642 
4643  /* TODO: implement SCIPintervalQuadLowerBound */
4644  {
4645  SCIP_INTERVAL minusb;
4647 
4648  quadlb = -SCIPintervalQuadUpperBound(SCIP_INTERVAL_INFINITY, -sqrcoef, minusb,
4649  SCIPexprGetActivity(qexpr));
4650  }
4651 
4652 #ifdef DEBUG_PROP
4653  SCIPinfoMessage(scip, NULL, "Computing activity for quadratic term %g <expr>^2 + [%g,%g] <expr>, where <expr> is: ", sqrcoef, b.inf, b.sup);
4654  SCIP_CALL( SCIPprintExpr(scip, qexpr, NULL) );
4655 #endif
4656  }
4657 #ifdef DEBUG_PROP
4658  SCIPinfoMessage(scip, NULL, " -> [%g, %g]\n", quadlb, quadub);
4659 #endif
4660 
4661  SCIPintervalSetBounds(&nlhdlrexprdata->quadactivities[i], quadlb, quadub);
4662  SCIPintervalAdd(SCIP_INTERVAL_INFINITY, &nlhdlrexprdata->quadactivity, nlhdlrexprdata->quadactivity, nlhdlrexprdata->quadactivities[i]);
4663 
4664  /* get number of +/-infinity contributions and compute finite activity */
4665  if( quadlb <= -SCIP_INTERVAL_INFINITY )
4666  nlhdlrexprdata->nneginfinityquadact++;
4667  else
4668  {
4669  SCIP_ROUNDMODE roundmode;
4670 
4671  roundmode = SCIPintervalGetRoundingMode();
4673 
4674  nlhdlrexprdata->minquadfiniteact += quadlb;
4675 
4676  SCIPintervalSetRoundingMode(roundmode);
4677  }
4678  if( quadub >= SCIP_INTERVAL_INFINITY )
4679  nlhdlrexprdata->nposinfinityquadact++;
4680  else
4681  {
4682  SCIP_ROUNDMODE roundmode;
4683 
4684  roundmode = SCIPintervalGetRoundingMode();
4686 
4687  nlhdlrexprdata->maxquadfiniteact += quadub;
4688 
4689  SCIPintervalSetRoundingMode(roundmode);
4690  }
4691  }
4692 
4693  SCIPdebugMsg(scip, "Activity of quadratic part is [%g, %g]\n", nlhdlrexprdata->quadactivity.inf, nlhdlrexprdata->quadactivity.sup);
4694  }
4695 
4696  /* interval evaluation is linear activity + quadactivity */
4697  SCIPintervalAdd(SCIP_INTERVAL_INFINITY, interval, nlhdlrexprdata->linactivity, nlhdlrexprdata->quadactivity);
4698 
4699  nlhdlrexprdata->activitiestag = SCIPgetCurBoundsTagNonlinear(SCIPfindConshdlr(scip, "nonlinear"));
4700 
4701  return SCIP_OKAY;
4702 }
4703 
4704 /** nonlinear handler reverse propagation callback
4705  *
4706  * @note the implemented technique is a proxy for solving the problem min/max{ x_i : quad expr in [quad expr] }
4707  * and as such can be improved.
4708  */
4709 static
4710 SCIP_DECL_NLHDLRREVERSEPROP(nlhdlrReversepropQuadratic)
4711 { /*lint --e{715}*/
4712  SCIP_EXPR** linexprs;
4713  SCIP_EXPR** bilinexprs; /* TODO: should this be stored in the nlhdlr expr data? */
4714  SCIP_Real* bilincoefs;
4715  SCIP_Real* lincoefs;
4716  SCIP_Real constant;
4717  int nquadexprs;
4718  int nlinexprs;
4719 
4720  SCIP_INTERVAL rhs;
4721  SCIP_INTERVAL quadactivity;
4722  int i;
4723 
4724  SCIPdebugMsg(scip, "Reverse propagation of quadratic expr given bounds = [%g,%g]\n", bounds.inf, bounds.sup);
4725 
4726  assert(scip != NULL);
4727  assert(expr != NULL);
4728  assert(infeasible != NULL);
4729  assert(nreductions != NULL);
4730  assert(nlhdlrexprdata != NULL);
4731  assert(nlhdlrexprdata->quadactivities != NULL);
4732  assert(nlhdlrexprdata->qexpr == expr);
4733 
4734  *nreductions = 0;
4735 
4736  /* not possible to conclude finite bounds if the interval of the expression is [-inf,inf] */
4738  {
4739  SCIPdebugMsg(scip, "expr's range is R -> cannot reverse propagate\n");
4740  return SCIP_OKAY;
4741  }
4742 
4743  /* ensure that partial activities as stored in nlhdlrexprdata are uptodate
4744  * if the activity stored in expr is more recent than the partial activities stored in this nlhdlrexprdata,
4745  * then we should reevaluate the partial activities
4746  */
4747  if( SCIPexprGetActivityTag(expr) > nlhdlrexprdata->activitiestag )
4748  {
4749  SCIP_CALL( nlhdlrIntevalQuadratic(scip, nlhdlr, expr, nlhdlrexprdata, &quadactivity, NULL, NULL) );
4750  }
4751 
4752  SCIPexprGetQuadraticData(expr, &constant, &nlinexprs, &linexprs, &lincoefs, &nquadexprs, NULL, NULL, NULL);
4753 
4754  /* propagate linear part in rhs = expr's interval - quadratic activity; first, reconstruct the quadratic activity */
4755  SCIPintervalSetBounds(&quadactivity,
4756  nlhdlrexprdata->nneginfinityquadact > 0 ? -SCIP_INTERVAL_INFINITY : nlhdlrexprdata->minquadfiniteact,
4757  nlhdlrexprdata->nposinfinityquadact > 0 ? SCIP_INTERVAL_INFINITY : nlhdlrexprdata->maxquadfiniteact);
4758 
4759  SCIPintervalSub(SCIP_INTERVAL_INFINITY, &rhs, bounds, quadactivity);
4760 
4761  SCIP_CALL( reversePropagateLinearExpr(scip, linexprs, nlinexprs, lincoefs, constant, rhs, infeasible, nreductions) );
4762 
4763  /* stop if we find infeasibility */
4764  if( *infeasible )
4765  return SCIP_OKAY;
4766 
4767  /* propagate quadratic part in expr's interval - linear activity, where linear activity was computed in INTEVAL.
4768  * The idea is basically to write interval quadratics for each expr and then solve for expr.
4769  *
4770  * One way of achieving this is:
4771  * - for each expression expr_i, write the quadratic expression as a_i expr^2_i + expr_i ( \sum_{j \in J_i} b_ij
4772  * expr_j + c_i ) + quadratic expression in expr_k for k \neq i
4773  * - 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
4774  * bilinear expression expr_i expr_j appears
4775  * - use some technique (like the one in nlhdlrIntevalQuadratic), to evaluate the activity of rest_i = [quadratic
4776  * expression in expr_k for k \neq i].
4777  * - solve a_i expr_i^2 + b expr_i \in rhs_i := [expr activity] - rest_i
4778  *
4779  * However, this might be expensive, especially computing rest_i. Hence, we implement a simpler version.
4780  * - we use the same partition as in nlhdlrIntevalQuadratic for the bilinear terms. This way, b = [\sum_{j \in P_i}
4781  * 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
4782  * - 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
4783  * \in P_k} b_jk expr_j + c_k] expr_k. The intervals [\min q_k, \max q_k] were already computed in
4784  * nlhdlrIntevalQuadratic, so we just reuse them.
4785  *
4786  * A downside of the above is that we might not deduce any bounds for variables that appear less often. For example,
4787  * 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
4788  * first parenthesis is interpreted as a function of x, while the second one as a function of z.
4789  * 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 -
4790  * x and propagate the y + z).
4791  * In general, after reverse propagating expr_i, we consider
4792  * \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,
4793  * compute an interval for the right hand side (see computeRangeForBilinearProp) and use that to propagate the
4794  * linear sum on the left hand side.
4795  *
4796  * Note: this last step generalizes a technique that appeared in the classic cons_quadratic.
4797  * The idea of that technique was to borrow a bilinear term expr_k expr_l when propagating expr_l and the quadratic
4798  * function for expr_k was simple enough.
4799  * Since in P_l we only consider the indices of expressions that appear multiplying expr_l as _second_ factor, we
4800  * would lose the bilinear terms expr_k * expr_l, which contributes to the dependency problem.
4801  * 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
4802  * rather quadactivities[k] (= max/min of a_k expr_k^2 + expr_k * [c_k + sum_i \in P_k b_ki expr_i]).
4803  * Thus, we _cannot_ just substract [b_kl * expr_k * expr_l] from rest_i.
4804  * But, if expr_k only appears as expr_k * expr_l, then quadactivities[k] = [b_kl * expr_k * expr_l]. So this
4805  * case was handled in old cons_quadratic.
4806  *
4807  *
4808  * TODO: handle simple cases
4809  * TODO: identify early when there is nothing to be gain
4810  */
4811  SCIPintervalSub(SCIP_INTERVAL_INFINITY, &rhs, bounds, nlhdlrexprdata->linactivity);
4812  SCIP_CALL( SCIPallocBufferArray(scip, &bilinexprs, nquadexprs) );
4813  SCIP_CALL( SCIPallocBufferArray(scip, &bilincoefs, nquadexprs) );
4814 
4815  for( i = 0; i < nquadexprs; ++i )
4816  {
4817  SCIP_INTERVAL rhs_i;
4818  SCIP_INTERVAL rest_i;
4819  SCIP_EXPR* qexpr;
4820  SCIP_Real lincoef;
4821  SCIP_Real sqrcoef;
4822  int nadjbilin;
4823  int* adjbilin;
4824  SCIP_EXPR* sqrexpr;
4825 
4826  SCIPexprGetQuadraticQuadTerm(expr, i, &qexpr, &lincoef, &sqrcoef, &nadjbilin, &adjbilin, &sqrexpr);
4827 
4828  /* rhs_i = rhs - rest_i.
4829  * to compute rest_i = [\sum_{k \neq i} q_k] we just have to substract
4830  * the activity of q_i from quadactivity; however, care must be taken about infinities;
4831  * if [q_i].sup = +infinity and there is = 1 contributing +infinity -> rest_i.sup = maxquadfiniteact
4832  * if [q_i].sup = +infinity and there is > 1 contributing +infinity -> rest_i.sup = +infinity
4833  * if [q_i].sup = finite and there is > 0 contributing +infinity -> rest_i.sup = +infinity
4834  * if [q_i].sup = finite and there is = 0 contributing +infinity -> rest_i.sup = maxquadfiniteact - [q_i].sup
4835  *
4836  * the same holds when replacing sup with inf, + with - and max(quadfiniteact) with min(...)
4837  */
4838  /* compute rest_i.sup */
4839  if( SCIPintervalGetSup(nlhdlrexprdata->quadactivities[i]) < SCIP_INTERVAL_INFINITY &&
4840  nlhdlrexprdata->nposinfinityquadact == 0 )
4841  {
4842  SCIP_ROUNDMODE roundmode;
4843 
4844  roundmode = SCIPintervalGetRoundingMode();
4846  rest_i.sup = nlhdlrexprdata->maxquadfiniteact - SCIPintervalGetSup(nlhdlrexprdata->quadactivities[i]);
4847 
4848  SCIPintervalSetRoundingMode(roundmode);
4849  }
4850  else if( SCIPintervalGetSup(nlhdlrexprdata->quadactivities[i]) >= SCIP_INTERVAL_INFINITY &&
4851  nlhdlrexprdata->nposinfinityquadact == 1 )
4852  rest_i.sup = nlhdlrexprdata->maxquadfiniteact;
4853  else
4854  rest_i.sup = SCIP_INTERVAL_INFINITY;
4855 
4856  /* compute rest_i.inf */
4857  if( SCIPintervalGetInf(nlhdlrexprdata->quadactivities[i]) > -SCIP_INTERVAL_INFINITY &&
4858  nlhdlrexprdata->nneginfinityquadact == 0 )
4859  {
4860  SCIP_ROUNDMODE roundmode;
4861 
4862  roundmode = SCIPintervalGetRoundingMode();
4864  rest_i.inf = nlhdlrexprdata->minquadfiniteact - SCIPintervalGetInf(nlhdlrexprdata->quadactivities[i]);
4865 
4866  SCIPintervalSetRoundingMode(roundmode);
4867  }
4868  else if( SCIPintervalGetInf(nlhdlrexprdata->quadactivities[i]) <= -SCIP_INTERVAL_INFINITY &&
4869  nlhdlrexprdata->nneginfinityquadact == 1 )
4870  rest_i.inf = nlhdlrexprdata->minquadfiniteact;
4871  else
4872  rest_i.inf = -SCIP_INTERVAL_INFINITY;
4873 
4874 #ifdef SCIP_DISABLED_CODE /* I (SV) added the following in cons_quadratic to fix/workaround some bug. Maybe we'll need this here, too? */
4875  /* FIXME in theory, rest_i should not be empty here
4876  * what we tried to do here is to remove the contribution of the i'th bilinear term (=bilinterm) to [minquadactivity,maxquadactivity] from rhs
4877  * however, quadactivity is computed differently (as x*(a1*y1+...+an*yn)) than q_i (a*ak*yk) and since interval arithmetics do overestimation,
4878  * it can happen that q_i is actually slightly larger than quadactivity, which results in rest_i being (slightly) empty
4879  * a proper fix could be to compute the quadactivity also as x*a1*y1+...+x*an*yn if sqrcoef=0, but due to taking
4880  * also infinite bounds into account, this complicates the code even further
4881  * instead, I'll just work around this by turning an empty rest_i into a small non-empty one
4882  */
4884  {
4885  assert(SCIPisSumRelEQ(scip, rest_i.inf, rest_i.sup));
4886  SCIPswapReals(&rest_i.inf, &rest_i.sup);
4887  }
4888 #endif
4889  assert(!SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, rest_i));
4890 
4891  /* compute rhs_i */
4892  SCIPintervalSub(SCIP_INTERVAL_INFINITY, &rhs_i, rhs, rest_i);
4893 
4895  continue;
4896 
4897  /* try to propagate */
4898  if( !isPropagableTerm(expr, i) )
4899  {
4900  assert(lincoef == 0.0);
4901 
4902  if( sqrcoef != 0.0 )
4903  {
4904  assert(sqrexpr != NULL);
4905  assert(nadjbilin == 0);
4906 
4907  /* solve sqrcoef sqrexpr in rhs_i */
4908  SCIP_CALL( propagateBoundsLinExpr(scip, sqrexpr, sqrcoef, rhs_i, infeasible, nreductions) );
4909  }
4910  else
4911  {
4912  /* qexpr only appears in a term of the form qexpr * other_expr (or other_expr * qexpr); we only care about
4913  * getting bounds for the product, thus we will compute these bounds when qexpr appears as qexpr *
4914  * other_expr; note that if it appears as other_expr * qexpr, then when we process other_expr bounds for the
4915  * product will be computed
4916  * TODO: we can actually avoid computing rhs_i in the case that qexpr is not propagable and it appears as
4917  * other_expr * qexpr
4918  */
4919  SCIP_EXPR* expr1;
4920  SCIP_EXPR* prodexpr;
4921  SCIP_Real prodcoef;
4922 
4923  assert(nadjbilin == 1);
4924  SCIPexprGetQuadraticBilinTerm(expr, adjbilin[0], &expr1, NULL, &prodcoef, NULL, &prodexpr);
4925 
4926  if( expr1 == qexpr )
4927  {
4928  /* solve prodcoef prodexpr in rhs_i */
4929  SCIP_CALL( propagateBoundsLinExpr(scip, prodexpr, prodcoef, rhs_i, infeasible, nreductions) );
4930  }
4931  }
4932  }
4933  else
4934  {
4935  SCIP_INTERVAL b;
4936  SCIP_EXPR* expr1 = NULL;
4937  SCIP_EXPR* expr2 = NULL;
4938  SCIP_Real bilincoef = 0.0;
4939  int nbilin = 0;
4940  int pos2 = 0;
4941  int j;
4942 
4943  /* set b to [c_l] */
4944  SCIPintervalSet(&b, lincoef);
4945 
4946  /* add [\sum_{j \in P_l} b_lj expr_j + c_l] into b */
4947  for( j = 0; j < nadjbilin; ++j )
4948  {
4949  SCIP_INTERVAL bterm;
4950  SCIP_INTERVAL expr2bounds;
4951 
4952  SCIPexprGetQuadraticBilinTerm(expr, adjbilin[j], &expr1, &expr2, &bilincoef, &pos2, NULL);
4953 
4954  if( expr1 != qexpr )
4955  continue;
4956 
4957  expr2bounds = SCIPgetExprBoundsNonlinear(scip, expr2);
4958  if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, expr2bounds) )
4959  {
4960  *infeasible = TRUE;
4961  break;
4962  }
4963 
4964  /* b += [b_lj * expr_j] for j \in P_l */
4965  SCIPintervalMulScalar(SCIP_INTERVAL_INFINITY, &bterm, expr2bounds, bilincoef);
4966  SCIPintervalAdd(SCIP_INTERVAL_INFINITY, &b, b, bterm);
4967 
4968  /* remember b_lj and expr_j to propagate them too */
4969  bilinexprs[nbilin] = expr2;
4970  bilincoefs[nbilin] = bilincoef;
4971  nbilin++;
4972  }
4973 
4974  if( !*infeasible )
4975  {
4976  /* solve a_i expr_i^2 + b expr_i in rhs_i */
4977  SCIP_CALL( propagateBoundsQuadExpr(scip, qexpr, sqrcoef, b, rhs_i, infeasible, nreductions) );
4978  }
4979 
4980  if( nbilin > 0 && !*infeasible )
4981  {
4982  /* if 0 is not in [expr_i], then propagate bilincoefs^T bilinexpr in rhs_i/expr_i - a_i expr_i - c_i */
4983  SCIP_INTERVAL bilinrhs;
4984  SCIP_INTERVAL qexprbounds;
4985 
4986  qexprbounds = SCIPgetExprBoundsNonlinear(scip, qexpr);
4987  if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, qexprbounds) )
4988  {
4989  *infeasible = TRUE;
4990  }
4991  else
4992  {
4993  /* compute bilinrhs := [rhs_i/expr_i - a_i expr_i] */
4994  computeRangeForBilinearProp(qexprbounds, sqrcoef, rhs_i, &bilinrhs);
4995 
4997  {
4998  int nreds;
4999 
5000  /* propagate \sum_{j \in P_i} b_ij expr_j + c_i in bilinrhs */
5001  SCIP_CALL( reversePropagateLinearExpr(scip, bilinexprs, nbilin, bilincoefs, lincoef, bilinrhs,
5002  infeasible, &nreds) );
5003 
5004  /* TODO FIXME: we are overestimating the number of reductions: an expr might be tightened many times! */
5005  *nreductions += nreds;
5006  }
5007  }
5008  }
5009  }
5010 
5011  /* stop if we find infeasibility */
5012  if( *infeasible )
5013  break;
5014  }
5015 
5016  SCIPfreeBufferArray(scip, &bilincoefs);
5017  SCIPfreeBufferArray(scip, &bilinexprs);
5018 
5019  return SCIP_OKAY;
5020 }
5021 
5022 /** callback to free data of handler */
5023 static
5024 SCIP_DECL_NLHDLRFREEHDLRDATA(nlhdlrFreehdlrdataQuadratic)
5025 { /*lint --e{715}*/
5026  assert(nlhdlrdata != NULL);
5027 
5028  SCIPfreeBlockMemory(scip, nlhdlrdata);
5029 
5030  return SCIP_OKAY;
5031 }
5032 
5033 /** nonlinear handler copy callback */
5034 static
5035 SCIP_DECL_NLHDLRCOPYHDLR(nlhdlrCopyhdlrQuadratic)
5036 { /*lint --e{715}*/
5037  assert(targetscip != NULL);
5038  assert(sourcenlhdlr != NULL);
5039  assert(strcmp(SCIPnlhdlrGetName(sourcenlhdlr), NLHDLR_NAME) == 0);
5040 
5041  SCIP_CALL( SCIPincludeNlhdlrQuadratic(targetscip) );
5042 
5043  return SCIP_OKAY;
5044 }
5045 
5046 /** includes quadratic nonlinear handler in nonlinear constraint handler */
5048  SCIP* scip /**< SCIP data structure */
5049  )
5050 {
5051  SCIP_NLHDLRDATA* nlhdlrdata;
5052  SCIP_NLHDLR* nlhdlr;
5053 
5054  assert(scip != NULL);
5055 
5056  /* create nonlinear handler specific data */
5057  SCIP_CALL( SCIPallocBlockMemory(scip, &nlhdlrdata) );
5058  BMSclearMemory(nlhdlrdata);
5059 
5061  NLHDLR_ENFOPRIORITY, nlhdlrDetectQuadratic, nlhdlrEvalauxQuadratic, nlhdlrdata) );
5062 
5063  SCIPnlhdlrSetCopyHdlr(nlhdlr, nlhdlrCopyhdlrQuadratic);
5064  SCIPnlhdlrSetFreeHdlrData(nlhdlr, nlhdlrFreehdlrdataQuadratic);
5065  SCIPnlhdlrSetFreeExprData(nlhdlr, nlhdlrFreeexprdataQuadratic);
5066  SCIPnlhdlrSetSepa(nlhdlr, NULL, nlhdlrEnfoQuadratic, NULL, NULL);
5067  SCIPnlhdlrSetProp(nlhdlr, nlhdlrIntevalQuadratic, nlhdlrReversepropQuadratic);
5068 
5069  /* parameters */
5070  SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" NLHDLR_NAME "/useintersectioncuts",
5071  "whether to use intersection cuts for quadratic constraints to separate",
5072  &nlhdlrdata->useintersectioncuts, FALSE, DEFAULT_USEINTERCUTS, NULL, NULL) );
5073 
5074  SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" NLHDLR_NAME "/usestrengthening",
5075  "whether the strengthening should be used",
5076  &nlhdlrdata->usestrengthening, FALSE, DEFAULT_USESTRENGTH, NULL, NULL) );
5077 
5078  SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" NLHDLR_NAME "/usemonoidal",
5079  "whether monoidal strengthening should be used",
5080  &nlhdlrdata->usemonoidal, FALSE, DEFAULT_USEMONOIDAL, NULL, NULL) );
5081 
5082  SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" NLHDLR_NAME "/useminrep",
5083  "whether the minimal representation of the S-free set should be used (instead of the gauge)",
5084  &nlhdlrdata->useminrep, FALSE, DEFAULT_USEMINREP, NULL, NULL) );
5085 
5086  SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" NLHDLR_NAME "/useboundsasrays",
5087  "use bounds of variables in quadratic as rays for intersection cuts",
5088  &nlhdlrdata->useboundsasrays, FALSE, DEFAULT_USEBOUNDS, NULL, NULL) );
5089 
5090  SCIP_CALL( SCIPaddIntParam(scip, "nlhdlr/" NLHDLR_NAME "/ncutslimit",
5091  "limit for number of cuts generated consecutively",
5092  &nlhdlrdata->ncutslimit, FALSE, DEFAULT_NCUTS, 0, INT_MAX, NULL, NULL) );
5093 
5094  SCIP_CALL( SCIPaddIntParam(scip, "nlhdlr/" NLHDLR_NAME "/ncutslimitroot",
5095  "limit for number of cuts generated at root node",
5096  &nlhdlrdata->ncutslimitroot, FALSE, DEFAULT_NCUTSROOT, 0, INT_MAX, NULL, NULL) );
5097 
5098  SCIP_CALL( SCIPaddIntParam(scip, "nlhdlr/" NLHDLR_NAME "/maxrank",
5099  "maximal rank a slackvar can have",
5100  &nlhdlrdata->maxrank, FALSE, INT_MAX, 0, INT_MAX, NULL, NULL) );
5101 
5102  SCIP_CALL( SCIPaddRealParam(scip, "nlhdlr/" NLHDLR_NAME "/mincutviolation",
5103  "minimal cut violation the generated cuts must fulfill to be added to the LP",
5104  &nlhdlrdata->mincutviolation, FALSE, 1e-4, 0.0, SCIPinfinity(scip), NULL, NULL) );
5105 
5106  SCIP_CALL( SCIPaddRealParam(scip, "nlhdlr/" NLHDLR_NAME "/minviolation",
5107  "minimal violation the constraint must fulfill such that a cut is generated",
5108  &nlhdlrdata->minviolation, FALSE, INTERCUTS_MINVIOL, 0.0, SCIPinfinity(scip), NULL, NULL) );
5109 
5110  SCIP_CALL( SCIPaddIntParam(scip, "nlhdlr/" NLHDLR_NAME "/atwhichnodes",
5111  "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",
5112  &nlhdlrdata->atwhichnodes, FALSE, 1, -1, INT_MAX, NULL, NULL) );
5113 
5114  SCIP_CALL( SCIPaddIntParam(scip, "nlhdlr/" NLHDLR_NAME "/nstrengthlimit",
5115  "limit for number of rays we do the strengthening for",
5116  &nlhdlrdata->nstrengthlimit, FALSE, INT_MAX, 0, INT_MAX, NULL, NULL) );
5117 
5118  SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" NLHDLR_NAME "/sparsifycuts",
5119  "should we try to sparisfy the intersection cut?",
5120  &nlhdlrdata->sparsifycuts, FALSE, FALSE, NULL, NULL) );
5121 
5122  SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" NLHDLR_NAME "/ignorebadrayrestriction",
5123  "should cut be generated even with bad numerics when restricting to ray?",
5124  &nlhdlrdata->ignorebadrayrestriction, FALSE, TRUE, NULL, NULL) );
5125 
5126  SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" NLHDLR_NAME "/ignorenhighre",
5127  "should cut be added even when range / efficacy is large?",
5128  &nlhdlrdata->ignorehighre, FALSE, TRUE, NULL, NULL) );
5129 
5130  SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" NLHDLR_NAME "/trackmore",
5131  "for monoidal strengthening, should we track more statistics (more expensive)?",
5132  &nlhdlrdata->trackmore, FALSE, FALSE, NULL, NULL) );
5133 
5134  /* statistic table */
5135  assert(SCIPfindTable(scip, TABLE_NAME_QUADRATIC) == NULL);
5137  NULL, NULL, NULL, NULL, NULL, NULL, tableOutputQuadratic,
5139  return SCIP_OKAY;
5140 }
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:7500
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:7490
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:453
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 sidefactor, SCIP_Real *a, SCIP_Real *b, SCIP_Real *c)
SCIP_VARTYPE SCIPvarGetType(SCIP_VAR *var)
Definition: var.c:17585
struct SCIP_NlhdlrData SCIP_NLHDLRDATA
Definition: type_nlhdlr.h:452
SCIP_Bool SCIPisZero(SCIP *scip, SCIP_Real val)
sum expression handler
int rayssize
int SCIPgetNLPCols(SCIP *scip)
Definition: scip_lp.c:527
SCIP_Real SCIPvarGetUbLocal(SCIP_VAR *var)
Definition: var.c:18145
SCIP_RETCODE SCIPaddRowprepTerm(SCIP *scip, SCIP_ROWPREP *rowprep, SCIP_VAR *var, SCIP_Real coef)
Definition: misc_rowprep.c:913
#define DEFAULT_USEMONOIDAL
SCIP_RETCODE SCIPcleanupRowprep(SCIP *scip, SCIP_ROWPREP *rowprep, SCIP_SOL *sol, SCIP_Real minviol, SCIP_Real *viol, SCIP_Bool *success)
SCIP_VAR * SCIPgetExprAuxVarNonlinear(SCIP_EXPR *expr)
#define BMSclearMemoryArray(ptr, num)
Definition: memory.h:130
SCIP_RETCODE SCIPgetLPRowsData(SCIP *scip, SCIP_ROW ***rows, int *nrows)
Definition: scip_lp.c:570
static void constructLPPos2ConsPosMap(SCIP_NLHDLREXPRDATA *nlhdlrexprdata, SCIP_VAR *auxvar, int *map)
static SCIP_RETCODE createBoundRays(SCIP *scip, RAYS **rays, int size)
#define BINSEARCH_MAXITERS
#define SCIPallocClearBlockMemory(scip, ptr)
Definition: scip_mem.h:91
static int countBasicVars(SCIP_NLHDLREXPRDATA *nlhdlrexprdata, SCIP_VAR *auxvar, SCIP_Bool *nozerostat)
public functions of nonlinear handlers of nonlinear constraints
#define DEFAULT_USEBOUNDS
SCIP_Longint SCIPgetNLPs(SCIP *scip)
SCIP_INTERVAL SCIPgetExprBoundsNonlinear(SCIP *scip, SCIP_EXPR *expr)
int SCIPcolGetLPPos(SCIP_COL *col)
Definition: lp.c:17093
SCIP_Real SCIPgetLhsNonlinear(SCIP_CONS *cons)
#define SCIP_NLHDLR_METHOD_SEPABOTH
Definition: type_nlhdlr.h:53
static SCIP_RETCODE computeRestrictionToLine(SCIP *scip, SCIP_NLHDLREXPRDATA *nlhdlrexprdata, SCIP_Real sidefactor, SCIP_Real *raycoefs, int *rayidx, int raynnonz, SCIP_Real *vb, SCIP_Real *vzlp, SCIP_Real kappa, SCIP_Real *apex, SCIP_Real *coefs2, SCIP_Bool *success)
nonlinear handler to handle quadratic expressions
static void computeVApexAndVRay(SCIP_NLHDLREXPRDATA *nlhdlrexprdata, SCIP_Real *apex, SCIP_Real *raycoefs, int *rayidx, int raynnonz, SCIP_Real *vapex, SCIP_Real *vray)
SCIP_Real SCIPgetSolVal(SCIP *scip, SCIP_SOL *sol, SCIP_VAR *var)
Definition: scip_sol.c:1217
static SCIP_Real isCase4a(SCIP_Real tsol, SCIP_Real *coefs4a, SCIP_Real *coefscondition)
#define SCIP_NLHDLR_METHOD_SEPAABOVE
Definition: type_nlhdlr.h:52
static SCIP_RETCODE propagateBoundsLinExpr(SCIP *scip, SCIP_EXPR *expr, SCIP_Real b, SCIP_INTERVAL rhs, SCIP_Bool *infeasible, int *nreductions)
SCIP_RETCODE SCIPaddRealParam(SCIP *scip, const char *name, const char *desc, SCIP_Real *valueptr, SCIP_Bool isadvanced, SCIP_Real defaultvalue, SCIP_Real minvalue, SCIP_Real maxvalue, SCIP_DECL_PARAMCHGD((*paramchgd)), SCIP_PARAMDATA *paramdata)
Definition: scip_param.c:139
#define ABS(x)
Definition: def.h:235
preparation of a linear inequality to become a SCIP_ROW
SCIP_RETCODE SCIPaddBoolParam(SCIP *scip, const char *name, const char *desc, SCIP_Bool *valueptr, SCIP_Bool isadvanced, SCIP_Bool defaultvalue, SCIP_DECL_PARAMCHGD((*paramchgd)), SCIP_PARAMDATA *paramdata)
Definition: scip_param.c:57
#define DEFAULT_NCUTS
SCIP_Real SCIPgetRowActivity(SCIP *scip, SCIP_ROW *row)
Definition: scip_lp.c:2104
int * raysbegin
SCIP_VAR ** SCIProwprepGetVars(SCIP_ROWPREP *rowprep)
Definition: misc_rowprep.c:639
SCIP_RETCODE SCIPcreateSol(SCIP *scip, SCIP_SOL **sol, SCIP_HEUR *heur)
Definition: scip_sol.c:184
#define SCIPreallocBufferArray(scip, ptr, num)
Definition: scip_mem.h:128
void SCIPnlhdlrSetCopyHdlr(SCIP_NLHDLR *nlhdlr, SCIP_DECL_NLHDLRCOPYHDLR((*copy)))
Definition: nlhdlr.c:76