Jpp 20.0.0-rc.2
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/JStats.hh"
15
16#include "Jeep/JPrint.hh"
17#include "Jeep/JParser.hh"
18#include "Jeep/JMessage.hh"
19
20namespace {
21
22 using JTOOLS::JStats;
23
24 struct JABC {
25
26 JABC(const char* title) :
27 f (title),
28 fp(title),
29 v (title)
30 {}
31
32 JStats f;
33 JStats fp;
34 JStats 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:1850
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.
Type definition of a spline interpolation method based on a JCollection with JResultPDF result type.
Auxiliary data structure for running average and standard deviation.
Definition JStats.hh:44
Auxiliary data structure for floating point format specification.
Definition JManip.hh:488