71 int MinimizerSimplex::mne_simplex_minimize(MatrixXf p,
75 float (*func)(
const VectorXf &x,
82 int (*report_func)(
int loop,
83 const VectorXf &fitpar,
89 float ytry,ysave,sum,rtol;
96 for (j = 0; j < ndim; j++) {
97 for (i = 0,sum = 0.0; i<mpts; i++)
101 if (report_func != NULL && report > 0)
102 (void)report_func(0,static_cast<VectorXf>(p.row(0)),ndim,-1.0);
104 for (;;count++,loop++) {
106 ihi = y[1]>y[2] ? (inhi = 2,1) : (inhi = 1,2);
107 for (i = 0; i < mpts; i++) {
113 }
else if (y[i] > y[inhi])
117 rtol = 2.0*fabs(y[ihi]-y[ilo])/(fabs(y[ihi])+fabs(y[ilo]));
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.";
131 if (neval >= max_eval) {
132 std::cout<<
"Maximum number of evaluations exceeded.";
136 ytry = tryit(p,y,psum,ndim,func,user_data,ihi,neval,-ALPHA);
138 ytry = tryit(p,y,psum,ndim,func,user_data,ihi,neval,GAMMA);
139 else if (ytry >= y[inhi]) {
141 ytry = tryit(p,y,psum,ndim,func,user_data,ihi,neval,BETA);
143 for (i = 0; i < mpts; i++) {
145 for (j = 0; j < ndim; j++) {
146 psum[j] = 0.5*(p(i,j)+p(ilo,j));
149 y[i] = (*func)(psum,ndim,user_data);
153 for (j = 0; j < ndim; j++) {
154 for (i = 0,sum = 0.0; i < mpts; i++)
168 float MinimizerSimplex::tryit(MatrixXf p,
172 float (*func)(
const VectorXf &x,
int npar,
void *user_data),
179 float fac1,fac2,ytry;
182 fac1 = (1.0-fac)/ndim;
185 for (j = 0; j < ndim; j++)
186 ptry[j] = psum[j]*fac1-p(ihi,j)*fac2;
188 ytry = (*func)(ptry,ndim,user_data);
194 for (j = 0; j < ndim; j++) {
195 psum[j] += ptry[j]-p(ihi,j);
MinimizerSimplex class declaration.