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 treemodel.c 26 * @brief Branching rules based on the Single-Variable-Branching (SVB) model 27 * @author Daniel Anderson 28 * @author Pierre Le Bodic 29 * 30 * The Single-Variable-Branching (SVB) model is a simplified model of 31 * Branch & Bound trees, from which several nontrivial variable selection 32 * rules arise. The Treemodel branching rule complements SCIP's hybrid 33 * branching by suggesting improved branching variables given the current 34 * pseudocosts and the current dual gap. 35 * 36 * Given a variable with dual bound changes (l, r) (both positive) 37 * and an absolute gap G, the SVB model describes the tree that needs to be 38 * built by branching on that same variable at every node until the value G 39 * is reached at every leaf, starting from 0 at the root node. 40 * If we do so for every variable, we can select the variable that produces 41 * the smallest tree. 42 * In the case where the gap is not known, then we can compute the growth rate 43 * of the tree, which we call the ratio. 44 * The ratio of a variable (l, r) is the factor by which the size of the tree 45 * built using (l, r) that closes a gap G must be multiplied by to close a gap 46 * G+1. This ratio is not constant for all gaps, but when G tends to infinity, 47 * it converges to a fixed value we can compute numerically using a root finding 48 * algorithm (e.g. Laguerre). 49 * The ratio is used when the gap is too large (e.g. no primal bound known) or 50 * to help approximate the size of the SVB tree for that variable. 51 * 52 * See the following publication for more detail: 53 * 54 * @par 55 * Pierre Le Bodic and George Nemhauser@n 56 * An abstract model for branching and its application to mixed integer programming@n 57 * Mathematical Programming, 2017@n 58 */ 59 60 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/ 61 62 #include "scip/treemodel.h" 63 64 #include "scip/history.h" 65 #include "scip/var.h" 66 67 #include <limits.h> 68 69 #define LAGUERRE_THRESHOLD 100 /**< Maximum value of r/l at which Laguerre is the prefered FP method */ 70 71 /* Default parameters for the Treemodel branching rules */ 72 #define DEFAULT_ENABLE FALSE /**< should candidate branching variables be scored using the Treemodel rule? */ 73 #define DEFAULT_HIGHRULE 'r' /**< scoring function to use at nodes predicted to be high in the tree. 74 * ('d'efault, 's'vts, 'r'atio, 't'ree sample) */ 75 #define DEFAULT_LOWRULE 'r' /**< scoring function to use at nodes predicted to be low in the tree 76 * ('d'efault, 's'vts, 'r'atio, 't'ree sample) */ 77 #define DEFAULT_HEIGHT 10 /**< estimated tree height at which we switch from using the low rule to 78 * the high rule */ 79 #define DEFAULT_FILTERHIGH 'a' /**< should dominated candidates be filtered before using the high scoring 80 * function? ('a'uto, 't'rue, 'f'alse) */ 81 #define DEFAULT_FILTERLOW 'a' /**< should dominated candidates be filtered before using the low scoring 82 * function? ('a'uto, 't'rue, 'f'alse) */ 83 #define DEFAULT_MAXFPITER 24 /**< maximum number of fixed-point iterations when computing the ratio */ 84 #define DEFAULT_MAXSVTSHEIGHT 100 /**< maximum height to compute the SVTS score exactly before approximating */ 85 #define DEFAULT_FALLBACKINF 'r' /**< which method should be used as a fallback if the tree size estimates are 86 * infinite? ('d'efault, 'r'atio) */ 87 #define DEFAULT_FALLBACKNOPRIM 'r' /**< which method should be used as a fallback if there is no primal bound 88 * available? ('d'efault, 'r'atio) */ 89 #define DEFAULT_SMALLPSCOST 0.1 /**< threshold at which pseudocosts are considered small, making hybrid scores 90 * more likely to be the deciding factor in branching */ 91 92 /** parameters required by the Treemodel branching rules */ 93 struct SCIP_Treemodel 94 { 95 SCIP_Bool enabled; /**< should candidate branching variables be scored using the Treemodel 96 * rule? */ 97 char highrule; /**< scoring function to use at nodes predicted to be high in the tree. 98 * ('d'efault, 's'vts, 'r'atio, 't'ree sample) */ 99 char lowrule; /**< scoring function to use at nodes predicted to be low in the tree 100 * ('d'efault, 's'vts, 'r'atio, 't'ree sample) */ 101 int height; /**< estimated tree height at which we switch from using the low rule to 102 * the high rule */ 103 char filterhigh; /**< should dominated candidates be filtered before using the high 104 * scoring function? ('a'uto, 't'rue, 'f'alse) [ADVANCED] */ 105 char filterlow; /**< should dominated candidates be filtered before using the low 106 * scoring function? ('a'uto, 't'rue, 'f'alse) [ADVANCED] */ 107 int maxfpiter; /**< maximum number of fixed-point iterations when computing the ratio 108 * [ADVANCED] */ 109 int maxsvtsheight; /**< maximum height to compute the SVTS score exactly before approximating 110 * [ADVANCED] */ 111 char fallbackinf; /**< which method should be used as a fallback if the tree size estimates 112 * are infinite? ('d'efault, 'r'atio) [ADVANCED] */ 113 char fallbacknoprim; /**< which method should be used as a fallback if there is no primal bound 114 * available? ('d'efault, 'r'atio) [ADVANCED] */ 115 SCIP_Real smallpscost; /**< threshold at which pseudocosts are considered small, making hybrid 116 * scores more likely to be the deciding factor in branching [ADVANCED] */ 117 }; 118 119 /** branching encoding of a variable's ratio 120 * A variable's ratio is defined based upon its left and right LP gains, as the unique root > 1 of 121 * the polynomial x^r - x^(r-l) -1, where l and r are the left and right LP gains. 122 * We store the root as upratio^(invleft), with invleft = 1/l. The value upratio is thus 123 * the ratio of the variable (1, r/l). 124 * An extra boolean stores whether the encoded ratio is valid, 125 * i.e. there were no numerical problems when computing it */ 126 struct SCIP_Ratio 127 { 128 SCIP_Real upratio; /**< "UnPowered" ratio, i.e. the ratio of the characteristic polynomial 129 * with gains (1, rightgain/leftgain) */ 130 SCIP_Real invleft; /**< "INVerse left gain, i.e. 1/leftgain */ 131 SCIP_Bool valid; /**< True iff the ratio computed is valid */ 132 }; 133 typedef struct SCIP_Ratio SCIP_RATIO; 134 135 /** a comparison method for the next method. It simply compares two SCIP_Real */ 136 static 137 SCIP_DECL_SORTINDCOMP(sciprealcomp) 138 { 139 SCIP_Real* value = (SCIP_Real*) dataptr; 140 SCIP_Real diffval; 141 142 assert(value != NULL); 143 assert(ind1 >= 0 && ind2 >= 0); 144 145 diffval = value[ind1] - value[ind2]; 146 if( diffval < 0.0 ) 147 return -1; 148 else if( diffval > 0.0) 149 return 1; 150 else 151 return 0; 152 } 153 154 /** given a pair of arrays of real non-negative values (a,b), with a <= b, computes 155 * the pairs that belong to the pareto front (with a tolerance). 156 * In other words, we are looking for non-dominated pairs of values. 157 * One value and one array are computed after this method. 158 * The value is the number of non-dominated elements. 159 * The array is a boolean array that indicates if an element is dominated. 160 * In case of a draw, only one variable is considered as non-dominated. 161 */ 162 static 163 SCIP_RETCODE findNonDominatedVars( 164 SCIP* scip, /**< SCIP data structure */ 165 SCIP_Real* a, /**< the first set of values */ 166 SCIP_Real* b, /**< the second set of values */ 167 int size, /**< the size of array a (and b) */ 168 int* ndominated, /**< returns the number of dominated elements */ 169 SCIP_Bool* dominated /**< returns the array of booleans that determine if an element is 170 * dominated */ 171 ) 172 { 173 SCIP_Real bestcurrenta; 174 SCIP_Real besta; 175 SCIP_Real currentb; 176 int* permb; 177 int* bestcurrents; 178 int nbestcurrent; 179 int indexinpermb; 180 int origindex; 181 int iterbestcurrent; 182 183 assert(scip != NULL); 184 assert(a != NULL); 185 assert(b != NULL); 186 assert(ndominated != NULL); 187 assert(dominated != NULL); 188 assert(size > 0); 189 190 SCIP_CALL( SCIPallocBufferArray(scip, &bestcurrents, size) ); 191 192 /* we first find the permutation of indices of array b that corresponds to 193 * the array of a non-increasing sort of its values */ 194 SCIP_CALL( SCIPallocBufferArray(scip, &permb, size) ); 195 for( origindex=0; origindex<size; ++origindex ) 196 permb[origindex] = origindex; 197 198 SCIPsortDownInd(permb, sciprealcomp, (void*)b, size); 199 200 *ndominated = 0; 201 /* Now we will traverse the pair of arrays a and b by non-decreasing order of values of b 202 * and mark the (non) dominated pairs */ 203 204 /* The current max value of a for all pairs that (almost) have the same value b */ 205 bestcurrenta = a[permb[0]]; 206 207 /* the current value b */ 208 currentb = b[permb[0]]; 209 /* the best pair(s) for the current value b */ 210 bestcurrents[0] = permb[0]; 211 nbestcurrent = 1; 212 /* the value a to "beat" to be non-dominated */ 213 besta = -1; 214 for( indexinpermb = 1; indexinpermb < size; ++indexinpermb ) 215 { 216 origindex = permb[indexinpermb]; 217 assert(b[origindex] <= currentb); 218 if( SCIPisLT(scip, b[origindex], currentb) ) 219 { 220 /* If the above is true, then we went through all the previous elements that had value currentb */ 221 /* Thus the best element for value currentb is non-dominated if its value bestcurrenta is better 222 * than the previous best besta */ 223 if( bestcurrenta > besta ) 224 { 225 for( iterbestcurrent=0; iterbestcurrent < nbestcurrent; ++iterbestcurrent ) 226 dominated[bestcurrents[iterbestcurrent]] = FALSE; 227 228 besta = bestcurrenta; 229 } 230 else 231 { 232 for( iterbestcurrent = 0; iterbestcurrent < nbestcurrent; ++iterbestcurrent ) 233 { 234 dominated[bestcurrents[iterbestcurrent]] = TRUE; 235 ++(*ndominated); 236 } 237 } 238 bestcurrenta = a[origindex]; 239 currentb = b[origindex]; 240 bestcurrents[0] = origindex; 241 nbestcurrent = 1; 242 } 243 else 244 { 245 /* Then the b values are (almost) equal and we need to compare values a */ 246 if( SCIPisGT(scip, a[origindex], bestcurrenta) ) 247 { 248 /* Then the new value is better than the old one(s) */ 249 for( iterbestcurrent = 0; iterbestcurrent < nbestcurrent; ++iterbestcurrent ) 250 { 251 dominated[bestcurrents[iterbestcurrent]] = TRUE; 252 ++(*ndominated); 253 } 254 255 bestcurrenta = a[origindex]; 256 bestcurrents[0] = origindex; 257 nbestcurrent = 1; 258 } 259 else 260 { 261 /* Then the new value is equal or dominated */ 262 if( SCIPisEQ(scip, a[origindex], bestcurrenta) ) 263 { 264 bestcurrents[nbestcurrent] = origindex; 265 ++nbestcurrent; 266 } 267 else 268 { 269 dominated[origindex] = TRUE; 270 ++(*ndominated); 271 } 272 } 273 } 274 } 275 /* Finally, we have to look at the last best variable */ 276 if( bestcurrenta > besta ) 277 { 278 for( iterbestcurrent = 0; iterbestcurrent < nbestcurrent; ++iterbestcurrent ) 279 dominated[bestcurrents[iterbestcurrent]] = FALSE; 280 } 281 else 282 { 283 for( iterbestcurrent = 0; iterbestcurrent < nbestcurrent; ++iterbestcurrent ) 284 { 285 dominated[bestcurrents[iterbestcurrent]] = TRUE; 286 ++(*ndominated); 287 } 288 } 289 290 SCIPfreeBufferArray(scip, &permb); 291 SCIPfreeBufferArray(scip, &bestcurrents); 292 return SCIP_OKAY; 293 } 294 295 /** returns true iff the variable with given gains has a ratio better (i.e smaller) than the given one */ 296 static 297 SCIP_Bool hasBetterRatio( 298 SCIP* scip, /**< SCIP data structure */ 299 SCIP_RATIO* branchratio, /**< The variable's ratio to compare against */ 300 SCIP_Real leftgain, /**< the left gain of a variable */ 301 SCIP_Real rightgain /**< the right gain of a variable */ 302 ) 303 { 304 SCIP_Real result; 305 306 assert(branchratio != NULL); 307 assert(branchratio->valid); 308 assert(SCIPisLE(scip, leftgain, rightgain)); 309 310 /* We evaluate the characteristic polynomial of the variable on the given ratio. */ 311 result = -1; 312 if( leftgain > 0.0 && rightgain > 0.0 ) 313 { 314 result = pow(branchratio->upratio, rightgain * branchratio->invleft) - pow(branchratio->upratio, (rightgain - leftgain) * branchratio->invleft) - 1; /*lint !e644*/ 315 } 316 317 /* If the result is positive, then it has a better ratio. */ 318 return (result > 0.0); 319 } 320 321 /** computes the variable ratio corresponding to the left and right gains */ 322 static 323 void computeVarRatio( 324 SCIP* scip, /**< SCIP data structure */ 325 SCIP_TREEMODEL* treemodel, /**< Treemodel parameter data structure */ 326 SCIP_VAR* var, /**< the candidate branching variable */ 327 SCIP_Real leftgain, /**< the left gain of the variable */ 328 SCIP_Real rightgain, /**< the right gain of the variable */ 329 SCIP_RATIO* branchratio /**< storage for the computed ratio */ 330 ) 331 { 332 SCIP_Real ratio; 333 SCIP_Real newratio; 334 SCIP_Real r; 335 int iters; 336 337 assert(SCIPisGE(scip, leftgain, 0.0)); 338 assert(SCIPisGE(scip, rightgain, leftgain)); 339 340 if( SCIPisZero(scip, leftgain) || SCIPisZero(scip, rightgain) ) 341 { 342 branchratio->valid = FALSE; 343 return; 344 } 345 346 /* We scale left and right gains by dividing both by left */ 347 r = rightgain / leftgain; 348 349 /* In the case where l and r are very close r may become < 1 */ 350 if( r <= 1 ) 351 { 352 branchratio->valid = TRUE; 353 branchratio->upratio = 2.0; 354 branchratio->invleft = 1.0 / leftgain; 355 return; 356 } 357 358 /* Check if this ratio has already been computed */ 359 if( SCIPhistoryIsRatioValid(var->history) && SCIPisEQ(scip, SCIPhistoryGetLastBalance(var->history), r) ) 360 { 361 branchratio->valid = TRUE; 362 branchratio->upratio = SCIPhistoryGetLastRatio(var->history); 363 branchratio->invleft = 1.0 / leftgain; 364 return; 365 } 366 367 /* Initialise the ratio at the previously computed ratio (if applicable) otherwise 368 * use the lower bound 2^(1/r) <= phi <= 2^(1/l). 369 * Note that we only use the previous ratio if the previous value of r/l was larger, 370 * ie. the previous ratio was smaller since we want to initialise at a lower bound. 371 */ 372 ratio = 1.0; 373 newratio = pow(2.0, 1.0/r); 374 if( SCIPhistoryIsRatioValid(var->history) && SCIPhistoryGetLastBalance(var->history) > r 375 && SCIPhistoryGetLastRatio(var->history) > newratio ) 376 newratio = SCIPhistoryGetLastRatio(var->history); 377 378 /* Depending on the value of rightgain/leftgain, we have two different methods to compute the ratio 379 * If this value is bigger than 100, we use a fixed-point method. Otherwise, we use Laguerre's method 380 * This is strictly for numerical efficiency and based on experiments. 381 */ 382 383 /* Use Laguerre's method */ 384 if( r <= LAGUERRE_THRESHOLD ) 385 { 386 /* We relax the epsilon after 5 iterations since we may not have enough precision to achieve any better 387 * convergence */ 388 for( iters = 0; ((iters <= 5 && !SCIPisEQ(scip, ratio, newratio)) || 389 (iters > 5 && !SCIPisSumEQ(scip, ratio, newratio))) 390 && iters < treemodel->maxfpiter && newratio > 1.0; iters++ ) 391 { 392 double G, H, a, p, p1, p2, phi_r; 393 394 ratio = newratio; 395 phi_r = pow(ratio, r); 396 p = phi_r - phi_r / ratio - 1.0; 397 if( p != 0 ) 398 { 399 p1 = (r * phi_r - (r - 1.0) * phi_r / ratio) / ratio; 400 p2 = (r * (r - 1.0) * phi_r - (r - 1.0) * (r - 2.0) * phi_r / ratio) / ratio / ratio; 401 G = p1 / p; 402 H = G * G - (p2 / p); 403 a = r / (G + (G >= 0 ? 1.0 : -1.0) * sqrt((r - 1.0) * (r * H - G * G))); 404 newratio = ratio - a; 405 } 406 } 407 } 408 /* Use fixed point method */ 409 else 410 { 411 /* We relax the epsilon after 10 iterations since we may not have enough precision to achieve any better 412 * convergence */ 413 for( iters = 0; ((iters <= 10 && !SCIPisEQ(scip, ratio, newratio)) || 414 (iters > 10 && !SCIPisSumEQ(scip, ratio, newratio))) 415 && iters < treemodel->maxfpiter && newratio > 1; iters++ ) 416 { 417 ratio = newratio; 418 newratio = pow(1.0-1.0/ratio, -1.0/r); 419 } 420 } 421 422 /* We think that everything worked. 423 * Note that the fixed point method is not guaranteed to converge due to numerical precision issues. 424 * In the case that the method fails to converge, a fallback strategy must be used. 425 * For instance, if this method is used for branching, then this variable can be ignored, 426 * or the scores of all variables could be recomputed using a different method. */ 427 if( iters < treemodel->maxfpiter && newratio > 1.0 ) 428 { 429 branchratio->valid = TRUE; 430 branchratio->upratio = (ratio + newratio) / 2.0; 431 branchratio->invleft = 1.0 / leftgain; 432 } 433 /* We (hopefully) make finding bugs easier by setting these values */ 434 else 435 { 436 branchratio->valid = FALSE; 437 branchratio->upratio = -1.0; 438 branchratio->invleft = -1.0; 439 } 440 441 /* Update the history */ 442 SCIPhistorySetRatioHistory(var->history, branchratio->valid, branchratio->upratio, r); 443 } 444 445 /** use the Ratio scoring function to select a branching candidate */ 446 static 447 SCIP_RETCODE selectCandidateUsingRatio( 448 SCIP* scip, /**< SCIP data structure */ 449 SCIP_TREEMODEL* treemodel, /**< Treemodel parameter data structure */ 450 SCIP_VAR** branchcands, /**< branching candidate storage */ 451 SCIP_Real* mingains, /**< minimum gain of rounding downwards or upwards */ 452 SCIP_Real* maxgains, /**< maximum gain of rounding downwards or upwards */ 453 SCIP_Bool filterdominated, /**< whether dominated variables have been filtered */ 454 SCIP_Bool* dominated, /**< whether each variable is dominated or not */ 455 int nbranchcands, /**< the number of branching candidates */ 456 int* bestcand /**< the best branching candidate found before the call, 457 and the best candidate after the call (possibly the same) */ 458 ) 459 { 460 SCIP_RATIO branchratio; 461 SCIP_RATIO bestbranchratio; 462 int c; 463 464 /* We initialize bestbranchratio at the default bestcand ratio, since it is likely to have 465 * a very good ratio and save evaluations of the ratio for many variables */ 466 int referencevar = *bestcand; 467 computeVarRatio(scip, treemodel, branchcands[referencevar], mingains[referencevar], maxgains[referencevar], &bestbranchratio); 468 469 for( c = 0; c < nbranchcands; ++c ) 470 { 471 if( (!filterdominated || !dominated[c]) && c != referencevar ) 472 { 473 if( !bestbranchratio.valid || hasBetterRatio(scip, &bestbranchratio, mingains[c], maxgains[c]) ) /*lint !e644*/ 474 { 475 computeVarRatio(scip, treemodel, branchcands[c], mingains[c], maxgains[c], &branchratio); 476 if( branchratio.valid ) /*lint !e644*/ 477 { 478 *bestcand = c; 479 bestbranchratio = branchratio; 480 } 481 } 482 } 483 } 484 485 return SCIP_OKAY; 486 } 487 488 /** Returns the predicted treesize for the gap and given up and down gains */ 489 static 490 SCIP_Real computeSVTS( 491 SCIP* scip, /**< SCIP data structure */ 492 SCIP_TREEMODEL* treemodel, /**< Treemodel parameter data structure */ 493 SCIP_VAR* var, /**< the candidate branching variable */ 494 SCIP_Real absgap, /**< the absolute gap to close (typically the local gap at the current node) */ 495 SCIP_Real mingain, /**< prediction of smaller objective gain of downwards/upwards */ 496 SCIP_Real maxgain /**< prediction of larger objective gain of downwards/upwards */ 497 ) 498 { 499 SCIP_Real prediction = SCIP_REAL_MAX; 500 501 if( SCIPisGT(scip, mingain, 0.0) && !SCIPisInfinity(scip, absgap) ) 502 { 503 SCIP_Real treesize; 504 SCIP_Real gaptoreach; 505 SCIP_Real scaledgap; 506 SCIP_Real scaledgain; 507 int mindepth; 508 int nr; 509 int ir; 510 511 /* We implicitly set the minimum gain to 1, and the maximum gain and gap accordingly, 512 * as the treesize does not change if we scale the gains and gap by a scalar */ 513 scaledgain = maxgain / mingain; 514 scaledgap = absgap / mingain; 515 516 mindepth = (int) SCIPceil(scip, scaledgap / scaledgain); 517 518 /* In the following case we compute the treesize for a smaller gap 519 * and we will deduce the treesize of the scaledgap using the ratio */ 520 if( mindepth > treemodel->maxsvtsheight ) 521 { 522 gaptoreach = scaledgap * (treemodel->maxsvtsheight - 1) / mindepth; 523 524 assert(!SCIPisInfinity(scip, gaptoreach)); 525 assert(gaptoreach > scaledgain); 526 } 527 else 528 { 529 gaptoreach = scaledgap; 530 } 531 532 mindepth = (int) ceil(gaptoreach / scaledgain); 533 assert(mindepth <= treemodel->maxsvtsheight); 534 treesize = 1; 535 536 /* nr is the number of times we turn right to reach a leaf */ 537 for( nr = 1; nr <= mindepth; ++nr ) 538 { 539 SCIP_Real binomcoeff = 1.0; 540 for( ir = 1; ir <= nr; ++ir ) 541 { 542 binomcoeff *= (nr + ceil((gaptoreach - (nr - 1) * scaledgain)) - ir) / ir; 543 } 544 treesize += binomcoeff; 545 } 546 547 treesize = 2.0 * treesize - 1.0; 548 549 assert(SCIPisGE(scip, treesize, 3.0)); 550 551 if( !SCIPisEQ(scip, scaledgap, gaptoreach) ) 552 { 553 /* If we have not computed the treesize for the scaled gap but for max gap instead, 554 * we use the ratio between two iterations to come up with an estimate of the treesize 555 * for the scaled gap */ 556 if( !SCIPisInfinity(scip,treesize) ) 557 { 558 SCIP_RATIO branchratio; 559 computeVarRatio(scip, treemodel, var, mingain, maxgain, &branchratio); 560 561 if( branchratio.valid ) /*lint !e644*/ 562 prediction = treesize * pow(branchratio.upratio, (scaledgap - gaptoreach) * branchratio.invleft); /*lint !e644*/ 563 } 564 } 565 else 566 { 567 prediction = treesize; 568 } 569 } 570 571 return prediction; 572 } 573 574 /** use the SVTS scoring function to select a branching candidate */ 575 static 576 SCIP_RETCODE selectCandidateUsingSVTS( 577 SCIP* scip, /**< SCIP data structure */ 578 SCIP_TREEMODEL* treemodel, /**< Treemodel parameter data structure */ 579 SCIP_VAR** branchcands, /**< branching candidate storage */ 580 SCIP_Real* mingains, /**< minimum gain of rounding downwards or upwards */ 581 SCIP_Real* maxgains, /**< maximum gain of rounding downwards or upwards */ 582 SCIP_Real* tiebreakerscore, /**< scores to use for tie breaking */ 583 SCIP_Real localabsgap, /**< The dual gap at the current node */ 584 SCIP_Bool filterdominated, /**< whether dominated variables have been filtered */ 585 SCIP_Bool* dominated, /**< whether each variable is dominated or not */ 586 int nbranchcands, /**< the number of branching candidates */ 587 int ndominated, /**< the number of dominated candidates */ 588 int* bestcand /**< the best branching candidate found before the call, 589 and the best candidate after the call (possibly the same) */ 590 ) 591 { 592 SCIP_Real* treesizes; 593 SCIP_Real referencetreesize; 594 SCIP_Real score; 595 SCIP_Real bestscore; 596 SCIP_Real avgtreesize; 597 int besttscand; 598 int referencevar; 599 int c; 600 601 /* We will first measure the treesize for scip's default variable. If it is infinite then we don't compute 602 * the treesize for other variables (even though it might be finite) and go directly to the fallback strategy */ 603 besttscand = *bestcand; 604 referencevar = *bestcand; 605 606 treesizes = NULL; 607 bestscore = 0.0; 608 avgtreesize = 0.0; 609 if( !SCIPisInfinity(scip, localabsgap) ) 610 { 611 referencetreesize = computeSVTS(scip, treemodel, branchcands[referencevar], localabsgap, mingains[referencevar], 612 maxgains[referencevar]); 613 if( !SCIPisInfinity(scip, referencetreesize) ) 614 { 615 SCIP_CALL( SCIPallocBufferArray(scip, &treesizes, nbranchcands) ); 616 treesizes[referencevar] = referencetreesize; 617 618 for( c = 0; c < nbranchcands; ++c ) 619 { 620 if( !filterdominated || !dominated[c] ) 621 { 622 if( c != referencevar ) 623 treesizes[c] = computeSVTS(scip, treemodel, branchcands[c], localabsgap, mingains[c], maxgains[c]); 624 else 625 treesizes[c] = referencetreesize; 626 627 avgtreesize += treesizes[c]; 628 } 629 else 630 treesizes[c] = SCIP_REAL_MAX; 631 } 632 avgtreesize = avgtreesize / (nbranchcands - ndominated); 633 634 for( c = 0; c < nbranchcands; ++c ) 635 { 636 score = (1.0 - 1.0 / (1.0 + avgtreesize / treesizes[c])) + 0.01 * tiebreakerscore[c]; 637 if(score > bestscore) 638 { 639 bestscore = score; 640 besttscand = c; 641 } 642 } 643 644 *bestcand = besttscand; 645 646 SCIPfreeBufferArray(scip, &treesizes); 647 } 648 /* Apply infinite treesize fallback strategy */ 649 else if( treemodel->fallbackinf == 'r' ) 650 { 651 SCIP_CALL( selectCandidateUsingRatio(scip, treemodel, branchcands, mingains, maxgains, filterdominated, dominated, 652 nbranchcands, bestcand) ); 653 } 654 } 655 /* Apply no primal bound fallback strategy */ 656 else if( treemodel->fallbacknoprim == 'r' ) 657 { 658 SCIP_CALL( selectCandidateUsingRatio(scip, treemodel, branchcands, mingains, maxgains, filterdominated, dominated, 659 nbranchcands, bestcand) ); 660 } 661 662 return SCIP_OKAY; 663 } 664 665 /** computes a^b for integer b */ 666 static 667 SCIP_Real integerpow( 668 SCIP_Real a, /**< the base */ 669 int b /**< the integer exponent */ 670 ) 671 { /*lint --e{644}*/ 672 SCIP_Real ans; 673 674 ans = 1.0; 675 for( ; b; b /= 2 ) 676 { 677 if( b & 1 ) 678 ans *= a; 679 a *= a; 680 } 681 return ans; 682 } 683 684 /** returns the sampled tree size for the given lp gains and dual gap */ 685 static 686 SCIP_Real computeSampleTreesize( 687 SCIP* scip, /**< SCIP data structure */ 688 SCIP_TREEMODEL* treemodel, /**< Treemodel parameter data structure */ 689 SCIP_VAR* var, /**< the candidate branching variable */ 690 SCIP_Real absgap, /**< the absolute gap to close (typically the local at the current node) */ 691 SCIP_Real leftgain, /**< The minimum gain from branching on this variable */ 692 SCIP_Real rightgain /**< The maximum gain from branching on this variable */ 693 ) 694 { 695 SCIP_RATIO branchratio; 696 SCIP_Real prediction; 697 SCIP_Real leftsize; 698 SCIP_Real rightsize; 699 SCIP_Real midsize; 700 701 computeVarRatio(scip, treemodel, var, leftgain, rightgain, &branchratio); 702 703 if( branchratio.valid ) /*lint !e644*/ 704 { /*lint --e{644}*/ 705 SCIP_Real phi_l = branchratio.upratio; 706 SCIP_Real phi_r = pow(branchratio.upratio, rightgain * branchratio.invleft); 707 int kl = (int)ceil(absgap / leftgain); 708 int kr = (int)ceil(absgap / rightgain); 709 int k = (int)ceil(absgap / (leftgain + rightgain)); 710 SCIP_Real phi_lr = phi_l * phi_r; 711 SCIP_Real phi_klr = integerpow(phi_lr, k); 712 713 /* We compute an estimate of the size of the tree using the left-most leaf, 714 * right-most leaf, and the leaf obtained from alternating left and right. */ 715 leftsize = (integerpow(phi_l, kl + 1) - 1.0) / (phi_l - 1.0); 716 rightsize = (integerpow(phi_r, kr + 1) - 1.0) / (phi_r - 1.0); 717 718 if( k * (leftgain + rightgain) < absgap + rightgain ) 719 midsize = (1.0 + phi_l) * (phi_klr * phi_lr - 1.0) / (phi_lr - 1.0) - phi_klr * phi_l; 720 else 721 midsize = (1.0 + phi_l) * (phi_klr - 1.0) / (phi_lr - 1.0); 722 723 prediction = (leftsize + rightsize + midsize) / 3.0; 724 } 725 else 726 { 727 prediction = SCIP_REAL_MAX; 728 } 729 730 return prediction; 731 } 732 733 /** use the Tree Sampling scoring function to select a branching candidate */ 734 static 735 SCIP_RETCODE selectCandidateUsingSampling( 736 SCIP* scip, /**< SCIP data structure */ 737 SCIP_TREEMODEL* treemodel, /**< Treemodel parameter data structure */ 738 SCIP_VAR** branchcands, /**< branching candidate storage */ 739 SCIP_Real* mingains, /**< minimum gain of rounding downwards or upwards */ 740 SCIP_Real* maxgains, /**< maximum gain of rounding downwards or upwards */ 741 SCIP_Real* tiebreakerscore, /**< scores to use for tie breaking */ 742 SCIP_Real localabsgap, /**< The dual gap at the current node */ 743 SCIP_Bool filterdominated, /**< whether dominated variables have been filtered */ 744 SCIP_Bool* dominated, /**< whether each variable is dominated or not */ 745 int nbranchcands, /**< the number of branching candidates */ 746 int ndominated, /**< the number of dominated candidates */ 747 int* bestcand /**< the best branching candidate found before the call, 748 and the best candidate after the call (possibly the same) */ 749 ) 750 { 751 SCIP_Real* treesizes; 752 SCIP_Real referencetreesize; 753 SCIP_Real score; 754 SCIP_Real bestscore; 755 SCIP_Real avgtreesize; 756 int besttscand; 757 int referencevar; 758 int c; 759 760 /* We will first measure the treesize for scip's default variable. If it is infinite then we don't compute 761 * the treesize for other variables (even though it might be finite) and go directly to the fallback strategy */ 762 besttscand = *bestcand; 763 referencevar = *bestcand; 764 765 treesizes = NULL; 766 bestscore = 0.0; 767 avgtreesize = 0.0; 768 if( !SCIPisInfinity(scip, localabsgap) ) 769 { 770 referencetreesize = computeSampleTreesize(scip, treemodel, branchcands[referencevar], localabsgap, mingains[referencevar], 771 maxgains[referencevar]); 772 773 if( !SCIPisInfinity(scip, referencetreesize) ) 774 { 775 SCIP_CALL( SCIPallocBufferArray(scip, &treesizes, nbranchcands) ); 776 treesizes[referencevar] = referencetreesize; 777 778 for( c = 0; c < nbranchcands; ++c ) 779 { 780 if( !filterdominated || !dominated[c] ) 781 { 782 if( c != referencevar ) 783 treesizes[c] = computeSampleTreesize(scip, treemodel, branchcands[c], localabsgap, mingains[c], maxgains[c]); 784 else 785 treesizes[c] = referencetreesize; 786 787 avgtreesize += treesizes[c]; 788 } 789 else 790 treesizes[c] = SCIP_REAL_MAX; 791 } 792 avgtreesize = avgtreesize / (nbranchcands - ndominated); 793 794 for( c = 0; c < nbranchcands; ++c ) 795 { 796 score = (1.0 - 1.0 / (1.0 + avgtreesize / treesizes[c])) + 0.01 * tiebreakerscore[c]; 797 if( score > bestscore ) 798 { 799 bestscore = score; 800 besttscand = c; 801 } 802 } 803 804 *bestcand = besttscand; 805 806 SCIPfreeBufferArray(scip, &treesizes); 807 } 808 /* Apply infinite treesize fallback strategy */ 809 else if( treemodel->fallbackinf == 'r' ) 810 { 811 SCIP_CALL( selectCandidateUsingRatio(scip, treemodel, branchcands, mingains, maxgains, filterdominated, dominated, 812 nbranchcands, bestcand) ); 813 } 814 } 815 /* Apply no primal bound fallback strategy */ 816 else if( treemodel->fallbacknoprim == 'r' ) 817 { 818 SCIP_CALL( selectCandidateUsingRatio(scip, treemodel, branchcands, mingains, maxgains, filterdominated, dominated, 819 nbranchcands, bestcand) ); 820 } 821 822 return SCIP_OKAY; 823 } 824 825 /** initialises the Treemodel parameter data structure */ 826 SCIP_RETCODE SCIPtreemodelInit( 827 SCIP* scip, /**< SCIP data structure */ 828 SCIP_TREEMODEL** treemodel /**< Treemodel parameter data structure */ 829 ) 830 { 831 assert(treemodel != NULL); 832 SCIP_CALL( SCIPallocBlockMemory(scip, treemodel) ); 833 assert(*treemodel != NULL); 834 835 SCIP_CALL( SCIPaddBoolParam(scip, "branching/treemodel/enable", 836 "should candidate branching variables be scored using the Treemodel branching rules?", 837 &(*treemodel)->enabled, FALSE, DEFAULT_ENABLE, 838 NULL, NULL) ); 839 SCIP_CALL( SCIPaddCharParam(scip, "branching/treemodel/highrule", 840 "scoring function to use at nodes predicted to be high in the tree ('d'efault, 's'vts, 'r'atio, 't'ree sample)", 841 &(*treemodel)->highrule, FALSE, DEFAULT_HIGHRULE, "dsrt", 842 NULL, NULL) ); 843 SCIP_CALL( SCIPaddCharParam(scip, "branching/treemodel/lowrule", 844 "scoring function to use at nodes predicted to be low in the tree ('d'efault, 's'vts, 'r'atio, 't'ree sample)", 845 &(*treemodel)->lowrule, FALSE, DEFAULT_LOWRULE, "dsrt", 846 NULL, NULL) ); 847 SCIP_CALL( SCIPaddIntParam(scip, "branching/treemodel/height", 848 "estimated tree height at which we switch from using the low rule to the high rule", 849 &(*treemodel)->height, FALSE, DEFAULT_HEIGHT, 0, INT_MAX, 850 NULL, NULL) ); 851 SCIP_CALL( SCIPaddCharParam(scip, "branching/treemodel/filterhigh", 852 "should dominated candidates be filtered before using the high scoring function? ('a'uto, 't'rue, 'f'alse)", 853 &(*treemodel)->filterhigh, TRUE, DEFAULT_FILTERHIGH, "atf", 854 NULL, NULL) ); 855 SCIP_CALL( SCIPaddCharParam(scip, "branching/treemodel/filterlow", 856 "should dominated candidates be filtered before using the low scoring function? ('a'uto, 't'rue, 'f'alse)", 857 &(*treemodel)->filterlow, TRUE, DEFAULT_FILTERLOW, "atf", 858 NULL, NULL) ); 859 SCIP_CALL( SCIPaddIntParam(scip, "branching/treemodel/maxfpiter", 860 "maximum number of fixed-point iterations when computing the ratio", 861 &(*treemodel)->maxfpiter, TRUE, DEFAULT_MAXFPITER, 1, INT_MAX, 862 NULL, NULL) ); 863 SCIP_CALL( SCIPaddIntParam(scip, "branching/treemodel/maxsvtsheight", 864 "maximum height to compute the SVTS score exactly before approximating", 865 &(*treemodel)->maxsvtsheight, TRUE, DEFAULT_MAXSVTSHEIGHT, 0, INT_MAX, 866 NULL, NULL) ); 867 SCIP_CALL( SCIPaddCharParam(scip, "branching/treemodel/fallbackinf", 868 "which method should be used as a fallback if the tree size estimates are infinite? ('d'efault, 'r'atio)", 869 &(*treemodel)->fallbackinf, TRUE, DEFAULT_FALLBACKINF, "dr", 870 NULL, NULL) ); 871 SCIP_CALL( SCIPaddCharParam(scip, "branching/treemodel/fallbacknoprim", 872 "which method should be used as a fallback if there is no primal bound available? ('d'efault, 'r'atio)", 873 &(*treemodel)->fallbacknoprim, TRUE, DEFAULT_FALLBACKNOPRIM, "dr", 874 NULL, NULL) ); 875 SCIP_CALL ( SCIPaddRealParam(scip, "branching/treemodel/smallpscost", 876 "threshold at which pseudocosts are considered small, making hybrid scores more likely to be the deciding factor in branching", 877 &(*treemodel)->smallpscost, TRUE, DEFAULT_SMALLPSCOST, 878 0.0, SCIP_REAL_MAX, NULL, NULL) ); 879 880 return SCIP_OKAY; 881 } 882 883 /** frees the Treemodel parameter data structure */ 884 SCIP_RETCODE SCIPtreemodelFree( 885 SCIP* scip, /**< SCIP data structure */ 886 SCIP_TREEMODEL** treemodel /**< Treemodel parameter data structure */ 887 ) 888 { 889 assert(treemodel != NULL); 890 assert(*treemodel != NULL); 891 892 SCIPfreeBlockMemory(scip, treemodel); 893 894 assert(*treemodel == NULL); 895 896 return SCIP_OKAY; 897 } 898 899 /** returns TRUE if the Treemodel branching rules are enabled */ 900 SCIP_Bool SCIPtreemodelIsEnabled( 901 SCIP* scip, /**< SCIP data structure */ 902 SCIP_TREEMODEL* treemodel /**< Treemodel parameter data structure */ 903 ) 904 { 905 assert(scip != NULL); 906 return treemodel->enabled; 907 } 908 909 /** apply the Treemodel branching rules to attempt to select a better 910 * branching candidate than the one selected by pseudocost branching 911 */ 912 SCIP_RETCODE SCIPtreemodelSelectCandidate( 913 SCIP* scip, /**< SCIP data structure */ 914 SCIP_TREEMODEL* treemodel, /**< Treemodel parameter data structure */ 915 SCIP_VAR** branchcands, /**< branching candidate storage */ 916 SCIP_Real* mingains, /**< minimum gain of rounding downwards or upwards */ 917 SCIP_Real* maxgains, /**< maximum gain of rounding downwards or upwards */ 918 SCIP_Real* tiebreakerscore, /**< scores to use for tie breaking */ 919 int nbranchcands, /**< the number of branching candidates */ 920 int* bestcand /**< the best branching candidate found before the call, 921 and the best candidate after the call (possibly the same) */ 922 ) 923 { 924 SCIP_Real localabsgap; /* The gap at the current node */ 925 int bestcandheight; /* The height of the best candidate according to SCIP */ 926 char scoringfunction; /* Scoring function to use (based on the estimated tree height) */ 927 char filtersetting; /* Whether we should apply filtering of dominated variables */ 928 929 assert(treemodel != NULL); 930 assert(treemodel->enabled); 931 assert(*bestcand >= 0); 932 933 /* Compute the dual gap at the current node */ 934 if( !SCIPisInfinity(scip, SCIPgetUpperbound(scip)) ) 935 localabsgap = SCIPgetUpperbound(scip) - SCIPgetNodeLowerbound(scip, SCIPgetCurrentNode(scip)); 936 else 937 localabsgap = SCIPinfinity(scip); 938 939 /* Compute an estimate of the height of the current node using the bestcand variable */ 940 if( !SCIPisInfinity(scip, localabsgap) && SCIPisGT(scip, mingains[*bestcand], 0.0) 941 && SCIPisLT(scip, localabsgap/mingains[*bestcand], 1.0 * INT_MAX)) 942 bestcandheight = (int)(localabsgap/mingains[*bestcand]); 943 else 944 bestcandheight = INT_MAX; 945 946 /* Decide which scoring function to use based on the estimated height of the tree */ 947 if( bestcandheight < treemodel->height ) 948 { 949 scoringfunction = treemodel->lowrule; 950 filtersetting = treemodel->filterlow; 951 } 952 else 953 { 954 scoringfunction = treemodel->highrule; 955 filtersetting = treemodel->filterhigh; 956 } 957 958 /* We are going to apply a Treemodel variable selection rule */ 959 if( scoringfunction != 'd' ) 960 { 961 SCIP_Bool* dominated; /* Whether variables are dominated */ 962 SCIP_Bool autofilter; /* If auto filtering is chosen, should variables be filtered? */ 963 SCIP_Bool filterdominated; /* Whether we should filter dominated variables */ 964 int ndominated; /* Number of dominated variables */ 965 966 /* Filtering dominated variables is suggested for SVTS and Tree Sampling rules */ 967 autofilter = (filtersetting == 'a' && (scoringfunction == 's' || scoringfunction == 't')); 968 filterdominated = (autofilter || filtersetting == 't'); 969 970 /* If selected, find the dominated variables */ 971 if( filterdominated ) 972 { 973 SCIP_CALL( SCIPallocBufferArray(scip, &dominated, nbranchcands) ); 974 SCIP_CALL( findNonDominatedVars(scip, mingains, maxgains, nbranchcands, &ndominated, dominated) ); 975 } 976 else 977 { 978 dominated = NULL; 979 ndominated = 0; 980 } 981 982 /* Invoke the selected scoring function */ 983 switch( scoringfunction ) 984 { 985 case 's': 986 SCIP_CALL( selectCandidateUsingSVTS(scip, treemodel, branchcands, mingains, maxgains, tiebreakerscore, 987 localabsgap, filterdominated, dominated, nbranchcands, ndominated, bestcand) ); 988 break; 989 case 'r': 990 SCIP_CALL( selectCandidateUsingRatio(scip, treemodel, branchcands, mingains, maxgains, filterdominated, 991 dominated, nbranchcands, bestcand) ); 992 break; 993 case 't': 994 SCIP_CALL( selectCandidateUsingSampling(scip, treemodel, branchcands, mingains, maxgains, tiebreakerscore, 995 localabsgap, filterdominated, dominated, nbranchcands, ndominated, bestcand) ); 996 break; 997 default: 998 return SCIP_PARAMETERWRONGVAL; 999 } 1000 1001 /* Free dominated variable buffer if it was used */ 1002 if( filterdominated ) 1003 { 1004 assert(dominated != NULL); 1005 SCIPfreeBufferArray(scip, &dominated); 1006 } 1007 } 1008 1009 return SCIP_OKAY; 1010 } 1011