1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 2 /* */ 3 /* This file is part of the class library */ 4 /* SoPlex --- the Sequential object-oriented simPlex. */ 5 /* */ 6 /* Copyright (c) 1996-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 SoPlex; see the file LICENSE. If not email to soplex@zib.de. */ 22 /* */ 23 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 24 25 26 #ifndef SOPLEX_WITH_PAPILO 27 28 namespace soplex 29 { 30 31 template <class R> class Presol : public SPxSimplifier<R> 32 { 33 34 public: 35 36 //------------------------------------ 37 ///@name Constructors / destructors 38 ///@{ 39 /// default constructor. 40 explicit Presol(Timer::TYPE ttype = Timer::USER_TIME) : SPxSimplifier<R>("PaPILO", ttype) 41 { ; }; 42 43 /// copy constructor. 44 Presol(const Presol& old) : SPxSimplifier<R>(old) { ; } 45 46 /// assignment operator 47 Presol& operator=(const Presol& rhs) 48 { 49 return *this; 50 } 51 52 /// destructor. 53 virtual ~Presol() 54 { 55 ; 56 } 57 58 SPxSimplifier<R>* clone() const 59 { 60 return new Presol(*this); 61 } 62 63 virtual typename SPxSimplifier<R>::Result simplify(SPxLPBase<R>& lp, 64 Real remainingTime, bool keepbounds, uint32_t seed) 65 { 66 assert(false); 67 return SPxSimplifier<R>::OKAY; 68 }; 69 70 virtual void unsimplify(const VectorBase<R>&, const VectorBase<R>&, 71 const VectorBase<R>&, const VectorBase<R>&, 72 const typename SPxSolverBase<R>::VarStatus[], 73 const typename SPxSolverBase<R>::VarStatus[], 74 bool isOptimal) 75 { 76 assert(false); 77 }; 78 79 /// returns result status of the simplification 80 virtual typename SPxSimplifier<R>::Result result() const 81 { 82 assert(false); 83 return SPxSimplifier<R>::OKAY; 84 } 85 86 /// specifies whether an optimal solution has already been unsimplified. 87 virtual bool isUnsimplified() const 88 { 89 assert(false); 90 return false; 91 } 92 93 /// returns a reference to the unsimplified primal solution. 94 virtual const VectorBase<R>& unsimplifiedPrimal() 95 { 96 assert(false); 97 static const VectorBase<R>& emptyVector = VectorBase<R>(); 98 return emptyVector; 99 } 100 101 /// returns a reference to the unsimplified dual solution. 102 virtual const VectorBase<R>& unsimplifiedDual() 103 { 104 assert(false); 105 static const VectorBase<R>& emptyVector = VectorBase<R>(); 106 return emptyVector; 107 } 108 109 /// returns a reference to the unsimplified slack values. 110 virtual const VectorBase<R>& unsimplifiedSlacks() 111 { 112 assert(false); 113 static const VectorBase<R>& emptyVector = VectorBase<R>(); 114 return emptyVector; 115 } 116 117 /// returns a reference to the unsimplified reduced costs. 118 virtual const VectorBase<R>& unsimplifiedRedCost() 119 { 120 assert(false); 121 static const VectorBase<R>& emptyVector = VectorBase<R>(); 122 return emptyVector; 123 } 124 125 /// gets basis status for a single row. 126 virtual typename SPxSolverBase<R>::VarStatus getBasisRowStatus(int i) const 127 { 128 assert(false); 129 return SPxSolverBase<R>::UNDEFINED; 130 } 131 132 /// gets basis status for a single column. 133 virtual typename SPxSolverBase<R>::VarStatus getBasisColStatus(int j) const 134 { 135 assert(false); 136 return SPxSolverBase<R>::UNDEFINED; 137 } 138 139 /// get optimal basis. 140 virtual void getBasis(typename SPxSolverBase<R>::VarStatus rows[], 141 typename SPxSolverBase<R>::VarStatus cols[], const int rowsSize, 142 const int colsSize) const 143 { 144 assert(false); 145 } 146 147 }; 148 149 } 150 151 152 #else 153 154 #include <memory> 155 156 #include "papilo/core/Presolve.hpp" 157 #include "papilo/core/ProblemBuilder.hpp" 158 #include "papilo/Config.hpp" 159 160 #include "soplex/spxsimplifier.h" 161 162 #include "soplex/spxdefines.h" 163 #include "soplex/spxsimplifier.h" 164 #include "soplex/array.h" 165 #include "soplex/exceptions.h" 166 #include "soplex/spxdefines.h" 167 168 169 namespace soplex 170 { 171 172 template<class R> 173 class Presol : public SPxSimplifier<R> 174 { 175 private: 176 177 #ifdef SOPLEX_DEBUG 178 const papilo::VerbosityLevel verbosityLevel = papilo::VerbosityLevel::kInfo; 179 #else 180 const papilo::VerbosityLevel verbosityLevel = papilo::VerbosityLevel::kQuiet; 181 #endif 182 183 VectorBase<R> m_prim; ///< unsimplified primal solution VectorBase<R>. 184 VectorBase<R> m_slack; ///< unsimplified slack VectorBase<R>. 185 VectorBase<R> m_dual; ///< unsimplified dual solution VectorBase<R>. 186 VectorBase<R> m_redCost; ///< unsimplified reduced cost VectorBase<R>. 187 DataArray<typename SPxSolverBase<R>::VarStatus> m_cBasisStat; ///< basis status of columns. 188 DataArray<typename SPxSolverBase<R>::VarStatus> m_rBasisStat; ///< basis status of rows. 189 190 191 papilo::PostsolveStorage<R> 192 postsolveStorage; ///< stored postsolve to recalculate the original solution 193 bool noChanges = false; ///< did PaPILO reduce the problem? 194 195 bool postsolved; ///< was the solution already postsolve? 196 bool vanished = false; 197 R modifyRowsFac; ///< 198 DataArray<int> m_stat; ///< preprocessing history. 199 typename SPxLPBase<R>::SPxSense m_thesense; ///< optimization sense. 200 201 bool enableSingletonCols; 202 bool enablePropagation; 203 bool enableParallelRows; 204 bool enableParallelCols; 205 bool enableSingletonStuffing; 206 bool enableDualFix; 207 bool enableFixContinuous; 208 bool enableDominatedCols; 209 210 // TODO: the following parameters were ignored? Maybe I don't exactly know what they suppose to be 211 bool m_keepbounds; ///< keep some bounds (for boundflipping) 212 typename SPxSimplifier<R>::Result m_result; ///< result of the simplification. 213 214 protected: 215 216 R epsZero() const 217 { 218 return this->tolerances()->epsilon(); 219 } 220 221 R feastol() const 222 { 223 return this->tolerances()->floatingPointFeastol(); 224 } 225 226 R opttol() const 227 { 228 return this->tolerances()->floatingPointOpttol(); 229 } 230 231 public: 232 233 //------------------------------------ 234 ///@name Constructors / destructors 235 ///@{ 236 /// default constructor. 237 explicit Presol(Timer::TYPE ttype = Timer::USER_TIME) 238 : SPxSimplifier<R>("PaPILO", ttype), postsolved(false), modifyRowsFac(1.0), 239 m_thesense(SPxLPBase<R>::MAXIMIZE), 240 m_keepbounds(false), m_result(this->OKAY) 241 { ; }; 242 243 /// copy constructor. 244 Presol(const Presol& old) 245 : SPxSimplifier<R>(old), m_prim(old.m_prim), m_slack(old.m_slack), m_dual(old.m_dual), 246 m_redCost(old.m_redCost), m_cBasisStat(old.m_cBasisStat), m_rBasisStat(old.m_rBasisStat), 247 postsolveStorage(old.postsolveStorage), postsolved(old.postsolved), 248 modifyRowsFac(old.modifyRowsFac), m_thesense(old.m_thesense), 249 m_keepbounds(old.m_keepbounds), m_result(old.m_result) 250 { 251 ; 252 } 253 254 /// assignment operator 255 Presol& operator=(const Presol& rhs) 256 { 257 if(this != &rhs) 258 { 259 SPxSimplifier<R>::operator=(rhs); 260 m_prim = rhs.m_prim; 261 m_slack = rhs.m_slack; 262 m_dual = rhs.m_dual; 263 m_redCost = rhs.m_redCost; 264 m_cBasisStat = rhs.m_cBasisStat; 265 m_rBasisStat = rhs.m_rBasisStat; 266 postsolved = rhs.postsolved; 267 m_thesense = rhs.m_thesense; 268 m_keepbounds = rhs.m_keepbounds; 269 m_result = rhs.m_result; 270 postsolveStorage = rhs.postsolveStorage; 271 modifyRowsFac = rhs.modifyRowsFac; 272 } 273 274 return *this; 275 } 276 277 /// destructor. 278 virtual ~Presol() 279 { 280 ; 281 } 282 283 SPxSimplifier<R>* clone() const 284 { 285 return new Presol(*this); 286 } 287 288 void 289 setModifyConsFrac(R value) 290 { 291 modifyRowsFac = value; 292 } 293 294 void 295 setEnableSingletonCols(bool value) 296 { 297 enableSingletonCols = value; 298 } 299 300 void 301 setEnablePropagation(bool value) 302 { 303 enablePropagation = value; 304 } 305 306 void 307 setEnableParallelRows(bool value) 308 { 309 enableParallelRows = value; 310 } 311 312 void 313 setEnableParallelCols(bool value) 314 { 315 enableParallelCols = value; 316 } 317 318 void 319 setEnableStuffing(bool value) 320 { 321 enableSingletonStuffing = value; 322 } 323 324 void 325 setEnableDualFix(bool value) 326 { 327 enableDualFix = value; 328 } 329 330 void 331 setEnableFixContinuous(bool value) 332 { 333 enableFixContinuous = value; 334 } 335 336 void 337 setEnableDomCols(bool value) 338 { 339 enableDominatedCols = value; 340 } 341 342 virtual typename SPxSimplifier<R>::Result simplify(SPxLPBase<R>& lp, 343 Real remainingTime, bool keepbounds, uint32_t seed); 344 345 virtual void unsimplify(const VectorBase<R>&, const VectorBase<R>&, const VectorBase<R>&, 346 const VectorBase<R>&, 347 const typename SPxSolverBase<R>::VarStatus[], 348 const typename SPxSolverBase<R>::VarStatus[], 349 bool isOptimal); 350 351 /// returns result status of the simplification 352 virtual typename SPxSimplifier<R>::Result result() const 353 { 354 return m_result; 355 } 356 357 /// specifies whether an optimal solution has already been unsimplified. 358 virtual bool isUnsimplified() const 359 { 360 return postsolved; 361 } 362 363 /// returns a reference to the unsimplified primal solution. 364 virtual const VectorBase<R>& unsimplifiedPrimal() 365 { 366 assert(postsolved); 367 return m_prim; 368 } 369 370 /// returns a reference to the unsimplified dual solution. 371 virtual const VectorBase<R>& unsimplifiedDual() 372 { 373 assert(postsolved); 374 return m_dual; 375 } 376 377 /// returns a reference to the unsimplified slack values. 378 virtual const VectorBase<R>& unsimplifiedSlacks() 379 { 380 assert(postsolved); 381 return m_slack; 382 } 383 384 /// returns a reference to the unsimplified reduced costs. 385 virtual const VectorBase<R>& unsimplifiedRedCost() 386 { 387 assert(postsolved); 388 return m_redCost; 389 } 390 391 /// gets basis status for a single row. 392 virtual typename SPxSolverBase<R>::VarStatus getBasisRowStatus(int i) const 393 { 394 assert(postsolved); 395 return m_rBasisStat[i]; 396 } 397 398 /// gets basis status for a single column. 399 virtual typename SPxSolverBase<R>::VarStatus getBasisColStatus(int j) const 400 { 401 assert(postsolved); 402 return m_cBasisStat[j]; 403 } 404 405 /// get optimal basis. 406 virtual void getBasis(typename SPxSolverBase<R>::VarStatus rows[], 407 typename SPxSolverBase<R>::VarStatus cols[], const int rowsSize = -1, 408 const int colsSize = -1) const 409 { 410 assert(postsolved); 411 assert(rowsSize < 0 || rowsSize >= m_rBasisStat.size()); 412 assert(colsSize < 0 || colsSize >= m_cBasisStat.size()); 413 414 for(int i = 0; i < m_rBasisStat.size(); ++i) 415 rows[i] = m_rBasisStat[i]; 416 417 for(int j = 0; j < m_cBasisStat.size(); ++j) 418 cols[j] = m_cBasisStat[j]; 419 } 420 421 private: 422 423 void initLocalVariables(const SPxLPBase <R>& lp); 424 425 void configurePapilo(papilo::Presolve<R>& presolve, R feasTolerance, R epsilon, uint32_t seed, 426 Real remainingTime) const; 427 428 void applyPresolveResultsToColumns(SPxLPBase <R>& lp, const papilo::Problem<R>& problem, 429 const papilo::PresolveResult<R>& res) const; 430 431 void applyPresolveResultsToRows(SPxLPBase <R>& lp, const papilo::Problem<R>& problem, 432 const papilo::PresolveResult<R>& res) const; 433 434 papilo::VarBasisStatus 435 convertToPapiloStatus(typename SPxSolverBase<R>::VarStatus status) const; 436 437 typename SPxSolverBase<R>::VarStatus 438 convertToSoplexStatus(papilo::VarBasisStatus status) const ; 439 }; 440 441 template <class R> 442 void Presol<R>::unsimplify(const VectorBase<R>& x, const VectorBase<R>& y, 443 const VectorBase<R>& s, const VectorBase<R>& r, 444 const typename SPxSolverBase<R>::VarStatus rows[], 445 const typename SPxSolverBase<R>::VarStatus cols[], 446 bool isOptimal) 447 { 448 449 SPX_MSG_INFO1((*this->spxout), (*this->spxout) 450 << " --- unsimplifying solution and basis" 451 << std::endl; 452 ) 453 454 assert(x.dim() <= m_prim.dim()); 455 assert(y.dim() <= m_dual.dim()); 456 assert(x.dim() == r.dim()); 457 assert(y.dim() == s.dim()); 458 459 //if presolving made no changes then copy the reduced solution to the original 460 if(noChanges) 461 { 462 for(int j = 0; j < x.dim(); ++j) 463 { 464 m_prim[j] = x[j]; 465 m_redCost[j] = r[j]; 466 m_cBasisStat[j] = cols[j]; 467 } 468 469 for(int i = 0; i < y.dim(); ++i) 470 { 471 m_dual[i] = y[i]; 472 m_slack[i] = s[i]; 473 m_rBasisStat[i] = rows[i]; 474 } 475 476 postsolved = true; 477 return; 478 } 479 480 int nColsReduced = (int)postsolveStorage.origcol_mapping.size(); 481 int nRowsReduced = (int)postsolveStorage.origrow_mapping.size(); 482 assert(x.dim() == (int)postsolveStorage.origcol_mapping.size() || vanished); 483 assert(y.dim() == (int)postsolveStorage.origrow_mapping.size() || vanished); 484 485 papilo::Solution<R> originalSolution{}; 486 papilo::Solution<R> reducedSolution{}; 487 reducedSolution.type = papilo::SolutionType::kPrimalDual; 488 reducedSolution.basisAvailabe = true; 489 490 reducedSolution.primal.clear(); 491 reducedSolution.reducedCosts.clear(); 492 reducedSolution.varBasisStatus.clear(); 493 reducedSolution.dual.clear(); 494 reducedSolution.rowBasisStatus.clear(); 495 496 reducedSolution.primal.resize(nColsReduced); 497 reducedSolution.reducedCosts.resize(nColsReduced); 498 reducedSolution.varBasisStatus.resize(nColsReduced); 499 reducedSolution.dual.resize(nRowsReduced); 500 reducedSolution.rowBasisStatus.resize(nRowsReduced); 501 502 postsolved = true; 503 504 // NOTE: for maximization problems, we have to switch signs of dual and 505 // reduced cost values, since simplifier assumes minimization problem 506 R switch_sign = m_thesense == SPxLPBase<R>::MAXIMIZE ? -1 : 1; 507 508 // assign values of variables in reduced LP 509 for(int j = 0; j < nColsReduced; ++j) 510 { 511 reducedSolution.primal[j] = isZero(x[j], this->epsZero()) ? 0.0 : x[j]; 512 reducedSolution.reducedCosts[j] = 513 isZero(r[j], this->epsZero()) ? 0.0 : switch_sign * r[j]; 514 reducedSolution.varBasisStatus[j] = convertToPapiloStatus(cols[j]); 515 } 516 517 for(int i = 0; i < nRowsReduced; ++i) 518 { 519 reducedSolution.dual[i] = isZero(y[i], this->epsZero()) ? 0.0 : switch_sign * y[i]; 520 reducedSolution.rowBasisStatus[i] = convertToPapiloStatus(rows[i]); 521 } 522 523 papilo::Num<R> num {}; 524 num.setEpsilon(this->epsZero()); 525 num.setFeasTol(this->feastol()); 526 /* since PaPILO verbosity is quiet it's irrelevant what the messenger is */ 527 papilo::Message msg{}; 528 msg.setVerbosityLevel(verbosityLevel); 529 530 papilo::Postsolve<R> postsolve {msg, num}; 531 auto status = postsolve.undo(reducedSolution, originalSolution, postsolveStorage, isOptimal); 532 533 if(status == PostsolveStatus::kFailed && isOptimal) 534 { 535 SPX_MSG_ERROR(std::cerr << "PaPILO did not pass validation" << std::endl;) 536 assert(false); 537 } 538 539 for(int j = 0; j < (int)postsolveStorage.nColsOriginal; ++j) 540 { 541 m_prim[j] = originalSolution.primal[j]; 542 m_redCost[j] = switch_sign * originalSolution.reducedCosts[j]; 543 m_cBasisStat[j] = convertToSoplexStatus(originalSolution.varBasisStatus[j]); 544 } 545 546 for(int i = 0; i < (int)postsolveStorage.nRowsOriginal; ++i) 547 { 548 m_dual[i] = switch_sign * originalSolution.dual[i]; 549 m_slack[i] = originalSolution.slack[i]; 550 m_rBasisStat[i] = convertToSoplexStatus(originalSolution.rowBasisStatus[i]); 551 } 552 553 } 554 555 template <class R> 556 papilo::VarBasisStatus 557 Presol<R>::convertToPapiloStatus(const typename SPxSolverBase<R>::VarStatus status) const 558 { 559 switch(status) 560 { 561 case SPxSolverBase<R>::ON_UPPER: 562 return papilo::VarBasisStatus::ON_UPPER; 563 564 case SPxSolverBase<R>::ON_LOWER: 565 return papilo::VarBasisStatus::ON_LOWER; 566 567 case SPxSolverBase<R>::FIXED: 568 return papilo::VarBasisStatus::FIXED; 569 570 case SPxSolverBase<R>::BASIC: 571 return papilo::VarBasisStatus::BASIC; 572 573 case SPxSolverBase<R>::UNDEFINED: 574 return papilo::VarBasisStatus::UNDEFINED; 575 576 case SPxSolverBase<R>::ZERO: 577 return papilo::VarBasisStatus::ZERO; 578 } 579 580 return papilo::VarBasisStatus::UNDEFINED; 581 } 582 583 template <class R> 584 typename SPxSolverBase<R>::VarStatus 585 Presol<R>::convertToSoplexStatus(papilo::VarBasisStatus status) const 586 { 587 switch(status) 588 { 589 case papilo::VarBasisStatus::ON_UPPER: 590 return SPxSolverBase<R>::ON_UPPER; 591 592 case papilo::VarBasisStatus::ON_LOWER: 593 return SPxSolverBase<R>::ON_LOWER; 594 595 case papilo::VarBasisStatus::ZERO: 596 return SPxSolverBase<R>::ZERO; 597 598 case papilo::VarBasisStatus::FIXED: 599 return SPxSolverBase<R>::FIXED; 600 601 case papilo::VarBasisStatus::UNDEFINED: 602 return SPxSolverBase<R>::UNDEFINED; 603 604 case papilo::VarBasisStatus::BASIC: 605 return SPxSolverBase<R>::BASIC; 606 } 607 608 return SPxSolverBase<R>::UNDEFINED; 609 } 610 611 612 template<class R> 613 papilo::Problem<R> buildProblem(SPxLPBase<R>& lp) 614 { 615 papilo::ProblemBuilder<R> builder; 616 617 /* build problem from matrix */ 618 int nnz = lp.nNzos(); 619 int ncols = lp.nCols(); 620 int nrows = lp.nRows(); 621 builder.reserve(nnz, nrows, ncols); 622 623 /* set up columns */ 624 builder.setNumCols(ncols); 625 626 R switch_sign = lp.spxSense() == SPxLPBase<R>::MAXIMIZE ? -1 : 1; 627 628 for(int i = 0; i < ncols; ++i) 629 { 630 R lowerbound = lp.lower(i); 631 R upperbound = lp.upper(i); 632 R objective = lp.obj(i); 633 builder.setColLb(i, lowerbound); 634 builder.setColUb(i, upperbound); 635 builder.setColLbInf(i, lowerbound <= -R(infinity)); 636 builder.setColUbInf(i, upperbound >= R(infinity)); 637 638 builder.setColIntegral(i, false); 639 builder.setObj(i, objective * switch_sign); 640 } 641 642 /* set up rows */ 643 builder.setNumRows(nrows); 644 645 for(int i = 0; i < nrows; ++i) 646 { 647 const SVectorBase<R> rowVector = lp.rowVector(i); 648 int rowlength = rowVector.size(); 649 int* indices = new int[rowlength]; 650 R* rowValues = new R[rowlength]; 651 652 for(int j = 0; j < rowlength; j++) 653 { 654 const Nonzero<R> element = rowVector.element(j); 655 indices[j] = element.idx; 656 rowValues[j] = element.val; 657 } 658 659 builder.addRowEntries(i, rowlength, indices, rowValues); 660 661 R lhs = lp.lhs(i); 662 R rhs = lp.rhs(i); 663 builder.setRowLhs(i, lhs); 664 builder.setRowRhs(i, rhs); 665 builder.setRowLhsInf(i, lhs <= -R(infinity)); 666 builder.setRowRhsInf(i, rhs >= R(infinity)); 667 } 668 669 return builder.build(); 670 } 671 672 673 template<class R> 674 typename SPxSimplifier<R>::Result 675 Presol<R>::simplify(SPxLPBase<R>& lp, Real remainingTime, bool keepbounds, uint32_t seed) 676 { 677 678 //TODO: how to use the keepbounds parameter? 679 m_keepbounds = keepbounds; 680 681 if(m_keepbounds) 682 SPX_MSG_WARNING((*this->spxout), 683 (*this->spxout) << "==== PaPILO doesn't handle parameter keepbounds" << 684 std::endl;) 685 686 initLocalVariables(lp); 687 688 papilo::Problem<R> problem = buildProblem(lp); 689 papilo::Presolve<R> presolve; 690 691 configurePapilo(presolve, this->tolerances()->floatingPointFeastol(), this->tolerances()->epsilon(), 692 seed, remainingTime); 693 SPX_MSG_INFO1((*this->spxout), (*this->spxout) 694 << " --- starting PaPILO" << std::endl; 695 ) 696 697 papilo::PresolveResult<R> res = presolve.apply(problem); 698 699 assert(res.postsolve.postsolveType == PostsolveType::kFull); 700 701 switch(res.status) 702 { 703 case papilo::PresolveStatus::kInfeasible: 704 m_result = SPxSimplifier<R>::INFEASIBLE; 705 SPX_MSG_INFO1((*this->spxout), (*this->spxout) 706 << " --- presolving detected infeasibility" << std::endl; 707 ) 708 return SPxSimplifier<R>::INFEASIBLE; 709 710 case papilo::PresolveStatus::kUnbndOrInfeas: 711 case papilo::PresolveStatus::kUnbounded: 712 m_result = SPxSimplifier<R>::UNBOUNDED; 713 SPX_MSG_INFO1((*this->spxout), (*this->spxout) << 714 "==== Presolving detected unboundedness of the problem" << std::endl; 715 ) 716 return SPxSimplifier<R>::UNBOUNDED; 717 718 case papilo::PresolveStatus::kUnchanged: 719 // since Soplex has no state unchanged store the value in a new variable 720 noChanges = true; 721 SPX_MSG_INFO1((*this->spxout), (*this->spxout) 722 << "==== Presolving found nothing " << std::endl; 723 ) 724 return SPxSimplifier<R>::OKAY; 725 726 case papilo::PresolveStatus::kReduced: 727 break; 728 } 729 730 731 int newNonzeros = problem.getConstraintMatrix().getNnz(); 732 733 if(newNonzeros == 0 || ((problem.getNRows() <= modifyRowsFac * lp.nRows() || 734 newNonzeros <= modifyRowsFac * lp.nNzos()))) 735 { 736 SPX_MSG_INFO1((*this->spxout), (*this->spxout) 737 << " --- presolved problem has " << problem.getNRows() << 738 " rows, " 739 << problem.getNCols() << " cols and " 740 << newNonzeros << " non-zeros and " 741 << presolve.getStatistics().nboundchgs << " boundchanges and " 742 << presolve.getStatistics().nsidechgs << " sidechanges" 743 << std::endl; 744 ) 745 postsolveStorage = res.postsolve; 746 747 // remove all constraints and variables 748 for(int j = lp.nCols() - 1; j >= 0; j--) 749 lp.removeCol(j); 750 751 for(int i = lp.nRows() - 1; i >= 0; i--) 752 lp.removeRow(i); 753 754 applyPresolveResultsToColumns(lp, problem, res); 755 applyPresolveResultsToRows(lp, problem, res); 756 assert(newNonzeros == lp.nNzos()); 757 } 758 else 759 { 760 noChanges = true; 761 SPX_MSG_INFO1((*this->spxout), 762 (*this->spxout) 763 764 << " --- presolve results smaller than the modifyconsfac" 765 << std::endl; 766 ) 767 } 768 769 if(newNonzeros == 0) 770 { 771 vanished = true; 772 m_result = SPxSimplifier<R>::VANISHED; 773 } 774 775 return m_result; 776 } 777 778 template<class R> 779 void Presol<R>::initLocalVariables(const SPxLPBase <R>& lp) 780 { 781 m_result = SPxSimplifier<R>::OKAY; 782 783 m_thesense = lp.spxSense(); 784 postsolved = false; 785 786 m_prim.reDim(lp.nCols()); 787 m_slack.reDim(lp.nRows()); 788 m_dual.reDim(lp.nRows()); 789 m_redCost.reDim(lp.nCols()); 790 m_cBasisStat.reSize(lp.nCols()); 791 m_rBasisStat.reSize(lp.nRows()); 792 793 this->m_timeUsed->reset(); 794 this->m_timeUsed->start(); 795 } 796 797 template<class R> 798 void Presol<R>::configurePapilo(papilo::Presolve<R>& presolve, R feasTolerance, R epsilon, 799 uint32_t seed, Real remainingTime) const 800 { 801 /* communicate the SOPLEX parameters to the presolve libary */ 802 803 /* communicate the random seed */ 804 presolve.getPresolveOptions().randomseed = (unsigned int) seed; 805 806 /* set number of threads to be used for presolve */ 807 /* TODO: set threads for PaPILO? Can Soplex be run with multiple threads?*/ 808 // presolve.getPresolveOptions().threads = data->threads; 809 810 presolve.getPresolveOptions().tlim = remainingTime; 811 presolve.getPresolveOptions().feastol = double(feasTolerance); 812 presolve.getPresolveOptions().epsilon = double(epsilon); 813 presolve.getPresolveOptions().detectlindep = 0; 814 presolve.getPresolveOptions().componentsmaxint = -1; 815 presolve.getPresolveOptions().calculate_basis_for_dual = true; 816 817 presolve.setVerbosityLevel(verbosityLevel); 818 819 /* enable lp presolvers with dual postsolve*/ 820 using uptr = std::unique_ptr<papilo::PresolveMethod<R>>; 821 822 /* fast presolvers*/ 823 if(enableSingletonCols) 824 presolve.addPresolveMethod(uptr(new papilo::SingletonCols<R>())); 825 826 if(enablePropagation) 827 presolve.addPresolveMethod(uptr(new papilo::ConstraintPropagation<R>())); 828 829 /* medium presolver */ 830 if(enableParallelRows) 831 presolve.addPresolveMethod(uptr(new papilo::ParallelRowDetection<R>())); 832 833 if(enableParallelCols) 834 presolve.addPresolveMethod(uptr(new papilo::ParallelColDetection<R>())); 835 836 if(enableSingletonStuffing) 837 presolve.addPresolveMethod(uptr(new papilo::SingletonStuffing<R>())); 838 839 if(enableDualFix) 840 presolve.addPresolveMethod(uptr(new papilo::DualFix<R>())); 841 842 if(enableFixContinuous) 843 presolve.addPresolveMethod(uptr(new papilo::FixContinuous<R>())); 844 845 /* exhaustive presolvers*/ 846 if(enableDominatedCols) 847 presolve.addPresolveMethod(uptr(new papilo::DominatedCols<R>())); 848 849 /** 850 * TODO: PaPILO doesn't support dualpostsolve for those presolvers 851 * presolve.addPresolveMethod(uptr(new papilo::SimpleSubstitution<R>())); 852 * presolve.addPresolveMethod(uptr(new papilo::DualInfer<R>())); 853 * presolve.addPresolveMethod(uptr(new papilo::Substitution<R>())); 854 * presolve.addPresolveMethod(uptr(new papilo::Sparsify<R>())); 855 * presolve.getPresolveOptions().removeslackvars = false; 856 * presolve.getPresolveOptions().maxfillinpersubstitution 857 * =data->maxfillinpersubstitution; 858 * presolve.getPresolveOptions().maxshiftperrow = data->maxshiftperrow; 859 */ 860 861 862 } 863 864 template<class R> 865 void Presol<R>::applyPresolveResultsToColumns(SPxLPBase <R>& lp, const papilo::Problem<R>& problem, 866 const papilo::PresolveResult<R>& res) const 867 { 868 869 const papilo::Objective<R>& objective = problem.getObjective(); 870 const papilo::Vec<R>& upperBounds = problem.getUpperBounds(); 871 const papilo::Vec<R>& lowerBounds = problem.getLowerBounds(); 872 const papilo::Vec<papilo::ColFlags>& colFlags = problem.getColFlags(); 873 874 R switch_sign = lp.spxSense() == SPxLPBase<R>::MAXIMIZE ? -1 : 1; 875 876 for(int col = 0; col < problem.getNCols(); col++) 877 { 878 DSVectorBase<R> emptyVector{0}; 879 R lb = lowerBounds[col]; 880 881 if(colFlags[col].test(papilo::ColFlag::kLbInf)) 882 lb = -R(infinity); 883 884 R ub = upperBounds[col]; 885 886 if(colFlags[col].test(papilo::ColFlag::kUbInf)) 887 ub = R(infinity); 888 889 LPColBase<R> column(objective.coefficients[col]* switch_sign, emptyVector, ub, lb); 890 lp.addCol(column); 891 assert(lp.lower(col) == lb); 892 assert(lp.upper(col) == ub); 893 } 894 895 lp.changeObjOffset(objective.offset); 896 897 assert(problem.getNCols() == lp.nCols()); 898 } 899 900 template<class R> 901 void Presol<R>::applyPresolveResultsToRows(SPxLPBase <R>& lp, const papilo::Problem<R>& problem, 902 const papilo::PresolveResult<R>& res) const 903 { 904 int size = res.postsolve.origrow_mapping.size(); 905 906 //add the adjusted constraints 907 for(int row = 0; row < size; row++) 908 { 909 R rhs = problem.getConstraintMatrix().getRightHandSides()[row]; 910 911 if(problem.getRowFlags()[row].test(papilo::RowFlag::kRhsInf)) 912 rhs = R(infinity); 913 914 R lhs = problem.getConstraintMatrix().getLeftHandSides()[row]; 915 916 if(problem.getRowFlags()[row].test(papilo::RowFlag::kLhsInf)) 917 lhs = -R(infinity); 918 919 const papilo::SparseVectorView<R> papiloRowVector = 920 problem.getConstraintMatrix().getRowCoefficients(row); 921 const int* indices = papiloRowVector.getIndices(); 922 const R* values = papiloRowVector.getValues(); 923 924 int length = papiloRowVector.getLength(); 925 DSVectorBase<R> soplexRowVector{length}; 926 927 for(int i = 0; i < length; i++) 928 { 929 soplexRowVector.add(indices[i], values[i]); 930 } 931 932 LPRowBase<R> lpRowBase(lhs, soplexRowVector, rhs); 933 lp.addRow(lpRowBase); 934 assert(lp.lhs(row) == lhs); 935 assert(lp.rhs(row) == rhs); 936 } 937 938 assert(problem.getNRows() == lp.nRows()); 939 } 940 941 } // namespace soplex 942 943 #endif 944