LOOS  v2.3.2
Simplex.hpp
1 /*
2  Nelder-Meade Simplex Optimizer
3  (c) 2009 Tod D. Romo, Grossfield Lab, URMC
4 
5  Based loosely on the NRC implementation (Numerical Recipes in C (1996):411)
6 */
7 
8 #if !defined(LOOS_SIMPLEX_HPP)
9 #define LOOS_SIMPLEX_HPP
10 
11 #include <iostream>
12 #include <vector>
13 
14 
16 
57 template<typename T = double>
58 class Simplex {
59  typedef std::vector< std::vector<T> > VVectors;
60 
61 
62  template<class C>
63  T modify(T factor, C& ftor) {
64  int j;
65  T f1, f2, val;
66 
67  f1 = (1.0 - factor) / ndim;
68  f2 = f1 - factor;
69 
70  // Compute the new simplex vertex
71 
72  for (j=0; j<ndim; ++j)
73  trial[j] = simpsum[j] * f1 - simplex[worst][j] * f2;
74 
75  val = ftor(trial);
76 
77  // Use this vertex in the simplex if it's better...
78 
79  if (val < values[worst]) {
80  values[worst] = val;
81  for (j=0; j<ndim; j++) {
82  simpsum[j] += trial[j] - simplex[worst][j];
83  simplex[worst][j] = trial[j];
84  }
85  }
86 
87  return(val);
88  }
89 
90 
91 
92  // The core of the simplex optimizer...
93  template<class C>
94  void core(C& ftor) {
95  int i, next_worst, j, mpts;
96  T sum, saved, val, num, den;
97 
98  mpts = ndim + 1;
99 
100  for (j=0; j<ndim; j++) {
101  for (sum = 0.0, i = 0; i<mpts; i++)
102  sum += simplex[i][j];
103  simpsum[j] = sum;
104  }
105 
106  n_evals = 0;
107 
108  while (n_evals <= maxiters) {
109  best = 1;
110  if (values[0] > values[1]) {
111  worst = 0;
112  next_worst = 1;
113  } else {
114  worst = 1;
115  next_worst = 0;
116  }
117 
118 
119  // Find candidate vertices in the simplex...
120 
121  for (i=0; i<mpts; i++) {
122  if (values[i] <= values[best])
123  best = i;
124  if (values[i] > values[worst]) {
125  next_worst = worst;
126  worst = i;
127  } else if (values[i] > values[next_worst] && i != worst)
128  next_worst = i;
129  }
130 
131  // Check for convergence...
132  // Some conditions may have the numerator and denominator equal (or both 0)
133  // which can cause problems below, hence the special test.
134 
135  num = fabs(values[worst] - values[best]);
136  den = fabs(values[worst]) + fabs(values[best]);
137  rtol = 2.0 * num / den;
138  if (rtol < tol || num == den)
139  return;
140 
141  // Now try reflecting, contracting, expanding the simplex...
142 
143  n_evals += 2;
144  val = modify(-1.0, ftor);
145  if (val <= values[best])
146  val = modify(2.0, ftor);
147  else if (val >= values[next_worst]) {
148 
149  saved = values[worst];
150  val = modify(0.5, ftor);
151 
152  if (val >= saved) {
153  for (i=0; i<mpts; i++) {
154  if (i != best) {
155  for (j=0; j<ndim; j++)
156  simplex[i][j] = simpsum[j] = 0.5 * (simplex[i][j] + simplex[best][j]);
157  values[i] = ftor(simpsum);
158  }
159  }
160  n_evals += ndim;
161 
162  for (j=0; j<ndim; j++) {
163  for (sum=0.0, i=0; i<mpts; i++)
164  sum += simplex[i][j];
165  simpsum[j] = sum;
166  }
167  }
168 
169  } else
170  --n_evals;
171 
172  }
173 
174  }
175 
176 
177  void allocateSpace(const int n) {
178  q = std::vector<T>(n+1, 0.0);
179  qq = std::vector<T>(n+1, 0.0);
180  simpsum = std::vector<T>(n, 0.0);
181  values = std::vector<T>(n+1, 0.0);
182  trial = std::vector<T>(n, 0.0);
183  simplex.clear();
184  for (int i=0; i<n+1; ++i)
185  simplex.push_back(std::vector<T>(n, 0.0));
186  }
187 
188 
189 
190 public:
191 
192  Simplex(const int n) : tol(1e-3), ndim(n), maxiters(2000), best(-1), worst(-1), n_evals(0) {
193  allocateSpace(n);
194  }
195 
197  void dim(const int n) { ndim = n; allocateSpace(n); }
198 
200  void seedLengths(const std::vector<T>& seeds) { characteristics = seeds; }
201 
203  void tolerance(const double d) { tol = d; }
204 
206  void maximumIterations(const int n) { maxiters = n; }
207 
208  int maximumIterations() const { return maxiters; }
209 
211  int numberOfIterations() const { return n_evals; }
212 
214  std::vector<T> finalParameters(void) const {
215  if (best < 0)
216  throw(std::logic_error("Simplex has not been optimized"));
217  return(simplex[best]);
218  }
219 
221  T finalValue(void) const { return(values[best]); }
222 
224  template<class C>
225  std::vector<T> optimize(std::vector<T>& f, C& ftor) {
226  int i, j, n;
227 
228  if (characteristics.size() != f.size())
229  throw(std::logic_error("Invalid seed"));
230 
231  n = ndim + 1;
232 
233  // Initial simplex is based on Nelder-Meade's construction...
234 
235  for (i=0; i<ndim; i++) {
236  q[i] = characteristics[i] * ( (sqrtf(n) + ndim) / (n * sqrt(2.0)) );
237  qq[i] = characteristics[i] * ( (sqrtf(n) - 1.0) / (n * sqrt(2.0)) );
238  }
239 
240  for (j=0; j<n; j++)
241  for (i=0; i<ndim; i++)
242  if (j == i+1)
243  simplex[j][i] = f[i] + qq[i];
244  else
245  simplex[j][i] = f[i] + q[i];
246 
247  for (j=0; j<n; ++j) {
248  values[j] = ftor(simplex[j]);
249  }
250 
251  core(ftor);
252  return(simplex[best]);
253  }
254 
255 private:
256  double tol;
257  int ndim, maxiters, best, worst;
258  double rtol;
259  int n_evals;
260 
261  std::vector<T> characteristics;
262  std::vector<T> simpsum;
263  std::vector<T> values;
264  std::vector<T> q;
265  std::vector<T> qq;
266  std::vector<T> trial;
267  VVectors simplex;
268 
269 };
270 
271 
272 
273 
274 #endif
void dim(const int n)
Set the number of dimensions.
Definition: Simplex.hpp:197
void maximumIterations(const int n)
Limit on the number of function evaluations to perform.
Definition: Simplex.hpp:206
std::vector< T > optimize(std::vector< T > &f, C &ftor)
Optimize via a functor.
Definition: Simplex.hpp:225
T finalValue(void) const
Final (best) value.
Definition: Simplex.hpp:221
void seedLengths(const std::vector< T > &seeds)
Initial guess.
Definition: Simplex.hpp:200
int numberOfIterations() const
Return the number of iterations that it actually went through.
Definition: Simplex.hpp:211
void tolerance(const double d)
Convergence criterion.
Definition: Simplex.hpp:203
std::vector< T > finalParameters(void) const
Retrieve the final (best fit) parameters.
Definition: Simplex.hpp:214
Nelder-Meade Simplex Optimizer (based loosely on the NRC (1996) implementation)
Definition: Simplex.hpp:58