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