nnpdf / yadism

Yet Another DIS Module

NNPDF: An open-source machine learning framework for global analyses of parton distributions

The NNPDF collaboration determines the structure of the proton using Machine Learning methods. This is the main repository of the fitting and analysis frameworks. In particular it contains all the necessary tools to reproduce the NNPDF4.0 PDF determinations.

Documentation

The documentation is available at https://docs.nnpdf.science/

Install

See the NNPDF installation guide for the conda package, and how to build from source.

Please note that the conda based workflow described in the documentation is the only supported one. While it may be possible to set up the code in different ways, we won't be able to provide any assistance.

We follow a rolling development model where the tip of the master branch is expected to be stable, tested and correct. For more information see our releases and compatibility policy.

Cite

This code is described in the following paper:

@article{NNPDF:2021uiq,
    author = "Ball, Richard D. and others",
    collaboration = "NNPDF",
    title = "{An open-source machine learning framework for global analyses of parton distributions}",
    eprint = "2109.02671",
    archivePrefix = "arXiv",
    primaryClass = "hep-ph",
    reportNumber = "Edinburgh 2021/13, Nikhef-2021-020, TIF-UNIMI-2021-12",
    doi = "10.1140/epjc/s10052-021-09747-9",
    journal = "Eur. Phys. J. C",
    volume = "81",
    number = "10",
    pages = "958",
    year = "2021"
}

If you use the code to produce new results in a scientific publication, please follow the Citation Policy, particularly in regards to the papers relevant for QCD NNLO and EW NLO calculations incorporated in the NNPDF dataset.

Contribute

We welcome bug reports or feature requests sent to the issue tracker. You may use the issue tracker for help and questions as well.

If you would like contribute to the code, please follow the Contribution Guidelines.

yadism's Issues

Join quad call

Numerics question: is it a good idea to join R+S+L under one integral?

because

• only then the "absoulte error" has its proper meaning
• it might be more acurate or even faster

in order to explore this, one would need to compute an exact convolution with say x^a(1-x)^b for various a and b and check the error and performance

a first but not followed attempt was made in #28

Import external formula for coefficient functions

In order to import external formulas, like ones in vogt, that are usually encoded in small fortran functions, it would be a good a idea to download and collect them in a single 'external' folder, and write a Makefile to compile this folder in a corresponding one importable from python.

Than Makefile should provide an install command to install the new python importable coefficient functions' expressions in yadism package.

In the package everything should be organized as it is, but the ESF implementations (F2light, F2charm, ...) instead of defining themselves the coefficient functions should only import them.

Refactor convolution

