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 sepa_eccuts.c 26 * @ingroup DEFPLUGINS_SEPA 27 * @brief edge concave cut separator 28 * @author Benjamin Mueller 29 */ 30 31 /**@todo only count number of fixed variables in the edge concave terms */ 32 /**@todo only add nonlinear row aggregations where at least ...% of the variables (bilinear terms?) are in edge concave 33 * terms */ 34 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/ 35 36 #include <assert.h> 37 #include <string.h> 38 39 #include "scip/scipdefplugins.h" 40 #include "scip/sepa_eccuts.h" 41 #include "scip/cons_xor.h" 42 #include "scip/nlp.h" 43 #include "tclique/tclique.h" 44 45 #define SEPA_NAME "eccuts" 46 #define SEPA_DESC "separator for edge-concave functions" 47 #define SEPA_PRIORITY -13000 48 #define SEPA_FREQ -1 49 #define SEPA_MAXBOUNDDIST 1.0 50 #define SEPA_USESSUBSCIP FALSE /**< does the separator use a secondary SCIP instance? */ 51 #define SEPA_DELAY FALSE /**< should separation method be delayed, if other separators found cuts? */ 52 53 #define CLIQUE_MAXFIRSTNODEWEIGHT 1000 /**< maximum weight of branching nodes in level 0; 0 if not used for cliques 54 * with at least one fractional node) */ 55 #define CLIQUE_MINWEIGHT 0 /**< lower bound for weight of generated cliques */ 56 #define CLIQUE_MAXNTREENODES 10000 /**< maximal number of nodes of b&b tree */ 57 #define CLIQUE_BACKTRACKFREQ 10000 /**< frequency to backtrack to first level of tree (0: no premature backtracking) */ 58 59 #define DEFAULT_DYNAMICCUTS TRUE /**< should generated cuts be removed from the LP if they are no longer tight? */ 60 #define DEFAULT_MAXROUNDS 10 /**< maximal number of separation rounds per node (-1: unlimited) */ 61 #define DEFAULT_MAXROUNDSROOT 250 /**< maximal number of separation rounds in the root node (-1: unlimited) */ 62 #define DEFAULT_MAXDEPTH -1 /**< maximal depth at which the separator is applied */ 63 #define DEFAULT_MAXSEPACUTS 10 /**< maximal number of e.c. cuts separated per separation round */ 64 #define DEFAULT_MAXSEPACUTSROOT 50 /**< maximal number of e.c. cuts separated per separation round in root node */ 65 #define DEFAULT_CUTMAXRANGE 1e+7 /**< maximal coefficient range of a cut (maximal coefficient divided by minimal 66 * coefficient) in order to be added to LP relaxation */ 67 #define DEFAULT_MINVIOLATION 0.3 /**< minimal violation of an e.c. cut to be separated */ 68 #define DEFAULT_MINAGGRSIZE 3 /**< search for e.c. aggregation of at least this size (has to be >= 3) */ 69 #define DEFAULT_MAXAGGRSIZE 4 /**< search for e.c. aggregation of at most this size (has to be >= minaggrsize) */ 70 #define DEFAULT_MAXBILINTERMS 500 /**< maximum number of bilinear terms allowed to be in a quadratic constraint */ 71 #define DEFAULT_MAXSTALLROUNDS 5 /**< maximum number of unsuccessful rounds in the e.c. aggregation search */ 72 73 #define SUBSCIP_NODELIMIT 100LL /**< node limit to solve the sub-SCIP */ 74 75 #define ADJUSTFACETTOL 1e-6 /**< adjust resulting facets in checkRikun() up to a violation of this value */ 76 #define USEDUALSIMPLEX TRUE /**< use dual or primal simplex algorithm? */ 77 78 /** first values for \f$2^n\f$ */ 79 static const int poweroftwo[] = { 1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096, 8192 }; 80 81 /* 82 * Data structures 83 */ 84 85 /** data to store a single edge-concave aggregations; an edge-concave aggregation of a quadratic constraint is a subset 86 * of nonconvex bilinear terms 87 */ 88 struct EcAggr 89 { 90 SCIP_VAR** vars; /**< variables */ 91 int nvars; /**< number of variables */ 92 int varsize; /**< size of vars array */ 93 94 SCIP_Real* termcoefs; /**< coefficients of bilinear terms */ 95 int* termvars1; /**< index of the first variable of each bilinear term */ 96 int* termvars2; /**< index of the second variable of each bilinear term*/ 97 int nterms; /**< number of bilinear terms in the aggregation */ 98 int termsize; /**< size of term{coefs,vars1,vars2} arrays */ 99 }; 100 typedef struct EcAggr SCIP_ECAGGR; 101 102 /** data to store all edge-concave aggregations and the remaining part of a nonlinear row of the form g(x) <= rhs */ 103 struct NlrowAggr 104 { 105 SCIP_NLROW* nlrow; /**< nonlinear row aggregation */ 106 SCIP_Bool rhsaggr; /**< consider nonlinear row aggregation for g(x) <= rhs (TRUE) or 107 * g(x) >= lhs (FALSE) */ 108 109 SCIP_ECAGGR** ecaggr; /**< array with all edge-concave aggregations */ 110 int necaggr; /**< number of edge-concave aggregation */ 111 112 SCIP_VAR** linvars; /**< linear variables */ 113 SCIP_Real* lincoefs; /**< linear coefficients */ 114 int nlinvars; /**< number of linear variables */ 115 int linvarssize; /**< size of linvars array */ 116 117 SCIP_VAR** quadvars; /**< quadratic variables */ 118 int* quadvar2aggr; /**< stores in which edge-concave aggregation the i-th quadratic variable 119 * is contained (< 0: in no edge-concave aggregation) */ 120 int nquadvars; /**< number of quadratic variables */ 121 int quadvarssize; /**< size of quadvars array */ 122 123 SCIP_VAR** remtermvars1; /**< first quadratic variable of remaining bilinear terms */ 124 SCIP_VAR** remtermvars2; /**< second quadratic variable of remaining bilinear terms */ 125 SCIP_Real* remtermcoefs; /**< coefficients for each remaining bilinear term */ 126 int nremterms; /**< number of remaining bilinear terms */ 127 int remtermsize; /**< size of remterm* arrays */ 128 129 SCIP_Real rhs; /**< rhs of the nonlinear row */ 130 SCIP_Real constant; /**< constant part of the nonlinear row */ 131 }; 132 typedef struct NlrowAggr SCIP_NLROWAGGR; 133 134 /** separator data */ 135 struct SCIP_SepaData 136 { 137 SCIP_NLROWAGGR** nlrowaggrs; /**< array containing all nonlinear row aggregations */ 138 int nnlrowaggrs; /**< number of nonlinear row aggregations */ 139 int nlrowaggrssize; /**< size of nlrowaggrs array */ 140 SCIP_Bool searchedforaggr; /**< flag if we already searched for nlrow aggregation candidates */ 141 int minaggrsize; /**< only search for e.c. aggregations of at least this size (has to be >= 3) */ 142 int maxaggrsize; /**< only search for e.c. aggregations of at most this size (has to be >= minaggrsize) */ 143 int maxecsize; /**< largest edge concave aggregation size */ 144 int maxbilinterms; /**< maximum number of bilinear terms allowed to be in a quadratic constraint */ 145 int maxstallrounds; /**< maximum number of unsuccessful rounds in the e.c. aggregation search */ 146 147 SCIP_LPI* lpi; /**< LP interface to solve the LPs to compute the facets of the convex envelopes */ 148 int lpisize; /**< maximum size of e.c. aggregations which can be handled by the LP interface */ 149 150 SCIP_Real cutmaxrange; /**< maximal coef range of a cut (maximal coefficient divided by minimal 151 * coefficient) in order to be added to LP relaxation */ 152 SCIP_Bool dynamiccuts; /**< should generated cuts be removed from the LP if they are no longer tight? */ 153 SCIP_Real minviolation; /**< minimal violation of an e.c. cut to be separated */ 154 155 int maxrounds; /**< maximal number of separation rounds per node (-1: unlimited) */ 156 int maxroundsroot; /**< maximal number of separation rounds in the root node (-1: unlimited) */ 157 int maxdepth; /**< maximal depth at which the separator is applied */ 158 int maxsepacuts; /**< maximal number of e.c. cuts separated per separation round */ 159 int maxsepacutsroot; /**< maximal number of e.c. cuts separated per separation round in root node */ 160 161 #ifdef SCIP_STATISTIC 162 SCIP_Real aggrsearchtime; /**< total time spent for searching edge concave aggregations */ 163 int nlhsnlrowaggrs; /**< number of found nonlinear row aggregations for SCIP_NLROWs of the form g(x) <= rhs */ 164 int nrhsnlrowaggrs; /**< number of found nonlinear row aggregations for SCIP_NLROWs of the form g(x) >= lhs */ 165 #endif 166 }; 167 168 169 /* 170 * Local methods 171 */ 172 173 /** creates an empty edge-concave aggregation (without bilinear terms) */ 174 static 175 SCIP_RETCODE ecaggrCreateEmpty( 176 SCIP* scip, /**< SCIP data structure */ 177 SCIP_ECAGGR** ecaggr, /**< pointer to store the edge-concave aggregation */ 178 int nquadvars, /**< number of quadratic variables */ 179 int nquadterms /**< number of bilinear terms */ 180 ) 181 { 182 assert(scip != NULL); 183 assert(ecaggr != NULL); 184 assert(nquadvars > 0); 185 assert(nquadterms >= nquadvars); 186 187 SCIP_CALL( SCIPallocBlockMemory(scip, ecaggr) ); 188 189 (*ecaggr)->nvars = 0; 190 (*ecaggr)->nterms = 0; 191 (*ecaggr)->varsize = nquadvars; 192 (*ecaggr)->termsize = nquadterms; 193 194 /* allocate enough memory for the quadratic variables and bilinear terms */ 195 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*ecaggr)->vars, nquadvars) ); 196 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*ecaggr)->termcoefs, nquadterms) ); 197 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*ecaggr)->termvars1, nquadterms) ); 198 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*ecaggr)->termvars2, nquadterms) ); 199 200 return SCIP_OKAY; 201 } 202 203 /** frees an edge-concave aggregation */ 204 static 205 SCIP_RETCODE ecaggrFree( 206 SCIP* scip, /**< SCIP data structure */ 207 SCIP_ECAGGR** ecaggr /**< pointer to store the edge-concave aggregation */ 208 ) 209 { 210 assert(scip != NULL); 211 assert(ecaggr != NULL); 212 213 SCIPfreeBlockMemoryArray(scip, &((*ecaggr)->termcoefs), (*ecaggr)->termsize); 214 SCIPfreeBlockMemoryArray(scip, &((*ecaggr)->termvars1), (*ecaggr)->termsize); 215 SCIPfreeBlockMemoryArray(scip, &((*ecaggr)->termvars2), (*ecaggr)->termsize); 216 SCIPfreeBlockMemoryArray(scip, &((*ecaggr)->vars), (*ecaggr)->varsize); 217 218 SCIPfreeBlockMemory(scip, ecaggr); 219 *ecaggr = NULL; 220 221 return SCIP_OKAY; 222 } 223 224 /** adds a quadratic variable to an edge-concave aggregation */ 225 static 226 SCIP_RETCODE ecaggrAddQuadvar( 227 SCIP_ECAGGR* ecaggr, /**< pointer to store the edge-concave aggregation */ 228 SCIP_VAR* x /**< first variable */ 229 ) 230 { 231 ecaggr->vars[ ecaggr->nvars++ ] = x; 232 return SCIP_OKAY; 233 } 234 235 /** adds a bilinear term to an edge-concave aggregation */ 236 static 237 SCIP_RETCODE ecaggrAddBilinTerm( 238 SCIP* scip, /**< SCIP data structure */ 239 SCIP_ECAGGR* ecaggr, /**< pointer to store the edge-concave aggregation */ 240 SCIP_VAR* x, /**< first variable */ 241 SCIP_VAR* y, /**< second variable */ 242 SCIP_Real coef /**< bilinear coefficient */ 243 ) 244 { 245 int idx1; 246 int idx2; 247 int i; 248 249 assert(x != NULL); 250 assert(y != NULL); 251 assert(ecaggr->nterms + 1 <= ((ecaggr->nvars + 1) * ecaggr->nvars) / 2); 252 assert(!SCIPisZero(scip, coef)); 253 254 idx1 = -1; 255 idx2 = -1; 256 257 /* search for the quadratic variables in the e.c. aggregation */ 258 for( i = 0; i < ecaggr->nvars && (idx1 == -1 || idx2 == -1); ++i ) 259 { 260 if( ecaggr->vars[i] == x ) 261 idx1 = i; 262 if( ecaggr->vars[i] == y ) 263 idx2 = i; 264 } 265 266 assert(idx1 != -1 && idx2 != -1); 267 268 ecaggr->termcoefs[ ecaggr->nterms ] = coef; 269 ecaggr->termvars1[ ecaggr->nterms ] = idx1; 270 ecaggr->termvars2[ ecaggr->nterms ] = idx2; 271 ++(ecaggr->nterms); 272 273 return SCIP_OKAY; 274 } 275 276 #ifdef SCIP_DEBUG 277 /** prints an edge-concave aggregation */ 278 static 279 void ecaggrPrint( 280 SCIP* scip, /**< SCIP data structure */ 281 SCIP_ECAGGR* ecaggr /**< pointer to store the edge-concave aggregation */ 282 ) 283 { 284 int i; 285 286 assert(scip != NULL); 287 assert(ecaggr != NULL); 288 289 SCIPdebugMsg(scip, " nvars = %d nterms = %d\n", ecaggr->nvars, ecaggr->nterms); 290 SCIPdebugMsg(scip, " vars: "); 291 for( i = 0; i < ecaggr->nvars; ++i ) 292 SCIPdebugMsgPrint(scip, "%s ", SCIPvarGetName(ecaggr->vars[i])); 293 SCIPdebugMsgPrint(scip, "\n"); 294 295 SCIPdebugMsg(scip, " terms: "); 296 for( i = 0; i < ecaggr->nterms; ++i ) 297 { 298 SCIP_VAR* x; 299 SCIP_VAR* y; 300 301 x = ecaggr->vars[ ecaggr->termvars1[i] ]; 302 y = ecaggr->vars[ ecaggr->termvars2[i] ]; 303 SCIPdebugMsgPrint(scip, "%e %s * %s ", ecaggr->termcoefs[i], SCIPvarGetName(x), SCIPvarGetName(y) ); 304 } 305 SCIPdebugMsgPrint(scip, "\n"); 306 } 307 #endif 308 309 /** stores linear terms in a given nonlinear row aggregation */ 310 static 311 SCIP_RETCODE nlrowaggrStoreLinearTerms( 312 SCIP* scip, /**< SCIP data structure */ 313 SCIP_NLROWAGGR* nlrowaggr, /**< nonlinear row aggregation */ 314 SCIP_VAR** linvars, /**< linear variables */ 315 SCIP_Real* lincoefs, /**< linear coefficients */ 316 int nlinvars /**< number of linear variables */ 317 ) 318 { 319 assert(scip != NULL); 320 assert(nlrowaggr != NULL); 321 assert(linvars != NULL || nlinvars == 0); 322 assert(lincoefs != NULL || nlinvars == 0); 323 assert(nlinvars >= 0); 324 325 nlrowaggr->nlinvars = 0; 326 nlrowaggr->linvarssize = 0; 327 nlrowaggr->linvars = NULL; 328 nlrowaggr->lincoefs = NULL; 329 330 if( nlinvars == 0 ) 331 return SCIP_OKAY; 332 333 SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &nlrowaggr->linvars, linvars, nlinvars) ); 334 SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &nlrowaggr->lincoefs, lincoefs, nlinvars) ); 335 nlrowaggr->nlinvars = nlinvars; 336 nlrowaggr->linvarssize = nlinvars; 337 338 /* if we have a nlrow of the form g(x) >= lhs, multiply every coefficient by -1 */ 339 if( !nlrowaggr->rhsaggr ) 340 { 341 int i; 342 343 for( i = 0; i < nlrowaggr->nlinvars; ++i ) 344 nlrowaggr->lincoefs[i] *= -1.0; 345 } 346 347 return SCIP_OKAY; 348 } 349 350 /** adds linear term to a given nonlinear row aggregation */ 351 static 352 SCIP_RETCODE nlrowaggrAddLinearTerm( 353 SCIP* scip, /**< SCIP data structure */ 354 SCIP_NLROWAGGR* nlrowaggr, /**< nonlinear row aggregation */ 355 SCIP_VAR* linvar, /**< linear variable */ 356 SCIP_Real lincoef /**< coefficient */ 357 ) 358 { 359 assert(scip != NULL); 360 assert(nlrowaggr != NULL); 361 assert(linvar != NULL); 362 363 if( nlrowaggr->nlinvars == nlrowaggr->linvarssize ) 364 { 365 int newsize = SCIPcalcMemGrowSize(scip, nlrowaggr->linvarssize+1); 366 SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &nlrowaggr->linvars, nlrowaggr->linvarssize, newsize) ); 367 SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &nlrowaggr->lincoefs, nlrowaggr->linvarssize, newsize) ); 368 nlrowaggr->linvarssize = newsize; 369 } 370 assert(nlrowaggr->linvarssize > nlrowaggr->nlinvars); 371 372 /* if we have a nlrow of the form g(x) >= lhs, multiply coefficient by -1 */ 373 if( !nlrowaggr->rhsaggr ) 374 lincoef = -lincoef; 375 376 nlrowaggr->linvars[nlrowaggr->nlinvars] = linvar; 377 nlrowaggr->lincoefs[nlrowaggr->nlinvars] = lincoef; 378 ++nlrowaggr->nlinvars; 379 380 return SCIP_OKAY; 381 } 382 383 /** adds quadratic variable to a given nonlinear row aggregation */ 384 static 385 SCIP_RETCODE nlrowaggrAddQuadraticVar( 386 SCIP* scip, /**< SCIP data structure */ 387 SCIP_NLROWAGGR* nlrowaggr, /**< nonlinear row aggregation */ 388 SCIP_VAR* quadvar /**< quadratic variable */ 389 ) 390 { 391 assert(scip != NULL); 392 assert(nlrowaggr != NULL); 393 assert(quadvar != NULL); 394 395 SCIP_CALL( SCIPensureBlockMemoryArray(scip, &nlrowaggr->quadvars, &nlrowaggr->quadvarssize, nlrowaggr->nquadvars+1) ); 396 assert(nlrowaggr->quadvarssize > nlrowaggr->nquadvars); 397 nlrowaggr->quadvars[nlrowaggr->nquadvars] = quadvar; 398 ++nlrowaggr->nquadvars; 399 400 return SCIP_OKAY; 401 } 402 403 /** adds a remaining bilinear term to a given nonlinear row aggregation */ 404 static 405 SCIP_RETCODE nlrowaggrAddRemBilinTerm( 406 SCIP_NLROWAGGR* nlrowaggr, /**< nonlinear row aggregation */ 407 SCIP_VAR* x, /**< first variable */ 408 SCIP_VAR* y, /**< second variable */ 409 SCIP_Real coef /**< bilinear coefficient */ 410 ) 411 { 412 assert(nlrowaggr != NULL); 413 assert(x != NULL); 414 assert(y != NULL); 415 assert(coef != 0.0); 416 assert(nlrowaggr->remtermcoefs != NULL); 417 assert(nlrowaggr->remtermvars1 != NULL); 418 assert(nlrowaggr->remtermvars2 != NULL); 419 420 nlrowaggr->remtermcoefs[ nlrowaggr->nremterms ] = coef; 421 nlrowaggr->remtermvars1[ nlrowaggr->nremterms ] = x; 422 nlrowaggr->remtermvars2[ nlrowaggr->nremterms ] = y; 423 ++(nlrowaggr->nremterms); 424 425 return SCIP_OKAY; 426 } 427 428 /** creates a nonlinear row aggregation */ 429 static 430 SCIP_RETCODE nlrowaggrCreate( 431 SCIP* scip, /**< SCIP data structure */ 432 SCIP_NLROW* nlrow, /**< nonlinear row */ 433 SCIP_NLROWAGGR** nlrowaggr, /**< pointer to store the nonlinear row aggregation */ 434 int* quadvar2aggr, /**< mapping between quadratic variables and edge-concave aggregation 435 * stores a negative value if the quadratic variables does not belong 436 * to any aggregation */ 437 int nfound, /**< number of edge-concave aggregations */ 438 SCIP_Bool rhsaggr /**< consider nonlinear row aggregation for g(x) <= rhs (TRUE) or 439 * lhs <= g(x) (FALSE) */ 440 ) 441 { 442 SCIP_EXPR* expr; 443 int* aggrnvars; /* count the number of variables in each e.c. aggregations */ 444 int* aggrnterms; /* count the number of bilinear terms in each e.c. aggregations */ 445 int nquadvars; 446 int nremterms; 447 int i; 448 449 assert(scip != NULL); 450 assert(nlrow != NULL); 451 assert(nlrowaggr != NULL); 452 assert(quadvar2aggr != NULL); 453 assert(nfound > 0); 454 455 expr = SCIPnlrowGetExpr(nlrow); 456 SCIPexprGetQuadraticData(expr, NULL, NULL, NULL, NULL, &nquadvars, NULL, NULL, NULL); 457 nremterms = 0; 458 459 SCIP_CALL( SCIPallocClearBufferArray(scip, &aggrnvars, nfound) ); 460 SCIP_CALL( SCIPallocClearBufferArray(scip, &aggrnterms, nfound) ); 461 462 /* create an empty nonlinear row aggregation */ 463 SCIP_CALL( SCIPallocBlockMemory(scip, nlrowaggr) ); 464 (*nlrowaggr)->nlrow = nlrow; 465 (*nlrowaggr)->rhsaggr = rhsaggr; 466 (*nlrowaggr)->rhs = rhsaggr ? SCIPnlrowGetRhs(nlrow) : -SCIPnlrowGetLhs(nlrow); 467 (*nlrowaggr)->constant = rhsaggr ? SCIPnlrowGetConstant(nlrow) : -SCIPnlrowGetConstant(nlrow); 468 469 (*nlrowaggr)->quadvars = NULL; 470 (*nlrowaggr)->nquadvars = 0; 471 (*nlrowaggr)->quadvarssize = 0; 472 (*nlrowaggr)->quadvar2aggr = NULL; 473 (*nlrowaggr)->remtermcoefs = NULL; 474 (*nlrowaggr)->remtermvars1 = NULL; 475 (*nlrowaggr)->remtermvars2 = NULL; 476 (*nlrowaggr)->nremterms = 0; 477 478 /* copy quadvar2aggr array */ 479 SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*nlrowaggr)->quadvar2aggr, quadvar2aggr, nquadvars) ); 480 481 /* store all linear terms */ 482 SCIP_CALL( nlrowaggrStoreLinearTerms(scip, *nlrowaggr, SCIPnlrowGetLinearVars(nlrow), SCIPnlrowGetLinearCoefs(nlrow), 483 SCIPnlrowGetNLinearVars(nlrow)) ); 484 485 /* store all quadratic variables and additional linear terms */ 486 /* count the number of variables in each e.c. aggregation */ 487 /* count the number of square and bilinear terms in each e.c. aggregation */ 488 for( i = 0; i < nquadvars; ++i ) 489 { 490 SCIP_EXPR* qterm; 491 SCIP_Real lincoef; 492 SCIP_Real sqrcoef; 493 int idx1; 494 int nadjbilin; 495 int* adjbilin; 496 int j; 497 498 SCIPexprGetQuadraticQuadTerm(expr, i, &qterm, &lincoef, &sqrcoef, &nadjbilin, &adjbilin, NULL); 499 assert(SCIPisExprVar(scip, qterm)); 500 501 SCIP_CALL( nlrowaggrAddQuadraticVar(scip, *nlrowaggr, SCIPgetVarExprVar(qterm)) ); 502 503 if( lincoef != 0.0 ) 504 { 505 SCIP_CALL( nlrowaggrAddLinearTerm(scip, *nlrowaggr, SCIPgetVarExprVar(qterm), lincoef) ); 506 } 507 508 if( quadvar2aggr[i] >= 0) 509 ++aggrnvars[ quadvar2aggr[i] ]; 510 511 idx1 = quadvar2aggr[i]; 512 if( rhsaggr ) 513 sqrcoef = -sqrcoef; 514 515 /* variable has to belong to an e.c. aggregation; square term has to be concave */ 516 if( idx1 >= 0 && SCIPisNegative(scip, sqrcoef) ) 517 ++aggrnterms[idx1]; 518 else 519 ++nremterms; 520 521 for( j = 0; j < nadjbilin; ++j ) 522 { 523 SCIP_EXPR* qterm1; 524 int pos2; 525 int idx2; 526 527 SCIPexprGetQuadraticBilinTerm(expr, adjbilin[j], &qterm1, NULL, NULL, &pos2, NULL); 528 529 /* only handle qterm1 == qterm here; the other case will be handled when its turn for qterm2 to be qterm */ 530 if( qterm1 != qterm ) 531 continue; 532 533 idx2 = quadvar2aggr[pos2]; 534 535 /* variables have to belong to the same e.c. aggregation; bilinear term has to be concave */ 536 if( idx1 >= 0 && idx2 >= 0 && idx1 == idx2 ) 537 ++aggrnterms[idx1]; 538 else 539 ++nremterms; 540 } 541 } 542 assert((*nlrowaggr)->nquadvars == nquadvars); 543 544 /* create all edge-concave aggregations (empty) and remaining terms */ 545 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*nlrowaggr)->ecaggr, nfound) ); 546 if( nremterms > 0 ) 547 { 548 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*nlrowaggr)->remtermcoefs, nremterms) ); 549 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*nlrowaggr)->remtermvars1, nremterms) ); 550 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*nlrowaggr)->remtermvars2, nremterms) ); 551 (*nlrowaggr)->remtermsize = nremterms; 552 } 553 (*nlrowaggr)->necaggr = nfound; 554 555 for( i = 0; i < nfound; ++i ) 556 { 557 SCIP_CALL( ecaggrCreateEmpty(scip, &(*nlrowaggr)->ecaggr[i], aggrnvars[i], aggrnterms[i]) ); 558 } 559 560 /* add quadratic variables to the edge-concave aggregations */ 561 for( i = 0; i < nquadvars; ++i ) 562 { 563 int idx; 564 565 idx = quadvar2aggr[i]; 566 567 if( idx >= 0) 568 { 569 SCIP_EXPR* qterm; 570 571 SCIPdebugMsg(scip, "add quadvar %d to aggr. %d\n", i, idx); 572 573 SCIPexprGetQuadraticQuadTerm(expr, i, &qterm, NULL, NULL, NULL, NULL, NULL); 574 assert(SCIPisExprVar(scip, qterm)); 575 576 SCIP_CALL( ecaggrAddQuadvar((*nlrowaggr)->ecaggr[idx], SCIPgetVarExprVar(qterm)) ); 577 } 578 } 579 580 /* add the bilinear/square terms to the edge-concave aggregations or in the remaining part */ 581 for( i = 0; i < nquadvars; ++i ) 582 { 583 SCIP_EXPR* qterm; 584 SCIP_VAR* x; 585 SCIP_Real coef; 586 int idx1; 587 int nadjbilin; 588 int* adjbilin; 589 int j; 590 591 SCIPexprGetQuadraticQuadTerm(expr, i, &qterm, NULL, &coef, &nadjbilin, &adjbilin, NULL); 592 593 x = SCIPgetVarExprVar(qterm); 594 595 idx1 = quadvar2aggr[i]; 596 if( rhsaggr ) 597 coef = -coef; 598 599 if( idx1 >= 0 && SCIPisNegative(scip, coef) ) 600 { 601 SCIP_CALL( ecaggrAddBilinTerm(scip, (*nlrowaggr)->ecaggr[idx1], x, x, coef) ); 602 SCIPdebugMsg(scip, "add term %e *%d^2 to aggr. %d\n", coef, i, idx1); 603 } 604 else 605 { 606 SCIP_CALL( nlrowaggrAddRemBilinTerm(*nlrowaggr, x, x, coef) ); 607 SCIPdebugMsg(scip, "add term %e *%d^2 to the remaining part\n", coef, idx1); 608 } 609 610 for( j = 0; j < nadjbilin; ++j ) 611 { 612 SCIP_EXPR* qterm1; 613 SCIP_EXPR* qterm2; 614 int pos2; 615 int idx2; 616 SCIP_VAR* y; 617 618 SCIPexprGetQuadraticBilinTerm(expr, adjbilin[j], &qterm1, &qterm2, &coef, &pos2, NULL); 619 620 /* only handle qterm1 == qterm here; the other case will be handled when its turn for qterm2 to be qterm */ 621 if( qterm1 != qterm ) 622 continue; 623 624 y = SCIPgetVarExprVar(qterm2); 625 626 idx2 = quadvar2aggr[pos2]; 627 if( rhsaggr ) 628 coef = -coef; 629 630 if( idx1 >= 0 && idx2 >= 0 && idx1 == idx2 ) 631 { 632 SCIP_CALL( ecaggrAddBilinTerm(scip, (*nlrowaggr)->ecaggr[idx1], x, y, coef) ); 633 SCIPdebugMsg(scip, "add term %e *%d*%d to aggr. %d\n", coef, i, pos2, idx1); 634 } 635 else 636 { 637 SCIP_CALL( nlrowaggrAddRemBilinTerm(*nlrowaggr, x, y, coef) ); 638 SCIPdebugMsg(scip, "add term %e *%d*%d to the remaining part\n", coef, i, pos2); 639 } 640 } 641 } 642 643 /* free allocated memory */ 644 SCIPfreeBufferArray(scip, &aggrnterms); 645 SCIPfreeBufferArray(scip, &aggrnvars); 646 647 return SCIP_OKAY; 648 } 649 650 /** frees a nonlinear row aggregation */ 651 static 652 SCIP_RETCODE nlrowaggrFree( 653 SCIP* scip, /**< SCIP data structure */ 654 SCIP_NLROWAGGR** nlrowaggr /**< pointer to free the nonlinear row aggregation */ 655 ) 656 { 657 int i; 658 659 assert(scip != NULL); 660 assert(nlrowaggr != NULL); 661 assert(*nlrowaggr != NULL); 662 (*nlrowaggr)->nlrow = NULL; 663 assert((*nlrowaggr)->quadvars != NULL); 664 assert((*nlrowaggr)->nquadvars > 0); 665 assert((*nlrowaggr)->nremterms >= 0); 666 667 /* free remaining part */ 668 SCIPfreeBlockMemoryArrayNull(scip, &(*nlrowaggr)->remtermcoefs, (*nlrowaggr)->remtermsize); 669 SCIPfreeBlockMemoryArrayNull(scip, &(*nlrowaggr)->remtermvars1, (*nlrowaggr)->remtermsize); 670 SCIPfreeBlockMemoryArrayNull(scip, &(*nlrowaggr)->remtermvars2, (*nlrowaggr)->remtermsize); 671 672 /* free quadratic variables */ 673 SCIPfreeBlockMemoryArray(scip, &(*nlrowaggr)->quadvars, (*nlrowaggr)->quadvarssize); 674 SCIPfreeBlockMemoryArray(scip, &(*nlrowaggr)->quadvar2aggr, (*nlrowaggr)->nquadvars); 675 676 /* free linear part */ 677 if( (*nlrowaggr)->nlinvars > 0 ) 678 { 679 SCIPfreeBlockMemoryArray(scip, &(*nlrowaggr)->linvars, (*nlrowaggr)->linvarssize); 680 SCIPfreeBlockMemoryArray(scip, &(*nlrowaggr)->lincoefs, (*nlrowaggr)->linvarssize); 681 } 682 683 /* free edge-concave aggregations */ 684 for( i = 0; i < (*nlrowaggr)->necaggr; ++i ) 685 { 686 SCIP_CALL( ecaggrFree(scip, &(*nlrowaggr)->ecaggr[i]) ); 687 } 688 SCIPfreeBlockMemoryArray(scip, &(*nlrowaggr)->ecaggr, (*nlrowaggr)->necaggr); 689 690 /* free nlrow aggregation */ 691 SCIPfreeBlockMemory(scip, nlrowaggr); 692 693 return SCIP_OKAY; 694 } 695 696 #ifdef SCIP_DEBUG 697 /** prints a nonlinear row aggregation */ 698 static 699 void nlrowaggrPrint( 700 SCIP* scip, /**< SCIP data structure */ 701 SCIP_NLROWAGGR* nlrowaggr /**< nonlinear row aggregation */ 702 ) 703 { 704 int i; 705 706 SCIPdebugMsg(scip, " nlrowaggr rhs = %e\n", nlrowaggr->rhs); 707 SCIPdebugMsg(scip, " #remaining terms = %d\n", nlrowaggr->nremterms); 708 709 SCIPdebugMsg(scip, "remaining terms: "); 710 for( i = 0; i < nlrowaggr->nremterms; ++i ) 711 SCIPdebugMsgPrint(scip, "%e %s * %s + ", nlrowaggr->remtermcoefs[i], SCIPvarGetName(nlrowaggr->remtermvars1[i]), 712 SCIPvarGetName(nlrowaggr->remtermvars2[i]) ); 713 for( i = 0; i < nlrowaggr->nlinvars; ++i ) 714 SCIPdebugMsgPrint(scip, "%e %s + ", nlrowaggr->lincoefs[i], SCIPvarGetName(nlrowaggr->linvars[i]) ); 715 SCIPdebugMsgPrint(scip, "\n"); 716 717 for( i = 0; i < nlrowaggr->necaggr; ++i ) 718 { 719 SCIPdebugMsg(scip, "print e.c. aggr %d\n", i); 720 ecaggrPrint(scip, nlrowaggr->ecaggr[i]); 721 } 722 return; 723 } 724 #endif 725 726 /** creates separator data */ 727 static 728 SCIP_RETCODE sepadataCreate( 729 SCIP* scip, /**< SCIP data structure */ 730 SCIP_SEPADATA** sepadata /**< pointer to store separator data */ 731 ) 732 { 733 assert(scip != NULL); 734 assert(sepadata != NULL); 735 736 SCIP_CALL( SCIPallocBlockMemory(scip, sepadata) ); 737 BMSclearMemory(*sepadata); 738 739 return SCIP_OKAY; 740 } 741 742 /** frees all nonlinear row aggregations */ 743 static 744 SCIP_RETCODE sepadataFreeNlrows( 745 SCIP* scip, /**< SCIP data structure */ 746 SCIP_SEPADATA* sepadata /**< pointer to store separator data */ 747 ) 748 { 749 assert(scip != NULL); 750 assert(sepadata != NULL); 751 752 /* free nonlinear row aggregations */ 753 if( sepadata->nlrowaggrs != NULL ) 754 { 755 int i; 756 757 for( i = sepadata->nnlrowaggrs - 1; i >= 0; --i ) 758 { 759 SCIP_CALL( nlrowaggrFree(scip, &sepadata->nlrowaggrs[i]) ); 760 } 761 762 SCIPfreeBlockMemoryArray(scip, &sepadata->nlrowaggrs, sepadata->nlrowaggrssize); 763 764 sepadata->nlrowaggrs = NULL; 765 sepadata->nnlrowaggrs = 0; 766 sepadata->nlrowaggrssize = 0; 767 } 768 769 return SCIP_OKAY; 770 } 771 772 /** frees separator data */ 773 static 774 SCIP_RETCODE sepadataFree( 775 SCIP* scip, /**< SCIP data structure */ 776 SCIP_SEPADATA** sepadata /**< pointer to store separator data */ 777 ) 778 { 779 assert(scip != NULL); 780 assert(sepadata != NULL); 781 assert(*sepadata != NULL); 782 783 /* free nonlinear row aggregations */ 784 SCIP_CALL( sepadataFreeNlrows(scip, *sepadata) ); 785 786 /* free LP interface */ 787 if( (*sepadata)->lpi != NULL ) 788 { 789 SCIP_CALL( SCIPlpiFree(&((*sepadata)->lpi)) ); 790 (*sepadata)->lpisize = 0; 791 } 792 793 SCIPfreeBlockMemory(scip, sepadata); 794 795 return SCIP_OKAY; 796 } 797 798 /** adds a nonlinear row aggregation to the separator data */ 799 static 800 SCIP_RETCODE sepadataAddNlrowaggr( 801 SCIP* scip, /**< SCIP data structure */ 802 SCIP_SEPADATA* sepadata, /**< separator data */ 803 SCIP_NLROWAGGR* nlrowaggr /**< non-linear row aggregation */ 804 ) 805 { 806 int i; 807 808 assert(scip != NULL); 809 assert(sepadata != NULL); 810 assert(nlrowaggr != NULL); 811 812 if( sepadata->nlrowaggrssize == 0 ) 813 { 814 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &sepadata->nlrowaggrs, 2) ); /*lint !e506*/ 815 sepadata->nlrowaggrssize = 2; 816 } 817 else if( sepadata->nlrowaggrssize < sepadata->nnlrowaggrs + 1 ) 818 { 819 SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &sepadata->nlrowaggrs, sepadata->nlrowaggrssize, 2 * sepadata->nlrowaggrssize) ); /*lint !e506 !e647*/ 820 sepadata->nlrowaggrssize *= 2; 821 assert(sepadata->nlrowaggrssize >= sepadata->nnlrowaggrs + 1); 822 } 823 824 sepadata->nlrowaggrs[ sepadata->nnlrowaggrs ] = nlrowaggr; 825 ++(sepadata->nnlrowaggrs); 826 827 /* update maximum e.c. aggregation size */ 828 for( i = 0; i < nlrowaggr->necaggr; ++i ) 829 sepadata->maxecsize = MAX(sepadata->maxecsize, nlrowaggr->ecaggr[i]->nvars); 830 831 #ifdef SCIP_STATISTIC 832 /* update statistics */ 833 if( nlrowaggr->rhsaggr ) 834 ++(sepadata->nrhsnlrowaggrs); 835 else 836 ++(sepadata->nlhsnlrowaggrs); 837 #endif 838 839 return SCIP_OKAY; 840 } 841 842 /** returns min{val-lb,ub-val} / (ub-lb) */ 843 static 844 SCIP_Real phi( 845 SCIP* scip, /**< SCIP data structure */ 846 SCIP_Real val, /**< solution value */ 847 SCIP_Real lb, /**< lower bound */ 848 SCIP_Real ub /**< upper bound */ 849 ) 850 { 851 if( SCIPisFeasEQ(scip, lb, ub) ) 852 return 0.0; 853 854 /* adjust */ 855 val = MAX(val, lb); 856 val = MIN(val, ub); 857 858 return MIN(ub - val, val - lb) / (ub - lb); 859 } 860 861 /** creates an MIP to search for cycles with an odd number of positive edges in the graph representation of a nonlinear row 862 * 863 * The model uses directed binary arc flow variables. 864 * We introduce for all quadratic elements a forward and backward edge. 865 * If the term is quadratic (e.g., loop in the graph) we fix the corresponding variables to zero. 866 * This leads to an easy mapping between quadratic elements and the variables of the MIP. 867 */ 868 static 869 SCIP_RETCODE createMIP( 870 SCIP* scip, /**< SCIP data structure */ 871 SCIP* subscip, /**< auxiliary SCIP to search aggregations */ 872 SCIP_SEPADATA* sepadata, /**< separator data */ 873 SCIP_NLROW* nlrow, /**< nonlinear row */ 874 SCIP_Bool rhsaggr, /**< consider nonlinear row aggregation for g(x) <= rhs (TRUE) or 875 * lhs <= g(x) (FALSE) */ 876 SCIP_VAR** forwardarcs, /**< array to store all forward arc variables */ 877 SCIP_VAR** backwardarcs, /**< array to store all backward arc variables */ 878 SCIP_Real* nodeweights, /**< weights for each node of the graph */ 879 int* nedges, /**< pointer to store the number of nonexcluded edges in the graph */ 880 int* narcs /**< pointer to store the number of created arc variables (number of square and bilinear terms) */ 881 ) 882 { 883 SCIP_VAR** oddcyclearcs; 884 SCIP_CONS** flowcons; 885 SCIP_CONS* cyclelengthcons; 886 SCIP_CONS* oddcyclecons; 887 char name[SCIP_MAXSTRLEN]; 888 SCIP_EXPR* expr; 889 int noddcyclearcs; 890 int nnodes; 891 int nquadexprs; 892 int nbilinexprs; 893 int i; 894 int arcidx; 895 896 assert(subscip != NULL); 897 assert(forwardarcs != NULL); 898 assert(backwardarcs != NULL); 899 assert(nedges != NULL); 900 assert(sepadata->minaggrsize <= sepadata->maxaggrsize); 901 902 expr = SCIPnlrowGetExpr(nlrow); 903 SCIPexprGetQuadraticData(expr, NULL, NULL, NULL, NULL, &nquadexprs, &nbilinexprs, NULL, NULL); 904 905 nnodes = nquadexprs; 906 *nedges = 0; 907 *narcs = 0; 908 909 assert(nnodes > 0); 910 911 noddcyclearcs = 0; 912 SCIP_CALL( SCIPallocBufferArray(subscip, &oddcyclearcs, 2*nbilinexprs) ); 913 914 /* create problem with default plug-ins */ 915 SCIP_CALL( SCIPcreateProbBasic(subscip, "E.C. aggregation MIP") ); 916 SCIP_CALL( SCIPsetObjsense(subscip, SCIP_OBJSENSE_MAXIMIZE) ); 917 SCIP_CALL( SCIPincludeDefaultPlugins(subscip) ); 918 919 /* create forward and backward arc variables */ 920 for( i = 0; i < nquadexprs; ++i ) 921 { 922 SCIP_EXPR* qterm; 923 SCIP_Real coef; 924 int nadjbilin; 925 int* adjbilin; 926 int j; 927 928 SCIPexprGetQuadraticQuadTerm(expr, i, &qterm, NULL, &coef, &nadjbilin, &adjbilin, NULL); 929 930 if( !SCIPisZero(scip, coef) ) 931 { 932 /* squares (loops) are fixed to zero */ 933 SCIPdebugMsg(scip, "edge {%d,%d} = {%s,%s} coeff=%e edgeweight=0\n", i, i, 934 SCIPvarGetName(SCIPgetVarExprVar(qterm)), SCIPvarGetName(SCIPgetVarExprVar(qterm)), 935 coef); 936 937 (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "x#%d#%d", i, i); 938 SCIP_CALL( SCIPcreateVarBasic(subscip, &forwardarcs[*narcs], name, 0.0, 0.0, 0.01, SCIP_VARTYPE_BINARY) ); 939 SCIP_CALL( SCIPaddVar(subscip, forwardarcs[*narcs]) ); 940 941 SCIP_CALL( SCIPcreateVarBasic(subscip, &backwardarcs[*narcs], name, 0.0, 0.0, 0.01, SCIP_VARTYPE_BINARY) ); 942 SCIP_CALL( SCIPaddVar(subscip, backwardarcs[*narcs]) ); 943 944 ++*narcs; 945 } 946 947 for( j = 0 ; j < nadjbilin; ++j ) 948 { 949 SCIP_EXPR* qterm1; 950 SCIP_EXPR* qterm2; 951 int pos2; 952 SCIP_Real edgeweight; 953 SCIP_CONS* noparallelcons; 954 955 SCIPexprGetQuadraticBilinTerm(expr, adjbilin[j], &qterm1, &qterm2, &coef, &pos2, NULL); 956 957 /* handle qterm == qterm2 later */ 958 if( qterm1 != qterm ) 959 continue; 960 961 edgeweight = nodeweights[i] + nodeweights[pos2]; 962 SCIPdebugMsg(scip, "edge {%d,%d} = {%s,%s} coeff=%e edgeweight=%e\n", i, pos2, 963 SCIPvarGetName(SCIPgetVarExprVar(qterm1)), SCIPvarGetName(SCIPgetVarExprVar(qterm2)), 964 coef, edgeweight); 965 966 (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "x#%d#%d", i, pos2); 967 SCIP_CALL( SCIPcreateVarBasic(subscip, &forwardarcs[*narcs], name, 0.0, 1.0, 0.01 + edgeweight, SCIP_VARTYPE_BINARY) ); 968 SCIP_CALL( SCIPaddVar(subscip, forwardarcs[*narcs]) ); 969 970 (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "x#%d#%d", i, pos2); 971 SCIP_CALL( SCIPcreateVarBasic(subscip, &backwardarcs[*narcs], name, 0.0, 1.0, 0.01 + edgeweight, SCIP_VARTYPE_BINARY) ); 972 SCIP_CALL( SCIPaddVar(subscip, backwardarcs[*narcs]) ); 973 974 ++(*nedges); 975 976 /* store all arcs which are important for the odd cycle property (no loops) */ 977 if( rhsaggr && SCIPisPositive(scip, coef) ) 978 { 979 assert(noddcyclearcs < 2*nbilinexprs-1); 980 oddcyclearcs[noddcyclearcs++] = forwardarcs[i]; 981 oddcyclearcs[noddcyclearcs++] = backwardarcs[i]; 982 } 983 984 if( !rhsaggr && SCIPisNegative(scip, coef) ) 985 { 986 assert(noddcyclearcs < 2*nbilinexprs-1); 987 oddcyclearcs[noddcyclearcs++] = forwardarcs[i]; 988 oddcyclearcs[noddcyclearcs++] = backwardarcs[i]; 989 } 990 991 /* add constraints to ensure no parallel edges */ 992 (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "cons_noparalleledges"); 993 SCIP_CALL( SCIPcreateConsBasicLinear(subscip, &noparallelcons, name, 0, NULL, NULL, 0.0, 1.0) ); 994 SCIP_CALL( SCIPaddCoefLinear(subscip, noparallelcons, forwardarcs[*narcs], 1.0) ); 995 SCIP_CALL( SCIPaddCoefLinear(subscip, noparallelcons, backwardarcs[*narcs], 1.0) ); 996 SCIP_CALL( SCIPaddCons(subscip, noparallelcons) ); 997 SCIP_CALL( SCIPreleaseCons(subscip, &noparallelcons) ); 998 999 ++*narcs; 1000 } 1001 } 1002 assert(*narcs > 0); 1003 1004 /* odd cycle property constraint */ 1005 (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "cons_oddcycle"); 1006 SCIP_CALL( SCIPcreateConsBasicXor(subscip, &oddcyclecons, name, TRUE, noddcyclearcs, oddcyclearcs) ); 1007 SCIP_CALL( SCIPaddCons(subscip, oddcyclecons) ); 1008 SCIP_CALL( SCIPreleaseCons(subscip, &oddcyclecons) ); 1009 SCIPfreeBufferArray(subscip, &oddcyclearcs); 1010 1011 /* cycle length constraint */ 1012 (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "cons_cyclelength"); 1013 SCIP_CALL( SCIPcreateConsBasicLinear(subscip, &cyclelengthcons, name, 0, NULL, NULL, 1014 (SCIP_Real) sepadata->minaggrsize, (SCIP_Real) sepadata->maxaggrsize) ); 1015 1016 for( i = 0; i < *narcs; ++i ) 1017 { 1018 SCIP_CALL( SCIPaddCoefLinear(subscip, cyclelengthcons, forwardarcs[i], 1.0) ); 1019 SCIP_CALL( SCIPaddCoefLinear(subscip, cyclelengthcons, backwardarcs[i], 1.0) ); 1020 } 1021 1022 SCIP_CALL( SCIPaddCons(subscip, cyclelengthcons) ); 1023 SCIP_CALL( SCIPreleaseCons(subscip, &cyclelengthcons) ); 1024 1025 /* create flow conservation constraints */ 1026 SCIP_CALL( SCIPallocBufferArray(subscip, &flowcons, nnodes) ); 1027 1028 for( i = 0; i < nnodes; ++i ) 1029 { 1030 (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "cons_flowconservation#%d", i); 1031 SCIP_CALL( SCIPcreateConsBasicLinear(subscip, &flowcons[i], name, 0, NULL, NULL, 0.0, 0.0) ); 1032 } 1033 1034 arcidx = 0; 1035 for( i = 0; i < nquadexprs; ++i ) 1036 { 1037 SCIP_EXPR* qterm; 1038 SCIP_Real coef; 1039 int nadjbilin; 1040 int* adjbilin; 1041 int j; 1042 1043 SCIPexprGetQuadraticQuadTerm(expr, i, &qterm, NULL, &coef, &nadjbilin, &adjbilin, NULL); 1044 1045 if( !SCIPisZero(scip, coef) ) 1046 ++arcidx; 1047 1048 for( j = 0 ; j < nadjbilin; ++j ) 1049 { 1050 SCIP_EXPR* qterm1; 1051 int pos2; 1052 1053 SCIPexprGetQuadraticBilinTerm(expr, adjbilin[j], &qterm1, NULL, NULL, &pos2, NULL); 1054 1055 /* handle qterm == qterm2 later */ 1056 if( qterm1 != qterm ) 1057 continue; 1058 1059 SCIP_CALL( SCIPaddCoefLinear(subscip, flowcons[i], forwardarcs[arcidx], 1.0) ); 1060 SCIP_CALL( SCIPaddCoefLinear(subscip, flowcons[i], backwardarcs[arcidx], -1.0) ); 1061 1062 SCIP_CALL( SCIPaddCoefLinear(subscip, flowcons[pos2], forwardarcs[arcidx], -1.0) ); 1063 SCIP_CALL( SCIPaddCoefLinear(subscip, flowcons[pos2], backwardarcs[arcidx], 1.0) ); 1064 1065 ++arcidx; 1066 } 1067 } 1068 assert(arcidx == *narcs); 1069 1070 for( i = 0; i < nnodes; ++i ) 1071 { 1072 SCIP_CALL( SCIPaddCons(subscip, flowcons[i]) ); 1073 SCIP_CALL( SCIPreleaseCons(subscip, &flowcons[i]) ); 1074 } 1075 1076 SCIPfreeBufferArray(subscip, &flowcons); 1077 1078 return SCIP_OKAY; 1079 } 1080 1081 /** fixed all arc variables (u,v) for which u or v is already in an edge-concave aggregation */ 1082 static 1083 SCIP_RETCODE updateMIP( 1084 SCIP* subscip, /**< auxiliary SCIP to search aggregations */ 1085 SCIP_NLROW* nlrow, /**< nonlinear row */ 1086 SCIP_VAR** forwardarcs, /**< forward arc variables */ 1087 SCIP_VAR** backwardarcs, /**< backward arc variables */ 1088 int* quadvar2aggr, /**< mapping of quadvars to e.c. aggr. index (< 0: in no aggr.) */ 1089 int* nedges /**< pointer to store the number of nonexcluded edges */ 1090 ) 1091 { 1092 SCIP_EXPR* expr; 1093 int nquadexprs; 1094 int arcidx; 1095 int i; 1096 1097 assert(subscip != NULL); 1098 assert(nlrow != NULL); 1099 assert(forwardarcs != NULL); 1100 assert(backwardarcs != NULL); 1101 assert(quadvar2aggr != NULL); 1102 assert(nedges != NULL); 1103 1104 SCIP_CALL( SCIPfreeTransform(subscip) ); 1105 1106 /* recompute the number of edges */ 1107 *nedges = 0; 1108 1109 expr = SCIPnlrowGetExpr(nlrow); 1110 SCIPexprGetQuadraticData(expr, NULL, NULL, NULL, NULL, &nquadexprs, NULL, NULL, NULL); 1111 1112 /* fix each arc to 0 if at least one of its nodes is contained in an e.c. aggregation */ 1113 arcidx = 0; 1114 for( i = 0; i < nquadexprs; ++i ) 1115 { 1116 SCIP_EXPR* qterm; 1117 SCIP_Real coef; 1118 int nadjbilin; 1119 int* adjbilin; 1120 int j; 1121 1122 SCIPexprGetQuadraticQuadTerm(expr, i, &qterm, NULL, &coef, &nadjbilin, &adjbilin, NULL); 1123 1124 if( !SCIPisZero(subscip, coef) ) 1125 { 1126 if( quadvar2aggr[i] != -1 ) 1127 { 1128 SCIP_CALL( SCIPchgVarUb(subscip, forwardarcs[arcidx], 0.0) ); 1129 SCIP_CALL( SCIPchgVarUb(subscip, backwardarcs[arcidx], 0.0) ); 1130 } 1131 ++arcidx; 1132 } 1133 1134 for( j = 0 ; j < nadjbilin; ++j ) 1135 { 1136 SCIP_EXPR* qterm1; 1137 int pos2; 1138 1139 SCIPexprGetQuadraticBilinTerm(expr, adjbilin[j], &qterm1, NULL, NULL, &pos2, NULL); 1140 1141 /* handle qterm == qterm2 later */ 1142 if( qterm1 != qterm ) 1143 continue; 1144 1145 if( quadvar2aggr[i] != -1 || quadvar2aggr[pos2] != -1 ) 1146 { 1147 SCIP_CALL( SCIPchgVarUb(subscip, forwardarcs[arcidx], 0.0) ); 1148 SCIP_CALL( SCIPchgVarUb(subscip, backwardarcs[arcidx], 0.0) ); 1149 } 1150 else 1151 ++*nedges; 1152 1153 ++arcidx; 1154 } 1155 } 1156 1157 return SCIP_OKAY; 1158 } 1159 1160 /** stores the best edge-concave aggregation found by the MIP model */ 1161 static 1162 SCIP_RETCODE storeAggrFromMIP( 1163 SCIP* subscip, /**< auxiliary SCIP to search aggregations */ 1164 SCIP_NLROW* nlrow, /**< nonlinear row */ 1165 SCIP_VAR** forwardarcs, /**< forward arc variables */ 1166 SCIP_VAR** backwardarcs, /**< backward arc variables */ 1167 int* quadvar2aggr, /**< mapping of quadvars to e.c. aggr. index (< 0: in no aggr.) */ 1168 int nfoundsofar /**< number of e.c. aggregation found so far */ 1169 ) 1170 { 1171 SCIP_SOL* sol; 1172 SCIP_EXPR* expr; 1173 int nquadexprs; 1174 int arcidx; 1175 int i; 1176 1177 assert(subscip != NULL); 1178 assert(nlrow != NULL); 1179 assert(forwardarcs != NULL); 1180 assert(backwardarcs != NULL); 1181 assert(quadvar2aggr != NULL); 1182 assert(nfoundsofar >= 0); 1183 assert(SCIPgetStatus(subscip) < SCIP_STATUS_INFEASIBLE); 1184 assert(SCIPgetNSols(subscip) > 0); 1185 1186 sol = SCIPgetBestSol(subscip); 1187 assert(sol != NULL); 1188 1189 expr = SCIPnlrowGetExpr(nlrow); 1190 SCIPexprGetQuadraticData(expr, NULL, NULL, NULL, NULL, &nquadexprs, NULL, NULL, NULL); 1191 1192 /* fix each arc to 0 if at least one of its nodes is contained in an e.c. aggregation */ 1193 arcidx = 0; 1194 for( i = 0; i < nquadexprs; ++i ) 1195 { 1196 SCIP_EXPR* qterm; 1197 SCIP_Real coef; 1198 int nadjbilin; 1199 int* adjbilin; 1200 int j; 1201 1202 SCIPexprGetQuadraticQuadTerm(expr, i, &qterm, NULL, &coef, &nadjbilin, &adjbilin, NULL); 1203 1204 if( !SCIPisZero(subscip, coef) ) 1205 { 1206 if( SCIPisGT(subscip, SCIPgetSolVal(subscip, sol, forwardarcs[arcidx]), 0.5) || 1207 SCIPisGT(subscip, SCIPgetSolVal(subscip, sol, backwardarcs[arcidx]), 0.5) ) 1208 { 1209 assert(quadvar2aggr[i] == -1 || quadvar2aggr[i] == nfoundsofar); 1210 quadvar2aggr[i] = nfoundsofar; 1211 } 1212 1213 ++arcidx; 1214 } 1215 1216 for( j = 0; j < nadjbilin; ++j ) 1217 { 1218 SCIP_EXPR* qterm1; 1219 int pos2; 1220 1221 SCIPexprGetQuadraticBilinTerm(expr, adjbilin[j], &qterm1, NULL, NULL, &pos2, NULL); 1222 1223 /* handle qterm == qterm2 later */ 1224 if( qterm1 != qterm ) 1225 continue; 1226 1227 if( SCIPisGT(subscip, SCIPgetSolVal(subscip, sol, forwardarcs[arcidx]), 0.5) || 1228 SCIPisGT(subscip, SCIPgetSolVal(subscip, sol, backwardarcs[arcidx]), 0.5) ) 1229 { 1230 assert(quadvar2aggr[i] == -1 || quadvar2aggr[i] == nfoundsofar); 1231 assert(quadvar2aggr[pos2] == -1 || quadvar2aggr[pos2] == nfoundsofar); 1232 1233 quadvar2aggr[i] = nfoundsofar; 1234 quadvar2aggr[pos2] = nfoundsofar; 1235 } 1236 1237 ++arcidx; 1238 } 1239 } 1240 1241 return SCIP_OKAY; 1242 } 1243 1244 /** searches for edge-concave aggregations with a MIP model based on binary flow variables */ 1245 static 1246 SCIP_RETCODE searchEcAggrWithMIP( 1247 SCIP* subscip, /**< SCIP data structure */ 1248 SCIP_Real timelimit, /**< time limit to solve the MIP */ 1249 int nedges, /**< number of nonexcluded undirected edges */ 1250 SCIP_Bool* aggrleft, /**< pointer to store if there might be a left aggregation */ 1251 SCIP_Bool* found /**< pointer to store if we have found an aggregation */ 1252 ) 1253 { 1254 assert(subscip != NULL); 1255 assert(aggrleft != NULL); 1256 assert(found != NULL); 1257 assert(nedges >= 0); 1258 1259 *aggrleft = TRUE; 1260 *found = FALSE; 1261 1262 if( SCIPisLE(subscip, timelimit, 0.0) ) 1263 return SCIP_OKAY; 1264 1265 /* set working limits */ 1266 SCIP_CALL( SCIPsetRealParam(subscip, "limits/time", timelimit) ); 1267 SCIP_CALL( SCIPsetLongintParam(subscip, "limits/totalnodes", SUBSCIP_NODELIMIT) ); 1268 1269 /* set heuristics to aggressive */ 1270 SCIP_CALL( SCIPsetHeuristics(subscip, SCIP_PARAMSETTING_AGGRESSIVE, TRUE) ); 1271 1272 /* disable output to console in optimized mode, enable in SCIP's debug mode */ 1273 #ifdef SCIP_DEBUG 1274 SCIP_CALL( SCIPsetIntParam(subscip, "display/verblevel", 5) ); 1275 SCIP_CALL( SCIPsetIntParam(subscip, "display/freq", 1) ); 1276 #else 1277 SCIP_CALL( SCIPsetIntParam(subscip, "display/verblevel", 0) ); 1278 #endif 1279 1280 SCIP_CALL( SCIPsolve(subscip) ); 1281 1282 /* no more aggregation left if the MIP is infeasible */ 1283 if( SCIPgetStatus(subscip) >= SCIP_STATUS_INFEASIBLE ) 1284 { 1285 *found = FALSE; 1286 *aggrleft = FALSE; 1287 return SCIP_OKAY; 1288 } 1289 1290 if( SCIPgetNSols(subscip) > 0 ) 1291 { 1292 *found = TRUE; 1293 *aggrleft = TRUE; 1294 1295 #ifdef SCIP_DEBUG 1296 if( SCIPgetNSols(subscip) > 0 ) 1297 { 1298 SCIP_CALL( SCIPprintSol(subscip, SCIPgetBestSol(subscip), NULL , FALSE) ); 1299 } 1300 #endif 1301 } 1302 1303 return SCIP_OKAY; 1304 } 1305 1306 /** creates a tclique graph from a given nonlinear row 1307 * 1308 * SCIP's clique code can only handle integer node weights; all node weights are scaled by a factor of 100; since the 1309 * clique code ignores nodes with weight of zero, we add an offset of 100 to each weight 1310 */ 1311 static 1312 SCIP_RETCODE createTcliqueGraph( 1313 SCIP_NLROW* nlrow, /**< nonlinear row */ 1314 TCLIQUE_GRAPH** graph, /**< TCLIQUE graph structure */ 1315 SCIP_Real* nodeweights /**< weights for each quadratic variable (nodes in the graph) */ 1316 ) 1317 { 1318 SCIP_EXPR* expr; 1319 int nquadexprs; 1320 int i; 1321 1322 assert(graph != NULL); 1323 assert(nlrow != NULL); 1324 1325 /* create the tclique graph */ 1326 if( !tcliqueCreate(graph) ) 1327 { 1328 SCIPerrorMessage("could not create clique graph\n"); 1329 return SCIP_ERROR; 1330 } 1331 1332 expr = SCIPnlrowGetExpr(nlrow); 1333 SCIPexprGetQuadraticData(expr, NULL, NULL, NULL, NULL, &nquadexprs, NULL, NULL, NULL); 1334 1335 /* add all nodes to the tclique graph */ 1336 for( i = 0; i < nquadexprs; ++i ) 1337 { 1338 int nodeweight; 1339 1340 /* note: clique code can only handle integer weights */ 1341 nodeweight = 100 + (int)(100 * nodeweights[i]); 1342 /* SCIPdebugMsg(scip, "%d (%s): nodeweight %d \n", i, SCIPvarGetName(SCIPnlrowGetQuadVars(nlrow)[i]), nodeweight); */ 1343 1344 if( !tcliqueAddNode(*graph, i, nodeweight) ) 1345 { 1346 SCIPerrorMessage("could not add node to clique graph\n"); 1347 return SCIP_ERROR; 1348 } 1349 } 1350 1351 /* add all edges */ 1352 for( i = 0; i < nquadexprs; ++i ) 1353 { 1354 SCIP_EXPR* qterm; 1355 int nadjbilin; 1356 int* adjbilin; 1357 int j; 1358 1359 SCIPexprGetQuadraticQuadTerm(expr, i, &qterm, NULL, NULL, &nadjbilin, &adjbilin, NULL); 1360 1361 for( j = 0; j < nadjbilin; ++j ) 1362 { 1363 SCIP_EXPR* qterm1; 1364 SCIP_EXPR* qterm2; 1365 int pos2; 1366 1367 SCIPexprGetQuadraticBilinTerm(expr, adjbilin[j], &qterm1, &qterm2, NULL, &pos2, NULL); 1368 1369 /* handle qterm == qterm2 later */ 1370 if( qterm1 != qterm ) 1371 continue; 1372 1373 #ifdef SCIP_DEBUG_DETAILED 1374 SCIPdebugMessage(" add edge (%d, %d) = (%s,%s) to tclique graph\n", 1375 SCIPvarGetIndex(SCIPgetVarExprVar(qterm1)), SCIPvarGetIndex(SCIPgetVarExprVar(qterm2)), 1376 SCIPvarGetName(SCIPgetVarExprVar(qterm1)), SCIPvarGetName(SCIPgetVarExprVar(qterm2))); 1377 #endif 1378 1379 if( !tcliqueAddEdge(*graph, i, pos2) ) 1380 { 1381 SCIPerrorMessage("could not add edge to clique graph\n"); 1382 return SCIP_ERROR; 1383 } 1384 } 1385 } 1386 1387 /* flush the clique graph */ 1388 if( !tcliqueFlush(*graph) ) 1389 { 1390 SCIPerrorMessage("could not flush the clique graph\n"); 1391 return SCIP_ERROR; 1392 } 1393 1394 return SCIP_OKAY; 1395 } 1396 1397 /** searches for edge-concave aggregations by computing cliques in the graph representation of a given nonlinear row 1398 * 1399 * update graph, compute clique, store clique; after computing a clique we heuristically check if the clique contains 1400 * at least one good cycle 1401 */ 1402 static 1403 SCIP_RETCODE searchEcAggrWithCliques( 1404 SCIP* scip, /**< SCIP data structure */ 1405 TCLIQUE_GRAPH* graph, /**< TCLIQUE graph structure */ 1406 SCIP_SEPADATA* sepadata, /**< separator data */ 1407 SCIP_NLROW* nlrow, /**< nonlinear row */ 1408 int* quadvar2aggr, /**< mapping of quadvars to e.c. aggr. index (< 0: in no aggr.) */ 1409 int nfoundsofar, /**< number of e.c. aggregation found so far */ 1410 SCIP_Bool rhsaggr, /**< consider nonlinear row aggregation for g(x) <= rhs (TRUE) or 1411 * lhs <= g(x) (FALSE) */ 1412 SCIP_Bool* foundaggr, /**< pointer to store if we have found an aggregation */ 1413 SCIP_Bool* foundclique /**< pointer to store if we have found a clique */ 1414 ) 1415 { 1416 SCIP_HASHMAP* cliquemap; 1417 TCLIQUE_STATUS status; 1418 SCIP_EXPR* expr; 1419 int nquadexprs; 1420 int* maxcliquenodes; 1421 int* degrees; 1422 int nmaxcliquenodes; 1423 int maxcliqueweight; 1424 int noddcycleedges; 1425 int ntwodegrees; 1426 int aggrsize; 1427 int i; 1428 1429 assert(graph != NULL); 1430 assert(nfoundsofar >= 0); 1431 assert(foundaggr != NULL); 1432 assert(foundclique != NULL); 1433 1434 cliquemap = NULL; 1435 *foundaggr = FALSE; 1436 *foundclique = FALSE; 1437 1438 expr = SCIPnlrowGetExpr(nlrow); 1439 SCIPexprGetQuadraticData(expr, NULL, NULL, NULL, NULL, &nquadexprs, NULL, NULL, NULL); 1440 assert(nquadexprs == tcliqueGetNNodes(graph)); 1441 1442 /* exclude all nodes which are already in an edge-concave aggregation (no flush is needed) */ 1443 for( i = 0; i < nquadexprs; ++i ) 1444 { 1445 if( quadvar2aggr[i] != -1 ) 1446 { 1447 SCIPdebugMsg(scip, "exclude node %d from clique graph\n", i); 1448 tcliqueChangeWeight(graph, i, 0); 1449 } 1450 } 1451 1452 SCIP_CALL( SCIPallocBufferArray(scip, &maxcliquenodes, nquadexprs) ); 1453 1454 /* solve clique problem */ 1455 tcliqueMaxClique(tcliqueGetNNodes, tcliqueGetWeights, tcliqueIsEdge, tcliqueSelectAdjnodes, graph, NULL, NULL, 1456 maxcliquenodes, &nmaxcliquenodes, &maxcliqueweight, CLIQUE_MAXFIRSTNODEWEIGHT, CLIQUE_MINWEIGHT, 1457 CLIQUE_MAXNTREENODES, CLIQUE_BACKTRACKFREQ, 0, -1, NULL, &status); 1458 1459 if( status != TCLIQUE_OPTIMAL || nmaxcliquenodes < sepadata->minaggrsize ) 1460 goto TERMINATE; 1461 1462 *foundclique = TRUE; 1463 aggrsize = MIN(sepadata->maxaggrsize, nmaxcliquenodes); 1464 SCIP_CALL( SCIPhashmapCreate(&cliquemap, SCIPblkmem(scip), aggrsize) ); 1465 1466 for( i = 0; i < aggrsize; ++i ) 1467 { 1468 SCIP_CALL( SCIPhashmapInsertInt(cliquemap, (void*) (size_t) maxcliquenodes[i], 0) ); /*lint !e571*/ 1469 } 1470 1471 /* count the degree of good cycle edges for each node in the clique */ 1472 SCIP_CALL( SCIPallocBufferArray(scip, °rees, aggrsize) ); 1473 BMSclearMemoryArray(degrees, aggrsize); 1474 ntwodegrees = 0; 1475 1476 /* count the number of positive or negative edges (depending on <= rhs or >= lhs) */ 1477 noddcycleedges = 0; 1478 for( i = 0; i < nquadexprs; ++i ) 1479 { 1480 SCIP_Bool isoddcycleedge; 1481 SCIP_EXPR* qterm; 1482 SCIP_Real coef; 1483 int nadjbilin; 1484 int* adjbilin; 1485 int j; 1486 1487 SCIPexprGetQuadraticQuadTerm(expr, i, &qterm, NULL, &coef, &nadjbilin, &adjbilin, NULL); 1488 1489 isoddcycleedge = (rhsaggr && SCIPisPositive(scip, coef)) || (!rhsaggr && SCIPisNegative(scip, coef)); 1490 1491 if( isoddcycleedge && SCIPhashmapExists(cliquemap, (void*) (size_t) i) ) 1492 { 1493 ++noddcycleedges; 1494 ++degrees[i]; 1495 } 1496 1497 for( j = 0; j < nadjbilin; ++j ) 1498 { 1499 SCIP_EXPR* qterm1; 1500 SCIP_EXPR* qterm2; 1501 int pos2; 1502 1503 SCIPexprGetQuadraticBilinTerm(expr, adjbilin[j], &qterm1, &qterm2, &coef, &pos2, NULL); 1504 1505 /* handle qterm == qterm2 later */ 1506 if( qterm1 != qterm ) 1507 continue; 1508 1509 isoddcycleedge = (rhsaggr && SCIPisPositive(scip, coef)) || (!rhsaggr && SCIPisNegative(scip, coef)); 1510 1511 if( isoddcycleedge 1512 && SCIPhashmapExists(cliquemap, (void*) (size_t) i) 1513 && SCIPhashmapExists(cliquemap, (void*) (size_t) pos2) ) 1514 { 1515 ++noddcycleedges; 1516 ++degrees[i]; 1517 ++degrees[pos2]; 1518 } 1519 } 1520 } 1521 1522 /* count the number of nodes with exactly two incident odd cycle edges */ 1523 for( i = 0; i < aggrsize; ++i ) 1524 if( degrees[i] == 2 ) 1525 ++ntwodegrees; 1526 1527 /* check cases for which we are sure that there are no good cycles in the clique */ 1528 if( noddcycleedges == 0 || (aggrsize == 3 && noddcycleedges == 2) || (aggrsize == 4 && ntwodegrees == 4) ) 1529 *foundaggr = FALSE; 1530 else 1531 *foundaggr = TRUE; 1532 1533 /* add the found clique as an edge-concave aggregation or exclude the nodes from the remaining search */ 1534 for( i = 0; i < aggrsize; ++i ) 1535 { 1536 quadvar2aggr[ maxcliquenodes[i] ] = *foundaggr ? nfoundsofar : -2; 1537 SCIPdebugMsg(scip, "%s %d\n", *foundaggr ? "aggregate node: " : "exclude node: ", maxcliquenodes[i]); 1538 } 1539 1540 SCIPfreeBufferArray(scip, °rees); 1541 1542 TERMINATE: 1543 if( cliquemap != NULL ) 1544 SCIPhashmapFree(&cliquemap); 1545 SCIPfreeBufferArray(scip, &maxcliquenodes); 1546 1547 return SCIP_OKAY; 1548 } 1549 1550 /** helper function for searchEcAggr() */ 1551 static 1552 SCIP_RETCODE doSeachEcAggr( 1553 SCIP* scip, /**< SCIP data structure */ 1554 SCIP* subscip, /**< sub-SCIP data structure */ 1555 SCIP_SEPADATA* sepadata, /**< separator data */ 1556 SCIP_NLROW* nlrow, /**< nonlinear row */ 1557 SCIP_SOL* sol, /**< current solution (might be NULL) */ 1558 SCIP_Bool rhsaggr, /**< consider nonlinear row aggregation for g(x) <= rhs (TRUE) or g(x) >= lhs (FALSE) */ 1559 int* quadvar2aggr, /**< array to store for each quadratic variable in which edge-concave 1560 * aggregation it is stored (< 0: in no aggregation); size has to be at 1561 * least SCIPnlrowGetNQuadVars(nlrow) */ 1562 int* nfound /**< pointer to store the number of found e.c. aggregations */ 1563 ) 1564 { 1565 TCLIQUE_GRAPH* graph = NULL; 1566 SCIP_EXPR* expr; 1567 SCIP_VAR** forwardarcs; 1568 SCIP_VAR** backwardarcs; 1569 SCIP_Real* nodeweights; 1570 SCIP_Real timelimit; 1571 SCIP_RETCODE retcode; 1572 int nunsucces = 0; 1573 int nedges = 0; 1574 int narcs; 1575 int nquadvars; 1576 int nbilinexprs; 1577 int i; 1578 1579 assert(subscip != NULL); 1580 assert(quadvar2aggr != NULL); 1581 assert(nfound != NULL); 1582 1583 expr = SCIPnlrowGetExpr(nlrow); 1584 SCIPexprGetQuadraticData(expr, NULL, NULL, NULL, NULL, &nquadvars, &nbilinexprs, NULL, NULL); 1585 1586 retcode = SCIP_OKAY; 1587 *nfound = 0; 1588 1589 /* arrays to store all arc variables of the MIP model; note that we introduce variables even for loops in the graph 1590 * to have an easy mapping from the edges of the graph to the quadratic elements 1591 * nquadvars + nbilinexprs is an upper bound on the actual number of square and bilinear terms 1592 */ 1593 SCIP_CALL( SCIPallocBufferArray(scip, &nodeweights, nquadvars) ); 1594 SCIP_CALL( SCIPallocBufferArray(scip, &forwardarcs, nquadvars + nbilinexprs) ); 1595 SCIP_CALL( SCIPallocBufferArray(scip, &backwardarcs, nquadvars + nbilinexprs) ); 1596 1597 /* initialize mapping from quadvars to e.c. aggregation index (-1: quadvar is in no aggregation); compute node 1598 * weights 1599 */ 1600 for( i = 0; i < nquadvars; ++i ) 1601 { 1602 SCIP_EXPR* qterm; 1603 SCIP_VAR* var; 1604 1605 SCIPexprGetQuadraticQuadTerm(expr, i, &qterm, NULL, NULL, NULL, NULL, NULL); 1606 assert(SCIPisExprVar(scip, qterm)); 1607 var = SCIPgetVarExprVar(qterm); 1608 1609 quadvar2aggr[i] = -1; 1610 nodeweights[i] = phi(scip, SCIPgetSolVal(scip, sol, var), SCIPvarGetLbLocal(var), SCIPvarGetUbLocal(var)); 1611 SCIPdebugMsg(scip, "%s = %e (%e in [%e, %e])\n", SCIPvarGetName(var), nodeweights[i], SCIPgetSolVal(scip, sol, var), 1612 SCIPvarGetLbLocal(var), SCIPvarGetUbLocal(var)); 1613 } 1614 1615 SCIP_CALL( createMIP(scip, subscip, sepadata, nlrow, rhsaggr, forwardarcs, backwardarcs, nodeweights, &nedges, &narcs) ); 1616 assert(nedges >= 0); 1617 assert(narcs > 0); 1618 SCIPdebugMsg(scip, "nedges (without loops) = %d\n", nedges); 1619 SCIPdebugMsg(scip, "narcs (number of quadratic terms) = %d\n", narcs); 1620 1621 SCIP_CALL( SCIPgetRealParam(scip, "limits/time", &timelimit) ); 1622 1623 /* main loop to search for edge-concave aggregations */ 1624 while( !SCIPisStopped(scip) ) 1625 { 1626 SCIP_Bool aggrleft; 1627 SCIP_Bool found; 1628 1629 SCIPdebugMsg(scip, "#remaining edges = %d\n", nedges); 1630 1631 /* not enough edges left */ 1632 if( nedges < sepadata->minaggrsize ) 1633 break; 1634 1635 /* check whether there is enough time left; update the remaining time */ 1636 if( !SCIPisInfinity(scip, timelimit) ) 1637 { 1638 timelimit -= SCIPgetSolvingTime(scip); 1639 if( timelimit <= 0.0 ) 1640 { 1641 SCIPdebugMsg(scip, "skip aggregation search since no time left\n"); 1642 goto TERMINATE; 1643 } 1644 } 1645 1646 /* 1.a - search for edge-concave aggregation with the help of the MIP model */ 1647 SCIP_CALL( searchEcAggrWithMIP(subscip, timelimit, nedges, &aggrleft, &found) ); 1648 1649 /* 1.b - there are no more edge-concave aggregations left */ 1650 if( !aggrleft ) 1651 { 1652 SCIPdebugMsg(scip, "no more aggregation left\n"); 1653 break; 1654 } 1655 1656 if( found ) 1657 { 1658 SCIP_CALL( storeAggrFromMIP(subscip, nlrow, forwardarcs, backwardarcs, quadvar2aggr, *nfound) ); 1659 ++(*nfound); 1660 nunsucces = 0; 1661 } 1662 /* try to find an edge-concave aggregation by computing cliques */ 1663 else 1664 { 1665 SCIP_Bool foundaggr; 1666 SCIP_Bool foundclique; 1667 1668 ++nunsucces; 1669 1670 /* create graph if necessary */ 1671 if( graph == NULL ) 1672 { 1673 SCIP_CALL_TERMINATE( retcode, createTcliqueGraph(nlrow, &graph, nodeweights), TERMINATE ); 1674 } 1675 1676 /* 2.a - search and store a single edge-concave aggregation by computing a clique with a good cycle */ 1677 SCIP_CALL_FINALLY( searchEcAggrWithCliques(scip, graph, sepadata, nlrow, quadvar2aggr, *nfound, rhsaggr, 1678 &foundaggr, &foundclique), tcliqueFree(&graph) ); 1679 1680 if( foundaggr ) 1681 { 1682 assert(foundclique); 1683 ++(*nfound); 1684 nunsucces = 0; 1685 } 1686 else 1687 ++nunsucces; 1688 1689 /* 2.b - no clique of at least minaggrsize size found */ 1690 if( !foundclique ) 1691 { 1692 assert(!foundaggr); 1693 SCIPdebugMsg(scip, "did not find a clique to exclude -> leave aggregation search\n"); 1694 break; 1695 } 1696 } 1697 1698 /* leave the algorithm if we did not find something for maxstallrounds many iterations */ 1699 if( nunsucces >= sepadata->maxstallrounds && *nfound == 0 ) 1700 { 1701 SCIPdebugMsg(scip, "did not find an e.c. aggregation for %d iterations\n", nunsucces); 1702 break; 1703 } 1704 1705 /* exclude all edges used in the last aggregation and nodes found in the clique solution */ 1706 SCIP_CALL_FINALLY( updateMIP(subscip, nlrow, forwardarcs, backwardarcs, quadvar2aggr, &nedges), tcliqueFree(&graph) ); 1707 } 1708 1709 TERMINATE: 1710 1711 #ifdef SCIP_DEBUG 1712 SCIPdebugMsg(scip, "aggregations found:\n"); 1713 for( i = 0; i < nquadvars; ++i ) 1714 { 1715 SCIPdebugMsg(scip, " %d in %d\n", i, quadvar2aggr[i]); 1716 } 1717 #endif 1718 1719 /* free clique graph */ 1720 if( graph != NULL ) 1721 tcliqueFree(&graph); 1722 1723 /* free sub-SCIP */ 1724 for( i = 0; i < narcs; ++i ) 1725 { 1726 SCIP_CALL( SCIPreleaseVar(subscip, &forwardarcs[i]) ); 1727 SCIP_CALL( SCIPreleaseVar(subscip, &backwardarcs[i]) ); 1728 } 1729 1730 SCIPfreeBufferArray(scip, &backwardarcs); 1731 SCIPfreeBufferArray(scip, &forwardarcs); 1732 SCIPfreeBufferArray(scip, &nodeweights); 1733 1734 return retcode; 1735 } 1736 1737 /** computes a partitioning into edge-concave aggregations for a given (quadratic) nonlinear row 1738 * 1739 * Each aggregation has to contain a cycle with an odd number of positive weighted edges (good cycles) in the corresponding graph representation. 1740 * For this we use the following algorithm: 1741 * -# use a MIP model based on binary flow variables to compute good cycles and store the implied subgraphs as an e.c. aggr. 1742 * -# if we find a good cycle, store the implied subgraph, delete it from the graph representation and go to 1) 1743 * -# if the MIP model is infeasible (there are no good cycles), STOP 1744 * -# we compute a large clique C if the MIP model fails (because of working limits, etc) 1745 * -# if we find a good cycle in C, store the implied subgraph of C, delete it from the graph representation and go to 1) 1746 * -# if C is not large enough, STOP 1747 */ 1748 static 1749 SCIP_RETCODE searchEcAggr( 1750 SCIP* scip, /**< SCIP data structure */ 1751 SCIP_SEPADATA* sepadata, /**< separator data */ 1752 SCIP_NLROW* nlrow, /**< nonlinear row */ 1753 SCIP_SOL* sol, /**< current solution (might be NULL) */ 1754 SCIP_Bool rhsaggr, /**< consider nonlinear row aggregation for g(x) <= rhs (TRUE) or g(x) >= lhs (FALSE) */ 1755 int* quadvar2aggr, /**< array to store for each quadratic variable in which edge-concave 1756 * aggregation it is stored (< 0: in no aggregation); size has to be at 1757 * least SCIPnlrowGetNQuadVars(nlrow) */ 1758 int* nfound /**< pointer to store the number of found e.c. aggregations */ 1759 ) 1760 { 1761 SCIP* subscip; 1762 SCIP_RETCODE retcode; 1763 1764 /* create and set up a sub-SCIP */ 1765 SCIP_CALL_FINALLY( SCIPcreate(&subscip), (void)SCIPfree(&subscip) ); 1766 1767 retcode = doSeachEcAggr(scip, subscip, sepadata, nlrow, sol, rhsaggr, quadvar2aggr, nfound); 1768 1769 SCIP_CALL( SCIPfree(&subscip) ); 1770 SCIP_CALL( retcode ); 1771 1772 return SCIP_OKAY; 1773 } 1774 1775 /** returns whether a given nonlinear row can be used to compute edge-concave aggregations for which their convex 1776 * envelope could dominate the termwise bilinear relaxation 1777 * 1778 * This is the case if there exists at least one cycle with 1779 * an odd number of positive edges in the corresponding graph representation of the nonlinear row. 1780 */ 1781 static 1782 SCIP_RETCODE isCandidate( 1783 SCIP* scip, /**< SCIP data structure */ 1784 SCIP_SEPADATA* sepadata, /**< separator data */ 1785 SCIP_NLROW* nlrow, /**< nonlinear row representation of a nonlinear constraint */ 1786 SCIP_Bool* rhscandidate, /**< pointer to store if we should compute edge-concave aggregations for 1787 * the <= rhs case */ 1788 SCIP_Bool* lhscandidate /**< pointer to store if we should compute edge-concave aggregations for 1789 * the >= lhs case */ 1790 ) 1791 { 1792 SCIP_EXPR* expr = NULL; 1793 SCIP_Bool takerow = FALSE; 1794 int nquadvars = 0; 1795 int* degrees; 1796 int ninterestingnodes; 1797 int nposedges; 1798 int nnegedges; 1799 int i; 1800 1801 assert(rhscandidate != NULL); 1802 assert(lhscandidate != NULL); 1803 1804 *rhscandidate = TRUE; 1805 *lhscandidate = TRUE; 1806 1807 /* check whether nlrow is in the NLP, is quadratic in variables, and there are enough quadratic variables */ 1808 if( SCIPnlrowIsInNLP(nlrow) && SCIPnlrowGetExpr(nlrow) != NULL ) 1809 { 1810 expr = SCIPnlrowGetExpr(nlrow); 1811 SCIP_CALL( SCIPcheckExprQuadratic(scip, expr, &takerow) ); 1812 } 1813 if( takerow ) 1814 takerow = SCIPexprAreQuadraticExprsVariables(expr); 1815 if( takerow ) 1816 { 1817 SCIPexprGetQuadraticData(expr, NULL, NULL, NULL, NULL, &nquadvars, NULL, NULL, NULL); 1818 takerow = nquadvars >= sepadata->minaggrsize; 1819 } 1820 if( !takerow ) 1821 { 1822 *rhscandidate = FALSE; 1823 *lhscandidate = FALSE; 1824 return SCIP_OKAY; 1825 } 1826 1827 /* check for infinite rhs or lhs */ 1828 if( SCIPisInfinity(scip, REALABS(SCIPnlrowGetRhs(nlrow))) ) 1829 *rhscandidate = FALSE; 1830 if( SCIPisInfinity(scip, REALABS(SCIPnlrowGetLhs(nlrow))) ) 1831 *lhscandidate = FALSE; 1832 1833 SCIP_CALL( SCIPallocClearBufferArray(scip, °rees, nquadvars) ); 1834 1835 ninterestingnodes = 0; 1836 nposedges = 0; 1837 nnegedges = 0; 1838 1839 for( i = 0; i < nquadvars; ++i ) 1840 { 1841 SCIP_EXPR* qterm; 1842 SCIP_VAR* var1; 1843 int nadjbilin; 1844 int* adjbilin; 1845 int j; 1846 1847 SCIPexprGetQuadraticQuadTerm(expr, i, &qterm, NULL, NULL, &nadjbilin, &adjbilin, NULL); 1848 assert(SCIPisExprVar(scip, qterm)); 1849 1850 var1 = SCIPgetVarExprVar(qterm); 1851 1852 /* do not consider global fixed variables */ 1853 if( SCIPisEQ(scip, SCIPvarGetLbGlobal(var1), SCIPvarGetUbGlobal(var1)) ) 1854 continue; 1855 1856 for( j = 0; j < nadjbilin; ++j ) 1857 { 1858 SCIP_EXPR* qterm1; 1859 SCIP_EXPR* qterm2; 1860 SCIP_VAR* var2; 1861 SCIP_Real coef; 1862 int pos2; 1863 1864 SCIPexprGetQuadraticBilinTerm(expr, adjbilin[j], &qterm1, &qterm2, &coef, &pos2, NULL); 1865 1866 if( qterm1 != qterm ) 1867 continue; 1868 1869 var2 = SCIPgetVarExprVar(qterm2); 1870 1871 /* do not consider loops or global fixed variables */ 1872 if( SCIPisEQ(scip, SCIPvarGetLbGlobal(var2), SCIPvarGetUbGlobal(var2)) ) 1873 continue; 1874 1875 ++degrees[i]; 1876 ++degrees[pos2]; 1877 1878 /* count the number of nodes with a degree of at least 2 */ 1879 if( degrees[i] == 2 ) 1880 ++ninterestingnodes; 1881 if( degrees[pos2] == 2 ) 1882 ++ninterestingnodes; 1883 1884 nposedges += SCIPisPositive(scip, coef) ? 1 : 0; 1885 nnegedges += SCIPisNegative(scip, coef) ? 1 : 0; 1886 } 1887 } 1888 1889 SCIPfreeBufferArray(scip, °rees); 1890 1891 SCIPdebugMsg(scip, "nlrow contains: %d edges\n", nposedges + nnegedges); 1892 1893 /* too many edges, too few edges, or to few nodes with degree at least 2 in the graph */ 1894 if( nposedges + nnegedges > sepadata->maxbilinterms || nposedges + nnegedges < sepadata->minaggrsize 1895 || ninterestingnodes < sepadata->minaggrsize ) 1896 { 1897 *rhscandidate = FALSE; 1898 *lhscandidate = FALSE; 1899 return SCIP_OKAY; 1900 } 1901 1902 /* check if there are enough positive/negative edges; for a 3-clique there has to be an odd number of those edges */ 1903 if( nposedges == 0 || (nposedges + nnegedges == 3 && (nposedges % 2) == 0) ) 1904 *rhscandidate = FALSE; 1905 if( nnegedges == 0 || (nposedges + nnegedges == 3 && (nnegedges % 2) == 0) ) 1906 *lhscandidate = FALSE; 1907 1908 return SCIP_OKAY; 1909 } 1910 1911 /** finds and stores edge-concave aggregations for a given nonlinear row */ 1912 static 1913 SCIP_RETCODE findAndStoreEcAggregations( 1914 SCIP* scip, /**< SCIP data structure */ 1915 SCIP_SEPADATA* sepadata, /**< separator data */ 1916 SCIP_NLROW* nlrow, /**< nonlinear row */ 1917 SCIP_SOL* sol /**< current solution (might be NULL) */ 1918 ) 1919 { 1920 int nquadvars; 1921 int* quadvar2aggr; 1922 SCIP_Bool rhscandidate; 1923 SCIP_Bool lhscandidate; 1924 1925 assert(scip != NULL); 1926 assert(nlrow != NULL); 1927 assert(sepadata != NULL); 1928 1929 #ifdef SCIP_DEBUG 1930 SCIPdebugMsg(scip, "search for edge-concave aggregation for the nonlinear row: \n"); 1931 SCIP_CALL( SCIPprintNlRow(scip, nlrow, NULL) ); 1932 #endif 1933 1934 /* check obvious conditions for existing cycles with an odd number of positive/negative edges */ 1935 SCIP_CALL( isCandidate(scip, sepadata, nlrow, &rhscandidate, &lhscandidate) ); 1936 SCIPdebugMsg(scip, "rhs candidate = %u lhs candidate = %u\n", rhscandidate, lhscandidate); 1937 1938 if( !rhscandidate && !lhscandidate ) 1939 return SCIP_OKAY; 1940 1941 SCIPexprGetQuadraticData(SCIPnlrowGetExpr(nlrow), NULL, NULL, NULL, NULL, &nquadvars, NULL, NULL, NULL); 1942 SCIP_CALL( SCIPallocBufferArray(scip, &quadvar2aggr, nquadvars) ); /*lint !e705*/ 1943 1944 /* search for edge-concave aggregations (consider <= rhs) */ 1945 if( rhscandidate ) 1946 { 1947 SCIP_NLROWAGGR* nlrowaggr; 1948 int nfound; 1949 1950 assert(!SCIPisInfinity(scip, REALABS(SCIPnlrowGetRhs(nlrow)))); 1951 1952 SCIPdebugMsg(scip, "consider <= rhs\n"); 1953 SCIP_CALL( searchEcAggr(scip, sepadata, nlrow, sol, TRUE, quadvar2aggr, &nfound) ); 1954 1955 if( nfound > 0 ) 1956 { 1957 SCIP_CALL( nlrowaggrCreate(scip, nlrow, &nlrowaggr, quadvar2aggr, nfound, TRUE) ); 1958 assert(nlrow != NULL); 1959 SCIPdebug(nlrowaggrPrint(scip, nlrowaggr)); 1960 SCIP_CALL( sepadataAddNlrowaggr(scip, sepadata, nlrowaggr) ); 1961 } 1962 } 1963 1964 /* search for edge-concave aggregations (consider <= lhs) */ 1965 if( lhscandidate ) 1966 { 1967 SCIP_NLROWAGGR* nlrowaggr; 1968 int nfound; 1969 1970 assert(!SCIPisInfinity(scip, REALABS(SCIPnlrowGetLhs(nlrow)))); 1971 1972 SCIPdebugMsg(scip, "consider >= lhs\n"); 1973 SCIP_CALL( searchEcAggr(scip, sepadata, nlrow, sol, FALSE, quadvar2aggr, &nfound) ); 1974 1975 if( nfound > 0 ) 1976 { 1977 SCIP_CALL( nlrowaggrCreate(scip, nlrow, &nlrowaggr, quadvar2aggr, nfound, FALSE) ); 1978 assert(nlrow != NULL); 1979 SCIPdebug(nlrowaggrPrint(scip, nlrowaggr)); 1980 SCIP_CALL( sepadataAddNlrowaggr(scip, sepadata, nlrowaggr) ); 1981 } 1982 } 1983 1984 SCIPfreeBufferArray(scip, &quadvar2aggr); 1985 return SCIP_OKAY; 1986 } 1987 1988 /* 1989 * methods to compute edge-concave cuts 1990 */ 1991 1992 #ifdef SCIP_DEBUG 1993 /** prints a given facet (candidate) */ 1994 static 1995 void printFacet( 1996 SCIP* scip, /**< SCIP data structure */ 1997 SCIP_VAR** vars, /**< variables contained in the edge-concave aggregation */ 1998 int nvars, /**< number of variables contained in the edge-concave aggregation */ 1999 SCIP_Real* facet, /**< current facet candidate */ 2000 SCIP_Real facetval /**< facet evaluated at the current solution */ 2001 ) 2002 { 2003 int i; 2004 2005 SCIPdebugMsg(scip, "print facet (val=%e): ", facetval); 2006 for( i = 0; i < nvars; ++i ) 2007 SCIPdebugMsgPrint(scip, "%e %s + ", facet[i], SCIPvarGetName(vars[i])); 2008 SCIPdebugMsgPrint(scip, "%e\n", facet[nvars]); 2009 } 2010 #endif 2011 2012 /** checks if a facet is really an underestimate for all corners of the domain [l,u] 2013 * 2014 * Because of numerics it can happen that a facet violates a corner of the domain. 2015 * To make the facet valid we subtract the maximum violation from the constant part of the facet. 2016 */ 2017 static 2018 SCIP_Bool checkRikun( 2019 SCIP* scip, /**< SCIP data structure */ 2020 SCIP_ECAGGR* ecaggr, /**< edge-concave aggregation data */ 2021 SCIP_Real* fvals, /**< array containing all corner values of the aggregation */ 2022 SCIP_Real* facet /**< current facet candidate (of dimension ecaggr->nvars + 1) */ 2023 ) 2024 { 2025 SCIP_Real maxviolation; 2026 SCIP_Real val; 2027 unsigned int i; 2028 unsigned int ncorner; 2029 unsigned int prev; 2030 2031 assert(scip != NULL); 2032 assert(ecaggr != NULL); 2033 assert(fvals != NULL); 2034 assert(facet != NULL); 2035 2036 ncorner = (unsigned int) poweroftwo[ecaggr->nvars]; 2037 maxviolation = 0.0; 2038 2039 /* check for the origin */ 2040 val = facet[ecaggr->nvars]; 2041 for( i = 0; i < (unsigned int) ecaggr->nvars; ++i ) 2042 val += facet[i] * SCIPvarGetLbLocal(ecaggr->vars[i]); 2043 2044 /* update maximum violation */ 2045 maxviolation = MAX(val - fvals[0], maxviolation); 2046 assert(SCIPisFeasEQ(scip, maxviolation, 0.0)); 2047 2048 prev = 0; 2049 for( i = 1; i < ncorner; ++i ) 2050 { 2051 unsigned int gray; 2052 unsigned int diff; 2053 unsigned int pos; 2054 2055 gray = i ^ (i >> 1); 2056 diff = gray ^ prev; 2057 2058 /* compute position of unique 1 of diff */ 2059 pos = 0; 2060 while( (diff >>= 1) != 0 ) 2061 ++pos; 2062 2063 if( gray > prev ) 2064 val += facet[pos] * (SCIPvarGetUbLocal(ecaggr->vars[pos]) - SCIPvarGetLbLocal(ecaggr->vars[pos])); 2065 else 2066 val -= facet[pos] * (SCIPvarGetUbLocal(ecaggr->vars[pos]) - SCIPvarGetLbLocal(ecaggr->vars[pos])); 2067 2068 /* update maximum violation */ 2069 maxviolation = MAX(val - fvals[gray], maxviolation); 2070 assert(SCIPisFeasEQ(scip, maxviolation, 0.0)); 2071 2072 prev = gray; 2073 } 2074 2075 SCIPdebugMsg(scip, "maximum violation of facet: %2.8e\n", maxviolation); 2076 2077 /* there seem to be numerical problems if the violation is too large; in this case we reject the facet */ 2078 if( maxviolation > ADJUSTFACETTOL ) 2079 return FALSE; 2080 2081 /* adjust constant part of the facet */ 2082 facet[ecaggr->nvars] -= maxviolation; 2083 2084 return TRUE; 2085 } 2086 2087 /** set up LP interface to solve LPs to compute the facet of the convex envelope */ 2088 static 2089 SCIP_RETCODE createLP( 2090 SCIP* scip, /**< SCIP data structure */ 2091 SCIP_SEPADATA* sepadata /**< separation data */ 2092 ) 2093 { 2094 SCIP_Real* obj; 2095 SCIP_Real* lb; 2096 SCIP_Real* ub; 2097 SCIP_Real* val; 2098 int* beg; 2099 int* ind; 2100 int nnonz; 2101 int ncols; 2102 int nrows; 2103 int i; 2104 int k; 2105 2106 assert(scip != NULL); 2107 assert(sepadata != NULL); 2108 assert(sepadata->nnlrowaggrs > 0); 2109 2110 /* LP interface has been already created with enough rows/columns*/ 2111 if( sepadata->lpi != NULL && sepadata->lpisize >= sepadata->maxecsize ) 2112 return SCIP_OKAY; 2113 2114 /* size of lpi is too small; reconstruct lpi */ 2115 if( sepadata->lpi != NULL ) 2116 { 2117 SCIP_CALL( SCIPlpiFree(&sepadata->lpi) ); 2118 sepadata->lpi = NULL; 2119 } 2120 2121 assert(sepadata->lpi == NULL); 2122 SCIP_CALL( SCIPlpiCreate(&(sepadata->lpi), SCIPgetMessagehdlr(scip), "e.c. LP", SCIP_OBJSEN_MINIMIZE) ); 2123 sepadata->lpisize = sepadata->maxecsize; 2124 2125 nrows = sepadata->maxecsize + 1; 2126 ncols = poweroftwo[nrows - 1]; 2127 nnonz = (ncols * (nrows + 1)) / 2; 2128 k = 0; 2129 2130 /* allocate necessary memory */ 2131 SCIP_CALL( SCIPallocBufferArray(scip, &obj, ncols) ); 2132 SCIP_CALL( SCIPallocBufferArray(scip, &lb, ncols) ); 2133 SCIP_CALL( SCIPallocBufferArray(scip, &ub, ncols) ); 2134 SCIP_CALL( SCIPallocBufferArray(scip, &beg, ncols) ); 2135 SCIP_CALL( SCIPallocBufferArray(scip, &val, nnonz) ); 2136 SCIP_CALL( SCIPallocBufferArray(scip, &ind, nnonz) ); 2137 2138 /* calculate nonzero entries in the LP; set obj, lb, and ub to zero */ 2139 for( i = 0; i < ncols; ++i ) 2140 { 2141 int row; 2142 int a; 2143 2144 obj[i] = 0.0; 2145 lb[i] = 0.0; 2146 ub[i] = 0.0; 2147 2148 SCIPdebugMsg(scip, "col %i starts at position %d\n", i, k); 2149 beg[i] = k; 2150 row = 0; 2151 a = 1; 2152 2153 /* iterate through the bit representation of i */ 2154 while( a <= i ) 2155 { 2156 if( (a & i) != 0 ) 2157 { 2158 val[k] = 1.0; 2159 ind[k] = row; 2160 2161 SCIPdebugMsg(scip, " val[%d][%d] = 1 (position %d)\n", row, i, k); 2162 2163 ++k; 2164 } 2165 2166 a <<= 1; /*lint !e701*/ 2167 ++row; 2168 assert(poweroftwo[row] == a); 2169 } 2170 2171 /* put 1 as a coefficient for sum_{i} \lambda_i = 1 row (last row) */ 2172 val[k] = 1.0; 2173 ind[k] = nrows - 1; 2174 ++k; 2175 SCIPdebugMsg(scip, " val[%d][%d] = 1 (position %d)\n", nrows - 1, i, k); 2176 } 2177 assert(k == nnonz); 2178 2179 /* 2180 * add all columns to the LP interface 2181 * CPLEX needs the row to exist before adding columns, so we create the rows with dummy sides 2182 * note that the assert is not needed once somebody fixes the LPI 2183 */ 2184 assert(nrows <= ncols); 2185 SCIP_CALL( SCIPlpiAddRows(sepadata->lpi, nrows, obj, obj, NULL, 0, NULL, NULL, NULL) ); 2186 SCIP_CALL( SCIPlpiAddCols(sepadata->lpi, ncols, obj, lb, ub, NULL, nnonz, beg, ind, val) ); 2187 2188 /* free allocated memory */ 2189 SCIPfreeBufferArray(scip, &ind); 2190 SCIPfreeBufferArray(scip, &val); 2191 SCIPfreeBufferArray(scip, &beg); 2192 SCIPfreeBufferArray(scip, &ub); 2193 SCIPfreeBufferArray(scip, &lb); 2194 SCIPfreeBufferArray(scip, &obj); 2195 2196 return SCIP_OKAY; 2197 } 2198 2199 /** evaluates an edge-concave aggregation at a corner of the domain [l,u] */ 2200 static 2201 SCIP_Real evalCorner( 2202 SCIP_ECAGGR* ecaggr, /**< edge-concave aggregation data */ 2203 int k /**< k-th corner */ 2204 ) 2205 { 2206 SCIP_Real val; 2207 int i; 2208 2209 assert(ecaggr != NULL); 2210 assert(k >= 0 && k < poweroftwo[ecaggr->nvars]); 2211 2212 val = 0.0; 2213 2214 for( i = 0; i < ecaggr->nterms; ++i ) 2215 { 2216 SCIP_Real coef; 2217 SCIP_Real bound1; 2218 SCIP_Real bound2; 2219 int idx1; 2220 int idx2; 2221 2222 idx1 = ecaggr->termvars1[i]; 2223 idx2 = ecaggr->termvars2[i]; 2224 coef = ecaggr->termcoefs[i]; 2225 assert(idx1 >= 0 && idx1 < ecaggr->nvars); 2226 assert(idx2 >= 0 && idx2 < ecaggr->nvars); 2227 2228 bound1 = ((poweroftwo[idx1]) & k) == 0 ? SCIPvarGetLbLocal(ecaggr->vars[idx1]) : SCIPvarGetUbLocal(ecaggr->vars[idx1]); /*lint !e661*/ 2229 bound2 = ((poweroftwo[idx2]) & k) == 0 ? SCIPvarGetLbLocal(ecaggr->vars[idx2]) : SCIPvarGetUbLocal(ecaggr->vars[idx2]); /*lint !e661*/ 2230 2231 val += coef * bound1 * bound2; 2232 } 2233 2234 return val; 2235 } 2236 2237 /** returns (val - lb) / (ub - lb) for a in [lb, ub] */ 2238 static 2239 SCIP_Real transformValue( 2240 SCIP* scip, /**< SCIP data structure */ 2241 SCIP_Real lb, /**< lower bound */ 2242 SCIP_Real ub, /**< upper bound */ 2243 SCIP_Real val /**< value in [lb,ub] */ 2244 ) 2245 { 2246 assert(scip != NULL); 2247 assert(!SCIPisInfinity(scip, -lb)); 2248 assert(!SCIPisInfinity(scip, ub)); 2249 assert(!SCIPisInfinity(scip, REALABS(val))); 2250 assert(!SCIPisFeasEQ(scip, ub - lb, 0.0)); /* this would mean that a variable has been fixed */ 2251 2252 /* adjust val */ 2253 val = MIN(val, ub); 2254 val = MAX(val, lb); 2255 2256 val = (val - lb) / (ub - lb); 2257 assert(val >= 0.0 && val <= 1.0); 2258 2259 return val; 2260 } 2261 2262 /** computes a facet of the convex envelope of an edge concave aggregation 2263 * 2264 * The algorithm solves the following LP: 2265 * \f{align}{ 2266 * \min & \sum_i \lambda_i f(v_i)\\ 2267 * s.t. & \sum_i \lambda_i v_i = x\\ 2268 * & \sum_i \lambda_i = 1\\ 2269 * & \lambda \geq 0 2270 * \f} 2271 * where \f$f\f$ is an edge concave function, \f$x\in [l,u]\f$ is a solution of the current relaxation, and \f$v_i\f$ are the vertices of \f$[l,u]\f$. 2272 * The method transforms the problem to the domain \f$[0,1]^n\f$, computes a facet, and transforms this facet to the 2273 * original space. The dual solution of the LP above are the coefficients of the facet. 2274 * 2275 * The complete algorithm works as follows: 2276 * -# compute \f$f(v_i)\f$ for each corner \f$v_i\f$ of \f$[l,u]\f$ 2277 * -# set up the described LP for the transformed space 2278 * -# solve the LP and store the resulting facet for the transformed space 2279 * -# transform the facet to original space 2280 * -# adjust and check facet with the algorithm of Rikun et al. 2281 */ 2282 static 2283 SCIP_RETCODE computeConvexEnvelopeFacet( 2284 SCIP* scip, /**< SCIP data structure */ 2285 SCIP_SEPADATA* sepadata, /**< separation data */ 2286 SCIP_SOL* sol, /**< solution (might be NULL) */ 2287 SCIP_ECAGGR* ecaggr, /**< edge-concave aggregation data */ 2288 SCIP_Real* facet, /**< array to store the coefficients of the resulting facet; size has to be at least (ecaggr->nvars + 1) */ 2289 SCIP_Real* facetval, /**< pointer to store the value of the facet evaluated at the current solution */ 2290 SCIP_Bool* success /**< pointer to store if we have found a facet */ 2291 ) 2292 { 2293 SCIP_Real* fvals; 2294 SCIP_Real* side; 2295 SCIP_Real* lb; 2296 SCIP_Real* ub; 2297 SCIP_Real perturbation; 2298 int* inds; 2299 int ncorner; 2300 int ncols; 2301 int nrows; 2302 int i; 2303 2304 assert(scip != NULL); 2305 assert(sepadata != NULL); 2306 assert(ecaggr != NULL); 2307 assert(facet != NULL); 2308 assert(facetval != NULL); 2309 assert(success != NULL); 2310 assert(ecaggr->nvars <= sepadata->maxecsize); 2311 2312 *facetval = -SCIPinfinity(scip); 2313 *success = FALSE; 2314 2315 /* create LP if this has not been done yet */ 2316 SCIP_CALL( createLP(scip, sepadata) ); 2317 2318 assert(sepadata->lpi != NULL); 2319 assert(sepadata->lpisize >= ecaggr->nvars); 2320 2321 SCIP_CALL( SCIPlpiGetNCols(sepadata->lpi, &ncols) ); 2322 SCIP_CALL( SCIPlpiGetNRows(sepadata->lpi, &nrows) ); 2323 ncorner = poweroftwo[ecaggr->nvars]; 2324 2325 assert(ncorner <= ncols); 2326 assert(ecaggr->nvars + 1 <= nrows); 2327 assert(nrows <= ncols); 2328 2329 /* allocate necessary memory */ 2330 SCIP_CALL( SCIPallocBufferArray(scip, &fvals, ncols) ); 2331 SCIP_CALL( SCIPallocBufferArray(scip, &inds, ncols) ); 2332 SCIP_CALL( SCIPallocBufferArray(scip, &lb, ncols) ); 2333 SCIP_CALL( SCIPallocBufferArray(scip, &ub, ncols) ); 2334 SCIP_CALL( SCIPallocBufferArray(scip, &side, ncols) ); 2335 2336 /* 2337 * 1. compute f(v_i) for each corner v_i of [l,u] 2338 * 2. set up the described LP for the transformed space 2339 */ 2340 for( i = 0; i < ncols; ++i ) 2341 { 2342 fvals[i] = i < ncorner ? evalCorner(ecaggr, i) : 0.0; 2343 inds[i] = i; 2344 2345 /* update bounds; fix variables to zero which are currently not in the LP */ 2346 lb[i] = 0.0; 2347 ub[i] = i < ncorner ? 1.0 : 0.0; 2348 SCIPdebugMsg(scip, "bounds of LP col %d = [%e, %e]; obj = %e\n", i, lb[i], ub[i], fvals[i]); 2349 } 2350 2351 /* update lhs and rhs */ 2352 perturbation = 0.001; 2353 for( i = 0; i < nrows; ++i ) 2354 { 2355 /* note that the last row corresponds to sum_{j} \lambda_j = 1 */ 2356 if( i < ecaggr->nvars ) 2357 { 2358 SCIP_VAR* x; 2359 2360 x = ecaggr->vars[i]; 2361 assert(x != NULL); 2362 2363 side[i] = transformValue(scip, SCIPvarGetLbLocal(x), SCIPvarGetUbLocal(x), SCIPgetSolVal(scip, sol, x)); 2364 2365 /* perturb point to enforce an LP solution with ecaggr->nvars + 1 nonzero */ 2366 side[i] += side[i] > perturbation ? -perturbation : perturbation; 2367 perturbation /= 1.2; 2368 } 2369 else 2370 { 2371 side[i] = (i == nrows - 1) ? 1.0 : 0.0; 2372 } 2373 2374 SCIPdebugMsg(scip, "LP row %d in [%e, %e]\n", i, side[i], side[i]); 2375 } 2376 2377 /* update LP */ 2378 SCIP_CALL( SCIPlpiChgObj(sepadata->lpi, ncols, inds, fvals) ); 2379 SCIP_CALL( SCIPlpiChgBounds(sepadata->lpi, ncols, inds, lb, ub) ); 2380 SCIP_CALL( SCIPlpiChgSides(sepadata->lpi, nrows, inds, side, side) ); 2381 2382 /* free memory used to build the LP */ 2383 SCIPfreeBufferArray(scip, &side); 2384 SCIPfreeBufferArray(scip, &ub); 2385 SCIPfreeBufferArray(scip, &lb); 2386 SCIPfreeBufferArray(scip, &inds); 2387 2388 /* 2389 * 3. solve the LP and store the resulting facet for the transformed space 2390 */ 2391 if( USEDUALSIMPLEX ) /*lint !e774 !e506*/ 2392 { 2393 SCIP_CALL( SCIPlpiSolveDual(sepadata->lpi) ); 2394 } 2395 else 2396 { 2397 SCIP_CALL( SCIPlpiSolvePrimal(sepadata->lpi) ); 2398 } 2399 2400 /* the dual solution corresponds to the coefficients of the facet in the transformed problem; note that it might be 2401 * the case that the dual solution has more components than the facet array 2402 */ 2403 if( ecaggr->nvars + 1 == ncols ) 2404 { 2405 SCIP_CALL( SCIPlpiGetSol(sepadata->lpi, NULL, NULL, facet, NULL, NULL) ); 2406 } 2407 else 2408 { 2409 SCIP_Real* dualsol; 2410 2411 SCIP_CALL( SCIPallocBufferArray(scip, &dualsol, nrows) ); 2412 2413 /* get the dual solution */ 2414 SCIP_CALL( SCIPlpiGetSol(sepadata->lpi, NULL, NULL, dualsol, NULL, NULL) ); 2415 2416 for( i = 0; i < ecaggr->nvars; ++i ) 2417 facet[i] = dualsol[i]; 2418 2419 /* constant part of the facet is the last component of the dual solution */ 2420 facet[ecaggr->nvars] = dualsol[nrows - 1]; 2421 2422 SCIPfreeBufferArray(scip, &dualsol); 2423 } 2424 2425 #ifdef SCIP_DEBUG 2426 SCIPdebugMsg(scip, "facet for the transformed problem: "); 2427 for( i = 0; i < ecaggr->nvars; ++i ) 2428 { 2429 SCIPdebugMsgPrint(scip, "%3.4e * %s + ", facet[i], SCIPvarGetName(ecaggr->vars[i])); 2430 } 2431 SCIPdebugMsgPrint(scip, "%3.4e\n", facet[ecaggr->nvars]); 2432 #endif 2433 2434 /* 2435 * 4. transform the facet to original space 2436 * we now have the linear underestimator L(x) = beta^T x + beta_0, which needs to be transform to the original space 2437 * the underestimator in the original space, G(x) = alpha^T x + alpha_0, is given by G(x) = L(T(x)), where T(.) is 2438 * the transformation applied in step 2; therefore, 2439 * alpha_i = beta_i/(ub_i - lb_i) 2440 * alpha_0 = beta_0 - sum_i lb_i * beta_i/(ub_i - lb_i) 2441 */ 2442 2443 SCIPdebugMsg(scip, "facet in orig. space: "); 2444 *facetval = 0.0; 2445 2446 for( i = 0; i < ecaggr->nvars; ++i ) 2447 { 2448 SCIP_Real varlb; 2449 SCIP_Real varub; 2450 2451 varlb = SCIPvarGetLbLocal(ecaggr->vars[i]); 2452 varub = SCIPvarGetUbLocal(ecaggr->vars[i]); 2453 assert(!SCIPisEQ(scip, varlb, varub)); 2454 2455 /* substract (\beta_i * lb_i) / (ub_i - lb_i) from current alpha_0 */ 2456 facet[ecaggr->nvars] -= (facet[i] * varlb) / (varub - varlb); 2457 2458 /* set \alpha_i := \beta_i / (ub_i - lb_i) */ 2459 facet[i] = facet[i] / (varub - varlb); 2460 *facetval += facet[i] * SCIPgetSolVal(scip, sol, ecaggr->vars[i]); 2461 2462 SCIPdebugMsgPrint(scip, "%3.4e * %s + ", facet[i], SCIPvarGetName(ecaggr->vars[i])); 2463 } 2464 2465 /* add constant part to the facet value */ 2466 *facetval += facet[ecaggr->nvars]; 2467 SCIPdebugMsgPrint(scip, "%3.4e\n", facet[ecaggr->nvars]); 2468 2469 /* 2470 * 5. adjust and check facet with the algorithm of Rikun et al. 2471 */ 2472 2473 if( checkRikun(scip, ecaggr, fvals, facet) ) 2474 { 2475 SCIPdebugMsg(scip, "facet pass the check of Rikun et al.\n"); 2476 *success = TRUE; 2477 } 2478 2479 /* free allocated memory */ 2480 SCIPfreeBufferArray(scip, &fvals); 2481 2482 return SCIP_OKAY; 2483 } 2484 2485 /* 2486 * miscellaneous methods 2487 */ 2488 2489 /** method to add a facet of the convex envelope of an edge-concave aggregation to a given cut */ 2490 static 2491 SCIP_RETCODE addFacetToCut( 2492 SCIP* scip, /**< SCIP data structure */ 2493 SCIP_SOL* sol, /**< current solution (might be NULL) */ 2494 SCIP_ROW* cut, /**< current cut (modifiable) */ 2495 SCIP_Real* facet, /**< coefficient of the facet (dimension nvars + 1) */ 2496 SCIP_VAR** vars, /**< variables of the facet */ 2497 int nvars, /**< number of variables in the facet */ 2498 SCIP_Real* cutconstant, /**< pointer to update the constant part of the facet */ 2499 SCIP_Real* cutactivity, /**< pointer to update the activity of the cut */ 2500 SCIP_Bool* success /**< pointer to store if everything went fine */ 2501 ) 2502 { 2503 int i; 2504 2505 assert(cut != NULL); 2506 assert(facet != NULL); 2507 assert(vars != NULL); 2508 assert(nvars > 0); 2509 assert(cutconstant != NULL); 2510 assert(cutactivity != NULL); 2511 assert(success != NULL); 2512 2513 *success = TRUE; 2514 2515 for( i = 0; i < nvars; ++i ) 2516 { 2517 if( SCIPisInfinity(scip, REALABS(facet[i])) ) 2518 { 2519 *success = FALSE; 2520 return SCIP_OKAY; 2521 } 2522 2523 if( !SCIPisZero(scip, facet[i]) ) 2524 { 2525 /* add only a constant if the variable has been fixed */ 2526 if( SCIPvarGetLbLocal(vars[i]) == SCIPvarGetUbLocal(vars[i]) ) /*lint !e777*/ 2527 { 2528 assert(SCIPisFeasEQ(scip, SCIPvarGetLbLocal(vars[i]), SCIPgetSolVal(scip, sol, vars[i]))); 2529 *cutconstant += facet[i] * SCIPgetSolVal(scip, sol, vars[i]); 2530 *cutactivity += facet[i] * SCIPgetSolVal(scip, sol, vars[i]); 2531 } 2532 else 2533 { 2534 *cutactivity += facet[i] * SCIPgetSolVal(scip, sol, vars[i]); 2535 SCIP_CALL( SCIPaddVarToRow(scip, cut, vars[i], facet[i]) ); 2536 } 2537 } 2538 } 2539 2540 /* add constant part of the facet */ 2541 *cutconstant += facet[nvars]; 2542 *cutactivity += facet[nvars]; 2543 2544 return SCIP_OKAY; 2545 } 2546 2547 /** method to add a linear term to a given cut */ 2548 static 2549 SCIP_RETCODE addLinearTermToCut( 2550 SCIP* scip, /**< SCIP data structure */ 2551 SCIP_SOL* sol, /**< current solution (might be NULL) */ 2552 SCIP_ROW* cut, /**< current cut (modifiable) */ 2553 SCIP_VAR* x, /**< linear variable */ 2554 SCIP_Real coeff, /**< coefficient */ 2555 SCIP_Real* cutconstant, /**< pointer to update the constant part of the facet */ 2556 SCIP_Real* cutactivity, /**< pointer to update the activity of the cut */ 2557 SCIP_Bool* success /**< pointer to store if everything went fine */ 2558 ) 2559 { 2560 SCIP_Real activity; 2561 2562 assert(cut != NULL); 2563 assert(x != NULL); 2564 assert(!SCIPisZero(scip, coeff)); 2565 assert(!SCIPisInfinity(scip, coeff)); 2566 assert(cutconstant != NULL); 2567 assert(cutactivity != NULL); 2568 assert(success != NULL); 2569 2570 *success = TRUE; 2571 activity = SCIPgetSolVal(scip, sol, x) * coeff; 2572 2573 /* do not add a term if the activity is -infinity */ 2574 if( SCIPisInfinity(scip, -1.0 * REALABS(activity)) ) 2575 { 2576 *success = FALSE; 2577 return SCIP_OKAY; 2578 } 2579 2580 /* add activity to the constant part if the variable has been fixed */ 2581 if( SCIPvarGetLbLocal(x) == SCIPvarGetUbLocal(x) ) /*lint !e777*/ 2582 { 2583 assert(SCIPisFeasEQ(scip, SCIPvarGetLbLocal(x), SCIPgetSolVal(scip, sol, x))); 2584 *cutconstant += activity; 2585 SCIPdebugMsg(scip, "add to cut: %e\n", activity); 2586 } 2587 else 2588 { 2589 SCIP_CALL( SCIPaddVarToRow(scip, cut, x, coeff) ); 2590 SCIPdebugMsg(scip, "add to cut: %e * %s\n", coeff, SCIPvarGetName(x)); 2591 } 2592 2593 *cutactivity += activity; 2594 2595 return SCIP_OKAY; 2596 } 2597 2598 /** method to add an underestimate of a bilinear term to a given cut */ 2599 static 2600 SCIP_RETCODE addBilinearTermToCut( 2601 SCIP* scip, /**< SCIP data structure */ 2602 SCIP_SOL* sol, /**< current solution (might be NULL) */ 2603 SCIP_ROW* cut, /**< current cut (modifiable) */ 2604 SCIP_VAR* x, /**< first bilinear variable */ 2605 SCIP_VAR* y, /**< seconds bilinear variable */ 2606 SCIP_Real coeff, /**< coefficient */ 2607 SCIP_Real* cutconstant, /**< pointer to update the constant part of the facet */ 2608 SCIP_Real* cutactivity, /**< pointer to update the activity of the cut */ 2609 SCIP_Bool* success /**< pointer to store if everything went fine */ 2610 ) 2611 { 2612 SCIP_Real activity; 2613 2614 assert(cut != NULL); 2615 assert(x != NULL); 2616 assert(y != NULL); 2617 assert(!SCIPisZero(scip, coeff)); 2618 assert(cutconstant != NULL); 2619 assert(cutactivity != NULL); 2620 assert(success != NULL); 2621 2622 *success = TRUE; 2623 activity = coeff * SCIPgetSolVal(scip, sol, x) * SCIPgetSolVal(scip, sol, y); 2624 2625 if( SCIPisInfinity(scip, REALABS(coeff)) ) 2626 { 2627 *success = FALSE; 2628 return SCIP_OKAY; 2629 } 2630 2631 /* do not add a term if the activity is -infinity */ 2632 if( SCIPisInfinity(scip, -1.0 * REALABS(activity)) ) 2633 { 2634 *success = FALSE; 2635 return SCIP_OKAY; 2636 } 2637 2638 /* quadratic case */ 2639 if( x == y ) 2640 { 2641 SCIP_Real refpoint; 2642 SCIP_Real lincoef; 2643 SCIP_Real linconst; 2644 2645 lincoef = 0.0; 2646 linconst = 0.0; 2647 refpoint = SCIPgetSolVal(scip, sol, x); 2648 2649 /* adjust the reference point */ 2650 refpoint = SCIPisLT(scip, refpoint, SCIPvarGetLbLocal(x)) ? SCIPvarGetLbLocal(x) : refpoint; 2651 refpoint = SCIPisGT(scip, refpoint, SCIPvarGetUbLocal(x)) ? SCIPvarGetUbLocal(x) : refpoint; 2652 assert(SCIPisLE(scip, refpoint, SCIPvarGetUbLocal(x)) && SCIPisGE(scip, refpoint, SCIPvarGetLbLocal(x))); 2653 2654 if( SCIPisPositive(scip, coeff) ) 2655 SCIPaddSquareLinearization(scip, coeff, refpoint, SCIPvarIsIntegral(x), &lincoef, &linconst, success); 2656 else 2657 SCIPaddSquareSecant(scip, coeff, SCIPvarGetLbLocal(x), SCIPvarGetUbLocal(x), &lincoef, &linconst, success); 2658 2659 *cutactivity += lincoef * refpoint + linconst; 2660 *cutconstant += linconst; 2661 2662 /* add underestimate to cut */ 2663 SCIP_CALL( SCIPaddVarToRow(scip, cut, x, lincoef) ); 2664 2665 SCIPdebugMsg(scip, "add to cut: %e * %s + %e\n", lincoef, SCIPvarGetName(x), linconst); 2666 } 2667 /* bilinear case */ 2668 else 2669 { 2670 SCIP_Real refpointx; 2671 SCIP_Real refpointy; 2672 SCIP_Real lincoefx; 2673 SCIP_Real lincoefy; 2674 SCIP_Real linconst; 2675 2676 lincoefx = 0.0; 2677 lincoefy = 0.0; 2678 linconst = 0.0; 2679 refpointx = SCIPgetSolVal(scip, sol, x); 2680 refpointy = SCIPgetSolVal(scip, sol, y); 2681 2682 /* adjust the reference points */ 2683 refpointx = SCIPisLT(scip, refpointx, SCIPvarGetLbLocal(x)) ? SCIPvarGetLbLocal(x) : refpointx; 2684 refpointx = SCIPisGT(scip, refpointx, SCIPvarGetUbLocal(x)) ? SCIPvarGetUbLocal(x) : refpointx; 2685 refpointy = SCIPisLT(scip, refpointy, SCIPvarGetLbLocal(y)) ? SCIPvarGetLbLocal(y) : refpointy; 2686 refpointy = SCIPisGT(scip, refpointy, SCIPvarGetUbLocal(y)) ? SCIPvarGetUbLocal(y) : refpointy; 2687 assert(SCIPisLE(scip, refpointx, SCIPvarGetUbLocal(x)) && SCIPisGE(scip, refpointx, SCIPvarGetLbLocal(x))); 2688 assert(SCIPisLE(scip, refpointy, SCIPvarGetUbLocal(y)) && SCIPisGE(scip, refpointy, SCIPvarGetLbLocal(y))); 2689 2690 SCIPaddBilinMcCormick(scip, coeff, SCIPvarGetLbLocal(x), SCIPvarGetUbLocal(x), refpointx, SCIPvarGetLbLocal(y), 2691 SCIPvarGetUbLocal(y), refpointy, FALSE, &lincoefx, &lincoefy, &linconst, success); 2692 2693 *cutactivity += lincoefx * refpointx + lincoefy * refpointy + linconst; 2694 *cutconstant += linconst; 2695 2696 /* add underestimate to cut */ 2697 SCIP_CALL( SCIPaddVarToRow(scip, cut, x, lincoefx) ); 2698 SCIP_CALL( SCIPaddVarToRow(scip, cut, y, lincoefy) ); 2699 2700 SCIPdebugMsg(scip, "add to cut: %e * %s + %e * %s + %e\n", lincoefx, SCIPvarGetName(x), lincoefy, 2701 SCIPvarGetName(y), linconst); 2702 } 2703 2704 return SCIP_OKAY; 2705 } 2706 2707 /** method to compute and add a cut for a nonlinear row aggregation and a given solution 2708 * 2709 * we compute for each edge concave aggregation one facet; 2710 * the remaining bilinear terms will be underestimated with McCormick, secants or linearizations; 2711 * constant and linear terms will be added to the cut directly 2712 */ 2713 static 2714 SCIP_RETCODE computeCut( 2715 SCIP* scip, /**< SCIP data structure */ 2716 SCIP_SEPA* sepa, /**< separator */ 2717 SCIP_SEPADATA* sepadata, /**< separator data */ 2718 SCIP_NLROWAGGR* nlrowaggr, /**< nonlinear row aggregation */ 2719 SCIP_SOL* sol, /**< current solution (might be NULL) */ 2720 SCIP_Bool* separated, /**< pointer to store if we could separate the current solution */ 2721 SCIP_Bool* cutoff /**< pointer to store if the current node gets cut off */ 2722 ) 2723 { 2724 SCIP_ROW* cut; 2725 SCIP_Real* bestfacet; 2726 SCIP_Real bestfacetval; 2727 SCIP_Real cutconstant; 2728 SCIP_Real cutactivity; 2729 int bestfacetsize; 2730 char cutname[SCIP_MAXSTRLEN]; 2731 SCIP_Bool found; 2732 SCIP_Bool islocalcut; 2733 int i; 2734 2735 assert(separated != NULL); 2736 assert(cutoff != NULL); 2737 assert(nlrowaggr->necaggr > 0); 2738 assert(nlrowaggr->nlrow != NULL); 2739 assert(SCIPnlrowIsInNLP(nlrowaggr->nlrow)); 2740 2741 *separated = FALSE; 2742 *cutoff = FALSE; 2743 /* we use SCIPgetDepth because we add the cut to the global cut pool if cut is globally valid */ 2744 islocalcut = SCIPgetDepth(scip) != 0; 2745 2746 /* create the cut */ 2747 (void) SCIPsnprintf(cutname, SCIP_MAXSTRLEN, "ec"); 2748 SCIP_CALL( SCIPcreateEmptyRowSepa(scip, &cut, sepa, cutname, -SCIPinfinity(scip), SCIPinfinity(scip), islocalcut, FALSE, 2749 sepadata->dynamiccuts) ); 2750 SCIP_CALL( SCIPcacheRowExtensions(scip, cut) ); 2751 2752 /* track rhs and activity of the cut */ 2753 cutconstant = nlrowaggr->constant; 2754 cutactivity = 0.0; 2755 2756 /* allocate necessary memory */ 2757 bestfacetsize = sepadata->maxaggrsize + 1; 2758 SCIP_CALL( SCIPallocBufferArray(scip, &bestfacet, bestfacetsize) ); 2759 2760 #ifdef SCIP_DEBUG 2761 SCIP_CALL( SCIPprintNlRow(scip, nlrowaggr->nlrow, NULL) ); 2762 2763 SCIPdebugMsg(scip, "current solution:\n"); 2764 for( i = 0; i < SCIPgetNVars(scip); ++i ) 2765 { 2766 SCIP_VAR* var = SCIPgetVars(scip)[i]; 2767 SCIPdebugMsg(scip, " %s = [%e, %e] solval = %e\n", SCIPvarGetName(var), SCIPvarGetLbLocal(var), 2768 SCIPvarGetUbLocal(var), SCIPgetSolVal(scip, sol, var)); 2769 } 2770 #endif 2771 2772 /* compute a facet for each edge-concave aggregation */ 2773 for( i = 0; i < nlrowaggr->necaggr; ++i ) 2774 { 2775 SCIP_ECAGGR* ecaggr; 2776 SCIP_Bool success; 2777 2778 ecaggr = nlrowaggr->ecaggr[i]; 2779 assert(ecaggr != NULL); 2780 2781 /* compute a facet of the convex envelope */ 2782 SCIP_CALL( computeConvexEnvelopeFacet(scip, sepadata, sol, ecaggr, bestfacet, &bestfacetval, &found) ); 2783 2784 SCIPdebugMsg(scip, "found facet for edge-concave aggregation %d/%d ? %s\n", i, nlrowaggr->necaggr, 2785 found ? "yes" : "no"); 2786 2787 #ifdef SCIP_DEBUG 2788 if( found ) 2789 printFacet(scip, ecaggr->vars, ecaggr->nvars, bestfacet, bestfacetval); 2790 #endif 2791 2792 /* do not add any cut because we did not found a facet for at least one edge-concave aggregation */ 2793 if( !found ) /*lint !e774*/ 2794 goto TERMINATE; 2795 2796 /* add facet to the cut and update the rhs and activity of the cut */ 2797 SCIP_CALL( addFacetToCut(scip, sol, cut, bestfacet, ecaggr->vars, ecaggr->nvars, &cutconstant, &cutactivity, 2798 &success) ); 2799 2800 if( !success ) 2801 goto TERMINATE; 2802 } 2803 2804 /* compute an underestimate for each bilinear term which is not in any edge-concave aggregation */ 2805 for( i = 0; i < nlrowaggr->nremterms; ++i ) 2806 { 2807 SCIP_VAR* x; 2808 SCIP_VAR* y; 2809 SCIP_Bool success; 2810 2811 x = nlrowaggr->remtermvars1[i]; 2812 y = nlrowaggr->remtermvars2[i]; 2813 assert(x != NULL); 2814 assert(y != NULL); 2815 2816 SCIP_CALL( addBilinearTermToCut(scip, sol, cut, x, y, nlrowaggr->remtermcoefs[i], &cutconstant, &cutactivity, 2817 &success) ); 2818 2819 if( !success ) 2820 goto TERMINATE; 2821 } 2822 2823 /* add all linear terms to the cut */ 2824 for( i = 0; i < nlrowaggr->nlinvars; ++i ) 2825 { 2826 SCIP_VAR* x; 2827 SCIP_Real coef; 2828 SCIP_Bool success; 2829 2830 x = nlrowaggr->linvars[i]; 2831 assert(x != NULL); 2832 2833 coef = nlrowaggr->lincoefs[i]; 2834 2835 SCIP_CALL( addLinearTermToCut(scip, sol, cut, x, coef, &cutconstant, &cutactivity, &success) ); 2836 2837 if( !success ) 2838 goto TERMINATE; 2839 } 2840 2841 SCIPdebugMsg(scip, "cut activity = %e rhs(nlrow) = %e\n", cutactivity, nlrowaggr->rhs); 2842 2843 /* set rhs of the cut (substract the constant part of the cut) */ 2844 SCIP_CALL( SCIPchgRowRhs(scip, cut, nlrowaggr->rhs - cutconstant) ); 2845 SCIP_CALL( SCIPflushRowExtensions(scip, cut) ); 2846 2847 /* check activity of the row; this assert can fail because of numerics */ 2848 /* assert(SCIPisFeasEQ(scip, cutactivity - cutconstant, SCIPgetRowSolActivity(scip, cut, sol)) ); */ 2849 2850 #ifdef SCIP_DEBUG 2851 SCIP_CALL( SCIPprintRow(scip, cut, NULL) ); 2852 #endif 2853 2854 SCIPdebugMsg(scip, "EC cut <%s>: act=%f eff=%f rank=%d range=%e\n", 2855 SCIProwGetName(cut), SCIPgetRowSolActivity(scip, cut, sol), SCIPgetCutEfficacy(scip, sol, cut), 2856 SCIProwGetRank(cut), SCIPgetRowMaxCoef(scip, cut) / SCIPgetRowMinCoef(scip, cut) ); 2857 2858 /* try to add the cut has a finite rhs, is efficacious, and does not exceed the maximum cut range */ 2859 if( !SCIPisInfinity(scip, nlrowaggr->rhs - cutconstant) && SCIPisCutEfficacious(scip, sol, cut) 2860 && SCIPgetRowMaxCoef(scip, cut) / SCIPgetRowMinCoef(scip, cut) < sepadata->cutmaxrange ) 2861 { 2862 /* add the cut if it is separating the given solution by at least minviolation */ 2863 if( SCIPisGE(scip, cutactivity - nlrowaggr->rhs, sepadata->minviolation) ) 2864 { 2865 SCIP_CALL( SCIPaddRow(scip, cut, FALSE, cutoff) ); 2866 *separated = TRUE; 2867 SCIPdebugMsg(scip, "added separating cut\n"); 2868 } 2869 2870 if( !(*cutoff) && !islocalcut ) 2871 { 2872 SCIP_CALL( SCIPaddPoolCut(scip, cut) ); 2873 SCIPdebugMsg(scip, "added cut to cut pool\n"); 2874 } 2875 } 2876 2877 TERMINATE: 2878 /* free allocated memory */ 2879 SCIPfreeBufferArray(scip, &bestfacet); 2880 2881 /* release the row */ 2882 SCIP_CALL( SCIPreleaseRow(scip, &cut) ); 2883 2884 return SCIP_OKAY; 2885 } 2886 2887 /** returns whether it is possible to compute a cut for a given nonlinear row aggregation */ 2888 static 2889 SCIP_Bool isPossibleToComputeCut( 2890 SCIP* scip, /**< SCIP data structure */ 2891 SCIP_SOL* sol, /**< current solution (might be NULL) */ 2892 SCIP_NLROWAGGR* nlrowaggr /**< nonlinear row aggregation */ 2893 ) 2894 { 2895 int i; 2896 2897 assert(scip != NULL); 2898 assert(nlrowaggr != NULL); 2899 2900 if( !SCIPnlrowIsInNLP(nlrowaggr->nlrow) ) 2901 { 2902 SCIPdebugMsg(scip, "nlrow is not in NLP anymore\n"); 2903 return FALSE; 2904 } 2905 2906 for( i = 0; i < nlrowaggr->nquadvars; ++i ) 2907 { 2908 SCIP_VAR* var = nlrowaggr->quadvars[i]; 2909 assert(var != NULL); 2910 2911 /* check whether the variable has infinite bounds */ 2912 if( SCIPisInfinity(scip, REALABS(SCIPvarGetLbLocal(var))) || SCIPisInfinity(scip, REALABS(SCIPvarGetUbLocal(var))) 2913 || SCIPisInfinity(scip, REALABS(SCIPgetSolVal(scip, sol, var))) ) 2914 { 2915 SCIPdebugMsg(scip, "nlrow aggregation contains unbounded variables\n"); 2916 return FALSE; 2917 } 2918 2919 /* check whether the variable has been fixed and is in one edge-concave aggregation */ 2920 if( nlrowaggr->quadvar2aggr[i] >= 0 && SCIPisFeasEQ(scip, SCIPvarGetLbLocal(var), SCIPvarGetUbLocal(var)) ) 2921 { 2922 SCIPdebugMsg(scip, "nlrow aggregation contains fixed variables in an e.c. aggregation\n"); 2923 return FALSE; 2924 } 2925 } 2926 2927 return TRUE; 2928 } 2929 2930 /** searches and tries to add edge-concave cuts */ 2931 static 2932 SCIP_RETCODE separateCuts( 2933 SCIP* scip, /**< SCIP data structure */ 2934 SCIP_SEPA* sepa, /**< separator */ 2935 SCIP_SEPADATA* sepadata, /**< separator data */ 2936 int depth, /**< current depth */ 2937 SCIP_SOL* sol, /**< current solution */ 2938 SCIP_RESULT* result /**< pointer to store the result of the separation call */ 2939 ) 2940 { 2941 int nmaxcuts; 2942 int ncuts; 2943 int i; 2944 2945 assert(*result == SCIP_DIDNOTRUN); 2946 2947 SCIPdebugMsg(scip, "separate cuts...\n"); 2948 2949 /* skip if there are no nonlinear row aggregations */ 2950 if( sepadata->nnlrowaggrs == 0 ) 2951 { 2952 SCIPdebugMsg(scip, "no aggregations exists -> skip call\n"); 2953 return SCIP_OKAY; 2954 } 2955 2956 /* get the maximal number of cuts allowed in a separation round */ 2957 nmaxcuts = depth == 0 ? sepadata->maxsepacutsroot : sepadata->maxsepacuts; 2958 ncuts = 0; 2959 2960 /* try to compute cuts for each nonlinear row independently */ 2961 for( i = 0; i < sepadata->nnlrowaggrs && ncuts < nmaxcuts && !SCIPisStopped(scip); ++i ) 2962 { 2963 SCIP_NLROWAGGR* nlrowaggr; 2964 SCIP_Bool separated; 2965 SCIP_Bool cutoff; 2966 2967 nlrowaggr = sepadata->nlrowaggrs[i]; 2968 assert(nlrowaggr != NULL); 2969 2970 /* skip nonlinear aggregations for which it is obviously not possible to compute a cut */ 2971 if( !isPossibleToComputeCut(scip, sol, nlrowaggr) ) 2972 return SCIP_OKAY; 2973 2974 *result = (*result == SCIP_DIDNOTRUN) ? SCIP_DIDNOTFIND : *result; 2975 2976 SCIPdebugMsg(scip, "try to compute a cut for nonlinear row aggregation %d\n", i); 2977 2978 /* compute and add cut */ 2979 SCIP_CALL( computeCut(scip, sepa, sepadata, nlrowaggr, sol, &separated, &cutoff) ); 2980 SCIPdebugMsg(scip, "found a cut: %s cutoff: %s\n", separated ? "yes" : "no", cutoff ? "yes" : "no"); 2981 2982 /* stop if the current node gets cut off */ 2983 if( cutoff ) 2984 { 2985 assert(separated); 2986 *result = SCIP_CUTOFF; 2987 return SCIP_OKAY; 2988 } 2989 2990 /* do not compute more cuts if we already separated the given solution */ 2991 if( separated ) 2992 { 2993 assert(!cutoff); 2994 *result = SCIP_SEPARATED; 2995 ++ncuts; 2996 } 2997 } 2998 2999 return SCIP_OKAY; 3000 } 3001 3002 /* 3003 * Callback methods of separator 3004 */ 3005 3006 /** copy method for separator plugins (called when SCIP copies plugins) */ 3007 static 3008 SCIP_DECL_SEPACOPY(sepaCopyEccuts) 3009 { /*lint --e{715}*/ 3010 assert(scip != NULL); 3011 assert(sepa != NULL); 3012 assert(strcmp(SCIPsepaGetName(sepa), SEPA_NAME) == 0); 3013 3014 /* call inclusion method of constraint handler */ 3015 SCIP_CALL( SCIPincludeSepaEccuts(scip) ); 3016 3017 return SCIP_OKAY; 3018 } 3019 3020 /** destructor of separator to free user data (called when SCIP is exiting) */ 3021 static 3022 SCIP_DECL_SEPAFREE(sepaFreeEccuts) 3023 { /*lint --e{715}*/ 3024 SCIP_SEPADATA* sepadata; 3025 3026 sepadata = SCIPsepaGetData(sepa); 3027 assert(sepadata != NULL); 3028 3029 SCIP_CALL( sepadataFree(scip, &sepadata) ); 3030 SCIPsepaSetData(sepa, NULL); 3031 3032 return SCIP_OKAY; 3033 } 3034 3035 /** solving process deinitialization method of separator (called before branch and bound process data is freed) */ 3036 static 3037 SCIP_DECL_SEPAEXITSOL(sepaExitsolEccuts) 3038 { /*lint --e{715}*/ 3039 SCIP_SEPADATA* sepadata; 3040 3041 sepadata = SCIPsepaGetData(sepa); 3042 assert(sepadata != NULL); 3043 3044 /* print statistics */ 3045 #ifdef SCIP_STATISTIC 3046 SCIPstatisticMessage("rhs-AGGR %d\n", sepadata->nrhsnlrowaggrs); 3047 SCIPstatisticMessage("lhs-AGGR %d\n", sepadata->nlhsnlrowaggrs); 3048 SCIPstatisticMessage("aggr. search time = %f\n", sepadata->aggrsearchtime); 3049 #endif 3050 3051 /* free nonlinear row aggregations */ 3052 SCIP_CALL( sepadataFreeNlrows(scip, sepadata) ); 3053 3054 /* mark that we should search again for nonlinear row aggregations */ 3055 sepadata->searchedforaggr = FALSE; 3056 3057 SCIPdebugMsg(scip, "exitsol\n"); 3058 3059 return SCIP_OKAY; 3060 } 3061 3062 /** LP solution separation method of separator */ 3063 static 3064 SCIP_DECL_SEPAEXECLP(sepaExeclpEccuts) 3065 { /*lint --e{715}*/ 3066 SCIP_SEPADATA* sepadata; 3067 int ncalls; 3068 3069 sepadata = SCIPsepaGetData(sepa); 3070 assert(sepadata != NULL); 3071 3072 *result = SCIP_DIDNOTRUN; 3073 3074 if( !allowlocal ) 3075 return SCIP_OKAY; 3076 3077 /* check min- and maximal aggregation size */ 3078 if( sepadata->maxaggrsize < sepadata->minaggrsize ) 3079 return SCIP_PARAMETERWRONGVAL; 3080 3081 /* only call separator, if we are not close to terminating */ 3082 if( SCIPisStopped(scip) ) 3083 return SCIP_OKAY; 3084 3085 /* skip if the LP is not constructed yet */ 3086 if( !SCIPisNLPConstructed(scip) ) 3087 { 3088 SCIPdebugMsg(scip, "Skip since NLP is not constructed yet.\n"); 3089 return SCIP_OKAY; 3090 } 3091 3092 /* only call separator up to a maximum depth */ 3093 if ( sepadata->maxdepth >= 0 && depth > sepadata->maxdepth ) 3094 return SCIP_OKAY; 3095 3096 /* only call separator a given number of times at each node */ 3097 ncalls = SCIPsepaGetNCallsAtNode(sepa); 3098 if ( (depth == 0 && sepadata->maxroundsroot >= 0 && ncalls >= sepadata->maxroundsroot) 3099 || (depth > 0 && sepadata->maxrounds >= 0 && ncalls >= sepadata->maxrounds) ) 3100 return SCIP_OKAY; 3101 3102 /* search for nonlinear row aggregations */ 3103 if( !sepadata->searchedforaggr ) 3104 { 3105 int i; 3106 3107 SCIPstatistic( sepadata->aggrsearchtime -= SCIPgetTotalTime(scip) ); 3108 3109 SCIPdebugMsg(scip, "search for nonlinear row aggregations\n"); 3110 for( i = 0; i < SCIPgetNNLPNlRows(scip) && !SCIPisStopped(scip); ++i ) 3111 { 3112 SCIP_NLROW* nlrow = SCIPgetNLPNlRows(scip)[i]; 3113 SCIP_CALL( findAndStoreEcAggregations(scip, sepadata, nlrow, NULL) ); 3114 } 3115 sepadata->searchedforaggr = TRUE; 3116 3117 SCIPstatistic( sepadata->aggrsearchtime += SCIPgetTotalTime(scip) ); 3118 } 3119 3120 /* search for edge-concave cuts */ 3121 SCIP_CALL( separateCuts(scip, sepa, sepadata, depth, NULL, result) ); 3122 3123 return SCIP_OKAY; 3124 } 3125 3126 /* 3127 * separator specific interface methods 3128 */ 3129 3130 /** creates the edge-concave separator and includes it in SCIP 3131 * 3132 * @ingroup SeparatorIncludes 3133 */ 3134 SCIP_RETCODE SCIPincludeSepaEccuts( 3135 SCIP* scip /**< SCIP data structure */ 3136 ) 3137 { 3138 SCIP_SEPADATA* sepadata; 3139 SCIP_SEPA* sepa; 3140 3141 /* create eccuts separator data */ 3142 SCIP_CALL( sepadataCreate(scip, &sepadata) ); 3143 3144 /* include separator */ 3145 SCIP_CALL( SCIPincludeSepaBasic(scip, &sepa, SEPA_NAME, SEPA_DESC, SEPA_PRIORITY, SEPA_FREQ, SEPA_MAXBOUNDDIST, 3146 SEPA_USESSUBSCIP, SEPA_DELAY, sepaExeclpEccuts, NULL, sepadata) ); 3147 3148 assert(sepa != NULL); 3149 3150 /* set non fundamental callbacks via setter functions */ 3151 SCIP_CALL( SCIPsetSepaCopy(scip, sepa, sepaCopyEccuts) ); 3152 SCIP_CALL( SCIPsetSepaFree(scip, sepa, sepaFreeEccuts) ); 3153 SCIP_CALL( SCIPsetSepaExitsol(scip, sepa, sepaExitsolEccuts) ); 3154 3155 /* add eccuts separator parameters */ 3156 SCIP_CALL( SCIPaddBoolParam(scip, 3157 "separating/" SEPA_NAME "/dynamiccuts", 3158 "should generated cuts be removed from the LP if they are no longer tight?", 3159 &sepadata->dynamiccuts, FALSE, DEFAULT_DYNAMICCUTS, NULL, NULL) ); 3160 3161 SCIP_CALL( SCIPaddIntParam(scip, 3162 "separating/" SEPA_NAME "/maxrounds", 3163 "maximal number of eccuts separation rounds per node (-1: unlimited)", 3164 &sepadata->maxrounds, FALSE, DEFAULT_MAXROUNDS, -1, INT_MAX, NULL, NULL) ); 3165 3166 SCIP_CALL( SCIPaddIntParam(scip, 3167 "separating/" SEPA_NAME "/maxroundsroot", 3168 "maximal number of eccuts separation rounds in the root node (-1: unlimited)", 3169 &sepadata->maxroundsroot, FALSE, DEFAULT_MAXROUNDSROOT, -1, INT_MAX, NULL, NULL) ); 3170 3171 SCIP_CALL( SCIPaddIntParam(scip, 3172 "separating/" SEPA_NAME "/maxdepth", 3173 "maximal depth at which the separator is applied (-1: unlimited)", 3174 &sepadata->maxdepth, FALSE, DEFAULT_MAXDEPTH, -1, INT_MAX, NULL, NULL) ); 3175 3176 SCIP_CALL( SCIPaddIntParam(scip, 3177 "separating/" SEPA_NAME "/maxsepacuts", 3178 "maximal number of edge-concave cuts separated per separation round", 3179 &sepadata->maxsepacuts, FALSE, DEFAULT_MAXSEPACUTS, 0, INT_MAX, NULL, NULL) ); 3180 3181 SCIP_CALL( SCIPaddIntParam(scip, 3182 "separating/" SEPA_NAME "/maxsepacutsroot", 3183 "maximal number of edge-concave cuts separated per separation round in the root node", 3184 &sepadata->maxsepacutsroot, FALSE, DEFAULT_MAXSEPACUTSROOT, 0, INT_MAX, NULL, NULL) ); 3185 3186 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/cutmaxrange", 3187 "maximal coef. range of a cut (max coef. divided by min coef.) in order to be added to LP relaxation", 3188 &sepadata->cutmaxrange, FALSE, DEFAULT_CUTMAXRANGE, 0.0, SCIPinfinity(scip), NULL, NULL) ); 3189 3190 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/minviolation", 3191 "minimal violation of an edge-concave cut to be separated", 3192 &sepadata->minviolation, FALSE, DEFAULT_MINVIOLATION, 0.0, 0.5, NULL, NULL) ); 3193 3194 SCIP_CALL( SCIPaddIntParam(scip, 3195 "separating/" SEPA_NAME "/minaggrsize", 3196 "search for edge-concave aggregations of at least this size", 3197 &sepadata->minaggrsize, TRUE, DEFAULT_MINAGGRSIZE, 3, 5, NULL, NULL) ); 3198 3199 SCIP_CALL( SCIPaddIntParam(scip, 3200 "separating/" SEPA_NAME "/maxaggrsize", 3201 "search for edge-concave aggregations of at most this size", 3202 &sepadata->maxaggrsize, TRUE, DEFAULT_MAXAGGRSIZE, 3, 5, NULL, NULL) ); 3203 3204 SCIP_CALL( SCIPaddIntParam(scip, 3205 "separating/" SEPA_NAME "/maxbilinterms", 3206 "maximum number of bilinear terms allowed to be in a quadratic constraint", 3207 &sepadata->maxbilinterms, TRUE, DEFAULT_MAXBILINTERMS, 0, INT_MAX, NULL, NULL) ); 3208 3209 SCIP_CALL( SCIPaddIntParam(scip, 3210 "separating/" SEPA_NAME "/maxstallrounds", 3211 "maximum number of unsuccessful rounds in the edge-concave aggregation search", 3212 &sepadata->maxstallrounds, TRUE, DEFAULT_MAXSTALLROUNDS, 0, INT_MAX, NULL, NULL) ); 3213 3214 return SCIP_OKAY; 3215 } 3216