ROOT  6.06/09
Reference Guide
ParallelTest.cxx
Go to the documentation of this file.
1 // @(#)root/minuit2:$Id$
2 // Authors: M. Winkler, F. James, L. Moneta, A. Zsenei 2003-2005
3 
4 /**********************************************************************
5  * *
6  * Copyright (c) 2005 LCG ROOT Math team, CERN/PH-SFT *
7  * *
8  **********************************************************************/
9 
10 
11 #include "GaussRandomGen.h"
14 #include "Minuit2/MnPrint.h"
15 #include "Minuit2/MnMigrad.h"
16 #include "Minuit2/MnMinos.h"
17 #include "Minuit2/MnPlot.h"
18 #include "Minuit2/MinosError.h"
19 #include "Minuit2/FCNBase.h"
20 #include <cmath>
21 #include <iostream>
22 
23 // example of a multi dimensional fit where parallelization can be used
24 // to speed up the result
25 // define the environment variable OMP_NUM_THREADS to the number of desired threads
26 // By default it will have thenumber of core of the machine
27 // The default number of dimension is 20 (fit in 40 parameters) on 1000 data events.
28 // One can change the dimension and the number of events by doing:
29 // ./test_Minuit2_Parallel ndim nevents
30 
31 using namespace ROOT::Minuit2;
32 
33 const int default_ndim = 20;
34 const int default_ndata = 1000;
35 
36 
37 double GaussPdf(double x, double x0, double sigma) {
38  double tmp = (x-x0)/sigma;
39  return (1.0/(std::sqrt(2 * M_PI) * std::fabs(sigma))) * std::exp(-tmp*tmp/2);
40 }
41 
42 double LogMultiGaussPdf(const std::vector<double> & x, const std::vector<double> & p) {
43  double f = 0;
44  int ndim = x.size();
45  for (int k = 0; k < ndim; ++k) {
46  double y = GaussPdf( x[k], p[ 2*k], p[2*k + 1] );
47  //std::cout << " k " << k << " " << y << " " << p[2*k] << " " << p[2*k+1] << std::endl;
48  y = std::max(y, 1.E-300);
49  double lgaus = std::log( y );
50  f += lgaus;
51  }
52  return f;
53 }
54 
55 typedef std::vector< std::vector< double> > Data;
56 
57 struct LogLikeFCN : public FCNBase {
58 
59  LogLikeFCN( const Data & data) : fData(data) {}
60 
61  double operator() (const std::vector<double> & p ) const {
62  double logl = 0;
63  int ndata = fData.size();
64  for (int i = 0; i < ndata; ++i) {
65  logl -= LogMultiGaussPdf(fData[i], p);
66  }
67  return logl;
68  }
69  double Up() const { return 0.5; }
70  const Data & fData;
71 };
72 
73 int doFit(int ndim, int ndata) {
74 
75  // generate the data (1000 data points) in 100 dimension
76 
77  Data data(ndata);
78  std::vector< double> event(ndim);
79 
80  std::vector<double> mean(ndim);
81  std::vector<double> sigma(ndim);
82  for (int k = 0; k < ndim; ++k) {
83  mean[k] = -double(ndim)/2 + k;
84  sigma[k] = 1. + 0.1*k;
85  }
86 
87  // generate the data
88  for (int i = 0; i < ndata; ++i) {
89  for (int k = 0; k < ndim;++k) {
90  GaussRandomGen rgaus(mean[k],sigma[k] );
91  event[k] = rgaus();
92  }
93  data[i] = event;
94  }
95 
96  // create FCN function
97  LogLikeFCN fcn(data);
98 
99  // create initial starting values for parameters
100  std::vector<double> init_par(2*ndim);
101  for (int k = 0; k < ndim; ++k) {
102  init_par[2*k] = 0;
103  init_par[2*k+1] = 1;
104  }
105 
106  std::vector<double> init_err(2*ndim);
107  for (int k = 0; k < 2*ndim; ++k) {
108  init_err[k] = 0.1;
109  }
110  // create minimizer (default constructor)
111  VariableMetricMinimizer fMinimizer;
112 
113  // Minimize
114  FunctionMinimum min = fMinimizer.Minimize(fcn, init_par, init_err);
115 
116  // output
117  std::cout<<"minimum: "<<min<<std::endl;
118 
119 
120 // // create MINOS Error factory
121 // MnMinos Minos(fFCN, min);
122 
123 // {
124 // // 3-sigma MINOS errors (minimal interface)
125 // fFCN.SetErrorDef(9.);
126 // std::pair<double,double> e0 = Minos(0);
127 // std::pair<double,double> e1 = Minos(1);
128 // std::pair<double,double> e2 = Minos(2);
129 
130 
131 // // output
132 // std::cout<<"3-sigma Minos errors with limits: "<<std::endl;
133 // std::cout.precision(16);
134 // std::cout<<"par0: "<<min.UserState().Value("mean")<<" "<<e0.first<<" "<<e0.second<<std::endl;
135 // std::cout<<"par1: "<<min.UserState().Value(1)<<" "<<e1.first<<" "<<e1.second<<std::endl;
136 // std::cout<<"par2: "<<min.UserState().Value("area")<<" "<<e2.first<<" "<<e2.second<<std::endl;
137 
138 
139 // }
140 
141 // }
142 
143 
144  return 0;
145 }
146 
147 int main(int argc, char **argv) {
148  int ndim = default_ndim;
149  int ndata = default_ndata;
150  if (argc > 1) {
151  ndim = atoi(argv[1] );
152  }
153  if (argc > 2) {
154  ndata = atoi(argv[2] );
155  }
156  std::cout << "do fit of " << ndim << " dimensional data on " << ndata << " events " << std::endl;
157  doFit(ndim,ndata);
158 }
const int ndata
static Vc_ALWAYS_INLINE int_v min(const int_v &x, const int_v &y)
Definition: vector.h:433
const int default_ndim
int doFit(int ndim, int ndata)
double LogMultiGaussPdf(const std::vector< double > &x, const std::vector< double > &p)
double sqrt(double)
Double_t x[n]
Definition: legend1.C:17
class holding the full result of the minimization; both internal and external (MnUserParameterState) ...
std::vector< std::vector< double > > Data
Interface (abstract class) defining the function to be minimized, which has to be implemented by the ...
Definition: FCNBase.h:47
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
virtual FunctionMinimum Minimize(const FCNBase &, const std::vector< double > &, const std::vector< double > &, unsigned int stra=1, unsigned int maxfcn=0, double toler=0.1) const
#define M_PI
Definition: Rotated.cxx:105
const int default_ndata
Double_t E()
Definition: TMath.h:54
Instantiates the SeedGenerator and MinimumBuilder for Variable Metric Minimization method...
double GaussPdf(double x, double x0, double sigma)
TRObject operator()(const T1 &t1) const
double f(double x)
Double_t y[n]
Definition: legend1.C:17
int main(int argc, char **argv)
static Vc_ALWAYS_INLINE int_v max(const int_v &x, const int_v &y)
Definition: vector.h:440
double exp(double)
double log(double)