`map_and_read_data` does not flatten coordinate array before passing it to C++

[r-walloner]

If I follow the example given in help(precice.Participant.map_and_read_data), map_and_read_data throws the following error.

Traceback (most recent call last):
  File "/home/robin/uni/ba-repo/investigate_relativeReadTime/B/B.py", line 24, in <module>
    data = participant.map_and_read_data("A-Mesh", "A-Data", coordinates, 0)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "cyprecice/cyprecice.pyx", line 1000, in cyprecice.Participant.map_and_read_data
  File "<stringsource>", line 48, in vector.from_py.__pyx_convert_vector_from_py_double
TypeError: only length-1 arrays can be converted to Python scalars

The problematic line seems to be this one.
It seems to assume the coordinates are already flattened: [x_0, y_0, x_1, y_1, ...], but they are not: [(x_0, y_0), (x_1, y_1), ...].

[r-walloner]

Minimal Example to reproduce

A.py

import precice

participant = precice.Participant("A", "../precice-config.xml", 0, 1)

mesh = [[0, 0]]
vertex_ids = participant.set_mesh_vertices("A-Mesh", mesh)
participant.initialize()

t = 0
dt = 1

while participant.is_coupling_ongoing():
    t = t + dt
    print("Timestep: ", t)

    data = t
    print("Writing data: ", data)
    participant.write_data("A-Mesh", "A-Data", vertex_ids, [data])

    participant.advance(dt)

participant.finalize()

B.py

import precice

participant = precice.Participant("B", "../precice-config.xml", 0, 1)

participant.set_mesh_access_region("A-Mesh", [-1, 1, -1, 1])
participant.initialize()

t = 0
dt = 1

while participant.is_coupling_ongoing():
    t = t + dt
    print("Timestep: ", t)

    relative_read_time = dt * 0
    print("Relative read time: ", relative_read_time)
    print("Absolute read time: ", t - dt + relative_read_time)

    coordinates = [(0.0, 0.0)]
    data = participant.map_and_read_data("A-Mesh", "A-Data", coordinates, relative_read_time)
    print("Reading data: ", data)

    participant.advance(dt)

participant.finalize()

precice-config.xml

<precice-configuration experimental="true">

  <data:scalar name="A-Data" />

  <mesh name="A-Mesh" dimensions="2">
    <use-data name="A-Data" />
  </mesh>

  <participant name="A">
    <provide-mesh name="A-Mesh" />
    <write-data name="A-Data" mesh="A-Mesh" />
  </participant>

  <participant name="B">
    <receive-mesh name="A-Mesh" from="A" api-access="true" />
    <read-data name="A-Data" mesh="A-Mesh" />
    <mapping:nearest-neighbor
      direction="read"
      from="A-Mesh"
      constraint="consistent">
    </mapping:nearest-neighbor>
  </participant>

  <m2n:sockets acceptor="A" connector="B" exchange-directory=".."/>

  <coupling-scheme:serial-explicit>
    <participants first="A" second="B" />
    <max-time value="4" />
    <time-window-size value="1" />
    <exchange data="A-Data" mesh="A-Mesh" from="A" to="B" />
  </coupling-scheme:serial-explicit>

</precice-configuration>

[BenjaminRodenberg]

Thanks for reporting this bug and also submitting a fix! Did you also face problems with the API function write_and_map_data? We just implemented these two functions for the latest preCICE release and I would not be surprised if the same bug as for map_and_read_data appears for write_and_map_data as well.

[r-walloner]

write_and_map_data seems to work fine.

But I find it slightly odd, that it takes the coordinates in a different format than map_and_read_data.
In write_and_map_data coordinates are a 2D numpy array np.array([[c1_x, c1_y], [c2_x, c2_y]]).
In map_and_read_data coordinates are a list of touples [(c1_x, c1_y), (c2_x, c2_y)].
Is this intended?

[davidscn]

I think this functionality was added by @MakisH, so maybe he knows best: was there any reason for the differences in the function interfaces?

[MakisH]

I think this functionality was added by @MakisH, so maybe he knows best: was there any reason for the differences in the function interfaces?

When implementing this with @BenjaminRodenberg in #231, we started by copying the existing write_data and read_data. This issue should also be present in those functions, correct?

