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   * @since 3.0
49   */
50  public abstract class AbstractSimplex implements OptimizationData {
51      /** Simplex. */
52      private PointValuePair[] simplex;
53      /** Start simplex configuration. */
54      private double[][] startConfiguration;
55      /** Simplex dimension (must be equal to {@code simplex.length - 1}). */
56      private final int dimension;
57  
58      /**
59       * Build a unit hypercube simplex.
60       *
61       * @param n Dimension of the simplex.
62       */
63      protected AbstractSimplex(int n) {
64          this(n, 1d);
65      }
66  
67      /**
68       * Build a hypercube simplex with the given side length.
69       *
70       * @param n Dimension of the simplex.
71       * @param sideLength Length of the sides of the hypercube.
72       */
73      protected AbstractSimplex(int n,
74                                double sideLength) {
75          this(createHypercubeSteps(n, sideLength));
76      }
77  
78      /**
79       * The start configuration for simplex is built from a box parallel to
80       * the canonical axes of the space. The simplex is the subset of vertices
81       * of a box parallel to the canonical axes. It is built as the path followed
82       * while traveling from one vertex of the box to the diagonally opposite
83       * vertex moving only along the box edges. The first vertex of the box will
84       * be located at the start point of the optimization.
85       * As an example, in dimension 3 a simplex has 4 vertices. Setting the
86       * steps to (1, 10, 2) and the start point to (1, 1, 1) would imply the
87       * start simplex would be: { (1, 1, 1), (2, 1, 1), (2, 11, 1), (2, 11, 3) }.
88       * The first vertex would be set to the start point at (1, 1, 1) and the
89       * last vertex would be set to the diagonally opposite vertex at (2, 11, 3).
90       *
91       * @param steps Steps along the canonical axes representing box edges. They
92       * may be negative but not zero.
93       * @throws NullArgumentException if {@code steps} is {@code null}.
94       * @throws ZeroException if one of the steps is zero.
95       */
96      protected AbstractSimplex(final double[] steps) {
97          if (steps == null) {
98              throw new NullArgumentException();
99          }
100         if (steps.length == 0) {
101             throw new ZeroException();
102         }
103         dimension = steps.length;
104 
105         // Only the relative position of the n final vertices with respect
106         // to the first one are stored.
107         startConfiguration = new double[dimension][dimension];
108         for (int i = 0; i < dimension; i++) {
109             final double[] vertexI = startConfiguration[i];
110             for (int j = 0; j < i + 1; j++) {
111                 if (steps[j] == 0) {
112                     throw new ZeroException(LocalizedFormats.EQUAL_VERTICES_IN_SIMPLEX);
113                 }
114                 System.arraycopy(steps, 0, vertexI, 0, j + 1);
115             }
116         }
117     }
118 
119     /**
120      * The real initial simplex will be set up by moving the reference
121      * simplex such that its first point is located at the start point of the
122      * optimization.
123      *
124      * @param referenceSimplex Reference simplex.
125      * @throws NotStrictlyPositiveException if the reference simplex does not
126      * contain at least one point.
127      * @throws DimensionMismatchException if there is a dimension mismatch
128      * in the reference simplex.
129      * @throws IllegalArgumentException if one of its vertices is duplicated.
130      */
131     protected AbstractSimplex(final double[][] referenceSimplex) {
132         if (referenceSimplex.length <= 0) {
133             throw new NotStrictlyPositiveException(LocalizedFormats.SIMPLEX_NEED_ONE_POINT,
134                                                    referenceSimplex.length);
135         }
136         dimension = referenceSimplex.length - 1;
137 
138         // Only the relative position of the n final vertices with respect
139         // to the first one are stored.
140         startConfiguration = new double[dimension][dimension];
141         final double[] ref0 = referenceSimplex[0];
142 
143         // Loop over vertices.
144         for (int i = 0; i < referenceSimplex.length; i++) {
145             final double[] refI = referenceSimplex[i];
146 
147             // Safety checks.
148             if (refI.length != dimension) {
149                 throw new DimensionMismatchException(refI.length, dimension);
150             }
151             for (int j = 0; j < i; j++) {
152                 final double[] refJ = referenceSimplex[j];
153                 boolean allEquals = true;
154                 for (int k = 0; k < dimension; k++) {
155                     if (refI[k] != refJ[k]) {
156                         allEquals = false;
157                         break;
158                     }
159                 }
160                 if (allEquals) {
161                     throw new MathIllegalArgumentException(LocalizedFormats.EQUAL_VERTICES_IN_SIMPLEX,
162                                                            i, j);
163                 }
164             }
165 
166             // Store vertex i position relative to vertex 0 position.
167             if (i > 0) {
168                 final double[] confI = startConfiguration[i - 1];
169                 for (int k = 0; k < dimension; k++) {
170                     confI[k] = refI[k] - ref0[k];
171                 }
172             }
173         }
174     }
175 
176     /**
177      * Get simplex dimension.
178      *
179      * @return the dimension of the simplex.
180      */
181     public int getDimension() {
182         return dimension;
183     }
184 
185     /**
186      * Get simplex size.
187      * After calling the {@link #build(double[]) build} method, this method will
188      * will be equivalent to {@code getDimension() + 1}.
189      *
190      * @return the size of the simplex.
191      */
192     public int getSize() {
193         return simplex.length;
194     }
195 
196     /**
197      * Compute the next simplex of the algorithm.
198      *
199      * @param evaluationFunction Evaluation function.
200      * @param comparator Comparator to use to sort simplex vertices from best
201      * to worst.
202      * @throws org.apache.commons.math3.exception.TooManyEvaluationsException
203      * if the algorithm fails to converge.
204      */
205     public abstract void iterate(final MultivariateFunction evaluationFunction,
206                                  final Comparator<PointValuePair> comparator);
207 
208     /**
209      * Build an initial simplex.
210      *
211      * @param startPoint First point of the simplex.
212      * @throws DimensionMismatchException if the start point does not match
213      * simplex dimension.
214      */
215     public void build(final double[] startPoint) {
216         if (dimension != startPoint.length) {
217             throw new DimensionMismatchException(dimension, startPoint.length);
218         }
219 
220         // Set first vertex.
221         simplex = new PointValuePair[dimension + 1];
222         simplex[0] = new PointValuePair(startPoint, Double.NaN);
223 
224         // Set remaining vertices.
225         for (int i = 0; i < dimension; i++) {
226             final double[] confI = startConfiguration[i];
227             final double[] vertexI = new double[dimension];
228             for (int k = 0; k < dimension; k++) {
229                 vertexI[k] = startPoint[k] + confI[k];
230             }
231             simplex[i + 1] = new PointValuePair(vertexI, Double.NaN);
232         }
233     }
234 
235     /**
236      * Evaluate all the non-evaluated points of the simplex.
237      *
238      * @param evaluationFunction Evaluation function.
239      * @param comparator Comparator to use to sort simplex vertices from best to worst.
240      * @throws org.apache.commons.math3.exception.TooManyEvaluationsException
241      * if the maximal number of evaluations is exceeded.
242      */
243     public void evaluate(final MultivariateFunction evaluationFunction,
244                          final Comparator<PointValuePair> comparator) {
245         // Evaluate the objective function at all non-evaluated simplex points.
246         for (int i = 0; i < simplex.length; i++) {
247             final PointValuePair vertex = simplex[i];
248             final double[] point = vertex.getPointRef();
249             if (Double.isNaN(vertex.getValue())) {
250                 simplex[i] = new PointValuePair(point, evaluationFunction.value(point), false);
251             }
252         }
253 
254         // Sort the simplex from best to worst.
255         Arrays.sort(simplex, comparator);
256     }
257 
258     /**
259      * Replace the worst point of the simplex by a new point.
260      *
261      * @param pointValuePair Point to insert.
262      * @param comparator Comparator to use for sorting the simplex vertices
263      * from best to worst.
264      */
265     protected void replaceWorstPoint(PointValuePair pointValuePair,
266                                      final Comparator<PointValuePair> comparator) {
267         for (int i = 0; i < dimension; i++) {
268             if (comparator.compare(simplex[i], pointValuePair) > 0) {
269                 PointValuePair tmp = simplex[i];
270                 simplex[i] = pointValuePair;
271                 pointValuePair = tmp;
272             }
273         }
274         simplex[dimension] = pointValuePair;
275     }
276 
277     /**
278      * Get the points of the simplex.
279      *
280      * @return all the simplex points.
281      */
282     public PointValuePair[] getPoints() {
283         final PointValuePair[] copy = new PointValuePair[simplex.length];
284         System.arraycopy(simplex, 0, copy, 0, simplex.length);
285         return copy;
286     }
287 
288     /**
289      * Get the simplex point stored at the requested {@code index}.
290      *
291      * @param index Location.
292      * @return the point at location {@code index}.
293      */
294     public PointValuePair getPoint(int index) {
295         if (index < 0 ||
296             index >= simplex.length) {
297             throw new OutOfRangeException(index, 0, simplex.length - 1);
298         }
299         return simplex[index];
300     }
301 
302     /**
303      * Store a new point at location {@code index}.
304      * Note that no deep-copy of {@code point} is performed.
305      *
306      * @param index Location.
307      * @param point New value.
308      */
309     protected void setPoint(int index, PointValuePair point) {
310         if (index < 0 ||
311             index >= simplex.length) {
312             throw new OutOfRangeException(index, 0, simplex.length - 1);
313         }
314         simplex[index] = point;
315     }
316 
317     /**
318      * Replace all points.
319      * Note that no deep-copy of {@code points} is performed.
320      *
321      * @param points New Points.
322      */
323     protected void setPoints(PointValuePair[] points) {
324         if (points.length != simplex.length) {
325             throw new DimensionMismatchException(points.length, simplex.length);
326         }
327         simplex = points;
328     }
329 
330     /**
331      * Create steps for a unit hypercube.
332      *
333      * @param n Dimension of the hypercube.
334      * @param sideLength Length of the sides of the hypercube.
335      * @return the steps.
336      */
337     private static double[] createHypercubeSteps(int n,
338                                                  double sideLength) {
339         final double[] steps = new double[n];
340         for (int i = 0; i < n; i++) {
341             steps[i] = sideLength;
342         }
343         return steps;
344     }
345 }