Currently, set_concentrations in our python interface does one pybind11 __setitem__ per species per grid cell in a Python loop. The same loop appears in set_user_defined_rate_parameters, get_concentrations, get_user_defined_rate_parameters, set_conditions, and get_conditions.
This is because when we expose a vector<double> with pybind11's bind_vector, we do not pass the buffer protocol. Therefore, each memory access is one Python operation per item, rather than vectorized.
This also loses writes. Without the buffer protocol, np.asarray(state.concentrations) copies through the sequence protocol, so np.asarray(state.concentrations)[0] = 99.0 is silently discarded.
We should fix this.
Acceptance criteria
- Update our state bindings so that each access in memory is not a
__setitem__.
- Add a test that writes through
np.asarray(...) and confirms the C++ object sees the value.
Ideas
-
consider defining our vector of double as py::bind_vector<std::vector<double>>(m, "VectorDouble", py::buffer_protocol());
-
Instead of exposing a vector of doubles, consider exposing a py::array_t similar to how miem does, which does support vectorization, and then do this everywhere we would return a vector. The shape comes with the array, so Python does no index math:
view = state.concentrations # zero-copy, shape (n_groups, n_species, L)
n_groups, _, L = view.shape
padded = np.zeros(n_groups * L)
padded[:n_cells] = values
view[:, i_species, :] = padded.reshape(n_groups, L)
-
note that GetConcentrationsStrides() is not enough to shape that array. For VectorMatrix, RowStride() is 1 and ColumnStride() is L, which only describe the layout inside one group. The n_species * L jump between groups is never reported, so a 2D (n_cells, n_species) view is wrong once n_cells > L. The general shape is 3D: (n_groups, n_species, L), with padding in the last group.
-
we could also consider exposing a buffer protocol directly on our micm matrices, but that breaks the separation of concerns and sort of negates the usefulness of the musica API
Currently, set_concentrations in our python interface does one pybind11
__setitem__per species per grid cell in a Python loop. The same loop appears in set_user_defined_rate_parameters, get_concentrations, get_user_defined_rate_parameters, set_conditions, and get_conditions.This is because when we expose a
vector<double>with pybind11'sbind_vector, we do not pass thebuffer protocol. Therefore, each memory access is one Python operation per item, rather than vectorized.This also loses writes. Without the buffer protocol,
np.asarray(state.concentrations)copies through the sequence protocol, sonp.asarray(state.concentrations)[0] = 99.0is silently discarded.We should fix this.
Acceptance criteria
__setitem__.np.asarray(...)and confirms the C++ object sees the value.Ideas
consider defining our vector of double as
py::bind_vector<std::vector<double>>(m, "VectorDouble", py::buffer_protocol());Instead of exposing a vector of doubles, consider exposing a
py::array_tsimilar to how miem does, which does support vectorization, and then do this everywhere we would return a vector. The shape comes with the array, so Python does no index math:note that
GetConcentrationsStrides()is not enough to shape that array. ForVectorMatrix,RowStride()is1andColumnStride()isL, which only describe the layout inside one group. Then_species * Ljump between groups is never reported, so a 2D(n_cells, n_species)view is wrong oncen_cells > L. The general shape is 3D:(n_groups, n_species, L), with padding in the last group.we could also consider exposing a buffer protocol directly on our micm matrices, but that breaks the separation of concerns and sort of negates the usefulness of the musica API