I remember a similar confusion while implementing, with the conclusion that "it is like this at the moment, and it would be a bigger task to change that". But I agree it is confusing and we should of course aim for consistency.

[davidscn]

But I find it slightly odd, that it takes the coordinates in a different format

the coordinates are only part of the newly added functions and not of read/write data. So where did you copy? Maybe a more directed question: how easy would it be to unify the interface and which of the two versions is more convenient? In my (non-python) opinion, the coordinates type should be the same type we pass to set_mesh_vertices. What does this function use?

[MakisH]

Looking at the documentation of the relevant functions, there is indeed an inconsistency on the type of coordinates (named positions in set_mesh_vertices), which is only demonstrated in the examples:

    def set_mesh_vertices (self, mesh_name, positions):
        """
        ...
        positions : array_like
            The coordinates of the vertices in a numpy array [N x D] where
            N = number of vertices and D = dimensions of geometry.

        Returns
        -------
        vertex_ids : numpy.ndarray
            IDs of the created vertices.

        Examples
        --------

        Set mesh vertices for a 3D problem with 5 mesh vertices.

        >>> positions = np.array([[1, 1, 1], [2, 2, 2], [3, 3, 3], [4, 4, 4], [5, 5, 5]])
        >>> positions.shape
        (5, 3)
        >>> mesh_name = "MeshOne"
        >>> vertex_ids = participant.set_mesh_vertices(mesh_name, positions)
        >>> vertex_ids.shape
        (5,)
        """

   def write_and_map_data (self, mesh_name, data_name, coordinates, values):
        """
        ...
        coordinates : array_like
            The coordinates of the vertices in a numpy array [N x D] where
            N = number of vertices and D = dimensions of geometry.

        Examples
        --------

        >>> coordinates = np.array([[c1_x, c1_y], [c2_x, c2_y], [c3_x, c3_y], [c4_x, c4_y], [c5_x, c5_y]])
        """

    def map_and_read_data (self, mesh_name, data_name, coordinates, relative_read_time):
        """

        coordinates : array_like
            Coordinates of the vertices.

        >>> coordinates = [(1.0, 1.0), (2.0, 2.0)]
        """

the coordinates are only part of the newly added functions and not of read/write data. So where did you copy? Maybe a more directed question: how easy would it be to unify the interface and which of the two versions is more convenient?

True. We did copy from the read_data and write_data, trying to stay consistent with the set_mesh_vertices. Apparently, at least in map_and_read_data, we did not, but there is a (not that good) reason.

Comparing the two implementations:

