2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
37 * \brief Implements routines in neldermead.h .
39 * \author Christian Blau <blau@kth.se>
44 #include "neldermead.h"
51 #include "gromacs/math/utilities.h"
52 #include "gromacs/utility/arrayref.h"
53 #include "gromacs/utility/gmxassert.h"
61 /*! \brief Evaluate the linear combination of two vectors a and b.
63 * \param[in] alpha scaling factor for a
64 * \param[in] a vector to be scaled
65 * \param[in] beta scaling factor for b
66 * \param[in] b vector to be scaled
68 * \returns alpha * a + beta * b.
70 std::vector<real> linearCombination(real alpha, ArrayRef<const real> a, real beta, ArrayRef<const real> b)
72 GMX_ASSERT(a.size() == b.size(),
73 "Input vectors have to have the same size to evaluate their linear combination.");
74 std::vector<real> result(a.size());
75 std::transform(std::begin(a), std::end(a), std::begin(b), std::begin(result),
76 [alpha, beta](auto a, auto b) { return alpha * a + beta * b; });
81 * \brief The parameters for a Nelder-Mead optimisation.
83 struct NelderMeadParameters
85 //! Factor to evaluate the reflection point
87 //! Factor to evaluate the expansion point
89 //! Factor to evaluate the contraction point
91 //! Factor to evaluate the simplex shrinkage
95 constexpr NelderMeadParameters defaultNelderMeadParameters = { 1, 2, 0.5, 0.5 };
100 NelderMeadSimplex::NelderMeadSimplex(const std::function<real(ArrayRef<const real>)>& f,
101 ArrayRef<const real> initalGuess)
103 // initial simplex contains the initally guessed vertex
104 std::vector<real> initalVertex = copyOf(initalGuess);
105 simplex_.push_back({ initalVertex, f(initalVertex) });
106 // create the missing verticies by moving 0.05 or 0.0025 if null
107 // from the initial vertex dimension after dimension
108 for (auto& v : initalVertex)
110 const auto oldValue = v;
119 simplex_.push_back({ initalVertex, f(initalVertex) });
122 simplex_.sort([](const RealFunctionvalueAtCoordinate& lhs,
123 const RealFunctionvalueAtCoordinate& rhs) { return lhs.value_ < rhs.value_; });
124 updateCentroidAndReflectionPoint();
127 RealFunctionvalueAtCoordinate
128 NelderMeadSimplex::evaluateReflectionPoint(const std::function<real(ArrayRef<const real>)>& f) const
130 return { reflectionPointCoordinates_, f(reflectionPointCoordinates_) };
133 const RealFunctionvalueAtCoordinate& NelderMeadSimplex::bestVertex() const
135 return simplex_.front();
138 const RealFunctionvalueAtCoordinate& NelderMeadSimplex::worstVertex() const
140 return simplex_.back();
143 real NelderMeadSimplex::secondWorstValue() const
145 // go backwards one step from the end of the sorted simplex list
146 // and look at the vertex value
147 return std::next(std::rbegin(simplex_))->value_;
150 RealFunctionvalueAtCoordinate
151 NelderMeadSimplex::evaluateExpansionPoint(const std::function<real(ArrayRef<const real>)>& f) const
153 const std::vector<real> expansionPointCoordinate =
154 linearCombination(1 - defaultNelderMeadParameters.gamma_, centroidWithoutWorstPoint_,
155 defaultNelderMeadParameters.gamma_, reflectionPointCoordinates_);
156 return { expansionPointCoordinate, f(expansionPointCoordinate) };
159 RealFunctionvalueAtCoordinate
160 NelderMeadSimplex::evaluateContractionPoint(const std::function<real(ArrayRef<const real>)>& f) const
162 std::vector<real> contractionPoint =
163 linearCombination(1 - defaultNelderMeadParameters.rho_, centroidWithoutWorstPoint_,
164 defaultNelderMeadParameters.rho_, worstVertex().coordinate_);
165 return { contractionPoint, f(contractionPoint) };
168 void NelderMeadSimplex::swapOutWorst(const RealFunctionvalueAtCoordinate& newVertex)
170 // drop the worst point - we know it's at the back of the simplex list, because
171 // we kept the list sorted
173 // find the point to insert the new vertex, so that the simplex vertices
174 // keep being sorted according to function value
175 const auto insertionPoint = std::lower_bound(
176 std::begin(simplex_), std::end(simplex_), newVertex.value_,
177 [](const RealFunctionvalueAtCoordinate& lhs, real value) { return lhs.value_ < value; });
178 simplex_.insert(insertionPoint, newVertex);
179 // now that the simplex has changed, it has a new centroid and reflection point
180 updateCentroidAndReflectionPoint();
183 void NelderMeadSimplex::shrinkSimplexPointsExceptBest(const std::function<real(ArrayRef<const real>)>& f)
185 std::vector<real> bestPointCoordinate = simplex_.front().coordinate_;
186 // skipping over the first simplex vertex, pull points closer to the best
188 std::transform(std::next(std::begin(simplex_)), std::end(simplex_), std::next(std::begin(simplex_)),
189 [bestPointCoordinate, f](const RealFunctionvalueAtCoordinate& d) -> RealFunctionvalueAtCoordinate {
190 const std::vector<real> shrinkPoint = linearCombination(
191 defaultNelderMeadParameters.sigma_, d.coordinate_,
192 1 - defaultNelderMeadParameters.sigma_, bestPointCoordinate);
193 return { shrinkPoint, f(shrinkPoint) };
196 simplex_.sort([](const RealFunctionvalueAtCoordinate& lhs,
197 const RealFunctionvalueAtCoordinate& rhs) { return lhs.value_ < rhs.value_; });
199 // now that the simplex has changed, it has a new centroid and reflection point
200 updateCentroidAndReflectionPoint();
203 real NelderMeadSimplex::orientedLength() const
206 const std::vector<real> firstSimplexVertexCoordinate = simplex_.front().coordinate_;
207 // find out which vertex coordinate has the largest distance to the first simplex vertex.
208 for (const auto& simplexVertex : simplex_)
210 const std::vector<real> differenceVector =
211 linearCombination(1, firstSimplexVertexCoordinate, -1, simplexVertex.coordinate_);
212 const real thisLength =
213 std::accumulate(std::begin(differenceVector), std::end(differenceVector), 0.,
214 [](real sum, real value) { return sum + value * value; });
215 result = std::max(result, thisLength);
220 void NelderMeadSimplex::updateCentroidAndReflectionPoint()
222 // intialize with first vertex, then add up all other vertex coordinates
224 centroidWithoutWorstPoint_ = std::accumulate(
225 std::next(std::begin(simplex_)), std::prev(std::end(simplex_)), simplex_.front().coordinate_,
226 [](std::vector<real> sum, const RealFunctionvalueAtCoordinate& x) {
227 std::transform(std::begin(sum), std::end(sum), std::begin(x.coordinate_),
228 std::begin(sum), std::plus<>());
232 // divide the summed up coordinates by N (the simplex has N+1 vertices)
233 std::transform(std::begin(centroidWithoutWorstPoint_), std::end(centroidWithoutWorstPoint_),
234 std::begin(centroidWithoutWorstPoint_),
235 [n = simplex_.size() - 1](const auto& x) { return x / n; });
237 // now, that we have evaluated the centroid, update the reflection points
238 reflectionPointCoordinates_ =
239 linearCombination(defaultNelderMeadParameters.alpha_ + 1, centroidWithoutWorstPoint_,
240 -1, worstVertex().coordinate_);