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 }