    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.optim.nonlinear.scalar.noderiv;
    019
    020    import java.util.Arrays;
    021    import java.util.Comparator;
    022
    023    import org.apache.commons.math3.analysis.MultivariateFunction;
    024    import org.apache.commons.math3.exception.NotStrictlyPositiveException;
    025    import org.apache.commons.math3.exception.DimensionMismatchException;
    026    import org.apache.commons.math3.exception.ZeroException;
    027    import org.apache.commons.math3.exception.OutOfRangeException;
    028    import org.apache.commons.math3.exception.NullArgumentException;
    029    import org.apache.commons.math3.exception.MathIllegalArgumentException;
    030    import org.apache.commons.math3.exception.util.LocalizedFormats;
    031    import org.apache.commons.math3.optim.PointValuePair;
    032    import org.apache.commons.math3.optim.OptimizationData;
    033
    034    /**
    035     * This class implements the simplex concept.
    036     * It is intended to be used in conjunction with {@link SimplexOptimizer}.
    037     * <br/>
    038     * The initial configuration of the simplex is set by the constructors
    039     * {@link #AbstractSimplex(double[])} or {@link #AbstractSimplex(double[][])}.
    040     * The other {@link #AbstractSimplex(int) constructor} will set all steps
    041     * to 1, thus building a default configuration from a unit hypercube.
    042     * <br/>
    043     * Users <em>must</em> call the {@link #build(double[]) build} method in order
    044     * to create the data structure that will be acted on by the other methods of
    045     * this class.
    046     *
    047     * @see SimplexOptimizer
    048     * @version $Id: AbstractSimplex.java 1397759 2012-10-13 01:12:58Z erans $
    049     * @since 3.0
    050     */
    051    public abstract class AbstractSimplex implements OptimizationData {
    052        /** Simplex. */
    053        private PointValuePair[] simplex;
    054        /** Start simplex configuration. */
    055        private double[][] startConfiguration;
    056        /** Simplex dimension (must be equal to {@code simplex.length - 1}). */
    057        private final int dimension;
    058
    059        /**
    060         * Build a unit hypercube simplex.
    061         *
    062         * @param n Dimension of the simplex.
    063         */
    064        protected AbstractSimplex(int n) {
    065            this(n, 1d);
    066        }
    067
    068        /**
    069         * Build a hypercube simplex with the given side length.
    070         *
    071         * @param n Dimension of the simplex.
    072         * @param sideLength Length of the sides of the hypercube.
    073         */
    074        protected AbstractSimplex(int n,
    075                                  double sideLength) {
    076            this(createHypercubeSteps(n, sideLength));
    077        }
    078
    079        /**
    080         * The start configuration for simplex is built from a box parallel to
    081         * the canonical axes of the space. The simplex is the subset of vertices
    082         * of a box parallel to the canonical axes. It is built as the path followed
    083         * while traveling from one vertex of the box to the diagonally opposite
    084         * vertex moving only along the box edges. The first vertex of the box will
    085         * be located at the start point of the optimization.
    086         * As an example, in dimension 3 a simplex has 4 vertices. Setting the
    087         * steps to (1, 10, 2) and the start point to (1, 1, 1) would imply the
    088         * start simplex would be: { (1, 1, 1), (2, 1, 1), (2, 11, 1), (2, 11, 3) }.
    089         * The first vertex would be set to the start point at (1, 1, 1) and the
    090         * last vertex would be set to the diagonally opposite vertex at (2, 11, 3).
    091         *
    092         * @param steps Steps along the canonical axes representing box edges. They
    093         * may be negative but not zero.
    094         * @throws NullArgumentException if {@code steps} is {@code null}.
    095         * @throws ZeroException if one of the steps is zero.
    096         */
    097        protected AbstractSimplex(final double[] steps) {
    098            if (steps == null) {
    099                throw new NullArgumentException();
    100            }
    101            if (steps.length == 0) {
    102                throw new ZeroException();
    103            }
    104            dimension = steps.length;
    105
    106            // Only the relative position of the n final vertices with respect
    107            // to the first one are stored.
    108            startConfiguration = new double[dimension][dimension];
    109            for (int i = 0; i < dimension; i++) {
    110                final double[] vertexI = startConfiguration[i];
    111                for (int j = 0; j < i + 1; j++) {
    112                    if (steps[j] == 0) {
    113                        throw new ZeroException(LocalizedFormats.EQUAL_VERTICES_IN_SIMPLEX);
    114                    }
    115                    System.arraycopy(steps, 0, vertexI, 0, j + 1);
    116                }
    117            }
    118        }
    119
    120        /**
    121         * The real initial simplex will be set up by moving the reference
    122         * simplex such that its first point is located at the start point of the
    123         * optimization.
    124         *
    125         * @param referenceSimplex Reference simplex.
    126         * @throws NotStrictlyPositiveException if the reference simplex does not
    127         * contain at least one point.
    128         * @throws DimensionMismatchException if there is a dimension mismatch
    129         * in the reference simplex.
    130         * @throws IllegalArgumentException if one of its vertices is duplicated.
    131         */
    132        protected AbstractSimplex(final double[][] referenceSimplex) {
    133            if (referenceSimplex.length <= 0) {
    134                throw new NotStrictlyPositiveException(LocalizedFormats.SIMPLEX_NEED_ONE_POINT,
    135                                                       referenceSimplex.length);
    136            }
    137            dimension = referenceSimplex.length - 1;
    138
    139            // Only the relative position of the n final vertices with respect
    140            // to the first one are stored.
    141            startConfiguration = new double[dimension][dimension];
    142            final double[] ref0 = referenceSimplex[0];
    143
    144            // Loop over vertices.
    145            for (int i = 0; i < referenceSimplex.length; i++) {
    146                final double[] refI = referenceSimplex[i];
    147
    148                // Safety checks.
    149                if (refI.length != dimension) {
    150                    throw new DimensionMismatchException(refI.length, dimension);
    151                }
    152                for (int j = 0; j < i; j++) {
    153                    final double[] refJ = referenceSimplex[j];
    154                    boolean allEquals = true;
    155                    for (int k = 0; k < dimension; k++) {
    156                        if (refI[k] != refJ[k]) {
    157                            allEquals = false;
    158                            break;
    159                        }
    160                    }
    161                    if (allEquals) {
    162                        throw new MathIllegalArgumentException(LocalizedFormats.EQUAL_VERTICES_IN_SIMPLEX,
    163                                                               i, j);
    164                    }
    165                }
    166
    167                // Store vertex i position relative to vertex 0 position.
    168                if (i > 0) {
    169                    final double[] confI = startConfiguration[i - 1];
    170                    for (int k = 0; k < dimension; k++) {
    171                        confI[k] = refI[k] - ref0[k];
    172                    }
    173                }
    174            }
    175        }
    176
    177        /**
    178         * Get simplex dimension.
    179         *
    180         * @return the dimension of the simplex.
    181         */
    182        public int getDimension() {
    183            return dimension;
    184        }
    185
    186        /**
    187         * Get simplex size.
    188         * After calling the {@link #build(double[]) build} method, this method will
    189         * will be equivalent to {@code getDimension() + 1}.
    190         *
    191         * @return the size of the simplex.
    192         */
    193        public int getSize() {
    194            return simplex.length;
    195        }
    196
    197        /**
    198         * Compute the next simplex of the algorithm.
    199         *
    200         * @param evaluationFunction Evaluation function.
    201         * @param comparator Comparator to use to sort simplex vertices from best
    202         * to worst.
    203         * @throws org.apache.commons.math3.exception.TooManyEvaluationsException
    204         * if the algorithm fails to converge.
    205         */
    206        public abstract void iterate(final MultivariateFunction evaluationFunction,
    207                                     final Comparator<PointValuePair> comparator);
    208
    209        /**
    210         * Build an initial simplex.
    211         *
    212         * @param startPoint First point of the simplex.
    213         * @throws DimensionMismatchException if the start point does not match
    214         * simplex dimension.
    215         */
    216        public void build(final double[] startPoint) {
    217            if (dimension != startPoint.length) {
    218                throw new DimensionMismatchException(dimension, startPoint.length);
    219            }
    220
    221            // Set first vertex.
    222            simplex = new PointValuePair[dimension + 1];
    223            simplex[0] = new PointValuePair(startPoint, Double.NaN);
    224
    225            // Set remaining vertices.
    226            for (int i = 0; i < dimension; i++) {
    227                final double[] confI = startConfiguration[i];
    228                final double[] vertexI = new double[dimension];
    229                for (int k = 0; k < dimension; k++) {
    230                    vertexI[k] = startPoint[k] + confI[k];
    231                }
    232                simplex[i + 1] = new PointValuePair(vertexI, Double.NaN);
    233            }
    234        }
    235
    236        /**
    237         * Evaluate all the non-evaluated points of the simplex.
    238         *
    239         * @param evaluationFunction Evaluation function.
    240         * @param comparator Comparator to use to sort simplex vertices from best to worst.
    241         * @throws org.apache.commons.math3.exception.TooManyEvaluationsException
    242         * if the maximal number of evaluations is exceeded.
    243         */
    244        public void evaluate(final MultivariateFunction evaluationFunction,
    245                             final Comparator<PointValuePair> comparator) {
    246            // Evaluate the objective function at all non-evaluated simplex points.
    247            for (int i = 0; i < simplex.length; i++) {
    248                final PointValuePair vertex = simplex[i];
    249                final double[] point = vertex.getPointRef();
    250                if (Double.isNaN(vertex.getValue())) {
    251                    simplex[i] = new PointValuePair(point, evaluationFunction.value(point), false);
    252                }
    253            }
    254
    255            // Sort the simplex from best to worst.
    256            Arrays.sort(simplex, comparator);
    257        }
    258
    259        /**
    260         * Replace the worst point of the simplex by a new point.
    261         *
    262         * @param pointValuePair Point to insert.
    263         * @param comparator Comparator to use for sorting the simplex vertices
    264         * from best to worst.
    265         */
    266        protected void replaceWorstPoint(PointValuePair pointValuePair,
    267                                         final Comparator<PointValuePair> comparator) {
    268            for (int i = 0; i < dimension; i++) {
    269                if (comparator.compare(simplex[i], pointValuePair) > 0) {
    270                    PointValuePair tmp = simplex[i];
    271                    simplex[i] = pointValuePair;
    272                    pointValuePair = tmp;
    273                }
    274            }
    275            simplex[dimension] = pointValuePair;
    276        }
    277
    278        /**
    279         * Get the points of the simplex.
    280         *
    281         * @return all the simplex points.
    282         */
    283        public PointValuePair[] getPoints() {
    284            final PointValuePair[] copy = new PointValuePair[simplex.length];
    285            System.arraycopy(simplex, 0, copy, 0, simplex.length);
    286            return copy;
    287        }
    288
    289        /**
    290         * Get the simplex point stored at the requested {@code index}.
    291         *
    292         * @param index Location.
    293         * @return the point at location {@code index}.
    294         */
    295        public PointValuePair getPoint(int index) {
    296            if (index < 0 ||
    297                index >= simplex.length) {
    298                throw new OutOfRangeException(index, 0, simplex.length - 1);
    299            }
    300            return simplex[index];
    301        }
    302
    303        /**
    304         * Store a new point at location {@code index}.
    305         * Note that no deep-copy of {@code point} is performed.
    306         *
    307         * @param index Location.
    308         * @param point New value.
    309         */
    310        protected void setPoint(int index, PointValuePair point) {
    311            if (index < 0 ||
    312                index >= simplex.length) {
    313                throw new OutOfRangeException(index, 0, simplex.length - 1);
    314            }
    315            simplex[index] = point;
    316        }
    317
    318        /**
    319         * Replace all points.
    320         * Note that no deep-copy of {@code points} is performed.
    321         *
    322         * @param points New Points.
    323         */
    324        protected void setPoints(PointValuePair[] points) {
    325            if (points.length != simplex.length) {
    326                throw new DimensionMismatchException(points.length, simplex.length);
    327            }
    328            simplex = points;
    329        }
    330
    331        /**
    332         * Create steps for a unit hypercube.
    333         *
    334         * @param n Dimension of the hypercube.
    335         * @param sideLength Length of the sides of the hypercube.
    336         * @return the steps.
    337         */
    338        private static double[] createHypercubeSteps(int n,
    339                                                     double sideLength) {
    340            final double[] steps = new double[n];
    341            for (int i = 0; i < n; i++) {
    342                steps[i] = sideLength;
    343            }
    344            return steps;
    345        }
    346    }