1 /*
2 * This file is part of the GROMACS molecular simulation package.
3 *
4 * Copyright (c) 2020, by the GROMACS development team, led by
5 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6 * and including many others, as listed in the AUTHORS file in the
7 * top-level source directory and at http://www.gromacs.org.
8 *
9 * GROMACS is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation; either version 2.1
12 * of the License, or (at your option) any later version.
13 *
14 * GROMACS is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Lesser General Public License for more details.
18 *
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with GROMACS; if not, see
21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
23 *
24 * If you want to redistribute modifications to GROMACS, please
25 * consider that scientific software is very special. Version
26 * control is crucial - bugs must be traceable. We will be happy to
27 * consider code for inclusion in the official distribution, but
28 * derived work must not be called official GROMACS. Details are found
29 * in the README & COPYING files - if they are missing, get the
30 * official version at http://www.gromacs.org.
31 *
32 * To help us fund GROMACS development, we humbly ask that you cite
33 * the research papers on the package. Check out http://www.gromacs.org.
34 */
35 /*! \internal \file
36 *
37 * \brief Implements routines in optimization.h .
38 *
39 * \author Christian Blau <blau@kth.se>
40 */
41
42 #include "gmxpre.h"
43
44 #include "optimization.h"
45
46 #include "neldermead.h"
47
48 namespace gmx
49 {
50
nelderMead(const std::function<real (ArrayRef<const real>)> & functionToMinimize,ArrayRef<const real> initalGuess,real minimumRelativeSimplexLength,int maxSteps)51 OptimisationResult nelderMead(const std::function<real(ArrayRef<const real>)>& functionToMinimize,
52 ArrayRef<const real> initalGuess,
53 real minimumRelativeSimplexLength,
54 int maxSteps)
55 {
56 // Set up the initial simplex, sorting vertices according to function value
57 NelderMeadSimplex nelderMeadSimplex(functionToMinimize, initalGuess);
58
59 // Run until maximum step size reached or algorithm is converged, e.g.,
60 // the oriented simplex length is smaller or equal a given number.
61 const real minimumSimplexLength = minimumRelativeSimplexLength * nelderMeadSimplex.orientedLength();
62 for (int currentStep = 0;
63 nelderMeadSimplex.orientedLength() > minimumSimplexLength && currentStep < maxSteps; ++currentStep)
64 {
65
66 // see if simplex can by improved by reflecing the worst vertex at the centroid
67 const RealFunctionvalueAtCoordinate& reflectionPoint =
68 nelderMeadSimplex.evaluateReflectionPoint(functionToMinimize);
69
70 // Reflection point is not better than best simplex vertex so far
71 // but better than second worst
72 if ((nelderMeadSimplex.bestVertex().value_ <= reflectionPoint.value_)
73 && (reflectionPoint.value_ < nelderMeadSimplex.secondWorstValue()))
74 {
75 nelderMeadSimplex.swapOutWorst(reflectionPoint);
76 continue;
77 }
78
79 // If the reflection point is better than the best one see if simplex
80 // can be further improved by continuing going in that direction
81 if (reflectionPoint.value_ < nelderMeadSimplex.bestVertex().value_)
82 {
83 RealFunctionvalueAtCoordinate expansionPoint =
84 nelderMeadSimplex.evaluateExpansionPoint(functionToMinimize);
85 if (expansionPoint.value_ < reflectionPoint.value_)
86 {
87 nelderMeadSimplex.swapOutWorst(expansionPoint);
88 }
89 else
90 {
91 nelderMeadSimplex.swapOutWorst(reflectionPoint);
92 }
93 continue;
94 }
95
96 // The reflection point was a poor choice, try contracting the
97 // worst point coordinates using the centroid instead
98 RealFunctionvalueAtCoordinate contractionPoint =
99 nelderMeadSimplex.evaluateContractionPoint(functionToMinimize);
100 if (contractionPoint.value_ < nelderMeadSimplex.worstVertex().value_)
101 {
102 nelderMeadSimplex.swapOutWorst(contractionPoint);
103 continue;
104 }
105
106 // If neither expansion nor contraction of the worst point give a
107 // good result shrink the whole simplex
108 nelderMeadSimplex.shrinkSimplexPointsExceptBest(functionToMinimize);
109 }
110
111 return { nelderMeadSimplex.bestVertex().coordinate_, nelderMeadSimplex.bestVertex().value_ };
112 }
113
114 } // namespace gmx
115