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    }