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 event_estim.c 26 * @brief event handler for tree size estimation and restarts 27 * 28 * This event handler plugin provides different methods for approximating the current fraction of the search 29 * that has already been completed and for estimating the total tree size at completion. 30 * It can trigger restarts of the current run if the current run seems hopeless. 31 * 32 * For details about the available approximations of search completion, please see 33 * 34 * Anderson, Hendel, Le Bodic, Pfetsch 35 * Estimating The Size of Branch-and-Bound Trees 36 * under preparation 37 * 38 * This code is a largely enriched version of a code that was used for clairvoyant restarts, see 39 * 40 * Anderson, Hendel, Le Bodic, Viernickel 41 * Clairvoyant Restarts in Branch-and-Bound Search Using Online Tree-Size Estimation 42 * AAAI-19: Proceedings of the Thirty-Third AAAI Conference on Artificial Intelligence, 2018 43 * 44 * @author Gregor Hendel 45 */ 46 47 /*--+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/ 48 49 #include <string.h> 50 #include "blockmemshell/memory.h" 51 #include "scip/event_estim.h" 52 #include "scip/prop_symmetry.h" 53 #include "scip/pub_disp.h" 54 #include "scip/pub_event.h" 55 #include "scip/pub_fileio.h" 56 #include "scip/pub_message.h" 57 #include "scip/pub_misc.h" 58 #include "scip/pub_tree.h" 59 #include "scip/scip_disp.h" 60 #include "scip/scip_event.h" 61 #include "scip/scip_general.h" 62 #include "scip/scip_mem.h" 63 #include "scip/scip_message.h" 64 #include "scip/scip_nlp.h" 65 #include "scip/scip_numerics.h" 66 #include "scip/scip_param.h" 67 #include "scip/scip_pricer.h" 68 #include "scip/scip_sol.h" 69 #include "scip/scip_solve.h" 70 #include "scip/scip_solvingstats.h" 71 #include "scip/scip_table.h" 72 #include "scip/scip_timing.h" 73 #include "scip/scip_tree.h" 74 #include "scip/type_disp.h" 75 #include "scip/type_event.h" 76 #include "scip/type_message.h" 77 #include "scip/type_misc.h" 78 #include "scip/type_retcode.h" 79 #include "scip/type_stat.h" 80 #include "scip/type_table.h" 81 82 #define EVENTHDLR_NAME "estim" 83 #define EVENTHDLR_DESC "event handler for tree size estimation and restarts" 84 #define EVENTTYPE_ESTIM (SCIP_EVENTTYPE_NODEDELETE | SCIP_EVENTTYPE_NODEBRANCHED) 85 86 /* 87 * Data structures 88 */ 89 90 /** enumerator for available restart policies */ 91 enum RestartPolicy 92 { 93 RESTARTPOLICY_NEVER = 0, /**< never restart (disable this event handler) */ 94 RESTARTPOLICY_ALWAYS = 1, /**< always restart (can be fine tuned by using minimum number of nodes and restart limit) */ 95 RESTARTPOLICY_ESTIMATION = 2, /**< base restart on the estimation method */ 96 RESTARTPOLICY_COMPLETION = 3 /**< trigger restart based on search completion approximation */ 97 }; 98 99 typedef enum RestartPolicy RESTARTPOLICY; 100 101 #define RESTARTPOLICY_CHAR_NEVER 'n' 102 #define RESTARTPOLICY_CHAR_ALWAYS 'a' 103 #define RESTARTPOLICY_CHAR_COMPLETION 'c' 104 #define RESTARTPOLICY_CHAR_ESTIMATION 'e' 105 106 #define DES_USETRENDINLEVEL TRUE /**< Should the trend be used in the level update? */ 107 108 /* constants for the table estimation */ 109 #define TABLE_NAME "estim" 110 #define TABLE_DESC "tree size estimations statistics table" 111 #define TABLE_POSITION 18500 /**< the position of the statistics table */ 112 #define TABLE_EARLIEST_STAGE SCIP_STAGE_INIT /**< output of the statistics table is only printed from this stage onwards */ 113 114 /* constants for the search completion display column */ 115 #define DISP_NAME "completed" 116 #define DISP_DESC "completion of search in percent (based on tree size estimation)" 117 #define DISP_HEADER "compl." 118 #define DISP_WIDTH 8 /**< the width of the display column */ 119 #define DISP_PRIORITY 110000 /**< the priority of the display column */ 120 #define DISP_POSITION 30100 /**< the relative position of the display column */ 121 #define DISP_STRIPLINE TRUE /**< the default for whether the display column should be separated 122 * with a line from its right neighbor */ 123 #define INITIALSIZE 100 124 #define SESCOEFF 0.75 /**< coefficient of single exponential smoothing of estimation */ 125 126 /* double exponential smoothing parameters for different time series */ 127 #define DES_ALPHA_TREEWEIGHT 0.65 128 #define DES_BETA_TREEWEIGHT 0.15 129 130 #define DES_ALPHA_GAP 0.6 131 #define DES_BETA_GAP 0.15 132 133 #define DES_ALPHA_LEAFFREQUENCY 0.3 134 #define DES_BETA_LEAFFREQUENCY 0.33 135 136 #define DES_ALPHA_SSG 0.6 137 #define DES_BETA_SSG 0.15 138 139 #define DES_ALPHA_OPENNODES 0.6 140 #define DES_BETA_OPENNODES 0.15 141 142 #define MAX_REGFORESTSIZE 10000000 /**< size limit (number of nodes) for regression forest */ 143 144 145 146 /* computation of search completion */ 147 #define COMPLETIONTYPE_AUTO 'a' /**< automatic (regression forest if available, else monotone regression on binary and SSG on nonbinary trees) */ 148 #define COMPLETIONTYPE_REGFOREST 'r' /**< regression forest (must be provided by user) */ 149 #define COMPLETIONTYPE_MONOREG 'm' /**< monotone regression (using tree weight and SSG) */ 150 #define COMPLETIONTYPE_TREEWEIGHT 'w' /**< use tree weight value as approximation of search tree completion */ 151 #define COMPLETIONTYPE_SSG 's' /**< use SSG value as approximation of search tree completion */ 152 #define COMPLETIONTYPE_GAP 'g' /**< use gap value as approximation of search tree completion */ 153 154 155 /* tree size estimation method */ 156 #define ESTIMMETHOD_COMPL 'c' /**< estimation based on projection of current search completion */ 157 #define ESTIMMETHOD_WBE 'b' /**< weighted backtrack estimation */ 158 #define ESTIMMETHOD_ENSMBL 'e' /**< estimation based on an ensemble of the individual estimations */ 159 #define ESTIMMETHOD_GAP 'g' /**< estimation based on double exponential smoothing for open nodes */ 160 #define ESTIMMETHOD_LFREQ 'l' /**< estimation based on double exponential smoothing for leaf frequency */ 161 #define ESTIMMETHOD_OPEN 'o' /**< estimation based on double exponential smoothing for open nodes */ 162 #define ESTIMMETHOD_SSG 's' /**< estimation based on double exponential smoothing for sum of subtree gaps */ 163 #define ESTIMMETHOD_TPROF 't' /**< estimation based on tree profile method */ 164 #define ESTIMMETHOD_TREEWEIGHT 'w' /**< estimation based on double exponential smoothing for tree weight */ 165 166 #define ESTIMMETHODS "bceglostw" 167 168 /* constants and default values for treeprofile parameters */ 169 #define TREEPROFILE_MINSIZE 512 /**< minimum size (depth) that tree profile can hold */ 170 #define SSG_STARTPRIMBOUND SCIP_INVALID /**< initial value of primal bound used within SSG */ 171 172 /** double exponential smoothing data structure */ 173 struct DoubleExpSmooth 174 { 175 SCIP_Real alpha; /**< level smoothing constant */ 176 SCIP_Real beta; /**< trend smoothing constant */ 177 SCIP_Real level; /**< estimation of the current level used for smoothing */ 178 SCIP_Real trend; /**< estimation of the current trend (slope) */ 179 SCIP_Real initialvalue; /**< the level value at 0 observations */ 180 SCIP_Bool usetrendinlevel; /**< Should the trend be used in the level update? */ 181 int n; /**< number of observations */ 182 }; 183 typedef struct DoubleExpSmooth DOUBLEEXPSMOOTH; 184 185 /** time series data structure for leaf time series 186 * 187 * These time series are the basic ingredient for tree size estimation via forecasting. 188 * 189 * This general class represents concrete time series such as the closed gap, tree weight, and leaf frequency. 190 * Through callbacks for data (de-)initialization and value queries, it provides a common interface 191 * to which double exponential smoothing or window forecasts can be applied. 192 */ 193 typedef struct TimeSeries TIMESERIES; 194 195 /** data structure for convenient access of tree information */ 196 typedef struct TreeData TREEDATA; 197 198 199 #define NTIMESERIES 5 200 201 /** time series position in event handler time series array */ 202 enum TsPos 203 { 204 TSPOS_NONE = -1, /**< invalid array position */ 205 TSPOS_GAP = 0, /**< time series position of gap */ 206 TSPOS_TREEWEIGHT = 1, /**< time series position of tree weight */ 207 TSPOS_LFREQ = 2, /**< time series position of leaf frequency */ 208 TSPOS_SSG = 3, /**< time series position of SSG */ 209 TSPOS_OPEN = 4 /**< time series position of open nodes */ 210 }; 211 212 typedef enum TsPos TSPOS; 213 214 /** regression forest data structure */ 215 typedef struct SCIP_RegForest SCIP_REGFOREST; 216 217 /** statistics collected from profile used for prediction */ 218 struct TreeProfileStats 219 { 220 int maxdepth; /**< maximum node depth encountered */ 221 int lastfulldepth; /**< deepest layer for which all nodes have been explored */ 222 int minwaistdepth; /**< minimum depth of the waist, i.e. the widest part of the tree */ 223 int maxwaistdepth; /**< maximum depth of the waist, i.e. the widest part of the tree */ 224 }; 225 226 typedef struct TreeProfileStats TREEPROFILESTATS; 227 228 229 /** profile data structure for tree */ 230 struct TreeProfile 231 { 232 SCIP_Longint* profile; /**< array to store the tree profile */ 233 int profilesize; /**< size of the profile array */ 234 TREEPROFILESTATS stats; /**< statistics collected from profile used for prediction */ 235 SCIP_Real lastestimate; /**< the last estimate predicted by predictTotalSizeTreeprofile() */ 236 TREEPROFILESTATS lastestimatestats; /**< tree profile statistics at last estimation */ 237 }; 238 239 typedef struct TreeProfile TREEPROFILE; 240 241 /* default values of user parameters */ 242 #define DEFAULT_USELEAFTS TRUE /**< Use leaf nodes as basic observations for time series, or all nodes? */ 243 #define DEFAULT_REPORTFREQ -1 /**< report frequency on estimation: -1: never, 0: always, k >= 1: k times evenly during search */ 244 #define DEFAULT_REGFORESTFILENAME "-" /**< default file name of user regression forest in RFCSV format */ 245 #define DEFAULT_COEFMONOWEIGHT 0.3667 /**< coefficient of tree weight in monotone approximation of search completion */ 246 #define DEFAULT_COEFMONOSSG 0.6333 /**< coefficient of 1 - SSG in monotone approximation of search completion */ 247 #define DEFAULT_COMPLETIONTYPE COMPLETIONTYPE_AUTO /**< default computation of search tree completion */ 248 #define DEFAULT_ESTIMMETHOD ESTIMMETHOD_TREEWEIGHT /**< default tree size estimation method: (c)ompletion, (e)nsemble, time series forecasts on either 249 * (g)ap, (l)eaf frequency, (o)open nodes, 250 * tree (w)eight, (s)sg, or (t)ree profile or w(b)e */ 251 #define DEFAULT_TREEPROFILE_ENABLED FALSE /**< Should the event handler collect data? */ 252 #define DEFAULT_TREEPROFILE_MINNODESPERDEPTH 20.0 /**< minimum average number of nodes at each depth before producing estimations */ 253 #define DEFAULT_RESTARTPOLICY 'e' /**< default restart policy: (a)lways, (c)ompletion, (e)stimation, (n)ever */ 254 #define DEFAULT_RESTARTLIMIT 1 /**< default restart limit */ 255 #define DEFAULT_MINNODES 1000L /**< minimum number of nodes before restart */ 256 #define DEFAULT_COUNTONLYLEAVES FALSE /**< should only leaves count for the minnodes parameter? */ 257 #define DEFAULT_RESTARTFACTOR 50.0 /**< factor by which the estimated number of nodes should exceed the current number of nodes */ 258 #define DEFAULT_RESTARTNONLINEAR FALSE /**< whether to apply a restart when nonlinear constraints are present */ 259 #define DEFAULT_RESTARTACTPRICERS FALSE /**< whether to apply a restart when active pricers are used */ 260 #define DEFAULT_HITCOUNTERLIM 50 /**< limit on the number of successive samples to really trigger a restart */ 261 #define DEFAULT_SSG_NMAXSUBTREES -1 /**< the maximum number of individual SSG subtrees; the old split is kept if 262 * a new split exceeds this number of subtrees ; -1: no limit */ 263 #define DEFAULT_SSG_NMINNODESLASTSPLIT 0L /**< minimum number of nodes to process between two consecutive SSG splits */ 264 #define DEFAULT_SHOWSTATS FALSE /**< should statistics be shown at the end? */ 265 266 /** event handler data */ 267 struct SCIP_EventhdlrData 268 { 269 SCIP_REGFOREST* regforest; /**< regression forest data structure */ 270 TIMESERIES* timeseries[NTIMESERIES]; /**< array of time series slots */ 271 TREEDATA* treedata; /**< tree data */ 272 TREEPROFILE* treeprofile; /**< tree profile data structure */ 273 char* regforestfilename; /**< file name of user regression forest in RFCSV format */ 274 SCIP_Real restartfactor; /**< factor by which the estimated number of nodes should exceed the current number of nodes */ 275 SCIP_Real weightlastreport; /**< tree weight at which last report was printed */ 276 SCIP_Real treeprofile_minnodesperdepth;/**< minimum average number of nodes at each depth before producing estimations */ 277 SCIP_Real coefmonoweight; /**< coefficient of tree weight in monotone approximation of search completion */ 278 SCIP_Real coefmonossg; /**< coefficient of 1 - SSG in monotone approximation of search completion */ 279 SCIP_Longint minnodes; /**< minimum number of nodes in a run before restart is triggered */ 280 int restartlimit; /**< How often should a restart be triggered? (-1 for no limit) */ 281 int nrestartsperformed; /**< number of restarts performed so far */ 282 int restarthitcounter; /**< the number of successive samples that would trigger a restart */ 283 int hitcounterlim; /**< limit on the number of successive samples to really trigger a restart */ 284 int nreports; /**< the number of reports already printed */ 285 int reportfreq; /**< report frequency on estimation: -1: never, 0:always, k >= 1: k times evenly during search */ 286 int lastrestartrun; /**< the last run at which this event handler triggered restart */ 287 char restartpolicyparam; /**< restart policy parameter */ 288 char estimmethod; /**< tree size estimation method: (c)ompletion, (e)nsemble, time series forecasts on either 289 * (g)ap, (l)eaf frequency, (o)open nodes, 290 * tree (w)eight, (s)sg, or (t)ree profile or w(b)e */ 291 char completiontypeparam;/**< approximation of search tree completion: 292 * (a)uto, (g)ap, tree (w)eight, (m)onotone regression, (r)egression forest, (s)sg */ 293 SCIP_Bool countonlyleaves; /**< Should only leaves count for the minnodes parameter? */ 294 SCIP_Bool useleafts; /**< Use leaf nodes as basic observations for time series, or all nodes? */ 295 SCIP_Bool treeprofile_enabled;/**< Should the event handler collect treeprofile data? */ 296 SCIP_Bool treeisbinary; /**< internal flag if all branching decisions produced 2 children */ 297 SCIP_Bool restartnonlinear; /**< whether to apply a restart when nonlinear constraints are present */ 298 SCIP_Bool restartactpricers; /**< whether to apply a restart when active pricers are used */ 299 SCIP_Bool showstats; /**< should statistics be shown at the end? */ 300 }; 301 302 typedef struct SubtreeSumGap SUBTREESUMGAP; 303 304 struct TreeData 305 { 306 SCIP_Longint nnodes; /**< the total number of nodes */ 307 SCIP_Longint nopen; /**< the current number of open nodes */ 308 SCIP_Longint ninner; /**< the number of inner nodes */ 309 SCIP_Longint nleaves; /**< the number of final leaf nodes */ 310 SCIP_Longint nvisited; /**< the number of visited nodes */ 311 long double weight; /**< the current tree weight (sum of leaf weights) */ 312 SUBTREESUMGAP* ssg; /**< subtree sum gap data structure */ 313 }; 314 315 struct SubtreeSumGap 316 { 317 SCIP_Real value; /**< the current subtree sum gap */ 318 SCIP_HASHMAP* nodes2info; /**< map between nodes and their subtree indices */ 319 SCIP_PQUEUE** subtreepqueues; /**< array of priority queues, one for each subtree */ 320 SCIP_Real scalingfactor; /**< the current scaling factor */ 321 SCIP_Real pblastsplit; /**< primal bound when last split occurred */ 322 SCIP_Longint nodelastsplit; /**< last node at which a subtree split occurred */ 323 SCIP_Longint nminnodeslastsplit; /**< minimum number of nodes to process between two consecutive SSG splits */ 324 int nmaxsubtrees; /**< the maximum number of individual SSG subtrees; the old split is kept if 325 * a new split exceeds this number of subtrees ; -1: no limit */ 326 int nsubtrees; /**< the current number n of subtrees labeled 0 .. n - 1 */ 327 }; 328 329 /** update callback of time series */ 330 #define DECL_TIMESERIESUPDATE(x) SCIP_RETCODE x (\ 331 SCIP* scip, \ 332 TIMESERIES* ts, \ 333 TREEDATA* treedata, \ 334 SCIP_Real* value \ 335 ) 336 337 /** time series data structure for leaf time series */ 338 struct TimeSeries 339 { 340 DOUBLEEXPSMOOTH des; /**< double exponential smoothing data structure */ 341 char* name; /**< name of this time series */ 342 SCIP_Real* vals; /**< value array of this time series */ 343 SCIP_Real* estimation; /**< array of estimations of this time series */ 344 SCIP_Real smoothestimation; /**< smoothened estimation value */ 345 SCIP_Real targetvalue; /**< target value of this time series */ 346 SCIP_Real currentvalue; /**< current value of time series */ 347 SCIP_Real initialvalue; /**< the initial value of time series */ 348 SCIP_Longint nobs; /**< total number of observations */ 349 int valssize; /**< size of value array */ 350 int nvals; /**< number of values */ 351 int resolution; /**< current (inverse of) resolution */ 352 SCIP_Bool useleafts; /**< Should this time series be recorded at leaf nodes, or at every node? */ 353 DECL_TIMESERIESUPDATE((*timeseriesupdate));/**< update callback at nodes */ 354 }; 355 356 /** extended node information for SSG priority queue */ 357 struct NodeInfo 358 { 359 SCIP_NODE* node; /**< search tree node */ 360 SCIP_Real lowerbound; /**< lower bound of the node at insertion into priority queue */ 361 int pos; /**< position of this node in priority queue */ 362 int subtreeidx; /**< subtree index of this node */ 363 }; 364 typedef struct NodeInfo NODEINFO; 365 366 struct SCIP_RegForest 367 { 368 int ntrees; /**< number of trees in this forest */ 369 int dim; /**< feature dimension */ 370 int* nbegin; /**< array of root node indices of each tree */ 371 int* child; /**< child index pair of each internal node, or (-1, -1) for leaves */ 372 int* splitidx; /**< data index for split at node, or -1 at a leaf */ 373 SCIP_Real* value; /**< split position at internal nodes, prediction at leaves */ 374 int size; /**< length of node arrays */ 375 }; 376 377 /* 378 * Local methods 379 */ 380 381 /** convert number to string and treat SCIP_INVALID as '-' */ 382 static 383 char* real2String( 384 SCIP_Real num, /**< number to convert to string */ 385 char* buf, /**< string buffer */ 386 int digits /**< number of decimal digits */ 387 ) 388 { 389 if( num == SCIP_INVALID )/*lint !e777*/ 390 (void) SCIPsnprintf(buf, 1, "-"); 391 else if( num >= 1e+20 ) /*lint !e777*/ 392 (void) SCIPsnprintf(buf, 3, "inf"); 393 else 394 (void) SCIPsnprintf(buf, SCIP_MAXSTRLEN, "%10.*f", digits, num); 395 396 return buf; 397 } 398 399 /** free a regression forest data structure */ 400 static 401 void SCIPregForestFree( 402 SCIP_REGFOREST** regforest /**< regression forest data structure */ 403 ) 404 { 405 SCIP_REGFOREST* regforestptr; 406 407 assert(regforest != NULL); 408 409 if( *regforest == NULL ) 410 return; 411 regforestptr = *regforest; 412 413 BMSfreeMemoryArrayNull(®forestptr->nbegin); 414 BMSfreeMemoryArrayNull(®forestptr->child); 415 BMSfreeMemoryArrayNull(®forestptr->splitidx); 416 BMSfreeMemoryArrayNull(®forestptr->value); 417 418 BMSfreeMemory(regforest); 419 } 420 421 /** make a prediction with regression forest */ 422 static 423 SCIP_Real SCIPregForestPredict( 424 SCIP_REGFOREST* regforest, /**< regression forest data structure */ 425 SCIP_Real* datapoint /**< a data point that matches the dimension of this regression forest */ 426 ) 427 { 428 int treeidx; 429 SCIP_Real value = 0.0; 430 431 assert(regforest != NULL); 432 assert(datapoint != NULL); 433 434 SCIPdebugMessage("Start prediction method of regression forest\n"); 435 436 /* loop through the trees */ 437 for( treeidx = 0; treeidx < regforest->ntrees; ++treeidx ) 438 { 439 int treepos = regforest->nbegin[treeidx]; 440 int* childtree = &(regforest->child[2 * treepos]); 441 int* splitidxtree = &(regforest->splitidx[treepos]); 442 int pos = 0; 443 SCIP_Real* valuetree = &(regforest->value[treepos]); 444 445 SCIPdebugMessage("Tree %d at position %d\n", treeidx, treepos); 446 447 /* find the correct leaf */ 448 while( splitidxtree[pos] != - 1 ) 449 { 450 int goright; 451 452 assert(splitidxtree[pos] < regforest->dim); 453 454 goright = (datapoint[splitidxtree[pos]] > valuetree[pos]) ? 1 : 0; 455 pos = childtree[2 * pos + goright]; 456 } 457 458 value += valuetree[pos]; 459 } 460 461 /* return the average value that the trees predict */ 462 return value / (SCIP_Real)(regforest->ntrees); 463 } 464 465 /** read a regression forest from an rfcsv file 466 * 467 * TODO improve this parser to better capture wrong user input, e.g., if the dimension is wrong 468 */ 469 static 470 SCIP_RETCODE SCIPregForestFromFile( 471 SCIP_REGFOREST** regforest, /**< regression forest data structure */ 472 const char* filename /**< name of file with the regression forest data */ 473 ) 474 { 475 SCIP_RETCODE retcode = SCIP_OKAY; 476 SCIP_FILE* file; 477 SCIP_REGFOREST* regforestptr; 478 char buffer[SCIP_MAXSTRLEN]; 479 char firstlineformat[SCIP_MAXSTRLEN]; 480 char dataformat[SCIP_MAXSTRLEN]; 481 char valuestr[SCIP_MAXSTRLEN]; 482 SCIP_Bool error = FALSE; 483 int ntrees; 484 int dim; 485 int size; 486 int sscanret; 487 int pos; 488 int treepos; 489 490 /* try to open file */ 491 file = SCIPfopen(filename, "r"); 492 493 if( file == NULL ) 494 return SCIP_NOFILE; 495 496 /* parse read the first line that contains the number of trees, feature dimension, and total number of nodes */ 497 (void) SCIPsnprintf(firstlineformat, SCIP_MAXSTRLEN, "### NTREES=%%10d FEATURE_DIM=%%10d LENGTH=%%10d\n"); 498 if( SCIPfgets(buffer, (int) sizeof(buffer), file) == NULL ) 499 { 500 error = TRUE; 501 SCIPerrorMessage("Could not read first line of regression file '%s'\n", filename); 502 goto CLOSEFILE; 503 } 504 505 /* coverity[secure_coding] */ 506 sscanret = sscanf(buffer, firstlineformat, &ntrees, &dim, &size); 507 508 if( sscanret != 3 ) 509 { 510 error = TRUE; 511 SCIPerrorMessage("Could not extract tree information from buffer line [%s]\n", buffer); 512 goto CLOSEFILE; 513 } 514 515 SCIPdebugMessage("Read ntrees=%d, dim=%d, size=%d (return value %d)\n", ntrees, dim, size, sscanret); 516 517 /* check if the tree is too big, or numbers are negative */ 518 if( size > MAX_REGFORESTSIZE ) 519 { 520 error = TRUE; 521 SCIPerrorMessage("Requested size %d exceeds size limit %d for regression trees", size, MAX_REGFORESTSIZE); 522 goto CLOSEFILE; 523 } 524 525 if( dim <= 0 || ntrees <= 0 || size <= 0 ) 526 { 527 error = TRUE; 528 SCIPerrorMessage("Cannot create regression tree with negative size, dimension, or number of trees\n"); 529 goto CLOSEFILE; 530 } 531 532 /* allocate memory in regression forest data structure */ 533 SCIP_ALLOC_TERMINATE( retcode, BMSallocMemory(regforest), FREEFOREST ); 534 BMSclearMemory(*regforest); 535 regforestptr = *regforest; 536 537 /* coverity[tainted_data] */ 538 SCIP_ALLOC_TERMINATE( retcode, BMSallocMemoryArray(®forestptr->nbegin, ntrees), FREEFOREST ); 539 SCIP_ALLOC_TERMINATE( retcode, BMSallocMemoryArray(®forestptr->child, 2 * size), FREEFOREST ); /*lint !e647*/ 540 SCIP_ALLOC_TERMINATE( retcode, BMSallocMemoryArray(®forestptr->splitidx, size), FREEFOREST ); 541 SCIP_ALLOC_TERMINATE( retcode, BMSallocMemoryArray(®forestptr->value, size), FREEFOREST ); 542 543 regforestptr->dim = dim; 544 regforestptr->size = size; 545 regforestptr->ntrees = ntrees; 546 547 SCIPdebugMessage("Random Forest allocated\n"); 548 549 /* loop through the rest of the file, which contains the comma separated node data */ 550 (void) SCIPsnprintf(dataformat, SCIP_MAXSTRLEN, "%%10d,%%10d,%%10d,%%10d,%%%ds\n", SCIP_MAXSTRLEN); 551 552 pos = 0; 553 treepos = 0; 554 while( !SCIPfeof(file) && !error ) 555 { 556 int node; 557 char* endptr; 558 559 /* get next line */ 560 if( SCIPfgets(buffer, (int) sizeof(buffer), file) == NULL ) 561 break; 562 563 sscanret = sscanf(buffer, dataformat, 564 &node, 565 ®forestptr->child[2 * pos], 566 ®forestptr->child[2 * pos + 1], 567 ®forestptr->splitidx[pos], 568 valuestr); 569 570 if( sscanret != 5 ) 571 { 572 SCIPerrorMessage("Something wrong with line %d '%s'", pos + 1, buffer); 573 error = TRUE; 574 } 575 576 (void)SCIPstrToRealValue(valuestr, ®forestptr->value[pos], &endptr); 577 578 /* new root node - increase the tree index position */ 579 if( node == 0 ) 580 { 581 assert(treepos < regforestptr->ntrees); 582 583 regforestptr->nbegin[treepos++] = pos; 584 } 585 586 ++pos; 587 } 588 589 goto CLOSEFILE; 590 591 /* insufficient memory for allocating regression forest */ 592 FREEFOREST: 593 assert(retcode == SCIP_NOMEMORY); 594 SCIPregForestFree(regforest); 595 596 CLOSEFILE: 597 SCIPfclose(file); 598 599 if( error ) 600 retcode = SCIP_INVALIDDATA; 601 602 return retcode; 603 } 604 605 /** compare two tree profile statistics for equality */ 606 static 607 SCIP_Bool isEqualTreeProfileStats( 608 TREEPROFILESTATS* stats, /**< first tree profile statistics */ 609 TREEPROFILESTATS* other /**< other tree profile statistics */ 610 ) 611 { 612 assert(stats != NULL); 613 assert(other != NULL); 614 615 return stats->maxdepth == other->maxdepth && 616 stats->lastfulldepth == other->lastfulldepth && 617 stats->minwaistdepth == other->minwaistdepth && 618 stats->maxwaistdepth == other->maxwaistdepth; 619 } 620 621 /** copy source tree profile into destination */ 622 static 623 void copyTreeProfileStats( 624 TREEPROFILESTATS* dest, /**< destination tree profile statistics */ 625 TREEPROFILESTATS* src /**< source tree profile statistics */ 626 ) 627 { 628 assert(dest != NULL); 629 assert(src != NULL); 630 631 dest->maxdepth = src->maxdepth; 632 dest->lastfulldepth = src->lastfulldepth; 633 dest->minwaistdepth = src->minwaistdepth; 634 dest->maxwaistdepth = src->maxwaistdepth; 635 } 636 637 /** reset tree profile statistics */ 638 static 639 void resetTreeProfileStats( 640 TREEPROFILESTATS* treeprofilestats /**< tree profile statistics */ 641 ) 642 { 643 assert(treeprofilestats != NULL); 644 645 BMSclearMemory(treeprofilestats); 646 } 647 648 649 /** extend tree profile to deeper tree */ 650 static 651 SCIP_RETCODE extendMemoryTreeProfile( 652 SCIP* scip, /**< SCIP data structure */ 653 TREEPROFILE* treeprofile, /**< tree profile data structure */ 654 int mindepth /**< minimum depth that the tree profile should hold */ 655 ) 656 { 657 if( mindepth < treeprofile->profilesize ) 658 return SCIP_OKAY; 659 660 if( treeprofile->profile == NULL ) 661 { 662 SCIP_CALL( SCIPallocClearMemoryArray(scip, &treeprofile->profile, mindepth) ); 663 treeprofile->profilesize = mindepth; 664 } 665 else 666 { 667 int newsize; 668 int nnewelems; 669 SCIP_Longint* newprofile; 670 671 newsize = SCIPcalcMemGrowSize(scip, mindepth + 1); 672 nnewelems = newsize - treeprofile->profilesize; 673 assert(newsize > treeprofile->profilesize); 674 675 SCIP_CALL( SCIPreallocMemoryArray(scip, &treeprofile->profile, newsize) ); 676 newprofile = &treeprofile->profile[treeprofile->profilesize]; 677 BMSclearMemoryArray(newprofile, nnewelems); 678 treeprofile->profilesize = newsize; 679 } 680 681 return SCIP_OKAY; 682 } 683 684 /** create a tree profile */ 685 static 686 SCIP_RETCODE createTreeProfile( 687 SCIP* scip, /**< SCIP data structure */ 688 TREEPROFILE** treeprofile /**< pointer to store tree profile data structure */ 689 ) 690 { 691 assert(scip != NULL); 692 assert(treeprofile != NULL); 693 694 SCIP_CALL( SCIPallocMemory(scip, treeprofile) ); 695 696 (*treeprofile)->profile = NULL; 697 (*treeprofile)->profilesize = 0; 698 SCIP_CALL( extendMemoryTreeProfile(scip, *treeprofile, TREEPROFILE_MINSIZE) ); 699 700 resetTreeProfileStats(&(*treeprofile)->stats); 701 resetTreeProfileStats(&(*treeprofile)->lastestimatestats); 702 703 (*treeprofile)->lastestimate = -1.0; 704 705 return SCIP_OKAY; 706 } 707 708 /** free a tree profile */ 709 static 710 void freeTreeProfile( 711 SCIP* scip, /**< SCIP data structure */ 712 TREEPROFILE** treeprofile /**< pointer to tree profile data structure */ 713 ) 714 { 715 assert(scip != NULL); 716 assert(treeprofile != NULL); 717 718 if( *treeprofile == NULL ) 719 return; 720 721 SCIPfreeMemoryArray(scip, &(*treeprofile)->profile); 722 723 SCIPfreeMemory(scip, treeprofile); 724 725 *treeprofile = NULL; 726 } 727 728 /** update tree profile */ 729 static 730 SCIP_RETCODE updateTreeProfile( 731 SCIP* scip, /**< SCIP data structure */ 732 TREEPROFILE* treeprofile, /**< tree profile data structure */ 733 SCIP_NODE* node /**< node that should be added to the profile */ 734 ) 735 { 736 int nodedepth; 737 unsigned long nbits; 738 SCIP_Longint nodedepthcnt; 739 SCIP_Longint maxnodes; 740 741 assert(scip != NULL); 742 assert(node != NULL); 743 744 if( treeprofile == NULL ) 745 return SCIP_OKAY; 746 747 nodedepth = SCIPnodeGetDepth(node); 748 assert(nodedepth >= 0); 749 maxnodes = treeprofile->profile[treeprofile->stats.minwaistdepth]; 750 assert(treeprofile->stats.minwaistdepth == treeprofile->stats.maxwaistdepth || 751 maxnodes == treeprofile->profile[treeprofile->stats.maxwaistdepth]); 752 753 /* ensure that the memory can hold at least this depth */ 754 SCIP_CALL( extendMemoryTreeProfile(scip, treeprofile, nodedepth) ); 755 756 nodedepthcnt = ++treeprofile->profile[nodedepth]; 757 758 /* Is this level fully explored? We assume binary branching. The first condition ensures that the bit shift operation 759 * of the second condition represents a feasible power of unsigned int. The largest power of 2 representable 760 * by unsigned int is 2^{8*sizeof(unsigned int) - 1}. */ 761 nbits = 8*sizeof(unsigned int); 762 /* coverity[overflow_before_widen] */ 763 if( (unsigned int)nodedepth < nbits && nodedepthcnt == (1U << nodedepth) )/*lint !e647*/ 764 { 765 SCIPdebugMsg(scip, "Level %d fully explored: %" SCIP_LONGINT_FORMAT " nodes\n", nodedepth, nodedepthcnt); 766 767 treeprofile->stats.lastfulldepth = nodedepth; 768 } 769 770 /* update maximum depth */ 771 if( treeprofile->stats.maxdepth < nodedepth ) 772 { 773 treeprofile->stats.maxdepth = nodedepth; 774 SCIPdebugMsg(scip, "Maximum depth increased to %d\n", treeprofile->stats.maxdepth); 775 } 776 777 /* minimum and maximum waist now coincide */ 778 if( nodedepthcnt > maxnodes ) 779 { 780 treeprofile->stats.minwaistdepth = treeprofile->stats.maxwaistdepth = nodedepth; 781 SCIPdebugMsg(scip, "Updating depth of tree waist: %d (%" SCIP_LONGINT_FORMAT " nodes)\n", 782 treeprofile->stats.minwaistdepth, nodedepthcnt); 783 } 784 else if( nodedepthcnt == maxnodes ) 785 { 786 /* enlarge the interval in which the waist lies */ 787 if( treeprofile->stats.minwaistdepth > nodedepth ) 788 treeprofile->stats.minwaistdepth = nodedepth; 789 else if( treeprofile->stats.maxwaistdepth < nodedepth ) 790 treeprofile->stats.maxwaistdepth = nodedepth; 791 } 792 assert(treeprofile->stats.minwaistdepth <= treeprofile->stats.maxwaistdepth); 793 794 return SCIP_OKAY; 795 } 796 797 /** make a prediction of the total tree size based on the current tree profile */ 798 static 799 SCIP_Real predictTotalSizeTreeProfile( 800 SCIP* scip, /**< SCIP data structure */ 801 TREEPROFILE* treeprofile, /**< tree profile data structure */ 802 SCIP_Real minnodesperdepth /**< minimum number of average nodes per depth to make a prediction */ 803 ) 804 { 805 SCIP_Real estimate; 806 SCIP_Real growthfac; 807 int d; 808 int waist; 809 810 /* prediction is disabled */ 811 if( treeprofile == NULL ) 812 return -1.0; 813 814 /* two few nodes to make a prediction */ 815 if( minnodesperdepth * treeprofile->stats.maxdepth > SCIPgetNNodes(scip) ) 816 return -1.0; 817 818 /* reuse previous estimation if tree profile hasn't changed */ 819 if( isEqualTreeProfileStats(&treeprofile->lastestimatestats, &treeprofile->stats) ) 820 { 821 SCIPdebugMsg(scip, "Reusing previous estimation result %g\n", treeprofile->lastestimate); 822 823 return treeprofile->lastestimate; 824 } 825 826 /* compute the (depth of the) waist as convex combination between the minimum and maximum waist depths */ 827 waist = (2 * treeprofile->stats.maxwaistdepth + treeprofile->stats.minwaistdepth) / 3; 828 829 growthfac = 2; 830 estimate = 1; 831 832 /* loop over all full levels */ 833 for( d = 1; d < treeprofile->stats.lastfulldepth; ++d ) 834 { 835 SCIP_Real gamma_d = 2.0; 836 837 estimate += growthfac; 838 growthfac *= gamma_d; 839 } 840 841 /* loop until the waist is reached */ 842 for( ; d < waist; ++d ) 843 { 844 SCIP_Real gamma_d = 2.0 - (d - treeprofile->stats.lastfulldepth + 1.0)/(waist - treeprofile->stats.lastfulldepth + 1.0); 845 846 assert(1.0 <= gamma_d && gamma_d <= 2.0); 847 estimate += growthfac; 848 growthfac *= gamma_d; 849 } 850 851 /* loop over the remaining levels */ 852 for( ; d <= treeprofile->stats.maxdepth; ++d ) 853 { 854 SCIP_Real gamma_d = (1.0 - (d - waist + 1.0)/(treeprofile->stats.maxdepth - waist + 1.0)); 855 assert(0.0 <= gamma_d && gamma_d <= 1.0); 856 857 estimate += growthfac; 858 growthfac *= gamma_d; 859 } 860 861 /* copy tree profile statistics */ 862 copyTreeProfileStats(&treeprofile->lastestimatestats, &treeprofile->stats); 863 864 treeprofile->lastestimate = estimate; 865 866 return estimate; 867 } 868 869 /** clean subtrees stored as priority queues */ 870 static 871 void subtreeSumGapDelSubtrees( 872 SCIP* scip, /**< SCIP data structure */ 873 SUBTREESUMGAP* ssg /**< subtree sum gap data structure */ 874 ) 875 { 876 assert(ssg->nsubtrees <= 1 || ssg->subtreepqueues != NULL); 877 878 /* free all previous priority queues */ 879 if( ssg->nsubtrees > 1 ) 880 { 881 int s; 882 883 for( s = 0; s < ssg->nsubtrees; ++s ) 884 { 885 int i; 886 SCIP_PQUEUE* pqueue = ssg->subtreepqueues[s]; 887 NODEINFO** nodeinfos; 888 889 assert(pqueue != NULL); 890 nodeinfos = (NODEINFO**)SCIPpqueueElems(pqueue); 891 892 /* free all remaining elements in reverse order */ 893 for( i = SCIPpqueueNElems(pqueue); --i >= 0; ) 894 { 895 NODEINFO* nodeinfo = nodeinfos[i]; 896 assert(nodeinfo != NULL); 897 SCIPfreeBlockMemory(scip, &nodeinfo); 898 } 899 900 SCIPpqueueFree(&pqueue); 901 } 902 903 SCIPfreeBlockMemoryArray(scip, &ssg->subtreepqueues, ssg->nsubtrees); 904 } 905 906 ssg->subtreepqueues = NULL; 907 } 908 909 /** reset subtree sum gap */ 910 static 911 SCIP_RETCODE subtreeSumGapReset( 912 SCIP* scip, /**< SCIP data structure */ 913 SUBTREESUMGAP* ssg /**< subtree sum gap data structure */ 914 ) 915 { 916 assert(ssg != NULL); 917 assert(ssg->nodes2info != NULL); 918 919 SCIP_CALL( SCIPhashmapRemoveAll(ssg->nodes2info) ); 920 921 subtreeSumGapDelSubtrees(scip, ssg); 922 923 ssg->value = 1.0; 924 ssg->scalingfactor = 1.0; 925 ssg->nsubtrees = 1; 926 ssg->subtreepqueues = NULL; 927 ssg->pblastsplit = SSG_STARTPRIMBOUND; 928 ssg->nodelastsplit = -1L; 929 930 return SCIP_OKAY; 931 } 932 933 /** create a subtree sum gap */ 934 static 935 SCIP_RETCODE subtreeSumGapCreate( 936 SCIP* scip, /**< SCIP data structure */ 937 SUBTREESUMGAP** ssg /**< pointer to store subtree sum gap data structure */ 938 ) 939 { 940 assert(scip != NULL); 941 assert(ssg != NULL); 942 943 /* allocate storage */ 944 SCIP_CALL( SCIPallocMemory(scip, ssg) ); 945 SCIP_CALL( SCIPhashmapCreate(&(*ssg)->nodes2info, SCIPblkmem(scip), INITIALSIZE) ); 946 947 /* explicitly set this to skip removal of subtrees during reset */ 948 (*ssg)->nsubtrees = 0; 949 950 /* reset ssg */ 951 SCIP_CALL( subtreeSumGapReset(scip, *ssg) ); 952 953 return SCIP_OKAY; 954 } 955 956 /** free a subtree sum gap */ 957 static 958 void subtreeSumGapFree( 959 SCIP* scip, /**< SCIP data structure */ 960 SUBTREESUMGAP** ssg /**< pointer to store subtree sum gap data structure */ 961 ) 962 { 963 assert(scip != NULL); 964 965 if( *ssg == NULL ) 966 return; 967 968 if( (*ssg)->nodes2info != NULL ) 969 { 970 SCIPhashmapFree(&(*ssg)->nodes2info); 971 } 972 973 /* delete all subtree data */ 974 subtreeSumGapDelSubtrees(scip, *ssg); 975 976 SCIPfreeMemory(scip, ssg); 977 } 978 979 /** compare two node infos by comparing their lower bound */ 980 static 981 SCIP_DECL_SORTPTRCOMP(compareNodeInfos) 982 { 983 NODEINFO* nodeinfo1 = (NODEINFO*)elem1; 984 NODEINFO* nodeinfo2 = (NODEINFO*)elem2; 985 986 if( nodeinfo1->lowerbound < nodeinfo2->lowerbound ) 987 return -1; 988 else if( nodeinfo1->lowerbound > nodeinfo2->lowerbound ) 989 return 1; 990 991 return 0; 992 } 993 994 /** position change callback of element in priority queue */ 995 static 996 SCIP_DECL_PQUEUEELEMCHGPOS(elemChgPosNodeInfo) 997 { 998 NODEINFO* nodeinfo = (NODEINFO*)elem; 999 1000 assert(oldpos == -1 || oldpos == nodeinfo->pos); 1001 nodeinfo->pos = newpos; 1002 } 1003 1004 /** store node in SSG data structure */ 1005 static 1006 SCIP_RETCODE subtreeSumGapStoreNode( 1007 SCIP* scip, /**< SCIP data structure */ 1008 SUBTREESUMGAP* ssg, /**< subtree sum gap data structure */ 1009 SCIP_NODE* node, /**< node that should be stored */ 1010 int subtreeidx /**< subtree index of that node */ 1011 ) 1012 { 1013 NODEINFO* nodeinfo; 1014 1015 assert(scip != NULL); 1016 assert(ssg != NULL); 1017 assert(node != NULL); 1018 1019 /* create a new node info */ 1020 SCIP_CALL( SCIPallocBlockMemory(scip, &nodeinfo) ); 1021 1022 /* store node information in data structure and insert into priority queue */ 1023 nodeinfo->node = node; 1024 nodeinfo->subtreeidx = subtreeidx; 1025 nodeinfo->pos = -1; 1026 nodeinfo->lowerbound = SCIPnodeGetLowerbound(node); 1027 1028 SCIPdebugMsg(scip, "Inserting label %d for node number %" SCIP_LONGINT_FORMAT " (%p)\n", 1029 subtreeidx, SCIPnodeGetNumber(node), (void*)node); 1030 1031 assert(!SCIPhashmapExists(ssg->nodes2info, (void*)node)); 1032 /* store node information in Hash Map */ 1033 SCIP_CALL( SCIPhashmapInsert(ssg->nodes2info, (void*)node, (void*)nodeinfo) ); 1034 1035 /* create the corresponding priority queue, if it does not exist yet */ 1036 assert(subtreeidx >= 0); 1037 assert(subtreeidx < ssg->nsubtrees); 1038 1039 if( ssg->subtreepqueues[subtreeidx] == NULL ) 1040 { 1041 SCIP_CALL( SCIPpqueueCreate(&ssg->subtreepqueues[subtreeidx], 5, 1.2, compareNodeInfos, elemChgPosNodeInfo) ); 1042 } 1043 1044 /* insert node and ensure that its position is up to date */ 1045 SCIP_CALL( SCIPpqueueInsert(ssg->subtreepqueues[subtreeidx], (void*)nodeinfo) ); 1046 assert(0 <= nodeinfo->pos); 1047 assert(SCIPpqueueNElems(ssg->subtreepqueues[subtreeidx]) > nodeinfo->pos); 1048 assert(SCIPpqueueElems(ssg->subtreepqueues[subtreeidx])[nodeinfo->pos] == (void*)nodeinfo); 1049 1050 return SCIP_OKAY; 1051 } 1052 1053 /** split the open nodes of the current tree */ 1054 static 1055 SCIP_RETCODE subtreeSumGapSplit( 1056 SCIP* scip, /**< SCIP data structure */ 1057 SUBTREESUMGAP* ssg, /**< subtree sum gap data structure */ 1058 SCIP_Bool addfocusnode /**< should the focus node be a subtree, too? */ 1059 ) 1060 { 1061 SCIP_NODE** opennodes[3]; 1062 int nopennodes[3]; 1063 int label; 1064 int t; 1065 int nnewsubtrees; 1066 1067 assert(scip != NULL); 1068 assert(ssg != NULL); 1069 1070 /* query the open nodes of SCIP */ 1071 SCIP_CALL( SCIPgetOpenNodesData(scip, &opennodes[0], &opennodes[1], &opennodes[2], &nopennodes[0], &nopennodes[1], &nopennodes[2]) ); 1072 1073 nnewsubtrees = nopennodes[0] + nopennodes[1] + nopennodes[2] + (addfocusnode ? 1 : 0); 1074 1075 /* clear hash map from entries */ 1076 SCIP_CALL( SCIPhashmapRemoveAll(ssg->nodes2info) ); 1077 1078 /* delete all subtrees */ 1079 subtreeSumGapDelSubtrees(scip, ssg); 1080 1081 ssg->nsubtrees = nnewsubtrees; 1082 SCIPdebugMsg(scip, "Splitting tree into %d subtrees\n", ssg->nsubtrees); 1083 1084 /* create priority queue array */ 1085 if( ssg->nsubtrees > 1 ) 1086 { 1087 SCIP_CALL( SCIPallocClearBlockMemoryArray(scip, &ssg->subtreepqueues, ssg->nsubtrees) ); 1088 } 1089 else 1090 { 1091 ssg->subtreepqueues = NULL; 1092 1093 return SCIP_OKAY; 1094 } 1095 1096 /* loop over node types (leaves, siblings, children) */ 1097 label = 0; 1098 for( t = 0; t < 3; ++t ) 1099 { 1100 SCIP_NODE** nodes = opennodes[t]; 1101 int nnodes = nopennodes[t]; 1102 int n; 1103 1104 /* label each open node as new, separate subtree */ 1105 for( n = 0; n < nnodes; ++n ) 1106 { 1107 SCIP_NODE* node = nodes[n]; 1108 SCIP_CALL( subtreeSumGapStoreNode(scip, ssg, node, label++) ); 1109 } 1110 } 1111 1112 if( addfocusnode ) 1113 { 1114 assert(SCIPgetFocusNode(scip) != NULL); 1115 SCIP_CALL( subtreeSumGapStoreNode(scip, ssg, SCIPgetFocusNode(scip), label) ); 1116 } 1117 1118 return SCIP_OKAY; 1119 } 1120 1121 /** compute a gap between a lower bound and the current upper bound */ 1122 static 1123 SCIP_Real calcGap( 1124 SCIP* scip, /**< SCIP data structure */ 1125 SCIP_Real lowerbound /**< lower bound value */ 1126 ) 1127 { 1128 SCIP_Real db; 1129 SCIP_Real pb; 1130 SCIP_Real abspb; 1131 SCIP_Real absdb; 1132 SCIP_Real gap; 1133 1134 if( SCIPisInfinity(scip, lowerbound) || lowerbound >= SCIPgetUpperbound(scip) ) 1135 return 0.0; 1136 1137 if( SCIPisInfinity(scip, SCIPgetUpperbound(scip)) ) 1138 return 1.0; 1139 1140 db = SCIPretransformObj(scip, lowerbound); 1141 pb = SCIPgetPrimalbound(scip); 1142 1143 if( SCIPisEQ(scip, db, pb) ) 1144 return 0.0; 1145 1146 abspb = REALABS(pb); 1147 absdb = REALABS(db); 1148 gap = REALABS(pb - db)/MAX(abspb,absdb); 1149 gap = MIN(gap, 1.0); 1150 1151 return gap; 1152 } 1153 1154 /** remove node from the subtree sum gap (because it has been solved by branching or is a leaf) */ 1155 static 1156 SCIP_RETCODE subtreeSumGapRemoveNode( 1157 SCIP* scip, /**< SCIP data structure */ 1158 SUBTREESUMGAP* ssg, /**< subtree sum gap data structure */ 1159 SCIP_NODE* node /**< node that should be removed */ 1160 ) 1161 { 1162 NODEINFO* nodeinfo; 1163 int subtreeidx; 1164 int pos; 1165 SCIP_PQUEUE* pqueue; 1166 1167 if( ssg->nsubtrees <= 1 ) 1168 return SCIP_OKAY; 1169 1170 nodeinfo = (NODEINFO*)SCIPhashmapGetImage(ssg->nodes2info, (void*)node); 1171 1172 /* it can happen that the node was not created via branching; search for the most recent ancestor in the queue */ 1173 if( nodeinfo == NULL ) 1174 { 1175 do 1176 { 1177 node = SCIPnodeGetParent(node); 1178 } while( node != NULL && (nodeinfo = (NODEINFO*)SCIPhashmapGetImage(ssg->nodes2info, (void*)node)) == NULL); 1179 1180 /* no ancestor found */ 1181 if( nodeinfo == NULL ) 1182 return SCIP_OKAY; 1183 } 1184 1185 /* get open nodes of this subtree stored as priority queue */ 1186 subtreeidx = nodeinfo->subtreeidx; 1187 pqueue = ssg->subtreepqueues[subtreeidx]; 1188 assert(pqueue != NULL); 1189 1190 /* delete the element from the priority queue */ 1191 pos = nodeinfo->pos; 1192 assert(pos >= 0); 1193 assert(pos < SCIPpqueueNElems(pqueue)); 1194 assert(SCIPpqueueElems(pqueue)[pos] == (void *)nodeinfo); 1195 SCIPpqueueDelPos(pqueue, pos); 1196 1197 /* update ssg if removed node was the lower bound defining node of its subtree */ 1198 if( pos == 0 ) 1199 { 1200 NODEINFO* nodeinfofirst; 1201 SCIP_Real oldgap; 1202 SCIP_Real newgap; 1203 1204 oldgap = calcGap(scip, nodeinfo->lowerbound); 1205 nodeinfofirst = (NODEINFO*)SCIPpqueueFirst(ssg->subtreepqueues[subtreeidx]); 1206 assert(nodeinfofirst == NULL || subtreeidx == nodeinfofirst->subtreeidx); 1207 newgap = calcGap(scip, nodeinfofirst != NULL ? nodeinfofirst->lowerbound : SCIPinfinity(scip) ); 1208 1209 assert(SCIPisLE(scip, newgap, oldgap)); 1210 1211 /* the SSG value is always up-to-date because it is recomputed when the primal bound changes */ 1212 ssg->value += ssg->scalingfactor * MIN(newgap - oldgap, 0.0); 1213 } 1214 1215 SCIP_CALL( SCIPhashmapRemove(ssg->nodes2info, (void*)node) ); 1216 1217 SCIPdebugMsg(scip, "Removed node %" SCIP_LONGINT_FORMAT " from open nodes of SSG\n", 1218 SCIPnodeGetNumber(node)); 1219 1220 SCIPfreeBlockMemory(scip, &nodeinfo); 1221 1222 return SCIP_OKAY; 1223 } 1224 1225 /** insert children into subtree sum gap */ 1226 static 1227 SCIP_RETCODE subtreeSumGapInsertChildren( 1228 SCIP* scip, /**< SCIP data structure */ 1229 SUBTREESUMGAP* ssg /**< subtree sum gap data structure */ 1230 ) 1231 { 1232 int nchildren; 1233 SCIP_NODE** children; 1234 SCIP_NODE* focusnode; 1235 SCIP_NODE* parentnode; 1236 NODEINFO* parentnodeinfo; 1237 int parentnodelabel; 1238 int n; 1239 1240 assert(scip != NULL); 1241 assert(ssg != NULL); 1242 1243 if( ssg->nsubtrees == 1 ) 1244 return SCIP_OKAY; 1245 1246 SCIP_CALL( SCIPgetChildren(scip, &children, &nchildren) ); 1247 1248 if( nchildren == 0 ) 1249 return SCIP_OKAY; 1250 1251 focusnode = SCIPgetFocusNode(scip); 1252 1253 /* a rare case: the search has been stopped at some point, and the current focus node is the only descendant 1254 * of its parent node 1255 */ 1256 if( !SCIPhashmapExists(ssg->nodes2info, (void*)focusnode) ) 1257 { 1258 parentnode = focusnode; 1259 do 1260 { 1261 parentnode = SCIPnodeGetParent(parentnode); 1262 } while( parentnode != NULL && !SCIPhashmapExists(ssg->nodes2info, (void *)parentnode)); 1263 1264 assert(parentnode != NULL && SCIPhashmapExists(ssg->nodes2info, (void *)parentnode)); 1265 } 1266 else 1267 parentnode = focusnode; 1268 1269 parentnodeinfo = (NODEINFO*)SCIPhashmapGetImage(ssg->nodes2info, (void *)parentnode); 1270 parentnodelabel = parentnodeinfo->subtreeidx; 1271 1272 /* loop over children and insert the focus node label */ 1273 for( n = 0; n < nchildren; ++n ) 1274 { 1275 assert(SCIPnodeGetParent(children[n]) == focusnode); 1276 1277 SCIPdebugMsg(scip, 1278 "Inserting label %d for node number %" SCIP_LONGINT_FORMAT " (parent %" SCIP_LONGINT_FORMAT ")\n", 1279 parentnodelabel, SCIPnodeGetNumber(children[n]), SCIPnodeGetNumber(parentnode)); 1280 1281 SCIP_CALL( subtreeSumGapStoreNode(scip, ssg, children[n], parentnodelabel) ); 1282 } 1283 1284 /* remove focus node from hash map */ 1285 SCIP_CALL( subtreeSumGapRemoveNode(scip, ssg, parentnode) ); 1286 1287 return SCIP_OKAY; 1288 } 1289 1290 /* this function is inefficient because it loops over all open nodes, but can be used for debugging */ 1291 #ifdef SCIP_DISABLED_CODE 1292 /** compute subtree sum gap from scratch (inefficiently because loop over all open nodes) */ 1293 static 1294 SCIP_RETCODE subtreesumgapComputeFromScratch( 1295 SCIP* scip, /**< SCIP data structure */ 1296 SUBTREESUMGAP* ssg, /**< subtree sum gap data structure */ 1297 SCIP_Bool updatescaling /**< should the scaling factor be updated? */ 1298 ) 1299 { 1300 SCIP_Real* lowerbounds; 1301 SCIP_NODE** opennodes[3]; 1302 SCIP_Real gapsum = 0; 1303 SCIP_Real pb; 1304 int nopennodes[3]; 1305 int l; 1306 int t; 1307 1308 /* treat trivial cases: only 1 subtree, no incumbent solution */ 1309 if( SCIPisInfinity(scip, SCIPgetUpperbound(scip)) ) 1310 { 1311 ssg->value = 1.0; 1312 1313 return SCIP_OKAY; 1314 } 1315 1316 /* simply use normal gap in trivial case */ 1317 if( ssg->nsubtrees == 1 ) 1318 { 1319 ssg->value = calcGap(scip, SCIPgetLowerbound(scip)); 1320 1321 return SCIP_OKAY; 1322 } 1323 1324 /* allocate temporary memory to store lower bound for every subtree */ 1325 SCIP_CALL( SCIPallocBufferArray(scip, &lowerbounds, ssg->nsubtrees) ); 1326 1327 /* initialize lower bounds as SCIPinfinity(scip) */ 1328 for( l = 0; l < ssg->nsubtrees; ++l ) 1329 lowerbounds[l] = SCIPinfinity(scip); 1330 1331 /* loop over children, siblings, and leaves to update subtree lower bounds */ 1332 SCIP_CALL( SCIPgetOpenNodesData(scip, &opennodes[0], &opennodes[1], &opennodes[2], &nopennodes[0], &nopennodes[1], &nopennodes[2]) ); 1333 1334 /* loop over the three types leaves, siblings, leaves */ 1335 for( t = 0; t < 3; ++t ) 1336 { 1337 int n; 1338 /* loop over nodes of this type */ 1339 for( n = 0; n < nopennodes[t]; ++n ) 1340 { 1341 SCIP_NODE* node = opennodes[t][n]; 1342 NODEINFO* nodeinfo; 1343 SCIP_Real lowerbound; 1344 int label; 1345 nodeinfo = (NODEINFO*)SCIPhashmapGetImage(ssg->nodes2info, (void *)node); 1346 label = nodeinfo->subtreeidx; 1347 lowerbound = nodeinfo->lowerbound; 1348 1349 assert(label >= 0 && label < ssg->nsubtrees); 1350 lowerbounds[label] = MIN(lowerbounds[label], lowerbound); 1351 } 1352 } 1353 1354 /* compute subtree gaps in original space; sum them up */ 1355 pb = SCIPgetPrimalbound(scip); 1356 for( l = 0; l < ssg->nsubtrees; ++l ) 1357 { 1358 SCIP_Real subtreedualbound; 1359 SCIP_Real subtreegap; 1360 /* skip subtrees with infinite lower bound; they are empty and contribute 0.0 to the gap sum term */ 1361 if( SCIPisInfinity(scip, lowerbounds[l]) ) 1362 continue; 1363 1364 subtreedualbound = SCIPretransformObj(scip, lowerbounds[l]); 1365 1366 if( SCIPisEQ(scip, subtreedualbound, pb) ) 1367 continue; 1368 1369 subtreegap = REALABS(pb - subtreedualbound)/MAX(REALABS(pb),REALABS(subtreedualbound)); 1370 subtreegap = MIN(subtreegap, 1.0); 1371 1372 gapsum += subtreegap; 1373 } 1374 1375 /* update the scaling factor by using the previous SSG value divided by the current gapsum */ 1376 if( updatescaling ) 1377 { 1378 ssg->scalingfactor = ssg->value / MAX(gapsum, 1e-6); 1379 } 1380 1381 /* update and store SSG value by considering scaling factor */ 1382 ssg->value = ssg->scalingfactor * gapsum; 1383 1384 SCIPfreeBufferArray(scip, &lowerbounds); 1385 1386 return SCIP_OKAY; 1387 } 1388 #endif 1389 1390 /** compute subtree sum gap from scratch efficiently (linear effort in the number of subtrees) */ 1391 static 1392 SCIP_RETCODE subtreeSumGapComputeFromScratchEfficiently( 1393 SCIP* scip, /**< SCIP data structure */ 1394 SUBTREESUMGAP* ssg, /**< subtree sum gap data structure */ 1395 SCIP_Bool updatescaling /**< should the scaling factor be updated? */ 1396 ) 1397 { 1398 SCIP_Real gapsum = 0.0; 1399 int l; 1400 1401 /* treat trivial cases: only 1 subtree, no incumbent solution */ 1402 if( SCIPisInfinity(scip, SCIPgetUpperbound(scip)) ) 1403 { 1404 ssg->value = 1.0; 1405 1406 return SCIP_OKAY; 1407 } 1408 1409 if( ssg->nsubtrees == 1 ) 1410 { 1411 ssg->value = calcGap(scip, SCIPgetLowerbound(scip)); 1412 1413 return SCIP_OKAY; 1414 } 1415 1416 /* compute subtree gaps in original space; sum them up */ 1417 for( l = 0; l < ssg->nsubtrees; ++l ) 1418 { 1419 SCIP_Real subtreegap; 1420 NODEINFO* nodeinfo; 1421 1422 assert(ssg->subtreepqueues[l] != NULL); 1423 1424 nodeinfo = (NODEINFO*)SCIPpqueueFirst(ssg->subtreepqueues[l]); 1425 1426 /* skip subtrees with infinite lower bound; they are empty and contribute 0.0 to the gap sum term */ 1427 if( nodeinfo == NULL || SCIPisInfinity(scip, nodeinfo->lowerbound) ) 1428 continue; 1429 1430 subtreegap = calcGap(scip, nodeinfo->lowerbound); 1431 1432 gapsum += subtreegap; 1433 } 1434 1435 /* update the scaling factor by using the previous SSG value divided by the current gapsum */ 1436 if( updatescaling ) 1437 { 1438 ssg->scalingfactor = ssg->value / MAX(gapsum, 1e-6); 1439 } 1440 1441 /* update and store SSG value by considering scaling factor */ 1442 ssg->value = ssg->scalingfactor * gapsum; 1443 1444 return SCIP_OKAY; 1445 } 1446 1447 /** update the subtree sum gap after a node event (branching or deletion of a node) */ 1448 static 1449 SCIP_RETCODE subtreeSumGapUpdate( 1450 SCIP* scip, /**< SCIP data structure */ 1451 SUBTREESUMGAP* ssg, /**< subtree sum gap data structure */ 1452 SCIP_NODE* node, /**< the corresponding node */ 1453 int nchildren, /**< number of children */ 1454 SCIP_Longint nsolvednodes /**< number of solved nodes so far, used as a time stamp */ 1455 ) 1456 { 1457 SCIP_Bool updatescaling = FALSE; 1458 SCIP_Bool insertchildren = (ssg->nsubtrees > 1 && nchildren > 0); 1459 1460 /* if the instance is solved or a node is cutoff at the initsolve stage or we are unbounded, the ssg is 0 */ 1461 if( SCIPgetStage(scip) == SCIP_STAGE_SOLVED || SCIPgetStage(scip) == SCIP_STAGE_INITSOLVE || SCIPisInfinity(scip, -SCIPgetUpperbound(scip)) ) 1462 { 1463 ssg->value = 0.0; 1464 1465 return SCIP_OKAY; 1466 } 1467 1468 /* make a new tree split if the primal bound has changed. */ 1469 if( ! SCIPisInfinity(scip, SCIPgetUpperbound(scip)) && ! SCIPisEQ(scip, SCIPgetPrimalbound(scip), ssg->pblastsplit) ) 1470 { 1471 int nnewsubtrees; 1472 SCIP_Bool addfocusnode; 1473 1474 addfocusnode = SCIPgetFocusNode(scip) != NULL && SCIPgetNChildren(scip) == 0 && !SCIPwasNodeLastBranchParent(scip, SCIPgetFocusNode(scip)); 1475 nnewsubtrees = SCIPgetNSiblings(scip) + SCIPgetNLeaves(scip) + SCIPgetNChildren(scip) + (addfocusnode ? 1 : 0); 1476 1477 /* check if number of new subtrees does not exceed maximum number of subtrees; always split if no split happened, yet */ 1478 if( ssg->nsubtrees <= 1 || 1479 ((ssg->nmaxsubtrees == -1 || nnewsubtrees <= ssg->nmaxsubtrees) && 1480 (nsolvednodes - ssg->nodelastsplit >= ssg->nminnodeslastsplit)) ) 1481 { 1482 SCIP_CALL( subtreeSumGapSplit(scip, ssg, addfocusnode) ); 1483 1484 /* remember time stamp */ 1485 ssg->nodelastsplit = nsolvednodes; 1486 } 1487 else 1488 { 1489 if( ssg->nmaxsubtrees != -1 && nnewsubtrees >= ssg->nmaxsubtrees ) 1490 { 1491 SCIPdebugMsg(scip, "Keep split into %d subtrees because new split into %d subtrees exceeds limit %d\n", 1492 ssg->nsubtrees, nnewsubtrees, ssg->nmaxsubtrees); 1493 } 1494 else 1495 { 1496 SCIPdebugMsg(scip, "Keep split into %d subtrees from %" SCIP_LONGINT_FORMAT " nodes ago\n", 1497 ssg->nsubtrees, nsolvednodes - ssg->nodelastsplit); 1498 } 1499 1500 /* no new split has happened; insert the new children to their SSG subtree */ 1501 if( insertchildren ) 1502 { 1503 SCIP_CALL( subtreeSumGapInsertChildren(scip, ssg) ); 1504 } 1505 } 1506 1507 ssg->pblastsplit = SCIPgetPrimalbound(scip); 1508 1509 updatescaling = TRUE; 1510 1511 /* compute the current SSG value from scratch */ 1512 SCIP_CALL( subtreeSumGapComputeFromScratchEfficiently(scip, ssg, updatescaling) ); 1513 } 1514 /* otherwise, if new children have been created, label them */ 1515 else if( insertchildren ) 1516 { 1517 SCIP_CALL( subtreeSumGapInsertChildren(scip, ssg) ); 1518 } 1519 1520 /* remove the node from the hash map if it is a leaf */ 1521 if( nchildren == 0 ) 1522 { 1523 SCIP_CALL( subtreeSumGapRemoveNode(scip, ssg, node) ); 1524 } 1525 1526 return SCIP_OKAY; 1527 } 1528 1529 /** reset tree data */ 1530 static 1531 SCIP_RETCODE resetTreeData( 1532 SCIP* scip, /**< SCIP data structure */ 1533 TREEDATA* treedata /**< tree data */ 1534 ) 1535 { 1536 /* simply set everything to 0 */ 1537 treedata->ninner = treedata->nleaves = treedata->nvisited = 0L; 1538 treedata->weight = 0.0; 1539 1540 /* set up root node */ 1541 treedata->nnodes = 1; 1542 treedata->nopen = 1; 1543 1544 SCIP_CALL( subtreeSumGapReset(scip, treedata->ssg) ); 1545 1546 return SCIP_OKAY; 1547 } 1548 1549 /** create tree data structure */ 1550 static 1551 SCIP_RETCODE createTreeData( 1552 SCIP* scip, /**< SCIP data structure */ 1553 TREEDATA** treedata /**< pointer to store tree data */ 1554 ) 1555 { 1556 assert(treedata != NULL); 1557 assert(scip != NULL); 1558 1559 SCIP_CALL( SCIPallocMemory(scip, treedata) ); 1560 1561 SCIP_CALL( subtreeSumGapCreate(scip, &(*treedata)->ssg) ); 1562 1563 SCIP_CALL( resetTreeData(scip, *treedata) ); 1564 1565 return SCIP_OKAY; 1566 } 1567 1568 /** free tree data structure */ 1569 static 1570 void freeTreeData( 1571 SCIP* scip, /**< SCIP data structure */ 1572 TREEDATA** treedata /**< pointer to tree data */ 1573 ) 1574 { 1575 assert(scip != NULL); 1576 1577 if( *treedata == NULL ) 1578 return; 1579 1580 subtreeSumGapFree(scip, &(*treedata)->ssg); 1581 1582 SCIPfreeMemory(scip, treedata); 1583 *treedata = NULL; 1584 } 1585 1586 /** update tree data structure after a node has been solved/is about to be deleted */ 1587 static 1588 SCIP_RETCODE updateTreeData( 1589 SCIP* scip, /**< SCIP data structure */ 1590 TREEDATA* treedata, /**< tree data */ 1591 SCIP_NODE* node, /**< the corresponding node */ 1592 int nchildren /**< the number of children */ 1593 ) 1594 { 1595 assert(node != NULL); 1596 1597 ++treedata->nvisited; 1598 treedata->nopen--; 1599 1600 if( nchildren == 0 ) 1601 { 1602 int depth = SCIPnodeGetDepth(node); 1603 treedata->nleaves++; 1604 treedata->weight += pow(0.5, (SCIP_Real)depth); 1605 } 1606 else 1607 { 1608 treedata->nnodes += nchildren; 1609 treedata->nopen += nchildren; 1610 ++treedata->ninner; 1611 } 1612 1613 /* update the subtree sum gap */ 1614 if( ! SCIPisInRestart(scip) ) 1615 { 1616 SCIP_CALL( subtreeSumGapUpdate(scip, treedata->ssg, node, nchildren, treedata->nvisited) ); 1617 } 1618 1619 return SCIP_OKAY; 1620 } 1621 1622 /** get weighted backtrack estimation from this tree data */ 1623 static 1624 SCIP_Real treeDataGetWbe( 1625 TREEDATA* treedata /**< tree data */ 1626 ) 1627 { 1628 if( treedata->weight <= 0.0 || treedata->nleaves == 0 ) 1629 return -1.0; 1630 1631 return 2.0 * treedata->nleaves / (SCIP_Real)treedata->weight - 1.0; 1632 } 1633 1634 #ifdef SCIP_DEBUG 1635 /* print method for tree data */ 1636 static 1637 char* treeDataPrint( 1638 TREEDATA* treedata, /**< tree data */ 1639 char* strbuf /**< string buffer */ 1640 ) 1641 { 1642 (void )SCIPsnprintf(strbuf, SCIP_MAXSTRLEN, 1643 "Tree Data: %" SCIP_LONGINT_FORMAT " nodes (" 1644 "%" SCIP_LONGINT_FORMAT " visited, " 1645 "%" SCIP_LONGINT_FORMAT " inner, " 1646 "%" SCIP_LONGINT_FORMAT " leaves, " 1647 "%" SCIP_LONGINT_FORMAT " open), " 1648 "weight: %.4Lf, ssg %.4f", 1649 treedata->nnodes, 1650 treedata->nvisited, 1651 treedata->ninner, 1652 treedata->nleaves, 1653 treedata->nopen, 1654 treedata->weight, 1655 treedata->ssg->value 1656 ); 1657 return strbuf; 1658 } 1659 #endif 1660 1661 /** reset double exponential smoothing */ 1662 static 1663 void doubleExpSmoothReset( 1664 DOUBLEEXPSMOOTH* des, /**< double exponential smoothing data structure */ 1665 SCIP_Real initialvalue /**< the initial value */ 1666 ) 1667 { 1668 des->n = 0; 1669 des->level = SCIP_INVALID; 1670 des->trend = SCIP_INVALID; 1671 des->initialvalue = initialvalue; 1672 } 1673 1674 /** initialize a double exponential smoothing data structure */ 1675 static 1676 void doubleExpSmoothInit( 1677 DOUBLEEXPSMOOTH* des, /**< double exponential smoothing data structure */ 1678 SCIP_Real x1 /**< the first sample value */ 1679 ) 1680 { 1681 assert(des != NULL); 1682 1683 des->n = 1; 1684 des->level = x1; 1685 des->trend = x1 - des->initialvalue; 1686 1687 des->usetrendinlevel = DES_USETRENDINLEVEL; 1688 1689 return; 1690 } 1691 1692 /** update a double exponential smoothing data structure */ 1693 static 1694 void doubleExpSmoothUpdate( 1695 DOUBLEEXPSMOOTH* des, /**< double exponential smoothing data structure */ 1696 SCIP_Real xnew /**< new sample value */ 1697 ) 1698 { 1699 if( des->n == 0 ) 1700 doubleExpSmoothInit(des, xnew); 1701 else 1702 { 1703 SCIP_Real newlevel; 1704 SCIP_Real newtrend; 1705 1706 newlevel = des->alpha * xnew + (1.0 - des->alpha) * (des->level + (des->usetrendinlevel ? des->trend : 0.0)); 1707 newtrend = des->beta * (newlevel - des->level) + (1.0 - des->beta) * des->trend; 1708 1709 des->level = newlevel; 1710 des->trend = newtrend; 1711 } 1712 } 1713 1714 /** get the current trend (slope) computed by this double exponential smoothing */ 1715 static 1716 SCIP_Real doubleExpSmoothGetTrend( 1717 DOUBLEEXPSMOOTH* des /**< double exponential smoothing data structure */ 1718 ) 1719 { 1720 assert(des != NULL); 1721 1722 if( des->n == 0 ) 1723 return SCIP_INVALID; 1724 1725 return des->trend; 1726 } 1727 1728 /** reset time series */ 1729 static 1730 void timeSeriesReset( 1731 TIMESERIES* timeseries /**< pointer to store time series */ 1732 ) 1733 { 1734 timeseries->resolution = 1; 1735 timeseries->nvals = 0; 1736 timeseries->nobs = 0L; 1737 timeseries->currentvalue = timeseries->initialvalue; 1738 timeseries->smoothestimation = SCIP_INVALID; 1739 1740 doubleExpSmoothReset(×eries->des, timeseries->initialvalue); 1741 } 1742 1743 /** create a time series object */ 1744 static 1745 SCIP_RETCODE timeSeriesCreate( 1746 SCIP* scip, /**< SCIP data structure */ 1747 TIMESERIES** timeseries, /**< pointer to store time series */ 1748 const char* name, /**< name of this time series */ 1749 SCIP_Real targetvalue, /**< target value of this time series */ 1750 SCIP_Real initialvalue, /**< the initial value of time series */ 1751 SCIP_Real alpha, /**< alpha parameter (level weight) for double exponential smoothing */ 1752 SCIP_Real beta, /**< beta parameter (level weight) for double exponential smoothing */ 1753 DECL_TIMESERIESUPDATE ((*timeseriesupdate)) /**< update callback at nodes, or NULL */ 1754 ) 1755 { 1756 TIMESERIES* timeseriesptr; 1757 assert(scip != NULL); 1758 assert(timeseries != NULL); 1759 assert(name != NULL); 1760 assert(alpha >= 0.0 && alpha <= 1); 1761 assert(beta >= 0.0 && beta <= 1); 1762 1763 SCIP_CALL( SCIPallocMemory(scip, timeseries) ); 1764 1765 timeseriesptr = *timeseries; 1766 assert(timeseriesptr != NULL); 1767 1768 /* copy name */ 1769 SCIP_ALLOC( BMSduplicateMemoryArray(×eriesptr->name, name, strlen(name)+1) ); 1770 1771 /* copy callbacks */ 1772 assert(timeseriesupdate != NULL); 1773 timeseriesptr->timeseriesupdate = timeseriesupdate; 1774 1775 timeseriesptr->targetvalue = targetvalue; 1776 timeseriesptr->valssize = 1024; 1777 timeseriesptr->initialvalue = initialvalue; 1778 1779 SCIP_CALL( SCIPallocMemoryArray(scip, ×eriesptr->vals, timeseriesptr->valssize) ); 1780 SCIP_CALL( SCIPallocMemoryArray(scip, ×eriesptr->estimation, timeseriesptr->valssize) ); 1781 1782 timeSeriesReset(timeseriesptr); 1783 1784 timeseriesptr->des.alpha = alpha; 1785 timeseriesptr->des.beta = beta; 1786 1787 SCIPdebugMsg(scip, "Finished creation of time series '%s'\n", timeseriesptr->name); 1788 1789 return SCIP_OKAY; 1790 } 1791 1792 /** free a time series */ 1793 static 1794 void timeSeriesFree( 1795 SCIP* scip, /**< SCIP data structure */ 1796 TIMESERIES** timeseries /**< pointer to time series */ 1797 ) 1798 { 1799 assert(scip != NULL); 1800 assert(timeseries != NULL); 1801 1802 BMSfreeMemoryArray(&(*timeseries)->name); 1803 1804 SCIPfreeMemoryArray(scip, &(*timeseries)->vals); 1805 SCIPfreeMemoryArray(scip, &(*timeseries)->estimation); 1806 1807 SCIPfreeMemory(scip, timeseries); 1808 1809 *timeseries = NULL; 1810 } 1811 1812 /** get current value of time series */ 1813 static 1814 SCIP_Real timeSeriesGetValue( 1815 TIMESERIES* timeseries /**< time series */ 1816 ) 1817 { 1818 assert(timeseries != NULL); 1819 1820 return timeseries->currentvalue; 1821 } 1822 1823 /** get target value (which this time series reaches at the end of the solution process) */ 1824 static 1825 SCIP_Real timeSeriesGetTargetValue( 1826 TIMESERIES* timeseries /**< time series */ 1827 ) 1828 { 1829 return timeseries->targetvalue; 1830 } 1831 1832 /** get resolution of time series */ 1833 static 1834 int timeSeriesGetResolution( 1835 TIMESERIES* timeseries /**< time series */ 1836 ) 1837 { 1838 return timeseries->resolution; 1839 } 1840 1841 /** estimate tree size at which time series reaches target value */ 1842 static 1843 SCIP_Real timeSeriesEstimate( 1844 TIMESERIES* timeseries, /**< time series */ 1845 TREEDATA* treedata /**< tree data for fallback estimation */ 1846 ) 1847 { 1848 SCIP_Real val; 1849 SCIP_Real targetval; 1850 SCIP_Real trend; 1851 SCIP_Real estimated; 1852 const SCIP_Real tolerance = 1e-6; 1853 1854 /* if no observations have been made yet, return infinity */ 1855 if( timeseries->nobs == 0L ) 1856 return -1.0; 1857 1858 val = timeSeriesGetValue(timeseries); 1859 targetval = timeSeriesGetTargetValue(timeseries); 1860 1861 /* if the value has reached the target value already, return the number of observations */ 1862 if( EPSZ(val - targetval, tolerance) ) 1863 return treedata->nnodes; 1864 1865 trend = doubleExpSmoothGetTrend(×eries->des); 1866 1867 /* Get current value and trend. The linear trend estimation may point into the wrong direction 1868 * In this case, we use the fallback mechanism that we will need twice as many nodes. 1869 */ 1870 if( (targetval > val && trend < tolerance) || (targetval < val && trend > -tolerance) ) 1871 { 1872 return 2.0 * treedata->nvisited; 1873 } 1874 1875 /* compute after how many additional steps the current trend reaches the target value; multiply by resolution */ 1876 estimated = timeSeriesGetResolution(timeseries) * (timeseries->nvals + (targetval - val) / (SCIP_Real)trend); 1877 return timeseries->useleafts ? 2.0 * estimated - 1.0 : estimated; 1878 } 1879 1880 /** update time series smoothened estimation */ 1881 static 1882 void timeSeriesUpdateSmoothEstimation( 1883 TIMESERIES* timeseries, /**< time series */ 1884 SCIP_Real estimation /**< estimation value */ 1885 ) 1886 { 1887 if( timeseries->smoothestimation == SCIP_INVALID )/*lint !e777*/ 1888 timeseries->smoothestimation = estimation; 1889 else 1890 { 1891 timeseries->smoothestimation *= (1.0 - SESCOEFF); 1892 timeseries->smoothestimation += SESCOEFF * estimation; 1893 } 1894 } 1895 1896 /** get smooth estimation of time series */ 1897 static 1898 SCIP_Real timeSeriesGetSmoothEstimation( 1899 TIMESERIES* timeseries /**< time series */ 1900 ) 1901 { 1902 return timeseries->smoothestimation; 1903 } 1904 1905 /** resample to lower resolution */ 1906 static 1907 void timeSeriesResample( 1908 TIMESERIES* timeseries /**< time series */ 1909 ) 1910 { 1911 DOUBLEEXPSMOOTH* des; 1912 int i; 1913 1914 assert(timeseries->nvals % 2 == 0); 1915 1916 des = ×eries->des; 1917 doubleExpSmoothReset(des, timeseries->initialvalue); 1918 1919 /* compress vals array to store only every second entry */ 1920 for( i = 0; i < timeseries->nvals / 2; ++i ) 1921 { 1922 timeseries->vals[i] = timeseries->vals[2 * i]; 1923 timeseries->estimation[i] = timeseries->estimation[2 * i]; 1924 doubleExpSmoothUpdate(des, timeseries->vals[i]); 1925 timeSeriesUpdateSmoothEstimation(timeseries, timeseries->estimation[i]); 1926 } 1927 1928 timeseries->resolution *= 2; 1929 timeseries->nvals = timeseries->nvals / 2; 1930 } 1931 1932 /** update time series */ 1933 static 1934 SCIP_RETCODE timeSeriesUpdate( 1935 SCIP* scip, /**< SCIP data structure */ 1936 TIMESERIES* timeseries, /**< time series */ 1937 TREEDATA* treedata, /**< tree data */ 1938 SCIP_Bool isleaf /**< are we at a leaf node? */ 1939 ) 1940 { 1941 SCIP_Real value; 1942 1943 assert(scip != NULL); 1944 assert(timeseries != NULL); 1945 assert(treedata != NULL); 1946 1947 /* call update callback */ 1948 assert(timeseries->timeseriesupdate != NULL); 1949 SCIP_CALL( timeseries->timeseriesupdate(scip, timeseries, treedata, &value) ); 1950 1951 /* store the value as current value */ 1952 timeseries->currentvalue = value; 1953 1954 if( timeseries->useleafts && ! isleaf ) 1955 return SCIP_OKAY; 1956 1957 timeseries->nobs++; 1958 1959 /* if this is a leaf that matches the time series resolution, store the value */ 1960 if( timeseries->nobs % timeseries->resolution == 0 ) 1961 { 1962 int tspos; 1963 SCIP_Real estimate; 1964 1965 assert(timeseries->nvals < timeseries->valssize); 1966 tspos = timeseries->nvals++; 1967 timeseries->vals[tspos] = value; 1968 doubleExpSmoothUpdate(×eries->des, value); 1969 estimate = timeSeriesEstimate(timeseries, treedata); 1970 timeseries->estimation[tspos] = estimate; 1971 timeSeriesUpdateSmoothEstimation(timeseries, estimate); 1972 } 1973 1974 /* if the time series has reached its capacity, resample and increase the resolution */ 1975 if( timeseries->nvals == timeseries->valssize ) 1976 timeSeriesResample(timeseries); 1977 1978 return SCIP_OKAY; 1979 } 1980 1981 /** get name of time series */ 1982 static 1983 char* timeSeriesGetName( 1984 TIMESERIES* timeseries /**< time series */ 1985 ) 1986 { 1987 return timeseries->name; 1988 } 1989 1990 /** reset all time series */ 1991 static 1992 void resetTimeSeries( 1993 SCIP_EVENTHDLRDATA* eventhdlrdata /**< event handler data */ 1994 ) 1995 { 1996 TIMESERIES** tss = eventhdlrdata->timeseries; 1997 int t; 1998 1999 /* loop over time series and reset them */ 2000 for( t = 0; t < NTIMESERIES; ++t ) 2001 { 2002 assert(tss[t] != NULL); 2003 timeSeriesReset(tss[t]); 2004 2005 tss[t]->useleafts = eventhdlrdata->useleafts; 2006 } 2007 } 2008 2009 /* 2010 * Callback methods of event handler 2011 */ 2012 2013 /** free all time series */ 2014 static 2015 void freeTimeSeries( 2016 SCIP* scip, /**< SCIP data structure */ 2017 SCIP_EVENTHDLRDATA* eventhdlrdata /**< event handler data */ 2018 ) 2019 { 2020 TIMESERIES** tss = eventhdlrdata->timeseries; 2021 int t; 2022 2023 /* loop over time series and reset them */ 2024 for( t = 0; t < NTIMESERIES; ++t ) 2025 { 2026 assert(tss[t] != NULL); 2027 timeSeriesFree(scip, &tss[t]); 2028 } 2029 } 2030 2031 /** get ensemble tree size estimation as a combination of the individual time series estimations 2032 * 2033 * the coefficients have been computed based on a nonlinear fit on a broad set of publicly available 2034 * MIP instances; please refer to the publication at the top of this file for further details. 2035 */ 2036 static 2037 SCIP_Real getEnsembleEstimation( 2038 SCIP_EVENTHDLRDATA* eventhdlrdata /**< event handler data */ 2039 ) 2040 { 2041 TREEDATA* treedata; 2042 SCIP_Real* coeffs; 2043 SCIP_Real estim; 2044 int t; 2045 2046 TSPOS tsposs[] = { 2047 TSPOS_GAP, 2048 TSPOS_TREEWEIGHT, 2049 TSPOS_LFREQ, 2050 TSPOS_SSG, 2051 TSPOS_OPEN 2052 }; 2053 2054 /* coefficients for the early stage (tree weight <= 0.3) */ 2055 SCIP_Real coeffs_early[] = { 2056 0.002, /* gap */ 2057 0.381, /* tree weight */ 2058 0.469, /* leaf-frequency */ 2059 0.292, /* SSG */ 2060 0.004 /* open-nodes */ 2061 }; 2062 2063 /* coefficients for the intermediate stage (0.3 < tree weight <= 0.6) */ 2064 SCIP_Real coeffs_intermediate[] = { 2065 0.011, /* gap */ 2066 0.193, /* tree weight */ 2067 0.351, /* leaf-frequency */ 2068 0.012, /* SSG */ 2069 0.051 /* open-nodes */ 2070 }; 2071 2072 /* coefficients for the late stage (tree weight > 0.6) */ 2073 SCIP_Real coeffs_late[] = { 2074 0.000, /* gap */ 2075 0.033, /* tree weight */ 2076 0.282, /* leaf-frequency */ 2077 0.003, /* SSG */ 2078 0.024 /* open-nodes */ 2079 }; 2080 2081 assert(eventhdlrdata != NULL); 2082 treedata = eventhdlrdata->treedata; 2083 2084 /* assign coeffs based on stage */ 2085 if( treedata->weight <= 0.3 ) 2086 { 2087 estim = 0.0; 2088 coeffs = coeffs_early; 2089 /* ensure that coeffs and time series are still aligned */ 2090 assert(sizeof(coeffs_early)/sizeof(SCIP_Real) == NTIMESERIES); /*lint !e506*/ 2091 } 2092 else if( treedata->weight <= 0.6 ) 2093 { 2094 coeffs = coeffs_intermediate; 2095 /* ensure that coeffs and time series are still aligned */ 2096 assert(sizeof(coeffs_intermediate)/sizeof(SCIP_Real) == NTIMESERIES); /*lint !e506*/ 2097 2098 /* initialize by intermediate WBE coefficient */ 2099 estim = 0.156 * treeDataGetWbe(treedata); 2100 } 2101 else 2102 { 2103 coeffs = coeffs_late; 2104 /* ensure that coeffs and time series are still aligned */ 2105 assert(sizeof(coeffs_late)/sizeof(SCIP_Real) == NTIMESERIES); /*lint !e506*/ 2106 2107 /* initialize by late WBE coefficient */ 2108 estim = 0.579 * treeDataGetWbe(treedata); 2109 } 2110 2111 /* combine estimation using the stage-dependent coefficients */ 2112 for( t = 0; t < NTIMESERIES; ++t ) 2113 { 2114 SCIP_Real testim; 2115 TSPOS tspos = tsposs[t]; 2116 testim = timeSeriesEstimate(eventhdlrdata->timeseries[tspos], treedata); 2117 2118 if( testim < 0.0 ) 2119 testim = treedata->nnodes; 2120 2121 estim += coeffs[t] * testim; 2122 } 2123 2124 if( estim < treedata->nnodes ) 2125 return (SCIP_Real)treedata->nnodes; 2126 else 2127 return estim; 2128 } 2129 2130 /** get approximation of search tree completion depending on the selected method */ 2131 static 2132 SCIP_RETCODE getSearchCompletion( 2133 SCIP_EVENTHDLRDATA* eventhdlrdata, /**< event handler data */ 2134 SCIP_Real* completed /**< pointer to store the search tree completion */ 2135 ) 2136 { 2137 SCIP_Real values[9]; 2138 TREEDATA* treedata; 2139 char completiontype; 2140 2141 assert(eventhdlrdata != NULL); 2142 treedata = eventhdlrdata->treedata; 2143 completiontype = eventhdlrdata->completiontypeparam; 2144 2145 /* infer automatic completion type 2146 * 2147 * use regression forest if available, 2148 * or 2149 * use monotone regression if both SSG and tree weight are meaningful; 2150 * or 2151 * use tree weight or SSG, depending which one is available, 2152 * or 2153 * use gap, which is always available 2154 */ 2155 if( completiontype == COMPLETIONTYPE_AUTO ) 2156 { 2157 SCIP_Bool useweight = eventhdlrdata->treeisbinary; 2158 SCIP_Bool usessg = treedata->ssg->pblastsplit != SSG_STARTPRIMBOUND;/*lint !e777*/ 2159 2160 if( eventhdlrdata->regforest != NULL ) 2161 completiontype = COMPLETIONTYPE_REGFOREST; 2162 else if( useweight && usessg ) 2163 completiontype = COMPLETIONTYPE_MONOREG; 2164 else if( useweight ) 2165 completiontype = COMPLETIONTYPE_TREEWEIGHT; 2166 else if( usessg ) 2167 completiontype = COMPLETIONTYPE_SSG; 2168 else 2169 completiontype = COMPLETIONTYPE_GAP; 2170 } 2171 2172 /* compute the search tree completion based on the selected method */ 2173 switch (completiontype) 2174 { 2175 /* use regression forest */ 2176 case COMPLETIONTYPE_REGFOREST: 2177 values[0] = timeSeriesGetValue(eventhdlrdata->timeseries[TSPOS_TREEWEIGHT]); 2178 values[1] = doubleExpSmoothGetTrend(&eventhdlrdata->timeseries[TSPOS_TREEWEIGHT]->des); 2179 values[2] = timeSeriesGetValue(eventhdlrdata->timeseries[TSPOS_SSG]); 2180 values[3] = doubleExpSmoothGetTrend(&eventhdlrdata->timeseries[TSPOS_SSG]->des); 2181 values[4] = timeSeriesGetValue(eventhdlrdata->timeseries[TSPOS_LFREQ]); 2182 values[5] = doubleExpSmoothGetTrend(&eventhdlrdata->timeseries[TSPOS_LFREQ]->des); 2183 values[6] = timeSeriesGetValue(eventhdlrdata->timeseries[TSPOS_GAP]); 2184 values[7] = doubleExpSmoothGetTrend(&eventhdlrdata->timeseries[TSPOS_GAP]->des); 2185 values[8] = doubleExpSmoothGetTrend(&eventhdlrdata->timeseries[TSPOS_OPEN]->des) < 0 ? 1.0 : 0.0; 2186 2187 *completed = SCIPregForestPredict(eventhdlrdata->regforest, values); 2188 break; 2189 2190 /* interpolate between ssg and tree weight */ 2191 case COMPLETIONTYPE_MONOREG: 2192 *completed = eventhdlrdata->coefmonoweight * (SCIP_Real)treedata->weight + 2193 eventhdlrdata->coefmonossg * (1.0 - treedata->ssg->value); 2194 break; 2195 2196 case COMPLETIONTYPE_TREEWEIGHT: 2197 *completed = (SCIP_Real)treedata->weight; 2198 break; 2199 2200 case COMPLETIONTYPE_GAP: 2201 *completed = timeSeriesGetValue(eventhdlrdata->timeseries[TSPOS_GAP]); /* gap is stored as 1 - gap */ 2202 break; 2203 2204 case COMPLETIONTYPE_SSG: 2205 *completed = 1.0 - treedata->ssg->value; /* ssg is decreasing */ 2206 break; 2207 2208 default: 2209 SCIPerrorMessage("Unsupported completion type '%c'\n", completiontype); 2210 SCIPABORT(); 2211 return SCIP_PARAMETERWRONGVAL; 2212 } 2213 return SCIP_OKAY; 2214 } 2215 2216 /** tree size estimation based on search tree completion */ 2217 static 2218 SCIP_RETCODE getEstimCompletion( 2219 SCIP* scip, /**< SCIP data structure */ 2220 SCIP_EVENTHDLRDATA* eventhdlrdata, /**< event handler data */ 2221 SCIP_Real* estim /**< pointer to store the estimation value */ 2222 ) 2223 { 2224 SCIP_Real completed; 2225 2226 *estim = -1.0; 2227 2228 SCIP_CALL( getSearchCompletion(eventhdlrdata, &completed) ); 2229 2230 completed = MIN(completed, 1.0); 2231 2232 if( completed > 0.0 ) 2233 *estim = SCIPgetNNodes(scip) / completed; 2234 2235 return SCIP_OKAY; 2236 } 2237 2238 /** update callback at nodes */ 2239 static 2240 DECL_TIMESERIESUPDATE(timeseriesUpdateGap) 2241 { /*lint --e{715}*/ 2242 SCIP_Real primalbound; 2243 SCIP_Real dualbound; 2244 2245 assert(scip != NULL); 2246 assert(ts != NULL); 2247 assert(value != NULL); 2248 2249 /* avoid to call SCIPgetDualbound during a restart where the queue is simply emptied */ 2250 if( SCIPisInRestart(scip) ) 2251 { 2252 *value = timeSeriesGetValue(ts); 2253 2254 return SCIP_OKAY; 2255 } 2256 2257 primalbound = SCIPgetPrimalbound(scip); 2258 dualbound = SCIPgetDualbound(scip); 2259 if( SCIPisInfinity(scip, REALABS(primalbound)) || SCIPisInfinity(scip, REALABS(dualbound)) ) 2260 *value = 0; 2261 else if( SCIPisEQ(scip, primalbound, dualbound) ) 2262 *value = 1.0; 2263 else 2264 { 2265 SCIP_Real abspb; 2266 SCIP_Real absdb; 2267 2268 abspb = REALABS(primalbound); 2269 absdb = REALABS(dualbound); 2270 *value = 1.0 - REALABS(primalbound - dualbound)/MAX(abspb, absdb); 2271 } 2272 2273 /* using this max, we set the closed gap to 0 in the case where the primal and dual bound differ in their sign */ 2274 *value = MAX(*value, 0.0); 2275 2276 return SCIP_OKAY; 2277 } 2278 2279 /** update callback at nodes */ 2280 static 2281 DECL_TIMESERIESUPDATE(timeseriesUpdateTreeWeight) 2282 { /*lint --e{715}*/ 2283 *value = (SCIP_Real)treedata->weight; 2284 2285 return SCIP_OKAY; 2286 } 2287 2288 /** update callback at nodes */ 2289 static 2290 DECL_TIMESERIESUPDATE(timeseriesUpdateLeafFreq) 2291 { /*lint --e{715}*/ 2292 if( treedata->nvisited == 0 ) 2293 *value = -0.5; 2294 else 2295 *value = (treedata->nleaves - 0.5)/(SCIP_Real)treedata->nvisited; 2296 2297 return SCIP_OKAY; 2298 } 2299 2300 /** update callback at nodes */ 2301 static 2302 DECL_TIMESERIESUPDATE(timeseriesUpdateSsg) 2303 { /*lint --e{715}*/ 2304 if( treedata->nvisited == 0 ) 2305 *value = 1.0; 2306 else 2307 *value = treedata->ssg->value; 2308 2309 return SCIP_OKAY; 2310 } 2311 2312 /** update callback at nodes */ 2313 static 2314 DECL_TIMESERIESUPDATE(timeseriesUpdateOpenNodes) 2315 { /*lint --e{715}*/ 2316 if( treedata->nvisited == 0 ) 2317 *value = 0.0; 2318 else 2319 *value = (SCIP_Real)treedata->nopen; 2320 2321 return SCIP_OKAY; 2322 } 2323 2324 /** include time series to forecast into event handler */ 2325 static 2326 SCIP_RETCODE includeTimeseries( 2327 SCIP* scip, /**< SCIP data structure */ 2328 SCIP_EVENTHDLRDATA* eventhdlrdata /**< event handler data */ 2329 ) 2330 { 2331 assert(scip != NULL); 2332 assert(eventhdlrdata != NULL); 2333 2334 /* include gap time series */ 2335 SCIP_CALL( timeSeriesCreate(scip, &eventhdlrdata->timeseries[TSPOS_GAP], "gap", 1.0, 0.0, 2336 DES_ALPHA_GAP, DES_BETA_GAP, timeseriesUpdateGap) ); 2337 2338 /* include tree weight time series */ 2339 SCIP_CALL( timeSeriesCreate(scip, &eventhdlrdata->timeseries[TSPOS_TREEWEIGHT], "tree-weight", 1.0, 0.0, 2340 DES_ALPHA_TREEWEIGHT, DES_BETA_TREEWEIGHT, timeseriesUpdateTreeWeight) ); 2341 2342 /* include leaf time series */ 2343 SCIP_CALL( timeSeriesCreate(scip, &eventhdlrdata->timeseries[TSPOS_LFREQ], "leaf-frequency", 0.5, -0.5, 2344 DES_ALPHA_LEAFFREQUENCY, DES_BETA_LEAFFREQUENCY, timeseriesUpdateLeafFreq) ); 2345 2346 /* include SSG time series */ 2347 SCIP_CALL( timeSeriesCreate(scip, &eventhdlrdata->timeseries[TSPOS_SSG], "ssg", 0.0, 1.0, 2348 DES_ALPHA_SSG, DES_BETA_SSG, timeseriesUpdateSsg) ); 2349 2350 /* include open nodes time series */ 2351 SCIP_CALL( timeSeriesCreate(scip, &eventhdlrdata->timeseries[TSPOS_OPEN], "open-nodes", 0.0, 0.0, 2352 DES_ALPHA_OPENNODES, DES_BETA_OPENNODES, timeseriesUpdateOpenNodes) ); 2353 2354 return SCIP_OKAY; 2355 } 2356 2357 /** get restartpolicy based on the value of the restart parameter */ 2358 static 2359 RESTARTPOLICY getRestartPolicy( 2360 SCIP_EVENTHDLRDATA* eventhdlrdata /**< event handler data */ 2361 ) 2362 { 2363 switch (eventhdlrdata->restartpolicyparam) 2364 { 2365 case RESTARTPOLICY_CHAR_ALWAYS: 2366 return RESTARTPOLICY_ALWAYS; 2367 case RESTARTPOLICY_CHAR_NEVER: 2368 return RESTARTPOLICY_NEVER; 2369 case RESTARTPOLICY_CHAR_COMPLETION: 2370 return RESTARTPOLICY_COMPLETION; 2371 case RESTARTPOLICY_CHAR_ESTIMATION: 2372 return RESTARTPOLICY_ESTIMATION; 2373 default: 2374 SCIPerrorMessage("Unknown restart policy %c\n", eventhdlrdata->restartpolicyparam); 2375 break; 2376 } 2377 2378 return RESTARTPOLICY_NEVER; 2379 } 2380 2381 /** check if a restart is applicable considering limit and threshold user parameters */ 2382 static 2383 SCIP_Bool isRestartApplicable( 2384 SCIP* scip, /**< SCIP data structure */ 2385 SCIP_EVENTHDLRDATA* eventhdlrdata /**< event handler data */ 2386 ) 2387 { 2388 SCIP_Longint nnodes; 2389 2390 /* check whether to apply restarts when there are active pricers available */ 2391 if( SCIPgetNActivePricers(scip) > 0 && ! eventhdlrdata->restartactpricers ) 2392 return FALSE; 2393 2394 /* check whether to apply a restart when nonlinear constraints are present */ 2395 if( SCIPisNLPConstructed(scip) && ! eventhdlrdata->restartnonlinear ) 2396 return FALSE; 2397 2398 /* check if max number of restarts has been reached */ 2399 if( eventhdlrdata->restartlimit != -1 && eventhdlrdata->nrestartsperformed >= eventhdlrdata->restartlimit ) 2400 return FALSE; 2401 2402 /* check if number of nodes exceeds the minimum number of nodes */ 2403 if( eventhdlrdata->countonlyleaves ) 2404 nnodes = eventhdlrdata->treedata->nleaves; 2405 else 2406 nnodes = eventhdlrdata->treedata->nvisited; 2407 2408 if( nnodes < eventhdlrdata->minnodes ) 2409 return FALSE; 2410 2411 return TRUE; 2412 } 2413 2414 /** should a restart be applied based on the value of the selected completion method? */ /*lint --e{715}*/ 2415 static 2416 SCIP_Bool shouldApplyRestartCompletion( 2417 SCIP* scip, /**< SCIP data structure */ 2418 SCIP_EVENTHDLRDATA* eventhdlrdata /**< event handler data */ 2419 ) 2420 { /*lint --e{715}*/ 2421 SCIP_Real completion; 2422 2423 SCIP_CALL_ABORT( getSearchCompletion(eventhdlrdata, &completion) ); 2424 2425 /* if the estimation exceeds the current number of nodes by a dramatic factor, restart */ 2426 if( completion < 1.0 / eventhdlrdata->restartfactor ) 2427 { 2428 SCIPstatistic( 2429 SCIPverbMessage(scip, SCIP_VERBLEVEL_FULL, NULL, 2430 "Completion %.5f less than restart threshold %.5f\n", 2431 completion, 1.0 / eventhdlrdata->restartfactor); 2432 ) 2433 2434 return TRUE; 2435 } 2436 2437 return FALSE; 2438 } 2439 2440 /** should a restart be applied based on the value of the selected completion method? */ 2441 static 2442 SCIP_Bool shouldApplyRestartEstimation( 2443 SCIP* scip, /**< SCIP data structure */ 2444 SCIP_EVENTHDLRDATA* eventhdlrdata /**< event handler data */ 2445 ) 2446 { 2447 SCIP_Real estimation; 2448 2449 estimation = SCIPgetTreesizeEstimation(scip); 2450 2451 if( estimation < 0.0 ) 2452 { 2453 SCIPstatistic( 2454 SCIPverbMessage(scip, SCIP_VERBLEVEL_FULL, NULL, 2455 "Estimation %g is still unavailable\n", 2456 estimation); 2457 ) 2458 2459 return TRUE; 2460 } 2461 2462 /* if the estimation exceeds the current number of nodes by a dramatic factor, restart */ 2463 if( estimation > eventhdlrdata->treedata->nnodes * eventhdlrdata->restartfactor ) 2464 { 2465 SCIPstatistic( 2466 SCIPverbMessage(scip, SCIP_VERBLEVEL_FULL, NULL, 2467 "Estimation %g exceeds number of estimation tree nodes %" SCIP_LONGINT_FORMAT " by a factor of %.1f\n", 2468 estimation, eventhdlrdata->treedata->nnodes, estimation / eventhdlrdata->treedata->nnodes); 2469 ) 2470 2471 return TRUE; 2472 } 2473 2474 return FALSE; 2475 } 2476 2477 /** check if a restart should be performed based on the given restart policy */ 2478 static 2479 SCIP_Bool shouldApplyRestart( 2480 SCIP* scip, /**< SCIP data structure */ 2481 SCIP_EVENTHDLRDATA* eventhdlrdata /**< event handler data */ 2482 ) 2483 { 2484 SCIP_Bool applyrestart = FALSE; 2485 2486 switch (getRestartPolicy(eventhdlrdata)) 2487 { 2488 case RESTARTPOLICY_ALWAYS: 2489 applyrestart = TRUE; 2490 break; 2491 case RESTARTPOLICY_NEVER: 2492 applyrestart = FALSE; 2493 break; 2494 case RESTARTPOLICY_COMPLETION: 2495 applyrestart = shouldApplyRestartCompletion(scip, eventhdlrdata); 2496 break; 2497 case RESTARTPOLICY_ESTIMATION: 2498 applyrestart = shouldApplyRestartEstimation(scip, eventhdlrdata); 2499 break; 2500 default: 2501 break; 2502 } 2503 2504 return applyrestart; 2505 } 2506 2507 /** update all time series */ 2508 static 2509 SCIP_RETCODE updateTimeseries( 2510 SCIP* scip, /**< SCIP data structure */ 2511 SCIP_EVENTHDLRDATA* eventhdlrdata, /**< event handler data */ 2512 TREEDATA* treedata, /**< tree data */ 2513 SCIP_Bool isleaf /**< are we at a leaf node? */ 2514 ) 2515 { 2516 TIMESERIES** tss = eventhdlrdata->timeseries; 2517 int t; 2518 2519 /* loop over time series */ 2520 for( t = 0; t < NTIMESERIES; ++t ) 2521 { 2522 assert(tss[t] != NULL); 2523 SCIP_CALL( timeSeriesUpdate(scip, tss[t], treedata, isleaf) ); 2524 2525 #ifdef SCIP_MORE_DEBUG 2526 SCIPdebugMsg(scip, 2527 "Update of time series '%s', current value %.4f (%" SCIP_LONGINT_FORMAT " observations)\n", 2528 timeSeriesGetName(tss[t]), timeSeriesGetValue(tss[t]), tss[t]->nobs); 2529 #endif 2530 } 2531 2532 return SCIP_OKAY; 2533 } 2534 2535 /** print a treesize estimation report into the string buffer */ 2536 static 2537 char* printReport( 2538 SCIP* scip, /**< SCIP data structure */ 2539 SCIP_EVENTHDLRDATA* eventhdlrdata, /**< event handler data */ 2540 char* strbuf, /**< string buffer */ 2541 int reportnum /**< report number, or 0 to omit number */ 2542 ) 2543 { 2544 TREEDATA* treedata = eventhdlrdata->treedata; 2545 char* ptr = strbuf; 2546 SCIP_Real completed; 2547 SCIP_Real wbeestim; 2548 char wbeestimstr[SCIP_MAXSTRLEN]; 2549 int t; 2550 2551 /* print report number */ 2552 if( reportnum > 0 ) 2553 ptr += SCIPsnprintf(ptr, SCIP_MAXSTRLEN, "Report %d\nTime Elapsed: %.2f\n", reportnum, SCIPgetSolvingTime(scip)); 2554 2555 ptr += SCIPsnprintf(ptr, SCIP_MAXSTRLEN, 2556 "Estim. Tree Size :%11" SCIP_LONGINT_FORMAT "\n", 2557 (SCIP_Longint)SCIPgetTreesizeEstimation(scip)); 2558 2559 SCIP_CALL_ABORT( getSearchCompletion(eventhdlrdata, &completed) ); 2560 2561 completed = MIN(1.0, completed); 2562 completed = MAX(0.0, completed); 2563 2564 /* print tree data */ 2565 ptr += SCIPsnprintf(ptr, SCIP_MAXSTRLEN, 2566 "%-19s: %" SCIP_LONGINT_FORMAT " nodes (" 2567 "%" SCIP_LONGINT_FORMAT " visited, " 2568 "%" SCIP_LONGINT_FORMAT " internal, " 2569 "%" SCIP_LONGINT_FORMAT " leaves, " 2570 "%" SCIP_LONGINT_FORMAT " open), " 2571 "weight: %.4Lf completed %.4f\n", 2572 "Estimation Tree", 2573 treedata->nnodes, 2574 treedata->nvisited, 2575 treedata->ninner, 2576 treedata->nleaves, 2577 treedata->nopen, 2578 treedata->weight, 2579 completed 2580 ); 2581 2582 /* print estimations */ 2583 ptr += SCIPsnprintf(ptr, SCIP_MAXSTRLEN, "Estimations : %10s %10s %10s %10s %10s", 2584 "estim", "value", "trend", "resolution", "smooth"); 2585 ptr += SCIPsnprintf(ptr, SCIP_MAXSTRLEN, "\n"); 2586 2587 wbeestim = treeDataGetWbe(eventhdlrdata->treedata); 2588 ptr += SCIPsnprintf(ptr, SCIP_MAXSTRLEN, " wbe : %10s %10s %10s %10s %10s\n", 2589 real2String(wbeestim, wbeestimstr, 0), "-", "-", "-", "-"); 2590 2591 ptr += SCIPsnprintf(ptr, SCIP_MAXSTRLEN, " tree-profile : %10.0f %10s %10s %10s %10s\n", 2592 predictTotalSizeTreeProfile(scip, eventhdlrdata->treeprofile, eventhdlrdata->treeprofile_minnodesperdepth), 2593 "-", "-", "-", "-"); 2594 2595 /* print time series forecasts */ 2596 for( t = 0; t < NTIMESERIES; ++t ) 2597 { 2598 SCIP_Real trend; 2599 SCIP_Real smoothestim; 2600 TIMESERIES* ts = eventhdlrdata->timeseries[t]; 2601 char trendstr[SCIP_MAXSTRLEN]; 2602 char smoothestimstr[SCIP_MAXSTRLEN]; 2603 2604 trend = doubleExpSmoothGetTrend(&ts->des); 2605 smoothestim = timeSeriesGetSmoothEstimation(ts); 2606 2607 ptr += SCIPsnprintf(ptr, SCIP_MAXSTRLEN, " %-17s: %10.0f %10.5f %10s %10d %10s\n", 2608 timeSeriesGetName(ts), 2609 timeSeriesEstimate(ts, eventhdlrdata->treedata), 2610 timeSeriesGetValue(ts), 2611 real2String(trend, trendstr, 5), 2612 timeSeriesGetResolution(ts), 2613 real2String(smoothestim, smoothestimstr, 0)); 2614 } 2615 2616 if( reportnum > 0 ) 2617 (void) SCIPsnprintf(ptr, SCIP_MAXSTRLEN, "End of Report %d\n", reportnum); 2618 2619 return strbuf; 2620 } 2621 2622 2623 /** copy method for event handler plugins (called when SCIP copies plugins) */ 2624 static 2625 SCIP_DECL_EVENTCOPY(eventCopyEstim) 2626 { /*lint --e{715}*/ 2627 assert(scip != NULL); 2628 2629 SCIP_CALL( SCIPincludeEventHdlrEstim(scip) ); 2630 2631 return SCIP_OKAY; 2632 } 2633 2634 /** destructor of event handler to free user data (called when SCIP is exiting) */ 2635 static 2636 SCIP_DECL_EVENTFREE(eventFreeEstim) 2637 { /*lint --e{715}*/ 2638 SCIP_EVENTHDLRDATA* eventhdlrdata; 2639 2640 eventhdlrdata = SCIPeventhdlrGetData(eventhdlr); 2641 assert(eventhdlrdata != NULL); 2642 2643 freeTreeData(scip, &eventhdlrdata->treedata); 2644 2645 freeTimeSeries(scip, eventhdlrdata); 2646 2647 SCIPfreeMemory(scip, &eventhdlrdata); 2648 2649 return SCIP_OKAY; 2650 } 2651 2652 /** initialization method of event handler (called after problem was transformed) */ 2653 static 2654 SCIP_DECL_EVENTINIT(eventInitEstim) 2655 { /*lint --e{715}*/ 2656 SCIP_EVENTHDLRDATA* eventhdlrdata; 2657 2658 eventhdlrdata = SCIPeventhdlrGetData(eventhdlr); 2659 assert(eventhdlrdata != NULL); 2660 2661 /* test if user specified a regression forest */ 2662 if( 0 != strncmp(eventhdlrdata->regforestfilename, DEFAULT_REGFORESTFILENAME, strlen(DEFAULT_REGFORESTFILENAME)) ) 2663 { 2664 SCIP_CALL( SCIPregForestFromFile(&eventhdlrdata->regforest, eventhdlrdata->regforestfilename) ); 2665 } 2666 2667 eventhdlrdata->lastrestartrun = 0; 2668 eventhdlrdata->nrestartsperformed = 0; 2669 2670 return SCIP_OKAY; 2671 } 2672 2673 /** deinitialization method of event handler (called before transformed problem is freed) */ 2674 static 2675 SCIP_DECL_EVENTEXIT(eventExitEstim) 2676 { /*lint --e{715}*/ 2677 SCIP_EVENTHDLRDATA* eventhdlrdata; 2678 2679 eventhdlrdata = SCIPeventhdlrGetData(eventhdlr); 2680 assert(eventhdlrdata != NULL); 2681 2682 SCIPregForestFree(&eventhdlrdata->regforest); 2683 2684 return SCIP_OKAY; 2685 } 2686 2687 /** solving process initialization method of event handler (called when branch and bound process is about to begin) */ 2688 static 2689 SCIP_DECL_EVENTINITSOL(eventInitsolEstim) 2690 { /*lint --e{715}*/ 2691 SCIP_EVENTHDLRDATA* eventhdlrdata; 2692 2693 eventhdlrdata = SCIPeventhdlrGetData(eventhdlr); 2694 assert(eventhdlrdata != NULL); 2695 2696 eventhdlrdata->restarthitcounter = 0; 2697 eventhdlrdata->weightlastreport = 0.0; 2698 eventhdlrdata->nreports = 0; 2699 2700 /* reset tree data */ 2701 SCIP_CALL( resetTreeData(scip, eventhdlrdata->treedata) ); 2702 2703 resetTimeSeries(eventhdlrdata); 2704 2705 SCIP_CALL( SCIPcatchEvent(scip, EVENTTYPE_ESTIM, eventhdlr, NULL, NULL) ); 2706 2707 if( eventhdlrdata->treeprofile_enabled ) 2708 { 2709 SCIP_CALL( createTreeProfile(scip, &eventhdlrdata->treeprofile) ); 2710 } 2711 2712 eventhdlrdata->treeisbinary = TRUE; 2713 2714 return SCIP_OKAY; 2715 } 2716 2717 /** solving process deinitialization method of event handler (called before branch and bound process data is freed) */ 2718 static 2719 SCIP_DECL_EVENTEXITSOL(eventExitsolEstim) 2720 { /*lint --e{715}*/ 2721 SCIP_EVENTHDLRDATA* eventhdlrdata; 2722 2723 eventhdlrdata = SCIPeventhdlrGetData(eventhdlr); 2724 assert(eventhdlrdata != NULL); 2725 2726 if( eventhdlrdata->treeprofile != NULL ) 2727 freeTreeProfile(scip, &eventhdlrdata->treeprofile); 2728 2729 SCIP_CALL( SCIPdropEvent(scip, EVENTTYPE_ESTIM, eventhdlr, NULL, -1) ); 2730 2731 return SCIP_OKAY; 2732 } 2733 2734 /** execution method of event handler */ 2735 static 2736 SCIP_DECL_EVENTEXEC(eventExecEstim) 2737 { /*lint --e{715}*/ 2738 SCIP_EVENTHDLRDATA* eventhdlrdata; 2739 SCIP_Bool isleaf; 2740 SCIP_EVENTTYPE eventtype; 2741 TREEDATA* treedata; 2742 char strbuf[SCIP_MAXSTRLEN]; 2743 2744 assert(scip != NULL); 2745 assert(eventhdlr != NULL); 2746 2747 eventhdlrdata = SCIPeventhdlrGetData(eventhdlr); 2748 assert(eventhdlrdata != NULL); 2749 eventtype = SCIPeventGetType(event); 2750 treedata = eventhdlrdata->treedata; 2751 2752 /* actual leaf nodes for our tree data are children/siblings/leaves or the focus node itself (deadend) 2753 * if it has not been branched on 2754 */ 2755 isleaf = (eventtype == SCIP_EVENTTYPE_NODEDELETE) && 2756 (SCIPnodeGetType(SCIPeventGetNode(event)) == SCIP_NODETYPE_CHILD || 2757 SCIPnodeGetType(SCIPeventGetNode(event)) == SCIP_NODETYPE_SIBLING || 2758 SCIPnodeGetType(SCIPeventGetNode(event)) == SCIP_NODETYPE_LEAF || 2759 (SCIPnodeGetType(SCIPeventGetNode(event)) == SCIP_NODETYPE_DEADEND && !SCIPwasNodeLastBranchParent(scip, SCIPeventGetNode(event)))); 2760 2761 if( eventtype == SCIP_EVENTTYPE_NODEBRANCHED || isleaf ) 2762 { 2763 SCIP_NODE* eventnode; 2764 int nchildren = 0; 2765 2766 if( eventtype == SCIP_EVENTTYPE_NODEBRANCHED ) 2767 { 2768 nchildren = SCIPgetNChildren(scip); 2769 2770 /* update whether the tree is still binary */ 2771 if( nchildren != 2 ) 2772 eventhdlrdata->treeisbinary = FALSE; 2773 } 2774 2775 eventnode = SCIPeventGetNode(event); 2776 SCIP_CALL( updateTreeData(scip, treedata, eventnode, nchildren) ); 2777 SCIP_CALL( updateTreeProfile(scip, eventhdlrdata->treeprofile, eventnode) ); 2778 2779 #ifdef SCIP_DEBUG 2780 SCIPdebugMsg(scip, "%s\n", treeDataPrint(treedata, strbuf)); 2781 #endif 2782 2783 SCIP_CALL( updateTimeseries(scip, eventhdlrdata, treedata, nchildren == 0) ); 2784 2785 /* should a new report be printed? */ 2786 if( eventhdlrdata->reportfreq >= 0 && SCIPgetStatus(scip) == SCIP_STATUS_UNKNOWN && 2787 (eventhdlrdata->reportfreq == 0 2788 || treedata->weight >= eventhdlrdata->weightlastreport + 1.0 / (SCIP_Real)eventhdlrdata->reportfreq) ) 2789 { 2790 SCIPverbMessage(scip, SCIP_VERBLEVEL_HIGH, NULL, "%s\n", printReport(scip, eventhdlrdata, strbuf, ++eventhdlrdata->nreports)); 2791 2792 if( eventhdlrdata->reportfreq > 0 ) 2793 eventhdlrdata->weightlastreport = 1 / (SCIP_Real)eventhdlrdata->reportfreq * SCIPfloor(scip, ((SCIP_Real)treedata->weight * eventhdlrdata->reportfreq)); 2794 else 2795 eventhdlrdata->weightlastreport = (SCIP_Real)treedata->weight; 2796 } 2797 } 2798 2799 /* if nodes have been pruned, things are progressing, don't restart right now */ 2800 if( isleaf ) 2801 return SCIP_OKAY; 2802 2803 /* check if all conditions are met such that the event handler should run */ 2804 if( ! isRestartApplicable(scip, eventhdlrdata) ) 2805 return SCIP_OKAY; 2806 2807 /* test if a restart should be applied */ 2808 if( shouldApplyRestart(scip, eventhdlrdata) ) 2809 { 2810 eventhdlrdata->restarthitcounter++; 2811 2812 if( eventhdlrdata->restarthitcounter >= eventhdlrdata->hitcounterlim ) 2813 { 2814 /* safe that we triggered a restart at this run */ 2815 if( SCIPgetNRuns(scip) > eventhdlrdata->lastrestartrun ) 2816 { 2817 eventhdlrdata->nrestartsperformed++; 2818 2819 SCIPverbMessage(scip, SCIP_VERBLEVEL_HIGH, NULL, 2820 "Restart triggered after %d consecutive estimations that the remaining tree will be large\n", 2821 eventhdlrdata->restarthitcounter); 2822 } 2823 2824 eventhdlrdata->lastrestartrun = SCIPgetNRuns(scip); 2825 2826 SCIP_CALL( SCIPrestartSolve(scip) ); 2827 } 2828 } 2829 else 2830 { 2831 eventhdlrdata->restarthitcounter = 0; 2832 } 2833 2834 return SCIP_OKAY; 2835 } 2836 2837 /** output method of statistics table to output file stream 'file' */ 2838 static 2839 SCIP_DECL_TABLEOUTPUT(tableOutputEstim) 2840 { /*lint --e{715}*/ 2841 SCIP_EVENTHDLR* eventhdlr; 2842 SCIP_EVENTHDLRDATA* eventhdlrdata; 2843 char strbuf[SCIP_MAXSTRLEN]; 2844 2845 eventhdlr = SCIPfindEventhdlr(scip, EVENTHDLR_NAME); 2846 assert(eventhdlr != NULL); 2847 2848 eventhdlrdata = SCIPeventhdlrGetData(eventhdlr); 2849 assert(eventhdlrdata != NULL); 2850 2851 if( eventhdlrdata->showstats ) 2852 SCIPinfoMessage(scip, file, "%s", printReport(scip, eventhdlrdata, strbuf, 0)); 2853 2854 return SCIP_OKAY; 2855 } 2856 2857 /** output method of search tree completion display column to output file stream 'file' */ 2858 static 2859 SCIP_DECL_DISPOUTPUT(dispOutputCompleted) 2860 { /*lint --e{715}*/ 2861 SCIP_EVENTHDLR* eventhdlr; 2862 SCIP_EVENTHDLRDATA* eventhdlrdata; 2863 TREEDATA* treedata; 2864 SCIP_Real completed; 2865 2866 assert(disp != NULL); 2867 assert(strcmp(SCIPdispGetName(disp), DISP_NAME) == 0); 2868 assert(scip != NULL); 2869 2870 eventhdlr = SCIPfindEventhdlr(scip, EVENTHDLR_NAME); 2871 assert(eventhdlr != NULL); 2872 eventhdlrdata = SCIPeventhdlrGetData(eventhdlr); 2873 assert(eventhdlrdata != NULL); 2874 treedata = eventhdlrdata->treedata; 2875 2876 SCIP_CALL( getSearchCompletion(eventhdlrdata, &completed) ); 2877 2878 completed = MIN(completed, 1.0); 2879 2880 if( treedata->weight >= 0.005 && completed > 0 ) 2881 SCIPinfoMessage(scip, file, "%7.2f%%", 100.0 * completed); 2882 else 2883 SCIPinfoMessage(scip, file, " unknown"); 2884 2885 return SCIP_OKAY; 2886 } 2887 2888 /** creates event handler for tree size estimation */ 2889 SCIP_RETCODE SCIPincludeEventHdlrEstim( 2890 SCIP* scip /**< SCIP data structure */ 2891 ) 2892 { 2893 SCIP_RETCODE retcode; 2894 SCIP_EVENTHDLRDATA* eventhdlrdata = NULL; 2895 SCIP_EVENTHDLR* eventhdlr = NULL; 2896 2897 /* create estim event handler data */ 2898 SCIP_CALL( SCIPallocMemory(scip, &eventhdlrdata) ); 2899 BMSclearMemory(eventhdlrdata); 2900 2901 SCIP_CALL_TERMINATE( retcode, createTreeData(scip, &eventhdlrdata->treedata), TERMINATE ); 2902 2903 SCIP_CALL( SCIPincludeEventhdlrBasic(scip, &eventhdlr, EVENTHDLR_NAME, EVENTHDLR_DESC, 2904 eventExecEstim, eventhdlrdata) ); 2905 assert(eventhdlr != NULL); 2906 2907 /* set non fundamental callbacks via setter functions */ 2908 SCIP_CALL( SCIPsetEventhdlrCopy(scip, eventhdlr, eventCopyEstim) ); 2909 SCIP_CALL( SCIPsetEventhdlrFree(scip, eventhdlr, eventFreeEstim) ); 2910 SCIP_CALL( SCIPsetEventhdlrInit(scip, eventhdlr, eventInitEstim) ); 2911 SCIP_CALL( SCIPsetEventhdlrExit(scip, eventhdlr, eventExitEstim) ); 2912 SCIP_CALL( SCIPsetEventhdlrInitsol(scip, eventhdlr, eventInitsolEstim) ); 2913 SCIP_CALL( SCIPsetEventhdlrExitsol(scip, eventhdlr, eventExitsolEstim) ); 2914 2915 /* add estimation event handler parameters */ 2916 SCIP_CALL( SCIPaddCharParam(scip, "estimation/restarts/restartpolicy", "restart policy: (a)lways, (c)ompletion, (e)stimation, (n)ever", 2917 &eventhdlrdata->restartpolicyparam, FALSE, DEFAULT_RESTARTPOLICY, "acen", NULL, NULL) ); 2918 2919 SCIP_CALL( SCIPaddCharParam(scip, "estimation/method", 2920 "tree size estimation method: (c)ompletion, (e)nsemble, " 2921 "time series forecasts on either (g)ap, (l)eaf frequency, (o)open nodes, tree (w)eight, (s)sg, " 2922 "or (t)ree profile or w(b)e", 2923 &eventhdlrdata->estimmethod, FALSE, DEFAULT_ESTIMMETHOD, ESTIMMETHODS, NULL, NULL) ); 2924 2925 SCIP_CALL( SCIPaddIntParam(scip, "estimation/restarts/restartlimit", "restart limit", 2926 &eventhdlrdata->restartlimit, FALSE, DEFAULT_RESTARTLIMIT, -1, INT_MAX, NULL, NULL) ); 2927 2928 SCIP_CALL( SCIPaddLongintParam(scip, "estimation/restarts/minnodes", "minimum number of nodes before restart", 2929 &eventhdlrdata->minnodes, FALSE, DEFAULT_MINNODES, -1L, SCIP_LONGINT_MAX, NULL, NULL) ); 2930 2931 SCIP_CALL( SCIPaddBoolParam(scip, "estimation/restarts/countonlyleaves", "should only leaves count for the minnodes parameter?", 2932 &eventhdlrdata->countonlyleaves, DEFAULT_COUNTONLYLEAVES, FALSE, NULL, NULL) ); 2933 2934 SCIP_CALL( SCIPaddRealParam(scip, "estimation/restarts/restartfactor", 2935 "factor by which the estimated number of nodes should exceed the current number of nodes", 2936 &eventhdlrdata->restartfactor, FALSE, DEFAULT_RESTARTFACTOR, 1.0, SCIP_REAL_MAX, NULL, NULL) ); 2937 2938 SCIP_CALL( SCIPaddBoolParam(scip, "estimation/restarts/restartnonlinear", 2939 "whether to apply a restart when nonlinear constraints are present", 2940 &eventhdlrdata->restartnonlinear, FALSE, DEFAULT_RESTARTNONLINEAR, NULL, NULL) ); 2941 2942 SCIP_CALL( SCIPaddBoolParam(scip, "estimation/restarts/restartactpricers", 2943 "whether to apply a restart when active pricers are used", 2944 &eventhdlrdata->restartactpricers, FALSE, DEFAULT_RESTARTACTPRICERS, NULL, NULL) ); 2945 2946 SCIP_CALL( SCIPaddRealParam(scip, "estimation/coefmonoweight", 2947 "coefficient of tree weight in monotone approximation of search completion", 2948 &eventhdlrdata->coefmonoweight, FALSE, DEFAULT_COEFMONOWEIGHT, 0.0, 1.0, NULL, NULL) ); 2949 2950 SCIP_CALL( SCIPaddRealParam(scip, "estimation/coefmonossg", 2951 "coefficient of 1 - SSG in monotone approximation of search completion", 2952 &eventhdlrdata->coefmonossg, FALSE, DEFAULT_COEFMONOSSG, 0.0, 1.0, NULL, NULL) ); 2953 2954 SCIP_CALL( SCIPaddIntParam(scip, "estimation/restarts/hitcounterlim", "limit on the number of successive samples to really trigger a restart", 2955 &eventhdlrdata->hitcounterlim, FALSE, DEFAULT_HITCOUNTERLIM, 1, INT_MAX, NULL, NULL) ); 2956 2957 SCIP_CALL( SCIPaddIntParam(scip, "estimation/reportfreq", 2958 "report frequency on estimation: -1: never, 0:always, k >= 1: k times evenly during search", 2959 &eventhdlrdata->reportfreq, TRUE, DEFAULT_REPORTFREQ, -1, INT_MAX / 2, NULL, NULL) ); 2960 2961 SCIP_CALL( SCIPaddStringParam(scip, "estimation/regforestfilename", "user regression forest in RFCSV format", 2962 &eventhdlrdata->regforestfilename, FALSE, DEFAULT_REGFORESTFILENAME, NULL, NULL) ); 2963 2964 SCIP_CALL( SCIPaddCharParam(scip, "estimation/completiontype", 2965 "approximation of search tree completion: (a)uto, (g)ap, tree (w)eight, (m)onotone regression, (r)egression forest, (s)sg", 2966 &eventhdlrdata->completiontypeparam, FALSE, DEFAULT_COMPLETIONTYPE, "agmrsw", NULL, NULL) ); 2967 2968 SCIP_CALL( SCIPaddBoolParam(scip, "estimation/treeprofile/enabled", 2969 "should the event handler collect data?", 2970 &eventhdlrdata->treeprofile_enabled, FALSE, DEFAULT_TREEPROFILE_ENABLED, NULL, NULL) ); 2971 2972 SCIP_CALL( SCIPaddRealParam(scip, "estimation/treeprofile/minnodesperdepth", 2973 "minimum average number of nodes at each depth before producing estimations", 2974 &eventhdlrdata->treeprofile_minnodesperdepth, FALSE, DEFAULT_TREEPROFILE_MINNODESPERDEPTH, 1.0, SCIP_REAL_MAX, NULL, NULL) ); 2975 2976 SCIP_CALL( SCIPaddBoolParam(scip, "estimation/useleafts", 2977 "use leaf nodes as basic observations for time series, or all nodes?", 2978 &eventhdlrdata->useleafts, TRUE, DEFAULT_USELEAFTS, NULL, NULL) ); 2979 2980 SCIP_CALL( SCIPaddBoolParam(scip, "estimation/showstats", 2981 "should statistics be shown at the end?", 2982 &eventhdlrdata->showstats, TRUE, DEFAULT_SHOWSTATS, NULL, NULL) ); 2983 2984 /* SSG parameters */ 2985 SCIP_CALL( SCIPaddIntParam(scip, "estimation/ssg/nmaxsubtrees", 2986 "the maximum number of individual SSG subtrees; -1: no limit", 2987 &eventhdlrdata->treedata->ssg->nmaxsubtrees, FALSE, DEFAULT_SSG_NMAXSUBTREES, -1, INT_MAX / 2, NULL, NULL) ); 2988 2989 SCIP_CALL( SCIPaddLongintParam(scip, "estimation/ssg/nminnodeslastsplit", 2990 "minimum number of nodes to process between two consecutive SSG splits", 2991 &eventhdlrdata->treedata->ssg->nminnodeslastsplit, FALSE, DEFAULT_SSG_NMINNODESLASTSPLIT, 0L, SCIP_LONGINT_MAX, NULL, NULL) ); 2992 2993 /* include statistics table */ 2994 SCIP_CALL( SCIPincludeTable(scip, TABLE_NAME, TABLE_DESC, TRUE, 2995 NULL, NULL, NULL, NULL, NULL, NULL, tableOutputEstim, 2996 NULL, TABLE_POSITION, TABLE_EARLIEST_STAGE) ); 2997 2998 /* include time series into event handler */ 2999 SCIP_CALL( includeTimeseries(scip, eventhdlrdata) ); 3000 3001 /* include display column */ 3002 SCIP_CALL( SCIPincludeDisp(scip, DISP_NAME, DISP_DESC, DISP_HEADER, SCIP_DISPSTATUS_AUTO, 3003 NULL, NULL, NULL, NULL, NULL, NULL, dispOutputCompleted, 3004 NULL, DISP_WIDTH, DISP_PRIORITY, DISP_POSITION, DISP_STRIPLINE) ); 3005 3006 TERMINATE: 3007 if( retcode != SCIP_OKAY ) 3008 { 3009 freeTreeData(scip, &eventhdlrdata->treedata); 3010 SCIPfreeMemory(scip, &eventhdlrdata); 3011 } 3012 3013 return retcode; 3014 } 3015 3016 /** return an estimation of the final tree size */ 3017 SCIP_Real SCIPgetTreesizeEstimation( 3018 SCIP* scip /**< SCIP data structure */ 3019 ) 3020 { 3021 SCIP_EVENTHDLR* eventhdlr; 3022 SCIP_EVENTHDLRDATA* eventhdlrdata; 3023 TSPOS tspos = TSPOS_NONE; 3024 SCIP_Real estim; 3025 3026 assert(scip != NULL); 3027 3028 eventhdlr = SCIPfindEventhdlr(scip, EVENTHDLR_NAME); 3029 if( eventhdlr == NULL ) 3030 { 3031 SCIPwarningMessage(scip, "SCIPgetTreesizeEstimation() called, but event handler " EVENTHDLR_NAME " is missing.\n"); 3032 return -1.0; 3033 } 3034 3035 eventhdlrdata = SCIPeventhdlrGetData(eventhdlr); 3036 assert(eventhdlrdata != NULL); 3037 3038 switch (eventhdlrdata->estimmethod) 3039 { 3040 case ESTIMMETHOD_COMPL: 3041 SCIP_CALL_ABORT( getEstimCompletion(scip, eventhdlrdata, &estim) ); 3042 return estim; 3043 3044 case ESTIMMETHOD_ENSMBL: 3045 return getEnsembleEstimation(eventhdlrdata); 3046 3047 /* for the requested time series methods, we specify the array position */ 3048 case ESTIMMETHOD_GAP: 3049 tspos = TSPOS_GAP; 3050 break; 3051 3052 case ESTIMMETHOD_LFREQ: 3053 tspos = TSPOS_LFREQ; 3054 break; 3055 3056 case ESTIMMETHOD_OPEN: 3057 tspos = TSPOS_OPEN; 3058 break; 3059 3060 case ESTIMMETHOD_TREEWEIGHT: 3061 tspos = TSPOS_TREEWEIGHT; 3062 break; 3063 3064 case ESTIMMETHOD_SSG: 3065 tspos = TSPOS_SSG; 3066 break; 3067 3068 /* tree profile estimation */ 3069 case ESTIMMETHOD_TPROF: 3070 return predictTotalSizeTreeProfile(scip, eventhdlrdata->treeprofile, eventhdlrdata->treeprofile_minnodesperdepth); 3071 3072 /* Weighted backtrack estimation */ 3073 case ESTIMMETHOD_WBE: 3074 return treeDataGetWbe(eventhdlrdata->treedata); 3075 3076 default: 3077 SCIPerrorMessage("Unknown estimation '%c' method specified, should be one of [%s]\n", 3078 eventhdlrdata->estimmethod, ESTIMMETHODS); 3079 SCIPABORT(); 3080 break; 3081 } 3082 3083 assert(tspos != TSPOS_NONE); 3084 return (tspos == TSPOS_NONE ? -1.0 : timeSeriesEstimate(eventhdlrdata->timeseries[tspos], eventhdlrdata->treedata)); 3085 } 3086