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