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_Directions.cpp

Go to the documentation of this file.
00001 // $Id: APPSPACK_Directions.cpp,v 1.90.2.1 2007/01/14 20:45:01 jgriffi Exp $ 
00002 // $Source: /space/CVS-Acro/acro/packages/appspack/appspack/src/APPSPACK_Directions.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_Directions.hpp"
00039 #include "APPSPACK_Print.hpp"
00040 #include "APPSPACK_Utils.hpp"
00041 #include "APPSPACK_Constraints_Linear.hpp"
00042 #include "APPSPACK_Float.hpp"
00043 #include "APPSPACK_CDDLIB_Interface.h"
00044 #include "APPSPACK_LAPACK_Wrappers.hpp"
00045 
00046 APPSPACK::Directions::Directions(Parameter::List& params, const Constraints::Linear& constraints_in) :
00047   constraints(constraints_in),
00048   nDimensions(constraints_in.getScaling().size()),
00049   zero(nDimensions, static_cast<double>(0)),
00050   stepTolerance(params.getParameter("Step Tolerance", 0.01)),
00051   minStep(params.getParameter("Minimum Step", 2 * stepTolerance)),
00052   theta(params.getParameter("Contraction Factor", 0.5)),
00053   epsilonMax(params.getParameter("Epsilon Max", 2 * stepTolerance)),
00054   withNormals(params.getParameter("Add Projected Normals", true)),
00055   withCompass(params.getParameter("Add Projected Compass", false)),
00056   nDirections(0),
00057   nCached(0),
00058   nLapack(0),
00059   nCddlib(0),
00060   nMaxDirections(0),
00061   nAppend(0)
00062 {
00063   // Check parameters
00064   if (stepTolerance <= 0)
00065   {
00066     cout << "APPSPACK::Directions::Directions - Error: \"Step Tolerance\" cannot be negative." << endl;
00067     throw "APPSACK Error";
00068   }
00069 
00070   if (minStep <= stepTolerance)
00071   {
00072     cout << "APPSPACK::Directions::Directions - Error: \"Minimum Step\" must be greater than \"Step Tolerance\"." << endl;
00073     throw "APPSACK Error";
00074   }
00075 
00076   if ((theta <= 0) || (theta >= 1))
00077   {
00078     cout << "APPSPACK::Directions::Directions - Error: \"Contraction Factor\" must be strictly between zero and one." << endl;
00079     throw "APPSACK Error";
00080   }
00081 
00082   epsilonMin = epsilonMax;
00083 }
00084 
00085 APPSPACK::Directions::~Directions()
00086 {
00087 }
00088 
00089 const APPSPACK::Vector& APPSPACK::Directions::getDirection(int i) const
00090 {
00091   return direction.getRow(i);
00092 }
00093 
00094 double APPSPACK::Directions::getStep(int i) const
00095 {
00096   return step[i];
00097 }
00098 
00099 void APPSPACK::Directions::setStepConverged(int i)
00100 {
00101   step[i] = stepTolerance/2;
00102 }
00103 
00104 void APPSPACK::Directions::getDirectionIndices(vector<int>& idvector) const
00105 {
00106   idvector.resize(0);
00107   for (int i = 0; i < nDirections; i ++)
00108     if ((step[i] >= stepTolerance) && (tag[i] == -1))
00109       idvector.push_back(i);
00110 }
00111 
00112 void APPSPACK::Directions::computeNewDirections(const Point& newPoint)
00113 {
00114   // Grab new point info.
00115   double newStep = max(newPoint.getStep(), minStep);
00116   const Vector& x = newPoint.getX();
00117   
00118   // Update distance vector for new point.
00119   constraints.formDistanceVector( x, xDistance );
00120 
00121   do 
00122   {
00123     // Update active state vector, member constraintState.
00124     if (updateConstraintState(newStep))
00125     {    
00126       // Generate search directions.
00127       generateForLinear(direction);
00128     }
00129 
00130     if (!direction.empty())
00131     {
00132       // Update step, trueStep, epsilonMin, and nDirections.
00133       updateDirectionInfo(newStep);
00134       break;
00135     }
00136 
00137     // If directions are empty, we will want a smaller value of newStep.
00138     newStep = theta * newStep;
00139   } while (newStep >= stepTolerance);
00140   
00141   if (direction.empty())
00142   {
00143     cout << "Error in APPSPACK::Directions::computeNewDirections().\n"
00144          << "Unable to generate generators for epsilon-tangent cone.\n"
00145          << "One of two possible problems may be the culprit:\n"
00146          << "(1) Parameter \"Step Tolerance\" is too large \n"
00147          << "(2) No feasible search directions exist at current point. \n"
00148          << endl;
00149     throw "APPSPACK Error";
00150   }
00151 }
00152 
00153 void APPSPACK::Directions::appendNewDirections()
00154 {
00155   double newEps = getSmallestStep();
00156   // If newEps >= epsilonMin active set remains the same.
00157   if (newEps >= epsilonMin)
00158     return;
00159 
00160   // If constraintState remains the same, do nothing.
00161   if (!updateConstraintState(newEps))
00162   {
00163     epsilonMin = newEps;
00164     return;
00165   }
00166 
00167   // Compute new directions.  Note that distance vector has not been updated
00168   // as the current point is assumed the same.
00169   Matrix newDirection;
00170   generateForLinear(newDirection);
00171 
00172   // Now add new directions to the current list.
00173   direction.addUniqueRows(newDirection, constraints.getEpsMach());
00174   
00175   // Update step, trueStep, espilonk, and nDirections.
00176   updateDirectionInfo(newEps, true);
00177 }
00178 
00179 void APPSPACK::Directions::setTrueStepAndTag(int i, double trueStep_in, int tag_in)
00180 {
00181   trueStep[i] = trueStep_in;
00182   tag[i] = tag_in;
00183 }
00184 
00185 
00186 void APPSPACK::Directions::print(const string label) const
00187 {
00188   if (!label.empty())
00189     cout << label << ":\n";
00190 
00191   for (int i = 0; i < nDirections; i ++)
00192   {
00193     cout << setw(4) << i << " : ";
00194 
00195     cout << "d = " << direction.getRow(i) << " ";
00196 
00197     cout<< "step = " << Print::formatPositiveDouble(step[i]) << " ";
00198 
00199     if (tag[i] != -1) 
00200     {
00201       cout << "tag = " << setw(6) << tag[i] << " ";
00202       cout << "trueStep = " <<  Print::formatPositiveDouble(trueStep[i]);
00203     }
00204 
00205     cout << "\n";
00206   }
00207   cout << "Number of times directions calculated by...\n"
00208        << "  LAPACK: "    << nLapack << endl
00209        << "  CDDLIB: "    << nCddlib << endl
00210        << "  Cached: " << nCached  << endl
00211        << "Max directions in single iteration : " << nMaxDirections << endl
00212        << "Number of times directions appended: " << nAppend << endl;
00213 }
00214 
00215 bool APPSPACK::Directions::isStepConverged() const
00216 {
00217   for (int i = 0; i < nDirections; i ++)
00218   {
00219     if (step[i] >= stepTolerance)
00220       return false;
00221   }
00222 
00223   return true;
00224 }
00225 
00226 bool APPSPACK::Directions::empty() const
00227 {
00228   return (direction.empty() == true);
00229 }
00230 
00231 void APPSPACK::Directions::reduceStep(int i)
00232 {
00233   double tmpStep = theta * step[i];
00234 
00235   step[i] = tmpStep;
00236   trueStep[i] = -1;
00237   tag[i] = -1;
00238 }
00239 
00240 //PRIVATE
00241 void APPSPACK::Directions::buildNormalCone(APPSPACK::Matrix& VpT,
00242                                                    APPSPACK::Matrix& VlT) const
00243 {
00244   // Always insert equality constraints in VlT.
00245   VlT.addMatrix(constraints.getAtildeEq());
00246   
00247   // Insert elements into VpT and VlT from aHat.
00248   const Matrix& aHat = constraints.getAhat();
00249   for (int i = 0; i < constraintState.size(); i++)
00250   {
00251     if (constraintState[i] == APPSPACK::Constraints::BothActive)
00252       VlT.addRow(aHat.getRow(i));
00253     else if (constraintState[i] == APPSPACK::Constraints::LowerActive)
00254       VpT.addRow(aHat.getRow(i), -1.0);
00255     else if (constraintState[i] == APPSPACK::Constraints::UpperActive)
00256       VpT.addRow(aHat.getRow(i));
00257   }
00258 }
00259 
00260 //PRIVATE
00261 void APPSPACK::Directions::buildTangentCone(const APPSPACK::Matrix& VpT, 
00262                                               const APPSPACK::Matrix& VlT,
00263                                               APPSPACK::Matrix& T)
00264 {
00265 #if !(defined HAVE_LAPACK) && !(defined HAVE_CDDLIB)
00266   buildWithNothing(T);
00267   return;
00268 #endif
00269   
00270   if ((VpT.empty()) && (VlT.empty()))
00271   {
00272     // Generate all compass directions. Problem looks locally unconstrained.
00273     generateUnconstrained(T);
00274     return;
00275   }
00276   
00277   if (buildWithLapack(VpT, VlT, T))
00278   {
00279     nLapack++;
00280     return;
00281   }
00282   
00283   // Either constraints degenerate or LAPACK unavailable.  Try cddlib.
00284   if (buildWithCddlib(VpT, VlT, T))
00285   {
00286     nCddlib++;
00287     return;
00288   }
00289   
00290   cout << "Error in APPSPACK::Directions::buildTangentCone().\n"
00291        << "Unable to generate generators for epsilon-tangent cone.\n"
00292        << "One of four possible problems may be the culprit:\n"
00293        << "(1) APPSPACK was configured without LAPACK\n"
00294        << "(2) APPSPACK was configured without CDDLIB and the constraints are degenerate\n"
00295        << "(3) Either CDDLIB or LAPACK capability disabled in .apps files\n"
00296        << "(4) CDDLIB failed to find generators (please submit bug report!)\n"
00297        << endl;
00298   throw "APPSPACK Error";
00299 }
00300 
00301 //PRIVATE
00302 void APPSPACK::Directions::buildWithNothing(APPSPACK::Matrix& D)                                             
00303 {
00304   D.clear();
00305   const Vector& scaling = constraints.getScaling();
00306   for (int i = 0; i < constraintState.size(); i++)
00307   {
00308     tmpVector = zero;
00309     if (constraintState[i] == APPSPACK::Constraints::NeitherActive)
00310     {
00311       // Add +e_i
00312       tmpVector[i] = scaling[i];
00313       D.addRow(tmpVector);
00314       // Add -e_i
00315       tmpVector[i] = -1 * scaling[i];
00316       D.addRow(tmpVector);
00317     }
00318     else if (constraintState[i] == APPSPACK::Constraints::LowerActive)
00319     {
00320       // Add +e_i
00321       tmpVector[i] = scaling[i];
00322       D.addRow(tmpVector);
00323     }
00324     else if (constraintState[i] == APPSPACK::Constraints::UpperActive)
00325     {
00326       // Add -e_i
00327       tmpVector[i] = -1 * scaling[i];
00328       D.addRow(tmpVector);
00329     }
00330   }
00331 }
00332 
00333 //PRIVATE
00334 bool APPSPACK::Directions::buildWithCddlib( const APPSPACK::Matrix& VpT,
00335                                                const APPSPACK::Matrix& VlT,
00336                                                APPSPACK::Matrix& Tcone)
00337 {
00338 #ifndef HAVE_CDDLIB
00339   return false;
00340 #endif
00341   Tcone.clear();  
00342   const Vector& scaling = constraints.getScaling(); 
00343   int nvar = scaling.size();
00344 
00345   // Get row pointers needed by call to compute_cone_generators.
00346   vector< double *> VpTptr;
00347   const_cast<Matrix &>(VpT).getRowPointers(VpTptr);
00348 
00349   vector< double *> VlTptr;
00350   const_cast<Matrix &>(VlT).getRowPointers(VlTptr);
00351   
00352   int num_pointy = 0;
00353   int num_lineality = 0;
00354   double **P;
00355   double **L;
00356 
00357   compute_cone_generators( &num_pointy, &P, &num_lineality, &L, nvar,
00358                            VlT.getNrows(), &VlTptr[0], VpT.getNrows(), &VpTptr[0], 0);
00359 
00360   // Add in lineality directions.
00361   Matrix Lmat(L, num_lineality, nvar);
00362   Lmat.normalize();
00363   Lmat.scale(scaling);
00364   Tcone.addMatrix(Lmat);
00365   Tcone.addMatrix(Lmat, -1.0);
00366 
00367   // Add in pointy directions.
00368   Matrix Pmat(P, num_pointy, nvar);
00369   Pmat.normalize();
00370   Pmat.scale(scaling);
00371   Tcone.addMatrix(Pmat);
00372   
00373   // Deallocate.
00374   for (int i = 0; i < num_pointy; i++)
00375     free(P[i]);
00376   free(P);
00377   for (int i = 0; i < num_lineality; i++)
00378     free(L[i]);
00379   free(L);
00380   
00381   return true;
00382 }
00383 
00384 //PRIVATE
00385 bool APPSPACK::Directions::buildWithLapack(const APPSPACK::Matrix& VpT, 
00386                                               const APPSPACK::Matrix& VlT,
00387                                               APPSPACK::Matrix& Tcone)
00388 { 
00389 #ifndef HAVE_LAPACK
00390   return false;
00391 #endif
00392   Tcone.clear();
00393   // Get orthonormal basis matrix Z for nullspace of VlT if Vlt is non-empty.
00394   // Otherwise Z becomes the idenity.
00395   Matrix ZT;
00396   if (!VlT.empty())
00397   {
00398     VlT.nullSpace(ZT, constraints.getEpsMach());
00399     if (ZT.empty())
00400       return true;
00401   }
00402   else
00403   {
00404     // The are no equality constraints.  null(VlT) = I (trivially).
00405     ZT.setToIdentity(VpT.getNcols());
00406   }
00407 
00408   // Only equality constraints are present, so +/- Z generates space.
00409   if (VpT.empty()) 
00410   {
00411     // ZTscaled = Z^T * S
00412     Matrix ZTscaled(ZT, constraints.getScaling());
00413     Tcone.addMatrix(ZTscaled);
00414     Tcone.addMatrix(ZTscaled, -1.0);
00415     return true;
00416   }
00417   
00418   // Compute product VpTZ = VpT * Z.
00419   Matrix VpTZ;
00420   VpT.multMat(ZT, VpTZ, Matrix::Transpose);
00421   
00422   // Compute matrices R and N such that VpTZ*R = I and VptZ*N = 0;
00423   Matrix RT, NT;
00424   if (!VpTZ.getRighInvAndNullBasis(RT, NT, constraints.getEpsMach()))
00425     return false;
00426 
00427   // Now form matrix Z*R.  Since we are working with transposes, we really form RT*ZT.
00428   Matrix ZRT(RT);
00429   ZRT.multMat(ZT);
00430   // Columns of Z*R are not necessarily unit length.
00431   ZRT.normalize();
00432   // Now scale columns and add to direction.
00433   ZRT.scale(constraints.getScaling());
00434   Tcone.addMatrix(ZRT, -1.0);
00435 
00436   // Now form Z*N.  Since we are working with transposes, we really form NT*ZT.
00437   // N may be empty if Z'*Vp is a square system, implying null(Z'*Vp) = {0}.
00438   if (!NT.empty())
00439   {
00440     Matrix ZNT(NT);
00441     ZNT.multMat(ZT); 
00442     // Normalizing unnecessary, N and Z both have orthonormal columns.
00443     // Now scale columns and add to direction.
00444     ZNT.scale(constraints.getScaling());
00445     Tcone.addMatrix(ZNT);
00446     Tcone.addMatrix(ZNT, -1.0);
00447   }
00448   
00449   return true;
00450 }
00451 
00452 //PRIVATE
00453 void APPSPACK::Directions::generateUnconstrained(APPSPACK::Matrix& D)                                                
00454 {
00455   D.clear();
00456   const Vector& scaling = constraints.getScaling();
00457 
00458   for (int i = 0; i < nDimensions; i ++)
00459   {
00460     tmpVector = zero;    
00461     // Add +e_i
00462     tmpVector[i] = scaling[i];
00463     D.addRow(tmpVector);
00464     // Add -e_i
00465     tmpVector[i] = -1 * scaling[i];
00466     D.addRow(tmpVector);
00467   }
00468 }
00469 
00470 //PRIVATE
00471 void APPSPACK::Directions::generateForLinear(APPSPACK::Matrix& D)
00472 {
00473   D.clear();
00474 
00475   // Check if directions are cached.
00476   cacheIter = directionCache.find(constraintState);
00477   if (cacheIter != directionCache.end())
00478   {
00479     nCached++;
00480     D = cacheIter->second;
00481     return;
00482   }
00483   
00484   // Directions are not cached.  Form the normal cone.
00485   Matrix VpT;
00486   Matrix VlT;
00487   buildNormalCone(VpT, VlT);
00488   buildTangentCone(VpT, VlT, D);
00489   
00490   // If there are no tangent direction, Deltak must be reduced.
00491   if (D.empty())
00492   {
00493     // We still want to cache direction here--the work to figure out D is 
00494     // empty may be nontrivial.
00495     directionCache[constraintState] = D;
00496     return;
00497   }
00498   
00499   // Add normals.
00500   if (withNormals)
00501     addNormalDirections(VpT, VlT, D);
00502 
00503   if (withCompass)
00504     addCompassDirections(VlT, D);
00505   
00506   // Cache newly computed directions.
00507   directionCache[constraintState] = D;
00508 }
00509 
00510 //PRIVATE
00511 void APPSPACK::Directions::addNormalDirections(const APPSPACK::Matrix& VpT, 
00512                                                const APPSPACK::Matrix& VlT,
00513                                                APPSPACK::Matrix& D)
00514 {
00515   // Add in normal directions, projecting with Z*ZT if necessary.
00516   if (VpT.empty())
00517     return;
00518   
00519   Matrix ND(VpT);
00520   if (!VlT.empty())
00521   {    
00522 #ifndef HAVE_LAPACK
00523     {
00524       cerr << "APPSPACK::Directions::addNormalDirections - Error: "
00525            << "Cannot add projected normals without LAPACK." << endl;
00526       throw "APPSACK Error";
00527     }
00528 #endif
00529     Matrix ZT;
00530     VlT.nullSpace(ZT, constraints.getEpsMach());
00531     if (ZT.empty()) // null(ZT) is empty, projecting to empty set.
00532       return;
00533     ND.multMat(ZT, Matrix::Transpose);
00534     ND.multMat(ZT);
00535   }
00536   
00537   ND.normalize();
00538   ND.scale(constraints.getScaling());
00539   D.addMatrix(ND);
00540 }
00541 
00542 //PRIVATE
00543 void APPSPACK::Directions::addCompassDirections(const APPSPACK::Matrix& VlT,
00544                                                 APPSPACK::Matrix& D)
00545 {
00546   //First form compass search directions.
00547   Matrix PCD;
00548   generateUnconstrained(D);
00549   //No project if necessary
00550   if (!VlT.empty())
00551   {    
00552 #ifndef HAVE_LAPACK
00553     {
00554       cerr << "APPSPACK::Directions::addNormalDirections - Error: "
00555            << "Cannot add projected normals without LAPACK." << endl;
00556       throw "APPSACK Error";
00557     }
00558 #endif
00559     Matrix ZT;
00560     VlT.nullSpace(ZT, constraints.getEpsMach());
00561     if (ZT.empty()) // null(ZT) is empty, projecting to empty set.
00562       return;
00563     PCD.multMat(ZT, Matrix::Transpose);
00564     PCD.multMat(ZT);
00565   }
00566   
00567   PCD.normalize();
00568   PCD.scale(constraints.getScaling());
00569   D.addUniqueRows(PCD, constraints.getEpsMach());
00570 }
00571 
00572 //PRIVATE
00573 bool APPSPACK::Directions::updateConstraintState(double newStep)                                             
00574 {
00575   
00576   // Ensure epsilon < epsilonMax.
00577   double epsilon = min(epsilonMax, newStep);
00578   // Determine indices of near active constraints wrt epsilon.
00579   vector<Constraints::ActiveType> newState;
00580   constraints.getActiveIndex(xDistance, epsilon, newState);
00581 
00582   // If constraintState remains the same and is not updated return false.
00583   if (newState == constraintState)
00584     return false;
00585   
00586   // Otherwise update constraintState to newState and return true.
00587   constraintState = newState;
00588   return true;
00589 }
00590 
00591 //PRIVATE
00592 void APPSPACK::Directions::updateDirectionInfo(double newStep, bool isAppend)
00593 {
00594   //Note: newStep \ge minStep if 1) isAppend = false, and 2) the corresponding
00595   //tangent cone is empty.  newStep is bumped up to minStep initially inside
00596   //computeNewDirections.  An empty tangent cone corresponds to a failed iteration
00597   //and hence we do not want to bumb up to minStep in this case.
00598   if (isAppend)
00599   {
00600     // Only update new directions.
00601     int numnew = direction.getNrows() - nDirections;
00602     if (numnew > 0)
00603       nAppend++;
00604     nDirections = direction.getNrows();
00605     step.append(numnew, newStep);
00606     trueStep.append(numnew, -1);
00607     tag.insert(tag.end(), numnew, -1);
00608   }
00609   else
00610   {
00611     // Update all directions.
00612     nDirections = direction.getNrows();
00613     step.assign(nDirections, newStep);
00614     trueStep.assign(nDirections, -1);
00615     tag.assign(nDirections, -1);
00616   }
00617   nMaxDirections=max(nMaxDirections, nDirections);
00618   // epsilonMin records smallest value of epsilon used so far.
00619   epsilonMin = min(epsilonMax, getSmallestStep());
00620 }
00621 
00622 double APPSPACK::Directions::getSmallestStep() const
00623 {
00624   double sstep = step.max();
00625   for (int i=0; i < step.size(); i++)
00626   {
00627     if (step[i] >= stepTolerance)
00628       sstep = min(sstep, step[i]);
00629   }
00630   return sstep;
00631 }
00632 
00633 

 

© Sandia Corporation | Site Contact | Privacy and Security

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