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