MNE-CPP  beta 1.0
minimizersimplex.cpp
1 //=============================================================================================================
37 //*************************************************************************************************************
38 //=============================================================================================================
39 // INCLUDES
40 //=============================================================================================================
41 
42 #include "minimizersimplex.h"
43 
44 
45 //*************************************************************************************************************
46 //=============================================================================================================
47 // Qt INCLUDES
48 //=============================================================================================================
49 
50 
51 //*************************************************************************************************************
52 //=============================================================================================================
53 // USED NAMESPACES
54 //=============================================================================================================
55 
56 using namespace UTILSLIB;
57 
58 
59 //*************************************************************************************************************
60 //=============================================================================================================
61 // DEFINE MEMBER METHODS
62 //=============================================================================================================
63 
65 {
66 }
67 
68 
69 //*************************************************************************************************************
70 
71 int MinimizerSimplex::mne_simplex_minimize(MatrixXf p,
72  VectorXf y,
73  int ndim,
74  float ftol,
75  float (*func)(const VectorXf &x,
76  int npar,
77  void *user_data),
78  void *user_data,
79  int max_eval,
80  int &neval,
81  int report,
82  int (*report_func)(int loop,
83  const VectorXf &fitpar,
84  int npar,
85  double fval))
86 {
87  int i,j,ilo,ihi,inhi;
88  int mpts = ndim+1;
89  float ytry,ysave,sum,rtol;
90  VectorXf psum(ndim);
91  int result = 0;
92  int count = 0;
93  int loop = 1;
94 
95  neval = 0;
96  for (j = 0; j < ndim; j++) {
97  for (i = 0,sum = 0.0; i<mpts; i++)
98  sum += p(i,j);
99  psum[j] = sum;
100  }
101  if (report_func != NULL && report > 0)
102  (void)report_func(0,static_cast<VectorXf>(p.row(0)),ndim,-1.0);
103 
104  for (;;count++,loop++) {
105  ilo = 1;
106  ihi = y[1]>y[2] ? (inhi = 2,1) : (inhi = 1,2);
107  for (i = 0; i < mpts; i++) {
108  if (y[i] < y[ilo])
109  ilo = i;
110  if (y[i] > y[ihi]) {
111  inhi = ihi;
112  ihi = i;
113  } else if (y[i] > y[inhi])
114  if (i != ihi)
115  inhi = i;
116  }
117  rtol = 2.0*fabs(y[ihi]-y[ilo])/(fabs(y[ihi])+fabs(y[ilo]));
118  /*
119  * Report that we are proceeding...
120  */
121  if (count == report && report_func != NULL) {
122  if (report_func(loop,static_cast<VectorXf>(p.row(ilo)),ndim,y[ilo])) {
123  std::cout<<"Interation interrupted.";
124  result = -1;
125  break;
126  }
127  count = 0;
128  }
129  if (rtol < ftol)
130  break;
131  if (neval >= max_eval) {
132  std::cout<<"Maximum number of evaluations exceeded.";
133  result = -1;
134  break;
135  }
136  ytry = tryit(p,y,psum,ndim,func,user_data,ihi,neval,-ALPHA);
137  if (ytry <= y[ilo])
138  ytry = tryit(p,y,psum,ndim,func,user_data,ihi,neval,GAMMA);
139  else if (ytry >= y[inhi]) {
140  ysave = y[ihi];
141  ytry = tryit(p,y,psum,ndim,func,user_data,ihi,neval,BETA);
142  if (ytry >= ysave) {
143  for (i = 0; i < mpts; i++) {
144  if (i != ilo) {
145  for (j = 0; j < ndim; j++) {
146  psum[j] = 0.5*(p(i,j)+p(ilo,j));
147  p(i,j) = psum[j];
148  }
149  y[i] = (*func)(psum,ndim,user_data);
150  }
151  }
152  neval += ndim;
153  for (j = 0; j < ndim; j++) {
154  for (i = 0,sum = 0.0; i < mpts; i++)
155  sum += p(i,j);
156  psum[j] = sum;
157  }
158  }
159  }
160  }
161 
162  return (result);
163 }
164 
165 
166 //*************************************************************************************************************
167 
168 float MinimizerSimplex::tryit(MatrixXf p,
169  VectorXf y,
170  VectorXf psum,
171  int ndim,
172  float (*func)(const VectorXf &x,int npar,void *user_data),
173  void *user_data,
174  int ihi,
175  int &neval,
176  float fac)
177 {
178  int j;
179  float fac1,fac2,ytry;
180  VectorXf ptry(ndim);
181 
182  fac1 = (1.0-fac)/ndim;
183  fac2 = fac1-fac;
184 
185  for (j = 0; j < ndim; j++)
186  ptry[j] = psum[j]*fac1-p(ihi,j)*fac2;
187 
188  ytry = (*func)(ptry,ndim,user_data);
189  ++(neval);
190 
191  if (ytry < y[ihi]) {
192  y[ihi] = ytry;
193 
194  for (j = 0; j < ndim; j++) {
195  psum[j] += ptry[j]-p(ihi,j);
196  p(ihi,j) = ptry[j];
197  }
198  }
199 
200  return ytry;
201 }
202 
203 #undef ALPHA
204 #undef BETA
205 #undef GAMMA
206 
MinimizerSimplex class declaration.