File indexing completed on 2024-04-06 12:24:23
0001 #include "PhysicsTools/Utilities/interface/Integral.h"
0002 #include <iostream>
0003 #include <cmath>
0004 #include <vector>
0005 #include <utility>
0006 #include <algorithm>
0007 #include <sys/time.h>
0008 #include "TCanvas.h"
0009 #include "TGraph.h"
0010 #include "TAxis.h"
0011 #include "TROOT.h"
0012 #include "TH2.h"
0013
0014 using namespace funct;
0015 using namespace std;
0016
0017 double getTime() {
0018 struct timeval t;
0019 if (gettimeofday(&t, 0) < 0)
0020 abort();
0021 return (double)t.tv_sec + (double(t.tv_usec) * 1E-6);
0022 }
0023
0024 struct gauss {
0025 static const double c;
0026 double operator()(double x) const { return c * exp(-x * x); }
0027 };
0028
0029 const double gauss::c = 2. / sqrt(M_PI);
0030
0031 struct gauss1 : public gauss {};
0032
0033 struct gauss2 : public gauss {};
0034
0035 struct gauss3 : public gauss {};
0036
0037 struct gauss4 : public gauss {};
0038
0039 struct gaussPrimitive {
0040 double operator()(double x) const { return erf(x); }
0041 };
0042
0043 NUMERICAL_FUNCT_INTEGRAL(gauss1, TrapezoidIntegrator);
0044
0045 NUMERICAL_FUNCT_INTEGRAL(gauss2, GaussLegendreIntegrator);
0046
0047 NUMERICAL_FUNCT_INTEGRAL(gauss3, GaussIntegrator);
0048
0049 NUMERICAL_FUNCT_INTEGRAL(gauss4, RootIntegrator);
0050
0051 template <typename G, typename I>
0052 pair<double, double> check(const G &g, const I &i) {
0053 gaussPrimitive pr;
0054 double xMax = 1;
0055 double i0 = pr(xMax);
0056 double t0, t1;
0057 t0 = getTime();
0058 double i1 = integral_f(g, 0, xMax, i);
0059 t1 = getTime();
0060 pair<double, double> p = make_pair(t1 - t0, fabs(i1 - i0));
0061 cout << ">>> time: " << p.first << ", accuracy: " << p.second << endl;
0062 return p;
0063 }
0064
0065 int main() {
0066 gauss1 g1;
0067 gauss2 g2;
0068 gauss3 g3;
0069 gauss4 g4;
0070
0071 cout << ">>> Trapezoidal integration" << endl;
0072 vector<pair<double, double> > t;
0073 for (size_t i = 2; i < 20; i += 1) {
0074 TrapezoidIntegrator i1(i);
0075 t.push_back(check(g1, i1));
0076 }
0077 for (size_t i = 20; i < 1000; i += 20) {
0078 TrapezoidIntegrator i1(i);
0079 t.push_back(check(g1, i1));
0080 }
0081 cout << ">>> Gauss Legendre integration" << endl;
0082 vector<pair<double, double> > l;
0083 for (size_t i = 1; i < 20; i += 1) {
0084 GaussLegendreIntegrator i2(i, 1.e-5);
0085 l.push_back(check(g2, i2));
0086 }
0087 for (size_t i = 20; i < 1000; i += 20) {
0088 GaussLegendreIntegrator i2(i, 1.e-5);
0089 l.push_back(check(g2, i2));
0090 }
0091 cout << ">>> Gauss integration" << endl;
0092 vector<pair<double, double> > g;
0093 for (double e = 100; e > 1.e-5; e /= 2) {
0094 GaussIntegrator i3(e);
0095 g.push_back(check(g3, i3));
0096 }
0097 cout << ">>> ROOT GSL integration" << endl;
0098 vector<pair<double, double> > r;
0099 for (double e = 100; e > 1.e-5; e /= 2) {
0100 RootIntegrator i4(ROOT::Math::IntegrationOneDim::kADAPTIVESINGULAR, e, e);
0101 r.push_back(check(g4, i4));
0102 }
0103 gROOT->SetStyle("Plain");
0104 TCanvas canvas;
0105 canvas.SetLogx();
0106 canvas.SetLogy();
0107 double xMin = 1e6, xMax = 0, yMin = 1e6, yMax = 0;
0108 size_t nt = t.size();
0109 double *xt = new double[nt], *yt = new double[nt];
0110 const double xAbsMin = 1e-10, yAbsMin = 1e-10;
0111 for (size_t i = 0; i < nt; ++i) {
0112 xt[i] = t[i].first;
0113 yt[i] = t[i].second;
0114 if (xt[i] < xAbsMin)
0115 xt[i] = xAbsMin;
0116 if (yt[i] < yAbsMin)
0117 yt[i] = yAbsMin;
0118 if (xMin > xt[i])
0119 xMin = xt[i];
0120 if (xMax < xt[i])
0121 xMax = xt[i];
0122 if (yMin > yt[i])
0123 yMin = yt[i];
0124 if (yMax < yt[i])
0125 yMax = yt[i];
0126 }
0127 size_t nl = l.size();
0128 double *xl = new double[nl], *yl = new double[nl];
0129 for (size_t i = 0; i < nl; ++i) {
0130 xl[i] = l[i].first;
0131 yl[i] = l[i].second;
0132 if (xl[i] < xAbsMin)
0133 xl[i] = xAbsMin;
0134 if (yl[i] < yAbsMin)
0135 yl[i] = yAbsMin;
0136 if (xMin > xl[i])
0137 xMin = xl[i];
0138 if (xMax < xl[i])
0139 xMax = xl[i];
0140 if (yMin > yl[i])
0141 yMin = yl[i];
0142 if (yMax < yl[i])
0143 yMax = yl[i];
0144 }
0145 size_t ng = g.size();
0146 double *xg = new double[ng], *yg = new double[ng];
0147 for (size_t i = 0; i < ng; ++i) {
0148 xg[i] = g[i].first;
0149 yg[i] = g[i].second;
0150 if (xg[i] < xAbsMin)
0151 xg[i] = xAbsMin;
0152 if (yg[i] < yAbsMin)
0153 yg[i] = yAbsMin;
0154 if (xMin > xg[i])
0155 xMin = xg[i];
0156 if (xMax < xg[i])
0157 xMax = xg[i];
0158 if (yMin > yg[i])
0159 yMin = yg[i];
0160 if (yMax < yg[i])
0161 yMax = yg[i];
0162 }
0163 size_t nr = r.size();
0164 double *xr = new double[nr], *yr = new double[nr];
0165 for (size_t i = 0; i < nr; ++i) {
0166 xr[i] = r[i].first;
0167 yr[i] = r[i].second;
0168 if (xr[i] < xAbsMin)
0169 xr[i] = xAbsMin;
0170 if (yr[i] < yAbsMin)
0171 yr[i] = yAbsMin;
0172 if (xMin > xr[i])
0173 xMin = xr[i];
0174 if (xMax < xr[i])
0175 xMax = xr[i];
0176 if (yMin > yr[i])
0177 yMin = yr[i];
0178 if (yMax < yr[i])
0179 yMax = yr[i];
0180 }
0181 TH2F frame("frame", "Red: T, Blue: G-L, Green: G, Black: GSL", 1, xMin, xMax, 1, yMin, yMax);
0182 frame.GetXaxis()->SetTitle("CPU time (sec)");
0183 frame.GetYaxis()->SetTitle("Accuracy");
0184 frame.Draw();
0185 TGraph gt(nt, xt, yt), gl(nl, xl, yl), gg(ng, xg, yg), gr(nr, xr, yr);
0186 gt.SetMarkerStyle(21);
0187 gt.SetLineWidth(3);
0188 gt.SetLineColor(kRed);
0189 gl.SetMarkerStyle(21);
0190 gl.SetLineWidth(3);
0191 gl.SetLineColor(kBlue);
0192 gg.SetMarkerStyle(21);
0193 gg.SetLineWidth(3);
0194 gg.SetLineColor(kGreen);
0195 gr.SetMarkerStyle(21);
0196 gr.SetLineWidth(3);
0197 gr.SetLineColor(kBlack);
0198 gt.Draw("PC");
0199 gl.Draw("PC");
0200 gg.Draw("PC");
0201 gr.Draw("PC");
0202 canvas.SaveAs("integralTiming.eps");
0203 delete[] xt;
0204 delete[] yt;
0205 delete[] xl;
0206 delete[] yl;
0207 delete[] xg;
0208 delete[] yg;
0209 delete[] xr;
0210 delete[] yr;
0211
0212 return 0;
0213 }
