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 }