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    }