[Numpy-discussion] Allow matrix multiplication in structured arrays

[abel . gutierrez]

I'm using structured arrays to store atoms data produced by LAMMPS (I'm using a 
structured array that follows its format). I need to rotate the positions:

```
import numpy as np

transform = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]], dtype=np.float64)
dtype = np.dtype([("x", np.float64), ("y", np.float64), ("z", np.float64)])  # 
real case with more fields, integers, bools, strings

atoms = np.array(
[
(0.0, 0.0, 0.0),
(1.0, 0.0, 0.0),
(0.0, 1.0, 0.0),
(1.0, 1.0, 1.0),
],
dtype=dtype,
)

atoms[["x", "y", "z"]] = atoms[["x", "y", "z"]] @ transform.T
```

But this produces:

```
Traceback (most recent call last):
  File "c:\Users\acgc99\Desktop\rotation.py", line 16, in 
atoms[["x", "y", "z"]] = atoms[["x", "y", "z"]] @ transform.T
 ~~~^
numpy._core._exceptions._UFuncNoLoopError: ufunc 'matmul' did not contain a 
loop with signature matching types (dtype([('x', ' None
```

I can convert to unstructured arrays, but I guess that doing that change 
multiple times is not efficient when working with tens of millions of atoms.
___
NumPy-Discussion mailing list -- numpy-discussion@python.org
To unsubscribe send an email to numpy-discussion-le...@python.org
https://mail.python.org/mailman3//lists/numpy-discussion.python.org
Member address: arch...@mail-archive.com

[Numpy-discussion] Re: Allow matrix multiplication in structured arrays

[abel . gutierrez]

Hi,

Thanks for replying. I was using:
```
from numpy.lib import recfunctions as rfn
xyz = rfn.structured_to_unstructured(atoms[["x", "y", "z"]])
xyz = xyz @ transform.T
atoms[["x", "y", "z"]] = rfn.unstructured_to_structured(
xyz, dtype=atoms[["x", "y", "z"]].dtype
)
```
But I think this is creating a copy not a view.

Anyway, I did the following:
```
fields = atoms.dtype.names
idxx, idxz = fields.index("x"), fields.index("z")
xyz_field = ("xyz", (np.float64, (3,)))
# create a new dtype with previous fields, xyz and next fields
new_fields = (
[(n, atoms.dtype[n]) for n in fields[:idxx]]
+ [xyz_field]
+ [(n, atoms.dtype[n]) for n in fields[idxz + 1 :]]
)  # usually x, y, z fields are one next to the other
new_dtype = np.dtype(new_fields)

xyz = atoms.view(new_dtype)["xyz"]
xyz[:] = xyz @ transform.T
```
And works fine. I guess there is no direct way computing the rotation without a 
view.

Thanks Marten,

Best,
Abel
___
NumPy-Discussion mailing list -- numpy-discussion@python.org
To unsubscribe send an email to numpy-discussion-le...@python.org
https://mail.python.org/mailman3//lists/numpy-discussion.python.org
Member address: arch...@mail-archive.com