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

Go to the documentation of this file.
00001 // $Id: APPSPACK_Solver.cpp,v 1.65.2.3 2007/02/03 23:37:05 jgriffi Exp $ 
00002 // $Source: /space/CVS-Acro/acro/packages/appspack/appspack/src/APPSPACK_Solver.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 <algorithm>
00039 #include "APPSPACK_Solver.hpp"
00040 #include "APPSPACK_Print.hpp"
00041 #include "APPSPACK_Constraints_Linear.hpp"
00042 #include "APPSPACK_Float.hpp"
00043 #include "APPSPACK_LAPACK_Wrappers.hpp"
00044 
00045 APPSPACK::Solver::Solver(const Parameter::List& params_in,
00046                          Executor::Interface& executor_in,
00047                          const Constraints::Linear& constraints_in,
00048                          Combiner::Generic& combiner_in) :
00049   constraints(constraints_in),
00050   params(params_in),
00051   print(params),                // modifies params and sets globals
00052   bestPointPtr(initializeBestPointPtr(params, constraints, combiner_in)), // modifies params
00053   directions(params, constraints), // modifies params
00054   conveyor(params, constraints.getScaling(), executor_in), // modifies params
00055   exchangeList(),
00056   state(Continue),
00057   isFunctionTolerance(false),
00058   isMaxEvaluations(false),
00059   useSnapTo(params.getParameter("Snap To Boundary",true))
00060 {
00061   setup(params_in, executor_in, constraints_in);
00062 }
00063 
00064 APPSPACK::Solver::Solver(const Parameter::List& params_in,
00065                          Executor::Interface& executor_in,
00066                          const Constraints::Linear& constraints_in) :
00067   constraints(constraints_in),
00068   params(params_in),
00069   print(params),                // modifies params and sets globals
00070   bestPointPtr(initializeBestPointPtr(params, constraints, combiner)), // modifies params
00071   directions(params, constraints), // modifies params
00072   conveyor(params, constraints.getScaling(), executor_in), // modifies params
00073   exchangeList(),
00074   state(Continue),
00075   isFunctionTolerance(false),
00076   isMaxEvaluations(false),
00077   useSnapTo(params.getParameter("Snap To Boundary",true))
00078 {
00079   setup(params_in, executor_in, constraints_in);
00080 }
00081 
00082 void APPSPACK::Solver::setup(const Parameter::List& params_in,
00083                              Executor::Interface& executor_in,
00084                              const Constraints::Linear& constraints_in)
00085 {
00086   useRandomOrder = params.getParameter("Use Random Order", false);
00087 
00088   cout << "\n"
00089        << "-----------------------------------------------------\n"
00090        << "APPSPACK: Asynchronous Parallel Pattern Search\n"
00091        << "Written by T. G. Kolda et al., Sandia National Labs\n"
00092        << "For more information visit \n"
00093        << "http://software.sandia.gov/appspack\n"
00094        << "-----------------------------------------------------\n"
00095        << "\n";
00096   
00097   // Stopping Criteria
00098   if (params.isParameterDouble("Function Tolerance"))
00099   {
00100     isFunctionTolerance = true;
00101     functionTolerance = params.getDoubleParameter("Function Tolerance");
00102   }
00103 
00104   if (params.isParameter("Maximum Evaluations"))
00105   {
00106     isMaxEvaluations = true;
00107     maxEvaluations = params.getParameter("Maximum Evaluations", 0);
00108   }
00109 
00110   boundsTolerance = params.getParameter("Bounds Tolerance", params.getDoubleParameter("Step Tolerance") / 2.0);
00111   if (boundsTolerance <= 0)
00112   {
00113     cout << "APPSPACK::Solver::Solver - Invalid non-positive value for \"Bounds Tolerance\"." << endl;
00114     throw "APPSPACK Error";
00115   }
00116 
00118   if (Print::doPrint(Print::InitialData))
00119     printInitializationInformation();
00120   
00121   processNewBestPoint();
00122 }
00123 
00124 APPSPACK::Solver::~Solver()
00125 {
00126   delete bestPointPtr;
00127 }
00128 
00129 const vector<double>& APPSPACK::Solver::getBestX() const
00130 {
00131   return bestPointPtr->getX().getStlVector();
00132 }
00133 
00134 double APPSPACK::Solver::getBestF() const
00135 {
00136   return bestPointPtr->getF(); 
00137 }
00138 
00139 const vector<double>& APPSPACK::Solver::getBestVecF() const
00140 {
00141   return bestPointPtr->getVecF().getStlVector(); 
00142 }
00143 
00144 APPSPACK::Point* APPSPACK::Solver::initializeBestPointPtr(Parameter::List& params, 
00145                                                           const Constraints::Linear& constraints,
00146                                                           Combiner::Generic& combiner)
00147 {
00148   // Set up the initial point
00149   Vector nominalX;
00150   if (!params.isParameter("Initial X"))
00151     constraints.getNominalX(nominalX);
00152   
00153   Vector initialX = params.getParameter("Initial X", nominalX);
00154   
00155   // Check that the initial point is the right size
00156   if (initialX.size() != constraints.getScaling().size())
00157   {
00158     cerr << "Error: Size mismatch\n";
00159     cerr << "The size of the initial X is not the same as the lower bounds.\n";
00160     throw "APPSPACK Error";
00161   }  
00162 
00163   // Try making feasible wrt to simple bounds.
00164   if (!constraints.isFeasible(initialX))
00165     constraints.makeBoundsFeasible(initialX);
00166   
00167 #ifdef HAVE_LAPACK
00168   double epsilonSnap = params.getParameter("Bounds Tolerance", params.getDoubleParameter("Step Tolerance") / 2.0);
00169   if (epsilonSnap <= 0)
00170   {
00171     cout << "APPSPACK::Solver::Solver - Invalid non-positive value for \"Bounds Tolerance\"." << endl;
00172     throw "APPSPACK Error";
00173   }
00174 
00175   // Try snapping to
00176   if (!constraints.isFeasible(initialX)) 
00177     constraints.snapToBoundary(initialX,epsilonSnap);
00178 #endif
00179 
00180   // Check that x is feasible
00181   constraints.assertFeasible(initialX);
00182 
00183   // Initial F
00184   Vector initialF;
00185   if (params.isParameterDouble("Initial F"))
00186   {
00187     initialF.assign(1,params.getDoubleParameter("Initial F"));
00188     params.setParameter("Initial F",initialF); // reset to be a vector!
00189   }
00190   else if (params.isParameterVector("Initial F"))
00191   {
00192     initialF = params.getVectorParameter("Initial F");
00193   }
00194 
00195   // Initial step
00196   double initialStep = params.getParameter("Initial Step", 1.0);
00197   
00198   // Point
00199   double alpha = params.getParameter("Sufficient Decrease Factor", 0.01);
00200   
00201   // Create first best point
00202   return new Point(initialX, initialF, initialStep, alpha, combiner);
00203 }
00204 
00205 APPSPACK::Solver::State APPSPACK::Solver::solve()
00206 {
00207   while (state == Continue)
00208   {
00209     iterate();
00210   }
00211 
00212   // Print the final solution
00213   if (Print::doPrint(Print::FinalSolution))
00214   {
00215     cout << "\nFinal State: " << state << "\n";
00216     printBestPoint("Final Min");
00217     conveyor.getCounter().print();
00218     directions.print("Final Directions");
00219     cout << "\n"
00220          << "-----------------------------------------------------\n"
00221          << "Thank you!\n"
00222          << "APPSPACK: Asynchronous Parallel Pattern Search\n"
00223          << "Written by T. G. Kolda et al., Sandia National Labs\n"
00224          << "For more information visit \n"
00225          << "http://software.sandia.gov/appspack\n"
00226          << "-----------------------------------------------------\n"
00227          << "\n";
00228   }
00229 
00230 
00231   
00232   // Print solution file if specified by user.
00233   string solfile;
00234   solfile = params.getParameter("Solution File", solfile);
00235   if (!solfile.empty())
00236     writeSolutionFile(solfile);
00237 
00238   return state;
00239 }
00240 
00241 APPSPACK::Solver::State APPSPACK::Solver::iterate()
00242 {
00243     generateTrialPoints();
00244     
00245     conveyor.exchange(exchangeList);
00246 
00247     processEvaluatedTrialPoints();
00248     
00249     return state;
00250 }
00251 
00252 // PRIVATE
00253 void APPSPACK::Solver::processNewBestPoint(Point* newBestPointPtr)
00254 {
00255   // Update the best point
00256   if (newBestPointPtr != NULL)
00257   {
00258     delete bestPointPtr;
00259     bestPointPtr = newBestPointPtr;
00260   }
00261 
00262   // Print
00263   if (Print::doPrint(Print::NewBestPoint))
00264     printBestPoint("New Min");
00265   
00266   // Check for convergence based on the function value
00267   if ((isFunctionTolerance) && (bestPointPtr->getF() < functionTolerance))
00268   {
00269     state = FunctionConverged;
00270     return;
00271   }
00272 
00273   // Update the (scaled) search directions
00274   directions.computeNewDirections(*bestPointPtr);
00275 
00276   // Print
00277   if (Print::doPrint(Print::NewBestDirections))
00278     directions.print("New Directions");
00279 }
00280 
00281 // PRIVATE
00282 void APPSPACK::Solver::generateTrialPoints()
00283 {
00284   // Local references
00285   const Vector& parentX = bestPointPtr->getX();
00286   
00287   Vector& x = tmpVector;
00288   int n = parentX.size();
00289   x.resize(n);
00290 
00291   vector<int> indices;
00292   directions.getDirectionIndices(indices);
00293 
00294   if (useRandomOrder == true)
00295   {
00296     // Randomize access order.    
00297     random_shuffle(indices.begin(), indices.end());
00298   }
00299   
00300   for (int i = 0; i < indices.size(); i ++)
00301   {
00302     int idx = indices[i];
00303     const Vector& direction = directions.getDirection(idx);
00304     const double step = directions.getStep(idx);
00305     double maxStep = constraints.maxStep(parentX, direction, step);
00306   
00307     // Only add nontrivial points to the exchange list.
00308     if (maxStep <= 0)
00309     {
00310       // Cannot move along this direction and remain feasible.
00311       directions.setStepConverged(idx);
00312     }
00313     else
00314     {
00315       for (int j = 0; j < n; j++)
00316         x[j] = parentX[j] + maxStep*direction[j];
00317       
00318 #ifdef HAVE_LAPACK
00319       if (useSnapTo)
00320       {
00321         Vector xsnap = x;
00322         constraints.snapToBoundary(xsnap,boundsTolerance);
00323         if (constraints.isFeasible(xsnap))
00324           x = xsnap;
00325       }
00326 #endif
00327       
00328       constraints.assertFeasible(x);
00329       
00330       // Create a new trial point
00331       Point* newPointPtr = new Point(x, step, *bestPointPtr, idx);
00332       
00333       // Save off trial point information
00334       directions.setTrueStepAndTag(idx, maxStep, newPointPtr->getTag());
00335       
00336       // Push this trial point onto the new trial point list.
00337       // Ownership of the pointer transfers to the list.
00338       exchangeList.push(newPointPtr);
00339     }
00340   }
00341   
00342   if (Print::doPrint(Print::Directions))
00343     directions.print("Directions After Trial Point Generation");
00344 
00345   if (Print::doPrint(Print::UnevaluatedPoints))
00346     exchangeList.print("Trial Points Before Conveyor");
00347 
00348 }
00349 
00350 void APPSPACK::Solver::processEvaluatedTrialPoints()
00351 {
00352   // Print
00353   if (Print::doPrint(Print::EvaluatedPoints))
00354     exchangeList.print("Trial Points After Conveyor");
00355 
00356   // Check for a new best point
00357   if (exchangeList.best() < *bestPointPtr)
00358   {
00359     processNewBestPoint(exchangeList.popBest());
00360     exchangeList.prune();
00361     conveyor.prune();
00362   }
00363   else 
00364   {
00365     // Otherwise, just process the list
00366     Point* ptr;
00367     bool stepReduced = false;
00368     while (!exchangeList.isEmpty())
00369     {
00370       ptr = exchangeList.pop();
00371       if (ptr->getParentTag() == bestPointPtr->getTag())
00372       {
00373         directions.reduceStep(ptr->getIndex());
00374         stepReduced = true;
00375       }
00376       delete ptr;
00377     }
00378     // Check for step length convergence.  If not converged, append new directions to list.
00379     if (directions.isStepConverged())
00380       state = StepConverged;
00381     else if (stepReduced)
00382       directions.appendNewDirections();
00383   }
00384 
00385   // Check for number of evaluations
00386   if ((state == Continue) && 
00387       (isMaxEvaluations) && 
00388       (conveyor.getCounter().getNumEvaluations() >= maxEvaluations))
00389     state = EvaluationsExhausted; 
00390 }
00391 
00392 void APPSPACK::Solver::printInitializationInformation() const
00393 {
00394   cout << "\n";
00395   cout << "###########################################" << "\n";
00396   cout << "###   APPSPACK Initialization Results   ###" << endl;
00397   cout << "\n";
00398   
00399   // Paramters
00400   cout << "*** Parameter List ***" << "\n";
00401   params.print();
00402 
00403   // Constraints
00404   cout << "\n" << "*** Constraints ***" << "\n";
00405   constraints.print();
00406   
00407   // Conveyor
00408   cout << "\n" << "*** Conveyor ***" << "\n";
00409   conveyor.print();
00410 
00411   cout << "\n";
00412   cout << "### End APPSPACK Initialization Results ###" << endl;
00413   cout << "###########################################" << "\n";
00414 }
00415 
00416 void APPSPACK::Solver::printBestPoint(const string label) const
00417 {
00418   cout << "\n" << label << ": " << *bestPointPtr << endl;
00419 }
00420 
00421 ostream& operator<<(ostream& stream, APPSPACK::Solver::State state) 
00422 {
00423   switch (state)
00424   {
00425   case APPSPACK::Solver::Continue:
00426     stream << "Continue";
00427     break;
00428   case APPSPACK::Solver::StepConverged:
00429     stream << "Step Converged";
00430     break;
00431   case APPSPACK::Solver::FunctionConverged:
00432     stream << "Function Converged";
00433     break;
00434   case APPSPACK::Solver::EvaluationsExhausted:
00435     stream << "Evaluations Exhausted";
00436     break;
00437   }
00438 
00439   return stream;
00440 }
00441 
00442 void APPSPACK::Solver::writeSolutionFile(const string filename) const
00443 {
00444   // Determine output precision.
00445   int precision = params.getParameter("Solution File Precision", 14);
00446     
00447   // Open file for append
00448   ofstream fout;
00449   fout.open(filename.c_str(), std::ios::out);
00450   if (!fout)
00451   {
00452     cerr << "APPSPACK::Solver::writeSolutionFile() - Unable to open solution file" << endl;
00453     return;
00454   }
00455 
00456   // Write best value to outname using cache file output style.
00457   fout << "f=[ ";
00458   (bestPointPtr->getVecF()).leftshift(fout,precision);
00459   fout << " ]";
00460   
00461   fout << " x=[ ";
00462   (bestPointPtr->getX()).leftshift(fout,precision);
00463   fout << " ]" << endl;
00464 }

 

© 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