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
018package org.apache.commons.math4.optim.nonlinear.scalar.noderiv;
019
020import java.util.Arrays;
021import java.util.Comparator;
022
023import org.apache.commons.math4.analysis.MultivariateFunction;
024import org.apache.commons.math4.exception.DimensionMismatchException;
025import org.apache.commons.math4.exception.MathIllegalArgumentException;
026import org.apache.commons.math4.exception.NotStrictlyPositiveException;
027import org.apache.commons.math4.exception.NullArgumentException;
028import org.apache.commons.math4.exception.OutOfRangeException;
029import org.apache.commons.math4.exception.ZeroException;
030import org.apache.commons.math4.exception.util.LocalizedFormats;
031import org.apache.commons.math4.optim.OptimizationData;
032import org.apache.commons.math4.optim.PointValuePair;
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 * @since 3.0
049 */
050public abstract class AbstractSimplex implements OptimizationData {
051    /** Simplex. */
052    private PointValuePair[] simplex;
053    /** Start simplex configuration. */
054    private double[][] startConfiguration;
055    /** Simplex dimension (must be equal to {@code simplex.length - 1}). */
056    private final int dimension;
057
058    /**
059     * Build a unit hypercube simplex.
060     *
061     * @param n Dimension of the simplex.
062     */
063    protected AbstractSimplex(int n) {
064        this(n, 1d);
065    }
066
067    /**
068     * Build a hypercube simplex with the given side length.
069     *
070     * @param n Dimension of the simplex.
071     * @param sideLength Length of the sides of the hypercube.
072     */
073    protected AbstractSimplex(int n,
074                              double sideLength) {
075        this(createHypercubeSteps(n, sideLength));
076    }
077
078    /**
079     * The start configuration for simplex is built from a box parallel to
080     * the canonical axes of the space. The simplex is the subset of vertices
081     * of a box parallel to the canonical axes. It is built as the path followed
082     * while traveling from one vertex of the box to the diagonally opposite
083     * vertex moving only along the box edges. The first vertex of the box will
084     * be located at the start point of the optimization.
085     * As an example, in dimension 3 a simplex has 4 vertices. Setting the
086     * steps to (1, 10, 2) and the start point to (1, 1, 1) would imply the
087     * start simplex would be: { (1, 1, 1), (2, 1, 1), (2, 11, 1), (2, 11, 3) }.
088     * The first vertex would be set to the start point at (1, 1, 1) and the
089     * last vertex would be set to the diagonally opposite vertex at (2, 11, 3).
090     *
091     * @param steps Steps along the canonical axes representing box edges. They
092     * may be negative but not zero.
093     * @throws NullArgumentException if {@code steps} is {@code null}.
094     * @throws ZeroException if one of the steps is zero.
095     */
096    protected AbstractSimplex(final double[] steps) {
097        if (steps == null) {
098            throw new NullArgumentException();
099        }
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.math4.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.math4.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}