Jpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JMultiFunction1D.cc
Go to the documentation of this file.
1 #include <string>
2 #include <iostream>
3 #include <iomanip>
4 #include <cmath>
5 
6 #include "TRandom3.h"
7 #include "TMath.h"
8 
9 #include "JTools/JMapList.hh"
10 #include "JTools/JFunction1D_t.hh"
13 #include "JTools/JMultiFunction.hh"
14 #include "JTools/JQuantile.hh"
15 #include "JTools/JToolsToolkit.hh"
16 
17 #include "Jeep/JTimer.hh"
18 #include "Jeep/JParser.hh"
19 #include "Jeep/JMessage.hh"
20 
21 namespace {
22 
23  using namespace JPP;
24 
25  /**
26  * 5D test function.
27  */
28  inline double f5(const double x0,
29  const double x1,
30  const double x2,
31  const double x3,
32  const double x4)
33  {
34  return TMath::Gaus(x4, x0+x1+x2+x3, 1.0, kTRUE);
35  }
36 }
37 
38 
39 /**
40  * \file
41  *
42  * Example program to create 1D function interpolator from multi-dimensional interpolator.
43  * \author mdejong
44  */
45 int main(int argc, char **argv)
46 {
47  using namespace std;
48  using namespace JPP;
49 
50  int numberOfEvents;
51  int numberOfBins;
52  double precision;
53  int debug;
54 
55  try {
56 
57  JParser<> zap("Example program to create 1D function interpolator from multi-dimensional interpolator.");
58 
59  zap['n'] = make_field(numberOfEvents) = 1000;
60  zap['N'] = make_field(numberOfBins) = 11;
61  zap['e'] = make_field(precision) = 1.0e-14;
62  zap['d'] = make_field(debug) = 3;
63 
64  zap(argc, argv);
65  }
66  catch(const exception &error) {
67  FATAL(error.what() << endl);
68  }
69 
70  ASSERT(numberOfEvents > 0);
71 
72  const double xmin = -1.0;
73  const double xmax = +1.0;
74  const double dx = (xmax - xmin) / (numberOfBins - 1);
75 
76  typedef JSplineFunction1D_t JFunction1D_t;
77  typedef JFunction1D_t::abscissa_type abscissa_type;
78  typedef JFunction1D_t::value_type value_type;
79 
81 
85  JPolint1FunctionalMap>::maplist JMaplist_t;
86 
89 
90  for (double x0 = xmin; x0 < xmax + 0.5*dx; x0 += dx) {
91  for (double x1 = xmin; x1 < xmax + 0.5*dx; x1 += dx) {
92  for (double x2 = xmin; x2 < xmax + 0.5*dx; x2 += dx) {
93  for (double x3 = xmin; x3 < xmax + 0.5*dx; x3 += dx) {
94  for (double x4 = xmin; x4 < xmax + 0.5*dx; x4 += dx) {
95  g5[x0][x1][x2][x3][x4] = f5(x0,x1,x2,x3,x4);
96  h5[x0][x1][x2][x3][x4] = f5(x0,x1,x2,x3,x4);
97  }
98  }
99  }
100  }
101  }
102 
103  g5.compile();
104 
105 
106  const double x0 = +0.15;
107  const double x1 = -0.25;
108  const double x2 = 0.25;
109  const double x3 = -0.15;
110 
111  JFunction1D_t g1;
112 
113  JTimer timer("4D interpolator");
114 
115  timer.start();
116 
117  for (int i = 0; i != 100; ++i) {
118  copy(h5(x0,x1,x2,x3), g1); // put interpolated data into 1D functional object
119  }
120 
121  timer.stop();
122  timer.print(cout, 1.0/100, micro_t);
123 
124  g1.compile();
125 
126  JTimer t1("1D interpolator");
127  JTimer t5("5D interpolator");
128 
129  JQuantile Q1("1D interpolator");
130  JQuantile Q5("5D interpolator");
131 
132  // mini buffer to reduce overhead in timer
133 
134  const int N = 100;
135 
136  double x[N];
137  double y[N];
138  double v[N];
139  double w[N];
140 
141  for (int i = 0; i != numberOfEvents; ++i) {
142 
143  for (int __i = 0; __i != N; ++__i) {
144 
145  const double x4 = gRandom->Uniform(xmin, xmax);
146 
147  x[__i] = x4;
148  y[__i] = f5(x0,x1,x2,x3,x4);
149  }
150 
151  t1.start();
152 
153  for (int __i = 0; __i != N; ++__i) {
154  v[__i] = g1(x[__i]);
155  }
156 
157  t1.stop();
158 
159  t5.start();
160 
161  for (int __i = 0; __i != N; ++__i) {
162  w[__i] = g5(x0,x1,x2,x3,x[__i]);
163  }
164 
165  t5.stop();
166 
167  for (int __i = 0; __i != N; ++__i) {
168  Q1.put(v[__i] - y[__i]);
169  Q5.put(w[__i] - y[__i]);
170  }
171  }
172 
173  if (debug >= notice_t) {
174 
175  Q1.print(cout);
176  Q5.print(cout, false);
177 
178  t1.print(cout, 1.0 / (N * numberOfEvents), nano_t);
179  t5.print(cout, 1.0 / (N * numberOfEvents), nano_t);
180  }
181 
182  ASSERT(fabs(Q1.getMean() - Q5.getMean()) <= precision);
183  ASSERT(fabs(Q1.getSTDev() - Q5.getSTDev()) <= precision);
184 
185  return 0;
186 }
Utility class to parse command line options.
Definition: JParser.hh:1410
micro
Definition: JScale.hh:29
void put(const double x, const double w=1.0)
Put value.
Definition: JQuantile.hh:77
void print(std::ostream &out, const JScale_t scale=milli_t) const
Print timer data.
Definition: JTimer.hh:161
notice
Definition: JMessage.hh:30
Quantile calculator.
Definition: JQuantile.hh:34
void stop()
Stop timer.
Definition: JTimer.hh:113
Type definition of a spline interpolation method based on a JCollection with double result type...
Various implementations of functional maps.
Template implementation of function object in one dimension returning a constant value.
#define ASSERT(A)
Assert macro.
Definition: JMessage.hh:72
double getSTDev() const
Get standard deviation.
Definition: JQuantile.hh:159
Multidimensional interpolation method.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1836
Type definition of a 1st degree polynomial interpolation based on a JMap implementation.
Auxiliary class for CPU timing and usage.
Definition: JTimer.hh:32
This include file contains various recursive methods to operate on multi-dimensional collections...
int debug
debug level
Definition: JSirene.cc:59
std::ostream & print(std::ostream &out, bool lpr=true) const
Print quantile.
Definition: JQuantile.hh:228
Auxiliary class for recursive map list generation.
Definition: JMapList.hh:108
General purpose messaging.
#define FATAL(A)
Definition: JMessage.hh:65
double getMean() const
Get mean value.
Definition: JQuantile.hh:131
nano
Definition: JScale.hh:28
Utility class to parse command line options.
double f5(const double x0, const double x1, const double x2, const double x3, const double x4)
5D function.
void copy(const Head &from, JHead &to)
Copy header from from to to.
Definition: JHead.cc:40
void start()
Start timer.
Definition: JTimer.hh:89
Double_t g1(const Double_t x)
Function.
Definition: JQuantiles.cc:25
int main(int argc, char *argv[])
Definition: Main.cpp:15