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 presol_milp.cpp 26 * @brief MILP presolver 27 * @author Leona Gottwald 28 * 29 * Calls the presolve library and communicates (multi-)aggregations, fixings, and bound 30 * changes to SCIP by utilizing the postsolve information. Constraint changes can currently 31 * only be communicated by deleting all constraints and adding new ones. 32 * 33 * @todo add infrastructure to SCIP for handling parallel columns 34 * @todo better communication of constraint changes by adding more information to the postsolve structure 35 * @todo allow to pass additional external locks to the presolve library that are considered when doing reductions 36 * 37 */ 38 39 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/ 40 41 #include "scip/presol_milp.h" 42 43 #ifndef SCIP_WITH_PAPILO 44 45 /** creates the MILP presolver and includes it in SCIP */ 46 SCIP_RETCODE SCIPincludePresolMILP( 47 SCIP* scip /**< SCIP data structure */ 48 ) 49 { 50 assert(scip != NULL); 51 return SCIP_OKAY; 52 } 53 54 #else 55 56 /* disable some warnings that come up in header files of PAPILOs dependencies */ 57 #ifdef __GNUC__ 58 #pragma GCC diagnostic ignored "-Wshadow" 59 #pragma GCC diagnostic ignored "-Wctor-dtor-privacy" 60 #pragma GCC diagnostic ignored "-Wredundant-decls" 61 62 /* disable false warning, !3076, https://gcc.gnu.org/bugzilla/show_bug.cgi?id=106199 */ 63 #if __GNUC__ == 12 && __GNUC__MINOR__ <= 2 64 #pragma GCC diagnostic ignored "-Wstringop-overflow" 65 #endif 66 #endif 67 68 #include <assert.h> 69 #include "scip/cons_linear.h" 70 #include "scip/pub_matrix.h" 71 #include "scip/pub_presol.h" 72 #include "scip/pub_var.h" 73 #include "scip/pub_cons.h" 74 #include "scip/pub_message.h" 75 #include "scip/scip_general.h" 76 #include "scip/scip_presol.h" 77 #include "scip/scip_var.h" 78 #include "scip/scip_mem.h" 79 #include "scip/scip_prob.h" 80 #include "scip/scip_param.h" 81 #include "scip/scip_cons.h" 82 #include "scip/scip_numerics.h" 83 #include "scip/scip_timing.h" 84 #include "scip/scip_message.h" 85 #include "scip/scip_randnumgen.h" 86 #include "papilo/core/Presolve.hpp" 87 #include "papilo/core/ProblemBuilder.hpp" 88 #include "papilo/Config.hpp" 89 90 #define PRESOL_NAME "milp" 91 #define PRESOL_DESC "MILP specific presolving methods" 92 #define PRESOL_PRIORITY 9999999 /**< priority of the presolver (>= 0: before, < 0: after constraint handlers); combined with propagators */ 93 #define PRESOL_MAXROUNDS -1 /**< maximal number of presolving rounds the presolver participates in (-1: no limit) */ 94 #define PRESOL_TIMING SCIP_PRESOLTIMING_MEDIUM /* timing of the presolver (fast, medium, or exhaustive) */ 95 96 /* default parameter values */ 97 #define DEFAULT_THREADS 1 /**< maximum number of threads presolving may use (0: automatic) */ 98 #define DEFAULT_MAXFILLINPERSUBST 3 /**< maximal possible fillin for substitutions to be considered */ 99 #define DEFAULT_MAXSHIFTPERROW 10 /**< maximal amount of nonzeros allowed to be shifted to make space for substitutions */ 100 #define DEFAULT_DETECTLINDEP 0 /**< should linear dependent equations and free columns be removed? (0: never, 1: for LPs, 2: always) */ 101 #define DEFAULT_MAXBADGESIZE_SEQ 15000 /**< the max badge size in Probing if PaPILO is executed in sequential mode*/ 102 #define DEFAULT_MAXBADGESIZE_PAR -1 /**< the max badge size in Probing if PaPILO is executed in parallel mode*/ 103 #define DEFAULT_RANDOMSEED 0 /**< the random seed used for randomization of tie breaking */ 104 #define DEFAULT_MODIFYCONSFAC 0.8 /**< modify SCIP constraints when the number of nonzeros or rows is at most this 105 * factor times the number of nonzeros or rows before presolving */ 106 #define DEFAULT_MARKOWITZTOLERANCE 0.01 /**< the markowitz tolerance used for substitutions */ 107 #define DEFAULT_VERBOSITY 0 108 #define DEFAULT_HUGEBOUND 1e8 /**< absolute bound value that is considered too huge for activitity based calculations */ 109 #define DEFAULT_ENABLEPARALLELROWS TRUE /**< should the parallel rows presolver be enabled within the presolve library? */ 110 #define DEFAULT_ENABLEDOMCOL TRUE /**< should the dominated column presolver be enabled within the presolve library? */ 111 #define DEFAULT_ENABLEDUALINFER TRUE /**< should the dualinfer presolver be enabled within the presolve library? */ 112 #define DEFAULT_ENABLEMULTIAGGR TRUE /**< should the multi-aggregation presolver be enabled within the presolve library? */ 113 #define DEFAULT_ENABLEPROBING TRUE /**< should the probing presolver be enabled within the presolve library? */ 114 #define DEFAULT_ENABLESPARSIFY FALSE /**< should the sparsify presolver be enabled within the presolve library? */ 115 #define DEFAULT_FILENAME_PROBLEM "-" /**< default filename to store the instance before presolving */ 116 117 /* 118 * Data structures 119 */ 120 121 /** presolver data */ 122 struct SCIP_PresolData 123 { 124 int lastncols; /**< the number of columns from the last call */ 125 int lastnrows; /**< the number of rows from the last call */ 126 int threads; /**< maximum number of threads presolving may use (0: automatic) */ 127 int maxfillinpersubstitution; /**< maximal possible fillin for substitutions to be considered */ 128 int maxbadgesizeseq; /**< the max badge size in Probing if PaPILO is called in sequential mode*/ 129 int maxbadgesizepar; /**< the max badge size in Probing if PaPILO is called in parallel mode */ 130 int maxshiftperrow; /**< maximal amount of nonzeros allowed to be shifted to make space for substitutions */ 131 int detectlineardependency; /**< should linear dependent equations and free columns be removed? (0: never, 1: for LPs, 2: always) */ 132 int randomseed; /**< the random seed used for randomization of tie breaking */ 133 int verbosity; 134 135 SCIP_Bool enablesparsify; /**< should the sparsify presolver be enabled within the presolve library? */ 136 SCIP_Bool enabledomcol; /**< should the dominated column presolver be enabled within the presolve library? */ 137 SCIP_Bool enableprobing; /**< should the probing presolver be enabled within the presolve library? */ 138 SCIP_Bool enabledualinfer; /**< should the dualinfer presolver be enabled within the presolve library? */ 139 SCIP_Bool enablemultiaggr; /**< should the multi-aggregation presolver be enabled within the presolve library? */ 140 SCIP_Bool enableparallelrows; /**< should the parallel rows presolver be enabled within the presolve library? */ 141 SCIP_Real modifyconsfac; /**< modify SCIP constraints when the number of nonzeros or rows is at most this 142 * factor times the number of nonzeros or rows before presolving */ 143 SCIP_Real markowitztolerance; /**< the markowitz tolerance used for substitutions */ 144 SCIP_Real hugebound; /**< absolute bound value that is considered too huge for activitity based calculations */ 145 146 char* filename = NULL; /**< filename to store the instance before presolving */ 147 }; 148 149 using namespace papilo; 150 151 /* 152 * Local methods 153 */ 154 155 /** builds the presolvelib problem datastructure from the matrix */ 156 static 157 Problem<SCIP_Real> buildProblem( 158 SCIP* scip, /**< SCIP data structure */ 159 SCIP_MATRIX* matrix /**< initialized SCIP_MATRIX data structure */ 160 ) 161 { 162 ProblemBuilder<SCIP_Real> builder; 163 164 /* build problem from matrix */ 165 int nnz = SCIPmatrixGetNNonzs(matrix); 166 int ncols = SCIPmatrixGetNColumns(matrix); 167 int nrows = SCIPmatrixGetNRows(matrix); 168 builder.reserve(nnz, nrows, ncols); 169 170 /* set up columns */ 171 builder.setNumCols(ncols); 172 for(int i = 0; i != ncols; ++i) 173 { 174 SCIP_VAR* var = SCIPmatrixGetVar(matrix, i); 175 SCIP_Real lb = SCIPvarGetLbGlobal(var); 176 SCIP_Real ub = SCIPvarGetUbGlobal(var); 177 builder.setColLb(i, lb); 178 builder.setColUb(i, ub); 179 builder.setColLbInf(i, SCIPisInfinity(scip, -lb)); 180 builder.setColUbInf(i, SCIPisInfinity(scip, ub)); 181 #if PAPILO_VERSION_MAJOR > 2 || (PAPILO_VERSION_MAJOR == 2 && PAPILO_VERSION_MINOR >= 1) 182 if ( SCIPvarGetType(var) == SCIP_VARTYPE_IMPLINT ) 183 builder.setColImplInt(i, TRUE); 184 else 185 builder.setColIntegral(i, SCIPvarIsIntegral(var)); 186 #else 187 builder.setColIntegral(i, SCIPvarIsIntegral(var)); 188 #endif 189 builder.setObj(i, SCIPvarGetObj(var)); 190 } 191 192 /* set up rows */ 193 builder.setNumRows(nrows); 194 for(int i = 0; i != nrows; ++i) 195 { 196 int* rowcols = SCIPmatrixGetRowIdxPtr(matrix, i); 197 SCIP_Real* rowvals = SCIPmatrixGetRowValPtr(matrix, i); 198 int rowlen = SCIPmatrixGetRowNNonzs(matrix, i); 199 builder.addRowEntries(i, rowlen, rowcols, rowvals); 200 201 SCIP_Real lhs = SCIPmatrixGetRowLhs(matrix, i); 202 SCIP_Real rhs = SCIPmatrixGetRowRhs(matrix, i); 203 builder.setRowLhs(i, lhs); 204 builder.setRowRhs(i, rhs); 205 builder.setRowLhsInf(i, SCIPisInfinity(scip, -lhs)); 206 builder.setRowRhsInf(i, SCIPisInfinity(scip, rhs)); 207 } 208 209 /* init objective offset - the value itself is irrelevant */ 210 builder.setObjOffset(0); 211 212 return builder.build(); 213 } 214 215 /* 216 * Callback methods of presolver 217 */ 218 219 /** copy method for constraint handler plugins (called when SCIP copies plugins) */ 220 static 221 SCIP_DECL_PRESOLCOPY(presolCopyMILP) 222 { /*lint --e{715}*/ 223 SCIP_CALL( SCIPincludePresolMILP(scip) ); 224 225 return SCIP_OKAY; 226 } 227 228 /** destructor of presolver to free user data (called when SCIP is exiting) */ 229 static 230 SCIP_DECL_PRESOLFREE(presolFreeMILP) 231 { /*lint --e{715}*/ 232 SCIP_PRESOLDATA* data = SCIPpresolGetData(presol); 233 assert(data != NULL); 234 235 SCIPpresolSetData(presol, NULL); 236 SCIPfreeBlockMemory(scip, &data); 237 return SCIP_OKAY; 238 } 239 240 /** initialization method of presolver (called after problem was transformed) */ 241 static 242 SCIP_DECL_PRESOLINIT(presolInitMILP) 243 { /*lint --e{715}*/ 244 SCIP_PRESOLDATA* data = SCIPpresolGetData(presol); 245 assert(data != NULL); 246 247 data->lastncols = -1; 248 data->lastnrows = -1; 249 250 return SCIP_OKAY; 251 } 252 253 /** execution method of presolver */ 254 static 255 SCIP_DECL_PRESOLEXEC(presolExecMILP) 256 { /*lint --e{715}*/ 257 SCIP_MATRIX* matrix; 258 SCIP_PRESOLDATA* data; 259 SCIP_Bool initialized; 260 SCIP_Bool complete; 261 SCIP_Bool infeasible; 262 SCIP_Real timelimit; 263 264 *result = SCIP_DIDNOTRUN; 265 266 data = SCIPpresolGetData(presol); 267 268 int nvars = SCIPgetNVars(scip); 269 int nconss = SCIPgetNConss(scip); 270 271 /* run only if the problem size reduced by some amount since the last call or if it is the first call */ 272 if( data->lastncols != -1 && data->lastnrows != -1 && 273 nvars > data->lastncols * 0.85 && 274 nconss > data->lastnrows * 0.85 ) 275 return SCIP_OKAY; 276 277 SCIP_CALL( SCIPmatrixCreate(scip, &matrix, TRUE, &initialized, &complete, &infeasible, 278 naddconss, ndelconss, nchgcoefs, nchgbds, nfixedvars) ); 279 280 /* if infeasibility was detected during matrix creation, return here */ 281 if( infeasible ) 282 { 283 if( initialized ) 284 SCIPmatrixFree(scip, &matrix); 285 286 *result = SCIP_CUTOFF; 287 return SCIP_OKAY; 288 } 289 290 /* we only work on pure MIPs, also disable to try building the matrix again if it failed once */ 291 if( !initialized || !complete ) 292 { 293 data->lastncols = 0; 294 data->lastnrows = 0; 295 296 if( initialized ) 297 SCIPmatrixFree(scip, &matrix); 298 299 return SCIP_OKAY; 300 } 301 302 /* only allow communication of constraint modifications by deleting all constraints when they have not been upgraded yet */ 303 SCIP_CONSHDLR* linconshdlr = SCIPfindConshdlr(scip, "linear"); 304 assert(linconshdlr != NULL); 305 bool allowconsmodification = (SCIPconshdlrGetNCheckConss(linconshdlr) == SCIPmatrixGetNRows(matrix)); 306 307 Problem<SCIP_Real> problem = buildProblem(scip, matrix); 308 Presolve<SCIP_Real> presolve; 309 310 /* store current numbers of aggregations, fixings, and changed bounds for statistics */ 311 int oldnaggrvars = *naggrvars; 312 int oldnfixedvars = *nfixedvars; 313 int oldnchgbds = *nchgbds; 314 315 /* important so that SCIP does not throw an error, e.g. when an integer variable is substituted 316 * into a knapsack constraint */ 317 presolve.getPresolveOptions().substitutebinarieswithints = false; 318 319 /* currently these changes cannot be communicated to SCIP correctly since a constraint needs 320 * to be modified in the cases where slackvariables are removed from constraints but for the 321 * presolve library those look like normal substitution on the postsolve stack */ 322 presolve.getPresolveOptions().removeslackvars = false; 323 324 /* communicate the SCIP parameters to the presolve library */ 325 presolve.getPresolveOptions().maxfillinpersubstitution = data->maxfillinpersubstitution; 326 presolve.getPresolveOptions().markowitz_tolerance = data->markowitztolerance; 327 presolve.getPresolveOptions().maxshiftperrow = data->maxshiftperrow; 328 presolve.getPresolveOptions().hugeval = data->hugebound; 329 330 /* removal of linear dependent equations has only an effect when constraint modifications are communicated */ 331 presolve.getPresolveOptions().detectlindep = allowconsmodification ? data->detectlineardependency : 0; 332 333 /* communicate the random seed */ 334 presolve.getPresolveOptions().randomseed = SCIPinitializeRandomSeed(scip, (unsigned int)data->randomseed); 335 336 #ifdef PAPILO_TBB 337 /* set number of threads to be used for presolve */ 338 presolve.getPresolveOptions().threads = data->threads; 339 #else 340 if (data->threads != DEFAULT_THREADS) 341 SCIPverbMessage(scip, SCIP_VERBLEVEL_HIGH, NULL, 342 "PaPILO can utilize only multiple threads if it is build with TBB.\n"); 343 presolve.getPresolveOptions().threads = 1; 344 #endif 345 346 347 /* disable dual reductions that are not permitted */ 348 if( !complete ) 349 presolve.getPresolveOptions().dualreds = 0; 350 else if( SCIPallowStrongDualReds(scip) ) 351 presolve.getPresolveOptions().dualreds = 2; 352 else if( SCIPallowWeakDualReds(scip) ) 353 presolve.getPresolveOptions().dualreds = 1; 354 else 355 presolve.getPresolveOptions().dualreds = 0; 356 357 /* set up the presolvers that shall participate */ 358 using uptr = std::unique_ptr<PresolveMethod<SCIP_Real>>; 359 360 /* fast presolvers*/ 361 presolve.addPresolveMethod( uptr( new SingletonCols<SCIP_Real>() ) ); 362 presolve.addPresolveMethod( uptr( new CoefficientStrengthening<SCIP_Real>() ) ); 363 presolve.addPresolveMethod( uptr( new ConstraintPropagation<SCIP_Real>() ) ); 364 365 /* medium presolver */ 366 presolve.addPresolveMethod( uptr( new SimpleProbing<SCIP_Real>() ) ); 367 if( data->enableparallelrows ) 368 presolve.addPresolveMethod( uptr( new ParallelRowDetection<SCIP_Real>() ) ); 369 /* todo: parallel cols cannot be handled by SCIP currently 370 * addPresolveMethod( uptr( new ParallelColDetection<SCIP_Real>() ) ); */ 371 presolve.addPresolveMethod( uptr( new SingletonStuffing<SCIP_Real>() ) ); 372 #if PAPILO_VERSION_MAJOR > 2 || (PAPILO_VERSION_MAJOR == 2 && PAPILO_VERSION_MINOR >= 1) 373 DualFix<SCIP_Real> *dualfix = new DualFix<SCIP_Real>(); 374 dualfix->set_fix_to_infinity_allowed(false); 375 presolve.addPresolveMethod( uptr( dualfix ) ); 376 #else 377 presolve.addPresolveMethod( uptr( new DualFix<SCIP_Real>() ) ); 378 #endif 379 presolve.addPresolveMethod( uptr( new FixContinuous<SCIP_Real>() ) ); 380 presolve.addPresolveMethod( uptr( new SimplifyInequalities<SCIP_Real>() ) ); 381 presolve.addPresolveMethod( uptr( new SimpleSubstitution<SCIP_Real>() ) ); 382 383 /* exhaustive presolvers*/ 384 presolve.addPresolveMethod( uptr( new ImplIntDetection<SCIP_Real>() ) ); 385 if( data->enabledualinfer ) 386 presolve.addPresolveMethod( uptr( new DualInfer<SCIP_Real>() ) ); 387 if( data->enableprobing ) 388 { 389 #if PAPILO_VERSION_MAJOR > 2 || (PAPILO_VERSION_MAJOR == 2 && PAPILO_VERSION_MINOR >= 1) 390 Probing<SCIP_Real> *probing = new Probing<SCIP_Real>(); 391 if( presolve.getPresolveOptions().runs_sequential() ) 392 { 393 probing->set_max_badge_size( data->maxbadgesizeseq ); 394 } 395 else 396 { 397 probing->set_max_badge_size( data->maxbadgesizepar ); 398 } 399 presolve.addPresolveMethod( uptr( probing ) ); 400 401 #else 402 presolve.addPresolveMethod( uptr( new Probing<SCIP_Real>() ) ); 403 if( data->maxbadgesizeseq != DEFAULT_MAXBADGESIZE_SEQ ) 404 SCIPverbMessage(scip, SCIP_VERBLEVEL_HIGH, NULL, 405 " The parameter 'presolving/milp/maxbadgesizeseq' can only be used with PaPILO 2.1.0 or later versions.\n"); 406 407 if( data->maxbadgesizepar != DEFAULT_MAXBADGESIZE_PAR ) 408 SCIPverbMessage(scip, SCIP_VERBLEVEL_HIGH, NULL, 409 " The parameter 'presolving/milp/maxbadgesizepar' can only be used with PaPILO 2.1.0 or later versions.\n"); 410 411 #endif 412 } 413 if( data->enabledomcol ) 414 presolve.addPresolveMethod( uptr( new DominatedCols<SCIP_Real>() ) ); 415 if( data->enablemultiaggr ) 416 presolve.addPresolveMethod( uptr( new Substitution<SCIP_Real>() ) ); 417 if( data->enablesparsify ) 418 presolve.addPresolveMethod( uptr( new Sparsify<SCIP_Real>() ) ); 419 420 /* set tolerances */ 421 presolve.getPresolveOptions().feastol = SCIPfeastol(scip); 422 presolve.getPresolveOptions().epsilon = SCIPepsilon(scip); 423 424 /* adjust output settings of presolve library */ 425 #ifdef SCIP_PRESOLLIB_ENABLE_OUTPUT 426 problem.setName(SCIPgetProbName(scip)); 427 #else 428 presolve.setVerbosityLevel((VerbosityLevel) data->verbosity); 429 #endif 430 431 /* communicate the time limit */ 432 SCIP_CALL( SCIPgetRealParam(scip, "limits/time", &timelimit) ); 433 if( !SCIPisInfinity(scip, timelimit) ) 434 presolve.getPresolveOptions().tlim = timelimit - SCIPgetSolvingTime(scip); 435 436 if( 0 != strncmp(data->filename, DEFAULT_FILENAME_PROBLEM, strlen(DEFAULT_FILENAME_PROBLEM)) ) 437 { 438 SCIP_CALL(SCIPwriteTransProblem(scip, data->filename, NULL, FALSE)); 439 } 440 441 /* call the presolving */ 442 SCIPverbMessage(scip, SCIP_VERBLEVEL_HIGH, NULL, 443 " (%.1fs) running MILP presolver\n", SCIPgetSolvingTime(scip)); 444 int oldnnz = problem.getConstraintMatrix().getNnz(); 445 446 /*call presolving without storing information for dual postsolve*/ 447 #if (PAPILO_VERSION_MAJOR >= 2) 448 PresolveResult<SCIP_Real> res = presolve.apply(problem, false); 449 #else 450 PresolveResult<SCIP_Real> res = presolve.apply(problem); 451 #endif 452 data->lastncols = problem.getNCols(); 453 data->lastnrows = problem.getNRows(); 454 455 /* evaluate the result */ 456 switch(res.status) 457 { 458 case PresolveStatus::kInfeasible: 459 *result = SCIP_CUTOFF; 460 SCIPverbMessage(scip, SCIP_VERBLEVEL_HIGH, NULL, 461 " (%.1fs) MILP presolver detected infeasibility\n", 462 SCIPgetSolvingTime(scip)); 463 SCIPmatrixFree(scip, &matrix); 464 return SCIP_OKAY; 465 case PresolveStatus::kUnbndOrInfeas: 466 case PresolveStatus::kUnbounded: 467 *result = SCIP_UNBOUNDED; 468 SCIPverbMessage(scip, SCIP_VERBLEVEL_HIGH, NULL, 469 " (%.1fs) MILP presolver detected unboundedness\n", 470 SCIPgetSolvingTime(scip)); 471 SCIPmatrixFree(scip, &matrix); 472 return SCIP_OKAY; 473 case PresolveStatus::kUnchanged: 474 *result = SCIP_DIDNOTFIND; 475 data->lastncols = nvars; 476 data->lastnrows = nconss; 477 SCIPverbMessage(scip, SCIP_VERBLEVEL_HIGH, NULL, 478 " (%.1fs) MILP presolver found nothing\n", 479 SCIPgetSolvingTime(scip)); 480 SCIPmatrixFree(scip, &matrix); 481 return SCIP_OKAY; 482 case PresolveStatus::kReduced: 483 data->lastncols = problem.getNCols(); 484 data->lastnrows = problem.getNRows(); 485 *result = SCIP_SUCCESS; 486 } 487 488 /* result indicated success, now populate the changes into the SCIP structures */ 489 std::vector<SCIP_VAR*> tmpvars; 490 std::vector<SCIP_Real> tmpvals; 491 492 /* if the number of nonzeros decreased by a sufficient factor, rather create all constraints from scratch */ 493 int newnnz = problem.getConstraintMatrix().getNnz(); 494 bool constraintsReplaced = false; 495 if( newnnz == 0 || (allowconsmodification && 496 (problem.getNRows() <= data->modifyconsfac * data->lastnrows || 497 newnnz <= data->modifyconsfac * oldnnz)) ) 498 { 499 int oldnrows = SCIPmatrixGetNRows(matrix); 500 int newnrows = problem.getNRows(); 501 502 constraintsReplaced = true; 503 504 /* capture constraints that are still present in the problem after presolve */ 505 for( int i = 0; i < newnrows; ++i ) 506 { 507 SCIP_CONS* c = SCIPmatrixGetCons(matrix, res.postsolve.origrow_mapping[i]); 508 SCIP_CALL( SCIPcaptureCons(scip, c) ); 509 } 510 511 /* delete all constraints */ 512 *ndelconss += oldnrows; 513 *naddconss += newnrows; 514 515 for( int i = 0; i < oldnrows; ++i ) 516 { 517 SCIP_CALL( SCIPdelCons(scip, SCIPmatrixGetCons(matrix, i)) ); 518 } 519 520 /* now loop over rows of presolved problem and create them as new linear constraints, 521 * then release the old constraint after its name was passed to the new constraint */ 522 const Vec<RowFlags>& rflags = problem.getRowFlags(); 523 const auto& consmatrix = problem.getConstraintMatrix(); 524 for( int i = 0; i < newnrows; ++i ) 525 { 526 auto rowvec = consmatrix.getRowCoefficients(i); 527 const int* rowcols = rowvec.getIndices(); 528 /* SCIPcreateConsBasicLinear() requires a non const pointer */ 529 SCIP_Real* rowvals = const_cast<SCIP_Real*>(rowvec.getValues()); 530 int rowlen = rowvec.getLength(); 531 532 /* retrieve SCIP compatible left and right hand sides */ 533 SCIP_Real lhs = rflags[i].test(RowFlag::kLhsInf) ? - SCIPinfinity(scip) : consmatrix.getLeftHandSides()[i]; 534 SCIP_Real rhs = rflags[i].test(RowFlag::kRhsInf) ? SCIPinfinity(scip) : consmatrix.getRightHandSides()[i]; 535 536 /* create variable array matching the value array */ 537 tmpvars.clear(); 538 tmpvars.reserve(rowlen); 539 for( int j = 0; j < rowlen; ++j ) 540 tmpvars.push_back(SCIPmatrixGetVar(matrix, res.postsolve.origcol_mapping[rowcols[j]])); 541 542 /* create and add new constraint with name of old constraint */ 543 SCIP_CONS* oldcons = SCIPmatrixGetCons(matrix, res.postsolve.origrow_mapping[i]); 544 SCIP_CONS* cons; 545 SCIP_CALL( SCIPcreateConsBasicLinear(scip, &cons, SCIPconsGetName(oldcons), rowlen, tmpvars.data(), rowvals, lhs, rhs) ); 546 SCIP_CALL( SCIPaddCons(scip, cons) ); 547 548 /* release old and new constraint */ 549 SCIP_CALL( SCIPreleaseCons(scip, &oldcons) ); 550 SCIP_CALL( SCIPreleaseCons(scip, &cons) ); 551 } 552 } 553 554 /* loop over res.postsolve and add all fixed variables and aggregations to scip */ 555 for( std::size_t i = 0; i != res.postsolve.types.size(); ++i ) 556 { 557 ReductionType type = res.postsolve.types[i]; 558 int first = res.postsolve.start[i]; 559 int last = res.postsolve.start[i + 1]; 560 561 switch( type ) 562 { 563 case ReductionType::kFixedCol: 564 { 565 SCIP_Bool infeas; 566 SCIP_Bool fixed; 567 int col = res.postsolve.indices[first]; 568 569 SCIP_VAR* var = SCIPmatrixGetVar(matrix, col); 570 571 SCIP_Real value = res.postsolve.values[first]; 572 573 SCIP_CALL( SCIPfixVar(scip, var, value, &infeas, &fixed) ); 574 *nfixedvars += 1; 575 576 assert(!infeas); 577 /* SCIP has different rules for aggregating variables than PaPILO: therefore the variable PaPILO 578 * tries to fix now may have been aggregated by SCIP before. Additionally, after aggregation SCIP 579 * sometimes performs bound tightening resulting in possible fixings. These cases need to be excluded. */ 580 assert(fixed || SCIPvarGetStatus(var) == SCIP_VARSTATUS_AGGREGATED || 581 SCIPvarGetStatus(var) == SCIP_VARSTATUS_FIXED || SCIPvarGetStatus(var) == SCIP_VARSTATUS_NEGATED); 582 break; 583 } 584 /* 585 * Dual-postsolving in PaPILO required introducing a postsolve-type for substitution with additional information. 586 * Further, the different Substitution-postsolving types store the required postsolving data differently (in different order) in the postsolving stack. 587 * Therefore, we need to distinguish how to parse the required data (rowLength, col, side, startRowCoefficients, lastRowCoefficients) from the postsolving stack. 588 * If these values are accessed, the procedure is the same for both. 589 */ 590 #if (PAPILO_VERSION_MAJOR >= 2) 591 case ReductionType::kSubstitutedColWithDual: 592 #endif 593 case ReductionType::kSubstitutedCol: 594 { 595 int col = 0; 596 SCIP_Real side = 0; 597 598 int rowlen = 0; 599 int startRowCoefficients = 0; 600 int lastRowCoefficients = 0; 601 602 if( type == ReductionType::kSubstitutedCol ) 603 { 604 rowlen = last - first - 1; 605 col = res.postsolve.indices[first]; 606 side = res.postsolve.values[first]; 607 608 startRowCoefficients = first + 1; 609 lastRowCoefficients = last; 610 } 611 #if (PAPILO_VERSION_MAJOR >= 2) 612 if( type == ReductionType::kSubstitutedColWithDual ) 613 { 614 rowlen = (int) res.postsolve.values[first]; 615 col = res.postsolve.indices[first + 3 + rowlen]; 616 side = res.postsolve.values[first + 1]; 617 618 startRowCoefficients = first + 3; 619 lastRowCoefficients = first + 3 + rowlen; 620 621 assert(side == res.postsolve.values[first + 2]); 622 assert(res.postsolve.indices[first + 1] == 0); 623 assert(res.postsolve.indices[first + 2] == 0); 624 } 625 assert( type == ReductionType::kSubstitutedCol || type == ReductionType::kSubstitutedColWithDual ); 626 #else 627 assert( type == ReductionType::kSubstitutedCol ); 628 #endif 629 SCIP_Bool infeas; 630 SCIP_Bool aggregated; 631 SCIP_Bool redundant = FALSE; 632 SCIP_Real constant = 0.0; 633 if( rowlen == 2 ) 634 { 635 SCIP_Real updatedSide; 636 SCIP_VAR* varx = SCIPmatrixGetVar(matrix, res.postsolve.indices[startRowCoefficients]); 637 SCIP_VAR* vary = SCIPmatrixGetVar(matrix, res.postsolve.indices[startRowCoefficients + 1]); 638 SCIP_Real scalarx = res.postsolve.values[startRowCoefficients]; 639 SCIP_Real scalary = res.postsolve.values[startRowCoefficients + 1]; 640 641 SCIP_CALL( SCIPgetProbvarSum(scip, &varx, &scalarx, &constant) ); 642 assert(SCIPvarGetStatus(varx) != SCIP_VARSTATUS_MULTAGGR); 643 644 SCIP_CALL( SCIPgetProbvarSum(scip, &vary, &scalary, &constant) ); 645 assert(SCIPvarGetStatus(vary) != SCIP_VARSTATUS_MULTAGGR); 646 647 updatedSide = side - constant; 648 649 SCIP_CALL( SCIPaggregateVars(scip, varx, vary, scalarx, scalary, updatedSide, &infeas, &redundant, &aggregated) ); 650 } 651 else 652 { 653 SCIP_Real colCoef = 0.0; 654 SCIP_Real updatedSide; 655 656 for( int j = startRowCoefficients; j < lastRowCoefficients; ++j ) 657 { 658 if( res.postsolve.indices[j] == col ) 659 { 660 colCoef = res.postsolve.values[j]; 661 break; 662 } 663 } 664 665 tmpvars.clear(); 666 tmpvals.clear(); 667 tmpvars.reserve(rowlen); 668 tmpvals.reserve(rowlen); 669 670 assert(colCoef != 0.0); 671 SCIP_VAR* aggrvar = SCIPmatrixGetVar(matrix, col); 672 673 SCIP_CALL( SCIPgetProbvarSum(scip, &aggrvar, &colCoef, &constant) ); 674 assert(SCIPvarGetStatus(aggrvar) != SCIP_VARSTATUS_MULTAGGR); 675 676 updatedSide = side - constant; 677 678 for( int j = startRowCoefficients; j < lastRowCoefficients; ++j ) 679 { 680 if( res.postsolve.indices[j] == col ) 681 continue; 682 683 tmpvars.push_back(SCIPmatrixGetVar(matrix, res.postsolve.indices[j])); 684 tmpvals.push_back(- res.postsolve.values[j] / colCoef); 685 } 686 687 SCIP_CALL( SCIPmultiaggregateVar(scip, aggrvar, tmpvars.size(), 688 tmpvars.data(), tmpvals.data(), updatedSide / colCoef, &infeas, &aggregated) ); 689 } 690 691 if( aggregated ) 692 *naggrvars += 1; 693 else if( constraintsReplaced && !redundant ) 694 { 695 /* if the constraints where replaced, we need to add the failed substitution as an equality to SCIP */ 696 tmpvars.clear(); 697 tmpvals.clear(); 698 for( int j = startRowCoefficients; j < lastRowCoefficients; ++j ) 699 { 700 tmpvars.push_back(SCIPmatrixGetVar(matrix, res.postsolve.indices[j])); 701 tmpvals.push_back(res.postsolve.values[j]); 702 } 703 704 SCIP_CONS* cons; 705 String name = fmt::format("{}_failed_aggregation_equality", SCIPvarGetName(SCIPmatrixGetVar(matrix, col))); 706 SCIP_CALL( SCIPcreateConsBasicLinear(scip, &cons, name.c_str(), 707 tmpvars.size(), tmpvars.data(), tmpvals.data(), side, side ) ); 708 SCIP_CALL( SCIPaddCons(scip, cons) ); 709 SCIP_CALL( SCIPreleaseCons(scip, &cons) ); 710 *naddconss += 1; 711 } 712 713 if( infeas ) 714 { 715 *result = SCIP_CUTOFF; 716 break; 717 } 718 719 break; 720 } 721 case ReductionType::kParallelCol: 722 return SCIP_INVALIDRESULT; 723 #if (PAPILO_VERSION_MAJOR <= 1 && PAPILO_VERSION_MINOR==0) 724 #else 725 case ReductionType::kFixedInfCol: { 726 /* todo: currently SCIP can not handle this kind of reduction (see issue #3391) */ 727 assert(false); 728 if(!constraintsReplaced) 729 continue; 730 SCIP_Bool infeas; 731 SCIP_Bool fixed; 732 SCIP_Real value = SCIPinfinity(scip); 733 734 int column = res.postsolve.indices[first]; 735 bool is_negative_infinity = res.postsolve.values[first] < 0; 736 SCIP_VAR* column_variable = SCIPmatrixGetVar(matrix, column); 737 738 if( is_negative_infinity ) 739 { 740 value = -SCIPinfinity(scip); 741 } 742 743 SCIP_CALL( SCIPfixVar(scip, column_variable, value, &infeas, &fixed) ); 744 *nfixedvars += 1; 745 746 assert(!infeas); 747 assert(fixed); 748 break; 749 } 750 #endif 751 #if (PAPILO_VERSION_MAJOR >= 2) 752 case ReductionType::kVarBoundChange : 753 case ReductionType::kRedundantRow : 754 case ReductionType::kRowBoundChange : 755 case ReductionType::kReasonForRowBoundChangeForcedByRow : 756 case ReductionType::kRowBoundChangeForcedByRow : 757 case ReductionType::kSaveRow : 758 case ReductionType::kReducedBoundsCost : 759 case ReductionType::kColumnDualValue : 760 case ReductionType::kRowDualValue : 761 case ReductionType::kCoefficientChange : 762 // dual ReductionTypes should be only calculated for dual reductions and should not appear for MIP 763 SCIPerrorMessage("PaPILO: PaPILO should not return dual postsolving reductions in SCIP!!\n"); 764 SCIPABORT(); /*lint --e{527}*/ 765 break; 766 #endif 767 default: 768 SCIPdebugMsg(scip, "PaPILO returned unknown data type: \n" ); 769 continue; 770 } 771 } 772 773 /* tighten bounds of variables that are still present after presolving */ 774 if( *result != SCIP_CUTOFF ) 775 { 776 VariableDomains<SCIP_Real>& varDomains = problem.getVariableDomains(); 777 for( int i = 0; i != problem.getNCols(); ++i ) 778 { 779 assert( ! varDomains.flags[i].test(ColFlag::kInactive) ); 780 SCIP_VAR* var = SCIPmatrixGetVar(matrix, res.postsolve.origcol_mapping[i]); 781 if( !varDomains.flags[i].test(ColFlag::kLbInf) ) 782 { 783 SCIP_Bool infeas; 784 SCIP_Bool tightened; 785 SCIP_CALL( SCIPtightenVarLb(scip, var, varDomains.lower_bounds[i], TRUE, &infeas, &tightened) ); 786 787 if( tightened ) 788 *nchgbds += 1; 789 790 if( infeas ) 791 { 792 *result = SCIP_CUTOFF; 793 break; 794 } 795 } 796 797 if( !varDomains.flags[i].test(ColFlag::kUbInf) ) 798 { 799 SCIP_Bool infeas; 800 SCIP_Bool tightened; 801 SCIP_CALL( SCIPtightenVarUb(scip, var, varDomains.upper_bounds[i], TRUE, &infeas, &tightened) ); 802 803 if( tightened ) 804 *nchgbds += 1; 805 806 if( infeas ) 807 { 808 *result = SCIP_CUTOFF; 809 break; 810 } 811 } 812 } 813 } 814 815 /* finish with a final verb message and return */ 816 SCIPverbMessage(scip, SCIP_VERBLEVEL_HIGH, NULL, 817 " (%.1fs) MILP presolver (%d rounds): %d aggregations, %d fixings, %d bound changes\n", 818 SCIPgetSolvingTime(scip), presolve.getStatistics().nrounds, *naggrvars - oldnaggrvars, 819 *nfixedvars - oldnfixedvars, *nchgbds - oldnchgbds); 820 821 /* free the matrix */ 822 assert(initialized); 823 SCIPmatrixFree(scip, &matrix); 824 825 return SCIP_OKAY; 826 } 827 828 829 /* 830 * presolver specific interface methods 831 */ 832 833 /** creates the xyz presolver and includes it in SCIP */ 834 SCIP_RETCODE SCIPincludePresolMILP( 835 SCIP* scip /**< SCIP data structure */ 836 ) 837 { 838 SCIP_PRESOLDATA* presoldata; 839 SCIP_PRESOL* presol; 840 841 #if defined(PAPILO_VERSION_TWEAK) && PAPILO_VERSION_TWEAK != 0 842 String name = fmt::format("PaPILO {}.{}.{}.{}", PAPILO_VERSION_MAJOR, PAPILO_VERSION_MINOR, PAPILO_VERSION_PATCH, PAPILO_VERSION_TWEAK); 843 #else 844 String name = fmt::format("PaPILO {}.{}.{}", PAPILO_VERSION_MAJOR, PAPILO_VERSION_MINOR, PAPILO_VERSION_PATCH); 845 #endif 846 847 #if defined(PAPILO_GITHASH_AVAILABLE) && defined(PAPILO_TBB) 848 String desc = fmt::format("parallel presolve for integer and linear optimization (github.com/scipopt/papilo) (built with TBB) [GitHash: {}]", PAPILO_GITHASH); 849 #elif !defined(PAPILO_GITHASH_AVAILABLE) && !defined(PAPILO_TBB) 850 String desc("parallel presolve for integer and linear optimization (github.com/scipopt/papilo)"); 851 #elif defined(PAPILO_GITHASH_AVAILABLE) && !defined(PAPILO_TBB) 852 String desc = fmt::format("parallel presolve for integer and linear optimization (github.com/scipopt/papilo) [GitHash: {}]", PAPILO_GITHASH); 853 #elif !defined(PAPILO_GITHASH_AVAILABLE) && defined(PAPILO_TBB) 854 String desc = fmt::format("parallel presolve for integer and linear optimization (github.com/scipopt/papilo) (built with TBB)"); 855 #endif 856 857 /* add external code info for the presolve library */ 858 SCIP_CALL( SCIPincludeExternalCodeInformation(scip, name.c_str(), desc.c_str()) ); 859 860 /* create MILP presolver data */ 861 presoldata = NULL; 862 SCIP_CALL( SCIPallocBlockMemory(scip, &presoldata) ); 863 BMSclearMemory(presoldata); 864 865 presol = NULL; 866 867 /* include presolver */ 868 SCIP_CALL( SCIPincludePresolBasic(scip, &presol, PRESOL_NAME, PRESOL_DESC, PRESOL_PRIORITY, PRESOL_MAXROUNDS, PRESOL_TIMING, 869 presolExecMILP, 870 presoldata) ); 871 872 assert(presol != NULL); 873 874 /* set non fundamental callbacks via setter functions */ 875 SCIP_CALL( SCIPsetPresolCopy(scip, presol, presolCopyMILP) ); 876 SCIP_CALL( SCIPsetPresolFree(scip, presol, presolFreeMILP) ); 877 SCIP_CALL( SCIPsetPresolInit(scip, presol, presolInitMILP) ); 878 879 /* add MILP presolver parameters */ 880 SCIP_CALL( SCIPaddIntParam(scip, 881 "presolving/" PRESOL_NAME "/threads", 882 "maximum number of threads presolving may use (0: automatic)", 883 &presoldata->threads, FALSE, DEFAULT_THREADS, 0, INT_MAX, NULL, NULL) ); 884 885 SCIP_CALL( SCIPaddIntParam(scip, 886 "presolving/" PRESOL_NAME "/maxfillinpersubstitution", 887 "maximal possible fillin for substitutions to be considered", 888 &presoldata->maxfillinpersubstitution, FALSE, DEFAULT_MAXFILLINPERSUBST, INT_MIN, INT_MAX, NULL, NULL) ); 889 890 SCIP_CALL( SCIPaddIntParam(scip, 891 "presolving/" PRESOL_NAME "/maxshiftperrow", 892 "maximal amount of nonzeros allowed to be shifted to make space for substitutions", 893 &presoldata->maxshiftperrow, TRUE, DEFAULT_MAXSHIFTPERROW, 0, INT_MAX, NULL, NULL) ); 894 895 SCIP_CALL( SCIPaddIntParam(scip, 896 "presolving/" PRESOL_NAME "/randomseed", 897 "the random seed used for randomization of tie breaking", 898 &presoldata->randomseed, FALSE, DEFAULT_RANDOMSEED, INT_MIN, INT_MAX, NULL, NULL) ); 899 900 if( DependentRows<double>::Enabled ) 901 { 902 SCIP_CALL( SCIPaddIntParam(scip, 903 "presolving/" PRESOL_NAME "/detectlineardependency", 904 "should linear dependent equations and free columns be removed? (0: never, 1: for LPs, 2: always)", 905 &presoldata->detectlineardependency, TRUE, DEFAULT_DETECTLINDEP, 0, 2, NULL, NULL) ); 906 } 907 else 908 presoldata->detectlineardependency = 0; 909 910 SCIP_CALL( SCIPaddRealParam(scip, 911 "presolving/" PRESOL_NAME "/modifyconsfac", 912 "modify SCIP constraints when the number of nonzeros or rows is at most this factor " 913 "times the number of nonzeros or rows before presolving", 914 &presoldata->modifyconsfac, FALSE, DEFAULT_MODIFYCONSFAC, 0.0, 1.0, NULL, NULL) ); 915 916 SCIP_CALL( SCIPaddRealParam(scip, 917 "presolving/" PRESOL_NAME "/markowitztolerance", 918 "the markowitz tolerance used for substitutions", 919 &presoldata->markowitztolerance, FALSE, DEFAULT_MARKOWITZTOLERANCE, 0.0, 1.0, NULL, NULL) ); 920 921 SCIP_CALL( SCIPaddRealParam(scip, 922 "presolving/" PRESOL_NAME "/hugebound", 923 "absolute bound value that is considered too huge for activitity based calculations", 924 &presoldata->hugebound, FALSE, DEFAULT_HUGEBOUND, 0.0, SCIP_REAL_MAX, NULL, NULL) ); 925 926 #if PAPILO_VERSION_MAJOR > 2 || (PAPILO_VERSION_MAJOR == 2 && PAPILO_VERSION_MINOR >= 1) 927 SCIP_CALL(SCIPaddIntParam(scip, "presolving/" PRESOL_NAME "/maxbadgesizeseq", 928 "maximal badge size in Probing in PaPILO if PaPILO is executed in sequential mode", 929 &presoldata->maxbadgesizeseq, FALSE, DEFAULT_MAXBADGESIZE_SEQ, -1, INT_MAX, NULL, NULL)); 930 931 SCIP_CALL(SCIPaddIntParam(scip, "presolving/" PRESOL_NAME "/maxbadgesizepar", 932 "maximal badge size in Probing in PaPILO if PaPILO is executed in parallel mode", 933 &presoldata->maxbadgesizepar, FALSE, DEFAULT_MAXBADGESIZE_PAR, -1, INT_MAX, NULL, NULL)); 934 #else 935 presoldata->maxbadgesizeseq = DEFAULT_MAXBADGESIZE_SEQ; 936 presoldata->maxbadgesizepar = DEFAULT_MAXBADGESIZE_PAR; 937 #endif 938 939 SCIP_CALL( SCIPaddBoolParam(scip, 940 "presolving/" PRESOL_NAME "/enableparallelrows", 941 "should the parallel rows presolver be enabled within the presolve library?", 942 &presoldata->enableparallelrows, TRUE, DEFAULT_ENABLEPARALLELROWS, NULL, NULL) ); 943 944 SCIP_CALL( SCIPaddBoolParam(scip, 945 "presolving/" PRESOL_NAME "/enabledomcol", 946 "should the dominated column presolver be enabled within the presolve library?", 947 &presoldata->enabledomcol, TRUE, DEFAULT_ENABLEDOMCOL, NULL, NULL) ); 948 949 SCIP_CALL( SCIPaddBoolParam(scip, 950 "presolving/" PRESOL_NAME "/enabledualinfer", 951 "should the dualinfer presolver be enabled within the presolve library?", 952 &presoldata->enabledualinfer, TRUE, DEFAULT_ENABLEDUALINFER, NULL, NULL) ); 953 954 SCIP_CALL( SCIPaddBoolParam(scip, 955 "presolving/" PRESOL_NAME "/enablemultiaggr", 956 "should the multi-aggregation presolver be enabled within the presolve library?", 957 &presoldata->enablemultiaggr, TRUE, DEFAULT_ENABLEMULTIAGGR, NULL, NULL) ); 958 959 SCIP_CALL( SCIPaddBoolParam(scip, 960 "presolving/" PRESOL_NAME "/enableprobing", 961 "should the probing presolver be enabled within the presolve library?", 962 &presoldata->enableprobing, TRUE, DEFAULT_ENABLEPROBING, NULL, NULL) ); 963 964 SCIP_CALL( SCIPaddBoolParam(scip, 965 "presolving/" PRESOL_NAME "/enablesparsify", 966 "should the sparsify presolver be enabled within the presolve library?", 967 &presoldata->enablesparsify, TRUE, DEFAULT_ENABLESPARSIFY, NULL, NULL) ); 968 969 SCIP_CALL( SCIPaddStringParam(scip, "presolving/" PRESOL_NAME "/probfilename", 970 "filename to store the problem before MILP presolving starts", 971 &presoldata->filename, TRUE, DEFAULT_FILENAME_PROBLEM, NULL, NULL) ); 972 973 SCIP_CALL(SCIPaddIntParam(scip, "presolving/" PRESOL_NAME "/verbosity", 974 "verbosity level of PaPILO (0: quiet, 1: errors, 2: warnings, 3: normal, 4: detailed)", 975 &presoldata->verbosity, FALSE, DEFAULT_VERBOSITY, 0, 4, NULL, NULL)); 976 977 return SCIP_OKAY; 978 } 979 980 #endif 981