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.Comparator;
021
022import org.apache.commons.math3.optimization.PointValuePair;
023import 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
033public 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}