Logo ROOT   6.08/07
Reference Guide
testMinimizer.cxx
Go to the documentation of this file.
1 // test of minimization using new minimizer classes
2 
3 #include "Math/Minimizer.h"
4 #include "Math/Factory.h"
5 #include "Math/Functor.h"
6 
7 #include "TVirtualFitter.h"
8 
9 #include "Math/IFunction.h"
10 #include "Math/Util.h"
11 #include <cmath>
12 #include <cassert>
13 
14 #include <string>
15 #include <iostream>
16 
17 #include "TStopwatch.h"
18 #include "TMatrixD.h"
19 #include "TVectorD.h"
20 #include "TRandom3.h"
21 #include "TMath.h"
22 
23 //#define DEBUG
24 
25 int gNCall = 0;
26 int gNCall2 = 0;
27 int gNmin = 1000;
28 int gVerbose = 0;
29 bool useGradient = false;
30 
31 bool minos = true;
32 
33 double gAbsTolerance = 0.005;
34 
35 // Rosenbrok function to be minimize
36 
37 typedef void (*FCN)(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag);
38 
39 
40 
41 // ROSENBROCK function
42 ////////////////////////////////////////////////////////////////////////////////
43 
44 void RosenBrock(Int_t &, Double_t *, Double_t &f, Double_t *par, Int_t /*iflag*/)
45 {
46  gNCall++;
47  const Double_t x = par[0];
48  const Double_t y = par[1];
49  const Double_t tmp1 = y-x*x;
50  const Double_t tmp2 = 1-x;
51  f = 100*tmp1*tmp1+tmp2*tmp2;
52 }
53 
54 
55 
56 class RosenBrockFunction : public ROOT::Math::IMultiGenFunction {
57 
58 public :
59 
60  virtual ~RosenBrockFunction() {}
61 
62  unsigned int NDim() const { return 2; }
63 
65  return new RosenBrockFunction();
66  }
67 
68  const double * TrueMinimum() const {
69  fTrueMin[0] = 1;
70  fTrueMin[1] = 1;
71  return fTrueMin;
72  }
73 
74  private:
75 
76  inline double DoEval (const double * x) const {
77 #ifdef USE_FREE_FUNC
78  double f = 0;
79  int ierr = 0;
80  int i = 0;
81  RosenBrock(i,0,f,const_cast<double *>(x),ierr);
82  return f;
83 #else
84  gNCall++;
85  const Double_t xx = x[0];
86  const Double_t yy = x[1];
87  const Double_t tmp1 = yy-xx*xx;
88  const Double_t tmp2 = 1-xx;
89  return 100*tmp1*tmp1+tmp2*tmp2;
90 #endif
91  }
92 
93 
94  mutable double fTrueMin[2];
95 };
96 
97 
98 // TRIGONOMETRIC FLETCHER FUNCTION
99 
100 class TrigoFletcherFunction : public ROOT::Math::IMultiGradFunction {
101 
102 public :
103 
104 
105  TrigoFletcherFunction(unsigned int dim) : fDim(dim) {
106  double seed = 3;
107  A.ResizeTo(dim,dim);
108  B.ResizeTo(dim,dim);
109  x0.ResizeTo(dim);
110  sx0.ResizeTo(dim);
111  cx0.ResizeTo(dim);
112  sx.ResizeTo(dim);
113  cx.ResizeTo(dim);
114  v0.ResizeTo(dim);
115  v.ResizeTo(dim);
116  r.ResizeTo(dim);
117  A.Randomize(-100.,100,seed);
118  B.Randomize(-100.,100,seed);
119  for (unsigned int i = 0; i < dim; i++) {
120  for (unsigned int j = 0; j < dim; j++) {
121  A(i,j) = int(A(i,j));
122  B(i,j) = int(B(i,j));
123  }
124  }
125  x0.Randomize(-TMath::Pi(),TMath::Pi(),seed);
126  // calculate vector Ei
127  for (unsigned int i = 0; i < fDim ; ++i) {
128  cx0[i] = std::cos(x0[i]);
129  sx0[i] = std::sin(x0[i]);
130  }
131  v0 = A*sx0+B*cx0;
132  }
133 
134 
135  unsigned int NDim() const { return fDim; }
136 
137  ROOT::Math::IMultiGenFunction * Clone() const {
138  TrigoFletcherFunction * f = new TrigoFletcherFunction(*this);
139 // std::cerr <<"cannot clone this function" << std::endl;
140 // assert(0);
141  return f;
142  }
143 
144 
145  void StartPoints(double * x, double * s) {
146  TRandom3 rndm;
147  const double stepSize = 0.01;
148  const double deltaAmp = 0.1;
149  const double pi = TMath::Pi();
150  for (unsigned int i = 0; i < fDim; ++i) {
151  double delta = rndm.Uniform(-deltaAmp*pi,deltaAmp*pi);
152  x[i] = x0(i) + 0.1*delta;
153  if (x[i] <= - pi) x[i] += 2.*pi;
154  if (x[i] > pi) x[i] -= 2.*pi;
155  s[i] = stepSize;
156  }
157  }
158 
159 
160  const double * TrueMinimum() const {
161  return x0.GetMatrixArray();
162  }
163 
164 
165  void Gradient (const double * x, double * g) const {
166  gNCall2++;
167 
168  for (unsigned int i = 0; i < fDim ; ++i) {
169  cx [i] = std::cos(x[i]);
170  sx [i] = std::sin(x[i]);
171  }
172 
173  v = A*sx +B*cx;
174  r = v0-v;
175 
176 
177  // calculate the grad components
178  for (unsigned int i = 0; i < fDim ; ++i) {
179  g[i] = 0;
180  for (unsigned int k = 0; k < fDim ; ++k) {
181  g[i] += 2. * r(k) * ( - A(k,i) * cx(i) + B(k,i) * sx(i) );
182  }
183  }
184 
185  }
186 
187 #ifdef USE_FDF
188  void FdF (const double * x, double & f, double * g) const {
189  gNCall++;
190 
191  for (unsigned int i = 0; i < fDim ; ++i) {
192  cx [i] = std::cos(x[i]);
193  sx [i] = std::sin(x[i]);
194  }
195 
196  v = A*sx +B*cx;
197  r = v0-v;
198 
199  f = r * r;
200 
201 
202  // calculate the grad components
203  for (unsigned int i = 0; i < fDim ; ++i) {
204  g[i] = 0;
205  for (unsigned int k = 0; k < fDim ; ++k) {
206  g[i] += 2. * r(k) * ( - A(k,i) * cx(i) + B(k,i) * sx(i) );
207  }
208  }
209  }
210 #endif
211 
212  private:
213 
214 // TrigoFletcherFunction(const TrigoFletcherFunction & ) {}
215 // TrigoFletcherFunction & operator=(const TrigoFletcherFunction &) { return *this; }
216 
217  double DoEval (const double * x) const {
218  gNCall++;
219 
220 
221  for (unsigned int i = 0; i < fDim ; ++i) {
222  cx [i] = std::cos(x[i]);
223  sx [i] = std::sin(x[i]);
224  }
225 
226  v = A*sx +B*cx;
227  r = v0-v;
228 
229  return r * r;
230  }
231 
232 
233  double DoDerivative (const double * x, unsigned int i ) const {
234  std::vector<double> g(fDim);
235  Gradient(x,&g[0]);
236  return g[i];
237  }
238 
239 private:
240 
241  unsigned int fDim;
242 
243  TMatrixD A;
244  TMatrixD B;
245  TVectorD x0;
246  mutable TVectorD sx0;
247  mutable TVectorD cx0;
248  mutable TVectorD sx;
249  mutable TVectorD cx;
250  mutable TVectorD v0;
251  mutable TVectorD v;
252  mutable TVectorD r;
253 
254 
255 };
256 
257 
258 // CHEBYQUAD FUNCTION
259 
260 class ChebyQuadFunction : public ROOT::Math::IMultiGradFunction {
261 
262 public :
263 
264  ChebyQuadFunction(unsigned int n) :
265  fDim(n),
266  fvec(std::vector<double>(n) ),
267  fTrueMin(std::vector<double>(n) )
268  {
269  }
270 
271  unsigned int NDim() const { return fDim; }
272 
273  ROOT::Math::IMultiGenFunction * Clone() const {
274  return new ChebyQuadFunction(*this);
275  }
276 
277  const double * TrueMinimum() const {
278  return &fTrueMin[0];
279  }
280 
281  // use equally spaced points
282  void StartPoints(double * x, double * s) {
283  for (unsigned int i = 0; i < fDim; ++i) {
284  s[i] = 0.01;
285  x[i] = double(i)/(double(fDim)+1.0);
286  }
287  }
288 
289  // compute gradient
290 
291  void Gradient(const double * x, double * g) const {
292  gNCall2++;
293  unsigned int n = fDim;
294  // estimate first the fvec
295  DoCalculatefi(x);
296 
297  for (unsigned int j = 0; j < n; ++j) {
298  g[j] = 0.0;
299  double t1 = 1.0;
300  double t2 = 2.0 * x[j] - 1.0;
301  double t = 2.0 * t2;
302  double s1 = 0.0;
303  double s2 = 2.0;
304  for (unsigned int i = 0; i < n; ++i) {
305  g[j] += fvec[i] * s2;
306  double th = 4.0 * t2 + t * s2 - s1;
307  s1 = s2;
308  s2 = th;
309  th = t * t2 - t1;
310  t1 = t2;
311  t2 = th;
312  }
313  g[j] = 2. * g[j] / double(n);
314  }
315 
316 
317  }
318 
319  private:
320 
321  double DoEval (const double * x) const {
322 
323  gNCall++;
324  DoCalculatefi(x);
325  double f = 0;
326  for (unsigned int i = 0; i < fDim; ++i)
327  f += fvec[i] * fvec[i];
328 
329  return f;
330 
331  }
332 
333  double DoDerivative (const double * x, unsigned int i ) const {
334  std::vector<double> g(fDim);
335  Gradient(x,&g[0]);
336  return g[i];
337  }
338 
339  void DoCalculatefi(const double * x) const {
340  // calculate the i- element ; F(X) = Sum {fi]
341  unsigned int n = fDim;
342  for (unsigned int i = 0; i < n; ++i)
343  fvec[i] = 0;
344 
345  for (unsigned int j = 0; j < n; ++j) {
346  double t1 = 1.0;
347  double t2 = 2.0 * x[j] - 1.0;
348  double t = 2.0 * t2;
349  for (unsigned int i = 0; i < n; ++i) {
350  fvec[i] += t2;
351  double th = t * t2 - t1;
352  t1 = t2;
353  t2 = th;
354  }
355  }
356 
357  // sum with the integral (integral is zero for odd Cheb polynomial and = 1/(i**2 -1) for the even ones
358  for (unsigned int i = 1; i <= n; ++i) {
359  int l = i-1;
360  fvec[l] /= double ( n );
361  if ( ( i % 2 ) == 0 ) {
362  fvec[l] += 1.0 / ( double ( i*i ) - 1.0 );
363  }
364  }
365  }
366 
367  unsigned int fDim;
368  mutable std::vector<double> fvec;
369  mutable std::vector<double> fTrueMin;
370 };
371 
372 //WOOD function (4 dim function)
373 
374 double WoodFunction(const double * par) {
375  gNCall++;
376 
377  const double w = par[0];
378  const double x = par[1];
379  const double y = par[2];
380  const double z = par[3];
381 
382  const double w1 = w-1;
383  const double x1 = x-1;
384  const double y1 = y-1;
385  const double z1 = z-1;
386  const double tmp1 = x-w*w;
387  const double tmp2 = z-y*y;
388 
389  double f = 100*tmp1*tmp1+w1*w1+90*tmp2*tmp2+y1*y1+10.1*(x1*x1+z1*z1)+19.8*x1*z1;
390  return f;
391 }
392 
393 //Powell Function (4 dim function)
394 double PowellFunction(const double * par) {
395  gNCall++;
396 
397  const double w = par[0];
398  const double x = par[1];
399  const double y = par[2];
400  const double z = par[3];
401 
402  const double tmp1 = w+10*x;
403  const double tmp2 = y-z;
404  const double tmp3 = x-2*y;
405  const double tmp4 = w-z;
406 
407  double f = tmp1*tmp1+5*tmp2*tmp2+tmp3*tmp3*tmp3*tmp3+10*tmp4*tmp4*tmp4*tmp4;
408  return f;
409 }
410 
411 double SimpleQuadFunction(const double * par) {
412  gNCall++;
413  double x = par[0];
414  double y = par[1];
415  double f = x * x + 2 * y * y - x*y;
416  std::cout << "Quadfunc " << gNCall << "\t" << x << " , " << y << " f = " << f << std::endl;
417  return f;
418 }
419 
421 
422  const RosenBrockFunction * fRB = dynamic_cast< const RosenBrockFunction *> (&func);
423  if (fRB != 0) return fRB->TrueMinimum();
424  const TrigoFletcherFunction * fTF = dynamic_cast< const TrigoFletcherFunction *> (&func);
425  if (fTF != 0) return fTF->TrueMinimum();
426 // const ChebyQuadFunction * fCQ = dynamic_cast< const ChebyQuadFunction *> (&func);
427 // if (fCQ != 0) return fCQ->TrueMinimum();
428  return 0;
429 }
430 
431 void printMinimum(const std::vector<double> & x) {
432  std::cout << "Minimum X values\n";
433  std::cout << "\t";
434  int pr = std::cout.precision(12);
435  unsigned int n = x.size();
436  for (unsigned int i = 0; i < n; ++i) {
437  std::cout << x[i];
438  if ( i != n-1 ) std::cout << " , ";
439  if ( i > 0 && i % 6 == 0 ) std::cout << "\n\t";
440  }
441  std::cout << std::endl;
442  std::cout.precision(pr);
443 }
444 
445 int DoNewMinimization( const ROOT::Math::IMultiGenFunction & func, const double * x0, const double * s0, ROOT::Math::Minimizer * min, double &minval, double &edm, double * minx) {
446 
447  int iret = 0;
448 
449  min->SetMaxFunctionCalls(1000000);
450  min->SetMaxIterations(100000);
452  if (func.NDim() >= 10) {
453  min->SetTolerance(0.01);
454 
455  }
456 
457  min->SetPrintLevel(gVerbose);
458  // check if func provides gradient
459  const ROOT::Math::IMultiGradFunction * gfunc = dynamic_cast<const ROOT::Math::IMultiGradFunction *>(&func);
460  if (gfunc != 0 && useGradient)
461  min->SetFunction(*gfunc);
462  else
463  min->SetFunction(func);
464 
465  for (unsigned int i = 0; i < func.NDim(); ++i) {
466 #ifdef DEBUG
467  std::cout << "set variable " << i << " to value " << x0[i] << std::endl;
468 #endif
469  min->SetVariable(i,"x" + ROOT::Math::Util::ToString(i),x0[i], s0[i]);
470  }
471 
472  bool ret = min->Minimize();
473  minval = min->MinValue();
474  edm = min->Edm();
475 
476  if (!ret) {
477  return -1;
478  }
479 
480  const double * xmin = min->X();
481 
482  bool ok = true;
483  const double * trueMin = TrueMinimum(func);
484  if (trueMin != 0) {
485  for (unsigned int i = 0; i < func.NDim(); ++i)
486  ok &= (std::fabs(xmin[i]-trueMin[i] ) < gAbsTolerance);
487  }
488 
489  if (!ok) iret = -2;
490  int ncall_migrad = gNCall;
491 
492  // test Minos (use the default up of 1)
493  if (minos) {
494 
495  double el,eu;
496  for (unsigned int i = 0; i < func.NDim(); ++i) {
497  ret = min->GetMinosError(i,el,eu);
498  std::cout << "ncalls " << gNCall << "\t";
499  if (ret) std::cout << "MINOS error for " << i << " = " << el << " " << eu << std::endl;
500  else std::cout << "MINOS failed for " << i << std::endl;
501  }
502  }
503 
504  std::cout << "\nMigrad Ncalls:\t " << ncall_migrad << std::endl;
505  std::cout << "Minos Ncalls :\t " << gNCall - ncall_migrad << std::endl;
506 
507 
508 // std::cout << "function at the minimum " << func(xmin) << std::endl;
509  std::copy(xmin,xmin+func.NDim(),minx);
510  min->Clear();
511 
512  return iret;
513 }
514 
515 int DoOldMinimization( FCN func, TVirtualFitter * min, double &minval, double &edm) {
516 
517  int iret = 0;
518 
519  assert(min != 0);
520  min->SetFCN( func );
521 
522  Double_t arglist[100];
523  arglist[0] = gVerbose-1;
524  min->ExecuteCommand("SET PRINT",arglist,1);
525 
526  min->SetParameter(0,"x",-1.2,0.01,0,0);
527  min->SetParameter(1,"y", 1.0,0.01,0,0);
528 
529  arglist[0] = 0;
530  min->ExecuteCommand("SET NOW",arglist,0);
531  arglist[0] = 1000; // number of function calls
532  arglist[1] = gAbsTolerance; // tolerance
533  //min->ExecuteCommand("MIGRAD",arglist,0);
534  min->ExecuteCommand("MIGRAD",arglist,2);
535 
536  if (minos) min->ExecuteCommand("MINOS",arglist,0);
537 
538  Double_t parx,pary;
539  Double_t we,al,bl;
540  Char_t parName[32];
541  min->GetParameter(0,parName,parx,we,al,bl);
542  min->GetParameter(1,parName,pary,we,al,bl);
543 
544  bool ok = ( TMath::Abs(parx-1.) < gAbsTolerance &&
545  TMath::Abs(pary-1.) < gAbsTolerance );
546 
547 
548  double errdef = 0;
549  int nvpar, nparx;
550  min->GetStats(minval,edm,errdef,nvpar,nparx);
551  if (!ok) iret = -2;
552 
553  min->Clear(); // for further use
554  return iret;
555 
556 }
557 
558 
559 int testNewMinimizer( const ROOT::Math::IMultiGenFunction & func, const double * x0, const double * s0, const std::string & minimizer, const std::string & algoType) {
560 
561  std::cout << "\n************************************************************\n";
562  std::cout << "\tTest new ROOT::Math::Minimizer\n";
563  std::cout << "\tMinimizer is " << minimizer << " " << algoType << std::endl;
564 
565  int iret = 0;
566  double minval = 0.,edm = 0.;
567  std::vector<double> xmin(func.NDim() );
568 
569  TStopwatch w;
570  w.Start();
571 
572  ROOT::Math::Minimizer * min = ROOT::Math::Factory::CreateMinimizer(minimizer, algoType);
573  if (min == 0) {
574  std::cout << "Error using minimizer " << minimizer << " " << algoType << std::endl;
575  return -1;
576  }
577 
578 
579  for (int i = 0; i < gNmin; ++i) {
580  gNCall = 0; gNCall2 = 0;
581  iret |= DoNewMinimization(func, x0, s0, min,minval,edm,&xmin[0]);
582  }
583 
584  w.Stop();
585  if (iret != 0) std::cout << "\n****** ERROR: Minimization FAILED ! \n";
586  int pr = std::cout.precision(18);
587  std::cout << "\nNCalls: \t" << gNCall << " , " << gNCall2
588  << "\tMinValue: \t" << minval << "\tEdm: \t" << edm; std::cout.precision(pr);
589  std::cout << "\nTime: \t" << w.RealTime() << " , " << w.CpuTime() << std::endl;
590  printMinimum(xmin );
591  std::cout << "\n************************************************************\n";
592 
593 #ifdef CHECK_WITHMINUIT
594  // do Minuit after BFGS
595  if (minimizer == "GSL_BFGS") {
596  std::cout << "DO Minuit2 from last point\n";
597  gNCall = 0;
598  iret |= DoNewMinimization(func, &xmin.front(), s0, "Minuit2","",minval,edm,&xmin[0]);
599  int pr = std::cout.precision(18);
600  std::cout << "\nNCalls: \t" << gNCall << "\tMinValue: \t" << minval << "\tEdm: \t" << edm; std::cout.precision(pr);
601  std::cout << std::endl;
602  }
603 #endif
604 
605  delete min;
606 
607  return iret;
608 }
609 
610 
611 int testOldMinimizer( FCN func, const std::string & fitter, int n=25) {
612 
613  std::cout << "\n************************************************************\n";
614  std::cout << "\tTest using TVirtualFitter\n";
615  std::cout << "\tFitter is " << fitter << std::endl;
616 
617  int iret = 0;
618  double minval = 0.,edm = 0.;
619 
620  TStopwatch w;
621  w.Start();
622 
623  TVirtualFitter::SetDefaultFitter(fitter.c_str());
624 
626 
627  //min->Dump();
628 
629  for (int i = 0; i < gNmin; ++i) {
630  gNCall = 0;
631  iret |= DoOldMinimization(func, min,minval,edm);
632  }
633 
634  w.Stop();
635  if (iret != 0) std::cout << "\n****** ERROR: Minimization FAILED ! \n";
636  int pr = std::cout.precision(18);
637  std::cout << "\nNCalls: \t" << gNCall << "\tMinValue: \t" << minval << "\tEdm: \t" << edm; std::cout.precision(pr);
638  std::cout << "\nTime: \t" << w.RealTime() << " , " << w.CpuTime() << std::endl;
639  std::cout << "\n************************************************************\n";
640 
641  return iret;
642 }
643 
645 
646  int iret = 0;
647 
648 
649  std::cout << "\n************************************************************\n";
650  std::cout << "\tROSENBROCK function test\n\n";
651 
652  double s0[2] = {0.01,0.01};
653 
654  // minimize using Rosenbrock Function
655 #ifndef DEBUG
656  gNmin = 2000;
657 #endif
658 
659  gNmin = 1;
660 
661 
662  RosenBrockFunction fRB;
663  double xRB[2] = { -1.,1.2};
664 
665  iret |= testNewMinimizer(fRB,xRB,s0,"Minuit","");
666  iret |= testNewMinimizer(fRB,xRB,s0,"Minuit2","");
667 
668 // iret |= testNewMinimizer(fRB,xRB,s0,"GSLMultiMin","ConjugateFR");
669 // iret |= testNewMinimizer(fRB,xRB,s0,"GSLMultiMin","ConjugatePR");
670 // iret |= testNewMinimizer(fRB,xRB,s0,"GSLMultiMin","BFGS");
671 // iret |= testNewMinimizer(fRB,xRB,s0,"GSLMultiMin","BFGS2");
672 
673 // iret |= testOldMinimizer(RosenBrock,"Minuit",2);
674 // iret |= testOldMinimizer(RosenBrock,"Minuit2",2);
675 
676 
677  return iret;
678 }
679 
680 
682 
683  int iret = 0;
684 
685  // test with ChebyQuad function
686 
687  const int n = 8;
688  ChebyQuadFunction func(n);
689 
690 #ifndef DEBUG
691  gNmin = std::max(1, int(20000/n/n) );
692 #endif
693  gNmin = 1;
694 
695 
696  double s0[n];
697  double x0[n];
698  func.StartPoints(x0,s0);
699 
700  std::cout << "\n************************************************************\n";
701  std::cout << "\tCHEBYQUAD function test , n = " << n << "\n\n";
702 
703 
704 // double x[8] = {0.043153E+00, 0.193091E+00, 0.266329E+00, 0.500000E+00,
705 // 0.500000E+00, 0.733671E+00, 0.806910E+00, 0.956847E+00 };
706 // double x[2] = {0.5, 0.5001};
707 // std::cout << "FUNC " << func(x) << std::endl;
708  double x1[100] = { 0.00712780070646 , 0.0123441993113 , 0.0195428378255 , 0.0283679084192 , 0.0385291131289 , 0.0492202424892 , 0.0591277976178 ,
709  0.0689433195252 , 0.0791293590525 , 0.088794974369 , 0.0988949579193 , 0.108607151294 , 0.118571075831 ,
710  0.128605446238 , 0.137918291068 , 0.149177761352 , 0.156665324587 , 0.170851061982 , 0.174688134016 ,
711  0.192838903364 , 0.193078270803 , 0.209255377225 , 0.217740096779 , 0.225888518345 , 0.241031047421 ,
712  0.244253844041 , 0.257830449676 , 0.269467652526 , 0.274286498012 , 0.288877029988 , 0.297549406597 ,
713  0.304950954529 , 0.319230811642 , 0.326387092206 , 0.335229058731 , 0.349178359226 , 0.355905988048 ,
714  0.365197862755 , 0.379068092603 , 0.385826036925 , 0.394978252826 , 0.408974425717 , 0.415968185065 ,
715  0.424621041584 , 0.438837361714 , 0.446214149031 , 0.454242324351 , 0.468614308013 , 0.476506553416 ,
716  0.483916944941 , 0.498229247409 , 0.506794629616 , 0.513736742474 , 0.527712475478 , 0.537073277673 ,
717  0.543731917673 , 0.557187513963 , 0.567346279639 , 0.57379846397 , 0.586691058785 , 0.597561941009 ,
718  0.60382873461 , 0.616316037506 , 0.627719652101 , 0.633760038662 , 0.646175283836 , 0.657809344891 ,
719  0.663569004722 , 0.676314563639 , 0.687674566849 , 0.69332205923 , 0.706839545953 , 0.716907408637 ,
720  0.723407327715 , 0.738019389561 , 0.744806584048 , 0.754657613362 , 0.769181875619 , 0.772250323489 ,
721  0.787104833193 , 0.795856360905 , 0.804099304478 , 0.82142178741 , 0.819589601284 , 0.839024540481 ,
722  0.842457233039 , 0.857393475964 , 0.86408033345 , 0.876329840525 , 0.884867318008 , 0.895744532071 ,
723  0.905113958163 , 0.915445338697 , 0.925148068352 , 0.935344457785 , 0.945127838313 , 0.955272197168 ,
724  0.965687518559 , 0.975129521484 , 0.982662007764 };
725 
726  std::cout << "FUNC " << func(x1) << std::endl;
727 
728 
729  iret |= testNewMinimizer(func,x0,s0, "Minuit","");
730  iret |= testNewMinimizer(func,x0,s0, "Minuit2","");
731 
732 // iret |= testNewMinimizer(func,x0,s0, "GSLMultiMin","ConjugateFR");
733 // iret |= testNewMinimizer(func,x0,s0, "GSLMultiMin","ConjugatePR");
734 // iret |= testNewMinimizer(func,x0,s0, "GSLMultiMin","BFGS");
735 
736  return iret;
737 }
738 
739 
741 
742  int iret = 0;
743 
744 
745  // test with fletcher trigonometric function
746 #ifndef DEBUG
747  gNmin = 2;
748 #endif
749  gNmin = 1;
750 
751  const int nT = 50;
752  TrigoFletcherFunction fTrigo(nT);
753  double sTrigo[nT];
754  double xTrigo[nT];
755  fTrigo.StartPoints(xTrigo,sTrigo);
756 
757  std::cout << "\n************************************************************\n";
758  std::cout << "\tTRIGONOMETRIC Fletcher function test , n = " << nT << "\n\n";
759 
760 
761  iret |= testNewMinimizer(fTrigo,xTrigo,sTrigo,"Minuit2","");
762  iret |= testNewMinimizer(fTrigo,xTrigo,sTrigo,"Minuit","");
763 
764 // iret |= testNewMinimizer(fTrigo,xTrigo,sTrigo,"GSLMultiMin","ConjugateFR");
765 // iret |= testNewMinimizer(fTrigo,xTrigo,sTrigo,"GSLMultiMin","ConjugatePR");
766 // iret |= testNewMinimizer(fTrigo,xTrigo,sTrigo,"GSLMultiMin","BFGS");
767 
768 
769  return iret;
770 }
771 
772 
773 int testWood() {
774 
775  int iret = 0;
776 
777 
778  // test with Wood function (4d)
779 // minimum : F(1,1,1,1) = 0.
780 
781 
782 #ifndef DEBUG
783  gNmin = 1000;
784 #endif
785  gNmin = 1;
786 
788 
789  double x0[4] = { -3, -1, -3, -1 };
790  double s0[4] = { 0.1, 0.1, 0.1, 0.1};
791 
792  std::cout << "\n************************************************************\n";
793  std::cout << "\tWOOD 4 function test \n\n";
794 
795 
796  iret |= testNewMinimizer(f, x0, s0,"Minuit","");
797  iret |= testNewMinimizer(f, x0, s0,"Minuit2","");
798 
799 
800  return iret;
801 }
802 
803 int testPowell() {
804 
805  int iret = 0;
806 
807 
808  // test with Powell function (4d)
809  // minimum is at F(0,0,0,0) = 0
810 #ifndef DEBUG
811  gNmin = 1000;
812 #endif
813  gNmin = 1;
814 
816 
817  double x0[4] = { -3, -1, 0, 1 };
818  double s0[4] = { 0.1, 0.1, 0.1, 0.1};
819 
820  std::cout << "\n************************************************************\n";
821  std::cout << "\tPOWELL function test \n\n";
822 
823 
824  iret |= testNewMinimizer(f, x0, s0,"Minuit","");
825  iret |= testNewMinimizer(f, x0, s0,"Minuit2","");
826 
827 
828  return iret;
829 }
830 
831 
833 
834  int iret = 0;
835 
836 
837  // test with a simple quadratic function 2d
838  // minimum is at F(0,0) = 0
839 #ifndef DEBUG
840  gNmin = 1000;
841 #endif
842  gNmin = 1;
843 
845 
846  double x0[4] = { -3, -3 };
847  double s0[4] = { 0.1, 0.1};
848 
849  std::cout << "\n************************************************************\n";
850  std::cout << "\tSIMPLE QUAD function test \n\n";
851 
852 
853  iret |= testNewMinimizer(f, x0, s0,"Minuit","");
854  iret |= testNewMinimizer(f, x0, s0,"Minuit2","");
855 
856 
857  return iret;
858 }
859 
860 
861 int main() {
862 
863  int iret = 0;
864 
865 #ifdef DEBUG
866  gVerbose = 4;
867  gNmin = 1;
868 #endif
869 
870 
871  iret |= testRosenBrock();
872  iret |= testChebyQuad();
873  iret |= testTrigoFletcher();
874 
875  iret |= testWood();
876  iret |= testPowell();
877 
878  iret |= testQuadFunc();
879 
880 
881 
882  if (iret != 0)
883  std::cerr << "testMinim :\t FAILED " << std::endl;
884  else
885  std::cerr << "testMinim :\t OK " << std::endl;
886  return iret;
887 
888 }
Interface (abstract class) for multi-dimensional functions providing a gradient calculation.
Definition: IFunction.h:322
double par[1]
Definition: unuranDistr.cxx:38
int gNCall
static double B[]
void SetMaxIterations(unsigned int maxiter)
set maximum iterations (one iteration can have many function calls)
Definition: Minimizer.h:459
float xmin
Definition: THbookFile.cxx:93
Random number generator class based on M.
Definition: TRandom3.h:29
Double_t RealTime()
Stop the stopwatch (if it is running) and return the realtime (in seconds) passed between the start a...
Definition: TStopwatch.cxx:110
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:58
const double pi
double PowellFunction(const double *par)
Documentation for class Functor class.
Definition: Functor.h:394
virtual Int_t SetParameter(Int_t ipar, const char *parname, Double_t value, Double_t verr, Double_t vlow, Double_t vhigh)=0
int testChebyQuad()
static void SetDefaultFitter(const char *name="")
static: set name of default fitter
int DoOldMinimization(FCN func, TVirtualFitter *min, double &minval, double &edm)
TVectorT.
Definition: TMatrixTBase.h:89
void printMinimum(const std::vector< double > &x)
int testRosenBrock()
virtual Double_t GetParameter(Int_t ipar) const =0
int testWood()
static ROOT::Math::Minimizer * CreateMinimizer(const std::string &minimizerType="", const std::string &algoType="")
static method to create the corrisponding Minimizer given the string Supported Minimizers types are: ...
Definition: Factory.cxx:63
int DoNewMinimization(const ROOT::Math::IMultiGenFunction &func, const double *x0, const double *s0, ROOT::Math::Minimizer *min, double &minval, double &edm, double *minx)
void RosenBrock(Int_t &, Double_t *, Double_t &f, Double_t *par, Int_t)
int gVerbose
Double_t CpuTime()
Stop the stopwatch (if it is running) and return the cputime (in seconds) passed between the start an...
Definition: TStopwatch.cxx:125
int Int_t
Definition: RtypesCore.h:41
int gNmin
TLatex * t1
Definition: textangle.C:20
static double A[]
Short_t Abs(Short_t d)
Definition: TMathBase.h:110
double cos(double)
TMatrixT.
Definition: TMatrixDfwd.h:24
virtual void Clear()
reset for consecutive minimizations - implement if needed
Definition: Minimizer.h:125
double gAbsTolerance
void Stop()
Stop the stopwatch.
Definition: TStopwatch.cxx:77
Double_t x[n]
Definition: legend1.C:17
Abstract Minimizer class, defining the interface for the various minimizer (like Minuit2, Minuit, GSL, etc..) Plug-in&#39;s exist in ROOT to be able to instantiate the derived classes like ROOT::Math::GSLMinimizer or ROOT::Math::Minuit2Minimizer via the plug-in manager.
Definition: Minimizer.h:86
virtual double MinValue() const =0
return minimum function value
virtual bool Minimize()=0
method to perform the minimization
int testTrigoFletcher()
virtual void SetFCN(void *fcn) R__DEPRECATED(6
To set the address of the minimization objective function.
virtual const double * X() const =0
return pointer to X values at the minimum
virtual void Clear(Option_t *option="")=0
Set name and title to empty strings ("").
int testPowell()
double sin(double)
int main()
void(* FCN)(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag)
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
virtual void SetFunction(const ROOT::Math::IMultiGenFunction &func)=0
set the function to minimize
const double * TrueMinimum(const ROOT::Math::IMultiGenFunction &func)
TRandom2 r(17)
SVector< double, 2 > v
Definition: Dict.h:5
virtual unsigned int NDim() const =0
Retrieve the dimension of the function.
virtual double Edm() const
return expected distance reached from the minimum (re-implement if minimizer provides it ...
Definition: Minimizer.h:262
int testOldMinimizer(FCN func, const std::string &fitter, int n=25)
double SimpleQuadFunction(const double *par)
TLine * l
Definition: textangle.C:4
virtual Int_t ExecuteCommand(const char *command, Double_t *args, Int_t nargs)=0
virtual Int_t GetStats(Double_t &amin, Double_t &edm, Double_t &errdef, Int_t &nvpar, Int_t &nparx) const =0
Double_t Pi()
Definition: TMath.h:44
void SetMaxFunctionCalls(unsigned int maxfcn)
set maximum of function calls
Definition: Minimizer.h:456
static const double x1[5]
double f(double x)
double Double_t
Definition: RtypesCore.h:55
void SetTolerance(double tol)
set the tolerance
Definition: Minimizer.h:462
virtual bool GetMinosError(unsigned int ivar, double &errLow, double &errUp, int option=0)
minos error for variable i, return false if Minos failed or not supported and the lower and upper err...
Definition: Minimizer.h:357
Double_t y[n]
Definition: legend1.C:17
double func(double *x, double *p)
Definition: stressTF1.cxx:213
virtual Double_t Uniform(Double_t x1=1)
Returns a uniform deviate on the interval (0, x1).
Definition: TRandom.cxx:606
static const double eu
Definition: Vavilov.cxx:52
bool useGradient
Abstract Base Class for Fitting.
you should not use this method at all Int_t Int_t z
Definition: TRolke.cxx:630
char Char_t
Definition: RtypesCore.h:29
typedef void((*Func_t)())
virtual bool SetVariable(unsigned int ivar, const std::string &name, double val, double step)=0
set a new free variable
static TVirtualFitter * Fitter(TObject *obj, Int_t maxpar=25)
Static function returning a pointer to the current fitter.
std::string ToString(const T &val)
Utility function for conversion to strings.
Definition: Util.h:42
int gNCall2
int testQuadFunc()
void SetPrintLevel(int level)
set print level
Definition: Minimizer.h:453
bool minos
Documentation for the abstract class IBaseFunctionMultiDim.
Definition: IFunction.h:63
const Int_t n
Definition: legend1.C:16
virtual double DoEval(const double *x) const =0
Implementation of the evaluation function.
int testNewMinimizer(const ROOT::Math::IMultiGenFunction &func, const double *x0, const double *s0, const std::string &minimizer, const std::string &algoType)
virtual IBaseFunctionMultiDim * Clone() const =0
Clone a function.
double WoodFunction(const double *par)
Stopwatch class.
Definition: TStopwatch.h:30