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