26#ifndef SOPLEX_WITH_PAPILO
163#include "papilo/core/Presolve.hpp"
164#include "papilo/core/ProblemBuilder.hpp"
165#include "papilo/Config.hpp"
180class Presol :
public SPxSimplifier<R>
185 const papilo::VerbosityLevel verbosityLevel = papilo::VerbosityLevel::kInfo;
187 const papilo::VerbosityLevel verbosityLevel = papilo::VerbosityLevel::kQuiet;
190 VectorBase<R> m_prim;
191 VectorBase<R> m_slack;
192 VectorBase<R> m_dual;
193 VectorBase<R> m_redCost;
194 DataArray<typename SPxSolverBase<R>::VarStatus> m_cBasisStat;
195 DataArray<typename SPxSolverBase<R>::VarStatus> m_rBasisStat;
198 papilo::PostsolveStorage<R>
200 bool noChanges =
false;
203 bool vanished =
false;
208 DataArray<int> m_stat;
241 m_thesense(SPxLPBase<R>::MAXIMIZE),
242 m_keepbounds(false), m_result(this->
OKAY)
247 :
SPxSimplifier<R>(old), m_prim(old.m_prim), m_slack(old.m_slack), m_dual(old.m_dual),
248 m_redCost(old.m_redCost), m_cBasisStat(old.m_cBasisStat), m_rBasisStat(old.m_rBasisStat),
249 postsolveStorage(old.postsolveStorage), postsolved(old.postsolved), m_epsilon(old.m_epsilon),
250 m_feastol(old.m_feastol), m_opttol(old.m_opttol),
251 modifyRowsFac(old.modifyRowsFac), m_thesense(old.m_thesense),
252 m_keepbounds(old.m_keepbounds), m_result(old.m_result)
264 m_slack = rhs.m_slack;
266 m_redCost = rhs.m_redCost;
267 m_cBasisStat = rhs.m_cBasisStat;
268 m_rBasisStat = rhs.m_rBasisStat;
269 postsolved = rhs.postsolved;
270 m_epsilon = rhs.m_epsilon;
271 m_feastol = rhs.m_feastol;
272 m_opttol = rhs.m_opttol;
273 m_thesense = rhs.m_thesense;
274 m_keepbounds = rhs.m_keepbounds;
275 m_result = rhs.m_result;
276 postsolveStorage = rhs.postsolveStorage;
277 modifyRowsFac = rhs.modifyRowsFac;
289 SPxSimplifier<R>*
clone()
const
295 setModifyConsFrac(R value)
297 modifyRowsFac = value;
303 return simplify(lp, eps, delta, delta, remainingTime,
false, 0);
308 bool keepbounds, uint32_t seed);
310 virtual void unsimplify(
const VectorBase<R>&,
const VectorBase<R>&,
const VectorBase<R>&,
311 const VectorBase<R>&,
360 return m_rBasisStat[i];
367 return m_cBasisStat[j];
373 const int colsSize = -1)
const
376 assert(rowsSize < 0 || rowsSize >= m_rBasisStat.size());
377 assert(colsSize < 0 || colsSize >= m_cBasisStat.size());
379 for(
int i = 0; i < m_rBasisStat.size(); ++i)
380 rows[i] = m_rBasisStat[i];
382 for(
int j = 0; j < m_cBasisStat.size(); ++j)
383 cols[j] = m_cBasisStat[j];
388 void initLocalVariables(
const SPxLPBase <R>& lp);
390 void configurePapilo(papilo::Presolve<R>& presolve, R feasTolerance, R epsilon, uint32_t seed,
391 Real remainingTime)
const;
393 void applyPresolveResultsToColumns(SPxLPBase <R>& lp,
const papilo::Problem<R>& problem,
394 const papilo::PresolveResult<R>& res)
const;
396 void applyPresolveResultsToRows(SPxLPBase <R>& lp,
const papilo::Problem<R>& problem,
397 const papilo::PresolveResult<R>& res)
const;
399 papilo::VarBasisStatus
403 convertToSoplexStatus(papilo::VarBasisStatus status)
const ;
407void Presol<R>::unsimplify(
const VectorBase<R>& x,
const VectorBase<R>& y,
408 const VectorBase<R>& s,
const VectorBase<R>& r,
409 const typename SPxSolverBase<R>::VarStatus rows[],
410 const typename SPxSolverBase<R>::VarStatus cols[],
414 MSG_INFO1((*this->spxout), (*this->spxout)
415 <<
" --- unsimplifying solution and basis"
419 assert(
x.dim() <= m_prim.dim());
420 assert(
y.dim() <= m_dual.dim());
427 for(
int j = 0;
j <
x.dim(); ++
j)
431 m_cBasisStat[
j] = cols[
j];
434 for(
int i = 0;
i <
y.dim(); ++
i)
438 m_rBasisStat[
i] = rows[
i];
488 papilo::Num<R> num {};
489 num.setEpsilon(m_epsilon);
490 num.setFeasTol(m_feastol);
492 papilo::Message msg{};
495 papilo::Postsolve<R>
postsolve {msg, num};
498 if(status == PostsolveStatus::kFailed &&
isOptimal)
500 MSG_ERROR(std::cerr <<
"PaPILO did not pass validation" << std::endl;)
521papilo::VarBasisStatus
522Presol<R>::convertToPapiloStatus(
const typename SPxSolverBase<R>::VarStatus status)
const
527 return papilo::VarBasisStatus::ON_UPPER;
530 return papilo::VarBasisStatus::ON_LOWER;
533 return papilo::VarBasisStatus::FIXED;
536 return papilo::VarBasisStatus::BASIC;
539 return papilo::VarBasisStatus::UNDEFINED;
542 return papilo::VarBasisStatus::ZERO;
545 return papilo::VarBasisStatus::UNDEFINED;
549typename SPxSolverBase<R>::VarStatus
550Presol<R>::convertToSoplexStatus(papilo::VarBasisStatus status)
const
554 case papilo::VarBasisStatus::ON_UPPER:
557 case papilo::VarBasisStatus::ON_LOWER:
560 case papilo::VarBasisStatus::ZERO:
563 case papilo::VarBasisStatus::FIXED:
566 case papilo::VarBasisStatus::UNDEFINED:
569 case papilo::VarBasisStatus::BASIC:
578papilo::Problem<R> buildProblem(SPxLPBase<R>& lp)
580 papilo::ProblemBuilder<R> builder;
583 int nnz = lp.nNzos();
584 int ncols = lp.nCols();
585 int nrows = lp.nRows();
586 builder.reserve(nnz, nrows, ncols);
589 builder.setNumCols(ncols);
591 R switch_sign = lp.spxSense() == SPxLPBase<R>::MAXIMIZE ? -1 : 1;
593 for(
int i = 0; i < ncols; ++i)
595 R lowerbound = lp.lower(i);
596 R upperbound = lp.upper(i);
597 R objective = lp.obj(i);
598 builder.setColLb(i, lowerbound);
599 builder.setColUb(i, upperbound);
600 builder.setColLbInf(i, lowerbound <= -R(infinity));
601 builder.setColUbInf(i, upperbound >= R(infinity));
603 builder.setColIntegral(i,
false);
604 builder.setObj(i, objective * switch_sign);
608 builder.setNumRows(nrows);
610 for(
int i = 0; i < nrows; ++i)
612 const SVectorBase<R> rowVector = lp.rowVector(i);
613 int rowlength = rowVector.size();
614 int* indices =
new int[rowlength];
615 R* rowValues =
new R[rowlength];
617 for(
int j = 0; j < rowlength; j++)
619 const Nonzero<R> element = rowVector.element(j);
620 indices[j] = element.idx;
621 rowValues[j] = element.val;
624 builder.addRowEntries(i, rowlength, indices, rowValues);
628 builder.setRowLhs(i, lhs);
629 builder.setRowRhs(i, rhs);
630 builder.setRowLhsInf(i, lhs <= -R(infinity));
631 builder.setRowRhsInf(i, rhs >= R(infinity));
634 return builder.build();
639typename SPxSimplifier<R>::Result
640Presol<R>::simplify(SPxLPBase<R>& lp, R eps, R ftol, R otol,
641 Real remainingTime,
bool keepbounds, uint32_t seed)
648 MSG_WARNING((*this->spxout), (*this->spxout) <<
"==== PaPILO doesn't handle parameter keepbounds" <<
657 MSG_INFO1((*this->spxout), (*this->spxout)
658 <<
" --- starting PaPILO" << std::endl;
665 case papilo::PresolveStatus::kInfeasible:
667 MSG_INFO1((*this->spxout), (*this->spxout)
668 <<
" --- presolving detected infeasibility" << std::endl;
672 case papilo::PresolveStatus::kUnbndOrInfeas:
673 case papilo::PresolveStatus::kUnbounded:
675 MSG_INFO1((*this->spxout), (*this->spxout) <<
676 "==== Presolving detected unboundedness of the problem" << std::endl;
680 case papilo::PresolveStatus::kUnchanged:
683 MSG_INFO1((*this->spxout), (*this->spxout)
684 <<
"==== Presolving found nothing " << std::endl;
688 case papilo::PresolveStatus::kReduced:
698 MSG_INFO1((*this->spxout), (*this->spxout)
699 <<
" --- presolved problem has " <<
problem.getNRows() <<
701 <<
problem.getNCols() <<
" cols and "
703 <<
presolve.getStatistics().nboundchgs <<
" boundchanges and "
704 <<
presolve.getStatistics().nsidechgs <<
" sidechanges"
710 for(
int j =
lp.nCols() - 1;
j >= 0;
j--)
713 for(
int i =
lp.nRows() - 1;
i >= 0;
i--)
726 <<
" --- presolve results smaller than the modifyconsfac"
741void Presol<R>::initLocalVariables(
const SPxLPBase <R>& lp)
745 m_thesense =
lp.spxSense();
748 m_prim.reDim(
lp.nCols());
749 m_slack.reDim(
lp.nRows());
750 m_dual.reDim(
lp.nRows());
751 m_redCost.reDim(
lp.nCols());
752 m_cBasisStat.reSize(
lp.nCols());
753 m_rBasisStat.reSize(
lp.nRows());
755 this->m_timeUsed->reset();
756 this->m_timeUsed->start();
760void Presol<R>::configurePapilo(papilo::Presolve<R>& presolve, R feasTolerance, R epsilon,
761 uint32_t seed, Real remainingTime)
const
775 presolve.getPresolveOptions().detectlindep = 0;
776 presolve.getPresolveOptions().componentsmaxint = -1;
777 presolve.getPresolveOptions().calculate_basis_for_dual =
true;
782 using uptr = std::unique_ptr<papilo::PresolveMethod<R>>;
785 presolve.addPresolveMethod(
uptr(
new papilo::SingletonCols<R>()));
786 presolve.addPresolveMethod(
uptr(
new papilo::ConstraintPropagation<R>()));
789 presolve.addPresolveMethod(
uptr(
new papilo::ParallelRowDetection<R>()));
790 presolve.addPresolveMethod(
uptr(
new papilo::ParallelColDetection<R>()));
791 presolve.addPresolveMethod(
uptr(
new papilo::SingletonStuffing<R>()));
792 presolve.addPresolveMethod(
uptr(
new papilo::DualFix<R>()));
793 presolve.addPresolveMethod(
uptr(
new papilo::FixContinuous<R>()));
796 presolve.addPresolveMethod(
uptr(
new papilo::DominatedCols<R>()));
814void Presol<R>::applyPresolveResultsToColumns(SPxLPBase <R>& lp,
const papilo::Problem<R>& problem,
815 const papilo::PresolveResult<R>& res)
const
825 for(
int col = 0; col <
problem.getNCols(); col++)
830 if(
colFlags[col].test(papilo::ColFlag::kLbInf))
835 if(
colFlags[col].test(papilo::ColFlag::kUbInf))
850void Presol<R>::applyPresolveResultsToRows(SPxLPBase <R>& lp,
const papilo::Problem<R>& problem,
851 const papilo::PresolveResult<R>& res)
const
856 for(
int row = 0; row <
size; row++)
858 R rhs =
problem.getConstraintMatrix().getRightHandSides()[row];
860 if(
problem.getRowFlags()[row].test(papilo::RowFlag::kRhsInf))
863 R lhs =
problem.getConstraintMatrix().getLeftHandSides()[row];
865 if(
problem.getRowFlags()[row].test(papilo::RowFlag::kLhsInf))
869 problem.getConstraintMatrix().getRowCoefficients(row);
876 for(
int i = 0;
i < length;
i++)
Save arrays of arbitrary types.
Safe arrays of data objects.
DataArray(const DataArray &old)
copy constructor
int size() const
return nr. of elements.
virtual SPxSimplifier< R >::Result simplify(SPxLPBase< R > &lp, R eps, R ftol, R otol, Real remainingTime, bool keepbounds, uint32_t seed)
virtual void unsimplify(const VectorBase< R > &, const VectorBase< R > &, const VectorBase< R > &, const VectorBase< R > &, const typename SPxSolverBase< R >::VarStatus[], const typename SPxSolverBase< R >::VarStatus[], bool isOptimal)
Presol(Timer::TYPE ttype=Timer::USER_TIME)
Presol & operator=(const Presol &rhs)
assignment operator
virtual const VectorBase< R > & unsimplifiedSlacks()
returns a reference to the unsimplified slack values.
virtual const VectorBase< R > & unsimplifiedRedCost()
returns a reference to the unsimplified reduced costs.
SPxSimplifier< R > * clone() const
virtual ~Presol()
destructor.
virtual SPxSolverBase< R >::VarStatus getBasisColStatus(int j) const
gets basis status for a single column.
virtual SPxSimplifier< R >::Result simplify(SPxLPBase< R > &lp, R eps, R delta, Real remainingTime)
virtual SPxSimplifier< R >::Result result() const
returns result status of the simplification
virtual const VectorBase< R > & unsimplifiedDual()
returns a reference to the unsimplified dual solution.
virtual SPxSolverBase< R >::VarStatus getBasisRowStatus(int i) const
gets basis status for a single row.
virtual void getBasis(typename SPxSolverBase< R >::VarStatus rows[], typename SPxSolverBase< R >::VarStatus cols[], const int rowsSize=-1, const int colsSize=-1) const
get optimal basis.
Presol(const Presol &old)
copy constructor.
virtual const VectorBase< R > & unsimplifiedPrimal()
returns a reference to the unsimplified primal solution.
virtual bool isUnsimplified() const
specifies whether an optimal solution has already been unsimplified.
SPxSense
Optimization sense.
LP simplification abstract base class.
Result
Result of the simplification.
@ INFEASIBLE
primal infeasibility was detected
@ UNBOUNDED
primal unboundedness was detected
@ OKAY
simplification could be done
@ VANISHED
the problem was so much simplified that it vanished
SPxSimplifier & operator=(const SPxSimplifier &rhs)
assignment operator
SPxSimplifier(const char *p_name, Timer::TYPE ttype=Timer::USER_TIME)
constructor
Sequential object-oriented SimPlex.
@ BASIC
variable is basic.
@ ON_LOWER
variable set to its lower bound.
@ ON_UPPER
variable set to its upper bound.
@ UNDEFINED
nothing known about basis status (possibly due to a singular basis in transformed problem)
@ FIXED
variable fixed to identical bounds.
@ ZERO
free variable fixed to zero.
Exception classes for SoPlex.
Everything should be within this namespace.
THREADLOCAL const Real infinity
Debugging, floating point type and parameter definitions.
#define MSG_WARNING(spxout, x)
Prints out message x if the verbosity level is at least SPxOut::WARNING.
#define DEFAULT_BND_VIOL
default allowed bound violation
#define DEFAULT_EPS_ZERO
default allowed additive zero: 1.0 + EPS_ZERO == 1.0
#define MSG_INFO1(spxout, x)
Prints out message x if the verbosity level is at least SPxOut::INFO1.
#define MSG_ERROR(x)
Prints out message x if the verbosity level is at least SPxOut::ERROR.
LP simplification base class.