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   	/**@file  spxleastsqsc.hpp
26   	 * @brief LP least squares scaling.
27   	 */
28   	
29   	#include <assert.h>
30   	#include <cmath>
31   	#include "soplex/spxout.h"
32   	#include "soplex/basevectors.h"
33   	#include "soplex/svsetbase.h"
34   	#include "soplex/svectorbase.h"
35   	#include "soplex/ssvectorbase.h"
36   	#include <array>
37   	
38   	namespace soplex
39   	{
40   	
41   	/* update scaling VectorBase<R> */
42   	template <class R>
43   	static void updateScale(
44   	   const SSVectorBase<R> vecnnzeroes,
45   	   const SSVectorBase<R> resnvec,
46   	   SSVectorBase<R>& tmpvec,
47   	   SSVectorBase<R>*& psccurr,
48   	   SSVectorBase<R>*& pscprev,
49   	   R qcurr,
50   	   R qprev,
51   	   R eprev1,
52   	   R eprev2,
53   	   R epsilon)
54   	{
55   	   assert(psccurr != NULL);
56   	   assert(pscprev != NULL);
57   	   assert(qcurr * qprev != 0.0);
58   	
59   	   R fac = -(eprev1 * eprev2);
60   	
61   	   SSVectorBase<R>* pssv;
62   	
63   	   *pscprev -= *psccurr;
64   	
65   	   if(isZero(fac, epsilon))
66   	      (*pscprev).clear();
67   	   else
68   	      *pscprev *= fac;
69   	
70   	   *pscprev += tmpvec.assignPWproduct4setup(resnvec, vecnnzeroes);
71   	
72   	   *pscprev *= 1.0 / (qcurr * qprev);
73   	   *pscprev += *psccurr;
74   	
75   	   /* swap pointers */
76   	   pssv = psccurr;
77   	   psccurr = pscprev;
78   	   pscprev = pssv;
79   	}
80   	
81   	/* update scaling VectorBase<R> after main loop */
82   	template <class R>
83   	static void updateScaleFinal(
84   	   const SSVectorBase<R> vecnnzeroes,
85   	   const SSVectorBase<R> resnvec,
86   	   SSVectorBase<R>& tmpvec,
87   	   SSVectorBase<R>*& psccurr,
88   	   SSVectorBase<R>*& pscprev,
89   	   R q,
90   	   R eprev1,
91   	   R eprev2,
92   	   R epsilon)
93   	{
94   	   assert(q != 0);
95   	   assert(psccurr != NULL);
96   	   assert(pscprev != NULL);
97   	
98   	   R fac = -(eprev1 * eprev2);
99   	
100  	   *pscprev -= *psccurr;
101  	
102  	   if(isZero(fac, epsilon))
103  	      (*pscprev).clear();
104  	   else
105  	      *pscprev *= fac;
106  	
107  	   *pscprev += tmpvec.assignPWproduct4setup(resnvec, vecnnzeroes);
108  	   *pscprev *= 1.0 / q;
109  	   *pscprev += *psccurr;
110  	
111  	   psccurr = pscprev;
112  	}
113  	
114  	/* update residual VectorBase<R> */
115  	template <class R>
116  	static inline void updateRes(
117  	   const SVSetBase<R> facset,
118  	   const SSVectorBase<R> resvecprev,
119  	   SSVectorBase<R>& resvec,
120  	   SSVectorBase<R>& tmpvec,
121  	   R eprev,
122  	   R qcurr,
123  	   R epsilon)
124  	{
125  	   assert(qcurr != 0.0);
126  	
127  	   if(isZero(eprev, epsilon))
128  	      resvec.clear();
129  	   else
130  	      resvec *= eprev;
131  	
132  	   int dummy1 = 0;
133  	   int dummy2 = 0;
134  	   tmpvec.assign2product4setup(facset, resvecprev, 0, 0, dummy1, dummy2);
135  	   tmpvec.setup();
136  	   resvec += tmpvec;
137  	
138  	   resvec *= (-1.0 / qcurr);
139  	   resvec.setup();
140  	}
141  	
142  	
143  	/* initialize constant vectors and matrices */
144  	template <class R>
145  	static void initConstVecs(
146  	   const SVSetBase<R>* vecset,
147  	   SVSetBase<R>& facset,
148  	   SSVectorBase<R>& veclogs,
149  	   SSVectorBase<R>& vecnnzinv,
150  	   R epsilon)
151  	{
152  	   assert(vecset != NULL);
153  	
154  	   const int nvec = vecset->num();
155  	
156  	   for(int k = 0; k < nvec; ++k)
157  	   {
158  	      R logsum = 0.0;
159  	      int nnz = 0;
160  	      // get kth row or column of LP
161  	      const SVectorBase<R>& lpvec = (*vecset)[k];
162  	      const int size = lpvec.size();
163  	
164  	      for(int i = 0; i < size; ++i)
165  	      {
166  	         const R a = lpvec.value(i);
167  	
168  	         if(!isZero(a, epsilon))
169  	         {
170  	            logsum += log2(double(spxAbs(a))); // todo spxLog2?
171  	            nnz++;
172  	         }
173  	      }
174  	
175  	      R nnzinv;
176  	
177  	      if(nnz > 0)
178  	      {
179  	         nnzinv = 1.0 / nnz;
180  	      }
181  	      else
182  	      {
183  	         /* all-0 entries */
184  	         logsum = 1.0;
185  	         nnzinv = 1.0;
186  	      }
187  	
188  	      veclogs.add(k, logsum);
189  	      vecnnzinv.add(k, nnzinv);
190  	
191  	      /* create new VectorBase<R> for facset */
192  	      SVectorBase<R>& vecnew = (*(facset.create(nnz)));
193  	
194  	      for(int i = 0; i < size; ++i)
195  	      {
196  	         if(!isZero(lpvec.value(i), epsilon))
197  	            vecnew.add(lpvec.index(i), nnzinv);
198  	      }
199  	
200  	      vecnew.sort();
201  	   }
202  	
203  	   assert(veclogs.isSetup());
204  	   assert(vecnnzinv.isSetup());
205  	}
206  	
207  	/* return name of scaler */
208  	static inline const char* makename()
209  	{
210  	   return "Least squares";
211  	}
212  	
213  	template <class R>
214  	SPxLeastSqSC<R>::SPxLeastSqSC()
215  	   : SPxScaler<R>(makename(), false, false)
216  	{}
217  	
218  	template <class R>
219  	SPxLeastSqSC<R>::SPxLeastSqSC(const SPxLeastSqSC<R>& old)
220  	   : SPxScaler<R>(old), acrcydivisor(old.acrcydivisor), maxrounds(old.maxrounds)
221  	{}
222  	
223  	template <class R>
224  	SPxLeastSqSC<R>& SPxLeastSqSC<R>::operator=(const SPxLeastSqSC<R>& rhs)
225  	{
226  	   if(this != &rhs)
227  	   {
228  	      SPxScaler<R>::operator=(rhs);
229  	   }
230  	
231  	   return *this;
232  	}
233  	
234  	template <class R>
235  	void SPxLeastSqSC<R>::setRealParam(R param, const char* name)
236  	{
237  	   assert(param >= 1.0);
238  	   acrcydivisor = param;
239  	}
240  	
241  	template <class R>
242  	void SPxLeastSqSC<R>::setIntParam(int param, const char* name)
243  	{
244  	   assert(param >= 0);
245  	   maxrounds = param;
246  	}
247  	
248  	// todo refactor this method. Has no abstraction and is too long
249  	template <class R>
250  	void SPxLeastSqSC<R>::scale(SPxLPBase<R>& lp,  bool persistent)
251  	{
252  	   SPX_MSG_INFO1((*this->spxout), (*this->spxout) << "Least squares LP scaling" <<
253  	                 (persistent ? " (persistent)" : "") << std::endl;)
254  	
255  	   this->setup(lp);
256  	
257  	   const int nrows = lp.nRows();
258  	   const int ncols = lp.nCols();
259  	   const int lpnnz = lp.nNzos();
260  	
261  	   // is contraints matrix empty?
262  	   //todo don't create the scaler in this case!
263  	   if(lpnnz == 0)
264  	   {
265  	      // to keep the invariants, we still need to call this method
266  	      this->applyScaling(lp);
267  	
268  	      return;
269  	   }
270  	
271  	   assert(nrows > 0 && ncols > 0 && lpnnz > 0);
272  	
273  	   /* constant factor matrices;
274  	    * in Curtis-Reid article
275  	    * facnrows equals E^T M^(-1)
276  	    * facncols equals E N^(-1)
277  	    * */
278  	   SVSetBase<R> facnrows(nrows, nrows, 1.1, 1.2);
279  	   SVSetBase<R> facncols(ncols, ncols, 1.1, 1.2);
280  	
281  	   /* column scaling factor vectors */
282  	   SSVectorBase<R> colscale1(ncols, this->_tolerances);
283  	   SSVectorBase<R> colscale2(ncols, this->_tolerances);
284  	
285  	   /* row scaling factor vectors */
286  	   SSVectorBase<R> rowscale1(nrows, this->_tolerances);
287  	   SSVectorBase<R> rowscale2(nrows, this->_tolerances);
288  	
289  	   /* residual vectors */
290  	   SSVectorBase<R> resnrows(nrows, this->_tolerances);
291  	   SSVectorBase<R> resncols(ncols, this->_tolerances);
292  	
293  	   /* vectors to store temporary values */
294  	   SSVectorBase<R> tmprows(nrows, this->_tolerances);
295  	   SSVectorBase<R> tmpcols(ncols, this->_tolerances);
296  	
297  	   /* vectors storing the row and column sums (respectively) of logarithms of
298  	    *(absolute values of) non-zero elements of left hand matrix of LP
299  	    */
300  	   SSVectorBase<R> rowlogs(nrows, this->_tolerances);
301  	   SSVectorBase<R> collogs(ncols, this->_tolerances);
302  	
303  	   /* vectors storing the inverted number of non-zeros in each row and column
304  	    *(respectively) of left hand matrix of LP
305  	    */
306  	   SSVectorBase<R> rownnzinv(nrows, this->_tolerances);
307  	   SSVectorBase<R> colnnzinv(ncols, this->_tolerances);
308  	
309  	   /* VectorBase<R> pointers */
310  	   SSVectorBase<R>* csccurr = &colscale1;
311  	   SSVectorBase<R>* cscprev = &colscale2;
312  	   SSVectorBase<R>* rsccurr = &rowscale1;
313  	   SSVectorBase<R>* rscprev = &rowscale2;
314  	
315  	   SPX_MSG_INFO2((*this->spxout), (*this->spxout) << "before scaling:"
316  	                 << " min= " << lp.minAbsNzo()
317  	                 << " max= " << lp.maxAbsNzo()
318  	                 << " col-ratio= " << this->maxColRatio(lp)
319  	                 << " row-ratio= " << this->maxRowRatio(lp)
320  	                 << std::endl;)
321  	
322  	   /* initialize scalars, vectors and matrices */
323  	
324  	   assert(acrcydivisor > 0.0);
325  	
326  	   const R smax = lpnnz / acrcydivisor;
327  	   R qcurr = 1.0;
328  	   R qprev = 0.0;
329  	
330  	   std::array<R, 3> eprev;
331  	   eprev.fill(0.0);
332  	
333  	   initConstVecs(lp.rowSet(), facnrows, rowlogs, rownnzinv, R(this->tolerances()->epsilon()));
334  	   initConstVecs(lp.colSet(), facncols, collogs, colnnzinv, R(this->tolerances()->epsilon()));
335  	
336  	   assert(tmprows.isSetup() && tmpcols.isSetup());
337  	   assert(rowscale1.isSetup() && rowscale2.isSetup());
338  	   assert(colscale1.isSetup() && colscale2.isSetup());
339  	
340  	   // compute first residual vector r0
341  	   int dummy1 = 0;
342  	   int dummy2 = 0;
343  	   resncols = collogs - tmpcols.assign2product4setup(facnrows, rowlogs, 0, 0, dummy1, dummy2);
344  	
345  	   resncols.setup();
346  	   resnrows.setup();
347  	
348  	   rowscale1.assignPWproduct4setup(rownnzinv, rowlogs);
349  	   rowscale2 = rowscale1;
350  	
351  	   R scurr = resncols * tmpcols.assignPWproduct4setup(colnnzinv, resncols);
352  	
353  	   int k;
354  	
355  	   /* conjugate gradient loop */
356  	   for(k = 0; k < maxrounds; ++k)
357  	   {
358  	      const R sprev = scurr;
359  	
360  	      // termination criterion met?
361  	      if(scurr < smax)
362  	         break;
363  	
364  	      // is k even?
365  	      if((k % 2) == 0)
366  	      {
367  	         // not in first iteration?
368  	         if(k != 0)   // true, then update row scaling factor vector
369  	            updateScale(rownnzinv, resnrows, tmprows, rsccurr, rscprev, qcurr, qprev, eprev[1], eprev[2],
370  	                        R(this->tolerances()->epsilon()));
371  	
372  	         updateRes(facncols, resncols, resnrows, tmprows, eprev[0], qcurr, R(this->tolerances()->epsilon()));
373  	         scurr = resnrows * tmprows.assignPWproduct4setup(resnrows, rownnzinv);
374  	      }
375  	      else // k is odd
376  	      {
377  	         // update column scaling factor vector
378  	         updateScale(colnnzinv, resncols, tmpcols, csccurr, cscprev, qcurr, qprev, eprev[1], eprev[2],
379  	                     R(this->tolerances()->epsilon()));
380  	
381  	         updateRes(facnrows, resnrows, resncols, tmpcols, eprev[0], qcurr, R(this->tolerances()->epsilon()));
382  	         scurr = resncols * tmpcols.assignPWproduct4setup(resncols, colnnzinv);
383  	      }
384  	
385  	      // shift eprev entries one to the right
386  	      for(unsigned l = 2; l > 0; --l)
387  	         eprev[l] = eprev[l - 1];
388  	
389  	      assert(isNotZero(sprev, R(this->tolerances()->epsilon())));
390  	
391  	      eprev[0] = (qcurr * scurr) / sprev;
392  	
393  	      const R tmp = qcurr;
394  	      qcurr = 1.0 - eprev[0];
395  	      qprev = tmp;
396  	   }
397  	
398  	   if(k > 0 && (k % 2) == 0)
399  	   {
400  	      // update column scaling factor vector
401  	      updateScaleFinal(colnnzinv, resncols, tmpcols, csccurr, cscprev, qprev, eprev[1], eprev[2],
402  	                       R(this->tolerances()->epsilon()));
403  	   }
404  	   else if(k > 0)
405  	   {
406  	      // update row scaling factor vector
407  	      updateScaleFinal(rownnzinv, resnrows, tmprows, rsccurr, rscprev, qprev, eprev[1], eprev[2],
408  	                       R(this->tolerances()->epsilon()));
409  	   }
410  	
411  	   /* compute actual scaling factors */
412  	
413  	   const SSVectorBase<R>& rowscale = *rsccurr;
414  	   const SSVectorBase<R>& colscale = *csccurr;
415  	
416  	   DataArray<int>& colscaleExp = *this->m_activeColscaleExp;
417  	   DataArray<int>& rowscaleExp = *this->m_activeRowscaleExp;
418  	
419  	   for(k = 0; k < nrows; ++k)
420  	      rowscaleExp[k] = -int(rowscale[k] + ((rowscale[k] >= 0.0) ? (+0.5) : (-0.5)));
421  	
422  	   for(k = 0; k < ncols; ++k)
423  	      colscaleExp[k] = -int(colscale[k] + ((colscale[k] >= 0.0) ? (+0.5) : (-0.5)));
424  	
425  	   // scale
426  	   this->applyScaling(lp);
427  	
428  	   SPX_MSG_INFO3((*this->spxout), (*this->spxout) << "Row scaling min= " << this->minAbsRowscale()
429  	                 << " max= " << this->maxAbsRowscale()
430  	                 << std::endl
431  	                 << "Col scaling min= " << this->minAbsColscale()
432  	                 << " max= " << this->maxAbsColscale()
433  	                 << std::endl;)
434  	
435  	   SPX_MSG_INFO2((*this->spxout), (*this->spxout) << "after scaling: "
436  	                 << " min= " << lp.minAbsNzo(false)
437  	                 << " max= " << lp.maxAbsNzo(false)
438  	                 << " col-ratio= " << this->maxColRatio(lp)
439  	                 << " row-ratio= " << this->maxRowRatio(lp)
440  	                 << std::endl;)
441  	}
442  	
443  	} // namespace soplex
444