• In the write_and_map_data, coordinates are given as an numpy.ndarray, flattened, and given to the C++ API, which expects a span<const double>. Same for the values, and both can be given in either format to the Python API and automatically converted.
• In the map_and_read_data, coordinates are expected as a flat array, which are passed directly to the C++ API (which again expects a span<const double> and then converted to a numpy.ndarray before being returned. A check here would be good, but the fix in Flatten coordinate array in map_and_read_data #237 is already one step ahead.

python-bindings/cyprecice/cyprecice.pyx

Lines 906 to 1012 in 43a82f4

	def write_and_map_data (self, mesh_name, data_name, coordinates, values):
	"""
	This function writes values at temporary locations to data of a mesh.
	As opposed to the writeData function using VertexIDs, this function allows to write data via coordinates,
	which don't have to be specified during the initialization. This is particularly useful for meshes, which
	vary over time. Note that using this function comes at a performance cost, since the specified mapping
	needs to be computed locally for the given locations, whereas the other variant (writeData) can typically
	exploit the static interface mesh and pre-compute data structures more efficiently.
	Values are passed identically to write_data.
	Parameters
	----------
	mesh_name : str
	name of the mesh to write to.
	data_name : str
	Data name to write to.
	coordinates : array_like
	The coordinates of the vertices in a numpy array [N x D] where
	N = number of vertices and D = dimensions of geometry.
	values : array_like
	Values of data
	Examples
	--------
	Write scalar data for a 2D problem with 5 vertices:
	>>> mesh_name = "MeshOne"
	>>> data_name = "DataOne"
	>>> coordinates = np.array([[c1_x, c1_y], [c2_x, c2_y], [c3_x, c3_y], [c4_x, c4_y], [c5_x, c5_y]])
	>>> values = np.array([v1, v2, v3, v4, v5])
	>>> participant.write_and_map_data(mesh_name, data_name, coordinates, values)
	"""
	check_array_like(coordinates, "coordinates", "write_and_map_data")
	check_array_like(values, "values", "write_and_map_data")
	if not isinstance(coordinates, np.ndarray):
	coordinates = np.asarray(coordinates)
	if not isinstance(values, np.ndarray):
	values = np.asarray(values)
	cdef vector[double] cpp_coordinates = coordinates.flatten()
	cdef vector[double] cpp_values = values.flatten()
	self.thisptr.writeAndMapData (convert(mesh_name), convert(data_name), cpp_coordinates, cpp_values)
	def map_and_read_data (self, mesh_name, data_name, coordinates, relative_read_time):
	"""
	This function reads values at temporary locations from data of a mesh.
	As opposed to the readData function using VertexIDs, this function allows reading data via coordinates,
	which don't have to be specified during the initialization. This is particularly useful for meshes, which
	vary over time. Note that using this function comes at a performance cost, since the specified mapping
	needs to be computed locally for the given locations, whereas the other variant (readData) can typically
	exploit the static interface mesh and pre-compute data structures more efficient.
	Values are read identically to read_data.
	Parameters
	----------
	mesh_name : str
	Name of the mesh to write to.
	data_name : str
	Name of the data to read from.
	coordinates : array_like
	Coordinates of the vertices.
	relative_read_time : double
	Point in time where data is read relative to the beginning of the current time step
	Returns
	-------
	values : numpy.ndarray
	Contains the read data.
	Examples
	--------
	Read scalar data for a 2D problem with 2 vertices:
	>>> mesh_name = "MeshOne"
	>>> data_name = "DataOne"
	>>> coordinates = [(1.0, 1.0), (2.0, 2.0)]
	>>> dt = 1.0
	>>> values = map_and_read_data(mesh_name, data_name, coordinates, dt)
	>>> values.shape
	>>> (2, )
	"""
	check_array_like(coordinates, "coordinates", "map_and_read_data")
	if not isinstance(coordinates, np.ndarray):
	coordinates = np.asarray(coordinates)
	size = coordinates.shape[0]
	dimensions = self.get_data_dimensions(mesh_name, data_name)
	cdef vector[double] cpp_coordinates = coordinates.flatten()
	cdef vector[double] cpp_values = [-1 for _ in range(size * dimensions)]
	self.thisptr.mapAndReadData (convert(mesh_name), convert(data_name), cpp_coordinates, relative_read_time, cpp_values)
	cdef np.ndarray[double, ndim=1] np_values = np.array(cpp_values, dtype=np.double)
	if len(coordinates) == 0:
	return np_values.reshape((size))
	elif self.get_data_dimensions(mesh_name, data_name) == 1:
	return np_values.reshape((size))
	else:
	return np_values.reshape((size, dimensions))

This implementation felt more natural from the C++ API requirements, but of course does not make much sense for the user.

Without being familiar enough with the Python bindings, I also think we should change both the write_and_map_data and map_and_read_data implementation to expect a numpy.ndarray.

@IshaanDesai would be the best to provide a recommendation on this.

[IshaanDesai]

All the functions in question have the following lines in them:

if not isinstance(coordinates, np.ndarray):
     coordinates = np.asarray(coordinates)

which allows the flexibility to pass a list of tuples instead of a NumPy array. The examples in the docstrings of each of the functions purposely show different variants, to highlight the flexibility. So it is in fact intended to have different types of inputs. Internally a NumPy ndarray is created and used.

If my explanation is satisfactory, can I close this issue?

[davidscn]

The examples in the docstrings of each of the functions purposely show different variants, to highlight the flexibility. So it is in fact intended to have different types of inputs. Internally a NumPy ndarray is created and used.

This means that the comment of @r-walloner

But I find it slightly odd, that it takes the coordinates in a different format

is misleading, because he could pass the argument in both functions as both variants, right? Would it make sense to mention this somewhere in the docs(trings)?

[IshaanDesai]

because he could pass the argument in both functions as both variants, right? Would it make sense to mention this somewhere in the docs(trings)?

That is correct. I will open a PR modifying the docstrings to clearly state the flexibility.