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 branch_distribution.c 26 * @ingroup DEFPLUGINS_BRANCH 27 * @ingroup BRANCHINGRULES 28 * @brief probability based branching rule based on an article by J. Pryor and J.W. Chinneck 29 * @author Gregor Hendel 30 * 31 * The distribution branching rule selects a variable based on its impact on row activity distributions. More formally, 32 * let \f$ a(x) = a_1 x_1 + \dots + a_n x_n \leq b \f$ be a row of the LP. Let further \f$ l_i, u_i \in R\f$ denote the 33 * (finite) lower and upper bound, respectively, of the \f$ i \f$-th variable \f$x_i\f$. 34 * Viewing every variable \f$x_i \f$ as (continuously) uniformly distributed within its bounds, we can approximately 35 * understand the row activity \f$a(x)\f$ as a Gaussian random variate with mean value \f$ \mu = E[a(x)] = \sum_i a_i\frac{l_i + u_i}{2}\f$ 36 * and variance \f$ \sigma^2 = \sum_i a_i^2 \sigma_i^2 \f$, with \f$ \sigma_i^2 = \frac{(u_i - l_i)^2}{12}\f$ for 37 * continuous and \f$ \sigma_i^2 = \frac{(u_i - l_i + 1)^2 - 1}{12}\f$ for discrete variables. 38 * With these two parameters, we can calculate the probability to satisfy the row in terms of the cumulative distribution 39 * of the standard normal distribution: \f$ P(a(x) \leq b) = \Phi(\frac{b - \mu}{\sigma})\f$. 40 * 41 * The impact of a variable on the probability to satisfy a constraint after branching can be estimated by altering 42 * the variable contribution to the sums described above. In order to keep the tree size small, 43 * variables are preferred which decrease the probability 44 * to satisfy a row because it is more likely to create infeasible subproblems. 45 * 46 * The selection of the branching variable is subject to the parameter @p scoreparam. For both branching directions, 47 * an individual score is calculated. Available options for scoring methods are: 48 * - @b d: select a branching variable with largest difference in satisfaction probability after branching 49 * - @b l: lowest cumulative probability amongst all variables on all rows (after branching). 50 * - @b h: highest cumulative probability amongst all variables on all rows (after branching). 51 * - @b v: highest number of votes for lowering row probability for all rows a variable appears in. 52 * - @b w: highest number of votes for increasing row probability for all rows a variable appears in. 53 * 54 * If the parameter @p usescipscore is set to @a TRUE, a single branching score is calculated from the respective 55 * up and down scores as defined by the SCIP branching score function (see advanced branching parameter @p scorefunc), 56 * otherwise, the variable with the single highest score is selected, and the maximizing direction is assigned 57 * higher branching priority. 58 * 59 * The original idea of probability based branching schemes appeared in: 60 * 61 * J. Pryor and J.W. Chinneck:@n 62 * Faster Integer-Feasibility in Mixed-Integer Linear Programs by Branching to Force Change@n 63 * Computers and Operations Research, vol. 38, 2011, p. 1143–1152@n 64 * (Paper: <a href="http://www.sce.carleton.ca/faculty/chinneck/docs/PryorChinneck.pdf">PDF</a>). 65 */ 66 67 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/ 68 69 #include "scip/branch_distribution.h" 70 #include "scip/pub_branch.h" 71 #include "scip/pub_event.h" 72 #include "scip/pub_lp.h" 73 #include "scip/pub_message.h" 74 #include "scip/pub_misc.h" 75 #include "scip/pub_var.h" 76 #include "scip/scip_branch.h" 77 #include "scip/scip_event.h" 78 #include "scip/scip_general.h" 79 #include "scip/scip_lp.h" 80 #include "scip/scip_message.h" 81 #include "scip/scip_mem.h" 82 #include "scip/scip_numerics.h" 83 #include "scip/scip_param.h" 84 #include "scip/scip_pricer.h" 85 #include "scip/scip_prob.h" 86 #include "scip/scip_probing.h" 87 #include <string.h> 88 89 90 #define BRANCHRULE_NAME "distribution" 91 #define BRANCHRULE_DESC "branching rule based on variable influence on cumulative normal distribution of row activities" 92 #define BRANCHRULE_PRIORITY 0 93 #define BRANCHRULE_MAXDEPTH -1 94 #define BRANCHRULE_MAXBOUNDDIST 1.0 95 96 #define SCOREPARAM_VALUES "dhlvw" 97 #define DEFAULT_SCOREPARAM 'v' 98 #define DEFAULT_PRIORITY 2.0 99 #define SQRTOFTWO 1.4142136 100 #define SQUARED(x) ((x) * (x)) 101 #define DEFAULT_ONLYACTIVEROWS FALSE /**< should only rows which are active at the current node be considered? */ 102 #define DEFAULT_USEWEIGHTEDSCORE FALSE /**< should the branching score weigh up- and down-scores of a variable */ 103 104 /* branching rule event handler to catch variable bound changes */ 105 #define EVENTHDLR_NAME "eventhdlr_distribution" /**< event handler name */ 106 #define EVENT_DISTRIBUTION SCIP_EVENTTYPE_BOUNDCHANGED /**< the event tyoe to be handled by this event handler */ 107 108 /* 109 * Data structures 110 */ 111 112 /** branching rule data */ 113 struct SCIP_BranchruleData 114 { 115 SCIP_EVENTHDLR* eventhdlr; /**< event handler pointer */ 116 SCIP_VAR** updatedvars; /**< variables to process bound change events for */ 117 SCIP_Real* rowmeans; /**< row activity mean values for all rows */ 118 SCIP_Real* rowvariances; /**< row activity variances for all rows */ 119 SCIP_Real* currentubs; /**< variable upper bounds as currently saved in the row activities */ 120 SCIP_Real* currentlbs; /**< variable lower bounds as currently saved in the row activities */ 121 int* rowinfinitiesdown; /**< count the number of variables with infinite bounds which allow for always 122 * repairing the constraint right hand side */ 123 int* rowinfinitiesup; /**< count the number of variables with infinite bounds which allow for always 124 * repairing the constraint left hand side */ 125 int* varposs; /**< array of variable positions in the updated variables array */ 126 int* varfilterposs; /**< array of event filter positions for variable events */ 127 int nupdatedvars; /**< the current number of variables with pending bound changes */ 128 int memsize; /**< memory size of current arrays, needed for dynamic reallocation */ 129 int varpossmemsize; /**< memory size of updated vars and varposs array */ 130 char scoreparam; /**< parameter how the branch score is calculated */ 131 SCIP_Bool onlyactiverows; /**< should only rows which are active at the current node be considered? */ 132 SCIP_Bool usescipscore; /**< should the branching use SCIP's branching score function */ 133 }; 134 135 struct SCIP_EventhdlrData 136 { 137 SCIP_BRANCHRULEDATA* branchruledata; /**< the branching rule data to access distribution arrays */ 138 }; 139 140 /* 141 * Local methods 142 */ 143 144 /** ensure that maxindex + 1 rows can be represented in data arrays; memory gets reallocated with 10% extra space 145 * to save some time for future allocations */ 146 static 147 SCIP_RETCODE branchruledataEnsureArraySize( 148 SCIP* scip, /**< SCIP data structure */ 149 SCIP_BRANCHRULEDATA* branchruledata, /**< branchruledata */ 150 int maxindex /**< row index at hand (size must be at least this large) */ 151 ) 152 { 153 int newsize; 154 int r; 155 156 /* maxindex fits in current array -> nothing to do */ 157 if( maxindex < branchruledata->memsize ) 158 return SCIP_OKAY; 159 160 /* new memory size is the max index + 1 plus 10% additional space */ 161 newsize = (int)SCIPfeasCeil(scip, (maxindex + 1) * 1.1); 162 assert(newsize > branchruledata->memsize); 163 assert(branchruledata->memsize >= 0); 164 165 /* alloc memory arrays for row information */ 166 if( branchruledata->memsize == 0 ) 167 { 168 SCIP_VAR** vars; 169 int v; 170 int nvars; 171 172 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &branchruledata->rowinfinitiesdown, newsize) ); 173 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &branchruledata->rowinfinitiesup, newsize) ); 174 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &branchruledata->rowmeans, newsize) ); 175 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &branchruledata->rowvariances, newsize) ); 176 177 assert(SCIPgetStage(scip) == SCIP_STAGE_SOLVING); 178 179 vars = SCIPgetVars(scip); 180 nvars = SCIPgetNVars(scip); 181 182 assert(nvars > 0); 183 184 /* allocate variable update event processing array storage */ 185 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &branchruledata->varfilterposs, nvars) ); 186 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &branchruledata->varposs, nvars) ); 187 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &branchruledata->updatedvars, nvars) ); 188 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &branchruledata->currentubs, nvars) ); 189 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &branchruledata->currentlbs, nvars) ); 190 191 branchruledata->varpossmemsize = nvars; 192 branchruledata->nupdatedvars = 0; 193 194 /* init variable event processing data */ 195 for( v = 0; v < nvars; ++v ) 196 { 197 assert(SCIPvarIsActive(vars[v])); 198 assert(SCIPvarGetProbindex(vars[v]) == v); 199 200 /* set up variable events to catch bound changes */ 201 SCIP_CALL( SCIPcatchVarEvent(scip, vars[v], EVENT_DISTRIBUTION, branchruledata->eventhdlr, NULL, &(branchruledata->varfilterposs[v])) ); 202 assert(branchruledata->varfilterposs[v] >= 0); 203 204 branchruledata->varposs[v] = -1; 205 branchruledata->updatedvars[v] = NULL; 206 branchruledata->currentlbs[v] = SCIP_INVALID; 207 branchruledata->currentubs[v] = SCIP_INVALID; 208 } 209 } 210 else 211 { 212 SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &branchruledata->rowinfinitiesdown, branchruledata->memsize, newsize) ); 213 SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &branchruledata->rowinfinitiesup, branchruledata->memsize, newsize) ); 214 SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &branchruledata->rowmeans, branchruledata->memsize, newsize) ); 215 SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &branchruledata->rowvariances, branchruledata->memsize, newsize) ); 216 } 217 218 /* loop over extended arrays and invalidate data to trigger initialization of this row when necessary */ 219 for( r = branchruledata->memsize; r < newsize; ++r ) 220 { 221 branchruledata->rowmeans[r] = SCIP_INVALID; 222 branchruledata->rowvariances[r] = SCIP_INVALID; 223 branchruledata->rowinfinitiesdown[r] = 0; 224 branchruledata->rowinfinitiesup[r] = 0; 225 } 226 227 /* adjust memsize */ 228 branchruledata->memsize = newsize; 229 230 return SCIP_OKAY; 231 } 232 233 /* update the variables current lower and upper bound */ 234 static 235 void branchruledataUpdateCurrentBounds( 236 SCIP* scip, /**< SCIP data structure */ 237 SCIP_BRANCHRULEDATA* branchruledata, /**< branchrule data */ 238 SCIP_VAR* var /**< the variable to update current bounds */ 239 ) 240 { 241 int varindex; 242 SCIP_Real lblocal; 243 SCIP_Real ublocal; 244 245 assert(var != NULL); 246 247 varindex = SCIPvarGetProbindex(var); 248 assert(0 <= varindex && varindex < branchruledata->varpossmemsize); 249 lblocal = SCIPvarGetLbLocal(var); 250 ublocal = SCIPvarGetUbLocal(var); 251 252 assert(SCIPisFeasLE(scip, lblocal, ublocal)); 253 254 branchruledata->currentlbs[varindex] = lblocal; 255 branchruledata->currentubs[varindex] = ublocal; 256 } 257 258 /** calculate the variable's distribution parameters (mean and variance) for the bounds specified in the arguments. 259 * special treatment of infinite bounds necessary */ 260 void SCIPvarCalcDistributionParameters( 261 SCIP* scip, /**< SCIP data structure */ 262 SCIP_Real varlb, /**< variable lower bound */ 263 SCIP_Real varub, /**< variable upper bound */ 264 SCIP_VARTYPE vartype, /**< type of the variable */ 265 SCIP_Real* mean, /**< pointer to store mean value */ 266 SCIP_Real* variance /**< pointer to store the variance of the variable uniform distribution */ 267 ) 268 { 269 assert(mean != NULL); 270 assert(variance != NULL); 271 272 /* calculate mean and variance of variable uniform distribution before and after branching */ 273 if( SCIPisInfinity(scip, varub) || SCIPisInfinity(scip, -varlb) ) 274 { 275 /* variables with infinite bounds are not kept in the row activity variance */ 276 *variance = 0.0; 277 278 if( !SCIPisInfinity(scip, varub) ) 279 *mean = varub; 280 else if( !SCIPisInfinity(scip, -varlb) ) 281 *mean = varlb; 282 else 283 *mean = 0.0; 284 } 285 else 286 { 287 /* if the variable is continuous, we assume a continuous uniform distribution, otherwise a discrete one */ 288 if( vartype == SCIP_VARTYPE_CONTINUOUS ) 289 *variance = SQUARED(varub - varlb); 290 else 291 *variance = SQUARED(varub - varlb + 1) - 1; 292 *variance /= 12.0; 293 *mean = (varub + varlb) * .5; 294 } 295 296 assert(!SCIPisNegative(scip, *variance)); 297 } 298 299 /** calculates the cumulative distribution P(-infinity <= x <= value) that a normally distributed 300 * random variable x takes a value between -infinity and parameter \p value. 301 * 302 * The distribution is given by the respective mean and deviation. This implementation 303 * uses the error function SCIPerf(). 304 */ 305 SCIP_Real SCIPcalcCumulativeDistribution( 306 SCIP* scip, /**< current SCIP */ 307 SCIP_Real mean, /**< the mean value of the distribution */ 308 SCIP_Real variance, /**< the square of the deviation of the distribution */ 309 SCIP_Real value /**< the upper limit of the calculated distribution integral */ 310 ) 311 { 312 SCIP_Real normvalue; 313 SCIP_Real std; 314 315 /* we need to calculate the standard deviation from the variance */ 316 assert(!SCIPisNegative(scip, variance)); 317 if( SCIPisFeasZero(scip, variance) ) 318 std = 0.0; 319 else 320 std = sqrt(variance); 321 322 /* special treatment for zero variance */ 323 if( SCIPisFeasZero(scip, std) ) 324 { 325 if( SCIPisFeasLE(scip, value, mean) ) 326 return 1.0; 327 else 328 return 0.0; 329 } 330 assert( std != 0.0 ); /* for lint */ 331 332 /* scale and translate to standard normal distribution. Factor sqrt(2) is needed for SCIPerf() function */ 333 normvalue = (value - mean)/(std * SQRTOFTWO); 334 335 SCIPdebugMsg(scip, " Normalized value %g = ( %g - %g ) / (%g * 1.4142136)\n", normvalue, value, mean, std); 336 337 /* calculate the cumulative distribution function for normvalue. For negative normvalues, we negate 338 * the normvalue and use the oddness of the SCIPerf()-function; special treatment for values close to zero. 339 */ 340 if( SCIPisFeasZero(scip, normvalue) ) 341 return .5; 342 else if( normvalue > 0 ) 343 { 344 SCIP_Real erfresult; 345 346 erfresult = SCIPerf(normvalue); 347 return erfresult / 2.0 + 0.5; 348 } 349 else 350 { 351 SCIP_Real erfresult; 352 353 erfresult = SCIPerf(-normvalue); 354 355 return 0.5 - erfresult / 2.0; 356 } 357 } 358 359 /** calculates the probability of satisfying an LP-row under the assumption 360 * of uniformly distributed variable values. 361 * 362 * For inequalities, we use the cumulative distribution function of the standard normal 363 * distribution PHI(rhs - mu/sqrt(sigma2)) to calculate the probability 364 * for a right hand side row with mean activity mu and variance sigma2 to be satisfied. 365 * Similarly, 1 - PHI(lhs - mu/sqrt(sigma2)) is the probability to satisfy a left hand side row. 366 * For equations (lhs==rhs), we use the centeredness measure p = min(PHI(lhs'), 1-PHI(lhs'))/max(PHI(lhs'), 1 - PHI(lhs')), 367 * where lhs' = lhs - mu / sqrt(sigma2). 368 */ 369 SCIP_Real SCIProwCalcProbability( 370 SCIP* scip, /**< current scip */ 371 SCIP_ROW* row, /**< the row */ 372 SCIP_Real mu, /**< the mean value of the row distribution */ 373 SCIP_Real sigma2, /**< the variance of the row distribution */ 374 int rowinfinitiesdown, /**< the number of variables with infinite bounds to DECREASE activity */ 375 int rowinfinitiesup /**< the number of variables with infinite bounds to INCREASE activity */ 376 ) 377 { 378 SCIP_Real rowprobability; 379 SCIP_Real lhs; 380 SCIP_Real rhs; 381 SCIP_Real lhsprob; 382 SCIP_Real rhsprob; 383 384 lhs = SCIProwGetLhs(row); 385 rhs = SCIProwGetRhs(row); 386 387 lhsprob = 1.0; 388 rhsprob = 1.0; 389 390 /* use the cumulative distribution if the row contains no variable to repair every infeasibility */ 391 if( !SCIPisInfinity(scip, rhs) && rowinfinitiesdown == 0 ) 392 rhsprob = SCIPcalcCumulativeDistribution(scip, mu, sigma2, rhs); 393 394 /* use the cumulative distribution if the row contains no variable to repair every infeasibility 395 * otherwise the row can always be made feasible by increasing activity far enough 396 */ 397 if( !SCIPisInfinity(scip, -lhs) && rowinfinitiesup == 0 ) 398 lhsprob = 1.0 - SCIPcalcCumulativeDistribution(scip, mu, sigma2, lhs); 399 400 assert(SCIPisFeasLE(scip, lhsprob, 1.0) && SCIPisFeasGE(scip, lhsprob, 0.0)); 401 assert(SCIPisFeasLE(scip, rhsprob, 1.0) && SCIPisFeasGE(scip, rhsprob, 0.0)); 402 403 /* use centeredness measure for equations; for inequalities, the minimum of the two probabilities is the 404 * probability to satisfy the row */ 405 if( SCIPisFeasEQ(scip, lhs, rhs) ) 406 { 407 SCIP_Real minprobability; 408 SCIP_Real maxprobability; 409 410 minprobability = MIN(rhsprob, lhsprob); 411 maxprobability = MAX(lhsprob, rhsprob); 412 rowprobability = minprobability / maxprobability; 413 } 414 else 415 rowprobability = MIN(rhsprob, lhsprob); 416 417 SCIPdebug( SCIPprintRow(scip, row, NULL) ); 418 SCIPdebugMsg(scip, " Row %s, mean %g, sigma2 %g, LHS %g, RHS %g has probability %g to be satisfied\n", 419 SCIProwGetName(row), mu, sigma2, lhs, rhs, rowprobability); 420 421 assert(SCIPisFeasGE(scip, rowprobability, 0.0) && SCIPisFeasLE(scip, rowprobability, 1.0)); 422 423 return rowprobability; 424 } 425 426 /** calculates the initial mean and variance of the row activity normal distribution. 427 * 428 * The mean value \f$ \mu \f$ is given by \f$ \mu = \sum_i=1^n c_i * (lb_i +ub_i) / 2 \f$ where 429 * \f$n \f$ is the number of variables, and \f$ c_i, lb_i, ub_i \f$ are the variable coefficient and 430 * bounds, respectively. With the same notation, the variance \f$ \sigma^2 \f$ is given by 431 * \f$ \sigma^2 = \sum_i=1^n c_i^2 * \sigma^2_i \f$, with the variance being 432 * \f$ \sigma^2_i = ((ub_i - lb_i + 1)^2 - 1) / 12 \f$ for integer variables and 433 * \f$ \sigma^2_i = (ub_i - lb_i)^2 / 12 \f$ for continuous variables. 434 */ 435 static 436 void rowCalculateGauss( 437 SCIP* scip, /**< SCIP data structure */ 438 SCIP_BRANCHRULEDATA* branchruledata, /**< the branching rule data */ 439 SCIP_ROW* row, /**< the row for which the gaussian normal distribution has to be calculated */ 440 SCIP_Real* mu, /**< pointer to store the mean value of the gaussian normal distribution */ 441 SCIP_Real* sigma2, /**< pointer to store the variance value of the gaussian normal distribution */ 442 int* rowinfinitiesdown, /**< pointer to store the number of variables with infinite bounds to DECREASE activity */ 443 int* rowinfinitiesup /**< pointer to store the number of variables with infinite bounds to INCREASE activity */ 444 ) 445 { 446 SCIP_COL** rowcols; 447 SCIP_Real* rowvals; 448 int nrowvals; 449 int c; 450 451 assert(scip != NULL); 452 assert(row != NULL); 453 assert(mu != NULL); 454 assert(sigma2 != NULL); 455 assert(rowinfinitiesup != NULL); 456 assert(rowinfinitiesdown != NULL); 457 458 rowcols = SCIProwGetCols(row); 459 rowvals = SCIProwGetVals(row); 460 nrowvals = SCIProwGetNNonz(row); 461 462 assert(nrowvals == 0 || rowcols != NULL); 463 assert(nrowvals == 0 || rowvals != NULL); 464 465 *mu = SCIProwGetConstant(row); 466 *sigma2 = 0.0; 467 *rowinfinitiesdown = 0; 468 *rowinfinitiesup = 0; 469 470 /* loop over nonzero row coefficients and sum up the variable contributions to mu and sigma2 */ 471 for( c = 0; c < nrowvals; ++c ) 472 { 473 SCIP_VAR* colvar; 474 SCIP_Real colval; 475 SCIP_Real colvarlb; 476 SCIP_Real colvarub; 477 SCIP_Real squarecoeff; 478 SCIP_Real varvariance; 479 SCIP_Real varmean; 480 int varindex; 481 482 assert(rowcols[c] != NULL); 483 colvar = SCIPcolGetVar(rowcols[c]); 484 assert(colvar != NULL); 485 486 colval = rowvals[c]; 487 colvarlb = SCIPvarGetLbLocal(colvar); 488 colvarub = SCIPvarGetUbLocal(colvar); 489 490 varmean = 0.0; 491 varvariance = 0.0; 492 varindex = SCIPvarGetProbindex(colvar); 493 assert((branchruledata->currentlbs[varindex] == SCIP_INVALID) == (branchruledata->currentubs[varindex] == SCIP_INVALID)); /*lint !e777*/ 494 495 /* variable bounds need to be watched from now on */ 496 if( branchruledata->currentlbs[varindex] == SCIP_INVALID ) /*lint !e777*/ 497 branchruledataUpdateCurrentBounds(scip, branchruledata, colvar); 498 499 assert(!SCIPisInfinity(scip, colvarlb)); 500 assert(!SCIPisInfinity(scip, -colvarub)); 501 assert(SCIPisFeasLE(scip, colvarlb, colvarub)); 502 503 /* variables with infinite bounds are skipped for the calculation of the variance; they need to 504 * be accounted for by the counters for infinite row activity decrease and increase and they 505 * are used to shift the row activity mean in case they have one nonzero, but finite bound */ 506 if( SCIPisInfinity(scip, -colvarlb) || SCIPisInfinity(scip, colvarub) ) 507 { 508 if( SCIPisInfinity(scip, colvarub) ) 509 { 510 /* an infinite upper bound gives the row an infinite maximum activity or minimum activity, if the coefficient is 511 * positive or negative, resp. 512 */ 513 if( colval < 0.0 ) 514 ++(*rowinfinitiesdown); 515 else 516 ++(*rowinfinitiesup); 517 } 518 519 /* an infinite lower bound gives the row an infinite maximum activity or minimum activity, if the coefficient is 520 * negative or positive, resp. 521 */ 522 if( SCIPisInfinity(scip, -colvarlb) ) 523 { 524 if( colval > 0.0 ) 525 ++(*rowinfinitiesdown); 526 else 527 ++(*rowinfinitiesup); 528 } 529 } 530 SCIPvarCalcDistributionParameters(scip, colvarlb, colvarub, SCIPvarGetType(colvar), &varmean, &varvariance); 531 532 /* actual values are updated; the contribution of the variable to mu is the arithmetic mean of its bounds */ 533 *mu += colval * varmean; 534 535 /* the variance contribution of a variable is c^2 * (u - l)^2 / 12.0 for continuous and c^2 * ((u - l + 1)^2 - 1) / 12.0 for integer */ 536 squarecoeff = SQUARED(colval); 537 *sigma2 += squarecoeff * varvariance; 538 539 assert(!SCIPisFeasNegative(scip, *sigma2)); 540 } 541 542 SCIPdebug( SCIPprintRow(scip, row, NULL) ); 543 SCIPdebugMsg(scip, " Row %s has a mean value of %g at a sigma2 of %g \n", SCIProwGetName(row), *mu, *sigma2); 544 } 545 546 /** update the up- and downscore of a single variable after calculating the impact of branching on a 547 * particular row, depending on the chosen score parameter 548 */ 549 SCIP_RETCODE SCIPupdateDistributionScore( 550 SCIP* scip, /**< current SCIP pointer */ 551 SCIP_Real currentprob, /**< the current probability */ 552 SCIP_Real newprobup, /**< the new probability if branched upwards */ 553 SCIP_Real newprobdown, /**< the new probability if branched downwards */ 554 SCIP_Real* upscore, /**< pointer to store the new score for branching up */ 555 SCIP_Real* downscore, /**< pointer to store the new score for branching down */ 556 char scoreparam /**< parameter to determine the way the score is calculated */ 557 ) 558 { 559 assert(scip != NULL); 560 assert(SCIPisFeasGE(scip, currentprob, 0.0) && SCIPisFeasLE(scip, currentprob, 1.0)); 561 assert(SCIPisFeasGE(scip, newprobup, 0.0) && SCIPisFeasLE(scip, newprobup, 1.0)); 562 assert(SCIPisFeasGE(scip, newprobdown, 0.0) && SCIPisFeasLE(scip, newprobdown, 1.0)); 563 assert(upscore != NULL); 564 assert(downscore != NULL); 565 566 /* update up and downscore depending on score parameter */ 567 switch( scoreparam ) 568 { 569 case 'l' : 570 /* 'l'owest cumulative probability */ 571 if( SCIPisGT(scip, 1.0 - newprobup, *upscore) ) 572 *upscore = 1.0 - newprobup; 573 if( SCIPisGT(scip, 1.0 - newprobdown, *downscore) ) 574 *downscore = 1.0 - newprobdown; 575 break; 576 577 case 'd' : 578 /* biggest 'd'ifference currentprob - newprob */ 579 if( SCIPisGT(scip, currentprob - newprobup, *upscore) ) 580 *upscore = currentprob - newprobup; 581 if( SCIPisGT(scip, currentprob - newprobdown, *downscore) ) 582 *downscore = currentprob - newprobdown; 583 break; 584 585 case 'h' : 586 /* 'h'ighest cumulative probability */ 587 if( SCIPisGT(scip, newprobup, *upscore) ) 588 *upscore = newprobup; 589 if( SCIPisGT(scip, newprobdown, *downscore) ) 590 *downscore = newprobdown; 591 break; 592 593 case 'v' : 594 /* 'v'otes lowest cumulative probability */ 595 if( SCIPisLT(scip, newprobup, newprobdown) ) 596 *upscore += 1.0; 597 else if( SCIPisGT(scip, newprobup, newprobdown) ) 598 *downscore += 1.0; 599 break; 600 601 case 'w' : 602 /* votes highest cumulative probability */ 603 if( SCIPisGT(scip, newprobup, newprobdown) ) 604 *upscore += 1.0; 605 else if( SCIPisLT(scip, newprobup, newprobdown) ) 606 *downscore += 1.0; 607 break; 608 609 default : 610 SCIPerrorMessage(" ERROR! No branching scheme selected! Exiting method.\n"); 611 return SCIP_INVALIDCALL; 612 } 613 614 return SCIP_OKAY; 615 } 616 617 /** calculate the branching score of a variable, depending on the chosen score parameter */ 618 static 619 SCIP_RETCODE calcBranchScore( 620 SCIP* scip, /**< current SCIP */ 621 SCIP_BRANCHRULEDATA* branchruledata, /**< branch rule data */ 622 SCIP_VAR* var, /**< candidate variable */ 623 SCIP_Real lpsolval, /**< current fractional LP-relaxation solution value */ 624 SCIP_Real* upscore, /**< pointer to store the variable score when branching on it in upward direction */ 625 SCIP_Real* downscore, /**< pointer to store the variable score when branching on it in downward direction */ 626 char scoreparam /**< the score parameter of this branching rule */ 627 ) 628 { 629 SCIP_COL* varcol; 630 SCIP_ROW** colrows; 631 SCIP_Real* rowvals; 632 SCIP_Real varlb; 633 SCIP_Real varub; 634 SCIP_Real squaredbounddiff; /* current squared difference of variable bounds (ub - lb)^2 */ 635 SCIP_Real newub; /* new upper bound if branching downwards */ 636 SCIP_Real newlb; /* new lower bound if branching upwards */ 637 SCIP_Real squaredbounddiffup; /* squared difference after branching upwards (ub - lb')^2 */ 638 SCIP_Real squaredbounddiffdown; /* squared difference after branching downwards (ub' - lb)^2 */ 639 SCIP_Real currentmean; /* current mean value of variable uniform distribution */ 640 SCIP_Real meanup; /* mean value of variable uniform distribution after branching up */ 641 SCIP_Real meandown; /* mean value of variable uniform distribution after branching down*/ 642 SCIP_VARTYPE vartype; 643 int ncolrows; 644 int i; 645 646 SCIP_Bool onlyactiverows; /* should only rows which are active at the current node be considered? */ 647 648 assert(scip != NULL); 649 assert(var != NULL); 650 assert(upscore != NULL); 651 assert(downscore != NULL); 652 assert(!SCIPisIntegral(scip, lpsolval)); 653 assert(SCIPvarGetStatus(var) == SCIP_VARSTATUS_COLUMN); 654 655 varcol = SCIPvarGetCol(var); 656 assert(varcol != NULL); 657 658 colrows = SCIPcolGetRows(varcol); 659 rowvals = SCIPcolGetVals(varcol); 660 ncolrows = SCIPcolGetNNonz(varcol); 661 varlb = SCIPvarGetLbLocal(var); 662 varub = SCIPvarGetUbLocal(var); 663 assert(SCIPisFeasLT(scip, varlb, varub)); 664 vartype = SCIPvarGetType(var); 665 666 /* calculate mean and variance of variable uniform distribution before and after branching */ 667 currentmean = 0.0; 668 squaredbounddiff = 0.0; 669 SCIPvarCalcDistributionParameters(scip, varlb, varub, vartype, ¤tmean, &squaredbounddiff); 670 671 newlb = SCIPfeasCeil(scip, lpsolval); 672 newub = SCIPfeasFloor(scip, lpsolval); 673 674 /* calculate the variable's uniform distribution after branching up and down, respectively. */ 675 squaredbounddiffup = 0.0; 676 meanup = 0.0; 677 SCIPvarCalcDistributionParameters(scip, newlb, varub, vartype, &meanup, &squaredbounddiffup); 678 679 /* calculate the distribution mean and variance for a variable with finite lower bound */ 680 squaredbounddiffdown = 0.0; 681 meandown = 0.0; 682 SCIPvarCalcDistributionParameters(scip, varlb, newub, vartype, &meandown, &squaredbounddiffdown); 683 684 /* initialize the variable's up and down score */ 685 *upscore = 0.0; 686 *downscore = 0.0; 687 688 onlyactiverows = branchruledata->onlyactiverows; 689 690 /* loop over the variable rows and calculate the up and down score */ 691 for( i = 0; i < ncolrows; ++i ) 692 { 693 SCIP_ROW* row; 694 SCIP_Real changedrowmean; 695 SCIP_Real rowmean; 696 SCIP_Real rowvariance; 697 SCIP_Real changedrowvariance; 698 SCIP_Real currentrowprob; 699 SCIP_Real newrowprobup; 700 SCIP_Real newrowprobdown; 701 SCIP_Real squaredcoeff; 702 SCIP_Real rowval; 703 int rowinfinitiesdown; 704 int rowinfinitiesup; 705 int rowpos; 706 707 row = colrows[i]; 708 rowval = rowvals[i]; 709 assert(row != NULL); 710 711 /* we access the rows by their index */ 712 rowpos = SCIProwGetIndex(row); 713 714 /* skip non-active rows if the user parameter was set this way */ 715 if( onlyactiverows && SCIPisSumPositive(scip, SCIPgetRowLPFeasibility(scip, row)) ) 716 continue; 717 718 /* call method to ensure sufficient data capacity */ 719 SCIP_CALL( branchruledataEnsureArraySize(scip, branchruledata, rowpos) ); 720 721 /* calculate row activity distribution if this is the first candidate to appear in this row */ 722 if( branchruledata->rowmeans[rowpos] == SCIP_INVALID ) /*lint !e777*/ 723 { 724 rowCalculateGauss(scip, branchruledata, row, &branchruledata->rowmeans[rowpos], &branchruledata->rowvariances[rowpos], 725 &branchruledata->rowinfinitiesdown[rowpos], &branchruledata->rowinfinitiesup[rowpos]); 726 } 727 728 /* retrieve the row distribution parameters from the branch rule data */ 729 rowmean = branchruledata->rowmeans[rowpos]; 730 rowvariance = branchruledata->rowvariances[rowpos]; 731 rowinfinitiesdown = branchruledata->rowinfinitiesdown[rowpos]; 732 rowinfinitiesup = branchruledata->rowinfinitiesdown[rowpos]; 733 assert(!SCIPisNegative(scip, rowvariance)); 734 735 currentrowprob = SCIProwCalcProbability(scip, row, rowmean, rowvariance, 736 rowinfinitiesdown, rowinfinitiesup); 737 738 /* get variable's current expected contribution to row activity */ 739 squaredcoeff = SQUARED(rowval); 740 741 /* first, get the probability change for the row if the variable is branched on upwards. The probability 742 * can only be affected if the variable upper bound is finite 743 */ 744 if( !SCIPisInfinity(scip, varub) ) 745 { 746 int rowinftiesdownafterbranch; 747 int rowinftiesupafterbranch; 748 749 /* calculate how branching would affect the row parameters */ 750 changedrowmean = rowmean + rowval * (meanup - currentmean); 751 changedrowvariance = rowvariance + squaredcoeff * (squaredbounddiffup - squaredbounddiff); 752 changedrowvariance = MAX(0.0, changedrowvariance); 753 754 rowinftiesdownafterbranch = rowinfinitiesdown; 755 rowinftiesupafterbranch = rowinfinitiesup; 756 757 /* account for changes of the row's infinite bound contributions */ 758 if( SCIPisInfinity(scip, -varlb) && rowval < 0.0 ) 759 rowinftiesupafterbranch--; 760 if( SCIPisInfinity(scip, -varlb) && rowval > 0.0 ) 761 rowinftiesdownafterbranch--; 762 763 assert(rowinftiesupafterbranch >= 0); 764 assert(rowinftiesdownafterbranch >= 0); 765 newrowprobup = SCIProwCalcProbability(scip, row, changedrowmean, changedrowvariance, rowinftiesdownafterbranch, 766 rowinftiesupafterbranch); 767 } 768 else 769 newrowprobup = currentrowprob; 770 771 /* do the same for the other branching direction */ 772 if( !SCIPisInfinity(scip, varlb) ) 773 { 774 int rowinftiesdownafterbranch; 775 int rowinftiesupafterbranch; 776 777 changedrowmean = rowmean + rowval * (meandown - currentmean); 778 changedrowvariance = rowvariance + squaredcoeff * (squaredbounddiffdown - squaredbounddiff); 779 changedrowvariance = MAX(0.0, changedrowvariance); 780 781 rowinftiesdownafterbranch = rowinfinitiesdown; 782 rowinftiesupafterbranch = rowinfinitiesup; 783 784 /* account for changes of the row's infinite bound contributions */ 785 if( SCIPisInfinity(scip, varub) && rowval > 0.0 ) 786 rowinftiesupafterbranch -= 1; 787 if( SCIPisInfinity(scip, varub) && rowval < 0.0 ) 788 rowinftiesdownafterbranch -= 1; 789 790 assert(rowinftiesdownafterbranch >= 0); 791 assert(rowinftiesupafterbranch >= 0); 792 newrowprobdown = SCIProwCalcProbability(scip, row, changedrowmean, changedrowvariance, rowinftiesdownafterbranch, 793 rowinftiesupafterbranch); 794 } 795 else 796 newrowprobdown = currentrowprob; 797 798 /* update the up and down score depending on the chosen scoring parameter */ 799 SCIP_CALL( SCIPupdateDistributionScore(scip, currentrowprob, newrowprobup, newrowprobdown, upscore, downscore, scoreparam) ); 800 801 SCIPdebugMsg(scip, " Variable %s changes probability of row %s from %g to %g (branch up) or %g;\n", 802 SCIPvarGetName(var), SCIProwGetName(row), currentrowprob, newrowprobup, newrowprobdown); 803 SCIPdebugMsg(scip, " --> new variable score: %g (for branching up), %g (for branching down)\n", 804 *upscore, *downscore); 805 } 806 807 return SCIP_OKAY; 808 } 809 810 /** free branchrule data */ 811 static 812 void branchruledataFreeArrays( 813 SCIP* scip, /**< SCIP data structure */ 814 SCIP_BRANCHRULEDATA* branchruledata /**< branching rule data */ 815 ) 816 { 817 assert(branchruledata->memsize == 0 || branchruledata->rowmeans != NULL); 818 assert(branchruledata->memsize >= 0); 819 820 if( branchruledata->memsize > 0 ) 821 { 822 SCIPfreeBlockMemoryArray(scip, &branchruledata->rowmeans, branchruledata->memsize); 823 SCIPfreeBlockMemoryArray(scip, &branchruledata->rowvariances, branchruledata->memsize); 824 SCIPfreeBlockMemoryArray(scip, &branchruledata->rowinfinitiesup, branchruledata->memsize); 825 SCIPfreeBlockMemoryArray(scip, &branchruledata->rowinfinitiesdown, branchruledata->memsize); 826 827 SCIPfreeBlockMemoryArray(scip, &branchruledata->varfilterposs, branchruledata->varpossmemsize); 828 SCIPfreeBlockMemoryArray(scip, &branchruledata->varposs, branchruledata->varpossmemsize); 829 SCIPfreeBlockMemoryArray(scip, &branchruledata->updatedvars, branchruledata->varpossmemsize); 830 SCIPfreeBlockMemoryArray(scip, &branchruledata->currentubs, branchruledata->varpossmemsize); 831 SCIPfreeBlockMemoryArray(scip, &branchruledata->currentlbs, branchruledata->varpossmemsize); 832 833 branchruledata->memsize = 0; 834 } 835 } 836 837 /** add variable to the bound change event queue; skipped if variable is already in there, or if variable has 838 * no row currently watched 839 */ 840 static 841 void branchruledataAddBoundChangeVar( 842 SCIP_BRANCHRULEDATA* branchruledata, /**< branchrule data */ 843 SCIP_VAR* var /**< the variable whose bound changes need to be processed */ 844 ) 845 { 846 int varindex; 847 int varpos; 848 849 assert(var != NULL); 850 851 varindex = SCIPvarGetProbindex(var); 852 assert(-1 <= varindex && varindex < branchruledata->varpossmemsize); 853 854 /* if variable is not active, it should not be watched */ 855 if( varindex == -1 ) 856 return; 857 varpos = branchruledata->varposs[varindex]; 858 assert(varpos < branchruledata->nupdatedvars); 859 860 /* nothing to do if variable is already in the queue */ 861 if( varpos >= 0 ) 862 { 863 assert(branchruledata->updatedvars[varpos] == var); 864 865 return; 866 } 867 868 /* if none of the variables rows was calculated yet, variable needs not to be watched */ 869 assert((branchruledata->currentlbs[varindex] == SCIP_INVALID) == (branchruledata->currentubs[varindex] == SCIP_INVALID)); /*lint !e777*/ 870 if( branchruledata->currentlbs[varindex] == SCIP_INVALID ) /*lint !e777*/ 871 return; 872 873 /* add the variable to the branch rule data of variables to process updates for */ 874 assert(branchruledata->varpossmemsize > branchruledata->nupdatedvars); 875 varpos = branchruledata->nupdatedvars; 876 branchruledata->updatedvars[varpos] = var; 877 branchruledata->varposs[varindex] = varpos; 878 ++branchruledata->nupdatedvars; 879 } 880 881 /** returns the next unprocessed variable (last in, first out) with pending bound changes, or NULL */ 882 static 883 SCIP_VAR* branchruledataPopBoundChangeVar( 884 SCIP_BRANCHRULEDATA* branchruledata /**< branchrule data */ 885 ) 886 { 887 SCIP_VAR* var; 888 int varpos; 889 int varindex; 890 891 assert(branchruledata->nupdatedvars >= 0); 892 893 /* return if no variable is currently pending */ 894 if( branchruledata->nupdatedvars == 0 ) 895 return NULL; 896 897 varpos = branchruledata->nupdatedvars - 1; 898 var = branchruledata->updatedvars[varpos]; 899 assert(var != NULL); 900 varindex = SCIPvarGetProbindex(var); 901 assert(0 <= varindex && varindex < branchruledata->varpossmemsize); 902 assert(varpos == branchruledata->varposs[varindex]); 903 904 branchruledata->varposs[varindex] = -1; 905 branchruledata->nupdatedvars--; 906 907 return var; 908 } 909 910 /** process a variable from the queue of changed variables */ 911 static 912 SCIP_RETCODE varProcessBoundChanges( 913 SCIP* scip, /**< SCIP data structure */ 914 SCIP_BRANCHRULEDATA* branchruledata, /**< branchrule data */ 915 SCIP_VAR* var /**< the variable whose bound changes need to be processed */ 916 ) 917 { 918 SCIP_ROW** colrows; 919 SCIP_COL* varcol; 920 SCIP_Real* colvals; 921 SCIP_Real oldmean; 922 SCIP_Real newmean; 923 SCIP_Real oldvariance; 924 SCIP_Real newvariance; 925 SCIP_Real oldlb; 926 SCIP_Real newlb; 927 SCIP_Real oldub; 928 SCIP_Real newub; 929 SCIP_VARTYPE vartype; 930 int ncolrows; 931 int r; 932 int varindex; 933 934 /* skip event execution if SCIP is in Probing mode because these bound changes will be undone anyway before branching 935 * rule is called again 936 */ 937 assert(!SCIPinProbing(scip)); 938 939 assert(var != NULL); 940 varcol = SCIPvarGetCol(var); 941 assert(varcol != NULL); 942 colrows = SCIPcolGetRows(varcol); 943 colvals = SCIPcolGetVals(varcol); 944 ncolrows = SCIPcolGetNNonz(varcol); 945 946 varindex = SCIPvarGetProbindex(var); 947 948 oldlb = branchruledata->currentlbs[varindex]; 949 oldub = branchruledata->currentubs[varindex]; 950 951 /* skip update if the variable has never been subject of previously calculated row activities */ 952 assert((oldlb == SCIP_INVALID) == (oldub == SCIP_INVALID)); /*lint !e777*/ 953 if( oldlb == SCIP_INVALID ) /*lint !e777*/ 954 return SCIP_OKAY; 955 956 newlb = SCIPvarGetLbLocal(var); 957 newub = SCIPvarGetUbLocal(var); 958 959 /* skip update if the bound change events have cancelled out */ 960 if( SCIPisFeasEQ(scip, oldlb, newlb) && SCIPisFeasEQ(scip, oldub, newub) ) 961 return SCIP_OKAY; 962 963 /* calculate old and new variable distribution mean and variance */ 964 oldvariance = 0.0; 965 newvariance = 0.0; 966 oldmean = 0.0; 967 newmean = 0.0; 968 vartype = SCIPvarGetType(var); 969 SCIPvarCalcDistributionParameters(scip, oldlb, oldub, vartype, &oldmean, &oldvariance); 970 SCIPvarCalcDistributionParameters(scip, newlb, newub, vartype, &newmean, &newvariance); 971 972 /* loop over all rows of this variable and update activity distribution */ 973 for( r = 0; r < ncolrows; ++r ) 974 { 975 int rowpos; 976 977 assert(colrows[r] != NULL); 978 rowpos = SCIProwGetIndex(colrows[r]); 979 assert(rowpos >= 0); 980 981 SCIP_CALL( branchruledataEnsureArraySize(scip, branchruledata, rowpos) ); 982 983 /* only consider rows for which activity distribution was already calculated */ 984 if( branchruledata->rowmeans[rowpos] != SCIP_INVALID ) /*lint !e777*/ 985 { 986 SCIP_Real coeff; 987 SCIP_Real coeffsquared; 988 /*lint -e777*/ 989 assert(branchruledata->rowvariances[rowpos] != SCIP_INVALID 990 && SCIPisFeasGE(scip, branchruledata->rowvariances[rowpos], 0.0)); 991 992 coeff = colvals[r]; 993 coeffsquared = SQUARED(coeff); 994 995 /* update variable contribution to row activity distribution */ 996 branchruledata->rowmeans[rowpos] += coeff * (newmean - oldmean); 997 branchruledata->rowvariances[rowpos] += coeffsquared * (newvariance - oldvariance); 998 branchruledata->rowvariances[rowpos] = MAX(0.0, branchruledata->rowvariances[rowpos]); 999 1000 /* account for changes of the infinite contributions to row activities */ 1001 if( coeff > 0.0 ) 1002 { 1003 /* if the coefficient is positive, upper bounds affect activity up */ 1004 if( SCIPisInfinity(scip, newub) && !SCIPisInfinity(scip, oldub) ) 1005 ++branchruledata->rowinfinitiesup[rowpos]; 1006 else if( !SCIPisInfinity(scip, newub) && SCIPisInfinity(scip, oldub) ) 1007 --branchruledata->rowinfinitiesup[rowpos]; 1008 1009 if( SCIPisInfinity(scip, newlb) && !SCIPisInfinity(scip, oldlb) ) 1010 ++branchruledata->rowinfinitiesdown[rowpos]; 1011 else if( !SCIPisInfinity(scip, newlb) && SCIPisInfinity(scip, oldlb) ) 1012 --branchruledata->rowinfinitiesdown[rowpos]; 1013 } 1014 else if( coeff < 0.0 ) 1015 { 1016 if( SCIPisInfinity(scip, newub) && !SCIPisInfinity(scip, oldub) ) 1017 ++branchruledata->rowinfinitiesdown[rowpos]; 1018 else if( !SCIPisInfinity(scip, newub) && SCIPisInfinity(scip, oldub) ) 1019 --branchruledata->rowinfinitiesdown[rowpos]; 1020 1021 if( SCIPisInfinity(scip, newlb) && !SCIPisInfinity(scip, oldlb) ) 1022 ++branchruledata->rowinfinitiesup[rowpos]; 1023 else if( !SCIPisInfinity(scip, newlb) && SCIPisInfinity(scip, oldlb) ) 1024 --branchruledata->rowinfinitiesup[rowpos]; 1025 } 1026 assert(branchruledata->rowinfinitiesdown[rowpos] >= 0); 1027 assert(branchruledata->rowinfinitiesup[rowpos] >= 0); 1028 } 1029 } 1030 1031 /* store the new local bounds in the data */ 1032 branchruledataUpdateCurrentBounds(scip, branchruledata, var); 1033 1034 return SCIP_OKAY; 1035 } 1036 1037 /** destructor of event handler to free user data (called when SCIP is exiting) */ 1038 static 1039 SCIP_DECL_EVENTFREE(eventFreeDistribution) 1040 { 1041 SCIP_EVENTHDLRDATA* eventhdlrdata; 1042 1043 eventhdlrdata = SCIPeventhdlrGetData(eventhdlr); 1044 assert(eventhdlrdata != NULL); 1045 1046 SCIPfreeBlockMemory(scip, &eventhdlrdata); 1047 SCIPeventhdlrSetData(eventhdlr, NULL); 1048 1049 return SCIP_OKAY; 1050 } 1051 1052 /* 1053 * Callback methods of branching rule 1054 */ 1055 1056 /** copy method for branchrule plugins (called when SCIP copies plugins) */ 1057 static 1058 SCIP_DECL_BRANCHCOPY(branchCopyDistribution) 1059 { /*lint --e{715}*/ 1060 assert(scip != NULL); 1061 1062 SCIP_CALL( SCIPincludeBranchruleDistribution(scip) ); 1063 1064 return SCIP_OKAY; 1065 } 1066 1067 /** solving process deinitialization method of branching rule (called before branch and bound process data is freed) */ 1068 static 1069 SCIP_DECL_BRANCHEXITSOL(branchExitsolDistribution) 1070 { 1071 SCIP_BRANCHRULEDATA* branchruledata; 1072 1073 assert(branchrule != NULL); 1074 assert(strcmp(SCIPbranchruleGetName(branchrule), BRANCHRULE_NAME) == 0); 1075 1076 branchruledata = SCIPbranchruleGetData(branchrule); 1077 assert(branchruledata != NULL); 1078 1079 /* drop variable events at the end of branch and bound process (cannot be used after restarts, anyway) */ 1080 if( branchruledata->varfilterposs != NULL) 1081 { 1082 SCIP_VAR** vars; 1083 int nvars; 1084 int v; 1085 1086 vars = SCIPgetVars(scip); 1087 nvars = SCIPgetNVars(scip); 1088 1089 assert(nvars > 0); 1090 for( v = 0; v < nvars; ++v ) 1091 { 1092 SCIP_CALL( SCIPdropVarEvent(scip, vars[v], EVENT_DISTRIBUTION, branchruledata->eventhdlr, NULL, branchruledata->varfilterposs[v]) ); 1093 } 1094 } 1095 1096 /* free row arrays when branch-and-bound data is freed */ 1097 branchruledataFreeArrays(scip, branchruledata); 1098 1099 return SCIP_OKAY; 1100 } 1101 1102 /** destructor of branching rule to free user data (called when SCIP is exiting) */ 1103 static 1104 SCIP_DECL_BRANCHFREE(branchFreeDistribution) 1105 { 1106 SCIP_BRANCHRULEDATA* branchruledata; 1107 1108 assert(branchrule != NULL); 1109 assert(strcmp(SCIPbranchruleGetName(branchrule), BRANCHRULE_NAME) == 0); 1110 1111 branchruledata = SCIPbranchruleGetData(branchrule); 1112 assert(branchruledata != NULL); 1113 1114 /* free internal arrays first */ 1115 branchruledataFreeArrays(scip, branchruledata); 1116 SCIPfreeBlockMemory(scip, &branchruledata); 1117 SCIPbranchruleSetData(branchrule, NULL); 1118 1119 return SCIP_OKAY; 1120 } 1121 1122 /** branching execution method for fractional LP solutions */ 1123 static 1124 SCIP_DECL_BRANCHEXECLP(branchExeclpDistribution) 1125 { /*lint --e{715}*/ 1126 SCIP_BRANCHRULEDATA* branchruledata; 1127 SCIP_VAR** lpcands; 1128 SCIP_VAR* bestcand; 1129 SCIP_NODE* downchild; 1130 SCIP_NODE* upchild; 1131 1132 SCIP_Real* lpcandssol; 1133 1134 SCIP_Real bestscore; 1135 SCIP_BRANCHDIR bestbranchdir; 1136 int nlpcands; 1137 int c; 1138 1139 assert(branchrule != NULL); 1140 assert(strcmp(SCIPbranchruleGetName(branchrule), BRANCHRULE_NAME) == 0); 1141 assert(scip != NULL); 1142 assert(result != NULL); 1143 1144 *result = SCIP_DIDNOTRUN; 1145 1146 SCIP_CALL( SCIPgetLPBranchCands(scip, &lpcands, &lpcandssol, NULL, NULL, &nlpcands, NULL) ); 1147 1148 if( nlpcands == 0 ) 1149 return SCIP_OKAY; 1150 1151 if( SCIPgetNActivePricers(scip) > 0 ) 1152 return SCIP_OKAY; 1153 1154 branchruledata = SCIPbranchruleGetData(branchrule); 1155 1156 /* if branching rule data arrays were not initialized before (usually the first call of the branching execution), 1157 * allocate arrays with an initial capacity of the number of LP rows */ 1158 if( branchruledata->memsize == 0 ) 1159 { 1160 int nlprows; 1161 1162 /* get LP rows data */ 1163 nlprows = SCIPgetNLPRows(scip); 1164 1165 /* initialize arrays only if there are actual LP rows to operate on */ 1166 if( nlprows > 0 ) 1167 { 1168 SCIP_CALL( branchruledataEnsureArraySize(scip, branchruledata, nlprows) ); 1169 } 1170 else /* if there are no LP rows, branching rule cannot be used */ 1171 return SCIP_OKAY; 1172 } 1173 1174 /* process pending bound change events */ 1175 while( branchruledata->nupdatedvars > 0 ) 1176 { 1177 SCIP_VAR* nextvar; 1178 1179 /* pop the next variable from the queue and process its bound changes */ 1180 nextvar = branchruledataPopBoundChangeVar(branchruledata); 1181 assert(nextvar != NULL); 1182 SCIP_CALL( varProcessBoundChanges(scip, branchruledata, nextvar) ); 1183 } 1184 1185 bestscore = -1; 1186 bestbranchdir = SCIP_BRANCHDIR_AUTO; 1187 bestcand = NULL; 1188 1189 /* invalidate the current row distribution data which is initialized on the fly when looping over the candidates */ 1190 1191 /* loop over candidate variables and calculate their score in changing the cumulative 1192 * probability of fulfilling each of their constraints */ 1193 for( c = 0; c < nlpcands; ++c ) 1194 { 1195 SCIP_Real upscore; 1196 SCIP_Real downscore; 1197 SCIP_VAR* lpcand; 1198 int varindex; 1199 1200 lpcand = lpcands[c]; 1201 assert(lpcand != NULL); 1202 1203 varindex = SCIPvarGetProbindex(lpcand); 1204 1205 /* in debug mode, ensure that all bound process events which occurred in the mean time have been captured 1206 * by the branching rule event system 1207 */ 1208 assert(SCIPisFeasLE(scip, SCIPvarGetLbLocal(lpcand), SCIPvarGetUbLocal(lpcand))); 1209 assert(0 <= varindex && varindex < branchruledata->varpossmemsize); 1210 1211 /*lint -e777*/ 1212 assert((branchruledata->currentlbs[varindex] == SCIP_INVALID) == (branchruledata->currentubs[varindex] == SCIP_INVALID)); 1213 assert((branchruledata->currentlbs[varindex] == SCIP_INVALID) 1214 || SCIPisFeasEQ(scip, SCIPvarGetLbLocal(lpcand), branchruledata->currentlbs[varindex])); 1215 assert((branchruledata->currentubs[varindex] == SCIP_INVALID) 1216 || SCIPisFeasEQ(scip, SCIPvarGetUbLocal(lpcand), branchruledata->currentubs[varindex])); 1217 1218 /* if the branching rule has not captured the variable bounds yet, this can be done now */ 1219 if( branchruledata->currentlbs[varindex] == SCIP_INVALID ) /*lint !e777*/ 1220 { 1221 branchruledataUpdateCurrentBounds(scip, branchruledata, lpcand); 1222 } 1223 1224 upscore = 0.0; 1225 downscore = 0.0; 1226 1227 /* loop over candidate rows and determine the candidate up- and down- branching score w.r.t. the score parameter */ 1228 SCIP_CALL( calcBranchScore(scip, branchruledata, lpcand, lpcandssol[c], 1229 &upscore, &downscore, branchruledata->scoreparam) ); 1230 1231 /* if weighted scoring is enabled, use the branching score method of SCIP to weigh up and down score */ 1232 if( branchruledata->usescipscore ) 1233 { 1234 SCIP_Real score; 1235 1236 score = SCIPgetBranchScore(scip, lpcand, downscore, upscore); 1237 1238 /* select the candidate with the highest branching score */ 1239 if( score > bestscore ) 1240 { 1241 bestscore = score; 1242 bestcand = lpcand; 1243 /* prioritize branching direction with the higher score */ 1244 if( upscore > downscore ) 1245 bestbranchdir = SCIP_BRANCHDIR_UPWARDS; 1246 else 1247 bestbranchdir = SCIP_BRANCHDIR_DOWNWARDS; 1248 } 1249 } 1250 else 1251 { 1252 /* no weighted score; keep candidate which has the single highest score in one direction */ 1253 if( upscore > bestscore && upscore > downscore ) 1254 { 1255 bestscore = upscore; 1256 bestbranchdir = SCIP_BRANCHDIR_UPWARDS; 1257 bestcand = lpcand; 1258 } 1259 else if( downscore > bestscore ) 1260 { 1261 bestscore = downscore; 1262 bestbranchdir = SCIP_BRANCHDIR_DOWNWARDS; 1263 bestcand = lpcand; 1264 } 1265 } 1266 1267 SCIPdebugMsg(scip, " Candidate %s has score down %g and up %g \n", SCIPvarGetName(lpcand), downscore, upscore); 1268 SCIPdebugMsg(scip, " Best candidate: %s, score %g, direction %d\n", SCIPvarGetName(bestcand), bestscore, bestbranchdir); 1269 } 1270 assert(!SCIPisFeasIntegral(scip, SCIPvarGetSol(bestcand, TRUE))); 1271 assert(bestbranchdir == SCIP_BRANCHDIR_DOWNWARDS || bestbranchdir == SCIP_BRANCHDIR_UPWARDS); 1272 assert(bestcand != NULL); 1273 1274 SCIPdebugMsg(scip, " Branching on variable %s with bounds [%g, %g] and solution value <%g>\n", SCIPvarGetName(bestcand), 1275 SCIPvarGetLbLocal(bestcand), SCIPvarGetUbLocal(bestcand), SCIPvarGetLPSol(bestcand)); 1276 1277 /* branch on the best candidate variable */ 1278 SCIP_CALL( SCIPbranchVar(scip, bestcand, &downchild, NULL, &upchild) ); 1279 1280 assert(downchild != NULL); 1281 assert(upchild != NULL); 1282 1283 if( bestbranchdir == SCIP_BRANCHDIR_UPWARDS ) 1284 { 1285 SCIP_CALL( SCIPchgChildPrio(scip, upchild, DEFAULT_PRIORITY) ); 1286 SCIPdebugMsg(scip, " Changing node priority of up-child.\n"); 1287 } 1288 else 1289 { 1290 assert(bestbranchdir == SCIP_BRANCHDIR_DOWNWARDS); 1291 SCIP_CALL( SCIPchgChildPrio(scip, downchild, DEFAULT_PRIORITY) ); 1292 SCIPdebugMsg(scip, " Changing node priority of down-child.\n"); 1293 } 1294 1295 *result = SCIP_BRANCHED; 1296 1297 return SCIP_OKAY; 1298 } 1299 1300 /** event execution method of distribution branching which handles bound change events of variables */ 1301 static 1302 SCIP_DECL_EVENTEXEC(eventExecDistribution) 1303 { /*lint --e{715}*/ 1304 SCIP_BRANCHRULEDATA* branchruledata; 1305 SCIP_EVENTHDLRDATA* eventhdlrdata; 1306 SCIP_VAR* var; 1307 1308 assert(eventhdlr != NULL); 1309 eventhdlrdata = SCIPeventhdlrGetData(eventhdlr); 1310 assert(eventhdlrdata != NULL); 1311 1312 branchruledata = eventhdlrdata->branchruledata; 1313 var = SCIPeventGetVar(event); 1314 1315 /* add the variable to the queue of unprocessed variables; method itself ensures that every variable is added at most once */ 1316 branchruledataAddBoundChangeVar(branchruledata, var); 1317 1318 return SCIP_OKAY; 1319 } 1320 1321 1322 /* 1323 * branching rule specific interface methods 1324 */ 1325 1326 /** creates the distribution branching rule and includes it in SCIP */ 1327 SCIP_RETCODE SCIPincludeBranchruleDistribution( 1328 SCIP* scip /**< SCIP data structure */ 1329 ) 1330 { 1331 SCIP_BRANCHRULE* branchrule = NULL; 1332 SCIP_BRANCHRULEDATA* branchruledata; 1333 SCIP_EVENTHDLRDATA* eventhdlrdata; 1334 1335 /* create distribution branching rule data */ 1336 branchruledata = NULL; 1337 SCIP_CALL( SCIPallocBlockMemory(scip, &branchruledata) ); 1338 1339 branchruledata->memsize = 0; 1340 branchruledata->rowmeans = NULL; 1341 branchruledata->rowvariances = NULL; 1342 branchruledata->rowinfinitiesdown = NULL; 1343 branchruledata->rowinfinitiesup = NULL; 1344 branchruledata->varfilterposs = NULL; 1345 branchruledata->currentlbs = NULL; 1346 branchruledata->currentubs = NULL; 1347 1348 /* create event handler first to finish branch rule data */ 1349 eventhdlrdata = NULL; 1350 SCIP_CALL( SCIPallocBlockMemory(scip, &eventhdlrdata) ); 1351 eventhdlrdata->branchruledata = branchruledata; 1352 1353 branchruledata->eventhdlr = NULL; 1354 SCIP_CALL( SCIPincludeEventhdlrBasic(scip, &branchruledata->eventhdlr, EVENTHDLR_NAME, 1355 "event handler for dynamic acitivity distribution updating", 1356 eventExecDistribution, eventhdlrdata) ); 1357 assert( branchruledata->eventhdlr != NULL); 1358 SCIP_CALL( SCIPsetEventhdlrFree(scip, branchruledata->eventhdlr, eventFreeDistribution) ); 1359 1360 /* include branching rule */ 1361 SCIP_CALL( SCIPincludeBranchruleBasic(scip, &branchrule, BRANCHRULE_NAME, BRANCHRULE_DESC, BRANCHRULE_PRIORITY, BRANCHRULE_MAXDEPTH, 1362 BRANCHRULE_MAXBOUNDDIST, branchruledata) ); 1363 1364 assert(branchrule != NULL); 1365 SCIP_CALL( SCIPsetBranchruleCopy(scip, branchrule, branchCopyDistribution) ); 1366 SCIP_CALL( SCIPsetBranchruleFree(scip, branchrule, branchFreeDistribution) ); 1367 SCIP_CALL( SCIPsetBranchruleExitsol(scip, branchrule, branchExitsolDistribution) ); 1368 SCIP_CALL( SCIPsetBranchruleExecLp(scip, branchrule, branchExeclpDistribution) ); 1369 1370 /* add distribution branching rule parameters */ 1371 SCIP_CALL( SCIPaddCharParam(scip, "branching/" BRANCHRULE_NAME "/scoreparam", 1372 "the score;largest 'd'ifference, 'l'owest cumulative probability,'h'ighest c.p., 'v'otes lowest c.p., votes highest c.p.('w') ", 1373 &branchruledata->scoreparam, TRUE, DEFAULT_SCOREPARAM, SCOREPARAM_VALUES, NULL, NULL) ); 1374 1375 SCIP_CALL( SCIPaddBoolParam(scip, "branching/" BRANCHRULE_NAME "/onlyactiverows", 1376 "should only rows which are active at the current node be considered?", 1377 &branchruledata->onlyactiverows, TRUE, DEFAULT_ONLYACTIVEROWS, NULL, NULL) ); 1378 1379 SCIP_CALL( SCIPaddBoolParam(scip, "branching/" BRANCHRULE_NAME "/weightedscore", 1380 "should the branching score weigh up- and down-scores of a variable", 1381 &branchruledata->usescipscore, TRUE, DEFAULT_USEWEIGHTEDSCORE, NULL, NULL) ); 1382 1383 return SCIP_OKAY; 1384 } 1385