PowellMultiOptimizer

/*
 * Copyright (C) 2012 Alberto Irurueta Carro (alberto@irurueta.com)
 *
 * Licensed 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 com.irurueta.numerical.optimization;

import com.irurueta.algebra.AlgebraException;
import com.irurueta.algebra.Matrix;
import com.irurueta.algebra.WrongSizeException;
import com.irurueta.numerical.EvaluationException;
import com.irurueta.numerical.LockedException;
import com.irurueta.numerical.MultiDimensionFunctionEvaluatorListener;
import com.irurueta.numerical.NotAvailableException;
import com.irurueta.numerical.NotReadyException;

/**
 * This class searches for a multi dimension function local minimum.
 * The local minimum is searched by starting the algorithm at a start point
 * and a given direction which should be close and point to the local minimum to
 * be found to achieve the best accuracy with the lowest number of iterations.
 * The implementation of this class is based on Numerical Recipes 3rd ed.
 * Section 10.7 page 509.
 */
public class PowellMultiOptimizer extends LineMultiOptimizer {
    /**
     * Constant defining default tolerance or accuracy to be achieved on the
     * minimum being estimated by this class.
     */
    public static final double DEFAULT_TOLERANCE = 3e-8;

    /**
     * Minimum allowed tolerance value.
     */
    public static final double MIN_TOLERANCE = 0.0;

    /**
     * Maximum allowed iterations.
     */
    public static final int ITMAX = 200;

    /**
     * A small number.
     */
    public static final double TINY = 1e-25;

    /**
     * The fractional tolerance in the function value such that failure to
     * decrease by more than this amount on one iteration signals doneness.
     */
    private double tolerance;

    /**
     * Member contains number of iterations that were needed to estimate a
     * minimum.
     */
    private int iter;

    /**
     * Value of the function at the minimum.
     */
    private double fret;

    /**
     * Set of directions.
     */
    private Matrix ximat;

    /**
     * Empty constructor.
     */
    public PowellMultiOptimizer() {
        super();
        tolerance = DEFAULT_TOLERANCE;
        iter = 0;
        fret = 0.0;
        ximat = null;
    }

    /**
     * Constructor.
     *
     * @param listener  Listener to evaluate a multidimensional function.
     * @param tolerance Tolerance or accuracy to be expected on estimated local
     *                  minimum.
     * @throws IllegalArgumentException Raised if tolerance is negative.
     */
    public PowellMultiOptimizer(final MultiDimensionFunctionEvaluatorListener listener, final double tolerance) {
        super(listener);
        internalSetTolerance(tolerance);
        iter = 0;
        fret = 0.0;
        ximat = null;
    }

    /**
     * Constructor.
     *
     * @param listener   Listener to evaluate a multidimensional function.
     * @param point      Start point where algorithm will be started. Start point
     *                   should be close to the local minimum to be found. Provided array must
     *                   have a length equal to the number of dimensions of the function being
     *                   evaluated, otherwise and exception will be raised when searching for the
     *                   minimum.
     * @param directions Set of directions to start looking for a minimum.
     *                   Provided matrix must have the same rows as the number of dimensions of
     *                   the function being evaluated. Each column will be a direction to search
     *                   for a minimum.
     * @param tolerance  Tolerance or accuracy to be expected on estimated local
     *                   minimum.
     * @throws IllegalArgumentException Raised if provided point and direction
     *                                  don't have the same length or if provided tolerance is negative.
     */
    public PowellMultiOptimizer(
            final MultiDimensionFunctionEvaluatorListener listener, final double[] point, final Matrix directions,
            final double tolerance) {
        super(listener);
        internalSetPointAndDirections(point, directions);
        internalSetTolerance(tolerance);
        iter = 0;
        fret = 0.0;
    }

    /**
     * Returns set of directions to start looking for a minimum.
     * Returned matrix will have the same number of rows as the number of
     * dimensions of the function (or the start point).
     * Each column of the matrix will represent a vector containing a direction
     * to search for a minimum.
     *
     * @return Set of directions.
     * @throws NotAvailableException Raised if this has not yet been provided or
     *                               computed.
     */
    public Matrix getDirections() throws NotAvailableException {
        if (!areDirectionsAvailable()) {
            throw new NotAvailableException();
        }
        return ximat;
    }

