Jpp  16.0.0-rc.2
the software that should make you happy
 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 
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  h5.compile();
105 
106 
107  const double x0 = +0.15;
108  const double x1 = -0.25;
109  const double x2 = 0.25;
110  const double x3 = -0.15;
111 
112  JFunction1D_t g1;
113 
114  JTimer timer("4D interpolator");
115 
116  timer.start();
117 
118  for (int i = 0; i != 100; ++i) {
119  copy(h5(x0,x1,x2,x3), g1); // put interpolated data into 1D functional object
120  }
121 
122  timer.stop();
123  timer.print(cout, 1.0/100, micro_t);
124 
125  g1.compile();
126 
127  JTimer t1("1D interpolator");
128  JTimer t5("5D interpolator");
129 
130  JQuantile Q1("1D interpolator");
131  JQuantile Q5("5D interpolator");
132 
133  // mini buffer to reduce overhead in timer
134 
135  const int N = 100;
136 
137  double x[N];
138  double y[N];
139  double v[N];
140  double w[N];
141 
142  for (int i = 0; i != numberOfEvents; ++i) {
143 
144  for (int __i = 0; __i != N; ++__i) {
145 
146  const double x4 = gRandom->Uniform(xmin, xmax);
147 
148  x[__i] = x4;
149  y[__i] = f5(x0,x1,x2,x3,x4);
150  }
151 
152  t1.start();
153 
154  for (int __i = 0; __i != N; ++__i) {
155  v[__i] = g1(x[__i]);
156  }
157 
158  t1.stop();
159 
160  t5.start();
161 
162  for (int __i = 0; __i != N; ++__i) {
163  w[__i] = g5(x0,x1,x2,x3,x[__i]);
164  }
165 
166  t5.stop();
167 
168  for (int __i = 0; __i != N; ++__i) {
169  Q1.put(v[__i] - y[__i]);
170  Q5.put(w[__i] - y[__i]);
171  }
172  }
173 
174  if (debug >= notice_t) {
175 
176  Q1.print(cout);
177  Q5.print(cout, false);
178 
179  t1.print(cout, 1.0 / (N * numberOfEvents), nano_t);
180  t5.print(cout, 1.0 / (N * numberOfEvents), nano_t);
181  }
182 
183  ASSERT(fabs(Q1.getMean() - Q5.getMean()) <= precision);
184  ASSERT(fabs(Q1.getSTDev() - Q5.getSTDev()) <= precision);
185 
186  return 0;
187 }
Utility class to parse command line options.
Definition: JParser.hh:1500
data_type w[N+1][M+1]
Definition: JPolint.hh:757
int main(int argc, char *argv[])
Definition: Main.cc:15
General purpose class for collection of elements, see: &lt;a href=&quot;JTools.PDF&quot;;&gt;Collection of elements...
Definition: JCollection.hh:73
double getSTDev() const
Get standard deviation.
Definition: JQuantile.hh:266
micro
Definition: JScale.hh:29
Auxiliary data structure for running average, standard deviation and quantiles.
Definition: JQuantile.hh:43
void print(std::ostream &out, const JScale_t scale=milli_t) const
Print timer data.
Definition: JTimer.hh:161
notice
Definition: JMessage.hh:32
then JShowerPostfit f $INPUT_FILE o $OUTPUT_FILE N
void stop()
Stop timer.
Definition: JTimer.hh:113
double getMean() const
Get mean value.
Definition: JQuantile.hh:252
Various implementations of functional maps.
Template class for spline interpolation in 1D.
Definition: JSpline.hh:666
Template implementation of function object in one dimension returning a constant value.
#define ASSERT(A,...)
Assert macro.
Definition: JMessage.hh:90
Multidimensional interpolation method.
std::ostream & print(std::ostream &out, bool lpr=true) const
Print quantile.
Definition: JQuantile.hh:335
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1961
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...
void put(const double x, const double w=1.0)
Put value.
Definition: JQuantile.hh:133
int debug
debug level
Definition: JSirene.cc:63
Auxiliary class for recursive map list generation.
Definition: JMapList.hh:108
General purpose messaging.
#define FATAL(A)
Definition: JMessage.hh:67
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:139
data_type v[N+1][M+1]
Definition: JPolint.hh:756
int numberOfBins
number of bins for average CDF integral of optical module
Definition: JSirene.cc:64
void start()
Start timer.
Definition: JTimer.hh:89
Double_t g1(const Double_t x)
Function.
Definition: JQuantiles.cc:25