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;
019    
020    import org.apache.commons.math.ConvergenceException;
021    import org.apache.commons.math.FunctionEvaluationException;
022    import org.apache.commons.math.MathRuntimeException;
023    import org.apache.commons.math.analysis.UnivariateRealFunction;
024    import org.apache.commons.math.random.RandomGenerator;
025    
026    /**
027     * Special implementation of the {@link UnivariateRealOptimizer} interface adding
028     * multi-start features to an existing optimizer.
029     * <p>
030     * This class wraps a classical optimizer to use it several times in
031     * turn with different starting points in order to avoid being trapped
032     * into a local extremum when looking for a global one.
033     * </p>
034     * @version $Revision: 811685 $ $Date: 2009-09-05 13:36:48 -0400 (Sat, 05 Sep 2009) $
035     * @since 2.0
036     */
037    public class MultiStartUnivariateRealOptimizer implements UnivariateRealOptimizer {
038    
039        /** Serializable version identifier. */
040        private static final long serialVersionUID = 5983375963110961019L;
041    
042        /** Underlying classical optimizer. */
043        private final UnivariateRealOptimizer optimizer;
044    
045        /** Maximal number of iterations allowed. */
046        private int maxIterations;
047    
048        /** Maximal number of evaluations allowed. */
049        private int maxEvaluations;
050    
051        /** Number of iterations already performed for all starts. */
052        private int totalIterations;
053    
054        /** Number of evaluations already performed for all starts. */
055        private int totalEvaluations;
056    
057        /** Number of starts to go. */
058        private int starts;
059    
060        /** Random generator for multi-start. */
061        private RandomGenerator generator;
062    
063        /** Found optima. */
064        private double[] optima;
065    
066        /** Found function values at optima. */
067        private double[] optimaValues;
068    
069        /**
070         * Create a multi-start optimizer from a single-start optimizer
071         * @param optimizer single-start optimizer to wrap
072         * @param starts number of starts to perform (including the
073         * first one), multi-start is disabled if value is less than or
074         * equal to 1
075         * @param generator random generator to use for restarts
076         */
077        public MultiStartUnivariateRealOptimizer(final UnivariateRealOptimizer optimizer,
078                                                 final int starts,
079                                                 final RandomGenerator generator) {
080            this.optimizer        = optimizer;
081            this.totalIterations  = 0;
082            this.starts           = starts;
083            this.generator        = generator;
084            this.optima           = null;
085            setMaximalIterationCount(Integer.MAX_VALUE);
086            setMaxEvaluations(Integer.MAX_VALUE);
087        }
088    
089        /** {@inheritDoc} */
090        public double getFunctionValue() {
091            return optimizer.getFunctionValue();
092        }
093    
094        /** {@inheritDoc} */
095        public double getResult() {
096            return optimizer.getResult();
097        }
098    
099        /** {@inheritDoc} */
100        public double getAbsoluteAccuracy() {
101            return optimizer.getAbsoluteAccuracy();
102        }
103    
104        /** {@inheritDoc} */
105        public int getIterationCount() {
106            return totalIterations;
107        }
108    
109        /** {@inheritDoc} */
110        public int getMaximalIterationCount() {
111            return maxIterations;
112        }
113    
114        /** {@inheritDoc} */
115        public int getMaxEvaluations() {
116            return maxEvaluations;
117        }
118    
119        /** {@inheritDoc} */
120        public int getEvaluations() {
121            return totalEvaluations;
122        }
123    
124        /** {@inheritDoc} */
125        public double getRelativeAccuracy() {
126            return optimizer.getRelativeAccuracy();
127        }
128    
129        /** {@inheritDoc} */
130        public void resetAbsoluteAccuracy() {
131            optimizer.resetAbsoluteAccuracy();
132        }
133    
134        /** {@inheritDoc} */
135        public void resetMaximalIterationCount() {
136            optimizer.resetMaximalIterationCount();
137        }
138    
139        /** {@inheritDoc} */
140        public void resetRelativeAccuracy() {
141            optimizer.resetRelativeAccuracy();
142        }
143    
144        /** {@inheritDoc} */
145        public void setAbsoluteAccuracy(double accuracy) {
146            optimizer.setAbsoluteAccuracy(accuracy);
147        }
148    
149        /** {@inheritDoc} */
150        public void setMaximalIterationCount(int count) {
151            this.maxIterations = count;
152        }
153    
154        /** {@inheritDoc} */
155        public void setMaxEvaluations(int maxEvaluations) {
156            this.maxEvaluations = maxEvaluations;
157        }
158    
159        /** {@inheritDoc} */
160        public void setRelativeAccuracy(double accuracy) {
161            optimizer.setRelativeAccuracy(accuracy);
162        }
163    
164        /** Get all the optima found during the last call to {@link
165         * #optimize(UnivariateRealFunction, GoalType, double, double) optimize}.
166         * <p>The optimizer stores all the optima found during a set of
167         * restarts. The {@link #optimize(UnivariateRealFunction, GoalType,
168         * double, double) optimize} method returns the best point only. This
169         * method returns all the points found at the end of each starts,
170         * including the best one already returned by the {@link
171         * #optimize(UnivariateRealFunction, GoalType, double, double) optimize}
172         * method.
173         * </p>
174         * <p>
175         * The returned array as one element for each start as specified
176         * in the constructor. It is ordered with the results from the
177         * runs that did converge first, sorted from best to worst
178         * objective value (i.e in ascending order if minimizing and in
179         * descending order if maximizing), followed by Double.NaN elements
180         * corresponding to the runs that did not converge. This means all
181         * elements will be NaN if the {@link #optimize(UnivariateRealFunction,
182         * GoalType, double, double) optimize} method did throw a {@link
183         * ConvergenceException ConvergenceException}). This also means that
184         * if the first element is not NaN, it is the best point found across
185         * all starts.</p>
186         * @return array containing the optima
187         * @exception IllegalStateException if {@link #optimize(UnivariateRealFunction,
188         * GoalType, double, double) optimize} has not been called
189         * @see #getOptimaValues()
190         */
191        public double[] getOptima() throws IllegalStateException {
192            if (optima == null) {
193                throw MathRuntimeException.createIllegalStateException("no optimum computed yet");
194            }
195            return optima.clone();
196        }
197    
198        /** Get all the function values at optima found during the last call to {@link
199         * #optimize(UnivariateRealFunction, GoalType, double, double) optimize}.
200         * <p>
201         * The returned array as one element for each start as specified
202         * in the constructor. It is ordered with the results from the
203         * runs that did converge first, sorted from best to worst
204         * objective value (i.e in ascending order if minimizing and in
205         * descending order if maximizing), followed by Double.NaN elements
206         * corresponding to the runs that did not converge. This means all
207         * elements will be NaN if the {@link #optimize(UnivariateRealFunction,
208         * GoalType, double, double) optimize} method did throw a {@link
209         * ConvergenceException ConvergenceException}). This also means that
210         * if the first element is not NaN, it is the best point found across
211         * all starts.</p>
212         * @return array containing the optima
213         * @exception IllegalStateException if {@link #optimize(UnivariateRealFunction,
214         * GoalType, double, double) optimize} has not been called
215         * @see #getOptima()
216         */
217        public double[] getOptimaValues() throws IllegalStateException {
218            if (optimaValues == null) {
219                throw MathRuntimeException.createIllegalStateException("no optimum computed yet");
220            }
221            return optimaValues.clone();
222        }
223    
224        /** {@inheritDoc} */
225        public double optimize(final UnivariateRealFunction f, final GoalType goalType,
226                               final double min, final double max)
227            throws ConvergenceException,
228                FunctionEvaluationException {
229    
230            optima           = new double[starts];
231            optimaValues     = new double[starts];
232            totalIterations  = 0;
233            totalEvaluations = 0;
234    
235            // multi-start loop
236            for (int i = 0; i < starts; ++i) {
237    
238                try {
239                    optimizer.setMaximalIterationCount(maxIterations - totalIterations);
240                    optimizer.setMaxEvaluations(maxEvaluations - totalEvaluations);
241                    final double bound1 = (i == 0) ? min : min + generator.nextDouble() * (max - min);
242                    final double bound2 = (i == 0) ? max : min + generator.nextDouble() * (max - min);
243                    optima[i]       = optimizer.optimize(f, goalType,
244                                                         Math.min(bound1, bound2),
245                                                         Math.max(bound1, bound2));
246                    optimaValues[i] = optimizer.getFunctionValue();
247                } catch (FunctionEvaluationException fee) {
248                    optima[i]       = Double.NaN;
249                    optimaValues[i] = Double.NaN;
250                } catch (ConvergenceException ce) {
251                    optima[i]       = Double.NaN;
252                    optimaValues[i] = Double.NaN;
253                }
254    
255                totalIterations  += optimizer.getIterationCount();
256                totalEvaluations += optimizer.getEvaluations();
257    
258            }
259    
260            // sort the optima from best to worst, followed by NaN elements
261            int lastNaN = optima.length;
262            for (int i = 0; i < lastNaN; ++i) {
263                if (Double.isNaN(optima[i])) {
264                    optima[i] = optima[--lastNaN];
265                    optima[lastNaN + 1] = Double.NaN;
266                    optimaValues[i] = optimaValues[--lastNaN];
267                    optimaValues[lastNaN + 1] = Double.NaN;
268                }
269            }
270    
271            double currX = optima[0];
272            double currY = optimaValues[0];
273            for (int j = 1; j < lastNaN; ++j) {
274                final double prevY = currY;
275                currX = optima[j];
276                currY = optimaValues[j];
277                if ((goalType == GoalType.MAXIMIZE) ^ (currY < prevY)) {
278                    // the current element should be inserted closer to the beginning
279                    int i = j - 1;
280                    double mIX = optima[i];
281                    double mIY = optimaValues[i];
282                    while ((i >= 0) && ((goalType == GoalType.MAXIMIZE) ^ (currY < mIY))) {
283                        optima[i + 1]       = mIX;
284                        optimaValues[i + 1] = mIY;
285                        if (i-- != 0) {
286                            mIX = optima[i];
287                            mIY = optimaValues[i];
288                        } else {
289                            mIX = Double.NaN;
290                            mIY = Double.NaN;
291                        }
292                    }
293                    optima[i + 1]       = currX;
294                    optimaValues[i + 1] = currY;
295                    currX = optima[j];
296                    currY = optimaValues[j];
297                }
298            }
299    
300            if (Double.isNaN(optima[0])) {
301                throw new OptimizationException(
302                        "none of the {0} start points lead to convergence",
303                        starts);
304            }
305    
306            // return the found point given the best objective function value
307            return optima[0];
308    
309        }
310    
311        /** {@inheritDoc} */
312        public double optimize(final UnivariateRealFunction f, final GoalType goalType,
313                               final double min, final double max, final double startValue)
314                throws ConvergenceException, FunctionEvaluationException {
315            return optimize(f, goalType, min, max);
316        }
317    
318    }