View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      https://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  
18  /*
19   * This is not the original file distributed by the Apache Software Foundation
20   * It has been modified by the Hipparchus project
21   */
22  package org.hipparchus.optim.nonlinear.scalar.noderiv;
23  
24  import java.util.Comparator;
25  
26  import org.hipparchus.analysis.MultivariateFunction;
27  import org.hipparchus.optim.PointValuePair;
28  
29  /**
30   * This class implements the Nelder-Mead simplex algorithm.
31   *
32   */
33  public class NelderMeadSimplex extends AbstractSimplex {
34      /** Default value for {@link #rho}: {@value}. */
35      private static final double DEFAULT_RHO = 1;
36      /** Default value for {@link #khi}: {@value}. */
37      private static final double DEFAULT_KHI = 2;
38      /** Default value for {@link #gamma}: {@value}. */
39      private static final double DEFAULT_GAMMA = 0.5;
40      /** Default value for {@link #sigma}: {@value}. */
41      private static final double DEFAULT_SIGMA = 0.5;
42      /** Reflection coefficient. */
43      private final double rho;
44      /** Expansion coefficient. */
45      private final double khi;
46      /** Contraction coefficient. */
47      private final double gamma;
48      /** Shrinkage coefficient. */
49      private final double sigma;
50  
51      /**
52       * Build a Nelder-Mead simplex with default coefficients.
53       * The default coefficients are 1.0 for rho, 2.0 for khi and 0.5
54       * for both gamma and sigma.
55       *
56       * @param n Dimension of the simplex.
57       */
58      public NelderMeadSimplex(final int n) {
59          this(n, 1d);
60      }
61  
62      /**
63       * Build a Nelder-Mead simplex with default coefficients.
64       * The default coefficients are 1.0 for rho, 2.0 for khi and 0.5
65       * for both gamma and sigma.
66       *
67       * @param n Dimension of the simplex.
68       * @param sideLength Length of the sides of the default (hypercube)
69       * simplex. See {@link AbstractSimplex#AbstractSimplex(int,double)}.
70       */
71      public NelderMeadSimplex(final int n, double sideLength) {
72          this(n, sideLength,
73               DEFAULT_RHO, DEFAULT_KHI, DEFAULT_GAMMA, DEFAULT_SIGMA);
74      }
75  
76      /**
77       * Build a Nelder-Mead simplex with specified coefficients.
78       *
79       * @param n Dimension of the simplex. See
80       * {@link AbstractSimplex#AbstractSimplex(int,double)}.
81       * @param sideLength Length of the sides of the default (hypercube)
82       * simplex. See {@link AbstractSimplex#AbstractSimplex(int,double)}.
83       * @param rho Reflection coefficient.
84       * @param khi Expansion coefficient.
85       * @param gamma Contraction coefficient.
86       * @param sigma Shrinkage coefficient.
87       */
88      public NelderMeadSimplex(final int n, double sideLength,
89                               final double rho, final double khi,
90                               final double gamma, final double sigma) {
91          super(n, sideLength);
92  
93          this.rho = rho;
94          this.khi = khi;
95          this.gamma = gamma;
96          this.sigma = sigma;
97      }
98  
99      /**
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.hipparchus.exception.MathIllegalArgumentException
172      * if the reference simplex does not contain at least one point.
173      * @throws org.hipparchus.exception.MathIllegalArgumentException
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 }