26 #ifndef SOPLEX_WITH_PAPILO 64 Real remainingTime,
bool keepbounds, uint32_t seed)
142 const int colsSize)
const 156 #include "papilo/core/Presolve.hpp" 157 #include "papilo/core/ProblemBuilder.hpp" 158 #include "papilo/Config.hpp" 178 const papilo::VerbosityLevel verbosityLevel = papilo::VerbosityLevel::kInfo;
180 const papilo::VerbosityLevel verbosityLevel = papilo::VerbosityLevel::kQuiet;
191 papilo::PostsolveStorage<R>
193 bool noChanges =
false;
196 bool vanished =
false;
201 bool enableSingletonCols;
202 bool enablePropagation;
203 bool enableParallelRows;
204 bool enableParallelCols;
205 bool enableSingletonStuffing;
207 bool enableFixContinuous;
208 bool enableDominatedCols;
223 return this->
tolerances()->floatingPointFeastol();
228 return this->
tolerances()->floatingPointOpttol();
240 m_keepbounds(
false), m_result(this->
OKAY)
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)
261 m_slack = rhs.m_slack;
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;
289 setModifyConsFrac(R value)
291 modifyRowsFac = value;
295 setEnableSingletonCols(
bool value)
297 enableSingletonCols = value;
301 setEnablePropagation(
bool value)
303 enablePropagation = value;
307 setEnableParallelRows(
bool value)
309 enableParallelRows = value;
313 setEnableParallelCols(
bool value)
315 enableParallelCols = value;
319 setEnableStuffing(
bool value)
321 enableSingletonStuffing = value;
325 setEnableDualFix(
bool value)
327 enableDualFix = value;
331 setEnableFixContinuous(
bool value)
333 enableFixContinuous = value;
337 setEnableDomCols(
bool value)
339 enableDominatedCols = value;
343 Real remainingTime,
bool keepbounds, uint32_t seed);
395 return m_rBasisStat[i];
402 return m_cBasisStat[j];
408 const int colsSize = -1)
const 411 assert(rowsSize < 0 || rowsSize >= m_rBasisStat.
size());
412 assert(colsSize < 0 || colsSize >= m_cBasisStat.
size());
414 for(
int i = 0; i < m_rBasisStat.
size(); ++i)
415 rows[i] = m_rBasisStat[i];
417 for(
int j = 0; j < m_cBasisStat.
size(); ++j)
418 cols[j] = m_cBasisStat[j];
425 void configurePapilo(papilo::Presolve<R>& presolve, R feasTolerance, R epsilon, uint32_t seed,
426 Real remainingTime)
const;
428 void applyPresolveResultsToColumns(
SPxLPBase <R>& lp,
const papilo::Problem<R>& problem,
429 const papilo::PresolveResult<R>& res)
const;
431 void applyPresolveResultsToRows(
SPxLPBase <R>& lp,
const papilo::Problem<R>& problem,
432 const papilo::PresolveResult<R>& res)
const;
434 papilo::VarBasisStatus
438 convertToSoplexStatus(papilo::VarBasisStatus status)
const ;
450 <<
" --- unsimplifying solution and basis" 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());
462 for(
int j = 0; j < x.
dim(); ++j)
466 m_cBasisStat[j] = cols[j];
469 for(
int i = 0; i < y.
dim(); ++i)
473 m_rBasisStat[i] = rows[i];
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);
485 papilo::Solution<R> originalSolution{};
486 papilo::Solution<R> reducedSolution{};
487 reducedSolution.type = papilo::SolutionType::kPrimalDual;
488 reducedSolution.basisAvailabe =
true;
490 reducedSolution.primal.clear();
491 reducedSolution.reducedCosts.clear();
492 reducedSolution.varBasisStatus.clear();
493 reducedSolution.dual.clear();
494 reducedSolution.rowBasisStatus.clear();
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);
509 for(
int j = 0; j < nColsReduced; ++j)
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]);
517 for(
int i = 0; i < nRowsReduced; ++i)
519 reducedSolution.dual[i] = isZero(y[i], this->epsZero()) ? 0.0 : switch_sign * y[i];
520 reducedSolution.rowBasisStatus[i] = convertToPapiloStatus(rows[i]);
523 papilo::Num<R> num {};
524 num.setEpsilon(this->epsZero());
525 num.setFeasTol(this->feastol());
527 papilo::Message msg{};
528 msg.setVerbosityLevel(verbosityLevel);
530 papilo::Postsolve<R> postsolve {msg, num};
531 auto status = postsolve.undo(reducedSolution, originalSolution, postsolveStorage, isOptimal);
533 if(status == PostsolveStatus::kFailed && isOptimal)
535 SPX_MSG_ERROR(std::cerr <<
"PaPILO did not pass validation" << std::endl;)
539 for(
int j = 0; j < (int)postsolveStorage.nColsOriginal; ++j)
541 m_prim[j] = originalSolution.primal[j];
542 m_redCost[j] = switch_sign * originalSolution.reducedCosts[j];
543 m_cBasisStat[j] = convertToSoplexStatus(originalSolution.varBasisStatus[j]);
546 for(
int i = 0; i < (int)postsolveStorage.nRowsOriginal; ++i)
548 m_dual[i] = switch_sign * originalSolution.dual[i];
549 m_slack[i] = originalSolution.slack[i];
550 m_rBasisStat[i] = convertToSoplexStatus(originalSolution.rowBasisStatus[i]);
556 papilo::VarBasisStatus
562 return papilo::VarBasisStatus::ON_UPPER;
565 return papilo::VarBasisStatus::ON_LOWER;
568 return papilo::VarBasisStatus::FIXED;
571 return papilo::VarBasisStatus::BASIC;
574 return papilo::VarBasisStatus::UNDEFINED;
577 return papilo::VarBasisStatus::ZERO;
580 return papilo::VarBasisStatus::UNDEFINED;
589 case papilo::VarBasisStatus::ON_UPPER:
592 case papilo::VarBasisStatus::ON_LOWER:
595 case papilo::VarBasisStatus::ZERO:
598 case papilo::VarBasisStatus::FIXED:
601 case papilo::VarBasisStatus::UNDEFINED:
604 case papilo::VarBasisStatus::BASIC:
615 papilo::ProblemBuilder<R> builder;
618 int nnz = lp.
nNzos();
619 int ncols = lp.
nCols();
620 int nrows = lp.
nRows();
621 builder.reserve(nnz, nrows, ncols);
624 builder.setNumCols(ncols);
628 for(
int i = 0; i < ncols; ++i)
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));
638 builder.setColIntegral(i,
false);
639 builder.setObj(i, objective * switch_sign);
643 builder.setNumRows(nrows);
645 for(
int i = 0; i < nrows; ++i)
648 int rowlength = rowVector.
size();
649 int* indices =
new int[rowlength];
650 R* rowValues =
new R[rowlength];
652 for(
int j = 0; j < rowlength; j++)
655 indices[j] = element.
idx;
656 rowValues[j] = element.
val;
659 builder.addRowEntries(i, rowlength, indices, rowValues);
663 builder.setRowLhs(i, lhs);
664 builder.setRowRhs(i, rhs);
665 builder.setRowLhsInf(i, lhs <= -R(
infinity));
666 builder.setRowRhsInf(i, rhs >= R(
infinity));
669 return builder.build();
679 m_keepbounds = keepbounds;
683 (*this->
spxout) <<
"==== PaPILO doesn't handle parameter keepbounds" <<
686 initLocalVariables(lp);
688 papilo::Problem<R> problem = buildProblem(lp);
689 papilo::Presolve<R> presolve;
691 configurePapilo(presolve, this->
tolerances()->floatingPointFeastol(), this->
tolerances()->epsilon(),
692 seed, remainingTime);
694 <<
" --- starting PaPILO" << std::endl;
697 papilo::PresolveResult<R> res = presolve.apply(problem);
699 assert(res.postsolve.postsolveType == PostsolveType::kFull);
703 case papilo::PresolveStatus::kInfeasible:
706 <<
" --- presolving detected infeasibility" << std::endl;
710 case papilo::PresolveStatus::kUnbndOrInfeas:
711 case papilo::PresolveStatus::kUnbounded:
714 "==== Presolving detected unboundedness of the problem" << std::endl;
718 case papilo::PresolveStatus::kUnchanged:
722 <<
"==== Presolving found nothing " << std::endl;
726 case papilo::PresolveStatus::kReduced:
731 int newNonzeros = problem.getConstraintMatrix().getNnz();
733 if(newNonzeros == 0 || ((problem.getNRows() <= modifyRowsFac * lp.
nRows() ||
734 newNonzeros <= modifyRowsFac * lp.
nNzos())))
737 <<
" --- presolved problem has " << problem.getNRows() <<
739 << problem.getNCols() <<
" cols and " 740 << newNonzeros <<
" non-zeros and " 741 << presolve.getStatistics().nboundchgs <<
" boundchanges and " 742 << presolve.getStatistics().nsidechgs <<
" sidechanges" 745 postsolveStorage = res.postsolve;
748 for(
int j = lp.
nCols() - 1; j >= 0; j--)
751 for(
int i = lp.
nRows() - 1; i >= 0; i--)
754 applyPresolveResultsToColumns(lp, problem, res);
755 applyPresolveResultsToRows(lp, problem, res);
756 assert(newNonzeros == lp.
nNzos());
764 <<
" --- presolve results smaller than the modifyconsfac" 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());
799 uint32_t seed,
Real remainingTime)
const 804 presolve.getPresolveOptions().randomseed = (
unsigned int) seed;
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;
817 presolve.setVerbosityLevel(verbosityLevel);
820 using uptr = std::unique_ptr<papilo::PresolveMethod<R>>;
823 if(enableSingletonCols)
824 presolve.addPresolveMethod(uptr(
new papilo::SingletonCols<R>()));
826 if(enablePropagation)
827 presolve.addPresolveMethod(uptr(
new papilo::ConstraintPropagation<R>()));
830 if(enableParallelRows)
831 presolve.addPresolveMethod(uptr(
new papilo::ParallelRowDetection<R>()));
833 if(enableParallelCols)
834 presolve.addPresolveMethod(uptr(
new papilo::ParallelColDetection<R>()));
836 if(enableSingletonStuffing)
837 presolve.addPresolveMethod(uptr(
new papilo::SingletonStuffing<R>()));
840 presolve.addPresolveMethod(uptr(
new papilo::DualFix<R>()));
842 if(enableFixContinuous)
843 presolve.addPresolveMethod(uptr(
new papilo::FixContinuous<R>()));
846 if(enableDominatedCols)
847 presolve.addPresolveMethod(uptr(
new papilo::DominatedCols<R>()));
866 const papilo::PresolveResult<R>& res)
const 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();
876 for(
int col = 0; col < problem.getNCols(); col++)
879 R lb = lowerBounds[col];
881 if(colFlags[col].test(papilo::ColFlag::kLbInf))
884 R ub = upperBounds[col];
886 if(colFlags[col].test(papilo::ColFlag::kUbInf))
889 LPColBase<R> column(objective.coefficients[col]* switch_sign, emptyVector, ub, lb);
891 assert(lp.
lower(col) == lb);
892 assert(lp.
upper(col) == ub);
897 assert(problem.getNCols() == lp.
nCols());
902 const papilo::PresolveResult<R>& res)
const 904 int size = res.postsolve.origrow_mapping.size();
907 for(
int row = 0; row < size; row++)
909 R rhs = problem.getConstraintMatrix().getRightHandSides()[row];
911 if(problem.getRowFlags()[row].test(papilo::RowFlag::kRhsInf))
914 R lhs = problem.getConstraintMatrix().getLeftHandSides()[row];
916 if(problem.getRowFlags()[row].test(papilo::RowFlag::kLhsInf))
919 const papilo::SparseVectorView<R> papiloRowVector =
920 problem.getConstraintMatrix().getRowCoefficients(row);
921 const int* indices = papiloRowVector.getIndices();
922 const R* values = papiloRowVector.getValues();
924 int length = papiloRowVector.getLength();
927 for(
int i = 0; i < length; i++)
929 soplexRowVector.
add(indices[i], values[i]);
934 assert(lp.
lhs(row) == lhs);
935 assert(lp.
rhs(row) == rhs);
938 assert(problem.getNRows() == lp.
nRows());
const VectorBase< R > & rhs() const
Returns right hand side vector.
virtual void removeRow(int i)
Removes i 'th row.
Safe arrays of data objects.Class DataArray provides safe arrays of Data Objects. For general C++ obj...
Sparse vector nonzero element.
#define SPX_MSG_WARNING(spxout, x)
Prints out message x if the verbosity level is at least SPxOut::WARNING.
const VectorBase< R > & upper() const
Returns upper bound vector.
Result
Result of the simplification.
virtual bool isUnsimplified() const
specifies whether an optimal solution has already been unsimplified.
Dense vector.Class VectorBase provides dense linear algebra vectors. Internally, VectorBase wraps std...
int size() const
Number of used indices.
Nonzero< R > & element(int n)
Reference to the n 'th nonzero element.
Sequential object-oriented SimPlex.SPxSolverBase is an LP solver class using the revised Simplex algo...
Exception classes for SoPlex.
Dynamic sparse vectors.Class DSVectorBase implements dynamic sparse vectors, i.e. SVectorBases with a...
Presol(Timer::TYPE ttype=Timer::USER_TIME)
#define SPX_MSG_ERROR(x)
Prints out message x if the verbosity level is at least SPxOut::ERROR.
virtual void removeCol(int i)
Removes i 'th column.
int idx
Index of nonzero element.
LP simplification base class.
virtual ~Presol()
destructor.
virtual const VectorBase< R > & unsimplifiedRedCost()
returns a reference to the unsimplified reduced costs.
int nRows() const
Returns number of rows in LP.
void add(const SVectorBase< S > &vec)
Append nonzeros of sv.
virtual void start()=0
start timer, resume accounting user, system and real time.
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)
simplification could be done
int nNzos() const
Returns number of nonzeros in LP.
Timer * m_timeUsed
user time used for simplification
SPxSimplifier & operator=(const SPxSimplifier &rhs)
assignment operator
SPxSense
Optimization sense.
const std::shared_ptr< Tolerances > tolerances() const
get the _tolerances member variable
SPxSense spxSense() const
Returns the optimization sense.
Presol(const Presol &old)
copy constructor.
LP simplification abstract base class.Instances of classes derived from SPxSimplifier may be loaded t...
int size() const
Returns number of elements.
R obj(int i) const
Returns objective value of column i.
const VectorBase< R > & lhs() const
Returns left hand side vector.
SPxSimplifier< R > * clone() const
virtual SPxSolverBase< R >::VarStatus getBasisColStatus(int j) const
gets basis status for a single column.
virtual void getBasis(typename SPxSolverBase< R >::VarStatus rows[], typename SPxSolverBase< R >::VarStatus cols[], const int rowsSize, const int colsSize) const
get optimal basis.
SOPLEX_THREADLOCAL const Real infinity
(In)equality for LPs.Class LPRowBase provides constraints for linear programs in the form where a is...
Debugging, floating point type and parameter definitions.
virtual SPxSimplifier< R >::Result simplify(SPxLPBase< R > &lp, Real remainingTime, bool keepbounds, uint32_t seed)
virtual void addCol(const LPColBase< R > &col, bool scale=false)
int dim() const
Dimension of vector.
virtual const VectorBase< R > & unsimplifiedSlacks()
returns a reference to the unsimplified slack values.
virtual const VectorBase< R > & unsimplifiedDual()
returns a reference to the unsimplified dual solution.
Everything should be within this namespace.
Saving LPs in a form suitable for SoPlex.Class SPxLPBase provides the data structures required for sa...
Presol & operator=(const Presol &rhs)
assignment operator
#define SPX_MSG_INFO1(spxout, x)
Prints out message x if the verbosity level is at least SPxOut::INFO1.
Save arrays of arbitrary types.
void changeObjOffset(const T &o)
const SVectorBase< R > & rowVector(int i) const
Gets row vector of row i.
virtual void addRow(const LPRowBase< R > &row, bool scale=false)
int size() const
return nr. of elements.
R val
Value of nonzero element.
int nCols() const
Returns number of columns in LP.
virtual SPxSimplifier< R >::Result result() const
returns result status of the simplification
Sparse vectors.Class SVectorBase provides packed sparse vectors. Such are a sparse vectors...
virtual SPxSolverBase< R >::VarStatus getBasisRowStatus(int i) const
gets basis status for a single row.
virtual void reset()=0
initialize timer, set timing accounts to zero.
SPxOut * spxout
message handler
virtual const VectorBase< R > & unsimplifiedPrimal()
returns a reference to the unsimplified primal solution.
const VectorBase< R > & lower() const
Returns (internal and possibly scaled) lower bound vector.
LP column.Class LPColBase provides a datatype for storing the column of an LP a the form similar to ...