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 */
017package org.apache.commons.math4.legacy.optim.nonlinear.scalar.noderiv;
018
019import java.util.Arrays;
020import org.apache.commons.math4.legacy.analysis.MultivariateFunction;
021import org.apache.commons.math4.legacy.exception.MathUnsupportedOperationException;
022import org.apache.commons.math4.legacy.exception.NotStrictlyPositiveException;
023import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException;
024import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
025import org.apache.commons.math4.legacy.optim.ConvergenceChecker;
026import org.apache.commons.math4.legacy.optim.PointValuePair;
027import org.apache.commons.math4.legacy.optim.nonlinear.scalar.GoalType;
028import org.apache.commons.math4.legacy.optim.nonlinear.scalar.MultivariateOptimizer;
029import org.apache.commons.math4.legacy.optim.univariate.UnivariatePointValuePair;
030import org.apache.commons.math4.core.jdkmath.JdkMath;
031
032/**
033 * Powell's algorithm.
034 * This code is translated and adapted from the Python version of this
035 * algorithm (as implemented in module {@code optimize.py} v0.5 of
036 * <em>SciPy</em>).
037 * <br>
038 * The default stopping criterion is based on the differences of the
039 * function value between two successive iterations. It is however possible
040 * to define a custom convergence checker that might terminate the algorithm
041 * earlier.
042 * <br>
043 * Constraints are not supported: the call to
044 * {@link #optimize(org.apache.commons.math4.legacy.optim.OptimizationData...)} will throw
045 * {@link MathUnsupportedOperationException} if bounds are passed to it.
046 * In order to impose simple constraints, the objective function must be
047 * wrapped in an adapter like
048 * {@link org.apache.commons.math4.legacy.optim.nonlinear.scalar.MultivariateFunctionMappingAdapter
049 * MultivariateFunctionMappingAdapter} or
050 * {@link org.apache.commons.math4.legacy.optim.nonlinear.scalar.MultivariateFunctionPenaltyAdapter
051 * MultivariateFunctionPenaltyAdapter}.
052 *
053 * @since 2.2
054 */
055public class PowellOptimizer
056    extends MultivariateOptimizer {
057    /**
058     * Minimum relative tolerance.
059     */
060    private static final double MIN_RELATIVE_TOLERANCE = 2 * JdkMath.ulp(1d);
061    /** Relative threshold. */
062    private final double relativeThreshold;
063    /** Absolute threshold. */
064    private final double absoluteThreshold;
065    /** Relative threshold. */
066    private final double lineSearchRelativeThreshold;
067    /** Absolute threshold. */
068    private final double lineSearchAbsoluteThreshold;
069
070    /**
071     * This constructor allows to specify a user-defined convergence checker,
072     * in addition to the parameters that control the default convergence
073     * checking procedure.
074     * <br>
075     * The internal line search tolerances are set to the square-root of their
076     * corresponding value in the multivariate optimizer.
077     *
078     * @param rel Relative threshold.
079     * @param abs Absolute threshold.
080     * @param checker Convergence checker.
081     * @throws NotStrictlyPositiveException if {@code abs <= 0}.
082     * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}.
083     */
084    public PowellOptimizer(double rel,
085                           double abs,
086                           ConvergenceChecker<PointValuePair> checker) {
087        this(rel, abs, JdkMath.sqrt(rel), JdkMath.sqrt(abs), checker);
088    }
089
090    /**
091     * This constructor allows to specify a user-defined convergence checker,
092     * in addition to the parameters that control the default convergence
093     * checking procedure and the line search tolerances.
094     *
095     * @param rel Relative threshold for this optimizer.
096     * @param abs Absolute threshold for this optimizer.
097     * @param lineRel Relative threshold for the internal line search optimizer.
098     * @param lineAbs Absolute threshold for the internal line search optimizer.
099     * @param checker Convergence checker.
100     * @throws NotStrictlyPositiveException if {@code abs <= 0}.
101     * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}.
102     */
103    public PowellOptimizer(double rel,
104                           double abs,
105                           double lineRel,
106                           double lineAbs,
107                           ConvergenceChecker<PointValuePair> checker) {
108        super(checker);
109
110        if (rel < MIN_RELATIVE_TOLERANCE) {
111            throw new NumberIsTooSmallException(rel, MIN_RELATIVE_TOLERANCE, true);
112        }
113        if (abs <= 0) {
114            throw new NotStrictlyPositiveException(abs);
115        }
116
117        relativeThreshold = rel;
118        absoluteThreshold = abs;
119        lineSearchRelativeThreshold = lineRel;
120        lineSearchAbsoluteThreshold = lineAbs;
121    }
122
123    /**
124     * The parameters control the default convergence checking procedure.
125     * <br>
126     * The internal line search tolerances are set to the square-root of their
127     * corresponding value in the multivariate optimizer.
128     *
129     * @param rel Relative threshold.
130     * @param abs Absolute threshold.
131     * @throws NotStrictlyPositiveException if {@code abs <= 0}.
132     * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}.
133     */
134    public PowellOptimizer(double rel,
135                           double abs) {
136        this(rel, abs, null);
137    }
138
139    /**
140     * Builds an instance with the default convergence checking procedure.
141     *
142     * @param rel Relative threshold.
143     * @param abs Absolute threshold.
144     * @param lineRel Relative threshold for the internal line search optimizer.
145     * @param lineAbs Absolute threshold for the internal line search optimizer.
146     * @throws NotStrictlyPositiveException if {@code abs <= 0}.
147     * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}.
148     */
149    public PowellOptimizer(double rel,
150                           double abs,
151                           double lineRel,
152                           double lineAbs) {
153        this(rel, abs, lineRel, lineAbs, null);
154    }
155
156    /** {@inheritDoc} */
157    @Override
158    protected PointValuePair doOptimize() {
159        checkParameters();
160
161        // Line search optimizer.
162        createLineSearch();
163
164        final GoalType goal = getGoalType();
165        final double[] guess = getStartPoint();
166        final MultivariateFunction func = getObjectiveFunction();
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 = func.value(x);
179        double[] x1 = x.clone();
180        while (true) {
181            incrementIterationCount();
182
183            double fX = fVal;
184            double fX2 = 0;
185            double delta = 0;
186            int bigInd = 0;
187            double alphaMin = 0;
188
189            for (int i = 0; i < n; i++) {
190                final double[] d = Arrays.copyOf(direc[i], direc[i].length);
191
192                fX2 = fVal;
193
194                final UnivariatePointValuePair optimum = lineSearch(x, d);
195                fVal = optimum.getValue();
196                alphaMin = optimum.getPoint();
197                final double[][] result = newPointAndDirection(x, d, alphaMin);
198                x = result[0];
199
200                if ((fX2 - fVal) > delta) {
201                    delta = fX2 - fVal;
202                    bigInd = i;
203                }
204            }
205
206            // Default convergence check.
207            boolean stop = 2 * (fX - fVal) <=
208                (relativeThreshold * (JdkMath.abs(fX) + JdkMath.abs(fVal)) +
209                 absoluteThreshold);
210
211            final PointValuePair previous = new PointValuePair(x1, fX);
212            final PointValuePair current = new PointValuePair(x, fVal);
213            if (!stop && checker != null) { // User-defined stopping criteria.
214                stop = checker.converged(getIterations(), previous, current);
215            }
216            if (stop) {
217                if (goal == GoalType.MINIMIZE) {
218                    return (fVal < fX) ? current : previous;
219                } else {
220                    return (fVal > fX) ? current : previous;
221                }
222            }
223
224            final double[] d = new double[n];
225            final double[] x2 = new double[n];
226            for (int i = 0; i < n; i++) {
227                d[i] = x[i] - x1[i];
228                x2[i] = 2 * x[i] - x1[i];
229            }
230
231            x1 = x.clone();
232            fX2 = func.value(x2);
233
234            if (fX > fX2) {
235                double t = 2 * (fX + fX2 - 2 * fVal);
236                double temp = fX - fVal - delta;
237                t *= temp * temp;
238                temp = fX - fX2;
239                t -= delta * temp * temp;
240
241                if (t < 0.0) {
242                    final UnivariatePointValuePair optimum = lineSearch(x, d);
243                    fVal = optimum.getValue();
244                    alphaMin = optimum.getPoint();
245                    final double[][] result = newPointAndDirection(x, d, alphaMin);
246                    x = result[0];
247
248                    final int lastInd = n - 1;
249                    direc[bigInd] = direc[lastInd];
250                    direc[lastInd] = result[1];
251                }
252            }
253        }
254    }
255
256    /**
257     * Compute a new point (in the original space) and a new direction
258     * vector, resulting from the line search.
259     *
260     * @param p Point used in the line search.
261     * @param d Direction used in the line search.
262     * @param optimum Optimum found by the line search.
263     * @return a 2-element array containing the new point (at index 0) and
264     * the new direction (at index 1).
265     */
266    private double[][] newPointAndDirection(double[] p,
267                                            double[] d,
268                                            double optimum) {
269        final int n = p.length;
270        final double[] nP = new double[n];
271        final double[] nD = new double[n];
272        for (int i = 0; i < n; i++) {
273            nD[i] = d[i] * optimum;
274            nP[i] = p[i] + nD[i];
275        }
276
277        final double[][] result = new double[2][];
278        result[0] = nP;
279        result[1] = nD;
280
281        return result;
282    }
283
284    /**
285     * @throws MathUnsupportedOperationException if bounds were passed to the
286     * {@link #optimize(org.apache.commons.math4.legacy.optim.OptimizationData[]) optimize} method.
287     */
288    private void checkParameters() {
289        if (getLowerBound() != null ||
290            getUpperBound() != null) {
291            throw new MathUnsupportedOperationException(LocalizedFormats.CONSTRAINT);
292        }
293    }
294}