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.optimization.direct;
019    
020    import java.util.Comparator;
021    
022    import org.apache.commons.math3.optimization.PointValuePair;
023    import org.apache.commons.math3.analysis.MultivariateFunction;
024    
025    /**
026     * This class implements the Nelder-Mead simplex algorithm.
027     *
028     * @version $Id: NelderMeadSimplex.java 1422230 2012-12-15 12:11:13Z erans $
029     * @deprecated As of 3.1 (to be removed in 4.0).
030     * @since 3.0
031     */
032    @Deprecated
033    public class NelderMeadSimplex extends AbstractSimplex {
034        /** Default value for {@link #rho}: {@value}. */
035        private static final double DEFAULT_RHO = 1;
036        /** Default value for {@link #khi}: {@value}. */
037        private static final double DEFAULT_KHI = 2;
038        /** Default value for {@link #gamma}: {@value}. */
039        private static final double DEFAULT_GAMMA = 0.5;
040        /** Default value for {@link #sigma}: {@value}. */
041        private static final double DEFAULT_SIGMA = 0.5;
042        /** Reflection coefficient. */
043        private final double rho;
044        /** Expansion coefficient. */
045        private final double khi;
046        /** Contraction coefficient. */
047        private final double gamma;
048        /** Shrinkage coefficient. */
049        private final double sigma;
050    
051        /**
052         * Build a Nelder-Mead simplex with default coefficients.
053         * The default coefficients are 1.0 for rho, 2.0 for khi and 0.5
054         * for both gamma and sigma.
055         *
056         * @param n Dimension of the simplex.
057         */
058        public NelderMeadSimplex(final int n) {
059            this(n, 1d);
060        }
061    
062        /**
063         * Build a Nelder-Mead simplex with default coefficients.
064         * The default coefficients are 1.0 for rho, 2.0 for khi and 0.5
065         * for both gamma and sigma.
066         *
067         * @param n Dimension of the simplex.
068         * @param sideLength Length of the sides of the default (hypercube)
069         * simplex. See {@link AbstractSimplex#AbstractSimplex(int,double)}.
070         */
071        public NelderMeadSimplex(final int n, double sideLength) {
072            this(n, sideLength,
073                 DEFAULT_RHO, DEFAULT_KHI, DEFAULT_GAMMA, DEFAULT_SIGMA);
074        }
075    
076        /**
077         * Build a Nelder-Mead simplex with specified coefficients.
078         *
079         * @param n Dimension of the simplex. See
080         * {@link AbstractSimplex#AbstractSimplex(int,double)}.
081         * @param sideLength Length of the sides of the default (hypercube)
082         * simplex. See {@link AbstractSimplex#AbstractSimplex(int,double)}.
083         * @param rho Reflection coefficient.
084         * @param khi Expansion coefficient.
085         * @param gamma Contraction coefficient.
086         * @param sigma Shrinkage coefficient.
087         */
088        public NelderMeadSimplex(final int n, double sideLength,
089                                 final double rho, final double khi,
090                                 final double gamma, final double sigma) {
091            super(n, sideLength);
092    
093            this.rho = rho;
094            this.khi = khi;
095            this.gamma = gamma;
096            this.sigma = sigma;
097        }
098    
099        /**
100         * Build a Nelder-Mead simplex with specified coefficients.
101         *
102         * @param n Dimension of the simplex. See
103         * {@link AbstractSimplex#AbstractSimplex(int)}.
104         * @param rho Reflection coefficient.
105         * @param khi Expansion coefficient.
106         * @param gamma Contraction coefficient.
107         * @param sigma Shrinkage coefficient.
108         */
109        public NelderMeadSimplex(final int n,
110                                 final double rho, final double khi,
111                                 final double gamma, final double sigma) {
112            this(n, 1d, rho, khi, gamma, sigma);
113        }
114    
115        /**
116         * Build a Nelder-Mead simplex with default coefficients.
117         * The default coefficients are 1.0 for rho, 2.0 for khi and 0.5
118         * for both gamma and sigma.
119         *
120         * @param steps Steps along the canonical axes representing box edges.
121         * They may be negative but not zero. See
122         */
123        public NelderMeadSimplex(final double[] steps) {
124            this(steps, DEFAULT_RHO, DEFAULT_KHI, DEFAULT_GAMMA, DEFAULT_SIGMA);
125        }
126    
127        /**
128         * Build a Nelder-Mead simplex with specified coefficients.
129         *
130         * @param steps Steps along the canonical axes representing box edges.
131         * They may be negative but not zero. See
132         * {@link AbstractSimplex#AbstractSimplex(double[])}.
133         * @param rho Reflection coefficient.
134         * @param khi Expansion coefficient.
135         * @param gamma Contraction coefficient.
136         * @param sigma Shrinkage coefficient.
137         * @throws IllegalArgumentException if one of the steps is zero.
138         */
139        public NelderMeadSimplex(final double[] steps,
140                                 final double rho, final double khi,
141                                 final double gamma, final double sigma) {
142            super(steps);
143    
144            this.rho = rho;
145            this.khi = khi;
146            this.gamma = gamma;
147            this.sigma = sigma;
148        }
149    
150        /**
151         * Build a Nelder-Mead simplex with default coefficients.
152         * The default coefficients are 1.0 for rho, 2.0 for khi and 0.5
153         * for both gamma and sigma.
154         *
155         * @param referenceSimplex Reference simplex. See
156         * {@link AbstractSimplex#AbstractSimplex(double[][])}.
157         */
158        public NelderMeadSimplex(final double[][] referenceSimplex) {
159            this(referenceSimplex, DEFAULT_RHO, DEFAULT_KHI, DEFAULT_GAMMA, DEFAULT_SIGMA);
160        }
161    
162        /**
163         * Build a Nelder-Mead simplex with specified coefficients.
164         *
165         * @param referenceSimplex Reference simplex. See
166         * {@link AbstractSimplex#AbstractSimplex(double[][])}.
167         * @param rho Reflection coefficient.
168         * @param khi Expansion coefficient.
169         * @param gamma Contraction coefficient.
170         * @param sigma Shrinkage coefficient.
171         * @throws org.apache.commons.math3.exception.NotStrictlyPositiveException
172         * if the reference simplex does not contain at least one point.
173         * @throws org.apache.commons.math3.exception.DimensionMismatchException
174         * if there is a dimension mismatch in the reference simplex.
175         */
176        public NelderMeadSimplex(final double[][] referenceSimplex,
177                                 final double rho, final double khi,
178                                 final double gamma, final double sigma) {
179            super(referenceSimplex);
180    
181            this.rho = rho;
182            this.khi = khi;
183            this.gamma = gamma;
184            this.sigma = sigma;
185        }
186    
187        /** {@inheritDoc} */
188        @Override
189        public void iterate(final MultivariateFunction evaluationFunction,
190                            final Comparator<PointValuePair> comparator) {
191            // The simplex has n + 1 points if dimension is n.
192            final int n = getDimension();
193    
194            // Interesting values.
195            final PointValuePair best = getPoint(0);
196            final PointValuePair secondBest = getPoint(n - 1);
197            final PointValuePair worst = getPoint(n);
198            final double[] xWorst = worst.getPointRef();
199    
200            // Compute the centroid of the best vertices (dismissing the worst
201            // point at index n).
202            final double[] centroid = new double[n];
203            for (int i = 0; i < n; i++) {
204                final double[] x = getPoint(i).getPointRef();
205                for (int j = 0; j < n; j++) {
206                    centroid[j] += x[j];
207                }
208            }
209            final double scaling = 1.0 / n;
210            for (int j = 0; j < n; j++) {
211                centroid[j] *= scaling;
212            }
213    
214            // compute the reflection point
215            final double[] xR = new double[n];
216            for (int j = 0; j < n; j++) {
217                xR[j] = centroid[j] + rho * (centroid[j] - xWorst[j]);
218            }
219            final PointValuePair reflected
220                = new PointValuePair(xR, evaluationFunction.value(xR), false);
221    
222            if (comparator.compare(best, reflected) <= 0 &&
223                comparator.compare(reflected, secondBest) < 0) {
224                // Accept the reflected point.
225                replaceWorstPoint(reflected, comparator);
226            } else if (comparator.compare(reflected, best) < 0) {
227                // Compute the expansion point.
228                final double[] xE = new double[n];
229                for (int j = 0; j < n; j++) {
230                    xE[j] = centroid[j] + khi * (xR[j] - centroid[j]);
231                }
232                final PointValuePair expanded
233                    = new PointValuePair(xE, evaluationFunction.value(xE), false);
234    
235                if (comparator.compare(expanded, reflected) < 0) {
236                    // Accept the expansion point.
237                    replaceWorstPoint(expanded, comparator);
238                } else {
239                    // Accept the reflected point.
240                    replaceWorstPoint(reflected, comparator);
241                }
242            } else {
243                if (comparator.compare(reflected, worst) < 0) {
244                    // Perform an outside contraction.
245                    final double[] xC = new double[n];
246                    for (int j = 0; j < n; j++) {
247                        xC[j] = centroid[j] + gamma * (xR[j] - centroid[j]);
248                    }
249                    final PointValuePair outContracted
250                        = new PointValuePair(xC, evaluationFunction.value(xC), false);
251                    if (comparator.compare(outContracted, reflected) <= 0) {
252                        // Accept the contraction point.
253                        replaceWorstPoint(outContracted, comparator);
254                        return;
255                    }
256                } else {
257                    // Perform an inside contraction.
258                    final double[] xC = new double[n];
259                    for (int j = 0; j < n; j++) {
260                        xC[j] = centroid[j] - gamma * (centroid[j] - xWorst[j]);
261                    }
262                    final PointValuePair inContracted
263                        = new PointValuePair(xC, evaluationFunction.value(xC), false);
264    
265                    if (comparator.compare(inContracted, worst) < 0) {
266                        // Accept the contraction point.
267                        replaceWorstPoint(inContracted, comparator);
268                        return;
269                    }
270                }
271    
272                // Perform a shrink.
273                final double[] xSmallest = getPoint(0).getPointRef();
274                for (int i = 1; i <= n; i++) {
275                    final double[] x = getPoint(i).getPoint();
276                    for (int j = 0; j < n; j++) {
277                        x[j] = xSmallest[j] + sigma * (x[j] - xSmallest[j]);
278                    }
279                    setPoint(i, new PointValuePair(x, Double.NaN, false));
280                }
281                evaluate(evaluationFunction, comparator);
282            }
283        }
284    }