Sandia Home Sandia Home
Main Page | Publications | Downloads | Configuration | Running the Code | Solver Parameters | FAQ | Namespace List | Class Hierarchy | Class List | File List | Namespace Members | Class Members | File Members | Related Pages

APPSPACK_Constraints_Linear.cpp

Go to the documentation of this file.
00001 // $Id: APPSPACK_Constraints_Linear.cpp,v 1.61.2.1 2007/01/14 21:51:27 jgriffi Exp $ 
00002 // $Source: /space/CVS-Acro/acro/packages/appspack/appspack/src/APPSPACK_Constraints_Linear.cpp,v $ 
00003 
00004 //@HEADER
00005 // ************************************************************************
00006 // 
00007 //          APPSPACK: Asynchronous Parallel Pattern Search
00008 //                 Copyright (2003) Sandia Corporation
00009 // 
00010 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00011 // license for use of this work by or on behalf of the U.S. Government.
00012 // 
00013 // This library is free software; you can redistribute it and/or modify
00014 // it under the terms of the GNU Lesser General Public License as
00015 // published by the Free Software Foundation; either version 2.1 of the
00016 // License, or (at your option) any later version.
00017 //  
00018 // This library is distributed in the hope that it will be useful, but
00019 // WITHOUT ANY WARRANTY; without even the implied warranty of
00020 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00021 // Lesser General Public License for more details.
00022 //                                                                                 
00023 // You should have received a copy of the GNU Lesser General Public
00024 // License along with this library; if not, write to the Free Software
00025 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00026 // USA.                                                                           .
00027 // 
00028 // Questions? Contact Tammy Kolda (tgkolda@sandia.gov) 
00029 // 
00030 // ************************************************************************
00031 //@HEADER
00032 
00038 #include "APPSPACK_Vector.hpp"
00039 #include "APPSPACK_Constraints_Linear.hpp"
00040 #include "APPSPACK_Print.hpp"
00041 #include "APPSPACK_Utils.hpp"
00042 #include "APPSPACK_Common.hpp"
00043 #include "APPSPACK_Float.hpp"
00044 #include "APPSPACK_LAPACK_Wrappers.hpp"
00045 #include "APPSPACK_CDDLIB_Interface.h"
00046 
00047 APPSPACK::Constraints::Linear::Linear(Parameter::List& params) :
00048   epsMach(params.getParameter("Machine Epsilon", 1.0e-12)),
00049   displayMatrix(params.getParameter("Display Linear Systems", true))
00050 {
00051   setup(params);
00052   errorCheck();
00053 }
00054 
00055 APPSPACK::Constraints::Linear::~Linear()
00056 {
00057 }
00058 
00059 void APPSPACK::Constraints::Linear::setup(Parameter::List& params)
00060 {  
00061   setupBounds(params);
00062   setupScaling(params);
00063   setupMatrix(params); 
00064   setupRhs(params);
00065   setupScaledSystem();
00066 }
00067 
00068 void APPSPACK::Constraints::Linear::setupBounds(APPSPACK::Parameter::List& params)
00069 { 
00070 
00071   vector<bool> isLower;
00072   vector<bool> isUpper;
00073 
00074   if (params.isParameterVector("Scaling"))
00075     scaling = params.getVectorParameter("Scaling");
00076   
00077   if (params.isParameterVector("Upper"))
00078     bUpper = params.getVectorParameter("Upper");
00079   
00080   if (params.isParameterVector("Lower"))
00081     bLower = params.getVectorParameter("Lower");
00082   
00083   if (params.isParameterVector("Is Lower"))
00084     params.getVectorParameter("Is Lower").copyToBool(isLower);
00085   
00086   if (params.isParameterVector("Is Upper"))
00087     params.getVectorParameter("Is Upper").copyToBool(isUpper);
00088   
00089   // Determine number of variables
00090   int nvar = max(scaling.size(), bLower.size());
00091   if (nvar == 0)
00092     error("setup", "Neither \"Scaling\" nor \"Lower\" is defined in the parameter list");
00093 
00094   // Fill lower (if necessary)
00095   if (bLower.empty())
00096   {
00097     bLower.assign(nvar, dne());
00098     params.getParameter("Lower", bLower); 
00099   }
00100 
00101   // Use isLower if it's defined (for backward compatibility)
00102   if (!isLower.empty())
00103   {
00104     if (isLower.size() != nvar)
00105       error("setupBounds", "\"Is Lower\" is the wrong size");
00106 
00107     for (int i = 0; i < nvar; i ++)
00108       if (isLower[i])
00109       {
00110         if (!exists(bLower[i]))
00111           error("setupBounds", "\"Is Lower\" has a true (nonzero) entry but \"Lower\" is not defined");    
00112       }
00113       else
00114       {
00115         bLower[i] = dne();
00116       }
00117   }
00118 
00119   // Fill upper 
00120   if (bUpper.empty())
00121   {
00122     bUpper.assign(nvar, dne());
00123     params.getParameter("Upper", bUpper);
00124   }
00125 
00126   // Use isUpper if it's defined (for backward compatibility)
00127   if (!isUpper.empty())
00128   {
00129     if (isUpper.size() != nvar)
00130       error("setupBounds", "\"Is Upper\" is the wrong size");
00131 
00132     for (int i = 0; i < nvar; i ++)
00133       if (isUpper[i])
00134       {
00135         if (!exists(bUpper[i]))
00136           error("setupBounds", "\"Is Upper\" has a true (nonzero) entry but \"Upper\" is not defined");    
00137       }
00138       else
00139       {
00140         bUpper[i] = dne();
00141       }
00142   }
00143 
00144   // Error check
00145   if (bLower.size() != nvar)
00146     error("setupBounds", "\"Lower\" vector has incorrect length");
00147   if (bUpper.size() != nvar)
00148     error("setupBounds", "\"Upper\" vector has incorrect length");
00149 }
00150 
00151 void APPSPACK::Constraints::Linear::setupScaling(APPSPACK::Parameter::List& params)
00152 {
00153   // Fix the scaling (if necessary)
00154   if (scaling.empty())
00155   {  
00156     // Check that all the bounds are specified
00157     for (int i = 0; i < bLower.size(); i++)
00158       if ((!exists(bLower[i])) || (!exists(bUpper[i])))
00159         error("setup", "\"Scaling\" must be specified if upper and lower bounds are not fully specified");
00160     
00161     for (int i = 0; i < bLower.size(); i++)
00162       scaling.push_back(bUpper[i] - bLower[i]);
00163     
00164     params.getParameter("Scaling", scaling);
00165   }
00166   
00167   // Error check
00168   if (scaling.size() != bLower.size())
00169    error("setupScaling", "\"Scaling\" vector has incorrect length");
00170 }
00171 
00172 void APPSPACK::Constraints::Linear::setupMatrix(APPSPACK::Parameter::List& params)
00173 {
00174   // Get inequality matrix
00175   if (params.isParameterMatrix("Inequality Matrix"))
00176   {
00177     aIneq = params.getMatrixParameter("Inequality Matrix");
00178 
00179     // Error check
00180     if ((!aIneq.empty()) && (aIneq.getNcols() != scaling.size()))
00181       error("setupMatrix", "\"Inequality Matrix\" matrix has incorrect size.");
00182   }
00183 
00184   // Get equality matrix
00185   if (params.isParameterMatrix("Equality Matrix"))
00186   {
00187     aEq = params.getMatrixParameter("Equality Matrix");
00188 
00189     // Error check
00190     if ((!aEq.empty()) && (aEq.getNcols() != scaling.size()))
00191         error("setupMatrix", "\"Equality Matrix\" matrix has incorrect size.");
00192   }
00193 }
00194 
00195 void APPSPACK::Constraints::Linear::setupRhs(APPSPACK::Parameter::List& params)
00196 {  
00197   // Get lower bounds on linear inequality constraints
00198   if (params.isParameterVector("Inequality Lower"))
00199     bIneqLower = params.getVectorParameter("Inequality Lower");
00200   else
00201     bIneqLower.assign( aIneq.getNrows(), dne());
00202 
00203   // Get upper bounds on linear inequality constraints
00204   if (params.isParameterVector("Inequality Upper"))
00205     bIneqUpper = params.getVectorParameter("Inequality Upper");
00206   else
00207     bIneqUpper.assign(aIneq.getNrows(), dne());
00208 
00209   // Get upper bounds on linear equality constraints
00210   if (params.isParameterVector("Equality Bound"))
00211     bEq = params.getVectorParameter("Equality Bound");
00212 
00213   // Error check
00214   if (bIneqLower.size() != aIneq.getNrows())
00215     error("setupRhs", "\"Inequality Lower\" vector has incorrect length");
00216 
00217   // Error check
00218   if (bIneqUpper.size() != aIneq.getNrows())
00219     error("setupRhs", "\"Inequality Upper\" vector has incorrect length");
00220 
00221   // Error check
00222   if (bEq.size() != aEq.getNrows())
00223     error("setupRhs", "\"Equality Bound\" vector has incorrect length");
00224 }
00225 
00226 void APPSPACK::Constraints::Linear::setupScaledSystem()
00227 {  
00228   // *** Form lHat ***
00229   for (int j = 0; j < scaling.size(); j++)
00230   {
00231     if (exists(bLower[j]))
00232       lHat.push_back(bLower[j]);
00233     else
00234       lHat.push_back(0);
00235   }
00236 
00237   // *** Form bHatLower and bHatUpper ***
00238 
00239   // Insert components corresponding to lower bounds
00240   for (int j = 0; j < scaling.size() ; j++)
00241   {
00242     if (exists(bLower[j]))
00243       bHatLower.push_back((bLower[j] - lHat[j]) / scaling[j]);
00244     else
00245       bHatLower.push_back(dne());
00246   }
00247 
00248   // Insert components corresponding to upper bounds
00249   for (int j = 0; j < scaling.size(); j++)
00250   {
00251     if (exists(bUpper[j]))
00252       bHatUpper.push_back((bUpper[j] - lHat[j]) / scaling[j]);
00253     else
00254       bHatUpper.push_back(dne());
00255   }
00256 
00257   // Insert component corresponding to general linear constraints
00258   if (!aIneq.empty())
00259   {
00260     // Compute ail = aIneq * lhat
00261     Vector ail(aIneq.getNrows());
00262     aIneq.multVec(lHat, ail);
00263 
00264     
00265     for (int j = 0; j < aIneq.getNrows(); j++)
00266     {
00267       // Form scaled left-hand side.
00268       if (exists(bIneqLower[j]))
00269         bHatLower.push_back(bIneqLower[j] - ail[j]);
00270       else
00271         bHatLower.push_back(dne());
00272 
00273       // Form scaled right-hand side.
00274       if (exists(bIneqUpper[j]))
00275         bHatUpper.push_back(bIneqUpper[j] - ail[j]);
00276       else
00277         bHatUpper.push_back(dne());
00278     }
00279   }
00280   
00281   // *** Form aHat ***
00282 
00283   // First add identity
00284   aHat.setToIdentity(scaling.size());
00285   // Now add in aTildeIneq
00286   aHat.addMatrix(aIneq, scaling);
00287   
00288   // *** Form aTildeEq and bTildeEq ***
00289   // Nullspace matrix, Z, of aTildeEq is needed to compute distance to
00290   // inequality constraints within nullspace.
00291   Matrix ZT;
00292   if (!aEq.empty())
00293   {
00294     // ael = aEq * lhat
00295     Vector ael(aEq.getNrows());
00296     aEq.multVec(lHat,ael);
00297 
00298     for (int i = 0; i < aEq.getNrows(); i++)
00299       bTildeEq.push_back(bEq[i] - ael[i]);
00300 
00301     aTildeEq.scale(aEq, scaling);
00302     aTildeEq.nullSpace(ZT, epsMach);
00303   }
00304 
00305   // Error check.
00306   if (bHatLower.size() != aIneq.getNrows() + scaling.size())
00307     error("setupScaledSystem", "\"bHatUpper\" vector has incorrect length");
00308   // Error check.
00309   if (bHatUpper.size() != aIneq.getNrows() + scaling.size())
00310     error("setupScaledSystem", "\"bHatUpper\" vector has incorrect length");
00311   // Error check.
00312   if ((aHat.getNrows() != aIneq.getNrows() + scaling.size()) || (aHat.getNcols() != scaling.size()))
00313     error("setupScaledSystem", "\"aHat\" matrix has incorrect size");
00314 
00315   // Store norms of rows for determining distances.
00316   aHatZNorm.resize(aHat.getNrows());
00317   aHatNorm.resize(aHat.getNrows());
00318 
00319   // Project matrix into nullspace of scaled equality constraints.
00320   Matrix AZ = aHat;
00321   if (!ZT.empty())
00322     AZ.multMat(ZT, Matrix::Transpose);
00323   for (int i = 0; i < aHat.getNrows(); i++)
00324   {
00325     aHatZNorm[i] = AZ.getRow(i).norm();
00326     aHatNorm[i] = aHat.getRow(i).norm();
00327   }
00328 }
00329 
00330 void APPSPACK::Constraints::Linear::errorCheck() const
00331 {
00332   int nvar = scaling.size();
00333 
00334   // Positive scaling values
00335   for (int i = 0; i < nvar; i ++)
00336     if (scaling[i] <= 0)
00337       error("errorCheck", "Negative \"Scaling\" value.");
00338 
00339   // Consistency of bounds
00340   for (int i = 0; i < nvar; i ++)
00341     if ((exists(bLower[i])) && (exists(bUpper[i])) && (bLower[i] > bUpper[i]))
00342       error("errorCheck", "Upper bound is less than lower bound.");
00343 
00344   // Consistency of constraints
00345   for (int i = 0; i < aIneq.getNrows(); i ++)
00346     if ((exists(bIneqLower[i])) && (exists(bIneqUpper[i])) && (bIneqLower[i] > bIneqUpper[i]))
00347       error("errorCheck", "Upper bound is less than lower bound for linear constraint.");
00348 
00349 #if !defined(HAVE_LAPACK) && !defined(HAVE_CDDLIB)
00350     // Error if trying to use general linear constraints without 
00351     // at least one of LAPACK or CDDLIB
00352     if (!(aEq.empty() && aIneq.empty()))
00353     {
00354       string errorMsg = "General linear constraint capability currently disabled.\n";
00355       errorMsg += "\nPlease reconfigure with LAPACK and/or CDDLIB.\n";
00356       error("errorCheck", errorMsg);
00357     }
00358 #endif
00359 }
00360 
00361 double APPSPACK::Constraints::Linear::getEpsMach() const
00362 {
00363   return epsMach;
00364 }
00365 
00366 const APPSPACK::Vector& APPSPACK::Constraints::Linear::getScaling() const
00367 {
00368   return scaling;
00369 }
00370 
00371 const APPSPACK::Vector& APPSPACK::Constraints::Linear::getLower() const
00372 {
00373   return bLower;
00374 }
00375 
00376 const APPSPACK::Vector& APPSPACK::Constraints::Linear::getUpper() const
00377 {
00378   return bUpper;
00379 }
00380 
00381 const APPSPACK::Vector& APPSPACK::Constraints::Linear::getBhatLower() const
00382 {
00383   return bHatLower;
00384 }
00385 
00386 const APPSPACK::Vector& APPSPACK::Constraints::Linear::getBhatUpper() const
00387 {
00388   return bHatUpper;
00389 }
00390 
00391 const APPSPACK::Vector& APPSPACK::Constraints::Linear::getBtildeEq() const
00392 {
00393   return bTildeEq;
00394 }
00395 
00396 const APPSPACK::Matrix& APPSPACK::Constraints::Linear::getAhat() const
00397 {
00398   return aHat;
00399 }
00400 
00401 const APPSPACK::Matrix& APPSPACK::Constraints::Linear::getAtildeEq() const
00402 {
00403   return aTildeEq;
00404 }
00405 
00406 void APPSPACK::Constraints::Linear::getNominalX(APPSPACK::Vector& x) const
00407 {
00408   x.resize(scaling.size());
00409   for (int i = 0; i < scaling.size(); i++)
00410   {
00411     if (exists(bUpper[i]) && exists(bLower[i]))
00412       x[i] = (bUpper[i] + bLower[i]) / 2.0 ;
00413     else if (exists(bUpper[i]))
00414       x[i] = bUpper[i];
00415     else if (exists(bLower[i]))
00416       x[i] = bLower[i];
00417     else
00418       x[i] = 0;
00419   }
00420 }
00421 
00422 void APPSPACK::Constraints::Linear::scale(APPSPACK::Vector& x) const
00423 { 
00424   if (x.size() != scaling.size())
00425     error("scale", "x vector has incorrect length");
00426 
00427   for (int i = 0; i < scaling.size(); i++)
00428     x[i] = (x[i] - lHat[i]) / scaling[i];
00429 }
00430 
00431 void APPSPACK::Constraints::Linear::unscale(APPSPACK::Vector& w) const
00432 { 
00433   if (w.size() != scaling.size())
00434     error("scale", "w vector has incorrect length");
00435 
00436   for (int i = 0; i < scaling.size(); i++)
00437     w[i] = (scaling[i] * w[i]) + lHat[i];
00438 }
00439 
00440 bool APPSPACK::Constraints::Linear::isFeasible(const APPSPACK::Vector& x) const
00441 {  
00442   int n = scaling.size();
00443 
00444   // Error check
00445   if (x.size() != n)
00446     error("isFeasible", "x vector has incorrect length");
00447 
00448   // Create scaled x
00449   Vector xTilde(x);
00450   scale(xTilde);
00451 
00452   if (!isInequalityFeasible(xTilde))
00453     return false;
00454 
00455   if (!isEqualityFeasible(xTilde))
00456     return false;
00457   
00458   return true;
00459 }
00460 
00461 void APPSPACK::Constraints::Linear::assertFeasible(const APPSPACK::Vector& x) const
00462 {  
00463   if (!isFeasible(x))
00464     error("assertFeasible","Infeasible point");
00465 }
00466 
00467 bool APPSPACK::Constraints::Linear::isBoundsOnly() const
00468 {  
00469   return ((aIneq.empty()) && (aEq.empty()));
00470 }
00471 
00472 double  APPSPACK::Constraints::Linear::maxStep(const APPSPACK::Vector& x,
00473                                                const APPSPACK::Vector& d,
00474                                                double step) const
00475 {
00476   double maxStep = step;
00477   int n = scaling.size();
00478 
00479   Vector xTilde(x);
00480   scale(xTilde);
00481 
00482   // Determine max step for bounds component
00483 
00484   for (int i = 0; i < scaling.size(); i++)
00485   {
00486     if (d[i] < -scaling[i] * epsMach) // Points to lower bound
00487       switch (getIneqState(i, LowerBound, xTilde))
00488       {
00489       case Active:           // Can't move if the constraint is active
00490         return 0;
00491       case Inactive:            
00492         maxStep = min(maxStep, (bLower[i] - x[i]) / d[i]);
00493       }
00494     else if (d[i] > scaling[i] * epsMach) // Points to upper bound
00495       switch (getIneqState(i, UpperBound, xTilde))
00496       {
00497       case Active:           // Can't move if the constraint is active
00498         return 0;
00499       case Inactive:
00500         maxStep = min(maxStep, (bUpper[i] - x[i]) / d[i]);
00501       }
00502   }
00503 
00504   // Determine max step for inequality component.
00505 
00506   if (!aIneq.empty())
00507   {
00508     int p = aIneq.getNrows();
00509 
00510     // aix = aIneq * x
00511     Vector aix(p);
00512     aIneq.multVec(x, aix);
00513 
00514     // aid = aIneq * d
00515     Vector aid(p);
00516     aIneq.multVec(d, aid);
00517 
00518     for (int j = 0; j < p; j++)
00519     {
00520       if (aid[j] < -epsMach) // Points to lower bound
00521         switch (getIneqState(n + j, LowerBound, xTilde))
00522         {
00523         case Active:         // Can't move if the constraint is active
00524           return 0;
00525         case Inactive:
00526           maxStep = min(maxStep, (bIneqLower[j] - aix[j]) / aid[j]);
00527         }
00528       else if (aid[j] > epsMach)  // Points to upper bound
00529         switch (getIneqState(n + j, UpperBound, xTilde))
00530         {
00531         case Active:         // Can't move if the constraint is active
00532           return 0;
00533         case Inactive:
00534           maxStep = min(maxStep, (bIneqUpper[j] - aix[j]) / aid[j]);
00535         }
00536     }
00537   }
00538 
00539   // Process equality constraints.
00540 
00541   if (!aEq.empty())
00542   {
00543     int p = aEq.getNrows();
00544 
00545     // z = aEq * d
00546     Vector z(p);
00547     aEq.multVec(d,z);
00548 
00549     // Check that direction stays on hyperplane defined by the
00550     // equality constraints
00551     for (int i = 0; i < p; i++)
00552       if (fabs(z[i]) > epsMach)
00553         return 0;
00554   }
00555 
00556   return maxStep;
00557 }
00558 
00559 void APPSPACK::Constraints::Linear::getActiveIndex(const APPSPACK::Vector& xdist, 
00560                                                    double epsilon,
00561                                                    vector<ActiveType>& index) const
00562 
00563 {
00564   int nrows = aHat.getNrows();
00565   // Size the index
00566   index.resize(nrows);
00567 
00568   // Determine which bounds are near.
00569   for (int i = 0; i < nrows; i++)
00570   {
00571     // Check if upper or lower bounds are near.
00572     bool nearlow = (xdist[i] < epsilon);
00573     bool nearupp = (xdist[i+nrows] < epsilon);
00574     
00575     // Update index.
00576     if (nearlow && nearupp)
00577       index[i] = BothActive;
00578     else if (nearlow)           // and not nearupp
00579       index[i] = LowerActive;
00580     else if (nearupp)           // and not nearlow
00581       index[i] = UpperActive;
00582     else
00583       index[i] = NeitherActive;
00584   }
00585 }
00586 
00587 void APPSPACK::Constraints::Linear::formDistanceVector(const Vector& x, Vector& xdist) const
00588 {
00589   // Transform x to scaled space..
00590   Vector xTilde = x;
00591   scale(xTilde);
00592   
00593   int nrows = aHat.getNrows();
00594 
00595   // z = aHat*xTilde.
00596   Vector z(nrows);
00597   aHat.multVec(xTilde, z);
00598   
00599   // xdist must be twice the size of z to hold upper and lower bound info.
00600   xdist.resize(2*nrows);
00601   for (int i = 0; i < nrows; i++)
00602   {
00603     // Process lower bound.
00604     if (exists(bHatLower[i]))
00605     {
00606       if (aHatZNorm[i] > epsMach)
00607         xdist[i] = fabs( z[i]  - bHatLower[i] ) / aHatZNorm[i];
00608       else if (fabs( z[i] - bHatLower[i] ) < epsMach)
00609         xdist[i] = 0;
00610       else
00611         xdist[i] = dne();
00612     }
00613     else
00614       xdist[i] = dne();
00615     
00616     // Process upper bound.
00617     if (exists(bHatUpper[i]))
00618     {
00619       if (aHatZNorm[i] > epsMach)       
00620         xdist[i+nrows] = fabs( z[i]  - bHatUpper[i] ) / aHatZNorm[i];
00621       else if (fabs( z[i] - bHatUpper[i] ) < epsMach)
00622         xdist[i+nrows] = 0;
00623       else
00624         xdist[i+nrows] = dne();
00625     }
00626     else
00627       xdist[i+nrows] = dne();
00628   }
00629 }
00630 
00631 void APPSPACK::Constraints::Linear::snapToBoundary(APPSPACK::Vector& x, 
00632                                                    double esnap) const
00633 {
00634   // Form scaled x.
00635   Vector xTilde = x;
00636   scale(xTilde);
00637   // Form snap-to system.
00638   Matrix Asnap;
00639   Vector bsnap;
00640   formSnapSystem(xTilde, esnap, Asnap, bsnap);
00641   // Find closest point satisfying snap constraints.
00642   Asnap.constrainedLSQR(xTilde,bsnap);
00643   unscale(xTilde);
00644   // Update x.
00645   x = xTilde;
00646 }
00647 
00648 void APPSPACK::Constraints::Linear::makeBoundsFeasible(APPSPACK::Vector& x) const
00649 {
00650   for (int i = 0; i < x.size(); i++)
00651   {
00652     if (exists(bLower[i]) && x[i] < bLower[i])
00653       x[i] = bLower[i];
00654     if (exists(bUpper[i]) && x[i] > bUpper[i])
00655       x[i] = bUpper[i];
00656   }
00657 }
00658 
00659 void APPSPACK::Constraints::Linear::formSnapSystem(const APPSPACK::Vector& xTilde, double esnap,
00660                                                    APPSPACK::Matrix& Asnap,
00661                                                    APPSPACK::Vector& bsnap) const
00662 {
00663   Asnap.clear();
00664   bsnap.resize(0);
00665   
00666   if (aHat.empty())
00667   {
00668     if (!isEqualityFeasible(xTilde))
00669     {
00670       Asnap = aTildeEq;
00671       bsnap = bTildeEq;
00672     }
00673     return;
00674   }
00675 
00676   int nrow = aHat.getNrows();
00677   int nvar = aHat.getNcols();
00678   // Form z = Ahat*xTilde;
00679   Vector z(nrow);
00680   aHat.multVec(xTilde, z);
00681   
00682   // Measure and sort distance in terms of scaled space.
00683   multimap<double,int> dmap;
00684   multimap<double,int>::iterator iter;
00685   for (int i = 0; i < nrow; i++)
00686   { 
00687     if (exists(bHatLower[i]))
00688     { 
00689       double distlow = fabs( z[i] - bHatLower[i] ) / aHatNorm[i];
00690      dmap.insert(std::pair<const double, int>(distlow, i));
00691     }
00692     if (exists(bHatUpper[i]))
00693     { // Indices for upper bounds in dmap denoted by i+nrow.  Symbolically stacking matrix.
00694       double distupp = fabs( z[i] - bHatUpper[i] ) / aHatNorm[i];
00695       dmap.insert(std::pair<const double, int>(distupp, i+nrow));
00696     }
00697   }
00698 
00699   // Erase constraints which are outside the epsilon ball.
00700   dmap.erase(dmap.upper_bound(esnap),dmap.end());
00701 
00702   if ((dmap.lower_bound(epsMach) == dmap.end()) && isEqualityFeasible(xTilde))
00703     return; // There are no constraints within distance esnap which are inactive.
00704 
00705   // Equality constraints are always added.
00706   Asnap = aTildeEq;
00707   bsnap = bTildeEq;
00708 
00709   // Now add inequality constraints within distance esnap.
00710   // Recall that we have symbolically stacked aHat, so we must mod out by nrow.
00711   for (iter=dmap.begin(); iter != dmap.end(); iter++)
00712   {
00713     if (bsnap.size() >= nvar)
00714       break;
00715     int irow = iter->second;
00716     if (irow < nrow)
00717     { // Snap to lower bound.
00718       Asnap.addRow(aHat.getRow(irow));
00719       bsnap.push_back(bHatLower[irow]);
00720     }
00721     else
00722     { // Snap to upper bound.
00723       Asnap.addRow(aHat.getRow(irow % nrow));
00724       bsnap.push_back(bHatUpper[irow % nrow]);
00725     }
00726   }
00727 
00728   // Now ensure that rank(Asnap) = Asnap.getNrows().
00729   Asnap.pruneDependentRows(bsnap,epsMach);
00730 }
00731 
00732 void APPSPACK::Constraints::Linear::error(const string& fname, const string& msg) const
00733 {
00734   cout << "APPSPACK::Constraints::" << fname << " - " << msg << endl;
00735   throw "APPSPACK Error";
00736 }
00737 
00738 void APPSPACK::Constraints::Linear::print() const
00739 {
00740   cout << "\n";
00741   cout << "Bound Constraints " << endl;
00742   cout << "lower = " << bLower << endl;
00743   cout << "upper = " << bUpper << endl;
00744   cout << "scaling = " << scaling << endl;
00745   
00746   if (!aIneq.empty() && displayMatrix)
00747   {
00748     cout << endl << "Inequality Constraints" << endl;
00749     cout << "Inequality matrix = " << aIneq << endl;
00750     cout << "Inequality lower = " << bIneqLower << endl;
00751     cout << "Inequality upper = " << bIneqUpper << endl;
00752   }
00753 
00754   if (!aEq.empty() && displayMatrix)
00755   {
00756     cout << endl << "Equality Constraints" << endl;
00757     cout << "Equality matrix = " << aEq << endl;
00758     cout << "Equality bound = " << bEq << endl;
00759   }
00760 
00761   // Determine total number of finite variable bounds.
00762   int varlow = 0;
00763   int varupp = 0;
00764   for (int i = 0; i < bLower.size(); i++)
00765   {
00766     if (exists(bLower[i]))
00767       varlow++;
00768     if (exists(bUpper[i]))
00769       varupp++;
00770   }
00771 
00772   // Determine total number of (finite) inequality bounds.
00773   int genlow = 0;
00774   int genupp = 0;
00775   for (int i = 0; i < bIneqLower.size(); i++)
00776   {
00777     if (exists(bIneqLower[i]))
00778       genlow++;
00779     if (exists(bIneqUpper[i]))
00780       genupp++;
00781   }
00782 
00783   cout << "Problem size summary:" << endl;
00784   cout << "Number of variables: " << scaling.size() << endl;
00785   cout << "Number of variable bounds: " << varlow + varupp 
00786        << " (" << varlow << " lower, " << varupp << " upper)" << endl;
00787   cout << "Number of inequality constraints: " << genlow + genupp
00788        << " (" << genlow << " lower, " << genupp << " upper)" << endl;
00789   cout << "Number of equality constraints: " << bEq.size() << endl;
00790 
00791 }
00792 
00793 // PRIVATE
00794 APPSPACK::Constraints::StateType 
00795 APPSPACK::Constraints::Linear::getIneqState(int i, BoundType bType, const Vector& xTilde, double epsilon) const
00796 {
00797   if (epsilon == -1)
00798     epsilon = epsMach;
00799 
00800   // a = constraint row
00801   const Vector& a = aHat.getRow(i);
00802   double anorm = aHatNorm[i];
00803 
00804   // b = lhs
00805   double b = (bType == LowerBound) ? bHatLower[i] : bHatUpper[i];
00806 
00807   // Check finiteness of the bound
00808   if (!exists(b))
00809     return DNE;
00810 
00811   // z = a' * xTilde
00812   double z;
00813   z = xTilde.dot(a);
00814 
00815   // Check if the constraint is epsilon-active
00816   if ( fabs(z - b) < (epsilon * anorm) )
00817     return Active;
00818 
00819   // Check if the constraint is otherwise satisfied
00820   if (((bType == LowerBound) && ( b <= z )) ||
00821       ((bType == UpperBound) && ( z <= b )))
00822     return Inactive;
00823 
00824   // Otherwise, it must be violated.
00825   return Violated;
00826 }
00827 
00828 // PRIVATE
00829 APPSPACK::Constraints::StateType 
00830 APPSPACK::Constraints::Linear::getEqState(int i, const APPSPACK::Vector& xTilde, double epsilon) const
00831 {
00832   if (epsilon == -1)
00833     epsilon = epsMach;
00834 
00835   // a = constraint row
00836   const Vector& a = aTildeEq.getRow(i);
00837   double anorm = a.norm();
00838 
00839   // b = rhs
00840   double b = bTildeEq[i];
00841 
00842   // z = a' * xTilde
00843   double z;
00844   z = xTilde.dot(a);
00845 
00846   // Check if the constraint is epsilon-active
00847   if ( fabs(z - b) < (epsilon * anorm) )
00848     return Active;
00849 
00850   // Otherwise, it must be violated.
00851   return Violated;
00852 }
00853 
00854 // PRIVATE
00855 bool APPSPACK::Constraints::Linear::isEqualityFeasible(const APPSPACK::Vector& xTilde) const
00856 {  
00857   // Check equality constraints
00858   for (int i = 0; i < aTildeEq.getNrows(); i ++)
00859   {
00860     if (Violated == getEqState(i, xTilde))
00861       return false;
00862   }
00863   
00864   return true;
00865 }
00866 
00867 // PRIVATE
00868 bool APPSPACK::Constraints::Linear::isInequalityFeasible(const APPSPACK::Vector& xTilde) const
00869 {  
00870   // Check inequality constraints
00871   for (int i = 0; i < aHat.getNrows(); i ++)
00872   {
00873     if (Violated == getIneqState(i, UpperBound, xTilde))
00874       return false;
00875     if (Violated == getIneqState(i, LowerBound, xTilde))
00876       return false;
00877   }
00878   
00879   return true;
00880 }

 

© Sandia Corporation | Site Contact | Privacy and Security

Generated on Fri Feb 16 10:33:34 2007 for APPSPACK 5.0.1 by doxygen 1.3.9.1 written by Dimitri van Heesch, © 1997-2002