Scippy

SoPlex

Sequential object-oriented simPlex

spxleastsqsc.cpp
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-2017 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SoPlex is distributed under the terms of the ZIB Academic Licence. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SoPlex; see the file COPYING. If not email to soplex@zib.de. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15 
16 /**@file spxleastsqsc.cpp
17  * @brief LP least squares scaling.
18  */
19 #include <cmath>
20 #include <assert.h>
21 #include "spxleastsqsc.h"
22 #include "spxout.h"
23 #include "basevectors.h"
24 #include "svsetbase.h"
25 #include "svectorbase.h"
26 #include "ssvectorbase.h"
27 #include <array>
28 
29 namespace soplex
30 {
31 
32 /* update scaling vector */
33 static void updateScale(
34  const SSVector vecnnzeroes,
35  const SSVector resnvec,
36  SSVector& tmpvec,
37  SSVector*& psccurr,
38  SSVector*& pscprev,
39  Real qcurr,
40  Real qprev,
41  Real eprev1,
42  Real eprev2)
43 {
44  assert(psccurr != NULL);
45  assert(pscprev != NULL);
46  assert(qcurr * qprev != 0.0);
47 
48  Real fac = -(eprev1 * eprev2);
49 
50  SSVector* pssv;
51 
52  *pscprev -= *psccurr;
53 
54  if( isZero(fac) )
55  (*pscprev).clear();
56  else
57  *pscprev *= fac;
58 
59  *pscprev += tmpvec.assignPWproduct4setup(resnvec, vecnnzeroes);
60 
61  *pscprev *= 1.0 / (qcurr * qprev);
62  *pscprev += *psccurr;
63 
64  /* swap pointers */
65  pssv = psccurr;
66  psccurr = pscprev;
67  pscprev = pssv;
68 }
69 
70 /* update scaling vector after main loop */
71 static void updateScaleFinal(
72  const SSVector vecnnzeroes,
73  const SSVector resnvec,
74  SSVector& tmpvec,
75  SSVector*& psccurr,
76  SSVector*& pscprev,
77  Real q,
78  Real eprev1,
79  Real eprev2)
80 {
81  assert(q != 0);
82  assert(psccurr != NULL);
83  assert(pscprev != NULL);
84 
85  Real fac = -(eprev1 * eprev2);
86 
87  *pscprev -= *psccurr;
88 
89  if( isZero(fac) )
90  (*pscprev).clear();
91  else
92  *pscprev *= fac;
93 
94  *pscprev += tmpvec.assignPWproduct4setup(resnvec, vecnnzeroes);
95  *pscprev *= 1.0 / q;
96  *pscprev += *psccurr;
97 
98  psccurr = pscprev;
99 }
100 
101 /* update residual vector */
102 static inline void updateRes(
103  const SVSet facset,
104  const SSVector resvecprev,
105  SSVector& resvec,
106  SSVector& tmpvec,
107  Real eprev,
108  Real qcurr)
109 {
110  assert(qcurr != 0.0);
111 
112  if( isZero(eprev) )
113  resvec.clear();
114  else
115  resvec *= eprev;
116 
117  tmpvec.assign2product4setup(facset, resvecprev);
118  tmpvec.setup();
119  resvec += tmpvec;
120 
121  resvec *= (-1.0 / qcurr);
122  resvec.setup();
123 }
124 
125 
126 /* initialize constant vectors and matrices */
127 static void initConstVecs(
128  const SVSet* vecset,
129  SVSet& facset,
130  SSVector& veclogs,
131  SSVector& vecnnzinv)
132 {
133  assert(vecset != NULL);
134 
135  const int nvec = vecset->num();
136 
137  for( int k = 0; k < nvec; ++k )
138  {
139  Real logsum = 0.0;
140  int nnz = 0;
141  // get kth row or column of LP
142  const SVector& lpvec = (*vecset)[k];
143  const int size = lpvec.size();
144 
145  for( int i = 0; i < size; ++i )
146  {
147  const Real a = lpvec.value(i);
148 
149  if( !isZero(a) )
150  {
151  logsum += log2(double(spxAbs(a))); // todo spxLog2?
152  nnz++;
153  }
154  }
155 
156  Real nnzinv;
157  if( nnz > 0)
158  {
159  nnzinv = 1.0 / nnz;
160  }
161  else
162  {
163  /* all-0 entries */
164  logsum = 1.0;
165  nnzinv = 1.0;
166  }
167 
168  veclogs.add(k, logsum);
169  vecnnzinv.add(k, nnzinv);
170 
171  /* create new vector for facset */
172  SVector& vecnew = (*(facset.create(nnz)));
173 
174  for( int i = 0; i < size; ++i )
175  {
176  if( !isZero(lpvec.value(i)) )
177  vecnew.add(lpvec.index(i), nnzinv);
178  }
179  vecnew.sort();
180  }
181 
182  assert(veclogs.isSetup());
183  assert(vecnnzinv.isSetup());
184 }
185 
186 /* return name of scaler */
187 static const char* makename()
188 {
189  return "Least squares";
190 }
191 
193  : SPxScaler(makename(), false, false)
194 {}
195 
198 {}
199 
201 {
202  if(this != &rhs)
203  {
205  }
206 
207  return *this;
208 }
209 
210 
211 void SPxLeastSqSC::setRealParam(Real param, const char* name)
212 {
213  assert(param >= 1.0);
214  acrcydivisor = param;
215 }
216 
217 void SPxLeastSqSC::setIntParam(int param, const char* name)
218 {
219  assert(param >= 0);
220  maxrounds = param;
221 }
222 
223 void SPxLeastSqSC::scale(SPxLP& lp, bool persistent)
224 {
225  MSG_INFO1( (*spxout), (*spxout) << "Least squares LP scaling" << (persistent ? " (persistent)" : "") << std::endl; )
226 
227  setup(lp);
228 
229  const int nrows = lp.nRows();
230  const int ncols = lp.nCols();
231  const int lpnnz = lp.nNzos();
232 
233  /* constant factor matrices;
234  * in Curtis-Reid article
235  * facnrows equals E^T M^(-1)
236  * facncols equals E N^(-1)
237  * */
238  SVSet facnrows(nrows, nrows, 1.1, 1.2);
239  SVSet facncols(ncols, ncols, 1.1, 1.2);
240 
241  /* column scaling factor vectors */
242  SSVector colscale1(ncols);
243  SSVector colscale2(ncols);
244 
245  /* row scaling factor vectors */
246  SSVector rowscale1(nrows);
247  SSVector rowscale2(nrows);
248 
249  /* residual vectors */
250  SSVector resnrows(nrows);
251  SSVector resncols(ncols);
252 
253  /* vectors to store temporary values */
254  SSVector tmprows(nrows);
255  SSVector tmpcols(ncols);
256 
257  /* vectors storing the row and column sums (respectively) of logarithms of
258  *(absolute values of) non-zero elements of left hand matrix of LP
259  */
260  SSVector rowlogs(nrows);
261  SSVector collogs(ncols);
262 
263  /* vectors storing the inverted number of non-zeros in each row and column
264  *(respectively) of left hand matrix of LP
265  */
266  SSVector rownnzinv(nrows);
267  SSVector colnnzinv(ncols);
268 
269  /* vector pointers */
270  SSVector* csccurr = &colscale1;
271  SSVector* cscprev = &colscale2;
272  SSVector* rsccurr = &rowscale1;
273  SSVector* rscprev = &rowscale2;
274 
275  MSG_INFO2( (*spxout), (*spxout) << "before scaling:"
276  << " min= " << lp.minAbsNzo()
277  << " max= " << lp.maxAbsNzo()
278  << " col-ratio= " << maxColRatio(lp)
279  << " row-ratio= " << maxRowRatio(lp)
280  << std::endl; )
281 
282  /* initialize scalars, vectors and matrices */
283 
284  assert(acrcydivisor > 0.0);
285 
286  const Real smax = lpnnz / acrcydivisor;
287  Real qcurr = 1.0;
288  Real qprev = 0.0;
289 
290  std::array<Real, 3> eprev;
291  eprev.fill(0.0);
292 
293  initConstVecs(lp.rowSet(), facnrows, rowlogs, rownnzinv);
294  initConstVecs(lp.colSet(), facncols, collogs, colnnzinv);
295 
296  assert(tmprows.isSetup());
297  assert(tmpcols.isSetup());
298  assert(rowscale1.isSetup());
299  assert(rowscale2.isSetup());
300  assert(colscale1.isSetup());
301  assert(colscale2.isSetup());
302 
303  // compute first residual vector r0
304  resncols = collogs - tmpcols.assign2product4setup(facnrows, rowlogs);
305 
306  resncols.setup();
307  resnrows.setup();
308 
309  rowscale1.assignPWproduct4setup(rownnzinv, rowlogs);
310  rowscale2 = rowscale1;
311 
312  Real scurr = resncols * tmpcols.assignPWproduct4setup(colnnzinv, resncols);
313 
314  int k;
315 
316  /* conjugate gradient loop */
317  for( k = 0; k < maxrounds; ++k )
318  {
319  const Real sprev = scurr;
320 
321  // is k even?
322  if( (k % 2) == 0 )
323  {
324  // not in first iteration?
325  if( k != 0 ) // true, then update row scaling factor vector
326  updateScale(rownnzinv, resnrows, tmprows, rsccurr, rscprev, qcurr, qprev, eprev[1], eprev[2]);
327 
328  updateRes(facncols, resncols, resnrows, tmprows, eprev[0], qcurr);
329  scurr = resnrows * tmprows.assignPWproduct4setup(resnrows, rownnzinv);
330  }
331  else // k is odd
332  {
333  // update column scaling factor vector
334  updateScale(colnnzinv, resncols, tmpcols, csccurr, cscprev, qcurr, qprev, eprev[1], eprev[2]);
335 
336  updateRes(facnrows, resnrows, resncols, tmpcols, eprev[0], qcurr);
337  scurr = resncols * tmpcols.assignPWproduct4setup(resncols, colnnzinv);
338  }
339 
340  // shift eprev entries one to the right
341  for( unsigned l = 2; l > 0; --l)
342  eprev[l] = eprev[l - 1];
343 
344  eprev[0] = (qcurr * scurr) / sprev;
345 
346  const Real tmp = qcurr;
347  qcurr = 1.0 - eprev[0];
348  qprev = tmp;
349 
350  // termination criterion met?
351  if( scurr < smax )
352  break;
353  }
354 
355  // is k even?
356  if( (k % 2) == 0 )
357  {
358  // update column scaling factor vector
359  updateScaleFinal(colnnzinv, resncols, tmpcols, csccurr, cscprev, qprev, eprev[1], eprev[2]);
360  }
361  else // k is odd
362  {
363  // update row scaling factor vector
364  updateScaleFinal(rownnzinv, resnrows, tmprows, rsccurr, rscprev, qprev, eprev[1], eprev[2]);
365  }
366 
367  /* compute actual scaling factors */
368 
369  const SSVector& rowscale = *rsccurr;
370  const SSVector& colscale = *csccurr;
371 
372  DataArray<int>& colscaleExp = *m_activeColscaleExp;
373  DataArray<int>& rowscaleExp = *m_activeRowscaleExp;
374 
375  for( k = 0; k < nrows; ++k )
376  rowscaleExp[k] = -int( rowscale[k] + ((rowscale[k] >= 0.0)? (+0.5) : (-0.5)) );
377 
378  for( k = 0; k < ncols; ++k )
379  colscaleExp[k] = -int( colscale[k] + ((colscale[k] >= 0.0)? (+0.5) : (-0.5)) );
380 
381  // scale
382  applyScaling(lp);
383 
384  MSG_INFO3( (*spxout), (*spxout) << "Row scaling min= " << minAbsRowscale()
385  << " max= " << maxAbsRowscale()
386  << std::endl
387  << "Col scaling min= " << minAbsColscale()
388  << " max= " << maxAbsColscale()
389  << std::endl; )
390 
391  MSG_INFO2( (*spxout), (*spxout) << "after scaling: "
392  << " min= " << lp.minAbsNzo(false)
393  << " max= " << lp.maxAbsNzo(false)
394  << " col-ratio= " << maxColRatio(lp)
395  << " row-ratio= " << maxRowRatio(lp)
396  << std::endl; )
397 }
398 
399 } // namespace soplex
Rational spxAbs(const Rational &r)
Absolute.
Definition: rational.cpp:3909
virtual void setup(SPxLPBase< Real > &lp)
clear and setup scaling arrays in the LP
Definition: spxscaler.cpp:124
bool isSetup() const
Returns setup status.
Definition: ssvectorbase.h:120
virtual R minAbsNzo(bool unscaled=true) const
Absolute smallest non-zero element in (possibly scaled) LP.
static const char * makename(bool doBoth)
Definition: spxequilisc.cpp:29
int size() const
Number of used indices.
Definition: svectorbase.h:152
static void updateRes(const SVSet facset, const SSVector resvecprev, SSVector &resvec, SSVector &tmpvec, Real eprev, Real qcurr)
Set of sparse vectors.
Least squares scaling.This SPxScaler implementation performs least squares scaling as suggested by Cu...
Definition: spxleastsqsc.h:36
virtual Real maxRowRatio(const SPxLPBase< Real > &lp) const
maximum ratio between absolute biggest and smallest element in any row.
Definition: spxscaler.cpp:855
SVectorBase< R > * create(int idxmax=0)
Creates new SVectorBase in set.
Definition: svsetbase.h:428
R & value(int n)
Reference to value of n &#39;th nonzero.
Definition: svectorbase.h:252
SPxLeastSqSC()
default constructor (this scaler makes no use of inherited member m_colFirst)
static void updateScale(const SSVector vecnnzeroes, const SSVector resnvec, SSVector &tmpvec, SSVector *&psccurr, SSVector *&pscprev, Real qcurr, Real qprev, Real eprev1, Real eprev2)
Wrapper for different output streams and verbosity levels.
int nRows() const
Returns number of rows in LP.
Definition: spxlpbase.h:151
DataArray< int > * m_activeColscaleExp
pointer to currently active column scaling factors
Definition: spxscaler.h:84
int nNzos() const
Returns number of nonzeros in LP.
Definition: spxlpbase.h:163
virtual void applyScaling(SPxLPBase< Real > &lp)
applies m_colscale and m_rowscale to the lp.
Definition: spxscaler.cpp:171
virtual void scale(SPxLP &lp, bool persistent=true)
Scale the loaded SPxLP.
static void initConstVecs(const SVSet *vecset, SVSet &facset, SSVector &veclogs, SSVector &vecnnzinv)
double Real
Definition: spxdefines.h:218
SPxOut * spxout
message handler
Definition: spxscaler.h:88
int & index(int n)
Reference to index of n &#39;th nonzero.
Definition: svectorbase.h:234
#define MSG_INFO2(spxout, x)
Prints out message x if the verbosity level is at least SPxOut::INFO2.
Definition: spxdefines.h:120
void clear()
Clears vector.
Definition: ssvectorbase.h:600
virtual Real maxAbsColscale() const
absolute biggest column scaling factor
Definition: spxscaler.cpp:773
virtual Real minAbsColscale() const
absolute smallest column scaling factor
Definition: spxscaler.cpp:760
virtual void setIntParam(int param, const char *name)
set int param (maximal conjugate gradient rounds)
#define MSG_INFO3(spxout, x)
Prints out message x if the verbosity level is at least SPxOut::INFO3.
Definition: spxdefines.h:122
DataArray< int > * m_activeRowscaleExp
pointer to currently active row scaling factors
Definition: spxscaler.h:85
void add(int i, const R &v)
Append one nonzero (i,v).
Definition: svectorbase.h:270
virtual Real minAbsRowscale() const
absolute smallest row scaling factor
Definition: spxscaler.cpp:787
LP least squares scaling.
Collection of dense, sparse, and semi-sparse vectors.
static void updateScaleFinal(const SSVector vecnnzeroes, const SSVector resnvec, SSVector &tmpvec, SSVector *&psccurr, SSVector *&pscprev, Real q, Real eprev1, Real eprev2)
Everything should be within this namespace.
SPxLeastSqSC & operator=(const SPxLeastSqSC &)
assignment operator
void setup()
Initializes nonzero indices for elements with absolute values above epsilon and sets all other elemen...
Definition: ssvectorbase.h:132
Sparse vectors.
virtual void setRealParam(Real param, const char *name)
set real param (conjugate gradient accuracy)
Semi sparse vector.
LP scaler abstract base class.Instances of classes derived from SPxScaler may be loaded to SoPlex in ...
Definition: spxscaler.h:76
void sort()
Sort nonzeros to increasing indices.
Definition: svectorbase.h:428
virtual Real maxAbsRowscale() const
absolute biggest row scaling factor
Definition: spxscaler.cpp:800
#define MSG_INFO1(spxout, x)
Prints out message x if the verbosity level is at least SPxOut::INFO1.
Definition: spxdefines.h:118
int nCols() const
Returns number of columns in LP.
Definition: spxlpbase.h:157
SSVectorBase< R > & assignPWproduct4setup(const SSVectorBase< S > &x, const SSVectorBase< T > &y)
Assigns pair wise vector product to SSVectorBase.
Definition: basevectors.h:511
const SVSetBase< R > * colSet() const
Returns the complete SVSetBase.
Definition: lpcolsetbase.h:68
bool isZero(Real a, Real eps=Param::epsilon())
returns true iff |a| <= eps
Definition: spxdefines.h:416
virtual Real maxColRatio(const SPxLPBase< Real > &lp) const
maximum ratio between absolute biggest and smallest element in any column.
Definition: spxscaler.cpp:817
const SVSetBase< R > * rowSet() const
Returns the complete SVSet.
Definition: lprowsetbase.h:69
void add(int i, R x)
Adds nonzero (i, x) to SSVectorBase.
Definition: ssvectorbase.h:208
int num() const
Current number of SVectorBases.
Definition: svsetbase.h:746
virtual R maxAbsNzo(bool unscaled=true) const
Absolute biggest non-zero element in (in rational case possibly scaled) LP.
SPxScaler & operator=(const SPxScaler &)
assignment operator
Definition: spxscaler.cpp:84
SSVectorBase< R > & assign2product4setup(const SVSetBase< S > &A, const SSVectorBase< T > &x)
Assigns SSVectorBase to for a setup x.
Definition: basevectors.h:602