001 /* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017 018 package org.apache.commons.math.optimization.direct; 019 020 import java.util.Arrays; 021 import java.util.Comparator; 022 023 import org.apache.commons.math.FunctionEvaluationException; 024 import org.apache.commons.math.MathRuntimeException; 025 import org.apache.commons.math.MaxEvaluationsExceededException; 026 import org.apache.commons.math.MaxIterationsExceededException; 027 import org.apache.commons.math.analysis.MultivariateRealFunction; 028 import org.apache.commons.math.optimization.GoalType; 029 import org.apache.commons.math.optimization.MultivariateRealOptimizer; 030 import org.apache.commons.math.optimization.OptimizationException; 031 import org.apache.commons.math.optimization.RealConvergenceChecker; 032 import org.apache.commons.math.optimization.RealPointValuePair; 033 import org.apache.commons.math.optimization.SimpleScalarValueChecker; 034 035 /** 036 * This class implements simplex-based direct search optimization 037 * algorithms. 038 * 039 * <p>Direct search methods only use objective function values, they don't 040 * need derivatives and don't either try to compute approximation of 041 * the derivatives. According to a 1996 paper by Margaret H. Wright 042 * (<a href="http://cm.bell-labs.com/cm/cs/doc/96/4-02.ps.gz">Direct 043 * Search Methods: Once Scorned, Now Respectable</a>), they are used 044 * when either the computation of the derivative is impossible (noisy 045 * functions, unpredictable discontinuities) or difficult (complexity, 046 * computation cost). In the first cases, rather than an optimum, a 047 * <em>not too bad</em> point is desired. In the latter cases, an 048 * optimum is desired but cannot be reasonably found. In all cases 049 * direct search methods can be useful.</p> 050 * 051 * <p>Simplex-based direct search methods are based on comparison of 052 * the objective function values at the vertices of a simplex (which is a 053 * set of n+1 points in dimension n) that is updated by the algorithms 054 * steps.<p> 055 * 056 * <p>The initial configuration of the simplex can be set using either 057 * {@link #setStartConfiguration(double[])} or {@link 058 * #setStartConfiguration(double[][])}. If neither method has been called 059 * before optimization is attempted, an explicit call to the first method 060 * with all steps set to +1 is triggered, thus building a default 061 * configuration from a unit hypercube. Each call to {@link 062 * #optimize(MultivariateRealFunction, GoalType, double[]) optimize} will reuse 063 * the current start configuration and move it such that its first vertex 064 * is at the provided start point of the optimization. If the same optimizer 065 * is used to solve different problems and the number of parameters change, 066 * the start configuration <em>must</em> be reset or a dimension mismatch 067 * will occur.</p> 068 * 069 * <p>If {@link #setConvergenceChecker(RealConvergenceChecker)} is not called, 070 * a default {@link SimpleScalarValueChecker} is used.</p> 071 * 072 * <p>Convergence is checked by providing the <em>worst</em> points of 073 * previous and current simplex to the convergence checker, not the best ones.</p> 074 * 075 * <p>This class is the base class performing the boilerplate simplex 076 * initialization and handling. The simplex update by itself is 077 * performed by the derived classes according to the implemented 078 * algorithms.</p> 079 * 080 * implements MultivariateRealOptimizer since 2.0 081 * 082 * @see MultivariateRealFunction 083 * @see NelderMead 084 * @see MultiDirectional 085 * @version $Revision: 885278 $ $Date: 2009-11-29 16:47:51 -0500 (Sun, 29 Nov 2009) $ 086 * @since 1.2 087 */ 088 public abstract class DirectSearchOptimizer implements MultivariateRealOptimizer { 089 090 /** Message for equal vertices. */ 091 private static final String EQUAL_VERTICES_MESSAGE = 092 "equal vertices {0} and {1} in simplex configuration"; 093 094 /** Message for dimension mismatch. */ 095 private static final String DIMENSION_MISMATCH_MESSAGE = 096 "dimension mismatch {0} != {1}"; 097 098 /** Simplex. */ 099 protected RealPointValuePair[] simplex; 100 101 /** Objective function. */ 102 private MultivariateRealFunction f; 103 104 /** Convergence checker. */ 105 private RealConvergenceChecker checker; 106 107 /** Maximal number of iterations allowed. */ 108 private int maxIterations; 109 110 /** Number of iterations already performed. */ 111 private int iterations; 112 113 /** Maximal number of evaluations allowed. */ 114 private int maxEvaluations; 115 116 /** Number of evaluations already performed. */ 117 private int evaluations; 118 119 /** Start simplex configuration. */ 120 private double[][] startConfiguration; 121 122 /** Simple constructor. 123 */ 124 protected DirectSearchOptimizer() { 125 setConvergenceChecker(new SimpleScalarValueChecker()); 126 setMaxIterations(Integer.MAX_VALUE); 127 setMaxEvaluations(Integer.MAX_VALUE); 128 } 129 130 /** Set start configuration for simplex. 131 * <p>The start configuration for simplex is built from a box parallel to 132 * the canonical axes of the space. The simplex is the subset of vertices 133 * of a box parallel to the canonical axes. It is built as the path followed 134 * while traveling from one vertex of the box to the diagonally opposite 135 * vertex moving only along the box edges. The first vertex of the box will 136 * be located at the start point of the optimization.</p> 137 * <p>As an example, in dimension 3 a simplex has 4 vertices. Setting the 138 * steps to (1, 10, 2) and the start point to (1, 1, 1) would imply the 139 * start simplex would be: { (1, 1, 1), (2, 1, 1), (2, 11, 1), (2, 11, 3) }. 140 * The first vertex would be set to the start point at (1, 1, 1) and the 141 * last vertex would be set to the diagonally opposite vertex at (2, 11, 3).</p> 142 * @param steps steps along the canonical axes representing box edges, 143 * they may be negative but not null 144 * @exception IllegalArgumentException if one step is null 145 */ 146 public void setStartConfiguration(final double[] steps) 147 throws IllegalArgumentException { 148 // only the relative position of the n final vertices with respect 149 // to the first one are stored 150 final int n = steps.length; 151 startConfiguration = new double[n][n]; 152 for (int i = 0; i < n; ++i) { 153 final double[] vertexI = startConfiguration[i]; 154 for (int j = 0; j < i + 1; ++j) { 155 if (steps[j] == 0.0) { 156 throw MathRuntimeException.createIllegalArgumentException( 157 EQUAL_VERTICES_MESSAGE, j, j + 1); 158 } 159 System.arraycopy(steps, 0, vertexI, 0, j + 1); 160 } 161 } 162 } 163 164 /** Set start configuration for simplex. 165 * <p>The real initial simplex will be set up by moving the reference 166 * simplex such that its first point is located at the start point of the 167 * optimization.</p> 168 * @param referenceSimplex reference simplex 169 * @exception IllegalArgumentException if the reference simplex does not 170 * contain at least one point, or if there is a dimension mismatch 171 * in the reference simplex or if one of its vertices is duplicated 172 */ 173 public void setStartConfiguration(final double[][] referenceSimplex) 174 throws IllegalArgumentException { 175 176 // only the relative position of the n final vertices with respect 177 // to the first one are stored 178 final int n = referenceSimplex.length - 1; 179 if (n < 0) { 180 throw MathRuntimeException.createIllegalArgumentException( 181 "simplex must contain at least one point"); 182 } 183 startConfiguration = new double[n][n]; 184 final double[] ref0 = referenceSimplex[0]; 185 186 // vertices loop 187 for (int i = 0; i < n + 1; ++i) { 188 189 final double[] refI = referenceSimplex[i]; 190 191 // safety checks 192 if (refI.length != n) { 193 throw MathRuntimeException.createIllegalArgumentException( 194 DIMENSION_MISMATCH_MESSAGE, refI.length, n); 195 } 196 for (int j = 0; j < i; ++j) { 197 final double[] refJ = referenceSimplex[j]; 198 boolean allEquals = true; 199 for (int k = 0; k < n; ++k) { 200 if (refI[k] != refJ[k]) { 201 allEquals = false; 202 break; 203 } 204 } 205 if (allEquals) { 206 throw MathRuntimeException.createIllegalArgumentException( 207 EQUAL_VERTICES_MESSAGE, i, j); 208 } 209 } 210 211 // store vertex i position relative to vertex 0 position 212 if (i > 0) { 213 final double[] confI = startConfiguration[i - 1]; 214 for (int k = 0; k < n; ++k) { 215 confI[k] = refI[k] - ref0[k]; 216 } 217 } 218 219 } 220 221 } 222 223 /** {@inheritDoc} */ 224 public void setMaxIterations(int maxIterations) { 225 this.maxIterations = maxIterations; 226 } 227 228 /** {@inheritDoc} */ 229 public int getMaxIterations() { 230 return maxIterations; 231 } 232 233 /** {@inheritDoc} */ 234 public void setMaxEvaluations(int maxEvaluations) { 235 this.maxEvaluations = maxEvaluations; 236 } 237 238 /** {@inheritDoc} */ 239 public int getMaxEvaluations() { 240 return maxEvaluations; 241 } 242 243 /** {@inheritDoc} */ 244 public int getIterations() { 245 return iterations; 246 } 247 248 /** {@inheritDoc} */ 249 public int getEvaluations() { 250 return evaluations; 251 } 252 253 /** {@inheritDoc} */ 254 public void setConvergenceChecker(RealConvergenceChecker convergenceChecker) { 255 this.checker = convergenceChecker; 256 } 257 258 /** {@inheritDoc} */ 259 public RealConvergenceChecker getConvergenceChecker() { 260 return checker; 261 } 262 263 /** {@inheritDoc} */ 264 public RealPointValuePair optimize(final MultivariateRealFunction function, 265 final GoalType goalType, 266 final double[] startPoint) 267 throws FunctionEvaluationException, OptimizationException, 268 IllegalArgumentException { 269 270 if (startConfiguration == null) { 271 // no initial configuration has been set up for simplex 272 // build a default one from a unit hypercube 273 final double[] unit = new double[startPoint.length]; 274 Arrays.fill(unit, 1.0); 275 setStartConfiguration(unit); 276 } 277 278 this.f = function; 279 final Comparator<RealPointValuePair> comparator = 280 new Comparator<RealPointValuePair>() { 281 public int compare(final RealPointValuePair o1, 282 final RealPointValuePair o2) { 283 final double v1 = o1.getValue(); 284 final double v2 = o2.getValue(); 285 return (goalType == GoalType.MINIMIZE) ? 286 Double.compare(v1, v2) : Double.compare(v2, v1); 287 } 288 }; 289 290 // initialize search 291 iterations = 0; 292 evaluations = 0; 293 buildSimplex(startPoint); 294 evaluateSimplex(comparator); 295 296 RealPointValuePair[] previous = new RealPointValuePair[simplex.length]; 297 while (true) { 298 299 if (iterations > 0) { 300 boolean converged = true; 301 for (int i = 0; i < simplex.length; ++i) { 302 converged &= checker.converged(iterations, previous[i], simplex[i]); 303 } 304 if (converged) { 305 // we have found an optimum 306 return simplex[0]; 307 } 308 } 309 310 // we still need to search 311 System.arraycopy(simplex, 0, previous, 0, simplex.length); 312 iterateSimplex(comparator); 313 314 } 315 316 } 317 318 /** Increment the iterations counter by 1. 319 * @exception OptimizationException if the maximal number 320 * of iterations is exceeded 321 */ 322 protected void incrementIterationsCounter() 323 throws OptimizationException { 324 if (++iterations > maxIterations) { 325 throw new OptimizationException(new MaxIterationsExceededException(maxIterations)); 326 } 327 } 328 329 /** Compute the next simplex of the algorithm. 330 * @param comparator comparator to use to sort simplex vertices from best to worst 331 * @exception FunctionEvaluationException if the function cannot be evaluated at 332 * some point 333 * @exception OptimizationException if the algorithm fails to converge 334 * @exception IllegalArgumentException if the start point dimension is wrong 335 */ 336 protected abstract void iterateSimplex(final Comparator<RealPointValuePair> comparator) 337 throws FunctionEvaluationException, OptimizationException, IllegalArgumentException; 338 339 /** Evaluate the objective function on one point. 340 * <p>A side effect of this method is to count the number of 341 * function evaluations</p> 342 * @param x point on which the objective function should be evaluated 343 * @return objective function value at the given point 344 * @exception FunctionEvaluationException if no value can be computed for the 345 * parameters or if the maximal number of evaluations is exceeded 346 * @exception IllegalArgumentException if the start point dimension is wrong 347 */ 348 protected double evaluate(final double[] x) 349 throws FunctionEvaluationException, IllegalArgumentException { 350 if (++evaluations > maxEvaluations) { 351 throw new FunctionEvaluationException(new MaxEvaluationsExceededException(maxEvaluations), 352 x); 353 } 354 return f.value(x); 355 } 356 357 /** Build an initial simplex. 358 * @param startPoint the start point for optimization 359 * @exception IllegalArgumentException if the start point does not match 360 * simplex dimension 361 */ 362 private void buildSimplex(final double[] startPoint) 363 throws IllegalArgumentException { 364 365 final int n = startPoint.length; 366 if (n != startConfiguration.length) { 367 throw MathRuntimeException.createIllegalArgumentException( 368 DIMENSION_MISMATCH_MESSAGE, n, startConfiguration.length); 369 } 370 371 // set first vertex 372 simplex = new RealPointValuePair[n + 1]; 373 simplex[0] = new RealPointValuePair(startPoint, Double.NaN); 374 375 // set remaining vertices 376 for (int i = 0; i < n; ++i) { 377 final double[] confI = startConfiguration[i]; 378 final double[] vertexI = new double[n]; 379 for (int k = 0; k < n; ++k) { 380 vertexI[k] = startPoint[k] + confI[k]; 381 } 382 simplex[i + 1] = new RealPointValuePair(vertexI, Double.NaN); 383 } 384 385 } 386 387 /** Evaluate all the non-evaluated points of the simplex. 388 * @param comparator comparator to use to sort simplex vertices from best to worst 389 * @exception FunctionEvaluationException if no value can be computed for the parameters 390 * @exception OptimizationException if the maximal number of evaluations is exceeded 391 */ 392 protected void evaluateSimplex(final Comparator<RealPointValuePair> comparator) 393 throws FunctionEvaluationException, OptimizationException { 394 395 // evaluate the objective function at all non-evaluated simplex points 396 for (int i = 0; i < simplex.length; ++i) { 397 final RealPointValuePair vertex = simplex[i]; 398 final double[] point = vertex.getPointRef(); 399 if (Double.isNaN(vertex.getValue())) { 400 simplex[i] = new RealPointValuePair(point, evaluate(point), false); 401 } 402 } 403 404 // sort the simplex from best to worst 405 Arrays.sort(simplex, comparator); 406 407 } 408 409 /** Replace the worst point of the simplex by a new point. 410 * @param pointValuePair point to insert 411 * @param comparator comparator to use to sort simplex vertices from best to worst 412 */ 413 protected void replaceWorstPoint(RealPointValuePair pointValuePair, 414 final Comparator<RealPointValuePair> comparator) { 415 int n = simplex.length - 1; 416 for (int i = 0; i < n; ++i) { 417 if (comparator.compare(simplex[i], pointValuePair) > 0) { 418 RealPointValuePair tmp = simplex[i]; 419 simplex[i] = pointValuePair; 420 pointValuePair = tmp; 421 } 422 } 423 simplex[n] = pointValuePair; 424 } 425 426 }