Since to convolute with a distribution a distribution is needed it makes sense to think the convolution function as a feature of a distribution itself (even more since the second operand can't be another distribution but should be a regular test function).

So:

• make convnd a method of DistributionVector
• define a custom __add__ method, able to sum dvec with dvec and dvec with regular functions
  • in the second case we assume that the distribution is simply the regular function itself, so has 0 component in the other directions
  • define also __radd__ and __iadd__
• define a custom __mul__ method, to multiply by a scalar
  • define also __rmul__ and __imul__
• define a suitable __iter__ method to be used in convnd function
  • get rid of __getitem__ once done
• in convnd: take care of integration limits properly

Github-workflow

There are still some stuffs that I would like to improve at some point:

• make setuptools responsible for eko installation
  • I prefer not to do it explicitly
  • maybe this would require to make setuptools in eko responsible for GSL installation (at least in principle)
• make the login 'bash -l {0}' and conda environment activation automatic
  • instead of cut and paste the code at each step that sounds always bad
• make actions/checkout@v2 responsible to checkout eko at the required commit
  • maybe this is not yet possible, and maybe it would be worth to open a feature-request on actions/checkout@v2
• remove my personal PAT (ALE_DIS_PAT) from the workflow
• extend the matrix
  • in principle it should work at least for python 3.8 too

Logging

Add log info (and verbose mode) with python logging library.

convolution

instead of making quad cluttering the output with integration errors logging them is the proper thing to do, so just use the full_output argument of quad, collect them and put in the global logs.

Numba optimization

Functions integrated by scipy.integrate.quad are currently composed by nesting multiple lambdas.

Calling this objects can be a bottleneck, so:

• we need to profile a full run of yadism while producing a "DIS operator"
• if needed: optimize the object passed to quad by constructing a proper numba object

More benchmarks

• Include benchmarks suggested by Juan Rojo, see the wiki.
• Include even more benchmarks, for example the code by Paolo Nason/Eric Laenen we discussed with Stefano C. (you can find it in https://github.com/nnpdf/external/tree/master/FONLLdis)
• xspace-bench

Code of conduct & co.

Just put:

• CONDUCT file
• issue template
• pull request template

Cut ESF out of PartonicChannel

The keyword is decoupling.

Since an ESF is a huge object collecting PartonicChannels and more (weights, is aware of computing methods) is a bad idea to pass to PartonicChannels a reference to all this mess.
Moreover ESF itself holds a reference to SF, and so literally to everything.

Q: What do we need for PartonicChannels from ESF?
A: cd cc; grep SF && cd nc; grep SF

Seriously there are only 4 objects needed:

1. nf, needed only for light
2. Q2, M2hq needed for both asy and heavy
3. x needed only for heavy (to set the threshold)

Just find a way to pass down this stuff by value.

Testing with limited resources

Separate apfel cache

Since apfel is static (from the point view of its implementation) as long as we keep the same input there is no reason to rerun separately in every workflow, so we can store the cache somewhere accessible from the workflow.
The following two strategies (and a mixed one) are available for generation:

1. keep apfel cache across runs (so push back the one generated by workflows)
2. generate apfel cache somewhere else (e.g.: Milan cluster)

Statistical approach

In order to keep the test size constant (or not growing exponentially with features), we can always run the test for isolated features, but selecting randomly only some combinations for each run, since there is no principle to decide a priori that some combination are more important than others.

• decide a fixed amount of slots
• decide some tests to be always run (e.g.: isolated features)
• attribute randomly the free slots
• keep logs of each random selection (for reproducibility and more)
• leave some hooks for adding more deterministic tests (e.g.: for running again failing ones)
• avoid most recent combinations in the random selection (for a quicker exploration)

NLO documentation

Update documentation

Document anything written up to this point.

yadism

• each module
• each function/method/class
• add something directly in the sphinx source (e.g. physics, formulas, ...)
  • or maybe not and write stuffs in the module and class headers

tests

Document also tests:

• what they test, why they are defined
• which kind of bug was found
• some structure in the module headers

Computational graph toolchain

Keeping investigating from #29 we decided to have a look into a proper way to replace and speed up the hard part of this library.

The goal

The situation is the following:

• we have a perturbative expansion f(x) = sum(alpha**k * f_k(x)), that involves multiplying and summing functions
• thus we need lazy evaluation
• at the end we will integrate the result

Proposals?

We are searching for a proper tool. Currently:

• lazy evaluation is implemented nesting lambdas, and the only efficiency measure is try to keep nesting as minimal as possible, based on a manual type check (the best we can achieve with lambdas only)

what we would like:

• a library for dealing automatically with computational graph; maybe we just need tensorflow, but maybe what we need is just the compile part, that maybe is xla (I'm still trying to understand the tf internals), and so jax would be enough (it just ship autograd, that we don't need, together with xla and nothing more, rather minimal if compared with tf)
• maybe also the integration should be managed on the same level, do we need something like tfp.math.ode.Solver to target performance?

Improve `against_apfel` tests

Cache

Compare with a cache of APFEL's results, if available

Split observables

Make a separate test for each observable: F2light, FLlight, F2charm, ...

Order by order

When introducing a new perturbative order the contribution of the new order should be compared with the correspondent from APFEL, instead of the full object.

Maybe the two things should be combined (if the correction is too irrelevant maybe it should not have to be too accurate)

Relative - absolute error

Probably no one of the two things make sense on its own, so maybe if the relative error is big enough also the absolute one should be checked against something else (what?), if it is too small maybe it's not really an error.

Candidates

• average over xs, Q2s
• some relevant standard point
• proposals?

Virtual puzzle

Some virtual contributions to DIS seems missing from calculations, there will be for sure a good reason but we haven't found it:

• photon - gluon @ NLO
• quark loop @ NNLO

In particular the first diagram will make a virtual correction proportional to nf at NLO, and both to nf^2 at NNLO, but it is easy to see that there is nothing like that around...

Inherit hadronic target to PineAPPL

in the PineAPPL generation we're hard coding at the moment the hadronic target: this might become a problem with NuTeV and CHORUS

https://github.com/N3PDF/yadism/blob/23c413de082d550b2161512b90a2157e2e69c3c2/src/yadism/output.py#L240

(the lepton is already correct:
https://github.com/N3PDF/yadism/blob/23c413de082d550b2161512b90a2157e2e69c3c2/src/yadism/output.py#L241
)

Reorganize Runner

• make it an object
• separate layers
  • input layer (inspector and so on)
  • output layer (an object that collects output and is able to format and write them)

Intrinsic charm

Literature:

• PhysLett.: https://inspirehep.net/literature/1395621
• EPJC: https://inspirehep.net/literature/1464831
• JHEP: https://inspirehep.net/literature/1397009

Issues:

• remember to update the weights for IC=0, current scenario is:
  • consistent IC=0 for NC at the level of the coefficient functions (LO not provided for heavy structure functions), not sure at the level of the weights (maybe they are provided also for heavy inputs)
  • consistent IC=1 (but also IB=1 and IT=1) for CC: all the weights are always implemented and coefficient functions as well, for both light and heavy inputs
• the correct setup for IntrinsicQuark (=IQ) specified in [m,n) range should be:
  • the weights should be implemented up to n (excluded), because all the heavy not intrinsic inputs should not be present as PDFs
  • the FNS scheme should be consistent, so only the flavors up to m (excluded) can be considered as light in the fitting region (than if it is a VFNS it can change at higher Q2)

Note that the heavy object here is referred to the input, not to the output/object coupling to the vector, so the heavy structure functions should be always implemented for both heavy and light inputs, and are the weights that controls if they are needed or not, not the coefficient functions (so the fact that is actually working for the NC is just a chance)

Structure Functions unit tests

Define unit tests for all ...StructureFunctions classes and logic:

• StructureFunction
• ESF
• ESFTMC
• ...

Useful ideas:

• factorize base mocking for all needed classes
• use pytest.monkeypatch to fake intermediate object difficult to access (you can register functions/objects to be used in place of something else, even without entering inside)

Implement travis pipeline

Steps:

• create a python setup.py
• create a conda-recipe
• enable travis for linux and osx
• enable tests

Alpha_s evolution

APFEL is doing something strange about alpha_s evolution order and FNS (doing one thing for "FFNS && NC" and another one for "everything else").

We should answer to following two questions:

• do we need to distinguish between two cases or alpha_s evolution should be defined in a unique way? (our temporary answer: unique, because we don't see nor found any place in which the opposite is argued)
• which one of the two is the correct one?
  1. lowering and raising perturbative order
  2. keeping the same perturbative order for DIS calculation and alpha_s evolution
     (our temporary answer: 1)

NLO precision

Both the light and heavy structure functions at NLO are still off by a few percent (even more in some cases).

This should be fixed, stripping the difference due to interpolation.

Split heavy quark masses

At the moment the theory runcard is defining three kinds of heavy quark masses:

• mc: the mass of the charm (e.g.)
• Qmc: the reference scale at which mc is evaluated as a Lagrangian parameter in a renormalization scheme (say MSbar)
• kcthr: so the actual new mass is given by kcthr * mc, and it is the threshold used by any VFNS for charm

For yadism they all coincide currently, because:

• mass schemes other than POLE are not implemented (so Qmc is irrelevant)
• kcthr is always considered one

Let's fix at least the second, and we will deal with the first as soon as HQMS will be implemented.

Implement entry point test

Steps:

• copy a LO F2 setup from NNPDF
• create a python test as dict
• send this dict to a run_dis function

Report parameters

TL;DR: default and/or fallbacks with report parameters or raising errors for anything, including missing arguments?

I'm wondering how to do this (actually if to do this).

In apfel once that defaults are applied, conflicts are solved and so on, it prints the parameters that is actually sing for both the evolution (not in this project) and the DIS observables (if you initialize also apfelDIS).

Discussing with @felixhekhorn this week he suggested that this feature may be not reimplemented, instead it's better to raise an error if there is any conflict at all.
I quite agree with Felix, errors are more pythonic than hidden fallback for sure, and I think that it is also better in general, but now I'm considering if including any default at all.

So the question is: if there is any input, any default of any kind, I think that would be better to report parameters used in some way (dumping on a file, printing as output, ...), otherwise we should raise an error also for any missing argument.

Heavy Quark Mass Scheme

The Heavy Quark Mass Scheme implementation it's not so hard, it will require:

• to get the pdf derivatives from eko
• to implement some new coefficient functions for the CC
• to figure out the interplay with IntrinsicCharm

Indeed at NLO there is nothing to do on NC (apart for the IntrinsicCharm case) because it is just LO Heavy Quark.

Reference: arXiv:1011.5790

Plus distributions - general form

We are wondering if it is possible on general grounds to reorganize something of this kind:

using constant coefficients for the distribution part:

If this is possible:

• we are thinking about rewriting our DistributionVec a little in order to collect all the singular bits in a single one (so we will have a fixed amount of components equal to 3 for all orders)
• just do this we only need to:
  • include in coefficient functions definition the log(1-z) ** k / (1-z)
  • multiply in the convolution just for the pdf difference pdf(x/z)/z - pdf(x), that is available since what we currently are doing is multiplying each one of these two terms by a factor that is constant in the second expression (\tilde{c}_k)

QED corrections

Add 'lhapdf' and 'eko' dependence in Travis

How can i add the dependence on

• lhapdf
• eko

in meta.yaml?

The second one prevent me to import interpolation.py (that I'm currently copying) and the first one is definitely causing the failure of the build.

features

see https://github.com/N3PDF/dis/wiki/List-of-features-to-be-implemented

Precompile nested functions

Since we are using a bunch of lambdas and functions calls in other functions (cf. all the full DistributionVec API) maybe at some point this can slow down the execution.

Instead of find a way to avoid this nesting and so on, rewriting stuff in a proper fashion with even more tricks, it would be nice to find a way to precompile the functions before using them (e.g. for integration in scipy.integrate.quad where a lot of functions' calls are involved).
Of course, if needed, this can involve rewriting the functions definition in a proper fashion, adding suitable decorator or whatsoever, but the goal is to avoid to change critically the structure because of this.

My idea of precompiling: of course I'm not an expert, otherwise I would have already figured out the proper way, but in the worst scenario means define a function that takes another one as input and go through the calls collecting everything in a single expanded expression, to be evaluated again as python code (hopefully an external library can do this for us, I don't know if numba or what else)

Modules

This code must contain at least 2 modules:

• the DIS observable implementation
  • by direct evaluation
  • through DIS operator (i.e. tensor that will be integrated with EKO)
• a simple API to retrieve both previous values.
  • we can think a similar approach to EKO, i.e. dict with setup and evaluation function which returns the DIS operator tensor.

In the first module we should restrict ourselves to C-style programming, this is due to the limitations imposed by numba.

FONLL-C Fcharm not working

Investigation work in progress.

QED corrections

• do we want to implement those?
• if yes: speak with Christopher

From the email exchange on 19.08.20 between SF, ERN, SC, CS, FH, AC we can already give two references

• "Factorized approach to radiative corrections for inelastic lepton-hadron collisions" https://arxiv.org/pdf/2008.02895.pdf
• https://wiki.bnl.gov/eic/index.php/DJANGOH

ESF output object

Make another object for managing ESF ouput structure, that for the moment consists in:

• x
• Q2
• q
• q_error
• g
• g_error
  but it is subject to change.

Furthermore it's quite uncomfortable not to be able to sum or multiply by a scalar output objects, since they are dictionaries, so it would be good to implement proper methods in the new class and use them for TMC calculation.

Implement automatic documentation

Steps:

• create a doc folder with sphinx setup
• add minimal layout for that
• integrate automatically on #3 .

Released based workflow

Set up released based workflow

Implement basic F2 calculation

Steps:

• encode the LO F2 observable
• benchmark results against APFEL

CC light NNLO: even and odd moments

Properly deal with them, even though at the moment it looks like that doing it improperly is not an issue for some observables...

Visualizations

some ideas of visualizations we may add to yadmark (as they do not belong into yadism itself)

1. for the true structure functions predicted by a certain PDF, we should show the traditional DIS plot (e.g. Figure 18.2 in PDG): y-axis: F, x-axis: x, Q2 by offset

2. for the DIS operators we could do smth similar to EKO:

• y-axis: position of basis function in grid
• x-axis: x
• make pixel with color = strength
• for fixed Q2 and PID

in LO this should generate smth like:

| o o x
| o x o
| x o o
-------

because this represents the delta-function

• However in NLO we would need to plot the pure NLO corrections maybe as everthing should be dominated by LO and I wonder if this would be visible in a combined plot
• in principle the operator should be physical and thus strictly positive, but in practice basis function can and will be negative, so the operator may go negative, so we would need to show positive and negative contributions

Output format revamped

Since the success in eko new output format implementation NNPDF/eko#76 I propose to follow the same path here, and so to split the operator from the metadata.
Main reason supporting this is that we'll save a lot in loading and dumping time.

Output structure

Here the typical structure of output is a bit different from those in eko, and a bit more complicate (but still perfectly suitable).

An Output object is still basically a dictionary, e.g. with the following keys:

[ins] In [19]: out.keys()
Out[19]: dict_keys(['XSHERACC', 'interpolation_is_log', 'interpolation_polynomial_degree', 'interpolation_xgrid', 'pids', 'projectilePID'])

At the moment we need an input runcard in order to be able to load the observable, because we don't know otherwise which are the observables (e.g. XSHERACC) names in the output object.

There are two possible solutions:

• we scope observables within a single entry obsvervables in the object instance
• we keep as it is, but we add an observables entry with a list of names of observables

While every other value is just a scalar (at most a list/array of scalars), each observable is a very nested object.

Observables -> ESF structure

Each observable is a list of ESF, potentially of different length. In turn, each ESF has the same structure.

[ins] In [14]: e = out["XSHERACC"][0]

[ins] In [15]: len(out["XSHERACC"])
Out[15]: 42

[ins] In [16]: np.array(e.get_raw()["orders"][0]["values"]).shape
Out[16]: (14, 50)

[ins] In [17]: e.get_raw()["orders"][0].keys()
Out[17]: dict_keys(['order', 'values', 'errors'])

[ins] In [18]: e.get_raw().keys()
Out[18]: dict_keys(['x', 'Q2', 'nf', 'orders'])

Thus we can make it an array, what we need is:

• one array per observable
• first dimension runs over ESFs (so we need to store x, Q2, and nf somewhere else, i..e in metadata)
• second dimension runs over orders; they are the same for each observable -> we need to store the list of orders
• two other dimensions run over flavor and interpolation grid (both already stored, respectively as pids and interpolation_xgrid)
• optionally one more dimension can run over [value, error]

Recap

We need:

• a tar archive containing
• metadata, i.e. all key-value pairs, but observables, to which we should add
  • the list of orders
  • x, Q2, and nf for each ESF in each observable (so one list per observable, made of 3-tuples)
• dump one array per observable; it's to dump mandatory separately because they might not be uniform

Benchmarks: Automatic detection of flavors cancellation

What

We noticed that there is a common source of differences (rel_err[%]) between us and APFEL coming from Q2 interpolation (or also anything APFEL-like, i.e. Q2 interpolating)

Proposed solution

• in yadmark we can get the numbers for the separate channels from yadism (the "DIS operator") simply calling runner.get_result() and only after calling .apply_pdf() on the result
• once the numbers for the separate channels are available we are able to detect the region, or better the exact ESF where a cancellation between flavors is happening (just looking for the maximum of each flavor channel and computing the absolute ratio with the summed result)
• once that a cancellation is detected raise the absolute error for APFEL comparison

Absolute error value problem

The actual number to choose for the absolute error is difficult to find, because the proper value should be the sum of absolute errors on the individual flavor channels, but in order to find this we should be able to break the APFEL's result as well on the individual flavors, but if these results where available we could have compared directly the channels instead of the recombined sum.

My personal guess for this number is a function of the sum/max(flavors) ratio, protecting for yadism recombined result going to 0 (and so ratio=0 -> rel_err=100%) and protect for the full yadism result where the cancellation is still present but it is more likely APFEL to go to 0 (1%<ratio<20% -> abs_err=yad_result).

Improve logs

Info level

Info printed during yadism runs are still partial, to be included

• FNS
• TMC (if enabled)

To be reorganized:

• couplings related: e.g. QCDNUM needs to load the coupling constants, and so it's doing many calls to get_weights, but it's reasonable that some values are constant per-run, like MZ and so on, and instead they are printed over and over
  • maybe we should make an object that loads some initial configs, and log only when instantiating the object, or changing top-level values (like MZ), while not logging for each and every

Debug level

• logging more in distribution_vec?

Differences to APFEL

At LO order APFEL is already dependent on Q2, while we are not: why?

Also:

• we have still ~0.05% error at LO, it's small but still not ~0
• Q2 dependence it's becoming relevant because it is probably enhanced by xiF, xiR != 1.0

Speed up

At some point we might deserve another round of analysis on:

• what is still not numbified? why? can we fix it?
• is it still to slow in some actual use case? can we quickly improve with a cheap parallelization?

Travelling

New observables

• heavy flavors
  • include properly F2-FL, and also the new ones
• decide something about / x * z

Multiple docs versions

Goals

• multiple versions of the docs deployed
• redirect to master/stable from the base domain (or write a home page)
  • when ready: delete all top level files, keep only .nojekyll, the redirecting, and the automatically deployed subfolders
• update the docs to show the correct version (below the logo)
  • instead of using the one manually update the file in the workflow with the content of $dest_dir (of course before building)
• add the switching version widget

Rationale

It seems useful to be able to look at the developing docs online, furthermore if there are multiple authors working on it, or if you want to present something to someone else.
On the other hand it seems a really bad idea to replace the master documentation with some unstable version, so it is important that the main (landing) one it is the one related to the current release (or master for the time being).

Number of active flavors

Which is the number of active flavors?

1. the Threshold object is the main responsible for the number of flavors always

• both in yadism
• and in eko

2. MaxNfPDFs should affect the number of flavors relative to the input

• both in yadism
• and in eko

3. MaxNfAlphas should affect the number of flavors in evolution

• `eko' only

4. the number of active flavors in the output is relevant only in the light structure functions, and in that case nf = 3, because light means u + d + s
5. which is the number of flavor to be used in the internal quark loops? the correct answer is "ask Vogt", or better look calculation details in proper papers

Deliver in the Flavor basis

In order to make things easy for the user the current setup is not the best one.

There are two basis involved:

1. the coefficient function basis: singlet, non-singlet, gluon
2. the pdf basis: flavor basis
   and the weights are a bridge between them.

Since the first one is needed only for internal purpose we should:

• apply weights internally (not showing them to the user)
• rotate the output in the flavor basis
  • maybe instead we should provide all the ingredients separately:
    • the coefficient functions vector (in their own basis: g, s, nsp, nsm)
    • the basis rotation (13x4 matrix)
    • the weights matrix (diagonal matrix of charges)

Further on benchmarks management

Split ParentTest class:

• conftest.py will contain everything needed at runtime, imported as a fixture
• something-else.py will contain all the management - analysis tools

Add _modified/_creation_time metadata to inputs:

• check also for cache timestamp to be consistent, otherwise recompute