Merge pull request #47 from N3PDF/fhe

my branch

.github/workflows/pythonapp.yml

@@ -6,13 +6,18 @@ jobs:
   build:
 
     runs-on: ubuntu-latest
+    strategy:
+#      max-parallel: 2
+      matrix:
+        python-version: [3.7, 3.8]
+      fail-fast: false
 
     steps:
     - uses: actions/checkout@v2
-    - name: Set up Python 3.8
+    - name: Set up Python ${{ matrix.python-version }} 🐍
       uses: actions/setup-python@v1
       with:
-        python-version: 3.8
+        python-version: ${{ matrix.python-version }}
     - name: Install dependencies
       run: |
         sudo apt-get install -y libgsl-dev
@@ -28,4 +33,7 @@ jobs:
     - name: Test with pytest
       run: |
         pip install pytest
+        # until https://github.com/numba/numba/pull/5660 is confirmed we need to deactivate numba prior running
+        export NUMBA_DISABLE_JIT=1
+        pip install pytest-cov
         pytest

.pylintrc

@@ -3,7 +3,7 @@
 # A comma-separated list of package or module names from where C extensions may
 # be loaded. Extensions are loading into the active Python interpreter and may
 # run arbitrary code
-extension-pkg-whitelist=numpy,lhapdf,_gsl_digamma
+extension-pkg-whitelist=numpy,lhapdf,_gsl_digamma,numba
 
 # Add files or directories to the blacklist. They should be base names, not
 # paths.
@@ -392,7 +392,7 @@ ignored-classes=optparse.Values,thread._local,_thread._local
 # (useful for modules/projects where namespaces are manipulated during runtime
 # and thus existing member attributes cannot be deduced by static analysis. It
 # supports qualified module names, as well as Unix pattern matching.
-ignored-modules=matplotlib.cm
+ignored-modules=matplotlib.cm,numba
 
 # Show a hint with possible names when a member name was not found. The aspect
 # of finding the hint is based on edit distance.

README.md

@@ -1,12 +1,12 @@
 <p align="center">
-  <a href="https://github.com/N3PDF/eko/actions?query=workflow%3A%22eko%22"><img alt="Tests" src="https://github.com/N3PDF/eko/workflows/eko/badge.svg" /></a> <a href="https://n3pdf.github.io/eko/"><img alt="Docs" src="https://github.com/N3PDF/eko/workflows/docs/badge.svg"></a>
+  <a href="https://n3pdf.github.io/eko/"><img alt="EKO" src="doc/source/img/Logo.svg" width=300></a>
 </p>
+<p align="center">
+  <a href="https://github.com/N3PDF/eko/actions?query=workflow%3A%22eko%22"><img alt="Tests" src="https://github.com/N3PDF/eko/workflows/eko/badge.svg" /></a> <a href="https://n3pdf.github.io/eko/"><img alt="Docs" src="https://github.com/N3PDF/eko/workflows/docs/badge.svg"></a> <a href="https://doi.org/10.5281/zenodo.3874238"><img src="https://zenodo.org/badge/DOI/10.5281/zenodo.3874238.svg" alt="DOI"></a>
 </p>
 
 # EKO
