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