On-device N-dimensional Linear Interpolation#
Introduction#
When there are cases where analytical function evaluations from particle properties are not tenable (eg. TRIM data), pre-calculated grids of function evaluations at a limited set of discrete points are usually the method used to approach the problem. Each dimensional component of the coordinates of these points corresponds to a value of a particle property, such as the locally sampled fluid density or locally sampled fluid temperature (“locally sampled” here just means sampled at the particle location). A pre-requisite for the pre-calculated grid is that each axis must be sorted to be either monotonically increasing or decreasing.
The method available in VANTAGE-Reactions to retrieve arbitrary function evaluations from a pre-calculated grid is to use a series of nested hypercube contractions via 1D linear interpolations. The general method is outlined in this paper Murman_S_apnum_jun13.
A pre-requisite for the axes of the grid is that they must be monotonic for the hypercube construction to work correctly.
This method can be extended to any arbitrary number of dimensions and scales (in terms of linear interpolation evaluations) as \(O(2^n - 1)\) where \(n\) is the number of dimensions. For further details on the underlying maths: Overview of the interpolation method
Limitations#
Whilst being quite powerful, this method does have some limitations. Chiefly, it’s quite sensitive to the sparseness of the pre-calculated grid. If the grid is too sparse then the linear interpolation will introduce inaccuracies that effectively come from the approximation of a constant linear gradient between points in the hypercube.
Grid construction#
There is an expectation that the ReactionDataBase-derived object that is an input of InterpolateData has the correct construction in the sense that the on-device calc_data() has the correct form of arguments and outputs. Specifically the form of calc_data() should follow the second definition in ReactionDataBaseOnDevice.
The calc_data() should provide a function evaluation at a given point, InterpolateData will use values from multiple invocations of calc_data() to construct the hypercube which will have function evaluations at each vertex.
Inputs for InterpolateData#
For flexibility, whilst there are guide-rails for what’s expected from the input ReactionDataBase, full manual configuration and usage is possible.
But for 2 common cases, there are some helper classes available, that aim to reduce the friction in constructing the input object.
Cartesian case#
In the case that the pre-calculated grid is a simple cartesian grid of size-1 function evaluations, CartesianGridData is available. An object of this type can be pre-constructed and then simply passed as an argument for the construction of a InterpolateData object.
If you already have the pre-requisite flattened vectors then it’s possible to directly construct CartesianGridData. But there’s an additional helper class that can aid in producing the flattened vectors. GridDescriptor only requires an array of vectors (individual vectors that contain coordinates for each axis) and a lambda that returns function evaluations given a set of coordinates. The constructed GridDescriptor object can then be used to construct CartesianGridData.
For further details: Cartesian case
TRIM case#
Similarly, if the pre-calculated grid is a cartesian grid of multi-dimensional function evaluations (specifically in the format of TRIM tables), then TrimEvalData is available and functions much the same way as CartesianGridData.
If the correct template arguments and standard arguments are given to GridDescriptor then it can also provide the necessary flattened vectors for TrimEvalData.
There are some further details (including the expected form of the lambda provided to GridDescriptor) in: TRIM case
Usage (size-1 function)#
The implementation of the interpolation is in InterpolateData which inherits from CompositeData. This is due to how the access to the pre-calculated grid is managed. For flexibility, rather than passing a full grid to the host-side InterpolateData object, a ReactionDataBase-derived object is passed, which will have a calc_data() that can retrieve values at coordinates that correspond to particle property values as described in Grid construction. This can be a simple retrieval from a look-up table or something more complicated like doing a limited set of calculations or even incorporating samples from a random distribution.
A few pre-requisites for the construction of InterpolateData:
output_ndim, asize_ttemplate parameter that specifies the number of outputs of the grid function. In the case of a single-valued function this would be 1.
interp_ndim, asize_ttemplate parameter that specifies how many dimensions out of the total number of dimensions of the grid to interpolate.
non_interp_ndim, asize_ttemplate parameter that is effectively the inverse ofinterp_ndim(default value is 0).
dims_vec, astd::vector<size_t>that contains the lengths of the extent of each dimension in the grid.
coords_vec, astd::vector<REAL>that is a vector containing a concatenated list of all of values of each dimension in the same order asdims_vec. For example ifdims_vec = {3, 2}thencoords_vec = {dim0_val_0, dim0_val_1, dim0_val_2, dim1_val_0, dim1_val_1}.
interp_indices, astd::array<REAL, interp_ndim>that just specifies which of the dimensions of the grid are to be interpolated.
extrapolation_type, an enum specifying the choice of how to handle extrapolation as explained in Extrapolation.
There are a few narrow interfaces for InterpolateData but in the example below, the widest interface is shown:
inline void interpolation_example(NP::ParticleGroupSharedPtr particle_group) {
// Number of dimensions of the pre-calculated grid
static constexpr int ndim = 3;
// Example coordinates for dimension 0.
std::vector<REAL> dim0_range = {1.0e+18, 2.0e+18, 3.0e+18, 4.0e+18,
5.0e+18, 6.0e+18, 7.0e+18, 8.0e+18};
// Example coordinates for dimension 1.
std::vector<REAL> dim1_range = {
1.00000000e+01, 2.78255940e+01, 7.74263683e+01, 2.15443469e+02,
5.99484250e+02, 1.66810054e+03, 4.64158883e+03, 1.29154967e+04,
3.59381366e+04, 1.00000000e+05};
// Example coordinates for dimension 2.
std::vector<REAL> dim2_range = {-98.5, -94., -86.5, -76., -62.5,
-46., -26.5, -4., 21.5, 50.,
81.5, 116., 153.5, 194., 237.5};
// Example lambda defining the values of the grid function at each coordinate.
static constexpr auto grid_func_lambda =
[](const std::array<REAL, ndim> &vals) {
return (vals[0] * vals[1] * vals[2]);
};
// Helper class to construct pre-requisites for CartesianGridData.
auto grid_descriptor = GridDescriptor<ndim>(
{dim0_range, dim1_range, dim2_range}, grid_func_lambda);
auto dims_vec = grid_descriptor.get_interp_dims();
auto coords_vec = grid_descriptor.get_flat_coords();
// Helper class that provides a calc_data() that retrieves values from the
// pre-calculated grid in grid_descriptor.
auto grid_func_data =
CartesianGridData<ndim>(grid_descriptor, particle_group->sycl_target);
// Construction of InterpolateData. There is a calc_data() within the
// on-device object associated with this class which may be used directly if
// the correct input_array is given. But it is meant to be used in a pipeline
// with results from another ReactionData-derived object acting as inputs for
// InterpolateData.
auto interpolate_data = InterpolateData<1, ndim, decltype(grid_func_data)>(
dims_vec, coords_vec, particle_group->sycl_target, grid_func_data);
// If a different extrapolation other than continue_linear is needed, then
// just add either:
// - ExtrapolationType::clamp_to_edge
// - ExtrapolationType::clamp_to_zero
// as the last argument when constructing InterpolateData.
// Example of construction of a ReactionData-derived object which has an
// output that matches the type of the expected input of the calc_data() from
// InterpolateDataOnDevice.
auto prop0_extract = extract<1>("PROP0");
auto prop1_extract = extract<1>("PROP1");
auto prop2_extract = extract<1>("PROP2");
auto concatenator =
ConcatenatorData(prop0_extract, prop1_extract, prop2_extract);
// Pipeline that handles the pass-through of values.
auto pipeline = pipe(concatenator, interpolate_data);
// Wrapping in a DataCalculator allows the extraction of values from
// interpolate_data via a NP::NDLocalArray buffer.
auto data_calc = DataCalculator(pipeline);
return;
}
Usage (multi-dimensional function)#
For multi-valued functions the process is similar but has a few key differences. Firstly, the helper functions are handled differently but also the constructor for InterpolateData requires an extra argument that specifies which dimensions of the grid that need to be interpolated. This is due to the fact that the “grid” for multi-valued functions is treated as having (grid_dimensions + n_function_outputs). For example, with TRIM data, if the tables were assigned to coordinates with dimensionality of 2 but the access convention for the table required 3 numbers (and also outputted 3 numbers) then the grid’s dimensionality is 5.
An example is shown here:
#include <random>
#include <tuple>
inline void
trim_interpolation_example(NP::ParticleGroupSharedPtr particle_group) {
// Number of dimensions of the pre-calculated grid
static constexpr int interp_ndim = 2;
// Number of dimensions associated with the TRIM tables
static constexpr int trim_ndim = 3;
static constexpr int trim_dim0 = 5;
static constexpr int trim_dim1 = 5;
static constexpr int trim_dim2 = 5;
static constexpr auto trim_dims_arr =
std::array<size_t, trim_ndim>{trim_dim0, trim_dim1, trim_dim2};
// Example coordinates for dimension 0.
std::vector<REAL> dim0_range = {1.0e+18, 2.0e+18, 3.0e+18, 4.0e+18,
5.0e+18, 6.0e+18, 7.0e+18, 8.0e+18};
// Example coordinates for dimension 1.
std::vector<REAL> dim1_range = {
1.00000000e+01, 2.78255940e+01, 7.74263683e+01, 2.15443469e+02,
5.99484250e+02, 1.66810054e+03, 4.64158883e+03, 1.29154967e+04,
3.59381366e+04, 1.00000000e+05};
// Example lambda for the TRIM tables.
static constexpr auto trim_grid_func_lambda =
[](const std::array<REAL, interp_ndim> &vals) {
// This calculates the correct length for the output.
// \sum_{i = 0}^{n}(\prod_{j=0}^{i} trim_dims_arr[j])
static constexpr size_t trim_table_length = [&] {
INT s = 0;
INT p = 1;
for (INT x : trim_dims_arr) {
p *= x;
s += p;
}
return s;
}();
std::array<REAL, trim_table_length> result;
for (size_t i = 0; i < trim_table_length; i++) {
result[i] = std::pow((i + 50), 3.0);
}
return result;
};
// Helper class to construct pre-requisites for TrimEvalData.
auto grid_descriptor = GridDescriptor<interp_ndim, trim_ndim>(
{dim0_range, dim1_range}, {trim_dim0, trim_dim1, trim_dim2},
trim_grid_func_lambda);
auto dims_vec = grid_descriptor.get_interp_dims();
auto coords_vec = grid_descriptor.get_flat_coords();
// Helper class that provides a calc_data() that retrieves values from the
// pre-calculated grid in grid_descriptor.
auto grid_func_data = TrimEvalData<interp_ndim + trim_ndim>(
grid_descriptor, particle_group->sycl_target);
// Given that the "grid" will have (interp_ndim + trim_ndim) dimensions, it's
// necessary to specify the indices of the dimensions that are to be
// interpolated.
std::array<size_t, interp_ndim> interp_indices = {0, 1};
// Construction of InterpolateData. There is a calc_data() within the
// on-device object associated with this class which may be used directly if
// the correct input_array is given. But it is meant to be used in a pipeline
// with results from another ReactionData-derived object acting as inputs for
// InterpolateData.
auto interpolate_data =
InterpolateData<trim_ndim, interp_ndim, decltype(grid_func_data),
trim_ndim>(dims_vec, coords_vec, interp_indices,
particle_group->sycl_target, grid_func_data);
// If a different extrapolation other than continue_linear is needed, then
// just add either:
// - ExtrapolationType::clamp_to_edge
// - ExtrapolationType::clamp_to_zero
// as the last argument when constructing InterpolateData.
// Example of construction of a ReactionData-derived object which has an
// output that matches the type of the expected input of the calc_data() from
// InterpolateDataOnDevice.
auto props_extract = extract<interp_ndim>("PROPS");
const int rank = particle_group->sycl_target->comm_pair.rank_parent;
auto rng = std::mt19937(52234126 + rank);
std::uniform_real_distribution<REAL> uniform_dist_2(0.0, 1.0);
auto rng_lambda = [&]() -> REAL {
REAL rng_sample = 0.0;
do {
rng_sample = uniform_dist_2(rng);
} while (rng_sample == 0.0);
return rng_sample;
};
auto trim_rng_kernel = NP::host_per_particle_block_rng<REAL>(rng_lambda, 1);
auto trim_sampler = SamplerData(trim_rng_kernel);
// This is hard-coded here to avoid bloated general implementation for
// arbitrary number of samplers. (quite easy in C++20)
auto trim_sampler_concat =
ConcatenatorData(trim_sampler, trim_sampler, trim_sampler);
auto concatenator = ConcatenatorData(props_extract, trim_sampler_concat);
// Pipeline that handles the pass-through of values.
auto pipeline = pipe(concatenator, interpolate_data);
// Wrapping in a DataCalculator allows the extraction of values from
// interpolate_data via a NP::NDLocalArray buffer.
auto data_calc = DataCalculator(pipeline);
return;
}
Extrapolation#
If the desired point has coordinates either partially or fully outside of the valid range of any/all dimensions of the pre-calculated grid then a choice must be made as to how to handle this. The options for this in VANTAGE-Reactions are:
default: Continue with the linear interpolation using the gradient and intercept calculation from the edge of the grid.
ExtrapolationType::clamp_to_zero: Clamp the function evaluation to zero if any dimensional components of the coordinates of the desired point are outside the grid.
ExtrapolationType::clamp_to_edge: Clamp the function evaluation to be as if it were calculated from the grid points at the edge of the grid (no gradient or intercept continuation).
Integration with Reactions#
The wrapped interpolation pipeline is effectively a composite ReactionData object whose generated values can be treated as either inputs for another ReactionData or as direct pre-requisite data for ReactionKernels. The TRIM data example, would pass the values directly to a ReactionKernels to decide velocities via scattering_kernel(). For more details on construction of compatible interpolation pipelines see: Interpolation method details.