-Evolution Kernel Operator
---
-Library for writing Evolution Kernel Operators to file in x-space.
+Python module to solve DGLAP equations in term of Evolution Kernel Operators in x-space.
 
 ## Installation
 ```
@@ -15,6 +15,3 @@ python setup.py install
 
 ## Documentation
 The documentation style of this code follows closely the [numpy documentation guide](https://numpydoc.readthedocs.io/en/latest/format.html).
-
-## Scope of the project
-The unique scope and aim of this code is (see above)

benchmarks/LHA-benchmark-paper.py

@@ -11,7 +11,6 @@
 from matplotlib.backends.backend_pdf import PdfPages
 
 from eko import run_dglap
-import eko.interpolation as interpolation
 
 from tools import plot_dist, save_all_operators_to_pdf, merge_dicts
 
@@ -205,12 +204,11 @@ def __init__(self, xgrid, polynomial_degree, flag):
             "PTO": 0,
             "alphas": 0.35,
             "Qref": np.sqrt(2),
-            "Qmc": np.sqrt(self._Q2init),
-            "Qmb": 4.5,
-            "Qmt": 175.0,
-            "xgrid_type": "custom",
-            "xgrid": xgrid,
-            "xgrid_polynom_rank": polynom_rank,
+            "mc": np.sqrt(self._Q2init),
+            "mb": 4.5,
+            "mt": 175.0,
+            "interpolation_xgrid": xgrid,
+            "interpolation_polynomial_degree": polynom_rank,
         }
 
     def _post_process(self, output, ref, tag):
@@ -237,7 +235,7 @@ def _post_process(self, output, ref, tag):
             save_all_operators_to_pdf(first_ops, assets_path / f"LHA-LO-{tag}-ops-{flag}.pdf")
         # dump operators to file
         if self.post_process_config["write_operator"]:
-            output.write_YAML_to_file(assets_path / f"LHA-LO-{tag}-ops-{flag}.yaml")
+            output.dump_yaml_to_file(assets_path / f"LHA-LO-{tag}-ops-{flag}.yaml")
 
     def _run_FFNS_raw(self, tag, Q2init, Q2final):
         """
@@ -295,19 +293,19 @@ def run_FFNS_twostep(self, Q2mid):
                 ret2t1 : dict
                     DGLAP result
         """
-        raise Exception("JCM: This cannot be done now")
-        # # suppress PDF plots for single steps
-        # plot_PDF = self.post_process_config["plot_PDF"]
-        # self.post_process_config["plot_PDF"] = False
-        # # step 1
-        # step1 = self._run_FFNS_raw("FFNS-twostep-step1", self._Q2init, Q2mid)
-        # # step 2
-        # step2 = self._run_FFNS_raw("FFNS-twostep-step2", Q2mid, self._Q2final)
-        # # join 2*1
-        # self.post_process_config["plot_PDF"] = plot_PDF
-        # ret2t1 = dglap.multiply_operators(step2, step1)
-        # self._post_process(ret2t1, LHA_final_dict_FFNS_ref, "FFNS-twostep-step2t1")
-        # return ret2t1
+        #raise Exception("JCM: This cannot be done now")
+        # suppress PDF plots for single steps
+        plot_PDF = self.post_process_config["plot_PDF"]
+        self.post_process_config["plot_PDF"] = False
+        # step 1
+        step1 = self._run_FFNS_raw("FFNS-twostep-step1", self._Q2init, Q2mid)
+        # step 2
+        step2 = self._run_FFNS_raw("FFNS-twostep-step2", Q2mid, self._Q2final)
+        # join 2*1
+        self.post_process_config["plot_PDF"] = plot_PDF
+        ret2t1 = step2.concat(step1)
+        self._post_process(ret2t1, LHA_final_dict_FFNS_ref, "FFNS-twostep-step2t1")
+        return ret2t1
 
     def run_ZMVFNS(self):
         """
@@ -346,7 +344,7 @@ def _save_final_scale_plots_to_pdf(self, path, output, ref):
         """
         pp = PdfPages(path)
         # get
-        pdf_grid = output.apply_PDF(LHA_init_pdfs, toy_xgrid)
+        pdf_grid = output.apply_pdf(LHA_init_pdfs, toy_xgrid)
         first_res = list(pdf_grid.values())[0]
         my_pdfs = first_res["pdfs"]
         my_pdf_errs = first_res["errors"]
@@ -412,10 +410,8 @@ def save_initial_scale_plots_to_pdf(path):
 
     # combine grid
     flag = f"l{n_low}m{n_mid}r{polynom_rank}"
-    xgrid_low = interpolation.get_xgrid_linear_at_log(
-        n_low, 1e-7, 1.0 if n_mid == 0 else 0.1
-    )
-    xgrid_mid = interpolation.get_xgrid_linear_at_id(n_mid, 0.1, 1.0)
+    xgrid_low = np.geomspace(1e-7, 1.0 if n_mid == 0 else 0.1,n_low)
+    xgrid_mid = np.linspace(0.1, 1.0, n_mid)
     xgrid_high = np.array([])
     xgrid = np.unique(np.concatenate((xgrid_low, xgrid_mid, xgrid_high)))
 
@@ -432,8 +428,8 @@ def save_initial_scale_plots_to_pdf(path):
     # check input scale
     #save_initial_scale_plots_to_pdf(assets_path / f"LHA-LO-FFNS-init-{flag}.pdf")
     # check fixed flavours
-    app.run_FFNS()
+    #app.run_FFNS()
     # do two steps
-    #app.run_FFNS_twostep(1e2)
+    app.run_FFNS_twostep(1e2)
     # check ZM-VFNS
-    # app.run_ZMVFNS()
+    #app.run_ZMVFNS()

benchmarks/LHAPDF-benchmarks.py

@@ -70,9 +70,8 @@ def __init__(self, pdfname, pdfmember, xgrid, polynomial_degree, flag):
             "Qref": float(self._pdfset.get_entry("MZ")),
             "FNS": "FFNS",
             "NfFF": int(self._pdfset.get_entry("NumFlavors")),
-            "xgrid_type": "custom",
-            "xgrid": xgrid,
-            "xgrid_polynom_rank": polynom_rank,
+            "interpolation_xgrid": xgrid,
+            "interpolation_polynomial_degree": polynom_rank,
         }
         self._Q2init = None
         self._Q2final = None
@@ -99,7 +98,7 @@ def _post_process(self, ret):
                 / f"LHAPDF-{self._pdfname}-{self._pdfmember}-ops-{self._flag}.pdf",
             )
         if self.post_process_config["write_operator"]:
-            dglap.write_YAML_to_file(ret, self._operator_path)
+            dglap.dump_yaml_to_file(ret, self._operator_path)
 
     def run_FFNS(self, Q2init, Q2final):
         """
@@ -236,10 +235,8 @@ def get_member(x):
 
     # combine grid
     flag = f"l{n_low}m{n_mid}r{polynom_rank}"
-    xgrid_low = interpolation.get_xgrid_linear_at_log(
-        n_low, 1e-8, 1.0 if n_mid == 0 else 0.1
-    )
-    xgrid_mid = interpolation.get_xgrid_linear_at_id(n_mid, 0.1, 1.0)
+    xgrid_low = np.geomspace(1e-8, 1.0 if n_mid == 0 else 0.1,n_low)
+    xgrid_mid = np.linspace(0.1, 1.0, n_mid)
     xgrid = np.unique(np.concatenate((xgrid_low, xgrid_mid)))
 
     # activate logging

benchmarks/tools.py

@@ -133,7 +133,7 @@ def plot_operator(ret, var_name, log_operator=True, abs_operator=False):
 
     ax = plt.subplot(1, 3, 3)
     plt.title("|error/value|")
-    err_to_val = np.abs(op_err / op)
+    err_to_val = np.abs(np.array(op_err) / np.array(op))
     im = plt.imshow(err_to_val, norm=LogNorm(), aspect="auto")
     plt.colorbar(im, ax=ax, fraction=0.034, pad=0.04)
     return fig

doc/source/Code/IO.txt

@@ -13,10 +13,9 @@ The input run card may contain the following key-value pairs:
 - ``PTO`` - *integer* : order of perturbation theory: ``0`` = LO, ...
 - ``alphas`` - *float* : reference value of the strong coupling :math:`\alpha_s(\mu_0^2)` (Note that we have to use
   :math:`\alpha_s(\mu_0^2)` here, instead of :math:`a_s(\mu_0^2)` for legacy reasons)
-- ``xgrid`` - *list(float)* : the interpolation grid
-- ``xgrid_type`` - *str* : generating function for the interpolation grid
-- ``xgrid_polynom_rank`` - *integer[4]* : polynomial degree of the interpolating function
-- ``log_interpol`` - *boolean [True]* : use logarithmic interpolation?
+- ``interpolation_xgrid`` - *list(float)* : the interpolation grid
+- ``interpolation_polynomial_degree`` - *integer[4]* : polynomial degree of the interpolating function
+- ``interpolation_is_log`` - *boolean [True]* : use logarithmic interpolation?
 - ``q2_grid`` - *list(float)* : the target grid in :math:`Q^2`
 
 

doc/source/Theory/pQCD.txt

@@ -12,6 +12,12 @@ We use perturbative QCD with the running coupling
       \frac{da_s}{d\ln\mu^2} = \beta(a_s) \
       = - \sum\limits_{j=0} \beta_n a_s^{j+2}
 
+Our main reference is :cite:`Herzog:2017ohr`.
+
+The matching conditions are given in :cite:`Chetyrkin:1997sg`.
+
+Implementation: :class:`~eko.strong_coupling.StrongCoupling`
+
 Splitting Functions
 -------------------
 

doc/source/index.txt

@@ -1,3 +1,8 @@
+.. image:: img/Logo.svg
+  :width: 300
+  :align: center
+  :alt: EKO logo
+
 EKO documentation
 =================
 

doc/source/refs.bib

@@ -122,4 +122,16 @@ @article{Vogt:2004ns
     title = "{Efficient evolution of unpolarized and polarized parton distributions with QCD-PEGASUS}",
     volume = "170",
     year = "2005"
+}
+@article{Chetyrkin:1997sg,
+    author = "Chetyrkin, K.G. and Kniehl, Bernd A. and Steinhauser, M.",
+    title = "{Strong coupling constant with flavor thresholds at four loops in the MS scheme}",
+    eprint = "hep-ph/9706430",
+    archivePrefix = "arXiv",
+    reportNumber = "MPI-PHT-97-025",
+    doi = "10.1103/PhysRevLett.79.2184",
+    journal = "Phys. Rev. Lett.",
+    volume = "79",
+    pages = "2184--2187",
+    year = "1997"
 }

refactor-bowler.py

@@ -0,0 +1,8 @@
+# -*- coding: utf-8 -*-
+from bowler import Query
+
+q = (
+    Query("src/eko/")
+    .select_class("Threshold")
+    .rename("ThresholdsConfig")
+)

setup.cfg

@@ -0,0 +1,9 @@
+[tool:pytest]
+norecursedirs = benchmarks/aux
+testpaths = tests benchmarks
+python_files = test_*.py benchmark_*.py regression.py
+python_classes = Test* Benchmark*
+python_functions = test_* benchmark_*
+
+addopts = --cov=eko --cov-report=html --cov-report=xml --strict-markers
+

setup.py

@@ -2,7 +2,7 @@
 from setuptools import setup, find_packages
 
 setup(name='eko',
-      version='0.0.1',
+      version='0.2.0',
       description='Evolution Kernel Operator',
       author = 'S.Carrazza, J.Cruz-Martinez, F. Hekhorn',
       author_email='stefano.carrazza@cern.ch',
@@ -36,5 +36,5 @@
           "cffi>1.0.0"
       ],
       cffi_modules=["src/cfunctions/digamma.py:ffibuilder"],
-      python_requires='>=3.6'
+      python_requires='>=3.7'
 )

src/cfunctions/digamma.py

@@ -1,5 +1,5 @@
-from cffi import FFI
-ffibuilder = FFI()
+import cffi
+ffibuilder = cffi.FFI()
 
 ffibuilder.set_source("_gsl_digamma",
    r""" // passed to the real C compiler,

src/eko/__init__.py

@@ -1,7 +1,8 @@
 # -*- coding: utf-8 -*-
 
 from eko._parameters import t_float, t_complex
-from eko.runner import Runner
+import eko.runner
+
 
 def run_dglap(setup):
     r"""
@@ -24,6 +25,6 @@ def run_dglap(setup):
             output : dict
                 output dictionary - see :doc:`/Code/IO`
     """
-    r = Runner(setup)
+    r = eko.runner.Runner(setup)
     output = r.get_output()
     return output

src/eko/constants.py

@@ -11,10 +11,11 @@
      - :math:`T_F` the normalization of fundamental generators - defaults to :math:`1/2`
      - :math:`C_A` second Casimir constant in the adjoint representation - defaults to :math:`N_C = 3`
      - :math:`C_F` second Casimir constant in the fundamental representation - defaults to :math:`\frac{N_C^2-1}{2N_C} = 4/3`
-"""#pylint: disable=line-too-long
+"""  # pylint: disable=line-too-long
 
 from eko import t_float
 
+
 class Constants:
     """
         Upon initialization all constants are set and a `lock` variable is put is place.

src/eko/ekomath.py

@@ -17,6 +17,7 @@
 cffi_utils.register_module(_gsl_digamma)
 c_digamma = _gsl_digamma.lib.digamma  # pylint: disable=no-member
 
+
 @nb.njit
 def gsl_digamma(N: t_complex):
     """
@@ -43,13 +44,11 @@ def gsl_digamma(N: t_complex):
 
         Note that, although not listed in the documentation, there **is** a
         `complex version <http://git.savannah.gnu.org/cgit/gsl.git/tree/specfunc/gsl_sf_psi.h#n76>`_
-    """#pylint: disable=line-too-long
+    """  # pylint: disable=line-too-long
     r = np.real(N)
     i = np.imag(N)
     out = np.empty(2)
-    c_digamma(
-        r, i, _gsl_digamma.ffi.from_buffer(out)
-    )
+    c_digamma(r, i, _gsl_digamma.ffi.from_buffer(out))
     result = np.complex(out[0], out[1])
     return result