    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.exception.util.LocalizedFormats;
    029    import org.apache.commons.math.optimization.GoalType;
    030    import org.apache.commons.math.optimization.MultivariateRealOptimizer;
    031    import org.apache.commons.math.optimization.OptimizationException;
    032    import org.apache.commons.math.optimization.RealConvergenceChecker;
    033    import org.apache.commons.math.optimization.RealPointValuePair;
    034    import org.apache.commons.math.optimization.SimpleScalarValueChecker;
    035
    036    /**
    037     * This class implements simplex-based direct search optimization
    038     * algorithms.
    039     *
    040     * <p>Direct search methods only use objective function values, they don't
    041     * need derivatives and don't either try to compute approximation of
    042     * the derivatives. According to a 1996 paper by Margaret H. Wright
    043     * (<a href="http://cm.bell-labs.com/cm/cs/doc/96/4-02.ps.gz">Direct
    044     * Search Methods: Once Scorned, Now Respectable</a>), they are used
    045     * when either the computation of the derivative is impossible (noisy
    046     * functions, unpredictable discontinuities) or difficult (complexity,
    047     * computation cost). In the first cases, rather than an optimum, a
    048     * <em>not too bad</em> point is desired. In the latter cases, an
    049     * optimum is desired but cannot be reasonably found. In all cases
    050     * direct search methods can be useful.</p>
    051     *
    052     * <p>Simplex-based direct search methods are based on comparison of
    053     * the objective function values at the vertices of a simplex (which is a
    054     * set of n+1 points in dimension n) that is updated by the algorithms
    055     * steps.<p>
    056     *
    057     * <p>The initial configuration of the simplex can be set using either
    058     * {@link #setStartConfiguration(double[])} or {@link
    059     * #setStartConfiguration(double[][])}. If neither method has been called
    060     * before optimization is attempted, an explicit call to the first method
    061     * with all steps set to +1 is triggered, thus building a default
    062     * configuration from a unit hypercube. Each call to {@link
    063     * #optimize(MultivariateRealFunction, GoalType, double[]) optimize} will reuse
    064     * the current start configuration and move it such that its first vertex
    065     * is at the provided start point of the optimization. If the {@code optimize}
    066     * method is called to solve a different problem and the number of parameters
    067     * change, the start configuration will be reset to a default one with the
    068     * appropriate dimensions.</p>
    069     *
    070     * <p>If {@link #setConvergenceChecker(RealConvergenceChecker)} is not called,
    071     * a default {@link SimpleScalarValueChecker} is used.</p>
    072     *
    073     * <p>Convergence is checked by providing the <em>worst</em> points of
    074     * previous and current simplex to the convergence checker, not the best ones.</p>
    075     *
    076     * <p>This class is the base class performing the boilerplate simplex
    077     * initialization and handling. The simplex update by itself is
    078     * performed by the derived classes according to the implemented
    079     * algorithms.</p>
    080     *
    081     * implements MultivariateRealOptimizer since 2.0
    082     *
    083     * @see MultivariateRealFunction
    084     * @see NelderMead
    085     * @see MultiDirectional
    086     * @version $Revision: 1070725 $ $Date: 2011-02-15 02:31:12 +0100 (mar. 15 f??vr. 2011) $
    087     * @since 1.2
    088     */
    089    public abstract class DirectSearchOptimizer implements MultivariateRealOptimizer {
    090
    091        /** Simplex. */
    092        protected RealPointValuePair[] simplex;
    093
    094        /** Objective function. */
    095        private MultivariateRealFunction f;
    096
    097        /** Convergence checker. */
    098        private RealConvergenceChecker checker;
    099
    100        /** Maximal number of iterations allowed. */
    101        private int maxIterations;
    102
    103        /** Number of iterations already performed. */
    104        private int iterations;
    105
    106        /** Maximal number of evaluations allowed. */
    107        private int maxEvaluations;
    108
    109        /** Number of evaluations already performed. */
    110        private int evaluations;
    111
    112        /** Start simplex configuration. */
    113        private double[][] startConfiguration;
    114
    115        /** Simple constructor.
    116         */
    117        protected DirectSearchOptimizer() {
    118            setConvergenceChecker(new SimpleScalarValueChecker());
    119            setMaxIterations(Integer.MAX_VALUE);
    120            setMaxEvaluations(Integer.MAX_VALUE);
    121        }
    122
    123        /** Set start configuration for simplex.
    124         * <p>The start configuration for simplex is built from a box parallel to
    125         * the canonical axes of the space. The simplex is the subset of vertices
    126         * of a box parallel to the canonical axes. It is built as the path followed
    127         * while traveling from one vertex of the box to the diagonally opposite
    128         * vertex moving only along the box edges. The first vertex of the box will
    129         * be located at the start point of the optimization.</p>
    130         * <p>As an example, in dimension 3 a simplex has 4 vertices. Setting the
    131         * steps to (1, 10, 2) and the start point to (1, 1, 1) would imply the
    132         * start simplex would be: { (1, 1, 1), (2, 1, 1), (2, 11, 1), (2, 11, 3) }.
    133         * The first vertex would be set to the start point at (1, 1, 1) and the
    134         * last vertex would be set to the diagonally opposite vertex at (2, 11, 3).</p>
    135         * @param steps steps along the canonical axes representing box edges,
    136         * they may be negative but not null
    137         * @exception IllegalArgumentException if one step is null
    138         */
    139        public void setStartConfiguration(final double[] steps)
    140            throws IllegalArgumentException {
    141            // only the relative position of the n final vertices with respect
    142            // to the first one are stored
    143            final int n = steps.length;
    144            startConfiguration = new double[n][n];
    145            for (int i = 0; i < n; ++i) {
    146                final double[] vertexI = startConfiguration[i];
    147                for (int j = 0; j < i + 1; ++j) {
    148                    if (steps[j] == 0.0) {
    149                        throw MathRuntimeException.createIllegalArgumentException(
    150                              LocalizedFormats.EQUAL_VERTICES_IN_SIMPLEX, 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                        LocalizedFormats.SIMPLEX_NEED_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                          LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, refI.length, n);
    188                }
    189                for (int j = 0; j < i; ++j) {
    190                    final double[] refJ = referenceSimplex[j];
    191                    boolean allEquals = true;
    192                    for (int k = 0; k < n; ++k) {
    193                        if (refI[k] != refJ[k]) {
    194                            allEquals = false;
    195                            break;
    196                        }
    197                    }
    198                    if (allEquals) {
    199                        throw MathRuntimeException.createIllegalArgumentException(
    200                              LocalizedFormats.EQUAL_VERTICES_IN_SIMPLEX, i, j);
    201                    }
    202                }
    203
    204                // store vertex i position relative to vertex 0 position
    205                if (i > 0) {
    206                    final double[] confI = startConfiguration[i - 1];
    207                    for (int k = 0; k < n; ++k) {
    208                        confI[k] = refI[k] - ref0[k];
    209                    }
    210                }
    211
    212            }
    213
    214        }
    215
    216        /** {@inheritDoc} */
    217        public void setMaxIterations(int maxIterations) {
    218            this.maxIterations = maxIterations;
    219        }
    220
    221        /** {@inheritDoc} */
    222        public int getMaxIterations() {
    223            return maxIterations;
    224        }
    225
    226        /** {@inheritDoc} */
    227        public void setMaxEvaluations(int maxEvaluations) {
    228            this.maxEvaluations = maxEvaluations;
    229        }
    230
    231        /** {@inheritDoc} */
    232        public int getMaxEvaluations() {
    233            return maxEvaluations;
    234        }
    235
    236        /** {@inheritDoc} */
    237        public int getIterations() {
    238            return iterations;
    239        }
    240
    241        /** {@inheritDoc} */
    242        public int getEvaluations() {
    243            return evaluations;
    244        }
    245
    246        /** {@inheritDoc} */
    247        public void setConvergenceChecker(RealConvergenceChecker convergenceChecker) {
    248            this.checker = convergenceChecker;
    249        }
    250
    251        /** {@inheritDoc} */
    252        public RealConvergenceChecker getConvergenceChecker() {
    253            return checker;
    254        }
    255
    256        /** {@inheritDoc} */
    257        public RealPointValuePair optimize(final MultivariateRealFunction function,
    258                                           final GoalType goalType,
    259                                           final double[] startPoint)
    260            throws FunctionEvaluationException, OptimizationException, IllegalArgumentException {
    261
    262            if ((startConfiguration == null) ||
    263                (startConfiguration.length != startPoint.length)) {
    264                // no initial configuration has been set up for simplex
    265                // build a default one from a unit hypercube
    266                final double[] unit = new double[startPoint.length];
    267                Arrays.fill(unit, 1.0);
    268                setStartConfiguration(unit);
    269            }
    270
    271            this.f = function;
    272            final Comparator<RealPointValuePair> comparator =
    273                new Comparator<RealPointValuePair>() {
    274                    public int compare(final RealPointValuePair o1,
    275                                       final RealPointValuePair o2) {
    276                        final double v1 = o1.getValue();
    277                        final double v2 = o2.getValue();
    278                        return (goalType == GoalType.MINIMIZE) ?
    279                                Double.compare(v1, v2) : Double.compare(v2, v1);
    280                    }
    281                };
    282
    283            // initialize search
    284            iterations  = 0;
    285            evaluations = 0;
    286            buildSimplex(startPoint);
    287            evaluateSimplex(comparator);
    288
    289            RealPointValuePair[] previous = new RealPointValuePair[simplex.length];
    290            while (true) {
    291
    292                if (iterations > 0) {
    293                    boolean converged = true;
    294                    for (int i = 0; i < simplex.length; ++i) {
    295                        converged &= checker.converged(iterations, previous[i], simplex[i]);
    296                    }
    297                    if (converged) {
    298                        // we have found an optimum
    299                        return simplex[0];
    300                    }
    301                }
    302
    303                // we still need to search
    304                System.arraycopy(simplex, 0, previous, 0, simplex.length);
    305                iterateSimplex(comparator);
    306
    307            }
    308
    309        }
    310
    311        /** Increment the iterations counter by 1.
    312         * @exception OptimizationException if the maximal number
    313         * of iterations is exceeded
    314         */
    315        protected void incrementIterationsCounter()
    316            throws OptimizationException {
    317            if (++iterations > maxIterations) {
    318                throw new OptimizationException(new MaxIterationsExceededException(maxIterations));
    319            }
    320        }
    321
    322        /** Compute the next simplex of the algorithm.
    323         * @param comparator comparator to use to sort simplex vertices from best to worst
    324         * @exception FunctionEvaluationException if the function cannot be evaluated at
    325         * some point
    326         * @exception OptimizationException if the algorithm fails to converge
    327         * @exception IllegalArgumentException if the start point dimension is wrong
    328         */
    329        protected abstract void iterateSimplex(final Comparator<RealPointValuePair> comparator)
    330            throws FunctionEvaluationException, OptimizationException, IllegalArgumentException;
    331
    332        /** Evaluate the objective function on one point.
    333         * <p>A side effect of this method is to count the number of
    334         * function evaluations</p>
    335         * @param x point on which the objective function should be evaluated
    336         * @return objective function value at the given point
    337         * @exception FunctionEvaluationException if no value can be computed for the
    338         * parameters or if the maximal number of evaluations is exceeded
    339         * @exception IllegalArgumentException if the start point dimension is wrong
    340         */
    341        protected double evaluate(final double[] x)
    342            throws FunctionEvaluationException, IllegalArgumentException {
    343            if (++evaluations > maxEvaluations) {
    344                throw new FunctionEvaluationException(new MaxEvaluationsExceededException(maxEvaluations), x);
    345            }
    346            return f.value(x);
    347        }
    348
    349        /** Build an initial simplex.
    350         * @param startPoint the start point for optimization
    351         * @exception IllegalArgumentException if the start point does not match
    352         * simplex dimension
    353         */
    354        private void buildSimplex(final double[] startPoint)
    355            throws IllegalArgumentException {
    356
    357            final int n = startPoint.length;
    358            if (n != startConfiguration.length) {
    359                throw MathRuntimeException.createIllegalArgumentException(
    360                      LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, n, startConfiguration.length);
    361            }
    362
    363            // set first vertex
    364            simplex = new RealPointValuePair[n + 1];
    365            simplex[0] = new RealPointValuePair(startPoint, Double.NaN);
    366
    367            // set remaining vertices
    368            for (int i = 0; i < n; ++i) {
    369                final double[] confI   = startConfiguration[i];
    370                final double[] vertexI = new double[n];
    371                for (int k = 0; k < n; ++k) {
    372                    vertexI[k] = startPoint[k] + confI[k];
    373                }
    374                simplex[i + 1] = new RealPointValuePair(vertexI, Double.NaN);
    375            }
    376
    377        }
    378
    379        /** Evaluate all the non-evaluated points of the simplex.
    380         * @param comparator comparator to use to sort simplex vertices from best to worst
    381         * @exception FunctionEvaluationException if no value can be computed for the parameters
    382         * @exception OptimizationException if the maximal number of evaluations is exceeded
    383         */
    384        protected void evaluateSimplex(final Comparator<RealPointValuePair> comparator)
    385            throws FunctionEvaluationException, OptimizationException {
    386
    387            // evaluate the objective function at all non-evaluated simplex points
    388            for (int i = 0; i < simplex.length; ++i) {
    389                final RealPointValuePair vertex = simplex[i];
    390                final double[] point = vertex.getPointRef();
    391                if (Double.isNaN(vertex.getValue())) {
    392                    simplex[i] = new RealPointValuePair(point, evaluate(point), false);
    393                }
    394            }
    395
    396            // sort the simplex from best to worst
    397            Arrays.sort(simplex, comparator);
    398
    399        }
    400
    401        /** Replace the worst point of the simplex by a new point.
    402         * @param pointValuePair point to insert
    403         * @param comparator comparator to use to sort simplex vertices from best to worst
    404         */
    405        protected void replaceWorstPoint(RealPointValuePair pointValuePair,
    406                                         final Comparator<RealPointValuePair> comparator) {
    407            int n = simplex.length - 1;
    408            for (int i = 0; i < n; ++i) {
    409                if (comparator.compare(simplex[i], pointValuePair) > 0) {
    410                    RealPointValuePair tmp = simplex[i];
    411                    simplex[i]         = pointValuePair;
    412                    pointValuePair     = tmp;
    413                }
    414            }
    415            simplex[n] = pointValuePair;
    416        }
    417
    418    }