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