    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 org.apache.commons.math.MaxIterationsExceededException;
    021    import org.apache.commons.math.analysis.UnivariateRealFunction;
    022    import org.apache.commons.math.FunctionEvaluationException;
    023    import org.apache.commons.math.optimization.GoalType;
    024    import org.apache.commons.math.optimization.OptimizationException;
    025    import org.apache.commons.math.optimization.RealPointValuePair;
    026    import org.apache.commons.math.optimization.general.AbstractScalarDifferentiableOptimizer;
    027    import org.apache.commons.math.optimization.univariate.AbstractUnivariateRealOptimizer;
    028    import org.apache.commons.math.optimization.univariate.BracketFinder;
    029    import org.apache.commons.math.optimization.univariate.BrentOptimizer;
    030
    031    /**
    032     * Powell algorithm.
    033     * This code is translated and adapted from the Python version of this
    034     * algorithm (as implemented in module {@code optimize.py} v0.5 of
    035     * <em>SciPy</em>).
    036     *
    037     * @version $Revision$ $Date$
    038     * @since 2.2
    039     */
    040    public class PowellOptimizer
    041        extends AbstractScalarDifferentiableOptimizer {
    042        /**
    043         * Default relative tolerance for line search ({@value}).
    044         */
    045        public static final double DEFAULT_LS_RELATIVE_TOLERANCE = 1e-7;
    046        /**
    047         * Default absolute tolerance for line search ({@value}).
    048         */
    049        public static final double DEFAULT_LS_ABSOLUTE_TOLERANCE = 1e-11;
    050        /**
    051         * Line search.
    052         */
    053        private final LineSearch line;
    054
    055        /**
    056         * Constructor with default line search tolerances (see the
    057         * {@link #PowellOptimizer(double,double) other constructor}).
    058         */
    059        public PowellOptimizer() {
    060            this(DEFAULT_LS_RELATIVE_TOLERANCE,
    061                 DEFAULT_LS_ABSOLUTE_TOLERANCE);
    062        }
    063
    064        /**
    065         * Constructor with default absolute line search tolerances (see
    066         * the {@link #PowellOptimizer(double,double) other constructor}).
    067         *
    068         * @param lsRelativeTolerance Relative error tolerance for
    069         * the line search algorithm ({@link BrentOptimizer}).
    070         */
    071        public PowellOptimizer(double lsRelativeTolerance) {
    072            this(lsRelativeTolerance,
    073                 DEFAULT_LS_ABSOLUTE_TOLERANCE);
    074        }
    075
    076        /**
    077         * @param lsRelativeTolerance Relative error tolerance for
    078         * the line search algorithm ({@link BrentOptimizer}).
    079         * @param lsAbsoluteTolerance Relative error tolerance for
    080         * the line search algorithm ({@link BrentOptimizer}).
    081         */
    082        public PowellOptimizer(double lsRelativeTolerance,
    083                               double lsAbsoluteTolerance) {
    084            line = new LineSearch(lsRelativeTolerance,
    085                                  lsAbsoluteTolerance);
    086        }
    087
    088        /** {@inheritDoc} */
    089        @Override
    090        protected RealPointValuePair doOptimize()
    091            throws FunctionEvaluationException,
    092                   OptimizationException {
    093            final double[] guess = point.clone();
    094            final int n = guess.length;
    095
    096            final double[][] direc = new double[n][n];
    097            for (int i = 0; i < n; i++) {
    098                direc[i][i] = 1;
    099            }
    100
    101            double[] x = guess;
    102            double fVal = computeObjectiveValue(x);
    103            double[] x1 = x.clone();
    104            while (true) {
    105                incrementIterationsCounter();
    106
    107                double fX = fVal;
    108                double fX2 = 0;
    109                double delta = 0;
    110                int bigInd = 0;
    111                double alphaMin = 0;
    112
    113                for (int i = 0; i < n; i++) {
    114                    final double[] d = /* Arrays.*/ copyOf(direc[i], n); // Java 1.5 does not support Arrays.copyOf()
    115
    116                    fX2 = fVal;
    117
    118                    line.search(x, d);
    119                    fVal = line.getValueAtOptimum();
    120                    alphaMin = line.getOptimum();
    121                    final double[][] result = newPointAndDirection(x, d, alphaMin);
    122                    x = result[0];
    123
    124                    if ((fX2 - fVal) > delta) {
    125                        delta = fX2 - fVal;
    126                        bigInd = i;
    127                    }
    128                }
    129
    130                final RealPointValuePair previous = new RealPointValuePair(x1, fX);
    131                final RealPointValuePair current = new RealPointValuePair(x, fVal);
    132                if (getConvergenceChecker().converged(getIterations(), previous, current)) {
    133                    if (goal == GoalType.MINIMIZE) {
    134                        return (fVal < fX) ? current : previous;
    135                    } else {
    136                        return (fVal > fX) ? current : previous;
    137                    }
    138                }
    139
    140                final double[] d = new double[n];
    141                final double[] x2 = new double[n];
    142                for (int i = 0; i < n; i++) {
    143                    d[i] = x[i] - x1[i];
    144                    x2[i] = 2 * x[i] - x1[i];
    145                }
    146
    147                x1 = x.clone();
    148                fX2 = computeObjectiveValue(x2);
    149
    150                if (fX > fX2) {
    151                    double t = 2 * (fX + fX2 - 2 * fVal);
    152                    double temp = fX - fVal - delta;
    153                    t *= temp * temp;
    154                    temp = fX - fX2;
    155                    t -= delta * temp * temp;
    156
    157                    if (t < 0.0) {
    158                        line.search(x, d);
    159                        fVal = line.getValueAtOptimum();
    160                        alphaMin = line.getOptimum();
    161                        final double[][] result = newPointAndDirection(x, d, alphaMin);
    162                        x = result[0];
    163
    164                        final int lastInd = n - 1;
    165                        direc[bigInd] = direc[lastInd];
    166                        direc[lastInd] = result[1];
    167                    }
    168                }
    169            }
    170        }
    171
    172        /**
    173         * Compute a new point (in the original space) and a new direction
    174         * vector, resulting from the line search.
    175         * The parameters {@code p} and {@code d} will be changed in-place.
    176         *
    177         * @param p Point used in the line search.
    178         * @param d Direction used in the line search.
    179         * @param optimum Optimum found by the line search.
    180         * @return a 2-element array containing the new point (at index 0) and
    181         * the new direction (at index 1).
    182         */
    183        private double[][] newPointAndDirection(double[] p,
    184                                                double[] d,
    185                                                double optimum) {
    186            final int n = p.length;
    187            final double[][] result = new double[2][n];
    188            final double[] nP = result[0];
    189            final double[] nD = result[1];
    190            for (int i = 0; i < n; i++) {
    191                nD[i] = d[i] * optimum;
    192                nP[i] = p[i] + nD[i];
    193            }
    194            return result;
    195        }
    196
    197        /**
    198         * Class for finding the minimum of the objective function along a given
    199         * direction.
    200         */
    201        private class LineSearch {
    202            /**
    203             * Optimizer.
    204             */
    205            private final AbstractUnivariateRealOptimizer optim = new BrentOptimizer();
    206            /**
    207             * Automatic bracketing.
    208             */
    209            private final BracketFinder bracket = new BracketFinder();
    210            /**
    211             * Value of the optimum.
    212             */
    213            private double optimum = Double.NaN;
    214            /**
    215             * Value of the objective function at the optimum.
    216             */
    217            private double valueAtOptimum = Double.NaN;
    218
    219            /**
    220             * @param relativeTolerance Relative tolerance.
    221             * @param absoluteTolerance Absolute tolerance.
    222             */
    223            public LineSearch(double relativeTolerance,
    224                              double absoluteTolerance) {
    225                optim.setRelativeAccuracy(relativeTolerance);
    226                optim.setAbsoluteAccuracy(absoluteTolerance);
    227            }
    228
    229            /**
    230             * Find the minimum of the function {@code f(p + alpha * d)}.
    231             *
    232             * @param p Starting point.
    233             * @param d Search direction.
    234             * @throws FunctionEvaluationException if function cannot be evaluated at some test point
    235             * @throws OptimizationException if algorithm fails to converge
    236             */
    237            public void search(final double[] p, final double[] d)
    238                throws OptimizationException, FunctionEvaluationException {
    239
    240                // Reset.
    241                optimum = Double.NaN;
    242                valueAtOptimum = Double.NaN;
    243
    244                try {
    245                    final int n = p.length;
    246                    final UnivariateRealFunction f = new UnivariateRealFunction() {
    247                            public double value(double alpha)
    248                                throws FunctionEvaluationException {
    249
    250                                final double[] x = new double[n];
    251                                for (int i = 0; i < n; i++) {
    252                                    x[i] = p[i] + alpha * d[i];
    253                                }
    254                                final double obj;
    255                                obj = computeObjectiveValue(x);
    256                                return obj;
    257                            }
    258                        };
    259
    260                    bracket.search(f, goal, 0, 1);
    261                    optimum = optim.optimize(f, goal,
    262                                             bracket.getLo(),
    263                                             bracket.getHi(),
    264                                             bracket.getMid());
    265                    valueAtOptimum = optim.getFunctionValue();
    266                } catch (MaxIterationsExceededException e) {
    267                    throw new OptimizationException(e);
    268                }
    269            }
    270
    271            /**
    272             * @return the optimum.
    273             */
    274            public double getOptimum() {
    275                return optimum;
    276            }
    277            /**
    278             * @return the value of the function at the optimum.
    279             */
    280            public double getValueAtOptimum() {
    281                return valueAtOptimum;
    282            }
    283        }
    284
    285        /**
    286         * Java 1.5 does not support Arrays.copyOf()
    287         *
    288         * @param source the array to be copied
    289         * @param newLen the length of the copy to be returned
    290         * @return the copied array, truncated or padded as necessary.
    291         */
    292         private double[] copyOf(double[] source, int newLen) {
    293             double[] output = new double[newLen];
    294             System.arraycopy(source, 0, output, 0, Math.min(source.length, newLen));
    295             return output;
    296         }
    297
    298    }