Main Page | Namespace List | Class Hierarchy | Alphabetical List | 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.30 2004/11/23 22:36:18 tgkolda 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
    00037 #include "APPSPACK_Solver.hpp"
    00038 #include "APPSPACK_Print.hpp"
    00039
    00040
    00041 APPSPACK::Solver::Solver(const Parameter::List& params_in,
    00042                          Executor::Interface& executor_in,
    00043                          const Constraints::Interface& constraints_in) :
    00044   constraints(constraints_in),
    00045   params(params_in),
    00046   print(params),                // modifies params and sets globals
    00047   bestPointPtr(initializeBestPointPtr(params, constraints)), // modifies params
    00048   directions(params, constraints), // modifies params
    00049   conveyor(params, constraints.getScaling(), executor_in), // modifies params
    00050   exchangeList(),
    00051   state(Continue),
    00052   isFunctionTolerance(false),
    00053   isMaxEvaluations(false)
    00054 {
    00055
    00056   cout << "\n"
    00057        << "-----------------------------------------------------\n"
    00058        << "APPSPACK: Asynchronous Parallel Pattern Search\n"
    00059        << "Written by T. G. Kolda et al., Sandia National Labs\n"
    00060        << "For more information visit \n"
    00061        << "http://software.sandia.gov/appspack\n"
    00062        << "-----------------------------------------------------\n"
    00063        << "\n";
    00064
    00065   // Stopping Criteria
    00066   if ((params.isParameterDouble("Function Tolerance")) || (params.isParameterValue("Function Tolerance")))
    00067   {
    00068     isFunctionTolerance = true;
    00069     if (params.isParameterDouble("Function Tolerance"))
    00070     {
    00071       functionTolerance = params.getDoubleParameter("Function Tolerance");
    00072     }
    00073     else if (params.isParameterValue("Function Tolerance"))
    00074     {
    00075       functionTolerance = params.getValueParameter("Function Tolerance");
    00076     }
    00077   }
    00078
    00079   if (params.isParameter("Maximum Evaluations"))
    00080   {
    00081     isMaxEvaluations = true;
    00082     maxEvaluations = params.getParameter("Maximum Evaluations", 0);
    00083   }
    00084
    00085   boundsTolerance = params.getParameter("Bounds Tolerance", params.getDoubleParameter("Step Tolerance") / 2.0);
    00086   if (boundsTolerance <= 0)
    00087   {
    00088     cout << "APPSPACK::Solver::Solver - Invalid non-positive value for \"Bounds Tolerance\"." << endl;
    00089     throw "APPSPACK Error";
    00090   }
    00091
    00092   // Machine Epsilon
    00093   epsMach = params.getParameter("Machine Epsilon", 1.0e-14);
    00094
    00096   if (Print::doPrint(Print::InitialData))
    00097     printInitializationInformation();
    00098
    00099   processNewBestPoint();
    00100 }
    00101
    00102 APPSPACK::Solver::~Solver()
    00103 {
    00104   delete bestPointPtr;
    00105 }
    00106
    00107 const APPSPACK::Vector& APPSPACK::Solver::getBestX() const
    00108 {
    00109   return bestPointPtr->getX();
    00110 }
    00111
    00112 bool APPSPACK::Solver::isBestF() const
    00113 {
    00114   return bestPointPtr->getF().getIsValue();
    00115 }
    00116
    00117 double APPSPACK::Solver::getBestF() const
    00118 {
    00119   return bestPointPtr->getF().getValue();
    00120 }
    00121
    00122
    00123 APPSPACK::Point* APPSPACK::Solver::initializeBestPointPtr(Parameter::List& params, const Constraints::Interface& constraints)
    00124 {
    00125   // Note that this needs to be called before the Cache constructor because of the Step Tolerance!
    00126
    00127   const Vector& lower = constraints.getLower();
    00128   const Vector& upper = constraints.getUpper();
    00129   const vector<bool>& isLower = constraints.getIsLower();
    00130   const vector<bool>& isUpper = constraints.getIsUpper();
    00131
    00132   // Initial point
    00133   Vector nominalX;
    00134   if (!params.isParameter("InitialX"))
    00135   {
    00136     nominalX.resize(constraints.getScaling().size());
    00137     for (int i = 0; i < nominalX.size(); i ++)
    00138     {
    00139       if ((isUpper[i]) && (isLower[i]))
    00140         nominalX[i] = ( upper[i] + lower[i] ) / 2.0 ;
    00141       else if (isUpper[i])
    00142         nominalX[i] = upper[i];
    00143       else if (isLower[i])
    00144         nominalX[i] = lower[i];
    00145       else
    00146         nominalX[i] = 0;
    00147     }
    00148   }
    00149   Vector initialX = params.getParameter("Initial X", nominalX);
    00150
    00151   // Check that the initial point is the right size
    00152   if (initialX.size() != lower.size())
    00153   {
    00154     cerr << "Error: Size mismatch\n";
    00155     cerr << "The size of the initial X is not the same as the lower bounds.\n";
    00156     throw "APPSPACK Error";
    00157   }
    00158
    00159   // Check that initial point is feasible
    00160   for (int i = 0; i < initialX.size(); i ++)
    00161   {
    00162     if (isUpper[i])
    00163       if (initialX[i] > upper[i])
    00164       {
    00165         cerr << "ERROR: Infeasible initial point\n";
    00166         cerr << "The " << i+1 << " component of the initial X (" << initialX[i]
    00167              << ") is greater than the upper bound (" << upper[i] << ")" << endl;
    00168         throw "APPSPACK Error";
    00169       }
    00170     if (isLower[i])
    00171       if (initialX[i] < lower[i])
    00172       {
    00173         cerr << "ERROR: Infeasible initial point\n";
    00174         cerr << "The " << i+1 << " component of the initial X (" << initialX[i]
    00175              << ") is lower than the lower bound (" << lower[i] << ")" << endl;
    00176         throw "APPSPACK Error";
    00177       }
    00178   }
    00179
    00180   // Initial F
    00181   Value initialF;
    00182   if (params.isParameterDouble("Initial F"))
    00183   {
    00184     initialF = params.getDoubleParameter("Initial F");
    00185   }
    00186   else if (params.isParameterValue("Initial F"))
    00187   {
    00188     initialF = params.getValueParameter("Initial F");
    00189   }
    00190
    00191   // Initial step
    00192   double initialStep = params.getParameter("Initial Step", 1.0);
    00193
    00194   // Point
    00195   double alpha = params.getParameter("Sufficient Decrease Factor", 0.01);
    00196
    00197   // Create first best point
    00198   return new Point(initialX, initialF, initialStep, alpha);
    00199 }
    00200
    00201 APPSPACK::Solver::State APPSPACK::Solver::solve()
    00202 {
    00203   while (state == Continue)
    00204   {
    00205     iterate();
    00206   }
    00207
    00208   // Print the final solution
    00209   if (Print::doPrint(Print::FinalSolution))
    00210   {
    00211     cout << "\nFinal State: " << state << "\n";
    00212     printBestPoint("Final Min");
    00213     directions.print("Final Directions");
    00214     conveyor.getCounter().print();
    00215   }
    00216
    00217   return state;
    00218 }
    00219
    00220 APPSPACK::Solver::State APPSPACK::Solver::iterate()
    00221 {
    00222     generateTrialPoints();
    00223
    00224     conveyor.exchange(exchangeList);
    00225
    00226     processEvaluatedTrialPoints();
    00227
    00228     return state;
    00229 }
    00230
    00231 // PRIVATE
    00232 void APPSPACK::Solver::processNewBestPoint(Point* newBestPointPtr)
    00233 {
    00234   // Update the best point
    00235   if (newBestPointPtr != NULL)
    00236   {
    00237     delete bestPointPtr;
    00238     bestPointPtr = newBestPointPtr;
    00239   }
    00240
    00241   // Print
    00242   if (Print::doPrint(Print::NewBestPoint))
    00243     printBestPoint("New Min");
    00244
    00245   // Check for convergence based on the function value
    00246   if ((isFunctionTolerance) && (bestPointPtr->getF() < functionTolerance))
    00247   {
    00248     state = FunctionConverged;
    00249     return;
    00250   }
    00251
    00252   // Update the (scaled) search directions
    00253   directions.computeNewDirections(*bestPointPtr);
    00254
    00255   // Print
    00256   if (Print::doPrint(Print::NewBestDirections))
    00257     directions.print("New Directions");
    00258 }
    00259
    00260 // PRIVATE
    00261 void APPSPACK::Solver::generateTrialPoints()
    00262 {
    00263   // Local references
    00264   const Vector& parentX = bestPointPtr->getX();
    00265   //const Vector& scaling = constraints.getScaling();
    00266   const Vector& lower = constraints.getLower();
    00267   const Vector& upper = constraints.getUpper();
    00268   const vector<bool>& isLower = constraints.getIsLower();
    00269   const vector<bool>& isUpper = constraints.getIsUpper();
    00270
    00271   Vector& x = tmpVector;
    00272   int n = parentX.size();
    00273   x.resize(n);
    00274
    00275   const vector<int> indices = directions.getDirectionIndices();
    00276
    00277   for (int i = 0; i < indices.size(); i ++)
    00278   {
    00279     int idx = indices[i];
    00280     const Vector& direction = directions.getDirection(idx);
    00281     const double step = directions.getStep(idx);
    00282
    00283     // Used in the calculation of the longest feasible step
    00284     double tmpStep = step;
    00285
    00286     // Calculate the new point
    00287     for (int k = 0; k < n; k ++)
    00288       x[k] = parentX[k] + tmpStep * direction[k];
    00289
    00290     // Correct violations of the constraints
    00291     for (int j = 0; j < n; j ++)
    00292     {
    00293       if ((isLower[j]) && ((lower[j] - x[j]) > 0))
    00294       {
    00295         if (lower[j] - x[j] > epsMach)
    00296         {
    00297           tmpStep = ( lower[j] - parentX[j] ) / direction[j];
    00298           for (int k = 0; k < n; k ++)
    00299             x[k] = parentX[k] + tmpStep * direction[k];
    00300         }
    00301         else
    00302           x[j] = lower[j];      // nudge onto exact boundary
    00303       }
    00304       else if ((isUpper[j]) && ((x[j] - upper[j]) > 0))
    00305       {
    00306         if ((x[j] - upper[j]) > 0)
    00307         {
    00308           tmpStep = ( upper[j] - parentX[j] ) / direction[j];
    00309           for (int k = 0; k < n; k ++)
    00310             x[k] = parentX[k] + tmpStep * direction[k];
    00311         }
    00312         else
    00313           x[j] = upper[j];      // nudge onto exact boundary
    00314       }
    00315     }
    00316
    00317     /*
    00318
    00319     // Snap nearby points to the constraints
    00320     for (int j = 0; j < n; j ++)
    00321     {
    00322       if ( (isLower[j]) &&
    00323            (abs(lower[j] - x[j]) < (boundsTolerance * scaling[j]) ) )
    00324       {
    00325         x[j] = lower[j];
    00326       }
    00327       else if ( (isUpper[j]) &&
    00328                 (abs(upper[j] - x[j]) < (boundsTolerance * scaling[j]) ) )
    00329       {
    00330         x[j] = upper[j];
    00331       }
    00332     }
    00333
    00334     */
    00335
    00336     // Create a new trial point
    00337     Point* newPointPtr = new Point(x, step, *bestPointPtr, idx);
    00338
    00339     // Save off trial point information
    00340     directions.setTrueStepAndTag(idx, tmpStep, newPointPtr->getTag());
    00341
    00342     // Push this trial point onto the new trial point list.
    00343     // Ownership of the pointer transfers to the list.
    00344     exchangeList.push(newPointPtr);
    00345
    00346   }
    00347
    00348   if (Print::doPrint(Print::Directions))
    00349     directions.print("Directions After Trial Point Generation");
    00350
    00351   if (Print::doPrint(Print::UnevaluatedPoints))
    00352     exchangeList.print("Trial Points Before Conveyor");
    00353
    00354 }
    00355
    00356 void APPSPACK::Solver::processEvaluatedTrialPoints()
    00357 {
    00358   // Print
    00359   if (Print::doPrint(Print::EvaluatedPoints))
    00360     exchangeList.print("Trial Points After Conveyor");
    00361
    00362   // Check for a new best point
    00363   if (exchangeList.best() < *bestPointPtr)
    00364   {
    00365     processNewBestPoint(exchangeList.popBest());
    00366     exchangeList.prune();
    00367     conveyor.prune();
    00368   }
    00369   else
    00370   {
    00371     // Otherwise, just process the list
    00372     Point* ptr;
    00373     while (!exchangeList.isEmpty())
    00374     {
    00375       ptr = exchangeList.pop();
    00376       if (ptr->getParentTag() == bestPointPtr->getTag())
    00377         directions.reduceStep(ptr->getIndex());
    00378       delete ptr;
    00379     }
    00380
    00381     // Check for step length convergence
    00382     if (directions.isStepConverged())
    00383       state = StepConverged;
    00384   }
    00385
    00386   // Check for number of evaluations
    00387   if ((state == Continue) &&
    00388       (isMaxEvaluations) &&
    00389       (conveyor.getCounter().getNumEvaluations() >= maxEvaluations))
    00390     state = EvaluationsExhausted;
    00391 }
    00392
    00393 void APPSPACK::Solver::printInitializationInformation() const
    00394 {
    00395   cout << "\n";
    00396   cout << "###########################################" << "\n";
    00397   cout << "###   APPSPACK Initialization Results   ###" << endl;
    00398   cout << "\n";
    00399
    00400   // Paramters
    00401   cout << "*** Parameter List ***" << "\n";
    00402   params.print();
    00403
    00404   // Constraints
    00405   cout << "\n" << "*** Constraints ***" << "\n";
    00406   constraints.print();
    00407
    00408   // Conveyor
    00409   cout << "\n" << "*** Conveyor ***" << "\n";
    00410   conveyor.print();
    00411
    00412   cout << "\n";
    00413   cout << "### End APPSPACK Initialization Results ###" << endl;
    00414   cout << "###########################################" << "\n";
    00415 }
    00416
    00417 void APPSPACK::Solver::printBestPoint(const string label) const
    00418 {
    00419   cout << "\n" << label << ": " << *bestPointPtr << endl;
    00420 }
    00421
    00422 ostream& operator<<(ostream& stream, APPSPACK::Solver::State state)
    00423 {
    00424   switch (state)
    00425   {
    00426   case APPSPACK::Solver::Continue:
    00427     stream << "Continue";
    00428     break;
    00429   case APPSPACK::Solver::StepConverged:
    00430     stream << "Step Converged";
    00431     break;
    00432   case APPSPACK::Solver::FunctionConverged:
    00433     stream << "Function Converged";
    00434     break;
    00435   case APPSPACK::Solver::EvaluationsExhausted:
    00436     stream << "Evaluations Exhausted";
    00437     break;
    00438   }
    00439
    00440   return stream;
    00441 }
    00442

© Sandia Corporation | Site Contact | Privacy and Security

Generated on Wed Dec 14 18:41:05 2005 for APPSPACK 4.0.2 by 1.3.8 written by Dimitri van Heesch, © 1997-2002