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