    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.math3.optimization.direct;
    019
    020    import org.apache.commons.math3.util.FastMath;
    021    import org.apache.commons.math3.util.MathArrays;
    022    import org.apache.commons.math3.analysis.UnivariateFunction;
    023    import org.apache.commons.math3.analysis.MultivariateFunction;
    024    import org.apache.commons.math3.exception.NumberIsTooSmallException;
    025    import org.apache.commons.math3.exception.NotStrictlyPositiveException;
    026    import org.apache.commons.math3.optimization.GoalType;
    027    import org.apache.commons.math3.optimization.PointValuePair;
    028    import org.apache.commons.math3.optimization.ConvergenceChecker;
    029    import org.apache.commons.math3.optimization.MultivariateOptimizer;
    030    import org.apache.commons.math3.optimization.univariate.BracketFinder;
    031    import org.apache.commons.math3.optimization.univariate.BrentOptimizer;
    032    import org.apache.commons.math3.optimization.univariate.UnivariatePointValuePair;
    033    import org.apache.commons.math3.optimization.univariate.SimpleUnivariateValueChecker;
    034
    035    /**
    036     * Powell algorithm.
    037     * This code is translated and adapted from the Python version of this
    038     * algorithm (as implemented in module {@code optimize.py} v0.5 of
    039     * <em>SciPy</em>).
    040     * <br/>
    041     * The default stopping criterion is based on the differences of the
    042     * function value between two successive iterations. It is however possible
    043     * to define a custom convergence checker that might terminate the algorithm
    044     * earlier.
    045     * <br/>
    046     * The internal line search optimizer is a {@link BrentOptimizer} with a
    047     * convergence checker set to {@link SimpleUnivariateValueChecker}.
    048     *
    049     * @version $Id: PowellOptimizer.java 1422313 2012-12-15 18:53:41Z psteitz $
    050     * @deprecated As of 3.1 (to be removed in 4.0).
    051     * @since 2.2
    052     */
    053    @Deprecated
    054    public class PowellOptimizer
    055        extends BaseAbstractMultivariateOptimizer<MultivariateFunction>
    056        implements MultivariateOptimizer {
    057        /**
    058         * Minimum relative tolerance.
    059         */
    060        private static final double MIN_RELATIVE_TOLERANCE = 2 * FastMath.ulp(1d);
    061        /**
    062         * Relative threshold.
    063         */
    064        private final double relativeThreshold;
    065        /**
    066         * Absolute threshold.
    067         */
    068        private final double absoluteThreshold;
    069        /**
    070         * Line search.
    071         */
    072        private final LineSearch line;
    073
    074        /**
    075         * This constructor allows to specify a user-defined convergence checker,
    076         * in addition to the parameters that control the default convergence
    077         * checking procedure.
    078         * <br/>
    079         * The internal line search tolerances are set to the square-root of their
    080         * corresponding value in the multivariate optimizer.
    081         *
    082         * @param rel Relative threshold.
    083         * @param abs Absolute threshold.
    084         * @param checker Convergence checker.
    085         * @throws NotStrictlyPositiveException if {@code abs <= 0}.
    086         * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}.
    087         */
    088        public PowellOptimizer(double rel,
    089                               double abs,
    090                               ConvergenceChecker<PointValuePair> checker) {
    091            this(rel, abs, FastMath.sqrt(rel), FastMath.sqrt(abs), checker);
    092        }
    093
    094        /**
    095         * This constructor allows to specify a user-defined convergence checker,
    096         * in addition to the parameters that control the default convergence
    097         * checking procedure and the line search tolerances.
    098         *
    099         * @param rel Relative threshold for this optimizer.
    100         * @param abs Absolute threshold for this optimizer.
    101         * @param lineRel Relative threshold for the internal line search optimizer.
    102         * @param lineAbs Absolute threshold for the internal line search optimizer.
    103         * @param checker Convergence checker.
    104         * @throws NotStrictlyPositiveException if {@code abs <= 0}.
    105         * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}.
    106         */
    107        public PowellOptimizer(double rel,
    108                               double abs,
    109                               double lineRel,
    110                               double lineAbs,
    111                               ConvergenceChecker<PointValuePair> checker) {
    112            super(checker);
    113
    114            if (rel < MIN_RELATIVE_TOLERANCE) {
    115                throw new NumberIsTooSmallException(rel, MIN_RELATIVE_TOLERANCE, true);
    116            }
    117            if (abs <= 0) {
    118                throw new NotStrictlyPositiveException(abs);
    119            }
    120            relativeThreshold = rel;
    121            absoluteThreshold = abs;
    122
    123            // Create the line search optimizer.
    124            line = new LineSearch(lineRel,
    125                                  lineAbs);
    126        }
    127
    128        /**
    129         * The parameters control the default convergence checking procedure.
    130         * <br/>
    131         * The internal line search tolerances are set to the square-root of their
    132         * corresponding value in the multivariate optimizer.
    133         *
    134         * @param rel Relative threshold.
    135         * @param abs Absolute threshold.
    136         * @throws NotStrictlyPositiveException if {@code abs <= 0}.
    137         * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}.
    138         */
    139        public PowellOptimizer(double rel,
    140                               double abs) {
    141            this(rel, abs, null);
    142        }
    143
    144        /**
    145         * Builds an instance with the default convergence checking procedure.
    146         *
    147         * @param rel Relative threshold.
    148         * @param abs Absolute threshold.
    149         * @param lineRel Relative threshold for the internal line search optimizer.
    150         * @param lineAbs Absolute threshold for the internal line search optimizer.
    151         * @throws NotStrictlyPositiveException if {@code abs <= 0}.
    152         * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}.
    153         * @since 3.1
    154         */
    155        public PowellOptimizer(double rel,
    156                               double abs,
    157                               double lineRel,
    158                               double lineAbs) {
    159            this(rel, abs, lineRel, lineAbs, null);
    160        }
    161
    162        /** {@inheritDoc} */
    163        @Override
    164        protected PointValuePair doOptimize() {
    165            final GoalType goal = getGoalType();
    166            final double[] guess = getStartPoint();
    167            final int n = guess.length;
    168
    169            final double[][] direc = new double[n][n];
    170            for (int i = 0; i < n; i++) {
    171                direc[i][i] = 1;
    172            }
    173
    174            final ConvergenceChecker<PointValuePair> checker
    175                = getConvergenceChecker();
    176
    177            double[] x = guess;
    178            double fVal = computeObjectiveValue(x);
    179            double[] x1 = x.clone();
    180            int iter = 0;
    181            while (true) {
    182                ++iter;
    183
    184                double fX = fVal;
    185                double fX2 = 0;
    186                double delta = 0;
    187                int bigInd = 0;
    188                double alphaMin = 0;
    189
    190                for (int i = 0; i < n; i++) {
    191                    final double[] d = MathArrays.copyOf(direc[i]);
    192
    193                    fX2 = fVal;
    194
    195                    final UnivariatePointValuePair optimum = line.search(x, d);
    196                    fVal = optimum.getValue();
    197                    alphaMin = optimum.getPoint();
    198                    final double[][] result = newPointAndDirection(x, d, alphaMin);
    199                    x = result[0];
    200
    201                    if ((fX2 - fVal) > delta) {
    202                        delta = fX2 - fVal;
    203                        bigInd = i;
    204                    }
    205                }
    206
    207                // Default convergence check.
    208                boolean stop = 2 * (fX - fVal) <=
    209                    (relativeThreshold * (FastMath.abs(fX) + FastMath.abs(fVal)) +
    210                     absoluteThreshold);
    211
    212                final PointValuePair previous = new PointValuePair(x1, fX);
    213                final PointValuePair current = new PointValuePair(x, fVal);
    214                if (!stop) { // User-defined stopping criteria.
    215                    if (checker != null) {
    216                        stop = checker.converged(iter, previous, current);
    217                    }
    218                }
    219                if (stop) {
    220                    if (goal == GoalType.MINIMIZE) {
    221                        return (fVal < fX) ? current : previous;
    222                    } else {
    223                        return (fVal > fX) ? current : previous;
    224                    }
    225                }
    226
    227                final double[] d = new double[n];
    228                final double[] x2 = new double[n];
    229                for (int i = 0; i < n; i++) {
    230                    d[i] = x[i] - x1[i];
    231                    x2[i] = 2 * x[i] - x1[i];
    232                }
    233
    234                x1 = x.clone();
    235                fX2 = computeObjectiveValue(x2);
    236
    237                if (fX > fX2) {
    238                    double t = 2 * (fX + fX2 - 2 * fVal);
    239                    double temp = fX - fVal - delta;
    240                    t *= temp * temp;
    241                    temp = fX - fX2;
    242                    t -= delta * temp * temp;
    243
    244                    if (t < 0.0) {
    245                        final UnivariatePointValuePair optimum = line.search(x, d);
    246                        fVal = optimum.getValue();
    247                        alphaMin = optimum.getPoint();
    248                        final double[][] result = newPointAndDirection(x, d, alphaMin);
    249                        x = result[0];
    250
    251                        final int lastInd = n - 1;
    252                        direc[bigInd] = direc[lastInd];
    253                        direc[lastInd] = result[1];
    254                    }
    255                }
    256            }
    257        }
    258
    259        /**
    260         * Compute a new point (in the original space) and a new direction
    261         * vector, resulting from the line search.
    262         *
    263         * @param p Point used in the line search.
    264         * @param d Direction used in the line search.
    265         * @param optimum Optimum found by the line search.
    266         * @return a 2-element array containing the new point (at index 0) and
    267         * the new direction (at index 1).
    268         */
    269        private double[][] newPointAndDirection(double[] p,
    270                                                double[] d,
    271                                                double optimum) {
    272            final int n = p.length;
    273            final double[] nP = new double[n];
    274            final double[] nD = new double[n];
    275            for (int i = 0; i < n; i++) {
    276                nD[i] = d[i] * optimum;
    277                nP[i] = p[i] + nD[i];
    278            }
    279
    280            final double[][] result = new double[2][];
    281            result[0] = nP;
    282            result[1] = nD;
    283
    284            return result;
    285        }
    286
    287        /**
    288         * Class for finding the minimum of the objective function along a given
    289         * direction.
    290         */
    291        private class LineSearch extends BrentOptimizer {
    292            /**
    293             * Value that will pass the precondition check for {@link BrentOptimizer}
    294             * but will not pass the convergence check, so that the custom checker
    295             * will always decide when to stop the line search.
    296             */
    297            private static final double REL_TOL_UNUSED = 1e-15;
    298            /**
    299             * Value that will pass the precondition check for {@link BrentOptimizer}
    300             * but will not pass the convergence check, so that the custom checker
    301             * will always decide when to stop the line search.
    302             */
    303            private static final double ABS_TOL_UNUSED = Double.MIN_VALUE;
    304            /**
    305             * Automatic bracketing.
    306             */
    307            private final BracketFinder bracket = new BracketFinder();
    308
    309            /**
    310             * The "BrentOptimizer" default stopping criterion uses the tolerances
    311             * to check the domain (point) values, not the function values.
    312             * We thus create a custom checker to use function values.
    313             *
    314             * @param rel Relative threshold.
    315             * @param abs Absolute threshold.
    316             */
    317            LineSearch(double rel,
    318                       double abs) {
    319                super(REL_TOL_UNUSED,
    320                      ABS_TOL_UNUSED,
    321                      new SimpleUnivariateValueChecker(rel, abs));
    322            }
    323
    324            /**
    325             * Find the minimum of the function {@code f(p + alpha * d)}.
    326             *
    327             * @param p Starting point.
    328             * @param d Search direction.
    329             * @return the optimum.
    330             * @throws org.apache.commons.math3.exception.TooManyEvaluationsException
    331             * if the number of evaluations is exceeded.
    332             */
    333            public UnivariatePointValuePair search(final double[] p, final double[] d) {
    334                final int n = p.length;
    335                final UnivariateFunction f = new UnivariateFunction() {
    336                        public double value(double alpha) {
    337                            final double[] x = new double[n];
    338                            for (int i = 0; i < n; i++) {
    339                                x[i] = p[i] + alpha * d[i];
    340                            }
    341                            final double obj = PowellOptimizer.this.computeObjectiveValue(x);
    342                            return obj;
    343                        }
    344                    };
    345
    346                final GoalType goal = PowellOptimizer.this.getGoalType();
    347                bracket.search(f, goal, 0, 1);
    348                // Passing "MAX_VALUE" as a dummy value because it is the enclosing
    349                // class that counts the number of evaluations (and will eventually
    350                // generate the exception).
    351                return optimize(Integer.MAX_VALUE, f, goal,
    352                                bracket.getLo(), bracket.getHi(), bracket.getMid());
    353            }
    354        }
    355    }