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