ROOT  6.07/01
Reference Guide
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
testMinim.cxx
Go to the documentation of this file.
1 // test of minimization usingnew minimizer classes
2 
3 #include "Math/Minimizer.h"
4 #include "Math/Factory.h"
5 
6 #include "TVirtualFitter.h"
7 
8 #include "Math/IFunction.h"
9 #include "Math/Util.h"
10 #include <cmath>
11 #include <cassert>
12 
13 #include <string>
14 #include <iostream>
15 
16 #include "TStopwatch.h"
17 #include "TMatrixD.h"
18 #include "TVectorD.h"
19 #include "TRandom3.h"
20 #include "TMath.h"
21 
22 #include "RVersion.h"
23 #include "RConfigure.h"
24 
25 //#define DEBUG
26 
27 int gNCall = 0;
28 int gNCall2 = 0;
29 int gNmin = 1;
30 int gVerbose = 1;
31 
32 
33 bool minos = false;
34 
35 double gAbsTolerance = 5.E-6; // otherwise gsl_conjugate_PR fails
36 
37 // Rosenbrok function to be minimize
38 
39 typedef void (*FCN)(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag);
40 
41 
42 
43 // ROSENBROCK function
44 ////////////////////////////////////////////////////////////////////////////////
45 
46 void RosenBrock(Int_t &, Double_t *, Double_t &f, Double_t *par, Int_t /*iflag*/)
47 {
48  gNCall++;
49  const Double_t x = par[0];
50  const Double_t y = par[1];
51  const Double_t tmp1 = y-x*x;
52  const Double_t tmp2 = 1-x;
53  f = 100*tmp1*tmp1+tmp2*tmp2;
54 }
55 
56 
57 
58 class RosenBrockFunction : public ROOT::Math::IMultiGenFunction {
59 
60 public :
61 
62 
63  unsigned int NDim() const { return 2; }
64 
66  return new RosenBrockFunction();
67  }
68 
69  const double * TrueXMinimum() const {
70  fTrueMin[0] = 1;
71  fTrueMin[1] = 1;
72  return fTrueMin;
73  }
74 
75  double TrueMinimum() const {
76  return 0;
77  }
78 
79  private:
80 
81  inline double DoEval (const double * x) const {
82 #ifdef USE_FREE_FUNC
83  double f = 0;
84  int ierr = 0;
85  int i = 0;
86  RosenBrock(i,0,f,const_cast<double *>(x),ierr);
87  return f;
88 #else
89  gNCall++;
90  const Double_t xx = x[0];
91  const Double_t yy = x[1];
92  const Double_t tmp1 = yy-xx*xx;
93  const Double_t tmp2 = 1-xx;
94  return 100*tmp1*tmp1+tmp2*tmp2;
95 #endif
96  }
97 
98 
99  mutable double fTrueMin[2];
100 };
101 
102 // TRIGONOMETRIC FLETCHER FUNCTION
103 
104 class TrigoFletcherFunction : public ROOT::Math::IMultiGradFunction {
105 
106 public :
107 
108 
109  TrigoFletcherFunction(unsigned int dim) : fDim(dim) {
110  double seed = 3;
111  A.ResizeTo(dim,dim);
112  B.ResizeTo(dim,dim);
113  x0.ResizeTo(dim);
114  sx0.ResizeTo(dim);
115  cx0.ResizeTo(dim);
116  sx.ResizeTo(dim);
117  cx.ResizeTo(dim);
118  v0.ResizeTo(dim);
119  v.ResizeTo(dim);
120  r.ResizeTo(dim);
121  A.Randomize(-100.,100,seed);
122  B.Randomize(-100.,100,seed);
123  for (unsigned int i = 0; i < dim; i++) {
124  for (unsigned int j = 0; j < dim; j++) {
125  A(i,j) = int(A(i,j));
126  B(i,j) = int(B(i,j));
127  }
128  }
129  x0.Randomize(-TMath::Pi(),TMath::Pi(),seed);
130  // calculate vector Ei
131  for (unsigned int i = 0; i < fDim ; ++i) {
132  cx0[i] = std::cos(x0[i]);
133  sx0[i] = std::sin(x0[i]);
134  }
135  v0 = A*sx0+B*cx0;
136  }
137 
138 
139  unsigned int NDim() const { return fDim; }
140 
142  TrigoFletcherFunction * f = new TrigoFletcherFunction(*this);
143 // std::cerr <<"cannot clone this function" << std::endl;
144 // assert(0);
145  return f;
146  }
147 
148 
149  void StartPoints(double * x, double * s) {
150  TRandom3 rndm;
151  const double stepSize = 0.01;
152  const double deltaAmp = 0.1;
153  const double pi = TMath::Pi();
154  for (unsigned int i = 0; i < fDim; ++i) {
155  double delta = rndm.Uniform(-deltaAmp*pi,deltaAmp*pi);
156  x[i] = x0(i) + 0.1*delta;
157  if (x[i] <= - pi) x[i] += 2.*pi;
158  if (x[i] > pi) x[i] -= 2.*pi;
159  s[i] = stepSize;
160  }
161  }
162 
163 
164  const double * TrueXMinimum() const {
165  return x0.GetMatrixArray();
166  }
167 
168  double TrueMinimum() const { return 0; }
169 
170  void Gradient (const double * x, double * g) const {
171  gNCall2++;
172 
173  for (unsigned int i = 0; i < fDim ; ++i) {
174  cx [i] = std::cos(x[i]);
175  sx [i] = std::sin(x[i]);
176  }
177 
178  v = A*sx +B*cx;
179  r = v0-v;
180 
181 
182  // calculate the grad components
183  for (unsigned int i = 0; i < fDim ; ++i) {
184  g[i] = 0;
185  for (unsigned int k = 0; k < fDim ; ++k) {
186  g[i] += 2. * r(k) * ( - A(k,i) * cx(i) + B(k,i) * sx(i) );
187  }
188  }
189 
190  }
191 
192 #ifdef USE_FDF
193  void FdF (const double * x, double & f, double * g) const {
194  gNCall++;
195 
196  for (unsigned int i = 0; i < fDim ; ++i) {
197  cx [i] = std::cos(x[i]);
198  sx [i] = std::sin(x[i]);
199  }
200 
201  v = A*sx +B*cx;
202  r = v0-v;
203 
204  f = r * r;
205 
206 
207  // calculate the grad components
208  for (unsigned int i = 0; i < fDim ; ++i) {
209  g[i] = 0;
210  for (unsigned int k = 0; k < fDim ; ++k) {
211  g[i] += 2. * r(k) * ( - A(k,i) * cx(i) + B(k,i) * sx(i) );
212  }
213  }
214  }
215 #endif
216 
217  private:
218 
219 // TrigoFletcherFunction(const TrigoFletcherFunction & ) {}
220 // TrigoFletcherFunction & operator=(const TrigoFletcherFunction &) { return *this; }
221 
222  double DoEval (const double * x) const {
223  gNCall++;
224 
225 
226  for (unsigned int i = 0; i < fDim ; ++i) {
227  cx [i] = std::cos(x[i]);
228  sx [i] = std::sin(x[i]);
229  }
230 
231  v = A*sx +B*cx;
232  r = v0-v;
233 
234  return r * r;
235  }
236 
237 
238  double DoDerivative (const double * x, unsigned int i ) const {
239  std::vector<double> g(fDim);
240  Gradient(x,&g[0]);
241  return g[i];
242  }
243 
244 private:
245 
246  unsigned int fDim;
247 
248  TMatrixD A;
249  TMatrixD B;
250  TVectorD x0;
251  mutable TVectorD sx0;
252  mutable TVectorD cx0;
253  mutable TVectorD sx;
254  mutable TVectorD cx;
255  mutable TVectorD v0;
256  mutable TVectorD v;
257  mutable TVectorD r;
258 
259 
260 };
261 
262 
263 // CHEBYQUAD FUNCTION
264 
265 class ChebyQuadFunction : public ROOT::Math::IMultiGradFunction {
266 
267 public :
268 
269  ChebyQuadFunction(unsigned int n) :
270  fDim(n),
271  fvec(std::vector<double>(n) ),
272  fTrueMin(std::vector<double>(n) )
273  {
274  }
275 
276  unsigned int NDim() const { return fDim; }
277 
279  return new ChebyQuadFunction(*this);
280  }
281 
282  const double * TrueXMinimum() const {
283  return &fTrueMin[0];
284  }
285 
286  double TrueMinimum() const { return 0; }
287 
288  // use equally spaced points
289  void StartPoints(double * x, double * s) {
290  for (unsigned int i = 0; i < fDim; ++i) {
291  s[i] = 0.01;
292  x[i] = double(i)/(double(fDim)+1.0);
293  }
294  }
295 
296  // compute gradient
297 
298  void Gradient(const double * x, double * g) const {
299  gNCall2++;
300  unsigned int n = fDim;
301  // estimate first the fvec
302  DoCalculatefi(x);
303 
304  for (unsigned int j = 0; j < n; ++j) {
305  g[j] = 0.0;
306  double t1 = 1.0;
307  double t2 = 2.0 * x[j] - 1.0;
308  double t = 2.0 * t2;
309  double s1 = 0.0;
310  double s2 = 2.0;
311  for (unsigned int i = 0; i < n; ++i) {
312  g[j] += fvec[i] * s2;
313  double th = 4.0 * t2 + t * s2 - s1;
314  s1 = s2;
315  s2 = th;
316  th = t * t2 - t1;
317  t1 = t2;
318  t2 = th;
319  }
320  g[j] = 2. * g[j] / double(n);
321  }
322 
323 
324  }
325 
326  private:
327 
328  double DoEval (const double * x) const {
329 
330  gNCall++;
331  DoCalculatefi(x);
332  double f = 0;
333  for (unsigned int i = 0; i < fDim; ++i)
334  f += fvec[i] * fvec[i];
335 
336  return f;
337 
338  }
339 
340  double DoDerivative (const double * x, unsigned int i ) const {
341  std::vector<double> g(fDim);
342  Gradient(x,&g[0]);
343  return g[i];
344  }
345 
346  void DoCalculatefi(const double * x) const {
347  // calculate the i- element ; F(X) = Sum {fi]
348  unsigned int n = fDim;
349  for (unsigned int i = 0; i < n; ++i)
350  fvec[i] = 0;
351 
352  for (unsigned int j = 0; j < n; ++j) {
353  double t1 = 1.0;
354  double t2 = 2.0 * x[j] - 1.0;
355  double t = 2.0 * t2;
356  for (unsigned int i = 0; i < n; ++i) {
357  fvec[i] += t2;
358  double th = t * t2 - t1;
359  t1 = t2;
360  t2 = th;
361  }
362  }
363 
364  // sum with the integral (integral is zero for odd Cheb polynomial and = 1/(i**2 -1) for the even ones
365  for (unsigned int i = 1; i <= n; ++i) {
366  int l = i-1;
367  fvec[l] /= double ( n );
368  if ( ( i % 2 ) == 0 ) {
369  fvec[l] += 1.0 / ( double ( i*i ) - 1.0 );
370  }
371  }
372  }
373 
374  unsigned int fDim;
375  mutable std::vector<double> fvec;
376  mutable std::vector<double> fTrueMin;
377 };
378 
379 
380 
382 
383  const RosenBrockFunction * fRB = dynamic_cast< const RosenBrockFunction *> (&func);
384  if (fRB != 0) return fRB->TrueXMinimum();
385  const TrigoFletcherFunction * fTF = dynamic_cast< const TrigoFletcherFunction *> (&func);
386  if (fTF != 0) return fTF->TrueXMinimum();
387 // const ChebyQuadFunction * fCQ = dynamic_cast< const ChebyQuadFunction *> (&func);
388 // if (fCQ != 0) return fCQ->TrueXMinimum();
389  return 0;
390 }
392 
393  const RosenBrockFunction * fRB = dynamic_cast< const RosenBrockFunction *> (&func);
394  if (fRB != 0) return fRB->TrueMinimum();
395  const TrigoFletcherFunction * fTF = dynamic_cast< const TrigoFletcherFunction *> (&func);
396  if (fTF != 0) return fTF->TrueMinimum();
397 // const ChebyQuadFunction * fCQ = dynamic_cast< const ChebyQuadFunction *> (&func);
398 // if (fCQ != 0) return fCQ->TrueXMinimum();
399  return 0;
400 }
401 
402 void printMinimum(const std::vector<double> & x) {
403  std::cout << "Minimum X values\n";
404  std::cout << "\t";
405  int pr = std::cout.precision(12);
406  unsigned int n = x.size();
407  for (unsigned int i = 0; i < n; ++i) {
408  std::cout << x[i];
409  if ( i != n-1 ) std::cout << " , ";
410  if ( i > 0 && i % 6 == 0 ) std::cout << "\n\t";
411  }
412  std::cout << std::endl;
413  std::cout.precision(pr);
414 }
415 
416 int DoNewMinimization( const ROOT::Math::IMultiGenFunction & func, const double * x0, const double * s0, ROOT::Math::Minimizer * min, double &minval, double &edm, double * minx) {
417 
418  int iret = 0;
419 
420  if (func.NDim() >= 10) {
421  min->SetMaxFunctionCalls(1000000);
422  min->SetMaxIterations(100000);
423  min->SetTolerance(0.01);
424  }
425  else
426  min->SetTolerance(0.001);
427 
428 
429  min->SetPrintLevel(gVerbose);
430  // check if func provides gradient
431  const ROOT::Math::IMultiGradFunction * gfunc = dynamic_cast<const ROOT::Math::IMultiGradFunction *>(&func);
432  if (gfunc != 0)
433  min->SetFunction(*gfunc);
434  else
435  min->SetFunction(func);
436 
437  for (unsigned int i = 0; i < func.NDim(); ++i) {
438  min->SetVariable(i,"x" + ROOT::Math::Util::ToString(i),x0[i], s0[i]);
439  }
440 
441  bool ret = min->Minimize();
442  minval = min->MinValue();
443  edm = min->Edm();
444 
445  if (!ret) {
446  delete min;
447  return -1;
448  }
449 
450  const double * xmin = min->X();
451 
452  bool ok = true;
453  const double * trueMin = TrueXMinimum(func);
454  if (trueMin != 0) {
455  ok &= (std::fabs(minval - TrueMinimum(func) ) < gAbsTolerance );
456  for (unsigned int i = 0; i < func.NDim(); ++i)
457  ok &= (std::fabs(xmin[i]-trueMin[i] ) < sqrt(gAbsTolerance));
458  }
459 
460  if (!ok) iret = -2;
461 
462  // test Minos (use the default up of 1)
463  if (minos) {
464 
465  double el,eu;
466  for (unsigned int i = 0; i < func.NDim(); ++i) {
467  ret = min->GetMinosError(i,el,eu);
468  if (ret) std::cout << "MINOS error for " << i << " = " << el << " " << eu << std::endl;
469  else std::cout << "MINOS failed for " << i << std::endl;
470  }
471  }
472 
473 #ifdef DEBUG
474  std::cout << "ncalls = " << min->NCalls() << std::endl;
475 #endif
476 
477 // std::cout << "function at the minimum " << func(xmin) << std::endl;
478  std::copy(xmin,xmin+func.NDim(),minx);
479  min->Clear();
480 
481  return iret;
482 }
483 
484 int DoOldMinimization( FCN func, TVirtualFitter * min, double &minval, double &edm) {
485 
486  int iret = 0;
487 
488  assert(min != 0);
489  min->SetFCN( func );
490 
491  Double_t arglist[100];
492  arglist[0] = gVerbose-1;
493  min->ExecuteCommand("SET PRINT",arglist,1);
494 
495  min->SetParameter(0,"x",-1.2,0.01,0,0);
496  min->SetParameter(1,"y", 1.0,0.01,0,0);
497 
498  arglist[0] = 0;
499  min->ExecuteCommand("SET NOW",arglist,0);
500  arglist[0] = 1000; // number of function calls
501  arglist[1] = 0.001; // tolerance
502  //min->ExecuteCommand("MIGRAD",arglist,0);
503  min->ExecuteCommand("MIGRAD",arglist,2);
504 
505  if (minos) min->ExecuteCommand("MINOS",arglist,0);
506 
507  Double_t parx,pary;
508  Double_t we,al,bl;
509  Char_t parName[32];
510  min->GetParameter(0,parName,parx,we,al,bl);
511  min->GetParameter(1,parName,pary,we,al,bl);
512 
513  bool ok = ( TMath::Abs(parx-1.) < sqrt(gAbsTolerance) &&
514  TMath::Abs(pary-1.) < sqrt(gAbsTolerance) );
515 
516 
517 
518  double errdef = 0;
519  int nvpar, nparx;
520  min->GetStats(minval,edm,errdef,nvpar,nparx);
521 
522  if (!ok) iret = -2;
523 
524  min->Clear(); // for further use
525  return iret;
526 
527 }
528 
529 
530 int testNewMinimizer( const ROOT::Math::IMultiGenFunction & func, const double * x0, const double * s0, const std::string & minimizer, const std::string & algoType) {
531 
532  std::cout << "\n************************************************************\n";
533  std::cout << "\tTest new ROOT::Math::Minimizer\n";
534  std::cout << "\tMinimizer is " << minimizer << " " << algoType << std::endl;
535 
536  int iret = 0;
537  double minval = 0., edm = 0.;
538  std::vector<double> xmin(func.NDim() );
539 
540  TStopwatch w;
541  w.Start();
542 
544  if (min == 0) {
545  std::cout << "Error using minimizer " << minimizer << " " << algoType << std::endl;
546  return -1;
547  }
548 
549  for (int i = 0; i < gNmin; ++i) {
550  gNCall = 0; gNCall2 = 0;
551  iret |= DoNewMinimization(func, x0, s0, min,minval,edm,&xmin[0]);
552  }
553 
554  w.Stop();
555  if (iret != 0) std::cout << "\n****** ERROR: Minimization FAILED ! \n";
556  int pr = std::cout.precision(18);
557  std::cout << "\nNCalls: \t" << gNCall << " , " << gNCall2
558  << "\tMinValue: \t" << minval << "\tEdm: \t" << edm; std::cout.precision(pr);
559  std::cout << "\nTime: \t" << w.RealTime() << " , " << w.CpuTime() << std::endl;
560  printMinimum(xmin );
561  std::cout << "\n************************************************************\n";
562 
563 #ifdef CHECK_WITHMINUIT
564  // do Minuit after BFGS
565  if (minimizer == "GSL_BFGS") {
566  std::cout << "DO Minuit2 from last point\n";
567  gNCall = 0;
568  iret |= DoNewMinimization(func, &xmin.front(), s0, "Minuit2","",minval,edm,&xmin[0]);
569  int pr = std::cout.precision(18);
570  std::cout << "\nNCalls: \t" << gNCall << "\tMinValue: \t" << minval << "\tEdm: \t" << edm; std::cout.precision(pr);
571  std::cout << std::endl;
572  }
573 #endif
574 
575  delete min;
576 
577  return iret;
578 }
579 
580 
581 int testOldMinimizer( FCN func, const std::string & fitter, int n=25) {
582 
583  std::cout << "\n************************************************************\n";
584  std::cout << "\tTest using TVirtualFitter\n";
585  std::cout << "\tFitter is " << fitter << std::endl;
586 
587  int iret = 0;
588  double minval = 0.,edm = 0.;
589 
590  TStopwatch w;
591  w.Start();
592 
593  TVirtualFitter::SetDefaultFitter(fitter.c_str());
594 
596 
597  //min->Dump();
598 
599  for (int i = 0; i < gNmin; ++i) {
600  gNCall = 0;
601  iret |= DoOldMinimization(func, min,minval,edm);
602  }
603 
604  w.Stop();
605  if (iret != 0) std::cout << "\n****** ERROR: Minimization FAILED ! \n";
606  int pr = std::cout.precision(18);
607  std::cout << "\nNCalls: \t" << gNCall << "\tMinValue: \t" << minval << "\tEdm: \t" << edm; std::cout.precision(pr);
608  std::cout << "\nTime: \t" << w.RealTime() << " , " << w.CpuTime() << std::endl;
609  std::cout << "\n************************************************************\n";
610 
611  return iret;
612 }
613 
615 
616  int iret = 0;
617 
618 
619  std::cout << "\n************************************************************\n";
620  std::cout << "\tROSENBROCK function test\n\n";
621 
622  double s0[2] = {0.01,0.01};
623 
624  // minimize using Rosenbrock Function
625 #ifndef DEBUG
626  gNmin = 1;
627 #endif
628 
629 
630 #if ROOT_VERSION_CODE < ROOT_VERSION(5,99,00)
631  iret |= testOldMinimizer(RosenBrock,"Minuit",2);
632  iret |= testOldMinimizer(RosenBrock,"Minuit2",2);
633 #endif
634 
635  RosenBrockFunction fRB;
636  double xRB[2] = { -1.,1.2};
637  iret |= testNewMinimizer(fRB,xRB,s0,"Minuit","");
638  iret |= testNewMinimizer(fRB,xRB,s0,"Minuit2","");
639 #ifdef R__HAS_MATHMORE
640  iret |= testNewMinimizer(fRB,xRB,s0,"GSLMultiMin","ConjugateFR");
641  iret |= testNewMinimizer(fRB,xRB,s0,"GSLMultiMin","ConjugatePR");
642  iret |= testNewMinimizer(fRB,xRB,s0,"GSLMultiMin","BFGS");
643  iret |= testNewMinimizer(fRB,xRB,s0,"GSLMultiMin","BFGS2");
644  //iret |= testNewMinimizer(fRB,xRB,s0,"Genetic","");
645 #endif
646 
647  return iret;
648 }
649 
650 
652 
653  int iret = 0;
654 
655 
656  // test with fletcher trigonometric function
657 #ifndef DEBUG
658  gNmin = 1;
659 #endif
660 
661  const int nT = 50;
662  TrigoFletcherFunction fTrigo(nT);
663  double sTrigo[nT];
664  double xTrigo[nT];
665  fTrigo.StartPoints(xTrigo,sTrigo);
666 
667  std::cout << "\n************************************************************\n";
668  std::cout << "\tTRIGONOMETRIC Fletcher function test , n = " << nT << "\n\n";
669 
670 
671  iret |= testNewMinimizer(fTrigo,xTrigo,sTrigo,"Minuit2","");
672  iret |= testNewMinimizer(fTrigo,xTrigo,sTrigo,"Minuit","");
673 #ifdef R__HAS_MATHMORE
674  iret |= testNewMinimizer(fTrigo,xTrigo,sTrigo,"GSLMultiMin","ConjugateFR");
675  iret |= testNewMinimizer(fTrigo,xTrigo,sTrigo,"GSLMultiMin","ConjugatePR");
676  iret |= testNewMinimizer(fTrigo,xTrigo,sTrigo,"GSLMultiMin","BFGS");
677 #endif
678 
679  return iret;
680 }
681 
683 
684  int iret = 0;
685 
686  // test with ChebyQuad function
687 
688  const int n = 8;
689  ChebyQuadFunction func(n);
690 
691 #ifndef DEBUG
692  gNmin = std::max(1, int(1000/n/n) );
693 #endif
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, "Minuit2","");
730  iret |= testNewMinimizer(func,x0,s0, "Minuit","");
731 #ifdef R__HAS_MATHMORE
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 #endif
736 
737  return iret;
738 }
739 
740 int main() {
741 
742  int iret = 0;
743 
744 #ifdef DEBUG
745  gVerbose = 3;
746  gNmin = 1;
747 #endif
748 
749  iret |= testRosenBrock();
750 // iret |= testChebyQuad();
751 // iret |= testTrigoFletcher();
752 
753 
754 
755  if (iret != 0)
756  std::cerr << "testMinim :\t FAILED " << std::endl;
757  else
758  std::cerr << "testMinim :\t OK " << std::endl;
759  return iret;
760 
761 }
int testOldMinimizer(FCN func, const std::string &fitter, int n=25)
Definition: testMinim.cxx:581
void RosenBrock(Int_t &, Double_t *, Double_t &f, Double_t *par, Int_t)
Definition: testMinim.cxx:46
Interface (abstract class) for multi-dimensional functions providing a gradient calculation.
Definition: IFunction.h:322
double par[1]
Definition: unuranDistr.cxx:38
static double B[]
int gVerbose
Definition: testMinim.cxx:30
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
static Vc_ALWAYS_INLINE int_v min(const int_v &x, const int_v &y)
Definition: vector.h:433
Double_t RealTime()
Stop the stopwatch (if it is running) and return the realtime (in seconds) passed between the start a...
Definition: TStopwatch.cxx:108
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:56
double TrueMinimum(const ROOT::Math::IMultiGenFunction &func)
Definition: testMinim.cxx:391
const double pi
virtual Int_t SetParameter(Int_t ipar, const char *parname, Double_t value, Double_t verr, Double_t vlow, Double_t vhigh)=0
#define assert(cond)
Definition: unittest.h:542
static void SetDefaultFitter(const char *name="")
static: set name of default fitter
virtual Double_t GetParameter(Int_t ipar) const =0
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
Double_t CpuTime()
Stop the stopwatch (if it is running) and return the cputime (in seconds) passed between the start an...
Definition: TStopwatch.cxx:123
void printMinimum(const std::vector< double > &x)
Definition: testMinim.cxx:402
int Int_t
Definition: RtypesCore.h:41
virtual void SetFCN(void *fcn)
To set the address of the minimization objective function.
int testRosenBrock()
Definition: testMinim.cxx:614
TLatex * t1
Definition: textangle.C:20
static double A[]
Short_t Abs(Short_t d)
Definition: TMathBase.h:110
double cos(double)
TFile * f
virtual void Clear()
reset for consecutive minimizations - implement if needed
Definition: Minimizer.h:125
TSocket * s1
Definition: hserv2.C:36
int DoOldMinimization(FCN func, TVirtualFitter *min, double &minval, double &edm)
Definition: testMinim.cxx:484
double sqrt(double)
void Stop()
Stop the stopwatch.
Definition: TStopwatch.cxx:75
Double_t x[n]
Definition: legend1.C:17
int testNewMinimizer(const ROOT::Math::IMultiGenFunction &func, const double *x0, const double *s0, const std::string &minimizer, const std::string &algoType)
Definition: testMinim.cxx:530
Abstract Minimizer class, defining the interface for the various minimizer (like Minuit2, Minuit, GSL, etc..) Plug-in'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
virtual void FdF(const double *x, double &f, double *df) const
Optimized method to evaluate at the same time the function value and derivative at a point x...
Definition: IFunction.h:358
int DoNewMinimization(const ROOT::Math::IMultiGenFunction &func, const double *x0, const double *s0, ROOT::Math::Minimizer *min, double &minval, double &edm, double *minx)
Definition: testMinim.cxx:416
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 ("").
double sin(double)
int main()
Definition: testMinim.cxx:740
double gAbsTolerance
Definition: testMinim.cxx:35
bool minos
Definition: testMinim.cxx:33
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
int testTrigoFletcher()
Definition: testMinim.cxx:651
ROOT::R::TRInterface & r
Definition: Object.C:4
int gNCall2
Definition: testMinim.cxx:28
SVector< double, 2 > v
Definition: Dict.h:5
int testChebyQuad()
Definition: testMinim.cxx:682
virtual unsigned int NDim() const =0
Retrieve the dimension of the function.
tuple w
Definition: qtexample.py:51
TLine * l
Definition: textangle.C:4
virtual Int_t ExecuteCommand(const char *command, Double_t *args, Int_t nargs)=0
void(* FCN)(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag)
Definition: testMinim.cxx:39
virtual Int_t GetStats(Double_t &amin, Double_t &edm, Double_t &errdef, Int_t &nvpar, Int_t &nparx) const =0
TText * t2
Definition: rootenv.C:28
Double_t Pi()
Definition: TMath.h:44
void SetMaxFunctionCalls(unsigned int maxfcn)
set maximum of function calls
Definition: Minimizer.h:456
TH1F * s2
Definition: threadsh2.C:15
static const double x1[5]
double Double_t
Definition: RtypesCore.h:55
void SetTolerance(double tol)
set the tolerance
Definition: Minimizer.h:462
const double * TrueXMinimum(const ROOT::Math::IMultiGenFunction &func)
Definition: testMinim.cxx:381
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
static Vc_ALWAYS_INLINE int_v max(const int_v &x, const int_v &y)
Definition: vector.h:440
Abstract Base Class for Fitting.
char Char_t
Definition: RtypesCore.h:29
typedef void((*Func_t)())
int gNmin
Definition: testMinim.cxx:29
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.
TRandom3 rndm
Definition: fit2dHist.C:44
std::string ToString(const T &val)
Utility function for conversion to strings.
Definition: Util.h:42
int gNCall
Definition: testMinim.cxx:27
virtual double Edm() const
return expected distance reached from the minimum (re-implement if minimizer provides it ...
Definition: Minimizer.h:262
TSocket * s0
Definition: hserv2.C:36
virtual unsigned int NCalls() const
number of function calls to reach the minimum
Definition: Minimizer.h:268
void SetPrintLevel(int level)
set print level
Definition: Minimizer.h:453
Documentation for the abstract class IBaseFunctionMultiDim.
Definition: IFunction.h:63
const Int_t n
Definition: legend1.C:16
virtual void Gradient(const double *x, double *grad) const
Evaluate all the vector of function derivatives (gradient) at a point x.
Definition: IFunction.h:342
virtual double DoEval(const double *x) const =0
Implementation of the evaluation function.
virtual IBaseFunctionMultiDim * Clone() const =0
Clone a function.
virtual double DoDerivative(const double *x, unsigned int icoord) const =0
function to evaluate the derivative with respect each coordinate.
Stopwatch class.
Definition: TStopwatch.h:30