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: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 2009) $ 086 * @since 1.2 087 */ 088 public abstract class DirectSearchOptimizer implements MultivariateRealOptimizer { 089 090 /** Simplex. */ 091 protected RealPointValuePair[] simplex; 092 093 /** Objective function. */ 094 private MultivariateRealFunction f; 095 096 /** Convergence checker. */ 097 private RealConvergenceChecker checker; 098 099 /** Maximal number of iterations allowed. */ 100 private int maxIterations; 101 102 /** Number of iterations already performed. */ 103 private int iterations; 104 105 /** Maximal number of evaluations allowed. */ 106 private int maxEvaluations; 107 108 /** Number of evaluations already performed. */ 109 private int evaluations; 110 111 /** Start simplex configuration. */ 112 private double[][] startConfiguration; 113 114 /** Simple constructor. 115 */ 116 protected DirectSearchOptimizer() { 117 setConvergenceChecker(new SimpleScalarValueChecker()); 118 setMaxIterations(Integer.MAX_VALUE); 119 setMaxEvaluations(Integer.MAX_VALUE); 120 } 121 122 /** Set start configuration for simplex. 123 * <p>The start configuration for simplex is built from a box parallel to 124 * the canonical axes of the space. The simplex is the subset of vertices 125 * of a box parallel to the canonical axes. It is built as the path followed 126 * while traveling from one vertex of the box to the diagonally opposite 127 * vertex moving only along the box edges. The first vertex of the box will 128 * be located at the start point of the optimization.</p> 129 * <p>As an example, in dimension 3 a simplex has 4 vertices. Setting the 130 * steps to (1, 10, 2) and the start point to (1, 1, 1) would imply the 131 * start simplex would be: { (1, 1, 1), (2, 1, 1), (2, 11, 1), (2, 11, 3) }. 132 * The first vertex would be set to the start point at (1, 1, 1) and the 133 * last vertex would be set to the diagonally opposite vertex at (2, 11, 3).</p> 134 * @param steps steps along the canonical axes representing box edges, 135 * they may be negative but not null 136 * @exception IllegalArgumentException if one step is null 137 */ 138 public void setStartConfiguration(final double[] steps) 139 throws IllegalArgumentException { 140 // only the relative position of the n final vertices with respect 141 // to the first one are stored 142 final int n = steps.length; 143 startConfiguration = new double[n][n]; 144 for (int i = 0; i < n; ++i) { 145 final double[] vertexI = startConfiguration[i]; 146 for (int j = 0; j < i + 1; ++j) { 147 if (steps[j] == 0.0) { 148 throw MathRuntimeException.createIllegalArgumentException( 149 "equals vertices {0} and {1} in simplex configuration", 150 j, j + 1); 151 } 152 System.arraycopy(steps, 0, vertexI, 0, j + 1); 153 } 154 } 155 } 156 157 /** Set start configuration for simplex. 158 * <p>The real initial simplex will be set up by moving the reference 159 * simplex such that its first point is located at the start point of the 160 * optimization.</p> 161 * @param referenceSimplex reference simplex 162 * @exception IllegalArgumentException if the reference simplex does not 163 * contain at least one point, or if there is a dimension mismatch 164 * in the reference simplex or if one of its vertices is duplicated 165 */ 166 public void setStartConfiguration(final double[][] referenceSimplex) 167 throws IllegalArgumentException { 168 169 // only the relative position of the n final vertices with respect 170 // to the first one are stored 171 final int n = referenceSimplex.length - 1; 172 if (n < 0) { 173 throw MathRuntimeException.createIllegalArgumentException( 174 "simplex must contain at least one point"); 175 } 176 startConfiguration = new double[n][n]; 177 final double[] ref0 = referenceSimplex[0]; 178 179 // vertices loop 180 for (int i = 0; i < n + 1; ++i) { 181 182 final double[] refI = referenceSimplex[i]; 183 184 // safety checks 185 if (refI.length != n) { 186 throw MathRuntimeException.createIllegalArgumentException( 187 "dimension mismatch {0} != {1}", 188 refI.length, n); 189 } 190 for (int j = 0; j < i; ++j) { 191 final double[] refJ = referenceSimplex[j]; 192 boolean allEquals = true; 193 for (int k = 0; k < n; ++k) { 194 if (refI[k] != refJ[k]) { 195 allEquals = false; 196 break; 197 } 198 } 199 if (allEquals) { 200 throw MathRuntimeException.createIllegalArgumentException( 201 "equals vertices {0} and {1} in simplex configuration", 202 i, j); 203 } 204 } 205 206 // store vertex i position relative to vertex 0 position 207 if (i > 0) { 208 final double[] confI = startConfiguration[i - 1]; 209 for (int k = 0; k < n; ++k) { 210 confI[k] = refI[k] - ref0[k]; 211 } 212 } 213 214 } 215 216 } 217 218 /** {@inheritDoc} */ 219 public void setMaxIterations(int maxIterations) { 220 this.maxIterations = maxIterations; 221 } 222 223 /** {@inheritDoc} */ 224 public int getMaxIterations() { 225 return maxIterations; 226 } 227 228 /** {@inheritDoc} */ 229 public void setMaxEvaluations(int maxEvaluations) { 230 this.maxEvaluations = maxEvaluations; 231 } 232 233 /** {@inheritDoc} */ 234 public int getMaxEvaluations() { 235 return maxEvaluations; 236 } 237 238 /** {@inheritDoc} */ 239 public int getIterations() { 240 return iterations; 241 } 242 243 /** {@inheritDoc} */ 244 public int getEvaluations() { 245 return evaluations; 246 } 247 248 /** {@inheritDoc} */ 249 public void setConvergenceChecker(RealConvergenceChecker checker) { 250 this.checker = checker; 251 } 252 253 /** {@inheritDoc} */ 254 public RealConvergenceChecker getConvergenceChecker() { 255 return checker; 256 } 257 258 /** {@inheritDoc} */ 259 public RealPointValuePair optimize(final MultivariateRealFunction f, 260 final GoalType goalType, 261 final double[] startPoint) 262 throws FunctionEvaluationException, OptimizationException, 263 IllegalArgumentException { 264 265 if (startConfiguration == null) { 266 // no initial configuration has been set up for simplex 267 // build a default one from a unit hypercube 268 final double[] unit = new double[startPoint.length]; 269 Arrays.fill(unit, 1.0); 270 setStartConfiguration(unit); 271 } 272 273 this.f = f; 274 final Comparator<RealPointValuePair> comparator = 275 new Comparator<RealPointValuePair>() { 276 public int compare(final RealPointValuePair o1, 277 final RealPointValuePair o2) { 278 final double v1 = o1.getValue(); 279 final double v2 = o2.getValue(); 280 return (goalType == GoalType.MINIMIZE) ? 281 Double.compare(v1, v2) : Double.compare(v2, v1); 282 } 283 }; 284 285 // initialize search 286 iterations = 0; 287 evaluations = 0; 288 buildSimplex(startPoint); 289 evaluateSimplex(comparator); 290 291 RealPointValuePair[] previous = new RealPointValuePair[simplex.length]; 292 while (true) { 293 294 if (iterations > 0) { 295 boolean converged = true; 296 for (int i = 0; i < simplex.length; ++i) { 297 converged &= checker.converged(iterations, previous[i], simplex[i]); 298 } 299 if (converged) { 300 // we have found an optimum 301 return simplex[0]; 302 } 303 } 304 305 // we still need to search 306 System.arraycopy(simplex, 0, previous, 0, simplex.length); 307 iterateSimplex(comparator); 308 309 } 310 311 } 312 313 /** Increment the iterations counter by 1. 314 * @exception OptimizationException if the maximal number 315 * of iterations is exceeded 316 */ 317 protected void incrementIterationsCounter() 318 throws OptimizationException { 319 if (++iterations > maxIterations) { 320 throw new OptimizationException(new MaxIterationsExceededException(maxIterations)); 321 } 322 } 323 324 /** Compute the next simplex of the algorithm. 325 * @param comparator comparator to use to sort simplex vertices from best to worst 326 * @exception FunctionEvaluationException if the function cannot be evaluated at 327 * some point 328 * @exception OptimizationException if the algorithm fails to converge 329 * @exception IllegalArgumentException if the start point dimension is wrong 330 */ 331 protected abstract void iterateSimplex(final Comparator<RealPointValuePair> comparator) 332 throws FunctionEvaluationException, OptimizationException, IllegalArgumentException; 333 334 /** Evaluate the objective function on one point. 335 * <p>A side effect of this method is to count the number of 336 * function evaluations</p> 337 * @param x point on which the objective function should be evaluated 338 * @return objective function value at the given point 339 * @exception FunctionEvaluationException if no value can be computed for the 340 * parameters or if the maximal number of evaluations is exceeded 341 * @exception IllegalArgumentException if the start point dimension is wrong 342 */ 343 protected double evaluate(final double[] x) 344 throws FunctionEvaluationException, IllegalArgumentException { 345 if (++evaluations > maxEvaluations) { 346 throw new FunctionEvaluationException(new MaxEvaluationsExceededException(maxEvaluations), 347 x); 348 } 349 return f.value(x); 350 } 351 352 /** Build an initial simplex. 353 * @param startPoint the start point for optimization 354 * @exception IllegalArgumentException if the start point does not match 355 * simplex dimension 356 */ 357 private void buildSimplex(final double[] startPoint) 358 throws IllegalArgumentException { 359 360 final int n = startPoint.length; 361 if (n != startConfiguration.length) { 362 throw MathRuntimeException.createIllegalArgumentException( 363 "dimension mismatch {0} != {1}", 364 n, startConfiguration.length); 365 } 366 367 // set first vertex 368 simplex = new RealPointValuePair[n + 1]; 369 simplex[0] = new RealPointValuePair(startPoint, Double.NaN); 370 371 // set remaining vertices 372 for (int i = 0; i < n; ++i) { 373 final double[] confI = startConfiguration[i]; 374 final double[] vertexI = new double[n]; 375 for (int k = 0; k < n; ++k) { 376 vertexI[k] = startPoint[k] + confI[k]; 377 } 378 simplex[i + 1] = new RealPointValuePair(vertexI, Double.NaN); 379 } 380 381 } 382 383 /** Evaluate all the non-evaluated points of the simplex. 384 * @param comparator comparator to use to sort simplex vertices from best to worst 385 * @exception FunctionEvaluationException if no value can be computed for the parameters 386 * @exception OptimizationException if the maximal number of evaluations is exceeded 387 */ 388 protected void evaluateSimplex(final Comparator<RealPointValuePair> comparator) 389 throws FunctionEvaluationException, OptimizationException { 390 391 // evaluate the objective function at all non-evaluated simplex points 392 for (int i = 0; i < simplex.length; ++i) { 393 final RealPointValuePair vertex = simplex[i]; 394 final double[] point = vertex.getPointRef(); 395 if (Double.isNaN(vertex.getValue())) { 396 simplex[i] = new RealPointValuePair(point, evaluate(point), false); 397 } 398 } 399 400 // sort the simplex from best to worst 401 Arrays.sort(simplex, comparator); 402 403 } 404 405 /** Replace the worst point of the simplex by a new point. 406 * @param pointValuePair point to insert 407 * @param comparator comparator to use to sort simplex vertices from best to worst 408 */ 409 protected void replaceWorstPoint(RealPointValuePair pointValuePair, 410 final Comparator<RealPointValuePair> comparator) { 411 int n = simplex.length - 1; 412 for (int i = 0; i < n; ++i) { 413 if (comparator.compare(simplex[i], pointValuePair) > 0) { 414 RealPointValuePair tmp = simplex[i]; 415 simplex[i] = pointValuePair; 416 pointValuePair = tmp; 417 } 418 } 419 simplex[n] = pointValuePair; 420 } 421 422 }