LOOS  v2.3.2
mc_fit.hpp
1 #if !defined(LOOS_MC_FIT_HPP)
2 #define LOOS_MC_FIT_HPP
3 
4 #include <iostream>
5 #include <vector>
6 #include <ctime>
7 #include <boost/random/linear_congruential.hpp>
8 #include <boost/random/uniform_real.hpp>
9 #include <boost/random/variate_generator.hpp>
10 
11 // @cond TOOLS_INTERNAL
12 
13 typedef boost::minstd_rand rand_num;
14 
15 struct ConstantAcceptor {
16  ConstantAcceptor() : val(0.25) { }
17  ConstantAcceptor(double d) : val(d) { }
18  double operator()(const uint iter) { return(val); }
19 
20  double val;
21 };
22 
23 
24 struct ExponentialAcceptor {
25  ExponentialAcceptor() : k(1.0) { }
26  ExponentialAcceptor(const double scale) : k(scale) { }
27  double operator()(const uint iter) { return(exp(-k * iter)); }
28 
29  double k;
30 };
31 
32 
33 template<typename T = double>
34 class mcoptimo {
35  typedef std::vector< std::vector<T> > VVectors;
36 
37 public:
38  // void reseed(const mcoptimo& previous, mcoptimo &current){
39  // current.setParams(randomize(previous.getParams()));
40  // // current.setParams() = randomize(previous.getParams());
41  // }
42 
43 
44  void setParams(const vector<T>& v) {
45  myparams = v;
46  }
47 
48  vector<T> getParams() const { return(myparams); }
49 
50 
51  // vector<T>& setParams() { return(myparams); }
52 
53 
54  //private:
55  // vector<T> myparams;
56 
57 
58 // void accept(&previous, const &current){
59 // if (current.getCov() > previous.getCov()){ //accept
60 // previous = current;
61 // // }else if (current.getCov() > e^(Energy/T_star)){ //accept
62 // // previous = current;
63 // } // else{ //reject
64 // // continue;
65 // // }
66 // }
67 
68  // void optimum(const &current, &best){
69  // if (current.getCov() > best.getCov()) best = current;
70  // }
71 
72  // Return a perturbed set of parameters
73  vector<T> randomize(const vector<T>& current, const vector<T>& sizes) {
74  //use the BOOST random generator per alan's suggestion
75  //using type pseudo-random number generator, based on the previous k's
76  //rand_num generator(time(0)); //<---Talk to Tod about other seeds!!!
77  // base_generator_type& generator = rng_singleton();
78 
79  vector<T> newParams;
80  for (int i = 0; i < previous.getParams->paramSize(); ++i){//i want the size of the vector myparams
81  //for the 'previous' instance of mcoptimize
82  T springval = 0;//add stuff here to grab that a particular param
83 
84  // Set range/bounds to the size here...
85  // T lowerbound = springval / 2;
86  // T upperbound = 3 * springval / 2;
87  // boost::uniform_real<> uni_dist(lowerbound , upperbound);
88  boost::uniform_real<> uni_dist(-1, 1);
89  boost::variate_generator<rand_num&, boost::uniform_real<> > uni(generator, uni_dist);
90  double scaled_random = (uni() + uni()) * current[i];
91  newParams.push_back(scaled_random + current[i]);
92  }
93  return(newParams);
94  }
95 
96  // template<class C>
97  T randomize(uint iter){
98  base_generator_type& single_random = rng_singleton();
99  T upperbound = acceptance(iter);
100  boost::uniform_real<> uni_dist(0, upperbound*2);
101  boost::variate_generator<rand_num&, boost::uniform_real<> > uni(single_random, uni_dist);
102  double my_random = uni();
103  return(my_random);
104  }
105 
106 // template<class C>
107 // T acceptance(uint iter, const uint stepSize, const uint absTemp){
108 // //maybe use scope to grab stepSize and absTemp
109 // //but what scope do they belong in??
110 
111 // //this should not be hard wired.
112 // //we want another ftor that pts
113 // //to an acceptance function
114 // T cutoff = absTemp - (stepSize * iter);
115 // return(cutoff);
116 // }
117 
118 
119  template<class C, class Acceptor = ConstantAcceptor>
120  vector<T> takeStep(const vector<T>& current, const vector<T>& sizes, C& ftor, Acceptor& acc, uint iter) {
121 
122  vector<T> newStep = randomize(current, sizes);
123  T prev = ftor(current);
124  T val = ftor(newStep);
125 
126  if (val < prev){
127  return(newStep);
128  }elseif (randomize(iter) < acc(iter)){//define both of these!!!
129  return(newStep);
130  }
131  return(current);
132  }
133 
134 
135 
136  // vector<double> takeStep(const vector<double>& current, const vector<double>& sizes, FitAggregator& ftor)
137  template<class C, class A = ConstantAcceptor>
138  vector<T> optimize(const vector<T>& current, C& ftor, A& ator) {
139 
140  vector<T> params(current);//<<----do i need this??
141  vector<T> sizes(initial_sizes);
142  uint iter = 0;
143 
144  while (!converged) {
145  vector<T> params = takeStep(params, sizes, ftor, ator, iter);
146  ++iter;
147 
148  return(params);
149  }
150 
151  void setSizes(const vector<T>& s) {
152  initial_sizes = s;
153  }
154 
155 
156 private:
157  //double cov;
158  //params K; //should this be a springFactory instance???
159  vector<T> initial_sizes;
160  vector<T> myparams;
161 
162 
164 
165 
166 
167  // The core of the optimizer...
168  template<class C>
169  void core(C& ftor) {
170  int i, next_worst, j, mpts;
171  T best, current, previous;//sum, saved, val, num, den;
172 
173  mpts = ndim + 1;
174 
175  for (j=0; j<ndim; j++) {
176  for (sum = 0.0, i = 0; i<mpts; i++)
177 
178  sum += simplex[i][j];
179  simpsum[j] = sum;
180  }
181 
182  int n_evals = 0;
183 
184  while (n_evals <= maxiters) {
185  best = 1;
186  if (values[0] > values[1]) {
187  worst = 0;
188  next_worst = 1;
189  } else {
190  worst = 1;
191  next_worst = 0;
192  }
193 
194 
195  // Find candidate vertices in the simplex...
196 
197  for (i=0; i<mpts; i++) {
198  if (values[i] <= values[best])
199  best = i;
200  if (values[i] > values[worst]) {
201  next_worst = worst;
202  worst = i;
203  } else if (values[i] > values[next_worst] && i != worst)
204  next_worst = i;
205  }
206 
207  // Check for convergence...
208  // Some conditions may have the numerator and denominator equal (or both 0)
209  // which can cause problems below, hence the special test.
210 
211  num = fabs(values[worst] - values[best]);
212  den = fabs(values[worst]) + fabs(values[best]);
213  rtol = 2.0 * num / den;
214  if (rtol < tol || num == den)
215  return;
216 
217  // Now try reflecting, contracting, expanding the simplex...
218 
219  n_evals += 2;
220  val = modify(-1.0, ftor);
221  if (val <= values[best])
222  val = modify(2.0, ftor);
223  else if (val >= values[next_worst]) {
224 
225  saved = values[worst];
226  val = modify(0.5, ftor);
227 
228  if (val >= saved) {
229  for (i=0; i<mpts; i++) {
230  if (i != best) {
231  for (j=0; j<ndim; j++)
232  simplex[i][j] = simpsum[j] = 0.5 * (simplex[i][j] + simplex[best][j]);
233  values[i] = ftor(simpsum);
234  }
235  }
236  n_evals += ndim;
237 
238  for (j=0; j<ndim; j++) {
239  for (sum=0.0, i=0; i<mpts; i++)
240  sum += simplex[i][j];
241  simpsum[j] = sum;
242  }
243  }
244 
245  } else
246  --n_evals;
247 
248  }
249 
250  }
251 
252 
253  void allocateSpace(const int n) {
254  q = std::vector<T>(n+1, 0.0);
255  qq = std::vector<T>(n+1, 0.0);
256  simpsum = std::vector<T>(n, 0.0);
257  values = std::vector<T>(n+1, 0.0);
258  trial = std::vector<T>(n, 0.0);
259  simplex.clear();
260  for (int i=0; i<n+1; ++i)
261  simplex.push_back(std::vector<T>(n, 0.0));
262  }
263 
264 
265 
266 public:
267 
268  Simplex(const int n) : tol(1e-3), ndim(n), maxiters(2000), best(-1), worst(-1) {
269  allocateSpace(n);
270  }
271 
273  void dim(const int n) { ndim = n; allocateSpace(n); }
274 
276  void seedLengths(const std::vector<T>& seeds) { characteristics = seeds; }
277 
279  void tolerance(const double d) { tol = d; }
280 
282  void maximumIterations(const int n) { maxiters = n; }
283 
285  std::vector<T> finalParameters(void) const {
286  if (best < 0)
287  throw(std::logic_error("Simplex has not been optimized"));
288  return(simplex[best]);
289  }
290 
292  T finalValue(void) const { return(values[best]); }
293 
295  template<class C>
296  std::vector<T> optimize(std::vector<T>& f, C& ftor) {
297  int i, j, n;
298 
299  if (characteristics.size() != f.size())
300  throw(std::logic_error("Invalid seed"));
301 
302  n = ndim + 1;
303 
304  // Initial simplex is based on Nelder-Meade's construction...
305 
306  for (i=0; i<ndim; i++) {
307  q[i] = characteristics[i] * ( (sqrtf(n) + ndim) / (n * sqrt(2.0)) );
308  qq[i] = characteristics[i] * ( (sqrtf(n) - 1.0) / (n * sqrt(2.0)) );
309  }
310 
311  for (j=0; j<n; j++)
312  for (i=0; i<ndim; i++)
313  if (j == i+1)
314  simplex[j][i] = f[i] + qq[i];
315  else
316  simplex[j][i] = f[i] + q[i];
317 
318  for (j=0; j<n; ++j) {
319  values[j] = ftor(simplex[j]);
320  }
321 
322  core(ftor);
323  return(simplex[best]);
324  }
325 
326 private:
327  double tol;
328  int ndim, maxiters, best, worst;
329  double rtol;
330 
331  std::vector<T> characteristics;
332  std::vector<T> simpsum;
333  std::vector<T> values;
334  std::vector<T> q;
335  std::vector<T> qq;
336  std::vector<T> trial;
337  VVectors simplex;
338 
339 };
340 
341 
342 // @endcond
343 
344 #endif
base_generator_type & rng_singleton(void)
Suite-wide random number generator singleton.
Nelder-Meade Simplex Optimizer (based loosely on the NRC (1996) implementation)
Definition: Simplex.hpp:58