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    }