From 3e4f3137bf6f20ad19171dc6d5d93fb7fd47b2b9 Mon Sep 17 00:00:00 2001
From: Frank Sauerburger <f.sauerburger@cern.ch>
Date: Thu, 1 Aug 2019 12:57:09 +0200
Subject: [PATCH] Implement PDFs
---
nnfwtbn/tests/test_toydata.py | 80 ++++++++++++++++++++++++++++++++++-
nnfwtbn/toydata.py | 36 ++++++++++++++++
2 files changed, 115 insertions(+), 1 deletion(-)
diff --git a/nnfwtbn/tests/test_toydata.py b/nnfwtbn/tests/test_toydata.py
index a67fecb..4ffa540 100644
--- a/nnfwtbn/tests/test_toydata.py
+++ b/nnfwtbn/tests/test_toydata.py
@@ -1,7 +1,7 @@
import unittest
import numpy as np
-from nnfwtbn.toydata import rand
+from nnfwtbn.toydata import rand, ipdf_gluon, ipdf_quark
class ToyDataTestBase(unittest.TestCase):
@@ -51,3 +51,81 @@ class ToyDataTestBase(unittest.TestCase):
self.assertAlmostEqual(a, 0.90141859)
self.assertAlmostEqual(b, 0.85225178)
self.assertAlmostEqual(c, 0.93632300)
+
+ def test_pdf_gluon_pos(self):
+ """
+ Check that the pdf is always positive.
+ """
+ N = 100
+ for i in range(N):
+ self.assertTrue(ipdf_gluon(i / N, 1 / N) >= 0)
+
+
+ def test_pdf_gluon_normalized(self):
+ """
+ Check that the pdf is normalized.
+ """
+ self.assertAlmostEqual(ipdf_gluon(0, 1), 1)
+
+ def test_pdf_gluon_max(self):
+ """
+ Check that the pdf has it's maximum at x=0.
+ """
+ N = 100
+ bins = np.empty(N)
+ for i in range(N):
+ bins[i] = ipdf_gluon(i / N, 1 / N)
+
+ self.assertEqual(np.argmax(bins), 0)
+
+
+ def test_pdf_gluon_at_1(self):
+ """
+ Check that the pdf vanishes at x=1.
+ """
+ d = 1e-4
+ self.assertAlmostEqual(ipdf_gluon(1 - d, d), 0)
+
+ def test_pdf_quark_pos(self):
+ """
+ Check that the pdf is always positive.
+ """
+ N = 100
+ for i in range(N):
+ self.assertTrue(ipdf_quark(i / N, 1 / N) >= 0)
+
+
+ def test_pdf_quark_normalized(self):
+ """
+ Check that the pdf is normalized.
+ """
+ self.assertAlmostEqual(ipdf_quark(0, 1), 1)
+
+ def test_pdf_quark_max(self):
+ """
+ Check that the pdf has it's maximum at x=0.
+ """
+ N = 100
+ bins = np.empty(N)
+ for i in range(N):
+ bins[i] = ipdf_quark(i / N, 1 / N)
+
+ self.assertEqual(np.argmax(bins), 0)
+
+
+ def test_pdf_quark_at_1(self):
+ """
+ Check that the pdf vanishes at x=1.
+ """
+ d = 1e-9
+ self.assertAlmostEqual(ipdf_quark(1 - d, d), 0)
+
+ def test_pdf_quark_peak(self):
+ """
+ Check that there is a peak at x=0.5.
+ """
+ d = 1e-4
+ self.assertGreater(
+ ipdf_quark(0.5, d),
+ ipdf_quark(0.3, d)
+ )
diff --git a/nnfwtbn/toydata.py b/nnfwtbn/toydata.py
index 2dffb59..068e48f 100644
--- a/nnfwtbn/toydata.py
+++ b/nnfwtbn/toydata.py
@@ -5,7 +5,9 @@ module does not rely on external random number generators (seeding numpy
might break user code).
"""
+import math
import numpy as np
+from scipy.integrate import quad
from numba import njit
@njit
@@ -31,3 +33,37 @@ def rand(size=1, seed=1991):
output[i] = (a * output[i - 1] + c) % m
return output / m
+
+def ipdf_gluon(x, dx):
+ """
+ Grossly simplification of the integrated parton distribution function of a
+ proton for gluons. The method returns the probability to find a gluon
+ within the proton carrying a momentum fraction between x and x + dx.
+
+ The pdf is normalized to unity.
+ >>> '%.4f' % ipdf_gluon(0, 1)
+ '1.0000'
+ """
+ def pdf(x):
+ return 7 * (x - 1)**6
+
+ return quad(pdf, x, x + dx)[0]
+
+def ipdf_quark(x, dx):
+ """
+ Grossly simplification of the integrated parton distribution function of a
+ proton for quarks (there is only a single quark flavor). The method
+ returns the probability to find a gluon within the proton carrying a
+ momentum fraction between x and x + dx.
+
+ The quark pdf has a bump at x=0.5.
+
+ The pdf is normalized to unity.
+ >>> '%.4f' % ipdf_quark(0, 1)
+ '1.0000'
+ """
+ def pdf(x):
+ return 0.6 * 7 * (x - 1)**6 \
+ + 0.4 / (20 * math.atan(5) * ((x - 0.5)**2 + 0.01))
+
+ return quad(pdf, x, x + dx)[0]
--
GitLab