    /**
     * Returns boolean indicating whether set of directions is available for
     * retrieval.
     *
     * @return True if available, false otherwise.
     */
    public boolean areDirectionsAvailable() {
        return ximat != null;
    }

    /**
     * Sets start point and set of directions to start looking for minimum.
     *
     * @param point      Start point where algorithm will be started. Start point
     *                   should be close to the local minimum to be found. Provided array must
     *                   have a length equal to the number of dimensions of the function being
     *                   evaluated, otherwise and exception will be raised when searching for the
     *                   minimum.
     * @param directions Set of directions to start looking for a minimum.
     *                   Provided matrix must have the same rows as the number of dimensions of
     *                   the function being evaluated. Each column will be a direction to search
     *                   for a minimum.
     * @throws LockedException          Raised if this instance is locked.
     * @throws IllegalArgumentException Raised if provided point and direction
     *                                  don't have the same length.
     */
    public void setPointAndDirections(final double[] point, final Matrix directions) throws LockedException {
        if (isLocked()) {
            throw new LockedException();
        }
        internalSetPointAndDirections(point, directions);
    }

    /**
     * Sets start point where local minimum is searched nearby.
     *
     * @param startPoint Start point to search for a local minimum
     * @throws LockedException Raised if this instance is locked.
     */
    public void setStartPoint(final double[] startPoint) throws LockedException {
        if (isLocked()) {
            throw new LockedException();
        }
        p = startPoint;
    }

    /**
     * Returns tolerance or accuracy to be expected on estimated local minimum.
     *
     * @return Tolerance or accuracy to be expected on estimated local minimum.
     */
    public double getTolerance() {
        return tolerance;
    }

    /**
     * Internal method to set tolerance or accuracy to be expected on estimated
     * local minimum.
     * This method does not check whether this instance is locked.
     *
     * @param tolerance Tolerance or accuracy to be expected on estimated local
     *                  minimum.
     * @throws IllegalArgumentException Raised if provided tolerance is
     *                                  negative.
     */
    private void internalSetTolerance(final double tolerance) {
        if (tolerance < MIN_TOLERANCE) {
            throw new IllegalArgumentException();
        }
        this.tolerance = tolerance;
    }

    /**
     * Sets tolerance or accuracy to be expected on estimated local minimum.
     *
     * @param tolerance Tolerance or accuracy to be expected on estimated local
     *                  minimum.
     * @throws LockedException          Raised if this instance is locked.
     * @throws IllegalArgumentException Raised if provided tolerance is
     *                                  negative.
     */
    public void setTolerance(final double tolerance) throws LockedException {
        if (isLocked()) {
            throw new LockedException();
        }
        internalSetTolerance(tolerance);
    }

