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