17 #if defined(DEBUG) || defined(WARNINGMSG) 28 class GradientCalculator;
36 std::cout <<
"Running Simplex with maxfcn = " << maxfcn <<
" minedm = " << minedm << std::endl;
43 unsigned int n = x.
size();
44 double wg = 1./double(n);
45 double alpha = 1.,
beta = 0.5,
gamma = 2., rhomin = 4., rhomax = 8.;
46 double rho1 = 1. + alpha;
49 double rho2 = 1. + alpha*
gamma;
52 std::vector<std::pair<double, MnAlgebraicVector> > simpl; simpl.reserve(n+1);
53 simpl.push_back(std::pair<double, MnAlgebraicVector>(seed.
Fval(),
x));
55 unsigned int jl = 0, jh = 0;
56 double amin = seed.
Fval(), aming = seed.
Fval();
58 for(
unsigned int i = 0; i <
n; i++) {
60 if(step(i) < dmin) step(i) = dmin;
71 simpl.push_back(std::pair<double, MnAlgebraicVector>(tmp, x));
77 std::cout <<
"simplex initial parameters - min " << jl <<
" " << amin <<
" max " << jh <<
" " << aming << std::endl;
78 for (
unsigned int i = 0; i < simplex.
Simplex().size(); ++i)
79 std::cout <<
" i = " << i <<
" x = " << simplex(i).second <<
" fval(x) = " << simplex(i).first << std::endl;
82 double edmPrev = simplex.
Edm();
87 amin = simplex(jl).first;
88 edmPrev = simplex.
Edm();
91 std::cout <<
"\n\nsimplex iteration: edm = " << simplex.
Edm()
92 <<
"\n--> Min Param is " << jl <<
" pmin " << simplex(jl).second <<
" f(pmin) " << amin
93 <<
"\n--> Max param is " << jh <<
" " << simplex(jh).first << std::endl;
107 for(
unsigned int i = 0; i < n+1; i++) {
108 if(i == jh)
continue;
109 pbar += (wg*simplex(i).second);
113 double ystar = mfcn(pstar);
116 std::cout <<
" pbar = " << pbar << std::endl;
117 std::cout <<
" pstar = " << pstar <<
" f(pstar) = " << ystar << std::endl;
121 if(ystar < simplex(jh).first) {
122 simplex.
Update(ystar, pstar);
123 if(jh != simplex.
Jh())
continue;
126 double ystst = mfcn(pstst);
128 std::cout <<
"Reduced simplex pstst = " << pstst <<
" f(pstst) = " << ystst << std::endl;
130 if(ystst > simplex(jh).first)
break;
131 simplex.
Update(ystst, pstst);
136 double ystst = mfcn(pstst);
138 std::cout <<
" pstst = " << pstst <<
" f(pstst) = " << ystst << std::endl;
141 double y1 = (ystar - simplex(jh).first)*rho2;
142 double y2 = (ystst - simplex(jh).first)*rho1;
143 double rho = 0.5*(rho2*y1 - rho1*y2)/(y1 - y2);
145 if(ystst < simplex(jl).first) simplex.
Update(ystst, pstst);
146 else simplex.
Update(ystar, pstar);
149 if(rho > rhomax) rho = rhomax;
151 double yrho = mfcn(prho);
153 std::cout <<
" prho = " << prho <<
" f(prho) = " << yrho << std::endl;
155 if(yrho < simplex(jl).first && yrho < ystst) {
156 simplex.
Update(yrho, prho);
159 if(ystst < simplex(jl).first) {
160 simplex.
Update(ystst, pstst);
163 if(yrho > simplex(jl).first) {
164 if(ystst < simplex(jl).first) simplex.
Update(ystst, pstst);
165 else simplex.
Update(ystar, pstar);
168 if(ystar > simplex(jh).first) {
169 pstst =
beta*simplex(jh).second + (1. -
beta)*pbar;
171 if(ystst > simplex(jh).first)
break;
172 simplex.
Update(ystst, pstst);
175 std::cout <<
"End loop : edm = " << simplex.
Edm() <<
" pstst = " << pstst <<
" f(pstst) = " << ystst << std::endl;
177 }
while( (simplex.
Edm() > minedm || edmPrev > minedm ) && mfcn.
NumOfCalls() < maxfcn);
181 amin = simplex(jl).first;
184 for(
unsigned int i = 0; i < n+1; i++) {
185 if(i == jh)
continue;
186 pbar += (wg*simplex(i).second);
188 double ybar = mfcn(pbar);
189 if(ybar < amin) simplex.
Update(ybar, pbar);
191 pbar = simplex(jl).second;
192 ybar = simplex(jl).first;
200 std::cout <<
"End simplex " << simplex.
Edm() <<
" pbar = " << pbar <<
" f(p) = " << ybar << std::endl;
212 MN_INFO_MSG(
"Simplex did not converge, #fcn calls exhausted.");
216 if(simplex.
Edm() > minedm) {
218 MN_INFO_MSG(
"Simplex did not converge, edm > minedm.");
const FunctionGradient & Gradient() const
This namespace contains pre-defined functions to be used in conjuction with TExecutor::Map and TExecu...
const MnAlgebraicVector & Vec() const
MinimumSeed contains the starting values for the minimization produced by the SeedGenerator.
const MnMachinePrecision & Precision() const
determines the relative floating point arithmetic precision.
double beta(double x, double y)
Calculates the beta function.
unsigned int NumOfCalls() const
const MinimumParameters & Parameters() const
class holding the full result of the minimization; both internal and external (MnUserParameterState) ...
static void PrintState(std::ostream &os, const MinimumState &state, const char *msg, int iter=-1)
Wrapper class to FCNBase interface used internally by Minuit.
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
void Update(double, const MnAlgebraicVector &)
unsigned int size() const
const MnAlgebraicVector & Gstep() const
class describing the simplex set of points (f(x), x ) which evolve during the minimization iteration ...
const std::vector< std::pair< double, MnAlgebraicVector > > & Simplex() const
virtual FunctionMinimum Minimum(const MnFcn &, const GradientCalculator &, const MinimumSeed &, const MnStrategy &, unsigned int, double) const
double Eps2() const
eps2 returns 2*sqrt(eps)
void TraceIteration(int iter, const MinimumState &state) const
MnAlgebraicVector Dirin() const
MinimumState keeps the information (position, Gradient, 2nd deriv, etc) after one minimization step (...
API class for defining three levels of strategies: low (0), medium (1), high (>=2); acts on: Migrad (...
interface class for gradient calculators