Scippy

SoPlex

Sequential object-oriented simPlex

spxpapilo.h
Go to the documentation of this file.
1/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2/* */
3/* This file is part of the class library */
4/* SoPlex --- the Sequential object-oriented simPlex. */
5/* */
6/* Copyright (c) 1996-2024 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
28namespace soplex
29{
30
31template <class R> class Presol : public SPxSimplifier<R>
32{
33
34public:
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
169namespace soplex
170{
171
172template<class R>
173class Presol : public SPxSimplifier<R>
174{
175private:
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
214protected:
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
231public:
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
421private:
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
441template <class R>
442void 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
555template <class R>
556papilo::VarBasisStatus
557Presol<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
583template <class R>
584typename SPxSolverBase<R>::VarStatus
585Presol<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
612template<class R>
613papilo::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
673template<class R>
674typename SPxSimplifier<R>::Result
675Presol<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
778template<class R>
779void 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
797template<class R>
798void 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
864template<class R>
865void 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
900template<class R>
901void 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
array.h
Save arrays of arbitrary types.
soplex::Presol::unsimplify
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)
Definition: spxpapilo.h:70
soplex::Presol::Presol
Presol(Timer::TYPE ttype=Timer::USER_TIME)
Definition: spxpapilo.h:40
soplex::Presol::operator=
Presol & operator=(const Presol &rhs)
assignment operator
Definition: spxpapilo.h:47
soplex::Presol::unsimplifiedSlacks
virtual const VectorBase< R > & unsimplifiedSlacks()
returns a reference to the unsimplified slack values.
Definition: spxpapilo.h:110
soplex::Presol::unsimplifiedRedCost
virtual const VectorBase< R > & unsimplifiedRedCost()
returns a reference to the unsimplified reduced costs.
Definition: spxpapilo.h:118
soplex::Presol::clone
SPxSimplifier< R > * clone() const
Definition: spxpapilo.h:58
soplex::Presol::getBasis
virtual void getBasis(typename SPxSolverBase< R >::VarStatus rows[], typename SPxSolverBase< R >::VarStatus cols[], const int rowsSize, const int colsSize) const
get optimal basis.
Definition: spxpapilo.h:140
soplex::Presol::~Presol
virtual ~Presol()
destructor.
Definition: spxpapilo.h:53
soplex::Presol::getBasisColStatus
virtual SPxSolverBase< R >::VarStatus getBasisColStatus(int j) const
gets basis status for a single column.
Definition: spxpapilo.h:133
soplex::Presol::result
virtual SPxSimplifier< R >::Result result() const
returns result status of the simplification
Definition: spxpapilo.h:80
soplex::Presol::unsimplifiedDual
virtual const VectorBase< R > & unsimplifiedDual()
returns a reference to the unsimplified dual solution.
Definition: spxpapilo.h:102
soplex::Presol::getBasisRowStatus
virtual SPxSolverBase< R >::VarStatus getBasisRowStatus(int i) const
gets basis status for a single row.
Definition: spxpapilo.h:126
soplex::Presol::Presol
Presol(const Presol &old)
copy constructor.
Definition: spxpapilo.h:44
soplex::Presol::simplify
virtual SPxSimplifier< R >::Result simplify(SPxLPBase< R > &lp, Real remainingTime, bool keepbounds, uint32_t seed)
Definition: spxpapilo.h:63
soplex::Presol::unsimplifiedPrimal
virtual const VectorBase< R > & unsimplifiedPrimal()
returns a reference to the unsimplified primal solution.
Definition: spxpapilo.h:94
soplex::Presol::isUnsimplified
virtual bool isUnsimplified() const
specifies whether an optimal solution has already been unsimplified.
Definition: spxpapilo.h:87
soplex::SPxLPBase
Saving LPs in a form suitable for SoPlex.
Definition: spxlpbase.h:108
soplex::SPxLPBase::SPxSense
SPxSense
Optimization sense.
Definition: spxlpbase.h:125
soplex::SPxSimplifier
LP simplification abstract base class.
Definition: spxsimplifier.h:52
soplex::SPxSimplifier::Result
Result
Result of the simplification.
Definition: spxsimplifier.h:93
soplex::SPxSimplifier::INFEASIBLE
@ INFEASIBLE
primal infeasibility was detected
Definition: spxsimplifier.h:95
soplex::SPxSimplifier::UNBOUNDED
@ UNBOUNDED
primal unboundedness was detected
Definition: spxsimplifier.h:97
soplex::SPxSimplifier::OKAY
@ OKAY
simplification could be done
Definition: spxsimplifier.h:94
soplex::SPxSimplifier::VANISHED
@ VANISHED
the problem was so much simplified that it vanished
Definition: spxsimplifier.h:98
soplex::SPxSimplifier::tolerances
const std::shared_ptr< Tolerances > tolerances() const
get the _tolerances member variable
Definition: spxsimplifier.h:273
soplex::SPxSimplifier::operator=
SPxSimplifier & operator=(const SPxSimplifier &rhs)
assignment operator
Definition: spxsimplifier.h:144
soplex::SPxSimplifier::SPxSimplifier
SPxSimplifier(const char *p_name, Timer::TYPE ttype=Timer::USER_TIME)
constructor
Definition: spxsimplifier.h:106
soplex::SPxSolverBase
Sequential object-oriented SimPlex.
Definition: spxsolver.h:104
soplex::SPxSolverBase::BASIC
@ BASIC
variable is basic.
Definition: spxsolver.h:201
soplex::SPxSolverBase::ON_LOWER
@ ON_LOWER
variable set to its lower bound.
Definition: spxsolver.h:198
soplex::SPxSolverBase::ON_UPPER
@ ON_UPPER
variable set to its upper bound.
Definition: spxsolver.h:197
soplex::SPxSolverBase::UNDEFINED
@ UNDEFINED
nothing known about basis status (possibly due to a singular basis in transformed problem)
Definition: spxsolver.h:202
soplex::SPxSolverBase::FIXED
@ FIXED
variable fixed to identical bounds.
Definition: spxsolver.h:199
soplex::SPxSolverBase::ZERO
@ ZERO
free variable fixed to zero.
Definition: spxsolver.h:200
soplex::Timer::TYPE
TYPE
types of timers
Definition: timer.h:109
soplex::VectorBase
Dense vector.
Definition: vectorbase.h:86
exceptions.h
Exception classes for SoPlex.
soplex
Everything should be within this namespace.
soplex::Real
double Real
Definition: spxdefines.h:269
soplex::infinity
SOPLEX_THREADLOCAL const Real infinity
Definition: spxdefines.cpp:41
spxdefines.h
Debugging, floating point type and parameter definitions.
SPX_MSG_ERROR
#define SPX_MSG_ERROR(x)
Prints out message x if the verbosity level is at least SPxOut::ERROR.
Definition: spxdefines.h:163
SPX_MSG_WARNING
#define SPX_MSG_WARNING(spxout, x)
Prints out message x if the verbosity level is at least SPxOut::WARNING.
Definition: spxdefines.h:165
SPX_MSG_INFO1
#define SPX_MSG_INFO1(spxout, x)
Prints out message x if the verbosity level is at least SPxOut::INFO1.
Definition: spxdefines.h:167
spxsimplifier.h
LP simplification base class.