1 // $Id: HOPSPACK_LinConstr.cpp 166 2010-03-22 19:58:07Z tplante $
2 // $URL: https://software.sandia.gov/svn/hopspack/trunk/src/src-shared/HOPSPACK_LinConstr.cpp $
3
4 //@HEADER
5 // ************************************************************************
6 //
7 // HOPSPACK: Hybrid Optimization Parallel Search Package
8 // Copyright 2009-2010 Sandia Corporation
9 //
10 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
11 // the U.S. Government retains certain rights in this software.
12 //
13 // This file is part of HOPSPACK.
14 //
15 // HOPSPACK is free software; you can redistribute it and/or modify
16 // it under the terms of the GNU Lesser General Public License as
17 // published by the Free Software Foundation; either version 2.1 of the
18 // License, or (at your option) any later version.
19 //
20 // This library is distributed in the hope that it will be useful, but
21 // WITHOUT ANY WARRANTY; without even the implied warranty of
22 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
23 // Lesser General Public License for more details.
24 //
25 // You should have received a copy of the GNU Lesser General Public
26 // License along with this library. If not, see http://www.gnu.org/licenses/.
27 //
28 // Questions? Contact Tammy Kolda (tgkolda@sandia.gov)
29 // or Todd Plantenga (tplante@sandia.gov)
30 //
31 // ************************************************************************
32 //@HEADER
33
34 /*!
35 @file HOPSPACK_LinConstr.cpp
36 @brief Implementation of HOPSPACK::LinConstr.
37 */
38
39 #include <iomanip>
40 #include <algorithm>
41 #include <math.h> //-- FOR fabs
42
43 #include "HOPSPACK_common.hpp"
44 #include "HOPSPACK_LinConstr.hpp"
45 #include "HOPSPACK_float.hpp"
46 #include "HOPSPACK_Matrix.hpp"
47 #include "HOPSPACK_ParameterList.hpp"
48 #include "HOPSPACK_ProblemDef.hpp"
49 #include "HOPSPACK_SolveLinConstrProj.hpp"
50 #include "HOPSPACK_Vector.hpp"
51
52
53 namespace HOPSPACK
54 {
55
LinConstr(const ProblemDef & probDef)56 LinConstr::LinConstr (const ProblemDef & probDef)
57 :
58 probDef (probDef),
59 scaling (probDef.getVarScaling())
60 {
61 //---- SET DEFAULTS IN CASE initialize IS NOT CALLED.
62
63 //---- THIS VALUE IS JUST SUFFICIENT FOR example2 ON ALL TEST MACHINES:
64 //---- _dActiveTol = 4.0 * HOPSPACK::getMachineEpsilon();
65 //---- HOWEVER, IT IS TOO TIGHT IF VARIABLES ARE SCALED TO THEIR LIMITS.
66 _dActiveTol = 1.0e-12;
67
68 _nDisplayFlag = 0;
69 return;
70 }
71
72
LinConstr(const LinConstr & cOther,const bool bDropEqs)73 LinConstr::LinConstr (const LinConstr & cOther,
74 const bool bDropEqs)
75 :
76 probDef (cOther.probDef),
77 _dActiveTol (cOther._dActiveTol),
78 _nDisplayFlag (cOther._nDisplayFlag),
79 scaling (cOther.probDef.getVarScaling()),
80 aIneq (cOther.aIneq),
81 bIneqLower (cOther.bIneqLower),
82 bIneqUpper (cOther.bIneqUpper)
83 {
84 if (bDropEqs == false)
85 {
86 aEq = cOther.aEq;
87 bEq = cOther.bEq;
88 }
89
90 if (setupScaledSystem() == false)
91 {
92 throwError ("constructor", "cannot set up scaled system");
93 }
94
95 return;
96 }
97
98
~LinConstr()99 LinConstr::~LinConstr()
100 {
101 }
102
initialize(const ParameterList & cLinConstrParams)103 bool LinConstr::initialize (const ParameterList & cLinConstrParams)
104 {
105 _dActiveTol = cLinConstrParams.getParameter ("Active Tolerance", _dActiveTol);
106
107 _nDisplayFlag = cLinConstrParams.getParameter ("Display", _nDisplayFlag);
108 if (_nDisplayFlag < 0)
109 _nDisplayFlag = 0;
110 if (_nDisplayFlag > 2)
111 _nDisplayFlag = 2;
112
113 if (setupMatrix (cLinConstrParams) == false)
114 return( false );
115 if (setupRhs (cLinConstrParams) == false)
116 return( false );
117 if (setupScaledSystem() == false)
118 return( false );
119
120 return( true );
121 }
122
123
124 // PRIVATE
setupMatrix(const ParameterList & params)125 bool LinConstr::setupMatrix(const ParameterList& params)
126 {
127 if (params.isParameterMatrix("Inequality Matrix"))
128 {
129 #if !defined(HAVE_LAPACK)
130 cerr << "ERROR: Cannot use linear constraints, no LAPACK in build" << endl;
131 throw INTERNAL_ERROR;
132 #endif
133
134 aIneq = params.getMatrixParameter("Inequality Matrix");
135 if ((aIneq.empty() == false) && (aIneq.getNcols() != scaling.size()))
136 {
137 cerr << "ERROR: Number of columns in 'Inequality Matrix' = "
138 << aIneq.getNcols() << " does not match number variables = "
139 << scaling.size() << endl;
140 return( false );
141 }
142
143 for (int i = 0; i < aIneq.getNrows(); i++)
144 {
145 const Vector v = aIneq.getRow (i);
146 for (int j = 0; j < v.size(); j++)
147 {
148 if (exists (v[j]) == false)
149 {
150 cerr << "ERROR: DNE value is not allowed in 'Inequality Matrix'"
151 << endl;
152 return( false );
153 }
154 }
155 }
156 }
157
158 if (params.isParameterMatrix("Equality Matrix"))
159 {
160 #if !defined(HAVE_LAPACK)
161 cerr << "ERROR: Cannot use linear constraints, no LAPACK in build" << endl;
162 throw INTERNAL_ERROR;
163 #endif
164
165 aEq = params.getMatrixParameter("Equality Matrix");
166 if ((aEq.empty() == false) && (aEq.getNcols() != scaling.size()))
167 {
168 cerr << "ERROR: Number of columns in 'Equality Matrix' = "
169 << aEq.getNcols() << " does not match number variables = "
170 << scaling.size() << endl;
171 return( false );
172 }
173
174 for (int i = 0; i < aEq.getNrows(); i++)
175 {
176 const Vector v = aEq.getRow (i);
177 for (int j = 0; j < v.size(); j++)
178 {
179 if (exists (v[j]) == false)
180 {
181 cerr << "ERROR: DNE value is not allowed in 'Equality Matrix'"
182 << endl;
183 return( false );
184 }
185 }
186 }
187 }
188
189 return( true );
190 }
191
192 // PRIVATE
setupRhs(const ParameterList & params)193 bool LinConstr::setupRhs(const ParameterList& params)
194 {
195 if (params.isParameterVector("Inequality Lower"))
196 bIneqLower = params.getVectorParameter("Inequality Lower");
197 else
198 bIneqLower.assign( aIneq.getNrows(), dne());
199 if (bIneqLower.size() != aIneq.getNrows())
200 {
201 cerr << "ERROR: Length of 'Inequality Lower' = " << bIneqLower.size()
202 << " does not match 'Inequality Matrix' = " << aIneq.getNrows()
203 << endl;
204 return( false );
205 }
206
207 if (params.isParameterVector("Inequality Upper"))
208 bIneqUpper = params.getVectorParameter("Inequality Upper");
209 else
210 bIneqUpper.assign(aIneq.getNrows(), dne());
211 if (bIneqUpper.size() != aIneq.getNrows())
212 {
213 cerr << "ERROR: Length of 'Inequality Upper' = " << bIneqUpper.size()
214 << " does not match 'Inequality Matrix' = " << aIneq.getNrows()
215 << endl;
216 return( false );
217 }
218
219 for (int i = 0; i < aIneq.getNrows(); i++)
220 {
221 if ((exists (bIneqLower[i]) == false) && (exists (bIneqUpper[i]) == false))
222 {
223 cerr << "ERROR: No bounds defined for inequality [" << (i + 1)
224 << "] in sublist 'Linear Constraints'" << endl;
225 return( false );
226 }
227 if ( (exists (bIneqLower[i])) && (exists (bIneqUpper[i]))
228 && (bIneqLower[i] > bIneqUpper[i]) )
229 {
230 cerr << "ERROR: Bounds are inconsistent for inequality [" << (i + 1)
231 << "] in sublist 'Linear Constraints'" << endl;
232 return( false );
233 }
234 }
235
236 if (params.isParameterVector("Equality Bounds"))
237 {
238 bEq = params.getVectorParameter("Equality Bounds");
239 if (bEq.size() != aEq.getNrows())
240 {
241 cerr << "ERROR: Length of 'Equality Bounds' = " << bEq.size()
242 << " does not match 'Equality Matrix' = " << aEq.getNrows()
243 << endl;
244 return( false );
245 }
246 for (int i = 0; i < bEq.size(); i++)
247 {
248 if (exists (bEq[i]) == false)
249 {
250 cerr << "ERROR: No bound defined for equality [" << (i + 1)
251 << "] in sublist 'Linear Constraints'" << endl;
252 return( false );
253 }
254 }
255 }
256 else if (aEq.empty() == false)
257 {
258 cerr << "ERROR: Need 'Equality Bounds' to go with 'Equality Matrix'" << endl;
259 return( false );
260 }
261
262 return( true );
263 }
264
265 // PRIVATE
setupScaledSystem()266 bool LinConstr::setupScaledSystem()
267 {
268 // *** Form lHat ***
269 const Vector & bLower = probDef.getLowerBnds();
270 for (int j = 0; j < scaling.size(); j++)
271 {
272 if (exists(bLower[j]))
273 lHat.push_back(bLower[j]);
274 else
275 lHat.push_back(0);
276 }
277
278 // *** Form bHatLower and bHatUpper ***
279
280 // Insert components corresponding to lower bounds
281 for (int j = 0; j < scaling.size() ; j++)
282 {
283 if (exists(bLower[j]))
284 bHatLower.push_back((bLower[j] - lHat[j]) / scaling[j]);
285 else
286 bHatLower.push_back(dne());
287 }
288
289 // Insert components corresponding to upper bounds
290 const Vector & bUpper = probDef.getUpperBnds();
291 for (int j = 0; j < scaling.size(); j++)
292 {
293 if (exists(bUpper[j]))
294 bHatUpper.push_back((bUpper[j] - lHat[j]) / scaling[j]);
295 else
296 bHatUpper.push_back(dne());
297 }
298
299 // Insert component corresponding to general linear constraints
300 if (!aIneq.empty())
301 {
302 // Compute ail = aIneq * lhat
303 Vector ail(aIneq.getNrows());
304 aIneq.multVec(lHat, ail);
305
306
307 for (int j = 0; j < aIneq.getNrows(); j++)
308 {
309 // Form scaled left-hand side.
310 if (exists(bIneqLower[j]))
311 bHatLower.push_back(bIneqLower[j] - ail[j]);
312 else
313 bHatLower.push_back(dne());
314
315 // Form scaled right-hand side.
316 if (exists(bIneqUpper[j]))
317 bHatUpper.push_back(bIneqUpper[j] - ail[j]);
318 else
319 bHatUpper.push_back(dne());
320 }
321 }
322
323
324 // *** Form aHat ***
325
326 // First add identity
327 aHat.setToIdentity(scaling.size());
328 // Now add in aTildeIneq
329 aHat.addMatrix(aIneq, scaling);
330
331 // *** Form aTildeEq and bTildeEq ***
332 // Nullspace matrix, Z, of aTildeEq is needed to compute distance to
333 // inequality constraints within nullspace.
334 Matrix ZT;
335 if (!aEq.empty())
336 {
337 // ael = aEq * lhat
338 Vector ael(aEq.getNrows());
339 aEq.multVec(lHat,ael);
340
341 for (int i = 0; i < aEq.getNrows(); i++)
342 bTildeEq.push_back(bEq[i] - ael[i]);
343
344 aTildeEq.scale(aEq, scaling);
345 aTildeEq.nullSpace(ZT, _dActiveTol);
346 }
347
348 // Error checks.
349 if (bHatLower.size() != aIneq.getNrows() + scaling.size())
350 {
351 cerr << "ERROR: Incorrect length for bHatLower <LinConstr.initialize()>"
352 << endl;
353 return( false );
354 }
355 if (bHatUpper.size() != aIneq.getNrows() + scaling.size())
356 {
357 cerr << "ERROR: Incorrect length for bHatUpper <LinConstr.initialize()>"
358 << endl;
359 return( false );
360 }
361 if ( (aHat.getNrows() != aIneq.getNrows() + scaling.size())
362 || (aHat.getNcols() != scaling.size()) )
363 {
364 cerr << "ERROR: Incorrect length for aHat <LinConstr.initialize()>"
365 << endl;
366 return( false );
367 }
368
369 // Store norms of rows for determining distances.
370 aHatZNorm.resize(aHat.getNrows());
371 // Project matrix into nullspace of scaled equality constraints.
372 Matrix AZ = aHat;
373 if (!ZT.empty())
374 AZ.multMat(ZT, Matrix::TRANSPOSE);
375 for (int i = 0; i < AZ.getNrows(); i++)
376 aHatZNorm[i] = AZ.getRow(i).norm();
377
378 return( true );
379 }
380
381
scale(Vector & x) const382 void LinConstr::scale(Vector& x) const
383 {
384 if (x.size() != scaling.size())
385 throwError("scale", "x vector has incorrect length");
386
387 for (int i = 0; i < scaling.size(); i++)
388 x[i] = (x[i] - lHat[i]) / scaling[i];
389 }
390
unscale(Vector & w) const391 void LinConstr::unscale(Vector& w) const
392 {
393 if (w.size() != scaling.size())
394 throwError("scale", "w vector has incorrect length");
395
396 for (int i = 0; i < scaling.size(); i++)
397 w[i] = (scaling[i] * w[i]) + lHat[i];
398 }
399
getActiveTol() const400 double LinConstr::getActiveTol() const
401 {
402 return _dActiveTol;
403 }
404
hasLinearConstraints() const405 bool LinConstr::hasLinearConstraints() const
406 {
407 if ((aIneq.getNrows() > 0) || (aEq.getNrows() > 0))
408 return( true );
409 return( false );
410 }
411
getAhat() const412 const Matrix & LinConstr::getAhat() const
413 {
414 return aHat;
415 }
416
getBhatLower() const417 const Vector & LinConstr::getBhatLower() const
418 {
419 return bHatLower;
420 }
421
getBhatUpper() const422 const Vector & LinConstr::getBhatUpper() const
423 {
424 return bHatUpper;
425 }
426
getAtildeEq() const427 const Matrix & LinConstr::getAtildeEq() const
428 {
429 return aTildeEq;
430 }
431
getBtildeEq() const432 const Vector & LinConstr::getBtildeEq() const
433 {
434 return bTildeEq;
435 }
436
437
isFeasible(const Vector & x,const bool bPrintViolationInfo) const438 bool LinConstr::isFeasible(const Vector& x,
439 const bool bPrintViolationInfo) const
440 {
441 int n = scaling.size();
442 // Error check
443 if (x.size() != n)
444 throwError("isFeasible", "x vector has incorrect length");
445
446 // Create scaled x
447 Vector xTilde(x);
448 scale(xTilde);
449
450 if (!isInequalityFeasible(xTilde, bPrintViolationInfo))
451 return false;
452
453 if (!isEqualityFeasible(xTilde, bPrintViolationInfo))
454 return false;
455
456 return true;
457 }
458
getL2Norm(const Vector & x) const459 double LinConstr::getL2Norm (const Vector & x) const
460 {
461 //---- COMPUTING WITH dnrm2 MAY BE MORE ACCURATE, BUT THEN
462 //---- THE METHOD REQUIRES LAPACK.
463
464 double dResult = 0.0;
465
466 //---- LOOP THRU VARIABLE BOUNDS.
467 const Vector & bLower = probDef.getLowerBnds();
468 const Vector & bUpper = probDef.getUpperBnds();
469 for (int i = 0; i < x.size(); i++)
470 {
471 if (exists (bLower[i]))
472 {
473 double dViolation = bLower[i] - x[i];
474 if (dViolation > dResult)
475 dResult += dViolation * dViolation;
476 }
477 if (exists (bUpper[i]))
478 {
479 double dViolation = x[i] - bUpper[i];
480 if (dViolation > dResult)
481 dResult += dViolation * dViolation;
482 }
483 }
484
485 //---- LOOP THRU EQUALITY CONSTRAINTS.
486 for (int i = 0; i < aEq.getNrows(); i++)
487 {
488 const Vector & nextCon = aEq.getRow (i);
489 double dConX = x.dot (nextCon);
490 double dViolation = fabs (dConX - bEq[i]);
491 dResult += dViolation * dViolation;
492 }
493
494 //---- LOOP THRU INEQUALITY CONSTRAINTS.
495 for (int i = 0; i < aIneq.getNrows(); i++)
496 {
497 const Vector & nextCon = aIneq.getRow(i);
498 double dConX = x.dot (nextCon);
499 double dViolation = 0.0;
500 if (exists (bIneqLower[i]))
501 {
502 if (dConX < bIneqLower[i])
503 dViolation = bIneqLower[i] - dConX;
504 }
505 if (exists (bIneqUpper[i]))
506 {
507 if (dConX > bIneqUpper[i])
508 dViolation = dConX - bIneqUpper[i];
509 }
510 dResult += dViolation * dViolation;
511 }
512
513 return( sqrt (dResult) );
514 }
515
getLInfNorm(const Vector & x) const516 double LinConstr::getLInfNorm (const Vector & x) const
517 {
518 double dResult = 0.0;
519
520 //---- LOOP THRU VARIABLE BOUNDS.
521 const Vector & bLower = probDef.getLowerBnds();
522 const Vector & bUpper = probDef.getUpperBnds();
523 for (int i = 0; i < x.size(); i++)
524 {
525 if (exists (bLower[i]))
526 {
527 double dViolation = bLower[i] - x[i];
528 if (dViolation > dResult)
529 dResult = dViolation;
530 }
531 if (exists (bUpper[i]))
532 {
533 double dViolation = x[i] - bUpper[i];
534 if (dViolation > dResult)
535 dResult = dViolation;
536 }
537 }
538
539 //---- LOOP THRU EQUALITY CONSTRAINTS.
540 for (int i = 0; i < aEq.getNrows(); i++)
541 {
542 const Vector & nextCon = aEq.getRow (i);
543 double dConX = x.dot (nextCon);
544 double dViolation = fabs (dConX - bEq[i]);
545 if (dViolation > dResult)
546 dResult = dViolation;
547 }
548
549 //---- LOOP THRU INEQUALITY CONSTRAINTS.
550 for (int i = 0; i < aIneq.getNrows(); i++)
551 {
552 const Vector & nextCon = aIneq.getRow(i);
553 double dConX = x.dot (nextCon);
554 double dViolation = 0.0;
555 if (exists (bIneqLower[i]))
556 {
557 dViolation = bIneqLower[i] - dConX;
558 if (dViolation > dResult)
559 dResult = dViolation;
560 }
561 if (exists (bIneqUpper[i]))
562 {
563 dViolation = dConX - bIneqUpper[i];
564 if (dViolation > dResult)
565 dResult = dViolation;
566 }
567 }
568
569 return( dResult );
570 }
571
572
maxStep(const Vector & x,const Vector & d,double maxLength) const573 double LinConstr::maxStep(const Vector& x,
574 const Vector& d,
575 double maxLength) const
576 {
577 double maxStep = maxLength;
578 int n = scaling.size();
579
580 Vector xTilde(x);
581 scale(xTilde);
582
583 const Vector & bLower = probDef.getLowerBnds();
584 const Vector & bUpper = probDef.getUpperBnds();
585
586 // Determine max step for bounds component
587
588 for (int i = 0; i < scaling.size(); i++)
589 {
590 if (d[i] < -scaling[i] * _dActiveTol) // Points to lower bound
591 switch (getIneqState(i, LOWER_BOUND, xTilde))
592 {
593 case ACTIVE: // Can't move if the constraint is active
594 return 0;
595 case VIOLATED: // Or violated
596 return 0;
597 case INACTIVE:
598 maxStep = min(maxStep, (bLower[i] - x[i]) / d[i]);
599 break;
600 case DNE: // This means there is no lower bound
601 break;
602 }
603 else if (d[i] > scaling[i] * _dActiveTol) // Points to upper bound
604 switch (getIneqState(i, UPPER_BOUND, xTilde))
605 {
606 case ACTIVE: // Can't move if the constraint is active
607 return 0;
608 case VIOLATED: // Or violated
609 return 0;
610 case INACTIVE:
611 maxStep = min(maxStep, (bUpper[i] - x[i]) / d[i]);
612 break;
613 case DNE: // This means there is no upper bound
614 break;
615 }
616 }
617
618 // Determine max step for inequality component.
619
620 if (!aIneq.empty())
621 {
622 int p = aIneq.getNrows();
623
624 // aix = aIneq * x
625 Vector aix(p);
626 aIneq.multVec(x, aix);
627
628 // aid = aIneq * d
629 Vector aid(p);
630 aIneq.multVec(d, aid);
631
632 for (int j = 0; j < p; j++)
633 {
634 if (aid[j] < -_dActiveTol) // Points to lower bound
635 switch (getIneqState(n + j, LOWER_BOUND, xTilde))
636 {
637 case ACTIVE: // Can't move if the constraint is active
638 return 0;
639 case VIOLATED: // Or violated
640 return 0;
641 case INACTIVE:
642 maxStep = min(maxStep, (bIneqLower[j] - aix[j]) / aid[j]);
643 break;
644 case DNE: // This means there is no lower bound
645 break;
646 }
647 else if (aid[j] > _dActiveTol) // Points to upper bound
648 switch (getIneqState(n + j, UPPER_BOUND, xTilde))
649 {
650 case ACTIVE: // Can't move if the constraint is active
651 return 0;
652 case VIOLATED: // Or violated
653 return 0;
654 case INACTIVE:
655 maxStep = min(maxStep, (bIneqUpper[j] - aix[j]) / aid[j]);
656 break;
657 case DNE: // This means there is no upper bound
658 break;
659 }
660 }
661 }
662
663 // Process equality constraints.
664
665 if (!aEq.empty())
666 {
667 int p = aEq.getNrows();
668
669 // z = aEq * d
670 Vector z(p);
671 aEq.multVec(d,z);
672
673 // Check that direction stays on hyperplane defined by the
674 // equality constraints
675 for (int i = 0; i < p; i++)
676 if (fabs(z[i]) > _dActiveTol)
677 return 0;
678 }
679
680 return maxStep;
681 }
682
getActiveIneqIndices(const Vector & xdist,double epsilon,vector<ActiveType> & index) const683 void LinConstr::getActiveIneqIndices(const Vector& xdist,
684 double epsilon,
685 vector<ActiveType>& index) const
686
687 {
688 int nrows = aHat.getNrows();
689 // Size the index
690 index.resize(nrows);
691
692 // Determine which bounds are near.
693 for (int i = 0; i < nrows; i++)
694 {
695 // Check if upper or lower bounds are near.
696 bool nearlow = (xdist[i] < epsilon);
697 bool nearupp = (xdist[i+nrows] < epsilon);
698
699 // Update index.
700 if (nearlow && nearupp)
701 index[i] = BOTH_ACTIVE;
702 else if (nearlow) // and not nearupp
703 index[i] = LOWER_ACTIVE;
704 else if (nearupp) // and not nearlow
705 index[i] = UPPER_ACTIVE;
706 else
707 index[i] = NEITHER_ACTIVE;
708 }
709 }
710
formDistanceVector(const Vector & x,Vector & xdist) const711 void LinConstr::formDistanceVector(const Vector& x, Vector& xdist) const
712 {
713 // Transform x to scaled space..
714 Vector xTilde = x;
715 scale(xTilde);
716
717 int nrows = aHat.getNrows();
718
719 // z = aHat*xTilde.
720 Vector z(nrows);
721 aHat.multVec(xTilde, z);
722
723 // xdist must be twice the size of z to hold upper and lower bound info.
724 xdist.resize(2*nrows);
725 for (int i = 0; i < nrows; i++)
726 {
727 // Process lower bound.
728 if (exists(bHatLower[i]))
729 {
730 if (aHatZNorm[i] > _dActiveTol)
731 xdist[i] = fabs( z[i] - bHatLower[i] ) / aHatZNorm[i];
732 else if (fabs( z[i] - bHatLower[i] ) < _dActiveTol)
733 xdist[i] = 0;
734 else
735 xdist[i] = dne();
736 }
737 else
738 xdist[i] = dne();
739
740 // Process upper bound.
741 if (exists(bHatUpper[i]))
742 {
743 if (aHatZNorm[i] > _dActiveTol)
744 xdist[i+nrows] = fabs( z[i] - bHatUpper[i] ) / aHatZNorm[i];
745 else if (fabs( z[i] - bHatUpper[i] ) < _dActiveTol)
746 xdist[i+nrows] = 0;
747 else
748 xdist[i+nrows] = dne();
749 }
750 else
751 xdist[i+nrows] = dne();
752 }
753 }
754
snapToBoundary(Vector & x,double esnap) const755 void LinConstr::snapToBoundary(Vector& x,
756 double esnap) const
757 {
758 // Form scaled x.
759 Vector xTilde = x;
760 scale(xTilde);
761
762 // Form snap-to system.
763 Matrix Asnap;
764 Vector bsnap;
765 formSnapSystem(xTilde, esnap, Asnap, bsnap);
766
767 // Find closest point satisfying snap constraints.
768 if (Asnap.specialConstrainedLSQR (xTilde, bsnap) == false)
769 return;
770
771 // Return the result.
772 unscale(xTilde);
773 x = xTilde;
774 return;
775 }
776
777 // PRIVATE
formSnapSystem(const Vector & xTilde,double esnap,Matrix & Asnap,Vector & bsnap) const778 void LinConstr::formSnapSystem(const Vector& xTilde, double esnap,
779 Matrix& Asnap,
780 Vector& bsnap) const
781 {
782 Asnap.clear();
783 bsnap.resize(0);
784
785 if (aHat.empty())
786 {
787 if (!isEqualityFeasible(xTilde))
788 {
789 Asnap = aTildeEq;
790 bsnap = bTildeEq;
791 }
792 return;
793 }
794
795 int nrow = aHat.getNrows();
796 int nvar = aHat.getNcols();
797 // Form z = Ahat*xTilde;
798 Vector z(nrow);
799 aHat.multVec(xTilde, z);
800
801 // Measure and sort distance in terms of scaled space.
802 multimap<double,int> dmap;
803 multimap<double,int>::iterator iter;
804 for (int i = 0; i < nrow; i++)
805 {
806 if (exists(bHatLower[i]))
807 {
808 double distlow = fabs( z[i] - bHatLower[i] ) / aHatZNorm[i];
809 dmap.insert(std::pair<const double, int>(distlow, i));
810 }
811 if (exists(bHatUpper[i]))
812 {
813 // Indices for upper bounds in dmap denoted by i+nrow.
814 // Symbolically stacking matrix.
815 double distupp = fabs( z[i] - bHatUpper[i] ) / aHatZNorm[i];
816 dmap.insert(std::pair<const double, int>(distupp, i+nrow));
817 }
818 }
819
820 // Erase constraints which are outside the esnap ball.
821 dmap.erase(dmap.upper_bound(esnap),dmap.end());
822
823 if ( (dmap.lower_bound(_dActiveTol) == dmap.end())
824 && isEqualityFeasible(xTilde) )
825 return; // There are no constraints within distance esnap which are inactive.
826
827 // Equality constraints are always added.
828 Asnap = aTildeEq;
829 bsnap = bTildeEq;
830
831 // Now add inequality constraints within distance esnap.
832 // Recall that we have symbolically stacked aHat, so we must mod out by nrow.
833 for (iter=dmap.begin(); iter != dmap.end(); iter++)
834 {
835 if (bsnap.size() >= nvar)
836 break;
837 int irow = iter->second;
838 if (irow < nrow)
839 { // Snap to lower bound.
840 Asnap.addRow(aHat.getRow(irow));
841 bsnap.push_back(bHatLower[irow]);
842 }
843 else
844 { // Snap to upper bound.
845 Asnap.addRow(aHat.getRow(irow % nrow));
846 bsnap.push_back(bHatUpper[irow % nrow]);
847 }
848 }
849
850 // Now ensure that rank(Asnap) = Asnap.getNrows().
851 Asnap.pruneDependentRows(bsnap,_dActiveTol);
852 }
853
854
projectToFeasibility(Vector & x) const855 bool LinConstr::projectToFeasibility (Vector & x) const
856 {
857 SolveLinConstrProj solver;
858 Vector xSolution;
859 if (solver.solve (probDef, *this, x, xSolution) == false)
860 return( false );
861
862 x = xSolution;
863 return( true );
864 }
865
866
printDefinition(const bool bDisplayFull) const867 void LinConstr::printDefinition (const bool bDisplayFull) const
868 {
869 if (_nDisplayFlag <= 0)
870 return;
871
872 if ((_nDisplayFlag < 2) || (bDisplayFull == false))
873 {
874 cout << "Linear Constraints" << endl;
875 }
876 else
877 {
878 cout << "Linear Constraints (full display)" << endl;
879 cout << " (Variable bounds are displayed in the Problem Definition)"
880 << endl;
881 }
882
883 printCounts_();
884 cout << " Tolerance for feasibility, active constraints = "
885 << setw(14) << setprecision(6)
886 << setiosflags (ios::scientific) << getActiveTol() << endl;
887
888 if ((_nDisplayFlag == 2) && (bDisplayFull == true))
889 {
890 //---- PRINT A LIST WITH EVERY INEQUALITY.
891 if (aIneq.empty() == false)
892 {
893 cout << " Inequality constraints:" << endl;
894 for (int i = 0; i < aIneq.getNrows(); i++)
895 {
896 cout << " ";
897
898 if (exists (bIneqLower[i]))
899 {
900 cout << setw(14) << setprecision(6)
901 << setiosflags (ios::scientific) << bIneqLower[i];
902 cout << " <= ";
903 }
904 else
905 cout << " " << " ";
906
907 printIneqName_ (i);
908
909 if (exists (bIneqUpper[i]))
910 {
911 cout << " <= ";
912 cout << setw(14) << setprecision(6)
913 << setiosflags (ios::scientific) << bIneqUpper[i];
914 }
915 else
916 cout << " " << " ";
917
918 cout << endl;
919 }
920 cout << " A_ineq = ";
921 aIneq.formattedPrint (" ", cout);
922 cout << endl;
923 }
924
925 //---- PRINT A LIST WITH EVERY EQUALITY.
926 if (bEq.empty() == false)
927 {
928 cout << " Equality constraints:" << endl;
929 for (int i = 0; i < bEq.size(); i++)
930 {
931 cout << " ";
932 printEqName_ (i);
933 cout << " = " << setw(14) << setprecision(6)
934 << setiosflags (ios::scientific) << bEq[i] << endl;
935 }
936 cout << " A_eq = ";
937 aEq.formattedPrint (" ", cout);
938 cout << endl;
939 }
940
941 cout << "End of Linear Constraints (full display)" << endl;
942 }
943
944 cout << endl;
945 return;
946 }
947
948 // PRIVATE
printCounts_(void) const949 void LinConstr::printCounts_ (void) const
950 {
951 int nNumIneqLower = 0;
952 int nNumIneqUpper = 0;
953 for (int i = 0; i < bIneqLower.size(); i++)
954 {
955 if (exists (bIneqLower[i]))
956 nNumIneqLower++;
957 if (exists(bIneqUpper[i]))
958 nNumIneqUpper++;
959 }
960
961 cout << " Constraint count summary:" << endl;
962 cout << " " << setw (5) << scaling.size() << " variables" << endl;
963 cout << " " << setw (5) << (nNumIneqLower + nNumIneqUpper)
964 << " inequality constraints" << endl;
965 cout << " " << setw (5) << bEq.size() << " equality constraints" << endl;
966 return;
967 }
968
969 // PRIVATE
printIneqName_(const int nIneqNum) const970 void LinConstr::printIneqName_ (const int nIneqNum) const
971 {
972 cout << "c_ineq[" << setw (3) << nIneqNum << "]";
973 return;
974 }
975
976 // PRIVATE
printEqName_(const int nEqNum) const977 void LinConstr::printEqName_ (const int nEqNum) const
978 {
979 cout << "c_eq[" << setw (3) << nEqNum << "]";
980 return;
981 }
982
983 // PRIVATE
throwError(const string & fname,const string & msg) const984 void LinConstr::throwError(const string& fname, const string& msg) const
985 {
986 cerr << "ERROR: " << msg << " <" << fname << ">" << endl;
987 throw INTERNAL_ERROR;
988 }
989
990 // PRIVATE
getIneqState(const int i,const BoundType bType,const Vector & xTilde,const bool bPrintViolationInfo) const991 LinConstr::StateType LinConstr::getIneqState
992 (const int i,
993 const BoundType bType,
994 const Vector& xTilde,
995 const bool bPrintViolationInfo) const
996 {
997 // a = constraint row
998 const Vector& a = aHat.getRow(i);
999 double anorm = aHatZNorm[i];
1000 // b = lhs
1001 double b = (bType == LOWER_BOUND) ? bHatLower[i] : bHatUpper[i];
1002
1003 // Check finiteness of the bound
1004 if (!exists(b))
1005 return DNE;
1006
1007 // z = a' * xTilde
1008 double z;
1009 z = xTilde.dot(a);
1010 double xnorm = xTilde.norm();
1011
1012 // Check if the constraint is epsilon-active
1013 if ( fabs(z - b) < (_dActiveTol * max (anorm, xnorm)) )
1014 return ACTIVE;
1015
1016 // Check if the constraint is otherwise satisfied
1017 if (((bType == LOWER_BOUND) && ( b <= z )) ||
1018 ((bType == UPPER_BOUND) && ( z <= b )))
1019 return INACTIVE;
1020
1021 // Otherwise, it must be violated.
1022 if (bPrintViolationInfo)
1023 {
1024 cout << " Inequality[" << i << "] violated by " << fabs(z - b)
1025 << " (tolerance = " << (_dActiveTol * max (anorm, xnorm)) << ")"
1026 << endl;
1027 }
1028 return VIOLATED;
1029 }
1030
1031 // PRIVATE
getEqState(const int i,const Vector & xTilde,const bool bPrintViolationInfo) const1032 LinConstr::StateType LinConstr::getEqState
1033 (const int i,
1034 const Vector& xTilde,
1035 const bool bPrintViolationInfo) const
1036 {
1037 // a = constraint row
1038 const Vector& a = aTildeEq.getRow(i);
1039 double anorm = a.norm();
1040
1041 // b = rhs
1042 double b = bTildeEq[i];
1043
1044 // z = a' * xTilde
1045 double z;
1046 z = xTilde.dot(a);
1047 double xnorm = xTilde.norm();
1048 // Check if the constraint is epsilon-active
1049 if ( fabs(z - b) < (_dActiveTol * max (anorm, xnorm)) )
1050 return ACTIVE;
1051
1052 // Otherwise, it must be violated.
1053 if (bPrintViolationInfo)
1054 {
1055 cout << " Equality[" << i << "] violated by " << fabs(z - b)
1056 << " (tolerance = " << (_dActiveTol * max (anorm, xnorm)) << ")"
1057 << endl;
1058 }
1059 return VIOLATED;
1060 }
1061
1062 // PRIVATE
isEqualityFeasible(const Vector & xTilde,const bool bPrintViolationInfo) const1063 bool LinConstr::isEqualityFeasible(const Vector& xTilde,
1064 const bool bPrintViolationInfo) const
1065 {
1066 // Check equality constraints
1067 for (int i = 0; i < aTildeEq.getNrows(); i ++)
1068 {
1069 if (VIOLATED == getEqState(i, xTilde, bPrintViolationInfo))
1070 return false;
1071 }
1072
1073 return true;
1074 }
1075
1076 // PRIVATE
isInequalityFeasible(const Vector & xTilde,const bool bPrintViolationInfo) const1077 bool LinConstr::isInequalityFeasible(const Vector& xTilde,
1078 const bool bPrintViolationInfo) const
1079 {
1080 // Check inequality constraints
1081 for (int i = 0; i < aHat.getNrows(); i ++)
1082 {
1083 if (VIOLATED == getIneqState(i, UPPER_BOUND, xTilde, bPrintViolationInfo))
1084 return false;
1085 if (VIOLATED == getIneqState(i, LOWER_BOUND, xTilde, bPrintViolationInfo))
1086 return false;
1087 }
1088
1089 return true;
1090 }
1091
1092 } //-- namespace HOPSPACK
1093