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