Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:33:33

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 }