Java tutorial
/* * Licensed to the Apache Software Foundation (ASF) under one or more * contributor license agreements. See the NOTICE file distributed with * this work for additional information regarding copyright ownership. * The ASF licenses this file to You under the Apache License, Version 2.0 * (the "License"); you may not use this file except in compliance with * the License. You may obtain a copy of the License at * * http://www.apache.org/licenses/LICENSE-2.0 * * Unless required by applicable law or agreed to in writing, software * distributed under the License is distributed on an "AS IS" BASIS, * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. * See the License for the specific language governing permissions and * limitations under the License. */ package org.apache.commons.math.optimization.direct; import org.apache.commons.math.MaxIterationsExceededException; import org.apache.commons.math.analysis.UnivariateRealFunction; import org.apache.commons.math.FunctionEvaluationException; import org.apache.commons.math.optimization.GoalType; import org.apache.commons.math.optimization.OptimizationException; import org.apache.commons.math.optimization.RealPointValuePair; import org.apache.commons.math.optimization.general.AbstractScalarDifferentiableOptimizer; import org.apache.commons.math.optimization.univariate.AbstractUnivariateRealOptimizer; import org.apache.commons.math.optimization.univariate.BracketFinder; import org.apache.commons.math.optimization.univariate.BrentOptimizer; /** * Powell algorithm. * This code is translated and adapted from the Python version of this * algorithm (as implemented in module {@code optimize.py} v0.5 of * <em>SciPy</em>). * * @version $Revision$ $Date$ * @since 2.2 */ public class PowellOptimizer extends AbstractScalarDifferentiableOptimizer { /** * Default relative tolerance for line search ({@value}). */ public static final double DEFAULT_LS_RELATIVE_TOLERANCE = 1e-7; /** * Default absolute tolerance for line search ({@value}). */ public static final double DEFAULT_LS_ABSOLUTE_TOLERANCE = 1e-11; /** * Line search. */ private final LineSearch line; /** * Constructor with default line search tolerances (see the * {@link #PowellOptimizer(double,double) other constructor}). */ public PowellOptimizer() { this(DEFAULT_LS_RELATIVE_TOLERANCE, DEFAULT_LS_ABSOLUTE_TOLERANCE); } /** * Constructor with default absolute line search tolerances (see * the {@link #PowellOptimizer(double,double) other constructor}). * * @param lsRelativeTolerance Relative error tolerance for * the line search algorithm ({@link BrentOptimizer}). */ public PowellOptimizer(double lsRelativeTolerance) { this(lsRelativeTolerance, DEFAULT_LS_ABSOLUTE_TOLERANCE); } /** * @param lsRelativeTolerance Relative error tolerance for * the line search algorithm ({@link BrentOptimizer}). * @param lsAbsoluteTolerance Relative error tolerance for * the line search algorithm ({@link BrentOptimizer}). */ public PowellOptimizer(double lsRelativeTolerance, double lsAbsoluteTolerance) { line = new LineSearch(lsRelativeTolerance, lsAbsoluteTolerance); } /** {@inheritDoc} */ @Override protected RealPointValuePair doOptimize() throws FunctionEvaluationException, OptimizationException { final double[] guess = point.clone(); final int n = guess.length; final double[][] direc = new double[n][n]; for (int i = 0; i < n; i++) { direc[i][i] = 1; } double[] x = guess; double fVal = computeObjectiveValue(x); double[] x1 = x.clone(); while (true) { incrementIterationsCounter(); double fX = fVal; double fX2 = 0; double delta = 0; int bigInd = 0; double alphaMin = 0; for (int i = 0; i < n; i++) { final double[] d = /* Arrays.*/ copyOf(direc[i], n); // Java 1.5 does not support Arrays.copyOf() fX2 = fVal; line.search(x, d); fVal = line.getValueAtOptimum(); alphaMin = line.getOptimum(); final double[][] result = newPointAndDirection(x, d, alphaMin); x = result[0]; if ((fX2 - fVal) > delta) { delta = fX2 - fVal; bigInd = i; } } final RealPointValuePair previous = new RealPointValuePair(x1, fX); final RealPointValuePair current = new RealPointValuePair(x, fVal); if (getConvergenceChecker().converged(getIterations(), previous, current)) { if (goal == GoalType.MINIMIZE) { return (fVal < fX) ? current : previous; } else { return (fVal > fX) ? current : previous; } } final double[] d = new double[n]; final double[] x2 = new double[n]; for (int i = 0; i < n; i++) { d[i] = x[i] - x1[i]; x2[i] = 2 * x[i] - x1[i]; } x1 = x.clone(); fX2 = computeObjectiveValue(x2); if (fX > fX2) { double t = 2 * (fX + fX2 - 2 * fVal); double temp = fX - fVal - delta; t *= temp * temp; temp = fX - fX2; t -= delta * temp * temp; if (t < 0.0) { line.search(x, d); fVal = line.getValueAtOptimum(); alphaMin = line.getOptimum(); final double[][] result = newPointAndDirection(x, d, alphaMin); x = result[0]; final int lastInd = n - 1; direc[bigInd] = direc[lastInd]; direc[lastInd] = result[1]; } } } } /** * Compute a new point (in the original space) and a new direction * vector, resulting from the line search. * The parameters {@code p} and {@code d} will be changed in-place. * * @param p Point used in the line search. * @param d Direction used in the line search. * @param optimum Optimum found by the line search. * @return a 2-element array containing the new point (at index 0) and * the new direction (at index 1). */ private double[][] newPointAndDirection(double[] p, double[] d, double optimum) { final int n = p.length; final double[][] result = new double[2][n]; final double[] nP = result[0]; final double[] nD = result[1]; for (int i = 0; i < n; i++) { nD[i] = d[i] * optimum; nP[i] = p[i] + nD[i]; } return result; } /** * Class for finding the minimum of the objective function along a given * direction. */ private class LineSearch { /** * Optimizer. */ private final AbstractUnivariateRealOptimizer optim = new BrentOptimizer(); /** * Automatic bracketing. */ private final BracketFinder bracket = new BracketFinder(); /** * Value of the optimum. */ private double optimum = Double.NaN; /** * Value of the objective function at the optimum. */ private double valueAtOptimum = Double.NaN; /** * @param relativeTolerance Relative tolerance. * @param absoluteTolerance Absolute tolerance. */ public LineSearch(double relativeTolerance, double absoluteTolerance) { optim.setRelativeAccuracy(relativeTolerance); optim.setAbsoluteAccuracy(absoluteTolerance); } /** * Find the minimum of the function {@code f(p + alpha * d)}. * * @param p Starting point. * @param d Search direction. * @throws FunctionEvaluationException if function cannot be evaluated at some test point * @throws OptimizationException if algorithm fails to converge */ public void search(final double[] p, final double[] d) throws OptimizationException, FunctionEvaluationException { // Reset. optimum = Double.NaN; valueAtOptimum = Double.NaN; try { final int n = p.length; final UnivariateRealFunction f = new UnivariateRealFunction() { public double value(double alpha) throws FunctionEvaluationException { final double[] x = new double[n]; for (int i = 0; i < n; i++) { x[i] = p[i] + alpha * d[i]; } final double obj; obj = computeObjectiveValue(x); return obj; } }; bracket.search(f, goal, 0, 1); optimum = optim.optimize(f, goal, bracket.getLo(), bracket.getHi(), bracket.getMid()); valueAtOptimum = optim.getFunctionValue(); } catch (MaxIterationsExceededException e) { throw new OptimizationException(e); } } /** * @return the optimum. */ public double getOptimum() { return optimum; } /** * @return the value of the function at the optimum. */ public double getValueAtOptimum() { return valueAtOptimum; } } /** * Java 1.5 does not support Arrays.copyOf() * * @param source the array to be copied * @param newLen the length of the copy to be returned * @return the copied array, truncated or padded as necessary. */ private double[] copyOf(double[] source, int newLen) { double[] output = new double[newLen]; System.arraycopy(source, 0, output, 0, Math.min(source.length, newLen)); return output; } }