View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  
18  package org.apache.commons.math3.optim.nonlinear.scalar.noderiv;
19  
20  import java.util.Arrays;
21  import java.util.Comparator;
22  
23  import org.apache.commons.math3.analysis.MultivariateFunction;
24  import org.apache.commons.math3.exception.NotStrictlyPositiveException;
25  import org.apache.commons.math3.exception.DimensionMismatchException;
26  import org.apache.commons.math3.exception.ZeroException;
27  import org.apache.commons.math3.exception.OutOfRangeException;
28  import org.apache.commons.math3.exception.NullArgumentException;
29  import org.apache.commons.math3.exception.MathIllegalArgumentException;
30  import org.apache.commons.math3.exception.util.LocalizedFormats;
31  import org.apache.commons.math3.optim.PointValuePair;
32  import org.apache.commons.math3.optim.OptimizationData;
33  
34  /**
35   * This class implements the simplex concept.
36   * It is intended to be used in conjunction with {@link SimplexOptimizer}.
37   * <br/>
38   * The initial configuration of the simplex is set by the constructors
39   * {@link #AbstractSimplex(double[])} or {@link #AbstractSimplex(double[][])}.
40   * The other {@link #AbstractSimplex(int) constructor} will set all steps
41   * to 1, thus building a default configuration from a unit hypercube.
42   * <br/>
43   * Users <em>must</em> call the {@link #build(double[]) build} method in order
44   * to create the data structure that will be acted on by the other methods of
45   * this class.
46   *
47   * @see SimplexOptimizer
48   * @version $Id: AbstractSimplex.java 1435539 2013-01-19 13:27:24Z tn $
49   * @since 3.0
50   */
51  public abstract class AbstractSimplex implements OptimizationData {
52      /** Simplex. */
53      private PointValuePair[] simplex;
54      /** Start simplex configuration. */
55      private double[][] startConfiguration;
56      /** Simplex dimension (must be equal to {@code simplex.length - 1}). */
57      private final int dimension;
58  
59      /**
60       * Build a unit hypercube simplex.
61       *
62       * @param n Dimension of the simplex.
63       */
64      protected AbstractSimplex(int n) {
65          this(n, 1d);
66      }
67  
68      /**
69       * Build a hypercube simplex with the given side length.
70       *
71       * @param n Dimension of the simplex.
72       * @param sideLength Length of the sides of the hypercube.
73       */
74      protected AbstractSimplex(int n,
75                                double sideLength) {
76          this(createHypercubeSteps(n, sideLength));
77      }
78  
79      /**
80       * The start configuration for simplex is built from a box parallel to
81       * the canonical axes of the space. The simplex is the subset of vertices
82       * of a box parallel to the canonical axes. It is built as the path followed
83       * while traveling from one vertex of the box to the diagonally opposite
84       * vertex moving only along the box edges. The first vertex of the box will
85       * be located at the start point of the optimization.
86       * As an example, in dimension 3 a simplex has 4 vertices. Setting the
87       * steps to (1, 10, 2) and the start point to (1, 1, 1) would imply the
88       * start simplex would be: { (1, 1, 1), (2, 1, 1), (2, 11, 1), (2, 11, 3) }.
89       * The first vertex would be set to the start point at (1, 1, 1) and the
90       * last vertex would be set to the diagonally opposite vertex at (2, 11, 3).
91       *
92       * @param steps Steps along the canonical axes representing box edges. They
93       * may be negative but not zero.
94       * @throws NullArgumentException if {@code steps} is {@code null}.
95       * @throws ZeroException if one of the steps is zero.
96       */
97      protected AbstractSimplex(final double[] steps) {
98          if (steps == null) {
99              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 }