1 #ifndef SimTK_SIMMATH_OPTIMIZER_H_ 2 #define SimTK_SIMMATH_OPTIMIZER_H_ 3 4 /* -------------------------------------------------------------------------- * 5 * Simbody(tm): SimTKmath * 6 * -------------------------------------------------------------------------- * 7 * This is part of the SimTK biosimulation toolkit originating from * 8 * Simbios, the NIH National Center for Physics-Based Simulation of * 9 * Biological Structures at Stanford, funded under the NIH Roadmap for * 10 * Medical Research, grant U54 GM072970. See https://simtk.org/home/simbody. * 11 * * 12 * Portions copyright (c) 2006-13 Stanford University and the Authors. * 13 * Authors: Jack Middleton * 14 * Contributors: Michael Sherman * 15 * * 16 * Licensed under the Apache License, Version 2.0 (the "License"); you may * 17 * not use this file except in compliance with the License. You may obtain a * 18 * copy of the License at http://www.apache.org/licenses/LICENSE-2.0. * 19 * * 20 * Unless required by applicable law or agreed to in writing, software * 21 * distributed under the License is distributed on an "AS IS" BASIS, * 22 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. * 23 * See the License for the specific language governing permissions and * 24 * limitations under the License. * 25 * -------------------------------------------------------------------------- */ 26 27 28 #include "SimTKcommon.h" 29 #include "simmath/internal/common.h" 30 #include "simmath/Differentiator.h" 31 32 namespace SimTK { 33 34 /** 35 * The available Optimizer algorithms. 36 * Gradient descent algorithms seek to find a local minimum, and are not 37 * guaranteed to find the global minimum. See the description of Optimizer for 38 * specific information about how to use the algorithms. 39 */ 40 enum OptimizerAlgorithm { 41 /// Simmath will select best Optimizer based on problem type. 42 BestAvailable = 0, 43 /// IpOpt algorithm (https://projects.coin-or.org/ipopt); 44 /// gradient descent. 45 InteriorPoint = 1, 46 /// Limited-memory Broyden-Fletcher-Goldfarb-Shanno algorithm; 47 /// gradient descent. 48 LBFGS = 2, 49 /// LBFGS with simple bound constraints; 50 /// gradient descent. 51 LBFGSB = 3, 52 /// C implementation of sequential quadratic programming 53 /// (requires external library: 54 /// ftp://frcatel.fri.uniza.sk/pub/soft/math/matprog/doc/fsqp.html); 55 /// gradient descent. 56 CFSQP = 4, 57 /// Covariance matrix adaptation, evolution strategy 58 /// (https://github.com/cma-es/c-cmaes); 59 /// this is a randomized algorithm that attempts to find a global minimum. 60 CMAES = 5, 61 UnknownOptimizerAlgorithm = 6, // the default impl. of getAlgorithm. 62 /// An algorithm that is implemented outside of Simmath. 63 UserSuppliedOptimizerAlgorithm = 7 64 }; 65 66 /** 67 * Abstract class which defines an objective/cost function which is optimized by 68 * and Optimizer object. The OptimizerSystem also defines any constraints which 69 * must be satisfied. 70 */ 71 class SimTK_SIMMATH_EXPORT OptimizerSystem { 72 public: OptimizerSystem()73 OptimizerSystem() : numParameters(0), 74 numEqualityConstraints(0), 75 numInequalityConstraints(0), 76 numLinearEqualityConstraints(0), 77 numLinearInequalityConstraints(0), 78 useLimits( false ), 79 lowerLimits(0), 80 upperLimits(0) { 81 } 82 OptimizerSystem(int nParameters)83 explicit OptimizerSystem(int nParameters ) { 84 new (this) OptimizerSystem(); // call the above constructor 85 setNumParameters(nParameters); 86 } 87 ~OptimizerSystem()88 virtual ~OptimizerSystem() { 89 if( useLimits ) { 90 delete lowerLimits; 91 delete upperLimits; 92 } 93 } 94 95 /// Objective/cost function which is to be optimized; return 0 when successful. 96 /// The value of f upon entry into the function is undefined. 97 /// This method must be supplied by concrete class. objectiveFunc(const Vector & parameters,bool new_parameters,Real & f)98 virtual int objectiveFunc ( const Vector& parameters, 99 bool new_parameters, Real& f ) const { 100 SimTK_THROW2(SimTK::Exception::UnimplementedVirtualMethod , "OptimizerSystem", "objectiveFunc" ); 101 return -1; } 102 103 /// Computes the gradient of the objective function; return 0 when successful. 104 /// This method does not have to be supplied if a numerical gradient is used. gradientFunc(const Vector & parameters,bool new_parameters,Vector & gradient)105 virtual int gradientFunc ( const Vector ¶meters, 106 bool new_parameters, Vector &gradient ) const { 107 SimTK_THROW2(SimTK::Exception::UnimplementedVirtualMethod , "OptimizerSystem", "gradientFunc" ); 108 return -1; } 109 /// Computes the value of the constraints; return 0 when successful. 110 /// This method must be supplied if the objective function has constraints. constraintFunc(const Vector & parameters,bool new_parameters,Vector & constraints)111 virtual int constraintFunc ( const Vector & parameters, 112 bool new_parameters, Vector & constraints ) const { 113 SimTK_THROW2(SimTK::Exception::UnimplementedVirtualMethod , "OptimizerSystem", "constraintFunc" ); 114 return -1; } 115 /// Computes Jacobian of the constraints; return 0 when successful. 116 /// This method does not have to be supplied if a numerical jacobian is used. constraintJacobian(const Vector & parameters,bool new_parameters,Matrix & jac)117 virtual int constraintJacobian ( const Vector& parameters, 118 bool new_parameters, Matrix& jac ) const { 119 SimTK_THROW2(SimTK::Exception::UnimplementedVirtualMethod , "OptimizerSystem", "constraintJacobian" ); 120 return -1; } 121 /// Computes Hessian of the objective function; return 0 when successful. 122 /// This method does not have to be supplied if limited memory is used. hessian(const Vector & parameters,bool new_parameters,Vector & gradient)123 virtual int hessian ( const Vector ¶meters, 124 bool new_parameters, Vector &gradient) const { 125 SimTK_THROW2(SimTK::Exception::UnimplementedVirtualMethod , "OptimizerSystem", "hessian" ); 126 return -1; } 127 128 /// Sets the number of parameters in the objective function. setNumParameters(const int nParameters)129 void setNumParameters( const int nParameters ) { 130 if( nParameters < 1 ) { 131 const char* where = " OptimizerSystem Constructor"; 132 const char* szName = "number of parameters"; 133 SimTK_THROW5(SimTK::Exception::ValueOutOfRange, szName, 1, nParameters, INT_MAX, where); 134 } else { 135 numParameters = nParameters; 136 } 137 } 138 /// Sets the number of equality constraints. setNumEqualityConstraints(const int n)139 void setNumEqualityConstraints( const int n ) { 140 if( n < 0 ) { 141 const char* where = " OptimizerSystem setNumEqualityConstraints"; 142 const char* szName = "number of equality constraints"; 143 SimTK_THROW3(SimTK::Exception::SizeWasNegative, szName, n, where); 144 } else { 145 numEqualityConstraints = n; 146 } 147 } 148 /// Sets the number of inequality constraints. setNumInequalityConstraints(const int n)149 void setNumInequalityConstraints( const int n ) { 150 if( n < 0 ) { 151 const char* where = " OptimizerSystem setNumInequalityConstraints"; 152 const char* szName = "number of inequality constraints"; 153 SimTK_THROW3(SimTK::Exception::SizeWasNegative, szName, n, where); 154 } else { 155 numInequalityConstraints = n; 156 } 157 } 158 /// Sets the number of lineaer equality constraints. setNumLinearEqualityConstraints(const int n)159 void setNumLinearEqualityConstraints( const int n ) { 160 if( n < 0 || n > numEqualityConstraints ) { 161 const char* where = " OptimizerSystem setNumLinearEqualityConstraints"; 162 const char* szName = "number of linear equality constraints"; 163 SimTK_THROW4(SimTK::Exception::SizeOutOfRange, szName, n, numEqualityConstraints, where); 164 } else { 165 numLinearEqualityConstraints = n; 166 } 167 } 168 /// Sets the number of lineaer inequality constraints. setNumLinearInequalityConstraints(const int n)169 void setNumLinearInequalityConstraints( const int n ) { 170 if( n < 0 || n > numInequalityConstraints ) { 171 const char* where = " OptimizerSystem setNumLinearInequalityConstraints"; 172 const char* szName = "number of linear inequality constraints"; 173 SimTK_THROW4(SimTK::Exception::SizeOutOfRange, szName, n, numInequalityConstraints, where); 174 } else { 175 numLinearInequalityConstraints = n; 176 } 177 } 178 /// Set the upper and lower bounds on the parameters. setParameterLimits(const Vector & lower,const Vector & upper)179 void setParameterLimits( const Vector& lower, const Vector& upper ) { 180 if( upper.size() != numParameters && upper.size() != 0) { 181 const char* where = " OptimizerSystem setParametersLimits"; 182 const char* szName = "upper limits length"; 183 SimTK_THROW5(Exception::IncorrectArrayLength, szName, upper.size(), "numParameters", numParameters, where); 184 } 185 if( lower.size() != numParameters && lower.size() != 0 ) { 186 const char* where = " OptimizerSystem setParametersLimits"; 187 const char* szName = "lower limits length"; 188 SimTK_THROW5(Exception::IncorrectArrayLength, szName, lower.size(), "numParameters", numParameters, where); 189 } 190 191 // set the upper and lower limits 192 if( useLimits ) { 193 delete lowerLimits; 194 delete upperLimits; 195 } 196 197 if( upper.size() == 0 ) { 198 useLimits = false; 199 } else { 200 lowerLimits = new Vector( lower ); 201 upperLimits = new Vector( upper ); 202 useLimits = true; 203 } 204 } 205 206 /// Returns the number of parameters, that is, the number of variables that 207 /// the Optimizer may adjust while searching for a solution. getNumParameters()208 int getNumParameters() const {return numParameters;} 209 /// Returns the total number of constraints. getNumConstraints()210 int getNumConstraints() const {return numEqualityConstraints+numInequalityConstraints;} 211 /// Returns the number of equality constraints. getNumEqualityConstraints()212 int getNumEqualityConstraints() const {return numEqualityConstraints;} 213 /// Returns the number of inequality constraints. getNumInequalityConstraints()214 int getNumInequalityConstraints() const {return numInequalityConstraints;} 215 /// Returns the number of linear equality constraints. getNumLinearEqualityConstraints()216 int getNumLinearEqualityConstraints() const {return numLinearEqualityConstraints;} 217 /// Returns the number of nonlinear equality constraints. getNumNonlinearEqualityConstraints()218 int getNumNonlinearEqualityConstraints() const {return numEqualityConstraints-numLinearEqualityConstraints;} 219 /// Returns the number of linear inequality constraints. getNumLinearInequalityConstraints()220 int getNumLinearInequalityConstraints() const {return numLinearInequalityConstraints;} 221 /// Returns the number of linear inequality constraints. getNumNonlinearInequalityConstraints()222 int getNumNonlinearInequalityConstraints() const {return numInequalityConstraints-numLinearInequalityConstraints;} 223 224 /// Returns true if there are limits on the parameters. getHasLimits()225 bool getHasLimits() const { return useLimits; } 226 /// Returns the limits on the allowed values of each parameter, as 227 /// an array of lower bounds and an array of upper bounds, with 228 /// assumed lengths matching the number of parameters. getParameterLimits(Real ** lower,Real ** upper)229 void getParameterLimits( Real **lower, Real **upper ) const { 230 *lower = &(*lowerLimits)[0]; 231 *upper = &(*upperLimits)[0]; 232 } 233 234 private: 235 int numParameters; 236 int numEqualityConstraints; 237 int numInequalityConstraints; 238 int numLinearEqualityConstraints; 239 int numLinearInequalityConstraints; 240 bool useLimits; 241 Vector* lowerLimits; 242 Vector* upperLimits; 243 244 }; // class OptimizerSystem 245 246 /** 247 * API for SimTK Simmath's optimizers. 248 * An optimizer finds a minimum to an objective function. Usually, this minimum 249 * is a local minimum. Some algorithms, like CMAES, are designed to find the 250 * global minumum. The optimizer can be constrained to search for a minimum 251 * within a feasible region. The feasible region is defined in two ways: via 252 * limits on the parameters of the objective function; and, for algorithms 253 * other than CMAES, by supplying constraint functions that must be satisfied. 254 * The optimizer starts searching for a minimum beginning at a user supplied 255 * initial value for the set of parameters. 256 * 257 * The objective function and constraints are specified by supplying the 258 * Optimizer with a concrete implemenation of an OptimizerSystem class. 259 * The OptimizerSystem can be passed to the Optimizer either through the 260 * Optimizer constructor or by calling the Optimizer::setOptimizerSystem 261 * method. The Optimizer class will select the best optimization algorithm to 262 * solve the problem based on the constraints supplied by the OptimizerSystem. 263 * A user can also override the optimization algorithm selected by the 264 * Optimizer by specifying the optimization algorithm. 265 * 266 * <h3> Optimization algorithms and advanced options </h3> 267 * 268 * See OptimizerAlgorithm for a brief description of the available algorithms. 269 * Most of these algorithms have options that are specific to the algorithm. 270 * These options are set via methods like Optimizer::setAdvancedStrOption. If 271 * you want to get going quickly, you can just use the default values of these 272 * options and ignore this section. As an example, an int option 273 * <b>popsize</b> would be set via: 274 * 275 * @code 276 * opt.setAdvancedIntOption("popsize", 5); 277 * @endcode 278 * 279 * For now, we only have detailed documentation for the CMAES algorithm. 280 * 281 * <h4> CMAES </h4> 282 * 283 * This is the c-cmaes algorithm written by Niko Hansen 284 * (https://github.com/cma-es/c-cmaes). 285 * 286 * Some notes: 287 * - This algorithm obeys parameter limits. 288 * - This is a derivative-free optimization algorithm, so methods like the 289 * following have no effect: 290 * - Optimizer::useNumericalGradient 291 * - Optimizer::setDifferentiatorMethod 292 * - Optimizer::setLimitedMemoryHistory 293 * - OptimizerSystem::gradientFunc 294 * - OptimizerSystem::hessian 295 * - This algorithm does not obey constraint functions, so methods like the 296 * following have no effect: 297 * - Optimizer::setConstraintTolerance 298 * - Optimizer::useNumericalJacobian 299 * - OptimizerSystem::constraintFunc 300 * - OptimizerSystem::constraintJacobian 301 * - OptimizerSystem::setNumEqualityConstraints 302 * - OptimizerSystem::setNumInequalityConstraints 303 * - OptimizerSystem::setNumLinearEqualityConstraints 304 * - OptimizerSystem::setNumLinearInequalityConstraints 305 * - The effect of the diagnostics level is as follows: 306 * - 0: minimal output to console (warnings, errors), some files are 307 * written to the current directory (errcmaes.err error log). 308 * - 1: additional output to console. 309 * - 2: all files are written to the current directory. 310 * - 3: output to console, and all files are written to the current directory. 311 * 312 * 313 * Encoding of Variables 314 * 315 * Inappropriate initialization of the algorithm may lead to resampling of the 316 * parameter distribution. Since some parameters may be defined in a bounded 317 * region [a, b]. The CMA algorithm involves sampling a random distribution for 318 * the variables (parameters) in the optimization problem. If the sampled values 319 * do not lie within the variables' bounds, CMA must resample the distribution 320 * until the variables lie within the bounds. This resampling is undesirable, 321 * and may prevent the algorithm from functioning properly. There are two ways 322 * to avoid excessive resampling:To overcome this the problem the user can 323 * formulate the problem as follows: 324 * 325 * 1) Either define the bounds of the parameter and choose the appropriate stddev 326 * and initial value for each parameter (see init_stepsize) 327 * 2) Or to reformulate the problem, by rescaling the parameter space so that 328 * each parameter map in a region between e.g. [0, 1] and each parameter has an 329 * initial value of 0.5 and stddev of 0.2 (see comments below). 330 * 331 * The specific formulation of a (real) optimization problem has a tremendous 332 * impact on the optimization performance. In particular, a reasonable parameter 333 * encoding is essential. All parameters should be rescaled such that they have 334 * presumably similar sensitivity (this makes the identity as initial covariance 335 * matrix the right choice). Usually, the best approach is to write a wrapper 336 * around the objective function that transforms the parameters before the actual 337 * function call. The wrapper scales, for example, in each parameter/coordinate 338 * the value [0; 10] into the typical actual domain of the parameter/coordinate. 339 * 340 * The natural encoding of (some of) the parameters can also be "logarithmic". 341 * That is, for a parameter that must always be positive, with a ratio between 342 * typical upper and lower value being larger than 100, we might use 10x instead 343 * of x to call the objective function. More specifically, to achieve the parameter 344 * range [10^–4,10^–1], we use 10^–4×10^3x/10 with x in [0; 10]. Again, the idea 345 * is to have similar sensitivity: this makes sense if we expect the change from 346 * 10^–4 to 10^–3 to have an impact similar to the change from 10^–2 to 10^–1. 347 * In order to avoid the problem that changes of very small values have too less 348 * an impact, an alternative is to choose 10^–1 × (x/10)2 >= 0. In the case where 349 * only a lower bound at zero is necessary, a simple and natural transformation is 350 * x2 × default_x, such that x=1 represents the default (or initial) value and x 351 * remains unbounded during optimization. 352 * 353 * In summary, to map the values [0;10] into [a;b] we have the alternative 354 * transformations a + (b-a) × x/10 or a + (b-a) × (x/10)2 >= a or a × (b/a)x/10 >= 0. 355 * 356 * For more details see: https://www.lri.fr/~hansen/cmaes_inmatlab.html 357 * 358 * Advanced options: 359 * 360 * The default values for options whose name begins with "stop" are specified 361 * at https://github.com/CMA-ES/c-cmaes/blob/master/cmaes_initials.par 362 * 363 * - <b>popsize</b> (int; default: depends on number of parameters) The 364 * population size (also known as lambda). 365 * - <b>init_stepsize</b> (Vector/Real; default: 0.3) Initial stddev; After the 366 * encoding of variables, the initial solution point x0 and the initial standard 367 * deviation (step_size) sigma0 must be chosen. In a practical application, one 368 * often wants to start by trying to improve a given solution locally. In this 369 * case we choose a rather small sigma0 (say in [0.001, 0.1], given the x-values 370 * "live" in [0,10]). Thereby we can also check whether the initial solution is 371 * possibly a local optimum. When a global optimum is sought-after on rugged or 372 * multimodal landscapes, sigma0 should be chosen such that the final desirable 373 * location (or at least some of its domain of attraction) is not far outside of 374 * x0 ± 2sigma0 in each coordinate. (Remark that in Rn, if each boundary domain is 375 * in distance sigma, then the boundary corner is sigma*sqrt(n) away, which poses 376 * a slight dilemma for larger n.) 377 * 378 * A warning is emitted if this is not set and default value is used for each 379 * variable. 380 * 381 * Example setting the init_stepsize: 382 * 383 * Vector initStepSize(N, 0.3); 384 * opt.setAdvancedVectorOption("init_stepsize", initStepSize); 385 * "or" 386 * opt.setAdvancedRealOption("init_stepsize", 0.3); 387 * - <b>seed</b> (int; default: 0, which uses clock time) Seed for the random 388 * number generator that is used to sample the population from a normal 389 * distribution. See note below. 390 * - <b>maxTimeFractionForEigendecomposition</b> (real; default: 0.2) 391 * Controls the amount of time spent generating eigensystem 392 * decompositions. 393 * - <b>stopMaxFunEvals</b> (int) Stop optimization after this 394 * number of evaluations of the objective function. 395 * - <b>stopFitness</b> (real) Stop if function value is smaller than 396 * stopFitness. 397 * - <b>stopTolFunHist</b> (real) Stop if function value differences of best 398 * values are smaller than stopTolFunHist. 399 * - <b>stopTolX</b> (real) Stop if step sizes are smaller than stopTolX. 400 * - <b>stopTolUpXFactor</b> (real) Stop if standard deviation increases 401 * by more than stopTolUpXFactor. 402 * - <b>parallel</b> (str) To run the optimization with multiple threads, set 403 * this to "multithreading". Only use this if your OptimizerSystem is 404 * threadsafe: you can't reliably modify any mutable variables in your 405 * OptimizerSystem::objectiveFun(). 406 * - <b>nthreads</b> (int) If the <b>parallel</b> option is set to 407 * "multithreading", this is the number of threads to use (by default, this 408 * is the number of processors/threads on the machine). 409 * 410 * If you want to generate identical results with repeated optimizations, 411 * you can set the <b>seed</b> option. In addition, you *must* set the 412 * <b>maxTimeFractionForEigendecomposition</b> option to be greater than or 413 * equal to 1.0. 414 * 415 * @code 416 * opt.setAdvancedIntOption("seed", 42); 417 * opt.setAdvancedRealOption("maxTimeFractionForEigendecomposition", 1); 418 * @endcode 419 * 420 */ 421 class SimTK_SIMMATH_EXPORT Optimizer { 422 public: 423 Optimizer(); 424 Optimizer( const OptimizerSystem& sys); 425 Optimizer( const OptimizerSystem& sys, OptimizerAlgorithm algorithm); 426 ~Optimizer(); 427 428 /// BestAvailable, UnknownAlgorithm, and UserSuppliedAlgorithm 429 /// are treated as never available. 430 static bool isAlgorithmAvailable(OptimizerAlgorithm algorithm); 431 432 /// Sets the relative accuracy used determine if the problem has converged. 433 void setConvergenceTolerance(Real accuracy ); 434 /// Sets the absolute tolerance used to determine whether constraint 435 /// violation is acceptable. 436 void setConstraintTolerance(Real tolerance); 437 438 439 /// Set the maximum number of iterations allowed of the optimization 440 /// method's outer stepping loop. Most optimizers also have an inner loop 441 /// ("line search") which is also iterative but is not affected by this 442 /// setting. Inner loop convergence is typically prescribed by theory, and 443 /// failure there is often an indication of an ill-formed problem. 444 void setMaxIterations( int iter ); 445 /// Set the maximum number of previous hessians used in a limited memory 446 /// hessian approximation. 447 void setLimitedMemoryHistory( int history ); 448 /// Set the level of debugging info displayed. 449 void setDiagnosticsLevel( int level ); 450 451 void setOptimizerSystem( const OptimizerSystem& sys ); 452 void setOptimizerSystem( const OptimizerSystem& sys, OptimizerAlgorithm algorithm ); 453 454 /// Set the value of an advanced option specified by a string. 455 bool setAdvancedStrOption( const char *option, const char *value ); 456 /// Set the value of an advanced option specified by a real value. 457 bool setAdvancedRealOption( const char *option, const Real value ); 458 /// Set the value of an advanced option specified by an integer value. 459 bool setAdvancedIntOption( const char *option, const int value ); 460 /// Set the value of an advanced option specified by an boolean value. 461 bool setAdvancedBoolOption( const char *option, const bool value ); 462 /// Set the value of an advanced option specified by an Vector value. 463 bool setAdvancedVectorOption( const char *option, const Vector value ); 464 465 466 /// Set which numerical differentiation algorithm is to be used for the next 467 /// useNumericalGradient() or useNumericalJacobian() call. Choices are 468 /// Differentiator::ForwardDifference (first order) or 469 /// Differentiator::CentralDifference (second order) with central the 470 /// default. 471 /// @warning This has no effect if you have already called 472 /// useNumericalGradient() or useNumericalJacobian(). Those must be called 473 /// \e after setDifferentiatorMethod(). 474 /// @see SimTK::Differentiator 475 void setDifferentiatorMethod(Differentiator::Method method); 476 /// Return the differentiation method last supplied in a call to 477 /// setDifferentiatorMethod(), \e not necessarily the method currently 478 /// in use. See setDifferentiatorMethod() for more information. 479 /// @see SimTK::Differentiator 480 Differentiator::Method getDifferentiatorMethod() const; 481 482 /// Return the algorithm used for the optimization. You may be interested 483 /// in this value if you didn't specify an algorithm, or specified for 484 /// Simbody to choose the BestAvailable algorithm. This method won't return 485 /// BestAvailable, even if it's the 'algorithm' that you chose. 486 OptimizerAlgorithm getAlgorithm() const; 487 488 /// Enable numerical calculation of gradient, with optional estimation of 489 /// the accuracy to which the objective function is calculated. For example, 490 /// if you are calculate about 6 significant digits, supply the estimated 491 /// accuracy as 1e-6. Providing the estimated accuracy improves the quality 492 /// of the calculated derivative. If no accuracy is provided we'll assume 493 /// the objective is calculated to near machine precision. The method used 494 /// for calculating the derivative will be whatever was \e previously 495 /// supplied in a call to setDifferentiatorMethod(), or the default which 496 /// is to use central differencing (two function evaluations per 497 /// gradient entry). See SimTK::Differentiator for more information. 498 /// @see setDifferentiatorMethod(), SimTK::Differentiator 499 void useNumericalGradient(bool flag, 500 Real estimatedAccuracyOfObjective = SignificantReal); 501 /// Enable numerical calculation of the constraint Jacobian, with optional 502 /// estimation of the accuracy to which the constraint functions are 503 /// calculated. For example, if you are calculating about 6 significant 504 /// digits, supply the estimated accuracy as 1e-6. Providing the estimated 505 /// accuracy improves the quality of the calculated derivative. If no 506 /// accuracy is provided we'll assume the constraints are calculated to near 507 /// machine precision. The method used for calculating the derivative will 508 /// be whatever was \e previously supplied in a call to 509 /// setDifferentiatorMethod(), or the default which is to use central 510 /// differencing (two function evaluations per Jacobian column. See 511 /// SimTK::Differentiator for more information. 512 /// @see setDifferentiatorMethod(), SimTK::Differentiator 513 void useNumericalJacobian(bool flag, 514 Real estimatedAccuracyOfConstraints = SignificantReal); 515 516 /// Compute optimization. 517 Real optimize(Vector&); 518 519 /// Return a reference to the OptimizerSystem currently associated with this Optimizer. 520 const OptimizerSystem& getOptimizerSystem() const; 521 522 /// Indicate whether the Optimizer is currently set to use a numerical gradient. 523 bool isUsingNumericalGradient() const; 524 /// Indicate whether the Optimizer is currently set to use a numerical Jacobian. 525 bool isUsingNumericalJacobian() const; 526 /// Return the estimated accuracy last specified in useNumericalGradient(). 527 Real getEstimatedAccuracyOfObjective() const; 528 /// Return the estimated accuracy last specified in useNumericalJacobian(). 529 Real getEstimatedAccuracyOfConstraints() const; 530 531 // This is a local class. 532 class OptimizerRep; 533 private: 534 Optimizer( const Optimizer& c ); 535 Optimizer& operator=(const Optimizer& rhs); 536 537 OptimizerRep* constructOptimizerRep(const OptimizerSystem&, OptimizerAlgorithm); getRep()538 const OptimizerRep& getRep() const {assert(rep); return *rep;} updRep()539 OptimizerRep& updRep() {assert(rep); return *rep;} 540 541 // Hidden implementation to preserve binary compatibility. 542 OptimizerRep* rep; 543 544 friend class OptimizerRep; 545 }; // class Optimizer 546 547 } // namespace SimTK 548 549 #endif //SimTK_SIMMATH_OPTIMIZER_H_ 550 551