1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 2 /* */ 3 /* This file is part of the program and library */ 4 /* SCIP --- Solving Constraint Integer Programs */ 5 /* */ 6 /* Copyright 2002-2022 Zuse Institute Berlin */ 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 conflict_dualproofanalysis.c 26 * @ingroup OTHER_CFILES 27 * @brief internal methods for dual proof conflict analysis 28 * @author Timo Berthold 29 * @author Jakob Witzig 30 * 31 * In dual proof analysis, an infeasible LP relaxation is analysed. 32 * Using the dual solution, a valid constraint is derived that is violated 33 * by all values in the domain. This constraint is added to the problem 34 * and can then be used for domain propagation. More details can be found in [1] 35 * 36 * [1] J. Witzig, T. Berthold, en S. Heinz, ‘Computational aspects of infeasibility analysis in mixed integer programming’, 37 * Math. Prog. Comp., mrt. 2021, doi: 10.1007/s12532-021-00202-0. 38 */ 39 40 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/ 41 42 #include "lpi/lpi.h" 43 #include "scip/clock.h" 44 #include "scip/conflict_general.h" 45 #include "scip/conflict_dualproofanalysis.h" 46 #include "scip/conflictstore.h" 47 #include "scip/cons.h" 48 #include "scip/cons_linear.h" 49 #include "scip/cuts.h" 50 #include "scip/history.h" 51 #include "scip/lp.h" 52 #include "scip/presolve.h" 53 #include "scip/prob.h" 54 #include "scip/prop.h" 55 #include "scip/pub_conflict.h" 56 #include "scip/pub_cons.h" 57 #include "scip/pub_lp.h" 58 #include "scip/pub_message.h" 59 #include "scip/pub_misc.h" 60 #include "scip/pub_misc_sort.h" 61 #include "scip/pub_paramset.h" 62 #include "scip/pub_prop.h" 63 #include "scip/pub_tree.h" 64 #include "scip/pub_var.h" 65 #include "scip/scip_conflict.h" 66 #include "scip/scip_cons.h" 67 #include "scip/scip_mem.h" 68 #include "scip/scip_sol.h" 69 #include "scip/scip_var.h" 70 #include "scip/set.h" 71 #include "scip/sol.h" 72 #include "scip/struct_conflict.h" 73 #include "scip/struct_lp.h" 74 #include "scip/struct_prob.h" 75 #include "scip/struct_set.h" 76 #include "scip/struct_stat.h" 77 #include "scip/struct_tree.h" 78 #include "scip/struct_var.h" 79 #include "scip/tree.h" 80 #include "scip/var.h" 81 #include "scip/visual.h" 82 83 #define BOUNDSWITCH 0.51 /**< threshold for bound switching - see cuts.c */ 84 #define POSTPROCESS FALSE /**< apply postprocessing to the cut - see cuts.c */ 85 #define USEVBDS FALSE /**< use variable bounds - see cuts.c */ 86 #define ALLOWLOCAL FALSE /**< allow to generate local cuts - see cuts. */ 87 #define MINFRAC 0.05 /**< minimal fractionality of floor(rhs) - see cuts.c */ 88 #define MAXFRAC 0.999 /**< maximal fractionality of floor(rhs) - see cuts.c */ 89 90 91 /* 92 * Proof Sets 93 */ 94 95 /** return the char associated with the type of the variable */ 96 static 97 char varGetChar( 98 SCIP_VAR* var /**< variable */ 99 ) 100 { 101 SCIP_VARTYPE vartype = SCIPvarGetType(var); 102 103 return (!SCIPvarIsIntegral(var) ? 'C' : 104 (vartype == SCIP_VARTYPE_BINARY ? 'B' : 105 (vartype == SCIP_VARTYPE_INTEGER ? 'I' : 'M'))); 106 } 107 108 /** resets the data structure of a proofset */ 109 static 110 void proofsetClear( 111 SCIP_PROOFSET* proofset /**< proof set */ 112 ) 113 { 114 assert(proofset != NULL); 115 116 proofset->nnz = 0; 117 proofset->rhs = 0.0; 118 proofset->validdepth = 0; 119 proofset->conflicttype = SCIP_CONFTYPE_UNKNOWN; 120 } 121 122 /** creates a proofset */ 123 static 124 SCIP_RETCODE proofsetCreate( 125 SCIP_PROOFSET** proofset, /**< proof set */ 126 BMS_BLKMEM* blkmem /**< block memory of transformed problem */ 127 ) 128 { 129 assert(proofset != NULL); 130 131 SCIP_ALLOC( BMSallocBlockMemory(blkmem, proofset) ); 132 (*proofset)->vals = NULL; 133 (*proofset)->inds = NULL; 134 (*proofset)->rhs = 0.0; 135 (*proofset)->nnz = 0; 136 (*proofset)->size = 0; 137 (*proofset)->validdepth = 0; 138 (*proofset)->conflicttype = SCIP_CONFTYPE_UNKNOWN; 139 140 return SCIP_OKAY; 141 } 142 143 /** frees a proofset */ 144 void SCIPproofsetFree( 145 SCIP_PROOFSET** proofset, /**< proof set */ 146 BMS_BLKMEM* blkmem /**< block memory */ 147 ) 148 { 149 assert(proofset != NULL); 150 assert(*proofset != NULL); 151 assert(blkmem != NULL); 152 153 BMSfreeBlockMemoryArrayNull(blkmem, &(*proofset)->vals, (*proofset)->size); 154 BMSfreeBlockMemoryArrayNull(blkmem, &(*proofset)->inds, (*proofset)->size); 155 BMSfreeBlockMemory(blkmem, proofset); 156 (*proofset) = NULL; 157 } 158 159 #ifdef SCIP_DEBUG 160 /** print a proof set */ 161 void proofsetPrint( 162 SCIP_PROOFSET* proofset, 163 SCIP_SET* set, 164 SCIP_PROB* transprob 165 ) 166 { 167 SCIP_VAR** vars; 168 int i; 169 170 assert(proofset != NULL); 171 172 vars = SCIPprobGetVars(transprob); 173 assert(vars != NULL); 174 175 printf("proofset: "); 176 for( i = 0; i < proofset->nnz; i++ ) 177 printf("%+.15g <%s> ", proofset->vals[i], SCIPvarGetName(vars[proofset->inds[i]])); 178 printf(" <= %.15g\n", proofset->rhs); 179 } 180 #endif 181 182 /** return the indices of variables in the proofset */ 183 static 184 int* proofsetGetInds( 185 SCIP_PROOFSET* proofset /**< proof set */ 186 ) 187 { 188 assert(proofset != NULL); 189 190 return proofset->inds; 191 } 192 193 /** return coefficient of variable in the proofset with given probindex */ 194 static 195 SCIP_Real* proofsetGetVals( 196 SCIP_PROOFSET* proofset /**< proof set */ 197 ) 198 { 199 assert(proofset != NULL); 200 201 return proofset->vals; 202 } 203 204 /** return the right-hand side if a proofset */ 205 static 206 SCIP_Real proofsetGetRhs( 207 SCIP_PROOFSET* proofset /**< proof set */ 208 ) 209 { 210 assert(proofset != NULL); 211 212 return proofset->rhs; 213 } 214 215 /** returns the number of variables in the proofset */ 216 int SCIPproofsetGetNVars( 217 SCIP_PROOFSET* proofset /**< proof set */ 218 ) 219 { 220 assert(proofset != NULL); 221 222 return proofset->nnz; 223 } 224 225 /** returns the number of variables in the proofset */ 226 static 227 SCIP_CONFTYPE proofsetGetConftype( 228 SCIP_PROOFSET* proofset /**< proof set */ 229 ) 230 { 231 assert(proofset != NULL); 232 233 return proofset->conflicttype; 234 } 235 236 /** adds given data as aggregation row to the proofset */ 237 static 238 SCIP_RETCODE proofsetAddSparseData( 239 SCIP_PROOFSET* proofset, /**< proof set */ 240 BMS_BLKMEM* blkmem, /**< block memory */ 241 SCIP_Real* vals, /**< variable coefficients */ 242 int* inds, /**< variable array */ 243 int nnz, /**< size of variable and coefficient array */ 244 SCIP_Real rhs /**< right-hand side of the aggregation row */ 245 ) 246 { 247 assert(proofset != NULL); 248 assert(blkmem != NULL); 249 250 if( proofset->size == 0 ) 251 { 252 assert(proofset->vals == NULL); 253 assert(proofset->inds == NULL); 254 255 SCIP_ALLOC( BMSduplicateBlockMemoryArray(blkmem, &proofset->vals, vals, nnz) ); 256 SCIP_ALLOC( BMSduplicateBlockMemoryArray(blkmem, &proofset->inds, inds, nnz) ); 257 258 proofset->size = nnz; 259 } 260 else 261 { 262 int i; 263 264 assert(proofset->vals != NULL); 265 assert(proofset->inds != NULL); 266 267 if( proofset->size < nnz ) 268 { 269 SCIP_ALLOC( BMSreallocBlockMemoryArray(blkmem, &proofset->vals, proofset->size, nnz) ); 270 SCIP_ALLOC( BMSreallocBlockMemoryArray(blkmem, &proofset->inds, proofset->size, nnz) ); 271 proofset->size = nnz; 272 } 273 274 for( i = 0; i < nnz; i++ ) 275 { 276 proofset->vals[i] = vals[i]; 277 proofset->inds[i] = inds[i]; 278 } 279 } 280 281 proofset->rhs = rhs; 282 proofset->nnz = nnz; 283 284 return SCIP_OKAY; 285 } 286 287 /** adds an aggregation row to the proofset */ 288 static 289 SCIP_RETCODE proofsetAddAggrrow( 290 SCIP_PROOFSET* proofset, /**< proof set */ 291 SCIP_SET* set, /**< global SCIP settings */ 292 BMS_BLKMEM* blkmem, /**< block memory */ 293 SCIP_AGGRROW* aggrrow /**< aggregation row to add */ 294 ) 295 { 296 SCIP_Real* vals; 297 int* inds; 298 int nnz; 299 int i; 300 301 assert(proofset != NULL); 302 assert(set != NULL); 303 304 inds = SCIPaggrRowGetInds(aggrrow); 305 assert(inds != NULL); 306 307 nnz = SCIPaggrRowGetNNz(aggrrow); 308 assert(nnz > 0); 309 310 SCIP_CALL( SCIPsetAllocBufferArray(set, &vals, nnz) ); 311 312 for( i = 0; i < nnz; i++ ) 313 { 314 vals[i] = SCIPaggrRowGetProbvarValue(aggrrow, inds[i]); 315 } 316 317 SCIP_CALL( proofsetAddSparseData(proofset, blkmem, vals, inds, nnz, SCIPaggrRowGetRhs(aggrrow)) ); 318 319 SCIPsetFreeBufferArray(set, &vals); 320 321 return SCIP_OKAY; 322 } 323 324 /** Removes a given variable @p var from position @p pos from the proofset and updates the right-hand side according 325 * to sign of the coefficient, i.e., rhs -= coef * bound, where bound = lb if coef >= 0 and bound = ub, otherwise. 326 * 327 * @note: The list of non-zero indices and coefficients will be updated by swapping the last non-zero index to @p pos. 328 */ 329 static 330 void proofsetCancelVarWithBound( 331 SCIP_PROOFSET* proofset, 332 SCIP_SET* set, 333 SCIP_VAR* var, 334 int pos, 335 SCIP_Bool* valid 336 ) 337 { 338 assert(proofset != NULL); 339 assert(var != NULL); 340 assert(pos >= 0 && pos < proofset->nnz); 341 assert(valid != NULL); 342 343 *valid = TRUE; 344 345 /* cancel with lower bound */ 346 if( proofset->vals[pos] > 0.0 ) 347 { 348 proofset->rhs -= proofset->vals[pos] * SCIPvarGetLbGlobal(var); 349 } 350 /* cancel with upper bound */ 351 else 352 { 353 assert(proofset->vals[pos] < 0.0); 354 proofset->rhs -= proofset->vals[pos] * SCIPvarGetUbGlobal(var); 355 } 356 357 --proofset->nnz; 358 359 proofset->vals[pos] = proofset->vals[proofset->nnz]; 360 proofset->inds[pos] = proofset->inds[proofset->nnz]; 361 proofset->vals[proofset->nnz] = 0.0; 362 proofset->inds[proofset->nnz] = 0; 363 364 if( SCIPsetIsInfinity(set, proofset->rhs) ) 365 *valid = FALSE; 366 } 367 368 /** creates and clears the proofset */ 369 SCIP_RETCODE SCIPconflictInitProofset( 370 SCIP_CONFLICT* conflict, /**< conflict analysis data */ 371 BMS_BLKMEM* blkmem /**< block memory of transformed problem */ 372 ) 373 { 374 assert(conflict != NULL); 375 assert(blkmem != NULL); 376 377 SCIP_CALL( proofsetCreate(&conflict->proofset, blkmem) ); 378 379 return SCIP_OKAY; 380 } 381 382 /** resizes proofsets array to be able to store at least num entries */ 383 static 384 SCIP_RETCODE conflictEnsureProofsetsMem( 385 SCIP_CONFLICT* conflict, /**< conflict analysis data */ 386 SCIP_SET* set, /**< global SCIP settings */ 387 int num /**< minimal number of slots in array */ 388 ) 389 { 390 assert(conflict != NULL); 391 assert(set != NULL); 392 393 if( num > conflict->proofsetssize ) 394 { 395 int newsize; 396 397 newsize = SCIPsetCalcMemGrowSize(set, num); 398 SCIP_ALLOC( BMSreallocMemoryArray(&conflict->proofsets, newsize) ); 399 conflict->proofsetssize = newsize; 400 } 401 assert(num <= conflict->proofsetssize); 402 403 return SCIP_OKAY; 404 } 405 406 /** add a proofset to the list of all proofsets */ 407 static 408 SCIP_RETCODE conflictInsertProofset( 409 SCIP_CONFLICT* conflict, /**< conflict analysis data */ 410 SCIP_SET* set, /**< global SCIP settings */ 411 SCIP_PROOFSET* proofset /**< proof set to add */ 412 ) 413 { 414 assert(conflict != NULL); 415 assert(proofset != NULL); 416 417 /* insert proofset into the sorted proofsets array */ 418 SCIP_CALL( conflictEnsureProofsetsMem(conflict, set, conflict->nproofsets + 1) ); 419 420 conflict->proofsets[conflict->nproofsets] = proofset; 421 ++conflict->nproofsets; 422 423 return SCIP_OKAY; 424 } 425 426 /** tighten the bound of a singleton variable in a constraint 427 * 428 * if the bound is contradicting with a global bound we cannot tighten the bound directly. 429 * in this case we need to create and add a constraint of size one such that propagating this constraint will 430 * enforce the infeasibility. 431 */ 432 static 433 SCIP_RETCODE tightenSingleVar( 434 SCIP_CONFLICT* conflict, /**< conflict analysis data */ 435 SCIP_SET* set, /**< global SCIP settings */ 436 SCIP_STAT* stat, /**< dynamic SCIP statistics */ 437 SCIP_TREE* tree, /**< tree data */ 438 BMS_BLKMEM* blkmem, /**< block memory */ 439 SCIP_PROB* origprob, /**< original problem */ 440 SCIP_PROB* transprob, /**< transformed problem */ 441 SCIP_REOPT* reopt, /**< reoptimization data */ 442 SCIP_LP* lp, /**< LP data */ 443 SCIP_BRANCHCAND* branchcand, /**< branching candidates */ 444 SCIP_EVENTQUEUE* eventqueue, /**< event queue */ 445 SCIP_CLIQUETABLE* cliquetable, /**< clique table */ 446 SCIP_VAR* var, /**< problem variable */ 447 SCIP_Real val, /**< coefficient of the variable */ 448 SCIP_Real rhs, /**< rhs of the constraint */ 449 SCIP_CONFTYPE prooftype, /**< type of the proof */ 450 int validdepth /**< depth where the bound change is valid */ 451 ) 452 { 453 SCIP_Real newbound; 454 SCIP_Bool applyglobal; 455 SCIP_BOUNDTYPE boundtype; 456 457 assert(tree != NULL); 458 assert(validdepth >= 0); 459 460 applyglobal = (validdepth <= SCIPtreeGetEffectiveRootDepth(tree)); 461 462 /* if variable and coefficient are integral the rhs can be rounded down */ 463 if( SCIPvarIsIntegral(var) && SCIPsetIsIntegral(set, val) ) 464 newbound = SCIPsetFeasFloor(set, rhs)/val; 465 else 466 newbound = rhs/val; 467 468 boundtype = (val > 0.0 ? SCIP_BOUNDTYPE_UPPER : SCIP_BOUNDTYPE_LOWER); 469 SCIPvarAdjustBd(var, set, boundtype, &newbound); 470 471 /* skip numerical unstable bound changes */ 472 if( applyglobal 473 && ((boundtype == SCIP_BOUNDTYPE_LOWER && SCIPsetIsLE(set, newbound, SCIPvarGetLbGlobal(var))) 474 || (boundtype == SCIP_BOUNDTYPE_UPPER && SCIPsetIsGE(set, newbound, SCIPvarGetUbGlobal(var)))) ) 475 { 476 return SCIP_OKAY; 477 } 478 479 /* the new bound contradicts a global bound, we can cutoff the root node immediately */ 480 if( applyglobal 481 && ((boundtype == SCIP_BOUNDTYPE_LOWER && SCIPsetIsGT(set, newbound, SCIPvarGetUbGlobal(var))) 482 || (boundtype == SCIP_BOUNDTYPE_UPPER && SCIPsetIsLT(set, newbound, SCIPvarGetLbGlobal(var)))) ) 483 { 484 SCIPsetDebugMsg(set, "detected global infeasibility at var <%s>: locdom=[%g,%g] glbdom=[%g,%g] new %s bound=%g\n", 485 SCIPvarGetName(var), SCIPvarGetLbLocal(var), 486 SCIPvarGetUbLocal(var), SCIPvarGetLbGlobal(var), SCIPvarGetUbGlobal(var), 487 (boundtype == SCIP_BOUNDTYPE_LOWER ? "lower" : "upper"), newbound); 488 SCIP_CALL( SCIPnodeCutoff(tree->path[0], set, stat, tree, transprob, origprob, reopt, lp, blkmem) ); 489 } 490 else 491 { 492 if( lp->strongbranching || !applyglobal ) 493 { 494 SCIP_CONS* cons; 495 SCIP_Real conslhs; 496 SCIP_Real consrhs; 497 char name[SCIP_MAXSTRLEN]; 498 499 SCIPsetDebugMsg(set, "add constraint <%s>[%c] %s %g to node #%lld in depth %d\n", 500 SCIPvarGetName(var), varGetChar(var), boundtype == SCIP_BOUNDTYPE_UPPER ? "<=" : ">=", newbound, 501 SCIPnodeGetNumber(tree->path[validdepth]), validdepth); 502 503 (void)SCIPsnprintf(name, SCIP_MAXSTRLEN, "pc_fix_%s", SCIPvarGetName(var)); 504 505 if( boundtype == SCIP_BOUNDTYPE_UPPER ) 506 { 507 conslhs = -SCIPsetInfinity(set); 508 consrhs = newbound; 509 } 510 else 511 { 512 conslhs = newbound; 513 consrhs = SCIPsetInfinity(set); 514 } 515 516 SCIP_CALL( SCIPcreateConsLinear(set->scip, &cons, name, 0, NULL, NULL, conslhs, consrhs, 517 FALSE, FALSE, FALSE, FALSE, TRUE, !applyglobal, FALSE, TRUE, TRUE, FALSE) ); 518 519 SCIP_CALL( SCIPaddCoefLinear(set->scip, cons, var, 1.0) ); 520 521 if( applyglobal ) 522 { 523 SCIP_CALL( SCIPprobAddCons(transprob, set, stat, cons) ); 524 } 525 else 526 { 527 SCIP_CALL( SCIPnodeAddCons(tree->path[validdepth], blkmem, set, stat, tree, cons) ); 528 } 529 530 SCIP_CALL( SCIPconsRelease(&cons, blkmem, set) ); 531 } 532 else 533 { 534 assert(applyglobal); 535 536 SCIPsetDebugMsg(set, "change global %s bound of <%s>[%c]: %g -> %g\n", 537 (boundtype == SCIP_BOUNDTYPE_LOWER ? "lower" : "upper"), 538 SCIPvarGetName(var), varGetChar(var), 539 (boundtype == SCIP_BOUNDTYPE_LOWER ? SCIPvarGetLbGlobal(var) : SCIPvarGetUbGlobal(var)), 540 newbound); 541 542 SCIP_CALL( SCIPnodeAddBoundchg(tree->path[0], blkmem, set, stat, transprob, origprob, tree, reopt, lp, branchcand, \ 543 eventqueue, cliquetable, var, newbound, boundtype, FALSE) ); 544 545 /* mark the node in the validdepth to be propagated again */ 546 SCIPnodePropagateAgain(tree->path[0], set, stat, tree); 547 } 548 } 549 550 if( applyglobal ) 551 ++conflict->nglbchgbds; 552 else 553 ++conflict->nlocchgbds; 554 555 if( prooftype == SCIP_CONFTYPE_INFEASLP || prooftype == SCIP_CONFTYPE_ALTINFPROOF ) 556 { 557 ++conflict->dualproofsinfnnonzeros; /* we count a global bound reduction as size 1 */ 558 ++conflict->ndualproofsinfsuccess; 559 ++conflict->ninflpsuccess; 560 561 if( applyglobal ) 562 ++conflict->ndualproofsinfglobal; 563 else 564 ++conflict->ndualproofsinflocal; 565 } 566 else 567 { 568 ++conflict->dualproofsbndnnonzeros; /* we count a global bound reduction as size 1 */ 569 ++conflict->ndualproofsbndsuccess; 570 ++conflict->nboundlpsuccess; 571 572 if( applyglobal ) 573 ++conflict->ndualproofsbndglobal; 574 else 575 ++conflict->ndualproofsbndlocal; 576 } 577 578 return SCIP_OKAY; 579 } 580 581 /** calculates the minimal activity of a given set of bounds and coefficients */ 582 static 583 SCIP_Real getMinActivity( 584 SCIP_SET* set, /**< global SCIP settings */ 585 SCIP_PROB* transprob, /**< transformed problem data */ 586 SCIP_Real* coefs, /**< coefficients in sparse representation */ 587 int* inds, /**< non-zero indices */ 588 int nnz, /**< number of non-zero indices */ 589 SCIP_Real* curvarlbs, /**< current lower bounds of active problem variables (or NULL for global bounds) */ 590 SCIP_Real* curvarubs /**< current upper bounds of active problem variables (or NULL for global bounds) */ 591 ) 592 { 593 SCIP_VAR** vars; 594 SCIP_Real QUAD(minact); 595 int i; 596 597 assert(coefs != NULL); 598 assert(inds != NULL); 599 600 vars = SCIPprobGetVars(transprob); 601 assert(vars != NULL); 602 603 QUAD_ASSIGN(minact, 0.0); 604 605 for( i = 0; i < nnz; i++ ) 606 { 607 SCIP_Real val; 608 SCIP_Real QUAD(delta); 609 int v = inds[i]; 610 611 assert(SCIPvarGetProbindex(vars[v]) == v); 612 613 val = coefs[i]; 614 615 if( val > 0.0 ) 616 { 617 SCIP_Real bnd; 618 619 assert(curvarlbs == NULL || !SCIPsetIsInfinity(set, -curvarlbs[v])); 620 621 bnd = (curvarlbs == NULL ? SCIPvarGetLbGlobal(vars[v]) : curvarlbs[v]); 622 623 if( SCIPsetIsInfinity(set, -bnd) ) 624 return -SCIPsetInfinity(set); 625 626 SCIPquadprecProdDD(delta, val, bnd); 627 } 628 else 629 { 630 SCIP_Real bnd; 631 632 assert(curvarubs == NULL || !SCIPsetIsInfinity(set, curvarubs[v])); 633 634 bnd = (curvarubs == NULL ? SCIPvarGetUbGlobal(vars[v]) : curvarubs[v]); 635 636 if( SCIPsetIsInfinity(set, bnd) ) 637 return -SCIPsetInfinity(set); 638 639 SCIPquadprecProdDD(delta, val, bnd); 640 } 641 642 /* update minimal activity */ 643 SCIPquadprecSumQQ(minact, minact, delta); 644 } 645 646 /* check whether the minmal activity is infinite */ 647 if( SCIPsetIsInfinity(set, QUAD_TO_DBL(minact)) ) 648 return SCIPsetInfinity(set); 649 if( SCIPsetIsInfinity(set, -QUAD_TO_DBL(minact)) ) 650 return -SCIPsetInfinity(set); 651 652 return QUAD_TO_DBL(minact); 653 } 654 655 /** calculates the minimal activity of a given set of bounds and coefficients */ 656 static 657 SCIP_Real getMaxActivity( 658 SCIP_SET* set, /**< global SCIP settings */ 659 SCIP_PROB* transprob, /**< transformed problem data */ 660 SCIP_Real* coefs, /**< coefficients in sparse representation */ 661 int* inds, /**< non-zero indices */ 662 int nnz, /**< number of non-zero indices */ 663 SCIP_Real* curvarlbs, /**< current lower bounds of active problem variables (or NULL for global bounds) */ 664 SCIP_Real* curvarubs /**< current upper bounds of active problem variables (or NULL for global bounds) */ 665 ) 666 { 667 SCIP_VAR** vars; 668 SCIP_Real QUAD(maxact); 669 int i; 670 671 assert(coefs != NULL); 672 assert(inds != NULL); 673 674 vars = SCIPprobGetVars(transprob); 675 assert(vars != NULL); 676 677 QUAD_ASSIGN(maxact, 0.0); 678 679 for( i = 0; i < nnz; i++ ) 680 { 681 SCIP_Real val; 682 SCIP_Real QUAD(delta); 683 int v = inds[i]; 684 685 assert(SCIPvarGetProbindex(vars[v]) == v); 686 687 val = coefs[i]; 688 689 if( val < 0.0 ) 690 { 691 SCIP_Real bnd; 692 693 assert(curvarlbs == NULL || !SCIPsetIsInfinity(set, -curvarlbs[v])); 694 695 bnd = (curvarlbs == NULL ? SCIPvarGetLbGlobal(vars[v]) : curvarlbs[v]); 696 697 if( SCIPsetIsInfinity(set, -bnd) ) 698 return SCIPsetInfinity(set); 699 700 SCIPquadprecProdDD(delta, val, bnd); 701 } 702 else 703 { 704 SCIP_Real bnd; 705 706 assert(curvarubs == NULL || !SCIPsetIsInfinity(set, curvarubs[v])); 707 708 bnd = (curvarubs == NULL ? SCIPvarGetUbGlobal(vars[v]) : curvarubs[v]); 709 710 if( SCIPsetIsInfinity(set, bnd) ) 711 return SCIPsetInfinity(set); 712 713 SCIPquadprecProdDD(delta, val, bnd); 714 } 715 716 /* update maximal activity */ 717 SCIPquadprecSumQQ(maxact, maxact, delta); 718 } 719 720 /* check whether the maximal activity got infinite */ 721 if( SCIPsetIsInfinity(set, QUAD_TO_DBL(maxact)) ) 722 return SCIPsetInfinity(set); 723 if( SCIPsetIsInfinity(set, -QUAD_TO_DBL(maxact)) ) 724 return -SCIPsetInfinity(set); 725 726 return QUAD_TO_DBL(maxact); 727 } 728 729 /** propagate a long proof */ 730 static 731 SCIP_RETCODE propagateLongProof( 732 SCIP_CONFLICT* conflict, /**< conflict analysis data */ 733 SCIP_SET* set, /**< global SCIP settings */ 734 SCIP_STAT* stat, /**< dynamic SCIP statistics */ 735 SCIP_REOPT* reopt, /**< reoptimization data */ 736 SCIP_TREE* tree, /**< tree data */ 737 BMS_BLKMEM* blkmem, /**< block memory */ 738 SCIP_PROB* origprob, /**< original problem */ 739 SCIP_PROB* transprob, /**< transformed problem */ 740 SCIP_LP* lp, /**< LP data */ 741 SCIP_BRANCHCAND* branchcand, /**< branching candidate storage */ 742 SCIP_EVENTQUEUE* eventqueue, /**< event queue */ 743 SCIP_CLIQUETABLE* cliquetable, /**< clique table data structure */ 744 SCIP_Real* coefs, /**< coefficients in sparse representation */ 745 int* inds, /**< non-zero indices */ 746 int nnz, /**< number of non-zero indices */ 747 SCIP_Real rhs, /**< right-hand side */ 748 SCIP_CONFTYPE conflicttype, /**< type of the conflict */ 749 int validdepth /**< depth where the proof is valid */ 750 ) 751 { 752 SCIP_VAR** vars; 753 SCIP_Real minact; 754 int i; 755 756 assert(coefs != NULL); 757 assert(inds != NULL); 758 assert(nnz >= 0); 759 760 vars = SCIPprobGetVars(transprob); 761 minact = getMinActivity(set, transprob, coefs, inds, nnz, NULL, NULL); 762 763 /* we cannot find global tightenings */ 764 if( SCIPsetIsInfinity(set, -minact) ) 765 return SCIP_OKAY; 766 767 for( i = 0; i < nnz; i++ ) 768 { 769 SCIP_VAR* var; 770 SCIP_Real val; 771 SCIP_Real resminact; 772 SCIP_Real lb; 773 SCIP_Real ub; 774 int pos; 775 776 pos = inds[i]; 777 val = coefs[i]; 778 var = vars[pos]; 779 lb = SCIPvarGetLbGlobal(var); 780 ub = SCIPvarGetUbGlobal(var); 781 782 assert(!SCIPsetIsZero(set, val)); 783 784 resminact = minact; 785 786 /* we got a potential new upper bound */ 787 if( val > 0.0 ) 788 { 789 SCIP_Real newub; 790 791 resminact -= (val * lb); 792 newub = (rhs - resminact)/val; 793 794 if( SCIPsetIsInfinity(set, newub) ) 795 continue; 796 797 /* we cannot tighten the upper bound */ 798 if( SCIPsetIsGE(set, newub, ub) ) 799 continue; 800 801 SCIP_CALL( tightenSingleVar(conflict, set, stat, tree, blkmem, origprob, transprob, reopt, lp, branchcand, \ 802 eventqueue, cliquetable, var, val, rhs-resminact, conflicttype, validdepth) ); 803 } 804 /* we got a potential new lower bound */ 805 else 806 { 807 SCIP_Real newlb; 808 809 resminact -= (val * ub); 810 newlb = (rhs - resminact)/val; 811 812 if( SCIPsetIsInfinity(set, -newlb) ) 813 continue; 814 815 /* we cannot tighten the lower bound */ 816 if( SCIPsetIsLE(set, newlb, lb) ) 817 continue; 818 819 SCIP_CALL( tightenSingleVar(conflict, set, stat, tree, blkmem, origprob, transprob, reopt, lp, branchcand, \ 820 eventqueue, cliquetable, var, val, rhs-resminact, conflicttype, validdepth) ); 821 } 822 823 /* the minimal activity should stay unchanged because we tightened the bound that doesn't contribute to the 824 * minimal activity 825 */ 826 assert(SCIPsetIsEQ(set, minact, getMinActivity(set, transprob, coefs, inds, nnz, NULL, NULL))); 827 } 828 829 return SCIP_OKAY; 830 } 831 832 /** creates a proof constraint and tries to add it to the storage */ 833 static 834 SCIP_RETCODE createAndAddProofcons( 835 SCIP_CONFLICT* conflict, /**< conflict analysis data */ 836 SCIP_CONFLICTSTORE* conflictstore, /**< conflict pool data */ 837 SCIP_PROOFSET* proofset, /**< proof set */ 838 SCIP_SET* set, /**< global SCIP settings */ 839 SCIP_STAT* stat, /**< dynamic SCIP statistics */ 840 SCIP_PROB* origprob, /**< original problem */ 841 SCIP_PROB* transprob, /**< transformed problem */ 842 SCIP_TREE* tree, /**< tree data */ 843 SCIP_REOPT* reopt, /**< reoptimization data */ 844 SCIP_LP* lp, /**< LP data */ 845 SCIP_BRANCHCAND* branchcand, /**< branching candidate storage */ 846 SCIP_EVENTQUEUE* eventqueue, /**< event queue */ 847 SCIP_CLIQUETABLE* cliquetable, /**< clique table data structure */ 848 BMS_BLKMEM* blkmem /**< block memory */ 849 ) 850 { 851 SCIP_CONS* cons; 852 SCIP_CONS* upgdcons; 853 SCIP_VAR** vars; 854 SCIP_Real* coefs; 855 int* inds; 856 SCIP_Real rhs; 857 SCIP_Real fillin; 858 SCIP_Real globalminactivity; 859 SCIP_Bool applyglobal; 860 SCIP_Bool toolong; 861 SCIP_Bool contonly; 862 SCIP_Bool hasrelaxvar; 863 SCIP_CONFTYPE conflicttype; 864 char name[SCIP_MAXSTRLEN]; 865 int nnz; 866 int i; 867 868 assert(conflict != NULL); 869 assert(conflictstore != NULL); 870 assert(proofset != NULL); 871 assert(proofset->validdepth == 0 || proofset->validdepth < SCIPtreeGetFocusDepth(tree)); 872 873 nnz = SCIPproofsetGetNVars(proofset); 874 875 if( nnz == 0 ) 876 return SCIP_OKAY; 877 878 vars = SCIPprobGetVars(transprob); 879 880 rhs = proofsetGetRhs(proofset); 881 assert(!SCIPsetIsInfinity(set, rhs)); 882 883 coefs = proofsetGetVals(proofset); 884 assert(coefs != NULL); 885 886 inds = proofsetGetInds(proofset); 887 assert(inds != NULL); 888 889 conflicttype = proofsetGetConftype(proofset); 890 891 applyglobal = (proofset->validdepth <= SCIPtreeGetEffectiveRootDepth(tree)); 892 893 if( applyglobal ) 894 { 895 SCIP_Real globalmaxactivity = getMaxActivity(set, transprob, coefs, inds, nnz, NULL, NULL); 896 897 /* check whether the alternative proof is redundant */ 898 if( SCIPsetIsLE(set, globalmaxactivity, rhs) ) 899 return SCIP_OKAY; 900 901 /* check whether the constraint proves global infeasibility */ 902 globalminactivity = getMinActivity(set, transprob, coefs, inds, nnz, NULL, NULL); 903 if( SCIPsetIsGT(set, globalminactivity, rhs) ) 904 { 905 SCIPsetDebugMsg(set, "detect global infeasibility: minactivity=%g, rhs=%g\n", globalminactivity, rhs); 906 907 SCIP_CALL( SCIPnodeCutoff(tree->path[proofset->validdepth], set, stat, tree, transprob, origprob, reopt, lp, blkmem) ); 908 909 goto UPDATESTATISTICS; 910 } 911 } 912 913 if( set->conf_minmaxvars >= nnz ) 914 toolong = FALSE; 915 else 916 { 917 SCIP_Real maxnnz; 918 919 if( transprob->startnconss < 100 ) 920 maxnnz = 0.85 * transprob->nvars; 921 else 922 maxnnz = (SCIP_Real)transprob->nvars; 923 924 fillin = nnz; 925 if( conflicttype == SCIP_CONFTYPE_INFEASLP || conflicttype == SCIP_CONFTYPE_ALTINFPROOF ) 926 { 927 fillin += SCIPconflictstoreGetNDualInfProofs(conflictstore) * SCIPconflictstoreGetAvgNnzDualInfProofs(conflictstore); 928 fillin /= (SCIPconflictstoreGetNDualInfProofs(conflictstore) + 1.0); 929 toolong = (fillin > MIN(2.0 * stat->avgnnz, maxnnz)); 930 } 931 else 932 { 933 assert(conflicttype == SCIP_CONFTYPE_BNDEXCEEDING || conflicttype == SCIP_CONFTYPE_ALTBNDPROOF); 934 935 fillin += SCIPconflictstoreGetNDualBndProofs(conflictstore) * SCIPconflictstoreGetAvgNnzDualBndProofs(conflictstore); 936 fillin /= (SCIPconflictstoreGetNDualBndProofs(conflictstore) + 1.0); 937 toolong = (fillin > MIN(1.5 * stat->avgnnz, maxnnz)); 938 } 939 940 toolong = (toolong && (nnz > set->conf_maxvarsfac * transprob->nvars)); 941 } 942 943 /* don't store global dual proofs that are too long / have too many non-zeros */ 944 if( toolong ) 945 { 946 if( applyglobal ) 947 { 948 SCIP_CALL( propagateLongProof(conflict, set, stat, reopt, tree, blkmem, origprob, transprob, lp, branchcand, 949 eventqueue, cliquetable, coefs, inds, nnz, rhs, conflicttype, proofset->validdepth) ); 950 } 951 return SCIP_OKAY; 952 } 953 954 /* check if conflict contains variables that are invalid after a restart to label it appropriately */ 955 hasrelaxvar = FALSE; 956 contonly = TRUE; 957 for( i = 0; i < nnz && (!hasrelaxvar || contonly); ++i ) 958 { 959 hasrelaxvar |= SCIPvarIsRelaxationOnly(vars[inds[i]]); 960 961 if( SCIPvarIsIntegral(vars[inds[i]]) ) 962 contonly = FALSE; 963 } 964 965 if( !applyglobal && contonly ) 966 return SCIP_OKAY; 967 968 if( conflicttype == SCIP_CONFTYPE_INFEASLP || conflicttype == SCIP_CONFTYPE_ALTINFPROOF ) 969 (void)SCIPsnprintf(name, SCIP_MAXSTRLEN, "dualproof_inf_%" SCIP_LONGINT_FORMAT, conflict->ndualproofsinfsuccess); 970 else if( conflicttype == SCIP_CONFTYPE_BNDEXCEEDING || conflicttype == SCIP_CONFTYPE_ALTBNDPROOF ) 971 (void)SCIPsnprintf(name, SCIP_MAXSTRLEN, "dualproof_bnd_%" SCIP_LONGINT_FORMAT, conflict->ndualproofsbndsuccess); 972 else 973 return SCIP_INVALIDCALL; 974 975 SCIP_CALL( SCIPcreateConsLinear(set->scip, &cons, name, 0, NULL, NULL, -SCIPsetInfinity(set), rhs, 976 FALSE, FALSE, FALSE, FALSE, TRUE, !applyglobal, 977 FALSE, TRUE, TRUE, FALSE) ); 978 979 for( i = 0; i < nnz; i++ ) 980 { 981 int v = inds[i]; 982 SCIP_CALL( SCIPaddCoefLinear(set->scip, cons, vars[v], coefs[i]) ); 983 } 984 985 /* do not upgrade linear constraints of size 1 */ 986 if( nnz > 1 ) 987 { 988 upgdcons = NULL; 989 /* try to automatically convert a linear constraint into a more specific and more specialized constraint */ 990 SCIP_CALL( SCIPupgradeConsLinear(set->scip, cons, &upgdcons) ); 991 if( upgdcons != NULL ) 992 { 993 SCIP_CALL( SCIPreleaseCons(set->scip, &cons) ); 994 cons = upgdcons; 995 996 if( conflicttype == SCIP_CONFTYPE_INFEASLP ) 997 conflicttype = SCIP_CONFTYPE_ALTINFPROOF; 998 else if( conflicttype == SCIP_CONFTYPE_BNDEXCEEDING ) 999 conflicttype = SCIP_CONFTYPE_ALTBNDPROOF; 1000 } 1001 } 1002 1003 /* mark constraint to be a conflict */ 1004 SCIPconsMarkConflict(cons); 1005 1006 /* add constraint to storage */ 1007 if( conflicttype == SCIP_CONFTYPE_INFEASLP || conflicttype == SCIP_CONFTYPE_ALTINFPROOF ) 1008 { 1009 /* add dual proof to storage */ 1010 SCIP_CALL( SCIPconflictstoreAddDualraycons(conflictstore, cons, blkmem, set, stat, transprob, reopt, hasrelaxvar) ); 1011 } 1012 else 1013 { 1014 SCIP_Real scale = 1.0; 1015 SCIP_Bool updateside = FALSE; 1016 1017 /* In some cases the constraint could not be updated to a more special type. However, it is possible that 1018 * constraint got scaled. Therefore, we need to be very careful when updating the lhs/rhs after the incumbent 1019 * solution has improved. 1020 */ 1021 if( conflicttype == SCIP_CONFTYPE_BNDEXCEEDING ) 1022 { 1023 SCIP_Real side; 1024 1025 #ifndef NDEBUG 1026 SCIP_CONSHDLR* conshdlr = SCIPconsGetHdlr(cons); 1027 1028 assert(conshdlr != NULL); 1029 assert(strcmp(SCIPconshdlrGetName(conshdlr), "linear") == 0); 1030 #endif 1031 side = SCIPgetLhsLinear(set->scip, cons); 1032 1033 if( !SCIPsetIsInfinity(set, -side) ) 1034 { 1035 if( SCIPsetIsZero(set, side) ) 1036 { 1037 scale = -1.0; 1038 } 1039 else 1040 { 1041 scale = proofsetGetRhs(proofset) / side; 1042 assert(SCIPsetIsNegative(set, scale)); 1043 } 1044 } 1045 else 1046 { 1047 side = SCIPgetRhsLinear(set->scip, cons); 1048 assert(!SCIPsetIsInfinity(set, side)); 1049 1050 if( SCIPsetIsZero(set, side) ) 1051 { 1052 scale = 1.0; 1053 } 1054 else 1055 { 1056 scale = proofsetGetRhs(proofset) / side; 1057 assert(SCIPsetIsPositive(set, scale)); 1058 } 1059 } 1060 updateside = TRUE; 1061 } 1062 1063 /* add dual proof to storage */ 1064 SCIP_CALL( SCIPconflictstoreAddDualsolcons(conflictstore, cons, blkmem, set, stat, transprob, reopt, scale, updateside, hasrelaxvar) ); 1065 } 1066 1067 if( applyglobal ) /*lint !e774*/ 1068 { 1069 /* add the constraint to the global problem */ 1070 SCIP_CALL( SCIPprobAddCons(transprob, set, stat, cons) ); 1071 } 1072 else 1073 { 1074 SCIP_CALL( SCIPnodeAddCons(tree->path[proofset->validdepth], blkmem, set, stat, tree, cons) ); 1075 } 1076 1077 SCIPsetDebugMsg(set, "added proof-constraint to node %p (#%lld) in depth %d (nproofconss %d)\n", 1078 (void*)tree->path[proofset->validdepth], SCIPnodeGetNumber(tree->path[proofset->validdepth]), 1079 proofset->validdepth, 1080 (conflicttype == SCIP_CONFTYPE_INFEASLP || conflicttype == SCIP_CONFTYPE_ALTINFPROOF) 1081 ? SCIPconflictstoreGetNDualInfProofs(conflictstore) : SCIPconflictstoreGetNDualBndProofs(conflictstore)); 1082 1083 /* release the constraint */ 1084 SCIP_CALL( SCIPreleaseCons(set->scip, &cons) ); 1085 1086 UPDATESTATISTICS: 1087 /* update statistics */ 1088 if( conflicttype == SCIP_CONFTYPE_INFEASLP || conflicttype == SCIP_CONFTYPE_ALTINFPROOF ) 1089 { 1090 conflict->dualproofsinfnnonzeros += nnz; 1091 if( applyglobal ) /*lint !e774*/ 1092 ++conflict->ndualproofsinfglobal; 1093 else 1094 ++conflict->ndualproofsinflocal; 1095 ++conflict->ndualproofsinfsuccess; 1096 } 1097 else 1098 { 1099 assert(conflicttype == SCIP_CONFTYPE_BNDEXCEEDING || conflicttype == SCIP_CONFTYPE_ALTBNDPROOF); 1100 conflict->dualproofsbndnnonzeros += nnz; 1101 if( applyglobal ) /*lint !e774*/ 1102 ++conflict->ndualproofsbndglobal; 1103 else 1104 ++conflict->ndualproofsbndlocal; 1105 1106 ++conflict->ndualproofsbndsuccess; 1107 } 1108 return SCIP_OKAY; 1109 } 1110 1111 /** create proof constraints out of proof sets */ 1112 SCIP_RETCODE SCIPconflictFlushProofset( 1113 SCIP_CONFLICT* conflict, /**< conflict analysis data */ 1114 SCIP_CONFLICTSTORE* conflictstore, /**< conflict store */ 1115 BMS_BLKMEM* blkmem, /**< block memory */ 1116 SCIP_SET* set, /**< global SCIP settings */ 1117 SCIP_STAT* stat, /**< dynamic problem statistics */ 1118 SCIP_PROB* transprob, /**< transformed problem after presolve */ 1119 SCIP_PROB* origprob, /**< original problem */ 1120 SCIP_TREE* tree, /**< branch and bound tree */ 1121 SCIP_REOPT* reopt, /**< reoptimization data structure */ 1122 SCIP_LP* lp, /**< current LP data */ 1123 SCIP_BRANCHCAND* branchcand, /**< branching candidate storage */ 1124 SCIP_EVENTQUEUE* eventqueue, /**< event queue */ 1125 SCIP_CLIQUETABLE* cliquetable /**< clique table data structure */ 1126 ) 1127 { 1128 assert(conflict != NULL); 1129 1130 if( proofsetGetConftype(conflict->proofset) != SCIP_CONFTYPE_UNKNOWN ) 1131 { 1132 /* only one variable has a coefficient different to zero, we add this bound change instead of a constraint */ 1133 if( SCIPproofsetGetNVars(conflict->proofset) == 1 ) 1134 { 1135 SCIP_VAR** vars; 1136 SCIP_Real* coefs; 1137 int* inds; 1138 SCIP_Real rhs; 1139 1140 vars = SCIPprobGetVars(transprob); 1141 1142 coefs = proofsetGetVals(conflict->proofset); 1143 inds = proofsetGetInds(conflict->proofset); 1144 rhs = proofsetGetRhs(conflict->proofset); 1145 1146 SCIP_CALL( tightenSingleVar(conflict, set, stat, tree, blkmem, origprob, transprob, reopt, lp, \ 1147 branchcand, eventqueue, cliquetable, vars[inds[0]], coefs[0], rhs, conflict->proofset->conflicttype, 1148 conflict->proofset->validdepth) ); 1149 } 1150 else 1151 { 1152 SCIP_Bool skipinitialproof = FALSE; 1153 1154 /* prefer an infeasibility proof 1155 * 1156 * todo: check whether this is really what we want 1157 */ 1158 if( set->conf_prefinfproof && proofsetGetConftype(conflict->proofset) == SCIP_CONFTYPE_BNDEXCEEDING ) 1159 { 1160 int i; 1161 1162 for( i = 0; i < conflict->nproofsets; i++ ) 1163 { 1164 if( proofsetGetConftype(conflict->proofsets[i]) == SCIP_CONFTYPE_INFEASLP ) 1165 { 1166 skipinitialproof = TRUE; 1167 break; 1168 } 1169 } 1170 } 1171 1172 if( !skipinitialproof ) 1173 { 1174 /* create and add the original proof */ 1175 SCIP_CALL( createAndAddProofcons(conflict, conflictstore, conflict->proofset, set, stat, origprob, transprob, \ 1176 tree, reopt, lp, branchcand, eventqueue, cliquetable, blkmem) ); 1177 } 1178 } 1179 1180 /* clear the proof set anyway */ 1181 proofsetClear(conflict->proofset); 1182 } 1183 1184 if( conflict->nproofsets > 0 ) 1185 { 1186 int i; 1187 1188 for( i = 0; i < conflict->nproofsets; i++ ) 1189 { 1190 assert(conflict->proofsets[i] != NULL); 1191 assert(proofsetGetConftype(conflict->proofsets[i]) != SCIP_CONFTYPE_UNKNOWN); 1192 1193 /* only one variable has a coefficient different to zero, we add this bound change instead of a constraint */ 1194 if( SCIPproofsetGetNVars(conflict->proofsets[i]) == 1 ) 1195 { 1196 SCIP_VAR** vars; 1197 SCIP_Real* coefs; 1198 int* inds; 1199 SCIP_Real rhs; 1200 1201 vars = SCIPprobGetVars(transprob); 1202 1203 coefs = proofsetGetVals(conflict->proofsets[i]); 1204 inds = proofsetGetInds(conflict->proofsets[i]); 1205 rhs = proofsetGetRhs(conflict->proofsets[i]); 1206 1207 SCIP_CALL( tightenSingleVar(conflict, set, stat, tree, blkmem, origprob, transprob, reopt, lp, 1208 branchcand, eventqueue, cliquetable, vars[inds[0]], coefs[0], rhs, 1209 conflict->proofsets[i]->conflicttype, conflict->proofsets[i]->validdepth) ); 1210 } 1211 else 1212 { 1213 /* create and add proof constraint */ 1214 SCIP_CALL( createAndAddProofcons(conflict, conflictstore, conflict->proofsets[i], set, stat, origprob, \ 1215 transprob, tree, reopt, lp, branchcand, eventqueue, cliquetable, blkmem) ); 1216 } 1217 } 1218 1219 /* free all proofsets */ 1220 for( i = 0; i < conflict->nproofsets; i++ ) 1221 SCIPproofsetFree(&conflict->proofsets[i], blkmem); 1222 1223 conflict->nproofsets = 0; 1224 } 1225 1226 return SCIP_OKAY; 1227 } 1228 1229 1230 1231 #ifdef SCIP_DEBUG 1232 /** print violation for debugging */ 1233 static 1234 void debugPrintViolationInfo( 1235 SCIP_SET* set, /**< global SCIP settings */ 1236 SCIP_Real minact, /**< min activity */ 1237 SCIP_Real rhs, /**< right hand side */ 1238 const char* infostr /**< additional info for this debug message, or NULL */ 1239 ) 1240 { 1241 SCIPsetDebugMsg(set, "-> %sminact=%.15g rhs=%.15g violation=%.15g\n",infostr != NULL ? infostr : "" , minact, rhs, minact - rhs); 1242 } 1243 #else 1244 #define debugPrintViolationInfo(...) /**/ 1245 #endif 1246 1247 /** apply coefficient tightening */ 1248 static 1249 void tightenCoefficients( 1250 SCIP_SET* set, /**< global SCIP settings */ 1251 SCIP_PROOFSET* proofset, /**< proof set */ 1252 int* nchgcoefs, /**< pointer to store number of changed coefficients */ 1253 SCIP_Bool* redundant /**< pointer to store whether the proof set is redundant */ 1254 ) 1255 { 1256 #ifdef SCIP_DEBUG 1257 SCIP_Real absmax = 0.0; 1258 SCIP_Real absmin = SCIPsetInfinity(set); 1259 int i; 1260 1261 for( i = 0; i < proofset->nnz; i++ ) 1262 { 1263 absmax = MAX(absmax, REALABS(proofset->vals[i])); 1264 absmin = MIN(absmin, REALABS(proofset->vals[i])); 1265 } 1266 #endif 1267 1268 (*redundant) = SCIPcutsTightenCoefficients(set->scip, FALSE, proofset->vals, &proofset->rhs, proofset->inds, &proofset->nnz, nchgcoefs); 1269 1270 #ifdef SCIP_DEBUG 1271 { 1272 SCIP_Real newabsmax = 0.0; 1273 SCIP_Real newabsmin = SCIPsetInfinity(set); 1274 1275 for( i = 0; i < proofset->nnz; i++ ) 1276 { 1277 newabsmax = MAX(newabsmax, REALABS(proofset->vals[i])); 1278 newabsmin = MIN(newabsmin, REALABS(proofset->vals[i])); 1279 } 1280 1281 SCIPsetDebugMsg(set, "coefficient tightening: [%.15g,%.15g] -> [%.15g,%.15g] (nnz: %d, nchg: %d rhs: %.15g)\n", 1282 absmin, absmax, newabsmin, newabsmax, proofsetGetNVars(proofset), *nchgcoefs, proofsetGetRhs(proofset)); 1283 printf("coefficient tightening: [%.15g,%.15g] -> [%.15g,%.15g] (nnz: %d, nchg: %d rhs: %.15g)\n", 1284 absmin, absmax, newabsmin, newabsmax, proofsetGetNVars(proofset), *nchgcoefs, proofsetGetRhs(proofset)); 1285 } 1286 #endif 1287 } 1288 1289 /** try to generate alternative proofs by applying subadditive functions */ 1290 static 1291 SCIP_RETCODE separateAlternativeProofs( 1292 SCIP_CONFLICT* conflict, /**< conflict analysis data */ 1293 SCIP_SET* set, /**< global SCIP settings */ 1294 SCIP_STAT* stat, /**< dynamic SCIP statistics */ 1295 SCIP_PROB* transprob, /**< transformed problem */ 1296 SCIP_TREE* tree, /**< tree data */ 1297 BMS_BLKMEM* blkmem, /**< block memory */ 1298 SCIP_AGGRROW* proofrow, /**< proof rows data */ 1299 SCIP_Real* curvarlbs, /**< current lower bounds of active problem variables */ 1300 SCIP_Real* curvarubs, /**< current upper bounds of active problem variables */ 1301 SCIP_CONFTYPE conflicttype /**< type of the conflict */ 1302 ) 1303 { 1304 SCIP_VAR** vars; 1305 SCIP_SOL* refsol; 1306 SCIP_Real* cutcoefs; 1307 SCIP_Real cutefficacy; 1308 SCIP_Real cutrhs; 1309 SCIP_Real proofefficacy; 1310 SCIP_Real efficacynorm; 1311 SCIP_Bool islocal; 1312 SCIP_Bool cutsuccess; 1313 SCIP_Bool success; 1314 SCIP_Bool infdelta; 1315 int* cutinds; 1316 int* inds; 1317 int cutnnz; 1318 int nnz; 1319 int nvars; 1320 int i; 1321 1322 vars = SCIPprobGetVars(transprob); 1323 nvars = SCIPprobGetNVars(transprob); 1324 1325 inds = SCIPaggrRowGetInds(proofrow); 1326 nnz = SCIPaggrRowGetNNz(proofrow); 1327 1328 proofefficacy = SCIPaggrRowGetMinActivity(set, transprob, proofrow, curvarlbs, curvarubs, &infdelta); 1329 1330 if( infdelta ) 1331 return SCIP_OKAY; 1332 1333 proofefficacy -= SCIPaggrRowGetRhs(proofrow); 1334 1335 efficacynorm = SCIPaggrRowCalcEfficacyNorm(set->scip, proofrow); 1336 proofefficacy /= MAX(1e-6, efficacynorm); 1337 1338 /* create reference solution */ 1339 SCIP_CALL( SCIPcreateSol(set->scip, &refsol, NULL) ); 1340 1341 /* initialize with average solution */ 1342 for( i = 0; i < nvars; i++ ) 1343 { 1344 SCIP_CALL( SCIPsolSetVal(refsol, set, stat, tree, vars[i], SCIPvarGetAvgSol(vars[i])) ); 1345 } 1346 1347 /* set all variables that are part of the proof to its active local bound */ 1348 for( i = 0; i < nnz; i++ ) 1349 { 1350 SCIP_Real val = SCIPaggrRowGetProbvarValue(proofrow, inds[i]); 1351 1352 if( val > 0.0 ) 1353 { 1354 SCIP_CALL( SCIPsolSetVal(refsol, set, stat, tree, vars[inds[i]], curvarubs[inds[i]]) ); 1355 } 1356 else 1357 { 1358 SCIP_CALL( SCIPsolSetVal(refsol, set, stat, tree, vars[inds[i]], curvarlbs[inds[i]]) ); 1359 } 1360 } 1361 1362 SCIP_CALL( SCIPsetAllocBufferArray(set, &cutcoefs, nvars) ); 1363 SCIP_CALL( SCIPsetAllocBufferArray(set, &cutinds, nvars) ); 1364 1365 cutnnz = 0; 1366 cutefficacy = -SCIPsetInfinity(set); 1367 1368 /* apply flow cover */ 1369 SCIP_CALL( SCIPcalcFlowCover(set->scip, refsol, POSTPROCESS, BOUNDSWITCH, ALLOWLOCAL, proofrow, \ 1370 cutcoefs, &cutrhs, cutinds, &cutnnz, &cutefficacy, NULL, &islocal, &cutsuccess) ); 1371 success = cutsuccess; 1372 1373 /* apply MIR */ 1374 SCIP_CALL( SCIPcutGenerationHeuristicCMIR(set->scip, refsol, POSTPROCESS, BOUNDSWITCH, USEVBDS, ALLOWLOCAL, INT_MAX, \ 1375 NULL, NULL, MINFRAC, MAXFRAC, proofrow, cutcoefs, &cutrhs, cutinds, &cutnnz, &cutefficacy, NULL, \ 1376 &islocal, &cutsuccess) ); 1377 success = (success || cutsuccess); 1378 1379 /* replace the current proof */ 1380 if( success && !islocal && SCIPsetIsPositive(set, cutefficacy) && cutefficacy * nnz > proofefficacy * cutnnz ) 1381 { 1382 SCIP_PROOFSET* alternativeproofset; 1383 SCIP_Bool redundant; 1384 int nchgcoefs; 1385 1386 SCIP_CALL( proofsetCreate(&alternativeproofset, blkmem) ); 1387 alternativeproofset->conflicttype = (conflicttype == SCIP_CONFTYPE_INFEASLP ? SCIP_CONFTYPE_ALTINFPROOF : SCIP_CONFTYPE_ALTBNDPROOF); 1388 1389 SCIP_CALL( proofsetAddSparseData(alternativeproofset, blkmem, cutcoefs, cutinds, cutnnz, cutrhs) ); 1390 1391 /* apply coefficient tightening */ 1392 tightenCoefficients(set, alternativeproofset, &nchgcoefs, &redundant); 1393 1394 if( !redundant ) 1395 { 1396 SCIP_CALL( conflictInsertProofset(conflict, set, alternativeproofset) ); 1397 } 1398 else 1399 { 1400 SCIPproofsetFree(&alternativeproofset, blkmem); 1401 } 1402 } /*lint !e438*/ 1403 1404 SCIPsetFreeBufferArray(set, &cutinds); 1405 SCIPsetFreeBufferArray(set, &cutcoefs); 1406 1407 SCIP_CALL( SCIPfreeSol(set->scip, &refsol) ); 1408 1409 return SCIP_OKAY; 1410 } 1411 1412 /** tighten a given infeasibility proof a^Tx <= b with minact > b w.r.t. local bounds 1413 * 1414 * 1) Apply cut generating functions 1415 * - c-MIR 1416 * - Flow-cover 1417 * - TODO: implement other subadditive functions 1418 * 2) Remove continuous variables contributing with its global bound 1419 * - TODO: implement a variant of non-zero-cancellation 1420 */ 1421 static 1422 SCIP_RETCODE tightenDualproof( 1423 SCIP_CONFLICT* conflict, /**< conflict analysis data */ 1424 SCIP_SET* set, /**< global SCIP settings */ 1425 SCIP_STAT* stat, /**< dynamic SCIP statistics */ 1426 BMS_BLKMEM* blkmem, /**< block memory */ 1427 SCIP_PROB* transprob, /**< transformed problem */ 1428 SCIP_TREE* tree, /**< tree data */ 1429 SCIP_AGGRROW* proofrow, /**< aggregated row representing the proof */ 1430 int validdepth, /**< depth where the proof is valid */ 1431 SCIP_Real* curvarlbs, /**< current lower bounds of active problem variables */ 1432 SCIP_Real* curvarubs, /**< current upper bounds of active problem variables */ 1433 SCIP_Bool initialproof /**< do we analyze the initial reason of infeasibility? */ 1434 ) 1435 { 1436 SCIP_VAR** vars; 1437 SCIP_Real* vals; 1438 int* inds; 1439 SCIP_PROOFSET* proofset; 1440 SCIP_Bool valid; 1441 SCIP_Bool redundant; 1442 int nnz; 1443 int nchgcoefs; 1444 int nbinvars; 1445 int ncontvars; 1446 int nintvars; 1447 int i; 1448 1449 assert(conflict->proofset != NULL); 1450 assert(curvarlbs != NULL); 1451 assert(curvarubs != NULL); 1452 1453 vars = SCIPprobGetVars(transprob); 1454 nbinvars = 0; 1455 nintvars = 0; 1456 ncontvars = 0; 1457 1458 inds = SCIPaggrRowGetInds(proofrow); 1459 nnz = SCIPaggrRowGetNNz(proofrow); 1460 1461 /* count number of binary, integer, and continuous variables */ 1462 for( i = 0; i < nnz; i++ ) 1463 { 1464 assert(SCIPvarGetProbindex(vars[inds[i]]) == inds[i]); 1465 1466 if( SCIPvarIsBinary(vars[inds[i]]) ) 1467 ++nbinvars; 1468 else if( SCIPvarIsIntegral(vars[inds[i]]) ) 1469 ++nintvars; 1470 else 1471 ++ncontvars; 1472 } 1473 1474 SCIPsetDebugMsg(set, "start dual proof tightening:\n"); 1475 SCIPsetDebugMsg(set, "-> tighten dual proof: nvars=%d (bin=%d, int=%d, cont=%d)\n", 1476 nnz, nbinvars, nintvars, ncontvars); 1477 debugPrintViolationInfo(set, aggrRowGetMinActivity(set, transprob, proofrow, curvarlbs, curvarubs, NULL), SCIPaggrRowGetRhs(proofrow), NULL); 1478 1479 /* try to find an alternative proof of local infeasibility that is stronger */ 1480 if( set->conf_sepaaltproofs ) 1481 { 1482 SCIP_CALL( separateAlternativeProofs(conflict, set, stat, transprob, tree, blkmem, proofrow, curvarlbs, curvarubs, 1483 conflict->conflictset->conflicttype) ); 1484 } 1485 1486 if( initialproof ) 1487 proofset = conflict->proofset; 1488 else 1489 { 1490 SCIP_CALL( proofsetCreate(&proofset, blkmem) ); 1491 } 1492 1493 /* start with a proofset containing all variables with a non-zero coefficient in the dual proof */ 1494 SCIP_CALL( proofsetAddAggrrow(proofset, set, blkmem, proofrow) ); 1495 proofset->conflicttype = conflict->conflictset->conflicttype; 1496 proofset->validdepth = validdepth; 1497 1498 /* get proof data */ 1499 vals = proofsetGetVals(proofset); 1500 inds = proofsetGetInds(proofset); 1501 nnz = SCIPproofsetGetNVars(proofset); 1502 1503 #ifndef NDEBUG 1504 for( i = 0; i < nnz; i++ ) 1505 { 1506 int idx = inds[i]; 1507 if( vals[i] > 0.0 ) 1508 assert(!SCIPsetIsInfinity(set, -curvarlbs[idx])); 1509 if( vals[i] < 0.0 ) 1510 assert(!SCIPsetIsInfinity(set, curvarubs[idx])); 1511 } 1512 #endif 1513 1514 /* remove continuous variable contributing with their global bound 1515 * 1516 * todo: check whether we also want to do that for bound exceeding proofs, but then we cannot update the 1517 * conflict anymore 1518 */ 1519 if( proofset->conflicttype == SCIP_CONFTYPE_INFEASLP ) 1520 { 1521 /* remove all continuous variables that have equal global and local bounds (ub or lb depend on the sign) 1522 * from the proof 1523 */ 1524 1525 for( i = 0; i < nnz && nnz > 1; ) 1526 { 1527 SCIP_Real val; 1528 int idx = inds[i]; 1529 1530 assert(vars[idx] != NULL); 1531 1532 val = vals[i]; 1533 assert(!SCIPsetIsZero(set, val)); 1534 1535 /* skip integral variables */ 1536 if( SCIPvarGetType(vars[idx]) != SCIP_VARTYPE_CONTINUOUS && SCIPvarGetType(vars[idx]) != SCIP_VARTYPE_IMPLINT ) 1537 { 1538 i++; 1539 continue; 1540 } 1541 else 1542 { 1543 SCIP_Real glbbd; 1544 SCIP_Real locbd; 1545 1546 /* get appropriate global and local bounds */ 1547 glbbd = (val < 0.0 ? SCIPvarGetUbGlobal(vars[idx]) : SCIPvarGetLbGlobal(vars[idx])); 1548 locbd = (val < 0.0 ? curvarubs[idx] : curvarlbs[idx]); 1549 1550 if( !SCIPsetIsEQ(set, glbbd, locbd) ) 1551 { 1552 i++; 1553 continue; 1554 } 1555 1556 SCIPsetDebugMsg(set, "-> remove continuous variable <%s>: glb=[%g,%g], loc=[%g,%g], val=%g\n", 1557 SCIPvarGetName(vars[idx]), SCIPvarGetLbGlobal(vars[idx]), SCIPvarGetUbGlobal(vars[idx]), 1558 curvarlbs[idx], curvarubs[idx], val); 1559 1560 proofsetCancelVarWithBound(proofset, set, vars[idx], i, &valid); 1561 assert(valid); /* this should be always fulfilled at this place */ 1562 1563 --nnz; 1564 } 1565 } 1566 } 1567 1568 /* apply coefficient tightening to initial proof */ 1569 tightenCoefficients(set, proofset, &nchgcoefs, &redundant); 1570 1571 /* it can happen that the constraints is almost globally redundant w.r.t to the maximal activity, 1572 * e.g., due to numerics. in this case, we want to discard the proof 1573 */ 1574 if( redundant ) 1575 { 1576 #ifndef NDEBUG 1577 SCIP_Real eps = MIN(0.01, 10.0*set->num_feastol); 1578 assert(proofset->rhs - getMaxActivity(set, transprob, proofset->vals, proofset->inds, proofset->nnz, NULL, NULL) < eps); 1579 #endif 1580 if( initialproof ) 1581 { 1582 proofsetClear(proofset); 1583 } 1584 else 1585 { 1586 SCIPproofsetFree(&proofset, blkmem); 1587 } 1588 } 1589 else 1590 { 1591 if( !initialproof ) 1592 { 1593 SCIP_CALL( conflictInsertProofset(conflict, set, proofset) ); 1594 } 1595 1596 if( nchgcoefs > 0 ) 1597 { 1598 if( proofset->conflicttype == SCIP_CONFTYPE_INFEASLP ) 1599 proofset->conflicttype = SCIP_CONFTYPE_ALTINFPROOF; 1600 else if( proofset->conflicttype == SCIP_CONFTYPE_BNDEXCEEDING ) 1601 proofset->conflicttype = SCIP_CONFTYPE_ALTBNDPROOF; 1602 } 1603 } 1604 1605 return SCIP_OKAY; 1606 } 1607 1608 /** perform conflict analysis based on a dual unbounded ray 1609 * 1610 * given an aggregation of rows lhs <= a^Tx such that lhs > maxactivity. if the constraint has size one we add a 1611 * bound change instead of the constraint. 1612 */ 1613 SCIP_RETCODE SCIPconflictAnalyzeDualProof( 1614 SCIP_CONFLICT* conflict, /**< conflict analysis data */ 1615 SCIP_SET* set, /**< global SCIP settings */ 1616 SCIP_STAT* stat, /**< dynamic SCIP statistics */ 1617 BMS_BLKMEM* blkmem, /**< block memory */ 1618 SCIP_PROB* origprob, /**< original problem */ 1619 SCIP_PROB* transprob, /**< transformed problem */ 1620 SCIP_TREE* tree, /**< tree data */ 1621 SCIP_REOPT* reopt, /**< reoptimization data */ 1622 SCIP_LP* lp, /**< LP data */ 1623 SCIP_AGGRROW* proofrow, /**< aggregated row representing the proof */ 1624 int validdepth, /**< valid depth of the dual proof */ 1625 SCIP_Real* curvarlbs, /**< current lower bounds of active problem variables */ 1626 SCIP_Real* curvarubs, /**< current upper bounds of active problem variables */ 1627 SCIP_Bool initialproof, /**< do we analyze the initial reason of infeasibility? */ 1628 SCIP_Bool* globalinfeasible, /**< pointer to store whether global infeasibility could be proven */ 1629 SCIP_Bool* success /**< pointer to store success result */ 1630 ) 1631 { 1632 SCIP_Real rhs; 1633 SCIP_Real minact; 1634 SCIP_Bool infdelta; 1635 int nnz; 1636 1637 assert(set != NULL); 1638 assert(transprob != NULL); 1639 assert(validdepth >= 0); 1640 assert(validdepth == 0 || validdepth < SCIPtreeGetFocusDepth(tree)); 1641 1642 /* get sparse data */ 1643 nnz = SCIPaggrRowGetNNz(proofrow); 1644 rhs = SCIPaggrRowGetRhs(proofrow); 1645 1646 *globalinfeasible = FALSE; 1647 *success = FALSE; 1648 1649 /* get minimal activity w.r.t. local bounds */ 1650 minact = SCIPaggrRowGetMinActivity(set, transprob, proofrow, curvarlbs, curvarubs, &infdelta); 1651 1652 if( infdelta ) 1653 return SCIP_OKAY; 1654 1655 /* only run is the proof proves local infeasibility */ 1656 if( SCIPsetIsFeasLE(set, minact, rhs) ) 1657 return SCIP_OKAY; 1658 1659 /* if the farkas-proof is empty, the node and its sub tree can be cut off completely */ 1660 if( nnz == 0 ) 1661 { 1662 SCIPsetDebugMsg(set, " -> empty farkas-proof in depth %d cuts off sub tree at depth %d\n", SCIPtreeGetFocusDepth(tree), validdepth); 1663 1664 SCIP_CALL( SCIPnodeCutoff(tree->path[validdepth], set, stat, tree, transprob, origprob, reopt, lp, blkmem) ); 1665 1666 *globalinfeasible = TRUE; 1667 *success = TRUE; 1668 1669 ++conflict->ndualproofsinfsuccess; 1670 1671 return SCIP_OKAY; 1672 } 1673 1674 /* try to enforce the constraint based on a dual ray */ 1675 SCIP_CALL( tightenDualproof(conflict, set, stat, blkmem, transprob, tree, proofrow, validdepth, 1676 curvarlbs, curvarubs, initialproof) ); 1677 1678 if( *globalinfeasible ) 1679 { 1680 SCIPsetDebugMsg(set, "detect global: cutoff root node\n"); 1681 SCIP_CALL( SCIPnodeCutoff(tree->path[0], set, stat, tree, transprob, origprob, reopt, lp, blkmem) ); 1682 *success = TRUE; 1683 1684 ++conflict->ndualproofsinfsuccess; 1685 } 1686 1687 return SCIP_OKAY; 1688 } 1689