1/*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements.  See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License.  You may obtain a copy of the License at
8 *
9 *      http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17
18package org.apache.commons.math.optimization.direct;
19
20import org.apache.commons.math.MaxIterationsExceededException;
21import org.apache.commons.math.analysis.UnivariateRealFunction;
22import org.apache.commons.math.FunctionEvaluationException;
23import org.apache.commons.math.optimization.GoalType;
24import org.apache.commons.math.optimization.OptimizationException;
25import org.apache.commons.math.optimization.RealPointValuePair;
26import org.apache.commons.math.optimization.general.AbstractScalarDifferentiableOptimizer;
27import org.apache.commons.math.optimization.univariate.AbstractUnivariateRealOptimizer;
28import org.apache.commons.math.optimization.univariate.BracketFinder;
29import org.apache.commons.math.optimization.univariate.BrentOptimizer;
30
31/**
32 * Powell algorithm.
33 * This code is translated and adapted from the Python version of this
34 * algorithm (as implemented in module {@code optimize.py} v0.5 of
35 * <em>SciPy</em>).
36 *
37 * @version $Revision$ $Date$
38 * @since 2.2
39 */
40public class PowellOptimizer
41    extends AbstractScalarDifferentiableOptimizer {
42    /**
43     * Default relative tolerance for line search ({@value}).
44     */
45    public static final double DEFAULT_LS_RELATIVE_TOLERANCE = 1e-7;
46    /**
47     * Default absolute tolerance for line search ({@value}).
48     */
49    public static final double DEFAULT_LS_ABSOLUTE_TOLERANCE = 1e-11;
50    /**
51     * Line search.
52     */
53    private final LineSearch line;
54
55    /**
56     * Constructor with default line search tolerances (see the
57     * {@link #PowellOptimizer(double,double) other constructor}).
58     */
59    public PowellOptimizer() {
60        this(DEFAULT_LS_RELATIVE_TOLERANCE,
61             DEFAULT_LS_ABSOLUTE_TOLERANCE);
62    }
63
64    /**
65     * Constructor with default absolute line search tolerances (see
66     * the {@link #PowellOptimizer(double,double) other constructor}).
67     *
68     * @param lsRelativeTolerance Relative error tolerance for
69     * the line search algorithm ({@link BrentOptimizer}).
70     */
71    public PowellOptimizer(double lsRelativeTolerance) {
72        this(lsRelativeTolerance,
73             DEFAULT_LS_ABSOLUTE_TOLERANCE);
74    }
75
76    /**
77     * @param lsRelativeTolerance Relative error tolerance for
78     * the line search algorithm ({@link BrentOptimizer}).
79     * @param lsAbsoluteTolerance Relative error tolerance for
80     * the line search algorithm ({@link BrentOptimizer}).
81     */
82    public PowellOptimizer(double lsRelativeTolerance,
83                           double lsAbsoluteTolerance) {
84        line = new LineSearch(lsRelativeTolerance,
85                              lsAbsoluteTolerance);
86    }
87
88    /** {@inheritDoc} */
89    @Override
90    protected RealPointValuePair doOptimize()
91        throws FunctionEvaluationException,
92               OptimizationException {
93        final double[] guess = point.clone();
94        final int n = guess.length;
95
96        final double[][] direc = new double[n][n];
97        for (int i = 0; i < n; i++) {
98            direc[i][i] = 1;
99        }
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}
299