diff --git a/.github/actions/cache-test-data/action.yml b/.github/actions/cache-test-data/action.yml new file mode 100644 index 00000000..be1504ca --- /dev/null +++ b/.github/actions/cache-test-data/action.yml @@ -0,0 +1,19 @@ +# .github/actions/cache-test-data/action.yml + +name: Cache Test Data +description: Caches and downloads test data +runs: + using: "composite" + steps: + - name: Cache test data + id: cache + uses: actions/cache@v4 + with: + path: test-data + key: test-data-v21 + - name: Download test data if cache miss + if: steps.cache.outputs.cache-hit != 'true' + run: | + cd maintainer + ./download-test-data.sh + shell: bash diff --git a/.github/workflows/capi.yaml b/.github/workflows/capi.yaml index 65f35070..daa60ecf 100644 --- a/.github/workflows/capi.yaml +++ b/.github/workflows/capi.yaml @@ -23,6 +23,9 @@ jobs: cargo cinstall --verbose --prefix=/usr/local/ --libdir=/usr/local/lib ldconfig + - name: Get test data + uses: ./.github/actions/cache-test-data + - name: Test C++ examples run: | cd examples/cpp diff --git a/.github/workflows/python.yml b/.github/workflows/python.yml index 5e4806c0..6df4bfc8 100644 --- a/.github/workflows/python.yml +++ b/.github/workflows/python.yml @@ -23,19 +23,8 @@ jobs: - name: Set up Rust run: | rustup component add llvm-tools - - name: Get test data - id: cache-test-data - uses: actions/cache@v4 - with: - path: test-data - key: test-data-v20 - - name: Download test data - if: steps.cache-test-data.outputs.cache-hit != 'true' - run: | - cd maintainer - ./download-test-data.sh - + uses: ./.github/actions/cache-test-data - name: Test env: # `-C link-dead-code` is needed to prevent 'warning: XX functions have mismatched data' warnings diff --git a/.github/workflows/rust.yml b/.github/workflows/rust.yml index 44953aa3..0f9e7528 100644 --- a/.github/workflows/rust.yml +++ b/.github/workflows/rust.yml @@ -21,16 +21,7 @@ jobs: - uses: actions/checkout@v4 - name: Get test data - id: cache-test-data - uses: actions/cache@v4 - with: - path: test-data - key: test-data-v20 - - name: Download test data - if: steps.cache-test-data.outputs.cache-hit != 'true' - run: | - cd maintainer - ./download-test-data.sh + uses: ./.github/actions/cache-test-data - name: Set RUSTDOCFLAGS run: | diff --git a/Cargo.lock b/Cargo.lock index 26ababf7..f45e9e33 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -990,6 +990,7 @@ name = "pineappl_capi" version = "1.0.0" dependencies = [ "itertools", + "ndarray", "pineappl 1.0.0", ] diff --git a/examples/cpp/Makefile b/examples/cpp/Makefile index aef4689b..b9f0406e 100644 --- a/examples/cpp/Makefile +++ b/examples/cpp/Makefile @@ -14,6 +14,8 @@ PROGRAMS = \ convolve-grid-deprecated \ convolve-grid \ get-subgrids \ + evolve-grid \ + evolve-grid-identity \ deprecated \ display-channels-deprecated \ display-channels \ @@ -51,6 +53,12 @@ deprecated: deprecated.cpp get-subgrids: get-subgrids.cpp $(CXX) $(CXXFLAGS) $< $(PINEAPPL_DEPS) -o $@ +evolve-grid: evolve-grid.cpp + $(CXX) $(CXXFLAGS) $< $(LHAPDF_DEPS) $(PINEAPPL_DEPS) -o $@ + +evolve-grid-identity: evolve-grid-identity.cpp + $(CXX) $(CXXFLAGS) $< $(LHAPDF_DEPS) $(PINEAPPL_DEPS) -o $@ + display-channels-deprecated: display-channels-deprecated.cpp $(CXX) $(CXXFLAGS) $< $(PINEAPPL_DEPS) -o $@ diff --git a/examples/cpp/evolve-grid-identity.cpp b/examples/cpp/evolve-grid-identity.cpp new file mode 100644 index 00000000..df22d424 --- /dev/null +++ b/examples/cpp/evolve-grid-identity.cpp @@ -0,0 +1,257 @@ +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +// NOTE: Uses the scale of the Grid as the starting scale such that we can use an IDENTITY EKO. +double FAC0 = 6456.44; + +/** @brief This struct can contain arbitrary parameters that need to be passed to Evolution + * Operator Callback (`generated_fake_ekos`). + */ +struct OperatorParams { + std::vector conv_types; +}; + +std::vector unravel_index(std::size_t flat_index, const std::vector& shape) { + std::size_t ndim = shape.size(); + std::vector coords(ndim); + + for (int i = ndim - 1; i >= 0; --i) { + coords[i] = flat_index % shape[i]; + flat_index /= shape[i]; + } + + return coords; +} + +extern "C" void generate_fake_ekos( + std::size_t op_index, + double /*fac1*/, + const int* /*pids_in*/, + const double* /*x_in*/, + const int* /*pids_out*/, + const double* /*x_out*/, + const std::size_t* eko_shape, + double* eko_buffer, + void* params_state +) { + // select the type of convolution based on the Operator index + OperatorParams* op_params = static_cast(params_state); + pineappl_conv_type _ = op_params->conv_types[op_index]; + + // NOTE: This has to work because the Evolution Operator is always 4D + std::vector shape(eko_shape, eko_shape + 4); + // Compute the length of the flattened shape by multiplying the entries + std::size_t flat_len = std::accumulate(shape.begin(), + shape.end(), 1, std::multiplies()); + for (std::size_t i = 0; i != flat_len; i++) { + std::vector coords = unravel_index(i, shape); + + double delta_ik = (coords[0] == coords[2]) ? 1.0 : 0.0; + double delta_jl = (coords[1] == coords[3]) ? 1.0 : 0.0; + + eko_buffer[i] = delta_ik * delta_jl; + } +} + +void print_results(std::vector dxsec_grid, std::vector dxsec_fktable) { + const int idx_width = 6; + const int num_width = 15; + const int dif_width = 15; + + // Print headers + std::cout << std::setw(idx_width) << "Bin" + << std::setw(num_width) << "Grid" + << std::setw(num_width) << "FkTable" + << std::setw(dif_width) << "reldiff" << std::endl; + + // Print dashed lines + std::cout << std::setw(idx_width) << std::string(idx_width - 2, '-') + << std::setw(num_width) << std::string(num_width - 2, '-') + << std::setw(num_width) << std::string(num_width - 2, '-') + << std::setw(dif_width) << std::string(dif_width - 2, '-') << std::endl; + + // Print the data + std::cout << std::scientific << std::setprecision(6); + for (size_t i = 0; i < dxsec_grid.size(); ++i) { + double reldiff = (dxsec_fktable[i] - dxsec_grid[i]) / dxsec_grid[i]; + std::cout << std::setw(idx_width) << i + << std::setw(num_width) << dxsec_grid[i] + << std::setw(num_width) << dxsec_fktable[i] + << std::setw(dif_width) << reldiff + << std::endl; + } +} + +int main() { + // TODO: How to get a Grid that can be evolved?? + std::string filename = "../../test-data/LHCB_WP_7TEV_opt.pineappl.lz4"; + + // disable LHAPDF banners to guarantee deterministic output + LHAPDF::setVerbosity(0); + std::string pdfset = "NNPDF31_nlo_as_0118_luxqed"; + auto pdf = std::unique_ptr(LHAPDF::mkPDF(pdfset, 0)); + + auto xfx = [](int32_t id, double x, double q2, void* pdf) { + return static_cast (pdf)->xfxQ2(id, x, q2); + }; + auto alphas = [](double q2, void* pdf) { + return static_cast (pdf)->alphasQ2(q2); + }; + + std::vector pdfs = {pdf.get(), pdf.get()}; + void** pdf_states = reinterpret_cast(pdfs.data()); + + // read the grid from a file + auto* grid = pineappl_grid_read(filename.c_str()); + + // Get the PID basis representation + pineappl_pid_basis pid_basis = pineappl_grid_pid_basis(grid); + assert(pid_basis == PINEAPPL_PID_BASIS_PDG); + + // Get the number of convolutions and their types + std::size_t n_convs = pineappl_grid_convolutions_len(grid); + std::vector conv_types(n_convs); + pineappl_grid_conv_types(grid, conv_types.data()); + + // Fill the vector of unique convolution types. If the Operators required for the Grid + // are the same, then it suffices to only pass ONE single Operator. + std::vector unique_convs; + for (std::size_t i = 0; i != n_convs; i++) { + pineappl_conv_type conv = conv_types[i]; + if (std::find(unique_convs.begin(), unique_convs.end(), conv) == unique_convs.end()) { + unique_convs.push_back(conv); + } + } + + // Get the shape of the evolve info objects + std::vector evinfo_shape(5); + // NOTE: The argument of `pineappl_grid_evolve_info_shape` must follow the following orders: + // - `grid`: PineAPPL Grid + // - `order_mask`: array of booleans to mask the order(s) to apply the Evolution to, + // `nullptr` selects all the orders + // - `evinfo_shape`: placeholder to store the shape of the Evolution Operator + pineappl_grid_evolve_info_shape(grid, nullptr, evinfo_shape.data()); + + // Get the values of the evolve info parameters. These contain, for example, the + // information on the `x`-grid and `PID` used to interpolate the Grid. + // NOTE: These are used to construct the Evolution Operator + std::vector fac1(evinfo_shape[0]); + std::vector frg1(evinfo_shape[1]); + std::vector pids_in(evinfo_shape[2]); + std::vector x_in(evinfo_shape[3]); + std::vector ren1(evinfo_shape[4]); + // NOTE: The argument of `pineappl_grid_evolve_info` must follow the following orders: + // - `grid`: PineAPPL Grid + // - `order_mask`: array of booleans to mask the order(s) to apply the Evolution to, + // `nullptr` selects all the orders + // The rest of the arguments are placeholders to store data + pineappl_grid_evolve_info(grid, nullptr, fac1.data(), + frg1.data(), pids_in.data(), x_in.data(), ren1.data()); + + // ------------------ Construct the Operator Info ------------------ + // The Operator Info is a vector with length `N_conv * N_Q2_slices` whose + // elements are `OperatorInfo` objects. + std::vector convtypes(unique_convs.size()); + std::vector opinfo_slices(unique_convs.size() * fac1.size()); + for (std::size_t i = 0; i != unique_convs.size(); i++) { + for (std::size_t j = 0; j != fac1.size(); j++) { + pineappl_operator_info opinfo = { + FAC0, // fac0 + fac1[j], // fac1 + pid_basis, + unique_convs[i], + }; + opinfo_slices[i * fac1.size() + j] = opinfo; + } + conv_types[i] = unique_convs[i]; + } + + // ------------------ Construct the Evolution Operator ------------------ + // Choose a different PID basis for the FK table + // std::vector pids_out = {-6, -5, -4, -3, -2, -1, 1, 2, 3, 4, 5, 6, 21, 22}; + std::vector pids_out = pids_in; + + // Construct the values of alphas table + std::vector alphas_table; + for (double q2 : ren1) { + double alpha = alphas(q2, pdf.get()); + alphas_table.push_back(alpha); + } + + // Construct the Parameters that will get passed to the Callback + OperatorParams* op_params = new OperatorParams; + op_params->conv_types = convtypes; + void* params = static_cast(op_params); + + std::vector 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 tensor_shape = {pids_in.size(), x_in.size(), pids_out.size(), x_in.size()}; + + // NOTE: The arguments of `pineappl_grid_evolve` must follow the following orders: + // - `grid`: PineAPPL Grid + // - `nb_slices`: the number of convolution(s)/Evolution Operator(s) required + // - `slices`: callback that returns the evolution operator(s) in slices + // - `operator_info`: operator info + // - `pids_in`: PIDs basis representation of the Grid + // - `x_in`: x-grid of the Grid + // - `pids_out`: PIDs basis representation of the FK table + // - `x_out`: x-grid of the FK table + // - `state`: parameters that get passed to `operator` + // - `order_mask`: array of booleans to mask the order(s) to apply the Evolution to, + // `nullptr` selects all the orders + // - `xi`: scale variation + // - `ren1`: values of the renormalization scales + // - `alphas_table`: values of alphas for each renormalization scales + // - `eko_shape`: shape of the evolution operators + pineappl_grid* fktable = pineappl_grid_evolve( + grid, // `grid` + unique_convs.size(), // `nb_slices` + generate_fake_ekos, // `slices` + opinfo_slices.data(), // `operator_info` + pids_in.data(), // `pids_in` + x_in.data(), // `x_in` + pids_out.data(), // `pids_out` + x_in.data(), // `x_out` + tensor_shape.data(), // `eko_shape` + params, // `state` + nullptr, // `order_mask` + xi.data(), // `xi` + ren1.data(), // `ren1` + alphas_table.data() // `alphas_table` + ); + + // ------------------ Compare Grid & FK after convolution ------------------ + // how many bins does this grid have? + std::size_t bins = pineappl_grid_bin_count(grid); + + // [ convolve the Grid ] + std::vector mu_scales = { 1.0, 1.0, 1.0 }; + std::vector dxsec_grid(bins); + pineappl_grid_convolve(grid, xfx, alphas, pdf_states, pdf.get(), + nullptr, nullptr, nullptr, 1, + mu_scales.data(), dxsec_grid.data()); + + // [ convolve the FK Table ] + std::vector dxsec_fktable(bins); + auto as_one = [](double /*q2*/, void* /*pdf*/) { return 1.0; }; + pineappl_grid_convolve(fktable, xfx, as_one, pdf_states, nullptr, + nullptr, nullptr, nullptr, 1, + mu_scales.data(), dxsec_fktable.data()); + + // Print the results + print_results(dxsec_grid, dxsec_fktable); + + pineappl_grid_write(fktable, "evolved-grid-identity.pineappl.lz4"); + + pineappl_grid_delete(grid); + pineappl_grid_delete(fktable); +} diff --git a/examples/cpp/evolve-grid.cpp b/examples/cpp/evolve-grid.cpp new file mode 100644 index 00000000..1ab863e5 --- /dev/null +++ b/examples/cpp/evolve-grid.cpp @@ -0,0 +1,313 @@ +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +// NOTE: Uses the scale of the Grid as the starting scale such that we can use an IDENTITY EKO. +double FAC0 = 2.7224999999999997; + +// x-grid nodes definition for the `in` and `out` +std::vector XGRID = { + 2.00000000000000e-07, + 3.03430476586795e-07, + 4.60350147489639e-07, + 6.98420853070036e-07, + 1.05960949591010e-06, + 1.60758549847081e-06, + 2.43894329289168e-06, + 3.70022720698550e-06, + 5.61375771693015e-06, + 8.51680667757335e-06, + 1.29210156907473e-05, + 1.96025050023917e-05, + 2.97384953722449e-05, + 4.51143839496404e-05, + 6.84374491896790e-05, + 1.03811729865769e-04, + 1.57456056008414e-04, + 2.38787829185619e-04, + 3.62054496381397e-04, + 5.48779532367080e-04, + 8.31406883648814e-04, + 1.25867971442728e-03, + 1.90346340228674e-03, + 2.87386758128175e-03, + 4.32850063882081e-03, + 6.49620619463380e-03, + 9.69915957404340e-03, + 1.43750685810901e-02, + 2.10891866837872e-02, + 3.05215840078289e-02, + 4.34149174170227e-02, + 6.04800287544474e-02, + 8.22812212620489e-02, + 1.09143757463307e-01, + 1.41120806444403e-01, + 1.78025660425694e-01, + 2.19504126500389e-01, + 2.65113704158282e-01, + 3.14387400769276e-01, + 3.66875318648224e-01, + 4.22166775358965e-01, + 4.79898902961025e-01, + 5.39757233788045e-01, + 6.01472197967335e-01, + 6.64813948247382e-01, + 7.29586844241431e-01, + 7.95624252292276e-01, + 8.62783932390611e-01, + 9.30944080871754e-01, + 1.00000000000000e+00 +}; + +// Particle PIDs for both `in` and `out` +std::vector PIDS = {- 22 , -6 , -5 , -4 , -3 , -2 , -1 , 21 , 1 , 2 , 3 , 4 , 5 , 6}; + +/** @brief This struct can contain arbitrary parameters that need to be passed to Evolution + * Operator Callback (`generated_fake_ekos`). + */ +struct OperatorParams { + std::vector conv_types; +}; + +std::vector unravel_index(std::size_t flat_index, const std::vector& shape) { + std::size_t ndim = shape.size(); + std::vector coords(ndim); + + for (int i = ndim - 1; i >= 0; --i) { + coords[i] = flat_index % shape[i]; + flat_index /= shape[i]; + } + + return coords; +} + +extern "C" void generate_fake_ekos( + std::size_t op_index, + double /*fac1*/, + const int* /*pids_in*/, + const double* /*x_in*/, + const int* /*pids_out*/, + const double* /*x_out*/, + const std::size_t* /*eko_shape*/, + double* eko_buffer, + void* params_state +) { + // select the type of convolution based on the Operator index + OperatorParams* op_params = static_cast(params_state); + pineappl_conv_type _ = op_params->conv_types[op_index]; + + std::ifstream input_file("../../test-data/EKO_LHCB_WP_7TEV.txt"); + double weight_value; + + std::size_t count = 0; + while (input_file >> weight_value) { + eko_buffer[count++] = weight_value; + } +} + +void print_results(std::vector dxsec_grid, std::vector dxsec_fktable) { + const int idx_width = 6; + const int num_width = 15; + const int dif_width = 15; + + // Print headers + std::cout << std::setw(idx_width) << "Bin" + << std::setw(num_width) << "Grid" + << std::setw(num_width) << "FkTable" + << std::setw(dif_width) << "reldiff" << std::endl; + + // Print dashed lines + std::cout << std::setw(idx_width) << std::string(idx_width - 2, '-') + << std::setw(num_width) << std::string(num_width - 2, '-') + << std::setw(num_width) << std::string(num_width - 2, '-') + << std::setw(dif_width) << std::string(dif_width - 2, '-') << std::endl; + + // Print the data + std::cout << std::scientific << std::setprecision(6); + for (size_t i = 0; i < dxsec_grid.size(); ++i) { + double reldiff = (dxsec_fktable[i] - dxsec_grid[i]) / dxsec_grid[i]; + std::cout << std::setw(idx_width) << i + << std::setw(num_width) << dxsec_grid[i] + << std::setw(num_width) << dxsec_fktable[i] + << std::setw(dif_width) << reldiff + << std::endl; + } +} + +int main() { + // TODO: How to get a Grid that can be evolved?? + std::string filename = "../../test-data/LHCB_WP_7TEV_opt.pineappl.lz4"; + + // disable LHAPDF banners to guarantee deterministic output + LHAPDF::setVerbosity(0); + std::string pdfset = "NNPDF31_nlo_as_0118_luxqed"; + auto pdf = std::unique_ptr(LHAPDF::mkPDF(pdfset, 0)); + + auto xfx = [](int32_t id, double x, double q2, void* pdf) { + return static_cast (pdf)->xfxQ2(id, x, q2); + }; + auto alphas = [](double q2, void* pdf) { + return static_cast (pdf)->alphasQ2(q2); + }; + + std::vector pdfs = {pdf.get(), pdf.get()}; + void** pdf_states = reinterpret_cast(pdfs.data()); + + // read the grid from a file + auto* grid = pineappl_grid_read(filename.c_str()); + + // Get the PID basis representation + pineappl_pid_basis pid_basis = pineappl_grid_pid_basis(grid); + assert(pid_basis == PINEAPPL_PID_BASIS_PDG); + + // Get the number of convolutions and their types + std::size_t n_convs = pineappl_grid_convolutions_len(grid); + std::vector conv_types(n_convs); + pineappl_grid_conv_types(grid, conv_types.data()); + + // Fill the vector of unique convolution types. If the Operators required for the Grid + // are the same, then it suffices to only pass ONE single Operator. + std::vector unique_convs; + for (std::size_t i = 0; i != n_convs; i++) { + pineappl_conv_type conv = conv_types[i]; + if (std::find(unique_convs.begin(), unique_convs.end(), conv) == unique_convs.end()) { + unique_convs.push_back(conv); + } + } + + // Get the shape of the evolve info objects + std::vector evinfo_shape(5); + // NOTE: The argument of `pineappl_grid_evolve_info_shape` must follow the following orders: + // - `grid`: PineAPPL Grid + // - `order_mask`: array of booleans to mask the order(s) to apply the Evolution to, + // `nullptr` selects all the orders + // - `evinfo_shape`: placeholder to store the shape of the Evolution Operator + pineappl_grid_evolve_info_shape(grid, nullptr, evinfo_shape.data()); + + // Get the values of the evolve info parameters. These contain, for example, the + // information on the `x`-grid and `PID` used to interpolate the Grid. + // NOTE: These are used to construct the Evolution Operator + std::vector fac1(evinfo_shape[0]); + std::vector frg1(evinfo_shape[1]); + std::vector _pids_in(evinfo_shape[2]); + std::vector _x_in(evinfo_shape[3]); + std::vector ren1(evinfo_shape[4]); + // NOTE: The argument of `pineappl_grid_evolve_info` must follow the following orders: + // - `grid`: PineAPPL Grid + // - `order_mask`: array of booleans to mask the order(s) to apply the Evolution to, + // `nullptr` selects all the orders + // The rest of the arguments are placeholders to store data + pineappl_grid_evolve_info(grid, nullptr, fac1.data(), + frg1.data(), _pids_in.data(), _x_in.data(), ren1.data()); + + // ------------------ Construct the Operator Info ------------------ + // The Operator Info is a vector with length `N_conv * N_Q2_slices` whose + // elements are `OperatorInfo` objects. + std::vector convtypes(unique_convs.size()); + std::vector opinfo_slices(unique_convs.size() * fac1.size()); + for (std::size_t i = 0; i != unique_convs.size(); i++) { + for (std::size_t j = 0; j != fac1.size(); j++) { + pineappl_operator_info opinfo = { + FAC0, // fac0 + fac1[j], // fac1 + pid_basis, + unique_convs[i], + }; + opinfo_slices[i * fac1.size() + j] = opinfo; + } + conv_types[i] = unique_convs[i]; + } + + // ------------------ Construct the Evolution Operator ------------------ + // Define the same PIDs for both `in` and `out` according to the EKO + std::vector pids_in = PIDS; + std::vector pids_out = PIDS; + + // Construct the values of alphas table + std::vector alphas_table; + for (double q2 : ren1) { + double alpha = alphas(q2, pdf.get()); + alphas_table.push_back(alpha); + } + + // Construct the Parameters that will get passed to the Callback + OperatorParams* op_params = new OperatorParams; + op_params->conv_types = convtypes; + void* params = static_cast(op_params); + + std::vector 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 tensor_shape = {pids_in.size(), XGRID.size(), pids_out.size(), XGRID.size()}; + + // NOTE: The arguments of `pineappl_grid_evolve` must follow the following orders: + // - `grid`: PineAPPL Grid + // - `nb_slices`: the number of convolution(s)/Evolution Operator(s) required + // - `slices`: callback that returns the evolution operator(s) in slices + // - `operator_info`: operator info + // - `pids_in`: PIDs basis representation of the Grid + // - `x_in`: x-grid of the Grid + // - `pids_out`: PIDs basis representation of the FK table + // - `x_out`: x-grid of the FK table + // - `state`: parameters that get passed to `operator` + // - `order_mask`: array of booleans to mask the order(s) to apply the Evolution to, + // `nullptr` selects all the orders + // - `xi`: scale variation + // - `ren1`: values of the renormalization scales + // - `alphas_table`: values of alphas for each renormalization scales + // - `eko_shape`: shape of the evolution operators + pineappl_grid* fktable = pineappl_grid_evolve( + grid, // `grid` + unique_convs.size(), // `nb_slices` + generate_fake_ekos, // `slices` + opinfo_slices.data(), // `operator_info` + pids_in.data(), // `pids_in` + XGRID.data(), // `x_in` + pids_out.data(), // `pids_out` + XGRID.data(), // `x_out` + tensor_shape.data(), // `eko_shape` + params, // `state` + nullptr, // `order_mask` + xi.data(), // `xi` + ren1.data(), // `ren1` + alphas_table.data() // `alphas_table` + ); + + // ------------------ Compare Grid & FK after convolution ------------------ + // how many bins does this grid have? + std::size_t bins = pineappl_grid_bin_count(grid); + + // [ convolve the Grid ] + std::vector mu_scales = { 1.0, 1.0, 1.0 }; + std::vector dxsec_grid(bins); + pineappl_grid_convolve(grid, xfx, alphas, pdf_states, pdf.get(), + nullptr, nullptr, nullptr, 1, + mu_scales.data(), dxsec_grid.data()); + + // [ convolve the FK Table ] + std::vector dxsec_fktable(bins); + auto as_one = [](double /*q2*/, void* /*pdf*/) { return 1.0; }; + pineappl_grid_convolve(fktable, xfx, as_one, pdf_states, nullptr, + nullptr, nullptr, nullptr, 1, + mu_scales.data(), dxsec_fktable.data()); + + // Print the results + print_results(dxsec_grid, dxsec_fktable); + + // write the unoptimised FK table into disk + pineappl_grid_write(fktable, "evolved-grid.pineappl.lz4"); + + // optimise the FK table and then write into disk + pineappl_fktable_optimize(fktable, PINEAPPL_FK_ASSUMPTIONS_NF3_SYM); + pineappl_grid_write(fktable, "evolved-grid-optimised.pineappl.lz4"); + + pineappl_grid_delete(grid); + pineappl_grid_delete(fktable); +} diff --git a/examples/cpp/get-subgrids.cpp b/examples/cpp/get-subgrids.cpp index b37f0db4..5cd1385d 100644 --- a/examples/cpp/get-subgrids.cpp +++ b/examples/cpp/get-subgrids.cpp @@ -5,11 +5,10 @@ #include #include #include +#include #include #include #include -#include - std::vector unravel_index(std::size_t flat_index, const std::vector& shape) { std::size_t ndim = shape.size(); @@ -23,13 +22,14 @@ std::vector unravel_index(std::size_t flat_index, const std::vector return coords; } -std::string coords_to_string(const std::vector& coords) { +template +std::string vector_to_string(const std::vector& coords) { std::ostringstream osstream; osstream << "("; for (std::size_t i = 0; i < coords.size(); ++i) { osstream << coords[i]; - if (i != coords.size() - 1) { + if (i != coords.size() - 1) { osstream << ", "; } } @@ -38,6 +38,41 @@ std::string coords_to_string(const std::vector& coords) { return osstream.str(); } +std::vector get_subgrid_array( + const pineappl_grid* grid, + std::vector subgrid_shape, + std::size_t bin, + std::size_t order, + std::size_t channel +) { + // Compute the length of the flattened shape by multiplying the entries + std::size_t flat_shape = std::accumulate(subgrid_shape.begin(), + subgrid_shape.end(), 1, std::multiplies()); + + // Extract the flattened subgrid values/weights + std::vector subgrid_array(flat_shape); + pineappl_grid_subgrid_array(grid, bin, order, channel, subgrid_array.data()); + + return subgrid_array; +} + +std::vector get_node_values( + const pineappl_grid* grid, + std::vector subgrid_shape, + std::size_t bin, + std::size_t order, + std::size_t channel +) { + // Compute the length of the flattend nodes + std::size_t nodes_size = std::accumulate(subgrid_shape.begin(), + subgrid_shape.end(), 0.0); + + // Extract the values of the nodes as a flattened array + std::vector node_values(nodes_size); + pineappl_grid_subgrid_node_values(grid, bin, order, channel, node_values.data()); + + return node_values; +} int main() { std::string filename = "drell-yan-rap-ll.pineappl.lz4"; @@ -54,35 +89,57 @@ int main() { std::size_t subgrid_dim = pineappl_grid_kinematics_len(grid); std::cout << std::right << std::setw(10) << "bin" << std::setw(10) << "sg idx" - << std::setw(6 * subgrid_dim) << "sg coordinates" << std::setw(16) - << "sg value" << "\n"; + << std::setw(6 * subgrid_dim) << "sg coordinates" + << std::setw(12 * subgrid_dim) << "node values" << std::setw(16) + << "weight value" << "\n"; std::cout << std::right << std::setw(10) << "---" << std::setw(10) << "------" - << std::setw(6 * subgrid_dim) << "--------------" << std::setw(16) - << "------------" << "\n"; + << std::setw(6 * subgrid_dim) << "--------------" + << std::setw(12 * subgrid_dim) << "--------------------------------" + << std::setw(16) << "------------" << "\n"; for (std::size_t b = 0; b < n_bins; ++b) { + // Extract the shape of the subgrids std::vector subgrid_shape(subgrid_dim); pineappl_grid_subgrid_shape(grid, b, order, channel, subgrid_shape.data()); // Check if the subgrid is not empty if (subgrid_shape[0] != 0) { - std::size_t flat_shape = std::accumulate(subgrid_shape.begin(), - subgrid_shape.end(), 1, std::multiplies()); - std::vector subgrid_array(flat_shape); - - pineappl_grid_subgrid_array(grid, b, order, channel, subgrid_array.data()); + std::vector subgrid_array = get_subgrid_array(grid, subgrid_shape, b, order, channel); + std::vector node_values = get_node_values(grid, subgrid_shape, b, order, channel); for (std::size_t index = 0; index < subgrid_array.size(); ++index) { if (subgrid_array[index] != 0) { // Unravel the index to recover the standard coordinates std::vector coords = unravel_index(index, subgrid_shape); - std::cout << std::right << std::setw(10) << b << std::setw(10) - << index << std::setw(6 * subgrid_dim) << coords_to_string(coords) - << std::setw(16) << subgrid_array[index] << "\n"; + // Store the values of the nodes in a vector. The vector therefore + // contains as elements {Scale, x1, x2, ..., xn} + std::vector node_values_index(coords.size()); + std::size_t start_index = 0; + for (std::size_t nd = 0; nd < coords.size(); ++nd) { + if (nd != 0) { start_index += subgrid_shape[nd - 1]; } + node_values_index[nd] = node_values[start_index + coords[nd]]; + } - // Compare to some reference value - if (b==0 && index==41020) { + std::cout << std::right << std::setw(10) << b << std::setw(10) + << index << std::setw(6 * subgrid_dim) + << vector_to_string(coords) << std::setw(12 * subgrid_dim) + << vector_to_string(node_values_index) << std::setw(16) + << subgrid_array[index] << "\n"; + + // Compare to some reference value. + if (b == 0 && index == 41020) { + // Check the unravelled index + assert(coords[0] == 16); + assert(coords[1] == 20); + assert(coords[2] == 20); + + // Check the values of the node entries. + assert(node_values_index[0] == 5442.30542919352900); // PyAPI: `subgrid.node_values[0][16]` + assert(node_values_index[1] == 0.03052158400782889); // PyAPI: `subgrid.node_values[1][20]` + assert(node_values_index[2] == 0.03052158400782889); // PyAPI: `subgrid.node_values[2][20]` + + // PyAPI: `grid.subgrid(0, 0, 0).to_array(subgrid.shape)[16][20][20]` assert(subgrid_array[index] == -4.936156925096015e-07); } break; diff --git a/examples/cpp/output b/examples/cpp/output index d5ad2d4d..5dfc4176 100644 --- a/examples/cpp/output +++ b/examples/cpp/output @@ -146,32 +146,52 @@ idx left right dsig/dx dx 21 2.1 2.2 3.507294e-02 0.1 22 2.2 2.3 1.976542e-02 0.1 23 2.3 2.4 5.590552e-03 0.1 - bin sg idx sg coordinates sg value - --- ------ -------------- ------------ - 0 41020 (16, 20, 20) -4.93616e-07 - 1 41020 (16, 20, 20) -4.5499e-08 - 2 40971 (16, 19, 21) -4.10008e-08 - 3 40971 (16, 19, 21) -2.78597e-07 - 4 40971 (16, 19, 21) -2.01924e-07 - 5 40922 (16, 18, 22) -3.39322e-10 - 6 40922 (16, 18, 22) -8.43358e-08 - 7 40922 (16, 18, 22) -2.99247e-07 - 8 40922 (16, 18, 22) -1.5117e-07 - 9 40872 (16, 17, 22) -9.95013e-11 - 10 40873 (16, 17, 23) -1.30578e-07 - 11 40873 (16, 17, 23) -2.55869e-07 - 12 40823 (16, 16, 23) -1.0823e-09 - 13 40823 (16, 16, 23) -4.22686e-09 - 14 40824 (16, 16, 24) -1.68149e-07 - 15 40774 (16, 15, 24) -7.58702e-10 - 16 40774 (16, 15, 24) -2.41887e-08 - 17 40774 (16, 15, 24) -8.83139e-09 - 18 40725 (16, 14, 25) -1.14932e-09 - 19 40725 (16, 14, 25) -3.03857e-08 - 20 40725 (16, 14, 25) -2.66414e-08 - 21 40675 (16, 13, 25) -3.60243e-10 - 22 40676 (16, 13, 26) -1.22785e-08 - 23 40626 (16, 12, 26) -7.04493e-11 + bin sg idx sg coordinates node values weight value + --- ------ -------------- -------------------------------- ------------ + 0 41020 (16, 20, 20) (5442.31, 0.0305216, 0.0305216) -4.93616e-07 + 1 41020 (16, 20, 20) (5442.31, 0.0305216, 0.0305216) -4.5499e-08 + 2 40971 (16, 19, 21) (5442.31, 0.0434149, 0.0210892) -4.10008e-08 + 3 40971 (16, 19, 21) (5442.31, 0.0434149, 0.0210892) -2.78597e-07 + 4 40971 (16, 19, 21) (5442.31, 0.0434149, 0.0210892) -2.01924e-07 + 5 40922 (16, 18, 22) (5442.31, 0.06048, 0.0143751) -3.39322e-10 + 6 40922 (16, 18, 22) (5442.31, 0.06048, 0.0143751) -8.43358e-08 + 7 40922 (16, 18, 22) (5442.31, 0.06048, 0.0143751) -2.99247e-07 + 8 40922 (16, 18, 22) (5442.31, 0.06048, 0.0143751) -1.5117e-07 + 9 40872 (16, 17, 22) (5442.31, 0.0822812, 0.0143751) -9.95013e-11 + 10 40873 (16, 17, 23) (5442.31, 0.0822812, 0.00969916) -1.30578e-07 + 11 40873 (16, 17, 23) (5442.31, 0.0822812, 0.00969916) -2.55869e-07 + 12 40823 (16, 16, 23) (5442.31, 0.109144, 0.00969916) -1.0823e-09 + 13 40823 (16, 16, 23) (5442.31, 0.109144, 0.00969916) -4.22686e-09 + 14 40824 (16, 16, 24) (5442.31, 0.109144, 0.00649621) -1.68149e-07 + 15 40774 (16, 15, 24) (5442.31, 0.141121, 0.00649621) -7.58702e-10 + 16 40774 (16, 15, 24) (5442.31, 0.141121, 0.00649621) -2.41887e-08 + 17 40774 (16, 15, 24) (5442.31, 0.141121, 0.00649621) -8.83139e-09 + 18 40725 (16, 14, 25) (5442.31, 0.178026, 0.0043285) -1.14932e-09 + 19 40725 (16, 14, 25) (5442.31, 0.178026, 0.0043285) -3.03857e-08 + 20 40725 (16, 14, 25) (5442.31, 0.178026, 0.0043285) -2.66414e-08 + 21 40675 (16, 13, 25) (5442.31, 0.219504, 0.0043285) -3.60243e-10 + 22 40676 (16, 13, 26) (5442.31, 0.219504, 0.00287387) -1.22785e-08 + 23 40626 (16, 12, 26) (5442.31, 0.265114, 0.00287387) -7.04493e-11 + Bin Grid FkTable reldiff + ---- ------------- ------------- ------------- + 0 7.545911e+02 7.559965e+02 1.862504e-03 + 1 6.902834e+02 6.921603e+02 2.718991e-03 + 2 6.002520e+02 6.024318e+02 3.631522e-03 + 3 4.855224e+02 4.877053e+02 4.496133e-03 + 4 3.619546e+02 3.638695e+02 5.290429e-03 + 5 2.458669e+02 2.473256e+02 5.933028e-03 + 6 1.158685e+02 1.166059e+02 6.363742e-03 + 7 2.751727e+01 2.765255e+01 4.916187e-03 + Bin Grid FkTable reldiff + ---- ------------- ------------- ------------- + 0 7.545911e+02 7.545911e+02 -3.762829e-08 + 1 6.902834e+02 6.902834e+02 -3.537025e-08 + 2 6.002520e+02 6.002520e+02 -3.293028e-08 + 3 4.855224e+02 4.855223e+02 -3.037598e-08 + 4 3.619546e+02 3.619545e+02 -2.782860e-08 + 5 2.458669e+02 2.458669e+02 -2.526536e-08 + 6 1.158685e+02 1.158685e+02 -2.182740e-08 + 7 2.751727e+01 2.751727e+01 -1.732084e-08 idx p-p c#0 l#0 p-p~ c#0 l# p-d c#0 l#0 p-d dx --- ------------ ----------- ------------- ------------ ------ 0 5.263109e-01 5.263109e-01 5.263109e-01 5.263109e-01 0.1 diff --git a/maintainer/download-test-data.sh b/maintainer/download-test-data.sh index 62c69348..fc036db1 100755 --- a/maintainer/download-test-data.sh +++ b/maintainer/download-test-data.sh @@ -35,6 +35,7 @@ files=( 'https://data.nnpdf.science/pineappl/test-data/LHCB_DY_8TEV.pineappl.lz4' 'https://data.nnpdf.science/pineappl/test-data/LHCB_DY_8TEV.tar' 'https://data.nnpdf.science/pineappl/test-data/LHCB_WP_7TEV.tar' + 'https://data.nnpdf.science/pineappl/test-data/EKO_LHCB_WP_7TEV.txt' 'https://data.nnpdf.science/pineappl/test-data/LHCB_WP_7TEV_old.pineappl.lz4' 'https://data.nnpdf.science/pineappl/test-data/LHCB_WP_7TEV_opt.pineappl.lz4' 'https://data.nnpdf.science/pineappl/test-data/LHCB_WP_7TEV_v2.tar' diff --git a/pineappl/src/fk_table.rs b/pineappl/src/fk_table.rs index 437fcc4f..4a9c6350 100644 --- a/pineappl/src/fk_table.rs +++ b/pineappl/src/fk_table.rs @@ -31,6 +31,7 @@ pub struct FkTable { /// tables are typically stored at very small `Q2 = Q0`, the PDFs `f(x,Q0)` of heavy quarks are /// typically set to zero at this scale or set to the same value as their anti-quark PDF. This is /// used to optimize the size of FK tables. +#[repr(C)] #[derive(Debug, Clone, Copy, Eq, PartialEq)] pub enum FkAssumptions { /// All quark PDFs are non-zero at the FK table scale and completely independent. diff --git a/pineappl_capi/Cargo.toml b/pineappl_capi/Cargo.toml index dfde3664..d0f2ff74 100644 --- a/pineappl_capi/Cargo.toml +++ b/pineappl_capi/Cargo.toml @@ -21,6 +21,7 @@ workspace = true [dependencies] pineappl = { path = "../pineappl", version = "=1.0.0" } itertools = "0.10.1" +ndarray = "0.15.4" [features] capi = [] diff --git a/pineappl_capi/cbindgen.toml b/pineappl_capi/cbindgen.toml index 7d3c6510..61ba6a88 100644 --- a/pineappl_capi/cbindgen.toml +++ b/pineappl_capi/cbindgen.toml @@ -38,6 +38,7 @@ rename_variants = "ScreamingSnakeCase" "Channels" = "pineappl_channels" "Conv" = "pineappl_conv" "ConvType" = "pineappl_conv_type" +"FkAssumptions" = "pineappl_fk_assumptions" "Grid" = "pineappl_grid" "GridOptFlags" = "pineappl_gof" "Interp" = "pineappl_interp" @@ -46,6 +47,8 @@ rename_variants = "ScreamingSnakeCase" "Kinematics" = "pineappl_kinematics" "Lumi" = "pineappl_lumi" "Map" = "pineappl_map" +"OperatorCallback" = "pineappl_operator_callback" +"OperatorInfo" = "pineappl_operator_info" "PidBasis" = "pineappl_pid_basis" "ReweightMeth" = "pineappl_reweight_meth" "ScaleFuncForm" = "pineappl_scale_func_form" diff --git a/pineappl_capi/src/lib.rs b/pineappl_capi/src/lib.rs index d36f9a88..ca82e8b5 100644 --- a/pineappl_capi/src/lib.rs +++ b/pineappl_capi/src/lib.rs @@ -55,9 +55,12 @@ //! //! [translation tables]: https://github.com/eqrion/cbindgen/blob/master/docs.md#std-types -use itertools::izip; +use itertools::{izip, Itertools}; +use ndarray::{Array4, CowArray}; use pineappl::boc::{Bin, BinsWithFillLimits, Channel, Kinematics, Order, ScaleFuncForm, Scales}; use pineappl::convolutions::{Conv, ConvType, ConvolutionCache}; +use pineappl::evolution::{AlphasTable, OperatorSliceInfo}; +use pineappl::fk_table::{FkAssumptions, FkTable}; use pineappl::grid::{Grid, GridOptFlags}; use pineappl::interpolation::{Interp as InterpMain, InterpMeth, Map, ReweightMeth}; use pineappl::packed_array::ravel_multi_index; @@ -1485,6 +1488,16 @@ pub struct Interp { pub interp_meth: InterpMeth, } +/// Type for defining the Operator info. +#[repr(C)] +#[derive(Debug)] +pub struct OperatorInfo { + fac0: f64, + fac1: f64, + pid_basis: PidBasis, + conv_type: ConvType, +} + /// An exact duplicate of `pineappl_lumi_new` to make naming (lumi -> channel) consistent. /// should be deleted using `pineappl_channels_delete`. #[no_mangle] @@ -1846,7 +1859,7 @@ pub unsafe extern "C" fn pineappl_grid_order_params2(grid: *const Grid, order_pa } } -/// A generalization of the convolution function. +/// A generalization of the convolution function for Grids. /// /// # Safety /// @@ -1926,6 +1939,24 @@ pub unsafe extern "C" fn pineappl_grid_convolve( )); } +/// Get the type of convolutions for this Grid. +/// +/// # Safety +/// +/// TODO +#[no_mangle] +pub unsafe extern "C" fn pineappl_grid_conv_types(grid: *const Grid, conv_types: *mut ConvType) { + let grid = unsafe { &*grid }; + let conv_types = unsafe { slice::from_raw_parts_mut(conv_types, grid.convolutions().len()) }; + let convs_array = grid + .convolutions() + .iter() + .map(Conv::conv_type) + .collect_vec(); + + conv_types.copy_from_slice(&convs_array); +} + /// Get the number of convolutions for a given Grid. /// /// # Safety @@ -1978,10 +2009,31 @@ pub unsafe extern "C" fn pineappl_grid_subgrid_shape( }; let shape = unsafe { slice::from_raw_parts_mut(shape, grid.kinematics().len()) }; - shape.copy_from_slice(&subgrid_shape); + shape.copy_from_slice(subgrid_shape); +} + +/// Get the nodes of the subgrid. +/// +/// # Safety +/// +#[no_mangle] +pub unsafe extern "C" fn pineappl_grid_subgrid_node_values( + grid: *const Grid, + bin: usize, + order: usize, + channel: usize, + node_values: *mut f64, +) { + let grid = unsafe { &*grid }; + let subgrid = &grid.subgrids()[[order, bin, channel]]; + + let flatten_subgrid: Vec = subgrid.node_values().into_iter().flatten().collect(); + let node_values = unsafe { slice::from_raw_parts_mut(node_values, flatten_subgrid.len()) }; + + node_values.copy_from_slice(flatten_subgrid.as_slice()); } -/// Get the subgrid for a given bin, channel, and order +/// Get the subgrid for a given bin, channel, and order. /// /// # Safety /// @@ -2007,8 +2059,247 @@ pub unsafe extern "C" fn pineappl_grid_subgrid_array( unsafe { slice::from_raw_parts_mut(subgrid_array, shape.iter().product()) }; for (index, value) in subgrid.indexed_iter() { - let ravel_index = ravel_multi_index(index.as_slice(), &shape); + let ravel_index = ravel_multi_index(index.as_slice(), shape); subgrid_array[ravel_index] = value; } } } + +/// Get the shape of the objects represented in the evolve info. +/// +/// # Safety +/// +/// If `grid` does not point to a valid `Grid` object, for example when `grid` is the null pointer, +/// this function is not safe to call. +#[no_mangle] +pub unsafe extern "C" fn pineappl_grid_evolve_info_shape( + grid: *const Grid, + order_mask: *const bool, + shape_info: *mut usize, +) { + let grid = unsafe { &*grid }; + + let order_mask = if order_mask.is_null() { + &[] + } else { + unsafe { slice::from_raw_parts(order_mask, grid.orders().len()) } + }; + + let shape_info = unsafe { slice::from_raw_parts_mut(shape_info, 5) }; + + let evolve_info = grid.evolve_info(order_mask); + let evinfo_shapes = [ + evolve_info.fac1.len(), + evolve_info.frg1.len(), + evolve_info.pids1.len(), + evolve_info.x1.len(), + evolve_info.ren1.len(), + ]; + + shape_info.copy_from_slice(&evinfo_shapes); +} + +/// Get the values of the parameters contained in evolve info. +/// +/// # Safety +/// +/// If `grid` does not point to a valid `Grid` object, for example when `grid` is the null pointer, +/// this function is not safe to call. +#[no_mangle] +pub unsafe extern "C" fn pineappl_grid_evolve_info( + grid: *const Grid, + order_mask: *const bool, + fac1: *mut f64, + frg1: *mut f64, + pids1: *mut i32, + x1: *mut f64, + ren1: *mut f64, +) { + let grid = unsafe { &*grid }; + + let order_mask = if order_mask.is_null() { + &[] + } else { + unsafe { slice::from_raw_parts(order_mask, grid.orders().len()) } + }; + + let evolve_info = grid.evolve_info(order_mask); + let fac1 = unsafe { slice::from_raw_parts_mut(fac1, evolve_info.fac1.len()) }; + let frg1 = unsafe { slice::from_raw_parts_mut(frg1, evolve_info.frg1.len()) }; + let pids1 = unsafe { slice::from_raw_parts_mut(pids1, evolve_info.pids1.len()) }; + let x1 = unsafe { slice::from_raw_parts_mut(x1, evolve_info.x1.len()) }; + let ren1 = unsafe { slice::from_raw_parts_mut(ren1, evolve_info.ren1.len()) }; + + fac1.copy_from_slice(&grid.evolve_info(order_mask).fac1); + frg1.copy_from_slice(&grid.evolve_info(order_mask).frg1); + pids1.copy_from_slice(&grid.evolve_info(order_mask).pids1); + x1.copy_from_slice(&grid.evolve_info(order_mask).x1); + ren1.copy_from_slice(&grid.evolve_info(order_mask).ren1); +} + +/// Type alias for the operator callback +pub type OperatorCallback = unsafe extern "C" fn( + usize, // index which selects Evolution parameters + f64, // fac1 + *const i32, // `pids_in` + *const f64, // `x_in` + *const i32, // `pids_out` + *const f64, // `x_out` + *const usize, // shape of the EKO + *mut f64, // Evolution Operator data buffer + *mut c_void, // Callable state of parameters +); + +/// Evolve a grid with an evolution operator and dump the resulting FK table. +/// +/// # Arguments +/// +/// * `grid` - A `Grid` object +/// * `op_info` - An array of `OperatorInfo` objects containing the information about the evolution. +/// * `operator` - A callack that returns the evolution operator. +/// * `max_orders` - The maximum QCD and EW orders `(αs, α)`. +/// * `params_state` - Parameters that get passed to `operator`. +/// * `x_in` - The x-grid that defines the Grid. +/// * `x_out` - The x-grid that will define the evolved Grid. +/// * `pids_in` - The list of PID values that defines the Grid. +/// * `pids_out` - The list of PID values that will define the evolved Grid. +/// * `eko_shape` - The shape of the evolution operator. +/// * `xi` - The values that defines that scale variations. +/// * `ren` - An array containing the values of the renormalization scale variation. +/// * `alphas` - An array containing the values of `αs`. It must have the same size as `ren1`. +/// +/// # Safety +/// +/// This function is not safe to call if: (a) the `grid` does not point to a valid `Grid` object or +/// is a null pointer; (b) the `op_info` and `operators` objects do not have the expected lengths, +/// (c) the shape of `eko_shape` is different from the actual size of `pids_out`, `x_out`, `pids_in`, +/// `x_in`. +/// +/// # Panics +/// +/// This function might panic if the either the `op_info` and/or `operators` are/is incompatible +/// with the `Grid`. +#[no_mangle] +pub unsafe extern "C" fn pineappl_grid_evolve( + grid: *const Grid, + nb_slices: usize, + slices: OperatorCallback, + operator_info: *const OperatorInfo, + pids_in: *const i32, + x_in: *const f64, + pids_out: *const i32, + x_out: *const f64, + eko_shape: *const usize, + state: *mut c_void, + order_mask: *const bool, + xi: *const f64, + ren1: *const f64, + alphas: *const f64, +) -> Box { + let grid = unsafe { &*grid }; + + let order_mask = if order_mask.is_null() { + &[] + } else { + unsafe { slice::from_raw_parts(order_mask, grid.orders().len()) } + }; + + let eko_shape = unsafe { slice::from_raw_parts(eko_shape, 4) }; + let pids_in = unsafe { slice::from_raw_parts(pids_in, eko_shape[0]) }; + let x_in = unsafe { slice::from_raw_parts(x_in, eko_shape[1]) }; + let pids_out = unsafe { slice::from_raw_parts(pids_out, eko_shape[2]) }; + let x_out = unsafe { slice::from_raw_parts(x_out, eko_shape[3]) }; + + let evolve_info = grid.evolve_info(order_mask); + let ren1_len = evolve_info.ren1.len(); + let ren1 = unsafe { slice::from_raw_parts(ren1, ren1_len) }; + let alphas = unsafe { slice::from_raw_parts(alphas, ren1_len) }; + let xi = unsafe { slice::from_raw_parts(xi, 3) }; + + let operator_info = unsafe { + slice::from_raw_parts(operator_info, nb_slices * evolve_info.fac1.len()) + .chunks_exact(evolve_info.fac1.len()) + }; + + let shape: (usize, usize, usize, usize) = <[usize; 4]>::try_from(eko_shape) + // UNWRAP: guaranteed to work since `eko_shape` is exactly 4 elements long + .unwrap() + .into(); + + let op_slices = operator_info + .map(|op_infos| { + op_infos.iter().enumerate().map(|(op_index, op_info)| { + let operator_slice_info = OperatorSliceInfo { + pid_basis: op_info.pid_basis, + fac0: op_info.fac0, + pids0: pids_out.to_vec(), + x0: x_out.to_vec(), + fac1: op_info.fac1, + pids1: pids_in.to_vec(), + x1: x_in.to_vec(), + conv_type: op_info.conv_type, + }; + + let mut array = CowArray::from(Array4::zeros(shape)); + + unsafe { + slices( + op_index, + op_info.fac1, + pids_in.as_ptr(), + x_in.as_ptr(), + pids_out.as_ptr(), + x_out.as_ptr(), + eko_shape.as_ptr(), + array + .as_slice_mut() + // UNWRAP: `array` is by construction contiguous and in standard order + .unwrap() + .as_mut_ptr(), + state, + ); + } + + // we specify an arbitrary error type since we don't return an error anywhere + Ok::<_, std::io::Error>((operator_slice_info, array)) + }) + }) + .collect(); + + let fk_table = grid + .evolve( + op_slices, + order_mask, + (xi[0], xi[1], xi[2]), + &AlphasTable { + ren1: ren1.to_vec(), + alphas: alphas.to_vec(), + }, + ) + // UNWRAP: error handling in the CAPI is to abort + .unwrap(); + + Box::new(fk_table.into_grid()) +} + +/// Optimize the size of this FK-table by removing heavy quark flavors assumed to be zero. +/// +/// # Safety +/// +/// If `grid` does not point to a valid `Grid` object, for example when `grid` is the null pointer, +/// this function is not safe to call. +/// +/// # Panics +/// +/// This function panics if `grid` is not an FK table-like object. +#[no_mangle] +pub unsafe extern "C" fn pineappl_fktable_optimize(grid: *mut Grid, assumptions: FkAssumptions) { + let grid = unsafe { &mut *grid }; + // SAFETY: this code has been copied from the `take_mut` crate + let read_grid = unsafe { std::ptr::read(grid) }; + let mut fktable = FkTable::try_from(read_grid) + // UNWRAP: error handling in the CAPI is to abort + .unwrap(); + fktable.optimize(assumptions); + unsafe { std::ptr::write(grid, fktable.into_grid()) }; +}