Scippy

SCIP

Solving Constraint Integer Programs

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