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
018package org.apache.commons.math3.optimization.direct;
019
020import org.apache.commons.math3.util.FastMath;
021import org.apache.commons.math3.util.MathArrays;
022import org.apache.commons.math3.analysis.UnivariateFunction;
023import org.apache.commons.math3.analysis.MultivariateFunction;
024import org.apache.commons.math3.exception.NumberIsTooSmallException;
025import org.apache.commons.math3.exception.NotStrictlyPositiveException;
026import org.apache.commons.math3.optimization.GoalType;
027import org.apache.commons.math3.optimization.PointValuePair;
028import org.apache.commons.math3.optimization.ConvergenceChecker;
029import org.apache.commons.math3.optimization.MultivariateOptimizer;
030import org.apache.commons.math3.optimization.univariate.BracketFinder;
031import org.apache.commons.math3.optimization.univariate.BrentOptimizer;
032import org.apache.commons.math3.optimization.univariate.UnivariatePointValuePair;
033import 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 1462503 2013-03-29 15:48:27Z luc $
050 * @deprecated As of 3.1 (to be removed in 4.0).
051 * @since 2.2
052 */
053@Deprecated
054public 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 && checker != null) {
215                stop = checker.converged(iter, previous, current);
216            }
217            if (stop) {
218                if (goal == GoalType.MINIMIZE) {
219                    return (fVal < fX) ? current : previous;
220                } else {
221                    return (fVal > fX) ? current : previous;
222                }
223            }
224
225            final double[] d = new double[n];
226            final double[] x2 = new double[n];
227            for (int i = 0; i < n; i++) {
228                d[i] = x[i] - x1[i];
229                x2[i] = 2 * x[i] - x1[i];
230            }
231
232            x1 = x.clone();
233            fX2 = computeObjectiveValue(x2);
234
235            if (fX > fX2) {
236                double t = 2 * (fX + fX2 - 2 * fVal);
237                double temp = fX - fVal - delta;
238                t *= temp * temp;
239                temp = fX - fX2;
240                t -= delta * temp * temp;
241
242                if (t < 0.0) {
243                    final UnivariatePointValuePair optimum = line.search(x, d);
244                    fVal = optimum.getValue();
245                    alphaMin = optimum.getPoint();
246                    final double[][] result = newPointAndDirection(x, d, alphaMin);
247                    x = result[0];
248
249                    final int lastInd = n - 1;
250                    direc[bigInd] = direc[lastInd];
251                    direc[lastInd] = result[1];
252                }
253            }
254        }
255    }
256
257    /**
258     * Compute a new point (in the original space) and a new direction
259     * vector, resulting from the line search.
260     *
261     * @param p Point used in the line search.
262     * @param d Direction used in the line search.
263     * @param optimum Optimum found by the line search.
264     * @return a 2-element array containing the new point (at index 0) and
265     * the new direction (at index 1).
266     */
267    private double[][] newPointAndDirection(double[] p,
268                                            double[] d,
269                                            double optimum) {
270        final int n = p.length;
271        final double[] nP = new double[n];
272        final double[] nD = new double[n];
273        for (int i = 0; i < n; i++) {
274            nD[i] = d[i] * optimum;
275            nP[i] = p[i] + nD[i];
276        }
277
278        final double[][] result = new double[2][];
279        result[0] = nP;
280        result[1] = nD;
281
282        return result;
283    }
284
285    /**
286     * Class for finding the minimum of the objective function along a given
287     * direction.
288     */
289    private class LineSearch extends BrentOptimizer {
290        /**
291         * Value that will pass the precondition check for {@link BrentOptimizer}
292         * but will not pass the convergence check, so that the custom checker
293         * will always decide when to stop the line search.
294         */
295        private static final double REL_TOL_UNUSED = 1e-15;
296        /**
297         * Value that will pass the precondition check for {@link BrentOptimizer}
298         * but will not pass the convergence check, so that the custom checker
299         * will always decide when to stop the line search.
300         */
301        private static final double ABS_TOL_UNUSED = Double.MIN_VALUE;
302        /**
303         * Automatic bracketing.
304         */
305        private final BracketFinder bracket = new BracketFinder();
306
307        /**
308         * The "BrentOptimizer" default stopping criterion uses the tolerances
309         * to check the domain (point) values, not the function values.
310         * We thus create a custom checker to use function values.
311         *
312         * @param rel Relative threshold.
313         * @param abs Absolute threshold.
314         */
315        LineSearch(double rel,
316                   double abs) {
317            super(REL_TOL_UNUSED,
318                  ABS_TOL_UNUSED,
319                  new SimpleUnivariateValueChecker(rel, abs));
320        }
321
322        /**
323         * Find the minimum of the function {@code f(p + alpha * d)}.
324         *
325         * @param p Starting point.
326         * @param d Search direction.
327         * @return the optimum.
328         * @throws org.apache.commons.math3.exception.TooManyEvaluationsException
329         * if the number of evaluations is exceeded.
330         */
331        public UnivariatePointValuePair search(final double[] p, final double[] d) {
332            final int n = p.length;
333            final UnivariateFunction f = new UnivariateFunction() {
334                    public double value(double alpha) {
335                        final double[] x = new double[n];
336                        for (int i = 0; i < n; i++) {
337                            x[i] = p[i] + alpha * d[i];
338                        }
339                        final double obj = PowellOptimizer.this.computeObjectiveValue(x);
340                        return obj;
341                    }
342                };
343
344            final GoalType goal = PowellOptimizer.this.getGoalType();
345            bracket.search(f, goal, 0, 1);
346            // Passing "MAX_VALUE" as a dummy value because it is the enclosing
347            // class that counts the number of evaluations (and will eventually
348            // generate the exception).
349            return optimize(Integer.MAX_VALUE, f, goal,
350                            bracket.getLo(), bracket.getHi(), bracket.getMid());
351        }
352    }
353}