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    }