Jpp 19.3.0-rc.3
the software that should make you happy
Loading...
Searching...
No Matches
JPolynome1D.cc
Go to the documentation of this file.
1#include <string>
2#include <iostream>
3#include <sstream>
4#include <iomanip>
5#include <map>
6
7#include "TRandom3.h"
8
9#include "JLang/JException.hh"
10#include "JMath/JPolynome.hh"
12#include "JTools/JGrid.hh"
13#include "JTools/JQuantile.hh"
15
16#include "Jeep/JPrint.hh"
17#include "Jeep/JParser.hh"
18#include "Jeep/JMessage.hh"
19
20namespace {
21
23
24 struct JABC {
25
26 JABC(const char* title) :
27 f (title),
28 fp(title),
29 v (title)
30 {}
31
32 JQuantile f;
33 JQuantile fp;
34 JQuantile v;
35 };
36
37
38 struct JPrecision_t {
39
40 JPrecision_t() :
41 f (0.0),
42 fp(0.0),
43 v (0.0)
44 {}
45
46 JPrecision_t(const double f,
47 const double fp,
48 const double v) :
49 f (f),
50 fp(fp),
51 v (v)
52 {}
53
54 friend std::istream& operator>>(std::istream& in, JPrecision_t& object)
55 {
56 return in >> object.f >> object.fp >> object.v;
57 }
58
59 friend std::ostream& operator<<(std::ostream& out, const JPrecision_t& object)
60 {
61 return out << object.f << ' ' << object.fp << ' ' << object.v;
62 }
63
64 double f;
65 double fp;
66 double v;
67 };
68}
69
70
71/**
72 * \file
73 *
74 * Example program to test 1D interpolation of a polynome.
75 * \author mdejong
76 */
77int main(int argc, char **argv)
78{
79 using namespace std;
80 using namespace JPP;
81
82 typedef map<string, JPrecision_t> map_type;
83
84 unsigned int numberOfEvents;
85 unsigned int numberOfBins;
86 JPolynome f1;
87 map_type precision;
88 int debug;
89
90 precision["spline"] = JPrecision_t(1.0e-3, 1.0e-1, 1.0e-4);
91 precision["hermite"] = JPrecision_t(1.0e-3, 1.0e-1, 1.0e-4);
92 precision["grid"] = JPrecision_t(1.0e-3, 1.0e-1, 1.0e-4);
93 precision["polint"] = JPrecision_t(1.0e-14, 1.0e-12, 1.0e-14);
94
95 try {
96
97 JParser<> zap("Example program to test 1D interpolation of a polynome.");
98
99 zap['n'] = make_field(numberOfEvents) = 0;
100 zap['N'] = make_field(numberOfBins) = 21;
101 zap['P'] = make_field(f1);
102 zap['e'] = make_field(precision) = JPARSER::initialised();
103 zap['d'] = make_field(debug) = 3;
104
105 zap(argc, argv);
106 }
107 catch(const exception &error) {
108 FATAL(error.what() << endl);
109 }
110
111
112 const int N = 3; // fixed degree of polynomial interpolation
113
114 JSplineFunction1S_t spline;
118 JPolintFunction1D_t<N+1> Polint; // integrand of polint
119
120 const double xmin = -1.0;
121 const double xmax = +1.0;
122
123 spline .configure(make_grid(numberOfBins, xmin, xmax), f1);
124 hermite.configure(make_grid(numberOfBins, xmin, xmax), f1);
125 grid .configure(make_grid(numberOfBins, xmin, xmax), f1);
126 polint .configure(make_grid(numberOfBins, xmin, xmax), f1);
127
128 spline .compile();
129 hermite.compile();
130 grid .compile();
131 polint .compile();
132
133 integrate(polint, Polint);
134
135 Polint.compile();
136
137 if (numberOfEvents != 0) {
138
139 JABC Qspline ("spline");
140 JABC Qhermite("hermite");
141 JABC Qgrid ("grid");
142 JABC Qpolint ("polint");
143
144 for (unsigned int i = 0; i != numberOfEvents; ++i) {
145
146 const double x = gRandom->Uniform(xmin, xmax);
147
148 const double y1 = f1.getValue(x);
149 const double fp = f1.getDerivative(x);
150 const double F1 = f1.getIntegral(x) - f1.getIntegral(xmin);
151
152 Qspline .f .put(y1 - spline (x).f);
153 Qhermite.f .put(y1 - hermite(x).f);
154 Qgrid .f .put(y1 - grid (x).f);
155 Qpolint .f .put(y1 - polint (x).f);
156
157 Qspline .fp.put(fp - spline (x).fp);
158 Qhermite.fp.put(fp - hermite(x).fp);
159 Qgrid .fp.put(fp - grid (x).fp);
160 Qpolint .fp.put(fp - polint (x).fp);
161
162 Qspline .v .put(F1 - spline (x).v);
163 Qhermite.v .put(F1 - hermite(x).v);
164 Qgrid .v .put(F1 - grid (x).v);
165 Qpolint .v .put(F1 - Polint (x));
166 }
167
168 if (debug >= debug_t) {
169
170 cout << "RMS: f f' F " << endl;
171 cout << "_________________________________________" << endl;
172
173 cout << "spline "
174 << SCIENTIFIC(10,3) << Qspline .f .getSTDev() << ' '
175 << SCIENTIFIC(10,3) << Qspline .fp.getSTDev() << ' '
176 << SCIENTIFIC(10,3) << Qspline .v .getSTDev() << endl;
177
178 cout << "hermite "
179 << SCIENTIFIC(10,3) << Qhermite.f .getSTDev() << ' '
180 << SCIENTIFIC(10,3) << Qhermite.fp.getSTDev() << ' '
181 << SCIENTIFIC(10,3) << Qhermite.v .getSTDev() << endl;
182
183 cout << "grid "
184 << SCIENTIFIC(10,3) << Qgrid .f .getSTDev() << ' '
185 << SCIENTIFIC(10,3) << Qgrid .fp.getSTDev() << ' '
186 << SCIENTIFIC(10,3) << Qgrid .v .getSTDev() << endl;
187
188 cout << "polint "
189 << SCIENTIFIC(10,3) << Qpolint .f .getSTDev() << ' '
190 << SCIENTIFIC(10,3) << Qpolint .fp.getSTDev() << ' '
191 << SCIENTIFIC(10,3) << Qpolint .v .getSTDev() << endl;
192 }
193
194 ASSERT(Qspline .f .getSTDev() <= precision["spline"].f);
195 ASSERT(Qspline .fp.getSTDev() <= precision["spline"].fp);
196 ASSERT(Qspline .v .getSTDev() <= precision["spline"].v);
197
198 ASSERT(Qhermite.f .getSTDev() <= precision["hermite"].f);
199 ASSERT(Qhermite.fp.getSTDev() <= precision["hermite"].fp);
200 ASSERT(Qhermite.v .getSTDev() <= precision["hermite"].v);
201
202 ASSERT(Qgrid .f .getSTDev() <= precision["grid"].f);
203 ASSERT(Qgrid .fp.getSTDev() <= precision["grid"].fp);
204 ASSERT(Qgrid .v .getSTDev() <= precision["grid"].v);
205
206 ASSERT(Qpolint .f .getSTDev() <= precision["polint"].f);
207 ASSERT(Qpolint .fp.getSTDev() <= precision["polint"].fp);
208 ASSERT(Qpolint .v .getSTDev() <= precision["polint"].v);
209
210 } else {
211
212 for ( ; ; ) {
213
214 string buffer;
215
216 cout << "> " << flush;
217
218 getline(cin, buffer);
219
220 if (buffer != "") {
221
222 try {
223
224 double x;
225
226 istringstream(buffer) >> x;
227
228 cout << "f1 " << f1 (x) << endl;
229 cout << "spline " << spline (x).f << endl;
230 cout << "hermite " << hermite(x).f << endl;
231 cout << "grid spline " << grid (x).f << endl;
232 cout << "polynomial " << polint (x).f << endl;
233 }
234 catch(const JException& exception) {
235 cout << exception.what() << endl;
236 }
237
238 } else {
239 break;
240 }
241 }
242 }
243
244 return 0;
245}
Exceptions.
std::istream & operator>>(std::istream &in, JAANET::JHead &header)
Read header from input.
Definition JHead.hh:1832
General purpose messaging.
#define ASSERT(A,...)
Assert macro.
Definition JMessage.hh:90
#define FATAL(A)
Definition JMessage.hh:67
int debug
debug level
Definition JSirene.cc:72
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
int main(int argc, char **argv)
I/O formatting auxiliaries.
int numberOfBins
number of bins for average CDF integral of optical module
Definition JSirene.cc:73
This include file contains various recursive methods to operate on multi-dimensional collections.
General exception.
Definition JException.hh:24
Utility class to parse command line options.
Definition JParser.hh:1698
std::ostream & operator<<(std::ostream &stream, const CLBCommonHeader &header)
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Recursive template class for polynomial function.
Definition JPolynome.hh:165
double getIntegral(const double x) const
Integral value.
Definition JPolynome.hh:276
double getValue(const double x) const
Function value.
Definition JMathlib.hh:1421
double getDerivative(const double x) const
Derivative value.
Definition JMathlib.hh:1433
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition JParser.hh:68
Type definition of a spline interpolation based on a JGridCollection with JResultPDF result type.
Type definition of a spline interpolation method based on a JCollection with JResultPDF result type.
Polynomial interpolation method with result type double.
Polynomial interpolation method with result type JResultDerivative.
Auxiliary data structure for running average, standard deviation and quantiles.
Definition JQuantile.hh:46
Type definition of a spline interpolation method based on a JCollection with JResultPDF result type.
Auxiliary data structure for floating point format specification.
Definition JManip.hh:488