Expose grid.evolve to the C-API
#344
Conversation
|
Hi @Radonirinaunimi, yes, it does make sense. |
|
Hi @vbertone, apologies for the late reply, I wanted to sort out the interface first before providing you with how the order of the flattening should look like. However, I had not anticipated how difficult the interface would be, and therefore this is not done yet. It might still takes some time. |
|
Hi @Radonirinaunimi, no problem at all. Take your time. |
Ok, this is not completely done (needs some clean up, fix some memory leak in the example, etc.) but I think we can start to test and experiment it. So there is a primary example in evolve-grid.cpp. The operators should be flattened following the std::size_t index = 0;
for a { for i { for b {for l { eko[index++] = ... } } } }The pineappl/examples/cpp/evolve-grid.cpp Lines 62 to 71 in 43b1619 Please let me know if something is wrong/doesn't make sense. |
|
This is now doing what we intended to do. @vbertone please test it and let me know. You can download the Grid for the test via: wget --no-verbose --no-clobber https://data.nnpdf.science/pineappl/test-data/LHCB_WP_7TEV_opt.pineappl.lz4@cschwan This is not super-clean yet and would appreciate if you can have a quick look (your feedback is always quite helpful). |
|
Hi @Radonirinaunimi, thanks a lot, I will test it as soon as possible! |
|
Hi @Radonirinaunimi, I just ran the code and it does work fine. Also, starting from here I should be able to follow the steps to eventually interface it to the evolution operator computed with APFEL++. I just have one question. Here is the output I get: So differences between grid and FK-table predictions are on the verge of 1%. Since the grid that is being used is (I guess) a realistic one and that the evolution is essentially the unity, this seems a purely numerical effect related, I suppose, to interpolation. The question is: is this accuracy satisfactory? And can in principle be improved? Thank you. |
|
Hey @vbertone and @felixhekhorn, there was a bug in the masking of the orders of the example (fixed in b2bf81b). Now the agreement is: |
|
Thank you, @Radonirinaunimi, this looks much better! |
|
Hi @Radonirinaunimi, I'm having some difficulties in understanding the behaviour of the code. Specifically, I would like to change the output PIDs because evolution can possibly generate channels that are not originally present. To this purpose, I made the following basic changes to the code The purpose is to allow for all possible channels. However, I must be doing something wrong because when I run the code I get this error that I cannot interpret: I also tried to play around with the order of variables here and there, but still get the same error (except that sometimes the two tensor shapes in the error message are swapped). Any idea of what I'm doing wrong? Thank you. |
|
Hi @vbertone, the error is because a mismatched shape between the EKO |
Hi @vbertone! My bad, my thinking, at least for that first proof of concept, was to expose as little objects as possible by abstracting away all the I believe that in 9efc9bc the exposed function now contains the least amount of needed arguments to construct pineappl/pineappl_capi/src/lib.rs Lines 2101 to 2115 in 9efc9bc Please let me know what you think. PS: @cschwan, in case you'd like to get access to the |
|
Hi @cschwan and @Radonirinaunimi, ok, understood, thanks. I will try to use this new interface to construct the evolution operator with |
|
@Radonirinaunimi, concerning the question of alpha_s, I don't think that's an issue. The function that returns alpha_s is fast enough that the overhead of computing it twice is most probably negligible. |
So, assuming that extern "C" void generate_fake_ekos(
const int* _pids_in,
const double* _x_in,
const int* pids_out,
const double* x_out,
double* eko_buffer,
void* pdf_state,
pineappl_conv_type conv_type,
double fac1,
std::size_t pids_in_len,
std::size_t x_in_len,
std::size_t pids_out_len,
std::size_t x_out_len
) {
// Ignore unused variables
(void) _pids_in;
(void) _x_in;
(void) pids_in_len;
(void) x_in_len;
// ------------------ Construct the Evolution Operator with APFEL++
// ------------------ Accumulate of thresholds and PIDs from LHAPDF set
std::vector<double> Thresholds; //{0, 0, 0, 1.3, 4.75};
std::vector<int> pids_in; //{-5, -4, -3, -2, -1, 21, 1, 2, 3, 4, 5};
const double mu0 = static_cast<LHAPDF::PDF*> (pdf_state)->qMin();
for (auto const& v : static_cast<LHAPDF::PDF*> (pdf_state)->flavors()) {
pids_in.push_back(v);
if (v > 0 && v < 7) {
Thresholds.push_back(static_cast<LHAPDF::PDF*> (pdf_state)->quarkThreshold(v));
}
}
// Get perturbative order from LHAPDF set
const int PerturbativeOrder = static_cast<LHAPDF::PDF*> (pdf_state)->qcdOrder();
// Construct running strong coupling
apfel::AlphaQCD a{static_cast<LHAPDF::PDF*> (pdf_state)->alphasQ(91.1876), 91.1876, Thresholds, PerturbativeOrder};
const apfel::TabulateObject<double> Alphas{a, 200, 0.9, 13001, 3};
const auto as = [&](double const& mu) -> double {
return Alphas.Evaluate(mu);
};
// Construct x-space grid for APFEL++ (can this be determined dynamically?)
// const apfel::Grid g{{apfel::SubGrid{100, 1e-5, 3}, apfel::SubGrid{100,
// 1e-1, 3}, apfel::SubGrid{100, 6e-1, 3}, apfel::SubGrid{80, 8.5e-1, 5}}};
const apfel::Grid g{
{apfel::SubGrid{80, 1e-5, 3}, apfel::SubGrid{40, 1e-1, 3}}};
// Get joint grid
std::vector<double> x_in(g.GetJointGrid().nx() + 1);
for (size_t i = 0; i < x_in.size(); i++)
x_in[i] = g.GetJointGrid().GetGrid()[i];
// The Evolution Operator is a vector with length `N_conv * N_Q2_slices * Σ
// product(OP shape)` Initial scale
const std::size_t flat_len = pids_out_len * x_out_len * pids_in.size() * x_in.size();
std::vector<std::size_t> shape = {pids_out_len, x_out_len, pids_in.size(), x_in.size()};
// Zero operator
const apfel::Operator ZeroOp{g, apfel::Null{}, apfel::eps5};
// Run over convolution types
int k = -1;
// Construct DGLAP objects for the relevant evolution and
// tabulate evolution operator
std::unique_ptr<apfel::Dglap<apfel::Operator>> EvDists;
if (conv_type == pineappl_conv_type::PINEAPPL_CONV_TYPE_UNPOL_PDF)
EvDists = BuildDglap(InitializeDglapObjectsQCDPhys(g, Thresholds, true), mu0, PerturbativeOrder, as);
else if (conv_type == pineappl_conv_type::PINEAPPL_CONV_TYPE_POL_PDF)
EvDists = BuildDglap(InitializeDglapObjectsQCDpolPhys(g, Thresholds, true), mu0, PerturbativeOrder, as);
else if (conv_type == pineappl_conv_type::PINEAPPL_CONV_TYPE_UNPOL_FF)
EvDists = BuildDglap(InitializeDglapObjectsQCDTPhys(g, Thresholds, true), mu0, PerturbativeOrder, as);
else
throw std::runtime_error("Unknow convolution type");
const apfel::TabulateObject<apfel::Set<apfel::Operator>> TabulatedOps{
*EvDists, 100, 1, 13000, 3};
// NOTE: The Loop over the scale `fac1` is also removed
....
}with the main changes that the loops on |
|
Hi @Radonirinaunimi, yes, this does make sense and makes it much clearer. Thank you. I still have a doubt though. In this example the evolution parameters are retrieved from the LHAPDF set, which in turn is passed to the function that returns the evolution operators through |
Hmhm, I understand you mean by generic interface. Let me think a bit more.
Could you perhaps provide an example of such a general case? In this way, the next iteration could be based upon that, and therefore also general xd. |
|
Hi @Radonirinaunimi, you mean an example of implementation? At the moment I don't have anything specific in mind. What I mean is: does the current signature of the function that returns the evolution operator accommodate a situation in which the For example, one may define a structure which encapsulates the evolution parameters and pass an instance of that structure through |
Yes, as long as it is not a complicated object (such as function, although function pointers might still work). See for example the patch below for evolve-grid-identity.cpp: diff --git a/examples/cpp/evolve-grid-identity.cpp b/examples/cpp/evolve-grid-identity.cpp
index f68a094d..42624e81 100644
--- a/examples/cpp/evolve-grid-identity.cpp
+++ b/examples/cpp/evolve-grid-identity.cpp
@@ -12,6 +12,11 @@
// NOTE: Uses the scale of the Grid as the starting scale such that we can use an IDENTITY EKO.
double FAC0 = 6456.44;
+struct ApfelxxParams {
+ int x_nodes;
+ double as_ref;
+};
+
std::vector<std::size_t> unravel_index(std::size_t flat_index, const std::vector<std::size_t>& shape) {
std::size_t ndim = shape.size();
std::vector<std::size_t> coords(ndim);
@@ -30,7 +35,7 @@ extern "C" void generate_fake_ekos(
const int* pids_out,
const double* x_out,
double* eko_buffer,
- void* pdf_state,
+ void* params,
pineappl_conv_type conv_type,
double fac1,
std::size_t pids_in_len,
@@ -46,9 +51,10 @@ extern "C" void generate_fake_ekos(
(void) conv_type;
(void) fac1;
- // Check to get the μ0 from the PDF
- const double mu0_scale = static_cast<LHAPDF::PDF*> (pdf_state)->q2Min();
- (void) mu0_scale; // Mark as unused variable
+ // Check external APFEL++ parameters
+ ApfelxxParams* op_params = static_cast<ApfelxxParams*>(params);
+ std::cout << op_params->as_ref << "\n";
+ std::cout << op_params->x_nodes << "\n";
std::size_t flat_len = pids_in_len * x_in_len * pids_out_len * x_out_len;
// NOTE: The EKO has to have as shape: (pids_in, x_in, pids_out, x_out)
@@ -175,6 +181,12 @@ int main() {
alphas_table.push_back(alpha);
}
+ // ApfelxxParams
+ ApfelxxParams* op_params = new ApfelxxParams;
+ op_params->x_nodes = 10;
+ op_params->as_ref = 0.118;
+ auto apfelxx_params = static_cast<void*>(op_params);
+
std::vector<double> xi = {1.0, 1.0, 1.0};
// NOTE: The EKO has to have as shape: (pids_in, x_in, pids_out, x_out)
std::vector<std::size_t> tensor_shape = {pids_in.size(), x_in.size(), pids_out.size(), x_in.size()};
@@ -193,7 +205,7 @@ int main() {
// - `ren1`: values of the renormalization scales
// - `alphas_table`: values of alphas for each renormalization scales
pineappl_fk_table* fktable = pineappl_grid_evolve(grid, opinfo_slices.data(),
- generate_fake_ekos, max_orders.data(), pdf.get(), x_in.data(),
+ generate_fake_ekos, max_orders.data(), apfelxx_params, x_in.data(),
x_in.data(), pids_in.data(), pids_out.data(),
tensor_shape.data(), xi.data(), ren1.data(), alphas_table.data());Would this be enough? PS: I should rename this to |
|
Hi @Radonirinaunimi, yes, that's precisely what I had in mind. I'll try that and let you know. Thank you. |
|
Hi @Radonirinaunimi and @cschwan, I managed to modify the code in Things seem to work fine. However, when computing the convolution between grids and evolution operators in the case of two convolution types, neither the memory footprint nor the performance seem to have got significantly better. I uploaded the corresponding code in the |
|
Hi @vbertone, out of curiosity, how are you estimating the memory usage? As for the runtime, I did not expect the speed to improve much to be honest. I am currently profiling the codes with the previous version and the current improvement but it is taking quite some time 😬 (I am using while previously it was of the order of Gbs (I am still checking the stack). |
|
Hi @Radonirinaunimi, I just used But if I use I finer grid (~100 nodes) I find: So now the memory usage is around 10G. Let me know if there is anything you would like to discuss. |
|
Stupid question: How many nodes are used in the current remote version of |
|
The current version uses the coarse grid. The grid is controlled through these lines: The first is perhaps excessively dense, the second is kind of fine, the third (which is the one commented you are probably using) is not accurate enough. |
|
Hi @vbertone, sorry the bit of a delay. Using the dense grid, I now get consistent benchmarks with yours: In hindsight, I'd say this is a fact of life. The optimization above was meant to not store all of the evolution operators in memory all at once but only one at time. While here, the problem is really with a single operator (maybe this is what you meant to tell me since the beginning but I missed it). I am not sure (and I don't think) there is really around this tbh (cc @cschwan) other than using reasonable number of nodes (in our case, we always use 50 for which the accuracy is acceptable). Nevertheless, I'd say that if 100 nodes can still be run on a personal computer, that is not too bad. |
|
Hi @Radonirinaunimi, ok, then that's what it is. Actually, I was wondering whether SIHP could be somehow optimised. If I remember correctly, you were mentioning that all partonic channels contribute to the cross section, which would make 13^3 (or 11^3) channels. If so, I'm pretty sure that not all of them are independent. Wouldn't that allow us to reduce the amount of computations do be done for the combination? Also, I launched the combination with around 100 nodes yesterday evening and it's still running. Hopefully, I should get it done by today. But still, the computation time it's of the order of one day, which is pretty heavy. |
The SIHP grids are already optimized so that non-contributing channels are removed. See here for example. But on that note, I am actually wondering whether the memory and speed issue is related to the explicit double-precision representation of the grid (?). But I am just thinking out loud... |
|
Hi @Radonirinaunimi, the code finally converged (~20 hours) and here is the output: The accuracy is not ideal but I suspect it has to do with the time-like evolution which does not precisely match that of the FF set that I'm using. But this is definitely encouraging. |
|
Hi @vbertone, that took quite some time 😅 As for the accuracy, that is indeed definitely because of mismatch in the evolution, which is also what I recall when I evolved these grids with The question is then: with a more reasonable number of grid nodes (say 40 or 50) would the accuracy still of the same order (without significant deterioration)? If so, then we can just run with coarse grid. |
|
Sure, I just launched a new run with a smaller grid with around 50 nodes. I'll let you know how long that takes and what the accuracy is. Tomorrow, I will also produce an LHAPDF grid with the exact same evolution as APFEL++ for the FFs. This will help us gauge any inaccuracies. |
|
Hi @Radonirinaunimi, the run with the coarser grid (~50 nodes) has come to an end and surprisingly it seems to have given rise to better accuracy than the finer grid: I wonder why. It might be accidental? |
Generally you want the number of channels to be as small as possible before the evolution. If you have the same hadronic initial states symmetrizing the channels is a quick option, but looking at @Radonirinaunimi's comment that seems to be done already. Furthermore, I think I even found an algorithm to minimize the number of channels further to the minimum possible. Furthermore, make sure you do not evolve and perform a basis change at the same time, because the basis change usually blows up the operators a lot. Instead, perform the evolution first and only then change the basis. On the technical side, I think we want to make sure that
|
Which orders are in your grid? If you have just LO + NLO and the NLO correction is very small, it could come from numeric cancellations. In any case, you can try evolving LO and the NLO correction separately and see how that changes the agreement. |
Hi @cschwan, could you please be more specific on this point? What the code does at the moment is to perform the evolution in a specific basis internal to
I'll leave @Radonirinaunimi to comment of these technical points. |
I think I'm actually doing LO + NLO. Specifically, I'm using: Ok, I can try the two orders separately, but I'm not quite sure how to do that. |
What I mean is: does the tensor you feed into I explained the problem a bit better here, NNPDF/pineko#183 (comment). From the numbers I reported there, I believe it's better to 1) evolve first without changing the basis, as the operators are smaller (22 GB vs. 45 Gb) and then 2) only after evolving perform the rotation using |
|
Hi @cschwan, ok, I think I see what you mean. But just to make sure I'm not mistaken, what you suggest is to combine hard matrix elements and evolution operators in the evolution basis, and only after the convolution rotate it into the physical basis. Did I understand correctly? If so, I can definitely provide the evolution operator in the evolution basis, but then I would not know how to tell |
|
@cschwan, we don't perform basis rotation exactly for this reason. So the FK table is computed in the same basis as the grid's (which is always the physical basis/PDG basis). What @vbertone mentioned about
I will have a closer look at this part which I haven't payed much attention when doing the benchmark. Perhaps there's something we can learn from there. |
Sorry, I was only a couple of minutes late with my comment 🙈. Based on #344 (comment), this would not be needed. |
My dad, it looks like I totally misunderstood what is going on under the hood 😄. |
Using
grid.evolveto evolve grid with APFEL++ evolutionIn addition to performing the evolution on the fly and outputting theory predictions, APFEL++ also wants the possibility to dump the evolved grids on disk for future usage (aka FK tables). This effectively implies
grid.evolveto be exposed at the C interface.@vbertone, the idea is then for you to provide an operator which is a rank-4 tensor. To be more illustrative, let's consider a simple DIS case. The convolution, for a fixed$Q^2$ , writes as:
where:
These form the dimensions of an EKO, ie$x$ -grid coordinates of the grid and FK table do not have to be same, and they do not have to be the same as what is actually inside the grid $-$ same is true for PID bases.
EKO[a][i][b][l]. Notice that theThe C++ function that consumes the operator from APFEL++ would therefore tentatively look as follows:! Deprecated, see evolve-grid.cpp.and this is the function that'll dump the evolved grid. I think this is the best way to do it as opposed to what we discussed during the meeting.
So the summarize, what APFEL++ needs to provide (for a given value of$Q^2$ ), is a rank-4 tensor representing the EKO. In addition, we need the $x$ and PID values (can be in any basis) from which the EKO was computed $-$ this way we don't rely on much information from the grids to construct the EKOs.
Does this make sense?
Extra-changes included in this PR (not fully related):
get-subgrids.cppexamplesTODO:
LHCB_WP_7TEV_opt.pineappl.lz4grid for the test which currently is available at:EKO_LHCB_WP_7TEV.txtwhich contains the flattened evolution operators