    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    }