    /**
     * This function estimates a function minimum.
     * Implementations of this class will usually search a local minimum close
     * to a start point and will start looking into provided start direction.
     * Minimization of a function f. Input consists of an initial starting
     * point p. The initial matrix ximat, whose columns contain the initial
     * set of directions, is set to the identity. Returned is the best point
     * found, at which point fret is the minimum function value and iter is
     * the number of iterations taken.
     *
     * @throws LockedException       Raised if this instance is locked, because
     *                               estimation is being computed.
     * @throws NotReadyException     Raised if this instance is not ready, because
     *                               a listener, a gradient listener and a start point haven't been provided.
     * @throws OptimizationException Raised if the algorithm failed because of
     *                               lack of convergence or because function couldn't be evaluated.
     */
    @Override
    public void minimize() throws LockedException, NotReadyException, OptimizationException {

        if (isLocked()) {
            throw new LockedException();
        }
        if (!isReady()) {
            throw new NotReadyException();
        }

        locked = true;

        buildDirections();

        // in case that set of directions have been directly provided, we check
        // for correctness
        final int n = p.length;
        if (ximat.getRows() != ximat.getColumns() || ximat.getRows() != n) {
            locked = false;
            throw new OptimizationException();
        }

        double fptt;
        final var pt = new double[n];
        final var ptt = new double[n];

        // set vector of directions
        if (!isDirectionAvailable() || xi.length != n) {
            xi = new double[n];
        }

        try {
            fret = listener.evaluate(p);

            // save the initial point
            System.arraycopy(p, 0, pt, 0, n);

            for (iter = 0; ; ++iter) {
                final var fp = fret;
                var ibig = 0;
                // Will be the biggest function decrease
                var del = 0.0;
                // In each iteration, loop over all directions in the set.
                for (var i = 0; i < n; i++) {
                    // Copy the direction contained in i-th column of ximat
                    ximat.getSubmatrixAsArray(0, i, n - 1, i, xi);
                    fptt = fret;

                    // minimize along it, and record it if it is the largest
                    // decrease so far.
                    fret = linmin();
                    if (fptt - fret > del) {
                        del = fptt - fret;
                        ibig = i + 1;
                    }
                }

                // Here comes the termination criterion.
                if (2.0 * (fp - fret) <= tolerance * (Math.abs(fp) +
                        Math.abs(fret)) + TINY) {
                    // minimum has been found
                    break;
                }
                if (iter == ITMAX) {
                    // too many iterations
                    locked = false;
                    throw new OptimizationException();
                }

                // Construct the extrapolated point and the average direction
                // moved. Save the old starting point
                for (var j = 0; j < n; j++) {
                    ptt[j] = 2.0 * p[j] - pt[j];
                    xi[j] = p[j] - pt[j];
                    pt[j] = p[j];
                }

                // Function value at extrapolated point
                fptt = listener.evaluate(ptt);

                if (fptt < fp) {
                    final var t = 2.0 * (fp - 2.0 * fret + fptt) * sqr(fp - fret - del) - del * sqr(fp - fptt);

                    if (t < 0.0) {
                        // Move to the minimum of the new direction, and save the
                        // new direction
                        fret = linmin();

                        ximat.setSubmatrix(0, ibig - 1,
                                n - 1, ibig - 1, ximat,
                                0, n - 1,
                                n - 1, n - 1);
                        ximat.setSubmatrix(0, n - 1,
                                n - 1, n - 1, xi);
                    }
                }

                if (iterationCompletedListener != null) {
                    iterationCompletedListener.onIterationCompleted(this, iter, ITMAX);
                }
            }
        } catch (final AlgebraException | EvaluationException e) {
            throw new OptimizationException(e);
        } finally {
            locked = false;
        }

        // set result
        xmin = p;
        resultAvailable = true;
        fmin = fret;
    }

    /**
     * Returns boolean indicating whether this instance is ready to start the
     * estimation of a local minimum.
     * An instance is ready once a listener and a start point are provided.
     *
     * @return True if this instance is ready, false otherwise.
     */
    @Override
    public boolean isReady() {
        return isListenerAvailable() && isStartPointAvailable();
    }

    /**
     * Internal method to set start point and set of directions to start looking
     * for minimum.
     * This method does not check whether this instance is locked.
     *
     * @param point      Start point where algorithm will be started. Start point
     *                   should be close to the local minimum to be found. Provided array must
     *                   have a length equal to the number of dimensions of the function being
     *                   evaluated, otherwise and exception will be raised when searching for the
     *                   minimum.
     * @param directions Set of directions to start looking for a minimum.
     *                   Provided matrix must have the same rows as the number of dimensions of
     *                   the function being evaluated. Each column will be a direction to search
     *                   for a minimum.
     * @throws IllegalArgumentException Raised if provided point and direction
     *                                  don't have the same length.
     */
    private void internalSetPointAndDirections(final double[] point, final Matrix directions) {
        if ((point.length != directions.getRows()) || (point.length != directions.getColumns())) {
            throw new IllegalArgumentException();
        }
        p = point;
        ximat = directions;
        xi = null;
    }

    /**
     * Internal method to build or rebuild the set of directions if needed.
     *
     * @throws NotReadyException Raised if no start point has yet been provided.
     */
    private void buildDirections() throws NotReadyException {
        if (!isStartPointAvailable()) {
            throw new NotReadyException();
        }

        final var n = p.length;

        if (areDirectionsAvailable() && (ximat.getRows() == n) && (ximat.getColumns() == n)) {
            return;
        }

        try {
            ximat = Matrix.identity(n, n);
        } catch (final WrongSizeException ignore) {
            // never happens
        }
    }

    /**
     * Computes the squared value of provided double.
     *
     * @param x Value to be squared.
     * @return Squared value.
     */
    private double sqr(double x) {
        return x * x;
    }
}