MNE-CPP  0.1.9
A Framework for Electrophysiology
simplex_algorithm.h
Go to the documentation of this file.
1 //=============================================================================================================
36 #ifndef SIMPLEXALGORITHM_H
37 #define SIMPLEXALGORITHM_H
38 
39 //=============================================================================================================
40 // INCLUDES
41 //=============================================================================================================
42 
43 #include "utils_global.h"
44 
45 //=============================================================================================================
46 // EIGEN INCLUDES
47 //=============================================================================================================
48 
49 #include <Eigen/Core>
50 
51 #define ALPHA 1.0
52 #define BETA 0.5
53 #define GAMMA 2.0
54 
55 //=============================================================================================================
56 // DEFINE NAMESPACE MNELIB
57 //=============================================================================================================
58 
59 namespace UTILSLIB
60 {
61 
62 //=============================================================================================================
69 {
70 
71 protected:
72  //=========================================================================================================
77 
78 public:
79  //=========================================================================================================
99  template <typename T>
100  static bool simplex_minimize(Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& p,
101  Eigen::Matrix<T,Eigen::Dynamic, 1>& y,
102  T ftol,
103  T (*func)(const Eigen::Matrix<T,Eigen::Dynamic, 1>& x, const void *user_data),
104  const void *user_data,
105  int max_eval,
106  int &neval,
107  int report,
108  bool (*report_func)(int loop, const Eigen::Matrix<T,Eigen::Dynamic, 1>& fitpar, double fval));
109 
110 private:
111 
112  template <typename T>
113  static T tryit(Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> &p,
114  Eigen::Matrix<T,Eigen::Dynamic, 1> &y,
115  Eigen::Matrix<T,Eigen::Dynamic, 1> &psum,
116  T (*func)( const Eigen::Matrix<T,Eigen::Dynamic, 1> &x,const void *user_data), /* The function to be evaluated */
117  const void *user_data, /* Data to be passed to the above function in each evaluation */
118  int ihi,
119  int &neval,
120  T fac);
121 };
122 
123 //=============================================================================================================
124 // DEFINE MEMBER METHODS
125 //=============================================================================================================
126 
127 template <typename T>
128 bool SimplexAlgorithm::simplex_minimize(Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& p,
129  Eigen::Matrix<T,Eigen::Dynamic, 1>& y,
130  T ftol,
131  T (*func)(const Eigen::Matrix<T,Eigen::Dynamic, 1>& x, const void *user_data),
132  const void *user_data,
133  int max_eval,
134  int &neval,
135  int report,
136  bool (*report_func)(int loop, const Eigen::Matrix<T,Eigen::Dynamic, 1>& fitpar, double fval))
137 {
138  int ndim = p.cols(); /* Number of variables */
139  int i,ilo,ihi,inhi;
140  int mpts = ndim+1;
141  T ytry,ysave,rtol;
142  Eigen::Matrix<T,Eigen::Dynamic, 1> psum(ndim);
143  bool result = true;
144  int count = 0;
145  int loop = 1;
146 
147  neval = 0;
148  psum = p.colwise().sum();
149 
150  if (report_func != NULL && report > 0) {
151  report_func(0,static_cast< Eigen::Matrix<T,Eigen::Dynamic, 1> >(p.row(0)),-1.0);
152  }
153 
154  for (;;count++,loop++) {
155  ilo = 1;
156  ihi = y[1]>y[2] ? (inhi = 2,1) : (inhi = 1,2);
157  for (i = 0; i < mpts; i++) {
158  if (y[i] < y[ilo])
159  ilo = i;
160  if (y[i] > y[ihi]) {
161  inhi = ihi;
162  ihi = i;
163  } else if (y[i] > y[inhi])
164  if (i != ihi)
165  inhi = i;
166  }
167  rtol = 2.0*std::fabs(y[ihi]-y[ilo])/(std::fabs(y[ihi])+std::fabs(y[ilo]));
168  /*
169  * Report that we are proceeding...
170  */
171  if (count == report && report_func != NULL) {
172  if (!report_func(loop,static_cast< Eigen::Matrix<T,Eigen::Dynamic, 1> >(p.row(ilo)),y[ilo])) {
173  qCritical("Interation interrupted.");
174  result = false;
175  break;
176  }
177  count = 0;
178  }
179  if (rtol < ftol) break;
180  if (neval >= max_eval) {
181  qCritical("Maximum number of evaluations exceeded.");
182  result = false;
183  break;
184  }
185  ytry = tryit<T>(p,y,psum,func,user_data,ihi,neval,-ALPHA);
186  if (ytry <= y[ilo])
187  tryit<T>(p,y,psum,func,user_data,ihi,neval,GAMMA);
188  else if (ytry >= y[inhi]) {
189  ysave = y[ihi];
190  ytry = tryit<T>(p,y,psum,func,user_data,ihi,neval,BETA);
191  if (ytry >= ysave) {
192  for (i = 0; i < mpts; i++) {
193  if (i != ilo) {
194  psum = 0.5 * ( p.row(i) + p.row(ilo) );
195  p.row(i) = psum;
196  y[i] = (*func)(psum,user_data);
197  }
198  }
199  neval += ndim;
200  psum = p.colwise().sum();
201  }
202  }
203  }
204 
205  return result;
206 }
207 
208 //=============================================================================================================
209 
210 template <typename T>
211 T SimplexAlgorithm::tryit(Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> &p,
212  Eigen::Matrix<T,Eigen::Dynamic, 1> &y,
213  Eigen::Matrix<T,Eigen::Dynamic, 1> &psum,
214  T (*func)( const Eigen::Matrix<T,Eigen::Dynamic, 1> &x,const void *user_data), /* The function to be evaluated */
215  const void *user_data, /* Data to be passed to the above function in each evaluation */
216  int ihi,
217  int &neval,
218  T fac)
219 {
220  int ndim = p.cols();
221  T fac1,fac2,ytry;
222 
223  Eigen::Matrix<T,Eigen::Dynamic, 1> ptry(ndim);
224 
225  fac1 = (1.0-fac)/ndim;
226  fac2 = fac1-fac;
227 
228  ptry = psum * fac1 - p.row(ihi).transpose() * fac2;
229 
230  ytry = (*func)(ptry,user_data);
231  ++neval;
232 
233  if (ytry < y[ihi]) {
234  y[ihi] = ytry;
235 
236  psum += ptry - p.row(ihi).transpose();
237  p.row(ihi) = ptry;
238  }
239 
240  return ytry;
241 }
242 } //NAMESPACE
243 
244 #undef ALPHA
245 #undef BETA
246 #undef GAMMA
247 
248 #endif // SIMPLEXALGORITHM_H
utils library export/import macros.
static bool simplex_minimize(Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &p, Eigen::Matrix< T, Eigen::Dynamic, 1 > &y, T ftol, T(*func)(const Eigen::Matrix< T, Eigen::Dynamic, 1 > &x, const void *user_data), const void *user_data, int max_eval, int &neval, int report, bool(*report_func)(int loop, const Eigen::Matrix< T, Eigen::Dynamic, 1 > &fitpar, double fval))
#define UTILSSHARED_EXPORT
Definition: utils_global.h:58