View Javadoc

1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  
18  package org.apache.commons.math.optimization.direct;
19  
20  import java.util.Arrays;
21  import java.util.Comparator;
22  
23  import org.apache.commons.math.FunctionEvaluationException;
24  import org.apache.commons.math.MathRuntimeException;
25  import org.apache.commons.math.MaxEvaluationsExceededException;
26  import org.apache.commons.math.MaxIterationsExceededException;
27  import org.apache.commons.math.analysis.MultivariateRealFunction;
28  import org.apache.commons.math.optimization.GoalType;
29  import org.apache.commons.math.optimization.MultivariateRealOptimizer;
30  import org.apache.commons.math.optimization.OptimizationException;
31  import org.apache.commons.math.optimization.RealConvergenceChecker;
32  import org.apache.commons.math.optimization.RealPointValuePair;
33  import org.apache.commons.math.optimization.SimpleScalarValueChecker;
34  
35  /** 
36   * This class implements simplex-based direct search optimization
37   * algorithms.
38   *
39   * <p>Direct search methods only use objective function values, they don't
40   * need derivatives and don't either try to compute approximation of
41   * the derivatives. According to a 1996 paper by Margaret H. Wright
42   * (<a href="http://cm.bell-labs.com/cm/cs/doc/96/4-02.ps.gz">Direct
43   * Search Methods: Once Scorned, Now Respectable</a>), they are used
44   * when either the computation of the derivative is impossible (noisy
45   * functions, unpredictable discontinuities) or difficult (complexity,
46   * computation cost). In the first cases, rather than an optimum, a
47   * <em>not too bad</em> point is desired. In the latter cases, an
48   * optimum is desired but cannot be reasonably found. In all cases
49   * direct search methods can be useful.</p>
50   *
51   * <p>Simplex-based direct search methods are based on comparison of
52   * the objective function values at the vertices of a simplex (which is a
53   * set of n+1 points in dimension n) that is updated by the algorithms
54   * steps.<p>
55   *
56   * <p>The initial configuration of the simplex can be set using either
57   * {@link #setStartConfiguration(double[])} or {@link
58   * #setStartConfiguration(double[][])}. If neither method has been called
59   * before optimization is attempted, an explicit call to the first method
60   * with all steps set to +1 is triggered, thus building a default
61   * configuration from a unit hypercube. Each call to {@link
62   * #optimize(MultivariateRealFunction, GoalType, double[]) optimize} will reuse
63   * the current start configuration and move it such that its first vertex
64   * is at the provided start point of the optimization. If the same optimizer
65   * is used to solve different problems and the number of parameters change,
66   * the start configuration <em>must</em> be reset or a dimension mismatch
67   * will occur.</p>
68   *
69   * <p>If {@link #setConvergenceChecker(RealConvergenceChecker)} is not called,
70   * a default {@link SimpleScalarValueChecker} is used.</p>
71   *
72   * <p>Convergence is checked by providing the <em>worst</em> points of
73   * previous and current simplex to the convergence checker, not the best ones.</p>
74   *
75   * <p>This class is the base class performing the boilerplate simplex
76   * initialization and handling. The simplex update by itself is
77   * performed by the derived classes according to the implemented
78   * algorithms.</p>
79   *
80   * implements MultivariateRealOptimizer since 2.0
81   * 
82   * @see MultivariateRealFunction
83   * @see NelderMead
84   * @see MultiDirectional
85   * @version $Revision: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 2009) $
86   * @since 1.2
87   */
88  public abstract class DirectSearchOptimizer implements MultivariateRealOptimizer {
89  
90      /** Simplex. */
91      protected RealPointValuePair[] simplex;
92  
93      /** Objective function. */
94      private MultivariateRealFunction f;
95  
96      /** Convergence checker. */
97      private RealConvergenceChecker checker;
98  
99      /** 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 }