Speedup Grid::convolute_eko

[cschwan]

The runtime for Grid::convolute_eko for double-hadronic grids is, unfortunately, quite long, of the order of 10-15 minutes for a typical case. I'm convinced that we can improve the runtime, and I have at least two ideas how:

1. Eliminate zeros of the EKO early. Commit 6aa3ba4 checks, for all source- and target-PID tuples, whether the operator is non-zero by going through all its elements once. For instance, in many of the grids we have a photon as an initial state, but often we ignore the PDF at the fitting scale and so we can skip everything that has a photon. This is also true for top and anti-top quarks, and possibly for other quark flavour as well.
2. The actual convolution can be rewritten as a loop over the different mu2 values, each performing the following operation in the general case where there are two hadronic initial states: (FK)_ij = sum_k O_il (G_lm)_k O_mj, so that we could leverage fast linear algebra libraries.

However, we need to test and benchmark this.

[codecov]

Codecov Report

Merging #103 (b7014fd) into master (e91a7df) will increase coverage by 0.06%.
The diff coverage is 88.00%.

@@            Coverage Diff             @@
##           master     #103      +/-   ##
==========================================
+ Coverage   88.43%   88.49%   +0.06%     
==========================================
  Files          31       31              
  Lines        3060     3085      +25     
==========================================
+ Hits         2706     2730      +24     
- Misses        354      355       +1     

Impacted Files	Coverage Δ
pineappl/src/grid.rs	88.37% <88.00%> (+0.12%)	⬆️
pineappl_cli/src/helpers.rs	89.89% <0.00%> (-1.12%)	⬇️
pineappl_cli/src/plot.rs	93.15% <0.00%> (+1.05%)	⬆️

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update e91a7df...b7014fd. Read the comment docs.

[felixhekhorn]

@cschwan I added an additional condition to the skipping to make it DIS compatible

[cschwan]

@felixhekhorn cheers, good catch!

[felixhekhorn]

On ATLAS_TTB_8TEV_TOT we went from 3:15 to 2:51, so not that much ... fine for me to merge

[cschwan]

@alecandido @felixhekhorn that was just 1), but the my hope is that 2) gives a much larger improvement (not yet implemented).

[cschwan]

I observed that the EKO for the default grid in Pineko is triangular in the source and target x, but the row with source x has 3 more entries than expected for a truly triangular matrix. That's probably because of the interpolation order (3?). However, this is different for points close to x = 1 (5 more entries). Did you increase the interpolation order in this region?

[felixhekhorn]

As explained in the docs we're taking the N_deg+1 (default: 4+1=5) nearest points - which on the border leads to additional offdiagonals

[felixhekhorn]

2. The actual convolution can be rewritten as a loop over the different mu2 values, each performing the following operation in the general case where there are two hadronic initial states: `(FK)_ij = sum_k O_il (G_lm)_k O_mj`, so that we could leverage fast linear algebra libraries.

However, we need to test and benchmark this.

Can we postpone 2. to a separate PR and close this one?

[alecandido]

Can we postpone 2. to a separate PR and close this one?

Maybe convert 2. in an issue first, since no one has time to dedicate to it right now.

[cschwan]

Or we simply leave it open because now this PR basically is 2 😄. Why do you want to close it?

[alecandido]

Or we simply leave it open because now this PR basically is 2 smile. Why do you want to close it?

Fine by me.
I believe @felixhekhorn wanted to start using 1., but if I remember correctly it wasn't a big change in performance.

[felixhekhorn]

this is getting stale (no activity the last half year) and the longer you wait the more difficult it will be to merge