From f26433ec09763b80f8b8d2af59a2646fc4488137 Mon Sep 17 00:00:00 2001 From: Adnan-Ali-Ahmad Date: Thu, 17 Feb 2022 21:06:34 +0100 Subject: [PATCH 1/8] added spatial.py & fixed slicing bug in array.py --- src/osyris/__init__.py | 1 + src/osyris/core/array.py | 22 ++++++- src/osyris/spatial/__init__.py | 6 ++ src/osyris/spatial/spatial.py | 104 +++++++++++++++++++++++++++++++++ 4 files changed, 132 insertions(+), 1 deletion(-) create mode 100644 src/osyris/spatial/__init__.py create mode 100644 src/osyris/spatial/spatial.py diff --git a/src/osyris/__init__.py b/src/osyris/__init__.py index 0e9f23fd..57e66ecc 100644 --- a/src/osyris/__init__.py +++ b/src/osyris/__init__.py @@ -6,3 +6,4 @@ from .config import config, units from .plot import histogram1d, histogram2d, plane, scatter, map, plot from .core import Array, Datagroup, Dataset, Plot +from .spatial import spatial diff --git a/src/osyris/core/array.py b/src/osyris/core/array.py index 1cda9e89..2783b820 100644 --- a/src/osyris/core/array.py +++ b/src/osyris/core/array.py @@ -47,7 +47,9 @@ def __getitem__(self, slice_): unit=self._unit, parent=self._parent, name=self._name) - if isinstance(slice_, int) and self.ndim > 1: + if self.shape[0] == 1: + return out + if (isinstance(slice_, int) or isinstance(slice_, np.integer)) and self.ndim > 1: return out.reshape(1, len(out)) return out @@ -253,6 +255,24 @@ def __mul__(self, other): unit=1.0 * result.units) return self.__class__(values=self._array * other, unit=self._unit) + + def __matmul__(self, other): + if isinstance(other, self.__class__): + scale_l = self._unit.to_base_units() + scale_r = other._unit.to_base_units() + result = scale_l * scale_r + lhs = self._array + rhs = other._array * result.magnitude + lhs, rhs = self._broadcast(lhs, rhs) + return self.__class__(values=np.transpose(lhs @ rhs.T), unit=1.0 * result.units) + if isinstance(other, Quantity): + scale_l = self._unit.to_base_units() + scale_r = other.to_base_units() + result = scale_l * scale_r + return self.__class__(values=self._array * result.magnitude, + unit=1.0 * result.units) + return self.__class__(values=self._array @ other, unit=self._unit) + def __imul__(self, other): if isinstance(other, self.__class__): scale_l = self._unit.to_base_units() diff --git a/src/osyris/spatial/__init__.py b/src/osyris/spatial/__init__.py new file mode 100644 index 00000000..7a6daeef --- /dev/null +++ b/src/osyris/spatial/__init__.py @@ -0,0 +1,6 @@ +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2022 Osyris contributors (https://github.com/nvaytet/osyris) + +# flake8: noqa + +from .spatial import cartesian_to_spherical, spherical_to_cartesian, rotation_matrix, rotation_matrix_axis_angle diff --git a/src/osyris/spatial/spatial.py b/src/osyris/spatial/spatial.py new file mode 100644 index 00000000..fc1a6c9a --- /dev/null +++ b/src/osyris/spatial/spatial.py @@ -0,0 +1,104 @@ +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2022 Osyris contributors (https://github.com/nvaytet/osyris) + +import numpy as np +from ..core import Array +from ..config import units + + +def rotation_matrix(alpha, beta, gamma): + """ + Returns the 3D rotation matrix of angles 'alpha' (yaw) around x axis + 'beta' (pitch) around y axis + 'gamma' (roll) around z axis + """ + Rx = np.array([[1, 0, 0], [0, np.cos(alpha), -np.sin(alpha)], + [0, np.sin(alpha), np.cos(alpha)]]) + Ry = np.array([[np.cos(beta), 0, np.sin(beta)], [0, 1, 0], + [-np.sin(beta), 0, np.cos(beta)]]) + Rz = np.array([[np.cos(gamma), -np.sin(gamma), 0], + [np.sin(gamma), np.cos(gamma), 0], [0, 0, 1]]) + + return Rz @ Ry @ Rx + + +def rotation_matrix_axis_angle(vec, angle): + """ + Returns 3D rotation matrix of angle 'angle' around rotation vector 'vec'. + """ + if isinstance(vec, list): + vec = np.array(vec) + R = np.cos(angle) * np.identity(3) + (np.sin(angle)) * np.cross( + vec, + np.identity(vec.shape[0]) * -1) + (1 - np.cos(angle)) * (np.outer(vec, vec)) + return R + + +def spherical_to_cartesian(r, theta, phi): + """ + Converts spherical components radius, colatitude and azimuth to x,y,z cartesian coordinates + Returns an osyris array of shape (len(r), 3) containing a stacked xyz + """ + x, y, z = (r * np.cos(phi) * np.sin(theta), r * np.sin(phi) * np.sin(theta), + r * np.cos(theta)) + xyz = np.vstack([x, y, z]) + return Array(values=np.transpose(xyz).values, unit=r.unit, name="position") + + +def cartesian_to_spherical(position, origin=[0, 0, 0]): + """ + Converts cartesian components in 'position' to spherical components + Returns a list of osyris Arrays containing radius, colatitude and azimuth + """ + if isinstance(origin, Array): + centered_pos = position - origin + else: + centered_pos = position - Array(values=origin, unit=centered_pos.unit) + + X = centered_pos[:, 0].values + Y = centered_pos[:, 1].values + Z = centered_pos[:, 2].values + + radius = np.linalg.norm([X, Y, Z], axis=0) + colatitude = np.arctan2(np.linalg.norm([X, Y], axis=0), Z) + + azimuth = np.arctan2(Y, X) + np.pi + + radius = Array(values=radius, unit=position.unit, name='radius') + colatitude = Array(values=colatitude, unit=units('rad'), name='colatitude') + azimuth = Array(values=azimuth, unit=units('rad'), name='azimuth') + + return radius, colatitude, azimuth + + +def cartesian_to_cylindrical(position, origin=[0, 0, 0]): + """ + Converts cartesian components in 'position' to cylindrical components + Returns a list of osyris Arrays containing radius, azimuth and elevation + """ + if isinstance(origin, Array): + centered_pos = position - origin + else: + centered_pos = position - Array(values=origin, unit=centered_pos.unit) + + elevation = centered_pos[:, 2].values + + radius = np.linalg.norm([centered_pos[:, 0].values, centered_pos[:, 1].values], + axis=0) + azimuth = np.arctan2(centered_pos[:, 1].values, centered_pos[:, 0].values) + + radius = Array(values=radius, unit=centered_pos.unit, name='radius') + azimuth = Array(values=azimuth, unit=units('rad'), name='azimuth') + elevation = Array(values=elevation, unit=centered_pos.unit, name='elevation') + + return radius, azimuth, elevation + + +def cylindrical_to_cartesian(r, azimuth, elevation): + """ + Converts cylindrical components radius, azimuth and elevation to x,y,z cartesian coordinates + Returns an osyris array of shape (len(r), 3) containing a stacked xyz + """ + x, y, z = (radius * np.cos(azimuth), radius * np.sin(azimuth), elevation) + xyz = np.vstack([x, y, z]) + return Array(values=np.transpose(xyz).values, unit=r.unit, name="position") From 8790f8e72908ab5ee4cd4dda633dab39cac6f3bb Mon Sep 17 00:00:00 2001 From: Adnan-Ali-Ahmad Date: Thu, 17 Feb 2022 21:29:58 +0100 Subject: [PATCH 2/8] forgot yapf formatting --- src/osyris/core/array.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/osyris/core/array.py b/src/osyris/core/array.py index 2783b820..23de304f 100644 --- a/src/osyris/core/array.py +++ b/src/osyris/core/array.py @@ -48,8 +48,9 @@ def __getitem__(self, slice_): parent=self._parent, name=self._name) if self.shape[0] == 1: - return out - if (isinstance(slice_, int) or isinstance(slice_, np.integer)) and self.ndim > 1: + return out + if (isinstance(slice_, int) + or isinstance(slice_, np.integer)) and self.ndim > 1: return out.reshape(1, len(out)) return out @@ -255,7 +256,6 @@ def __mul__(self, other): unit=1.0 * result.units) return self.__class__(values=self._array * other, unit=self._unit) - def __matmul__(self, other): if isinstance(other, self.__class__): scale_l = self._unit.to_base_units() @@ -264,7 +264,8 @@ def __matmul__(self, other): lhs = self._array rhs = other._array * result.magnitude lhs, rhs = self._broadcast(lhs, rhs) - return self.__class__(values=np.transpose(lhs @ rhs.T), unit=1.0 * result.units) + return self.__class__(values=np.transpose(lhs @ rhs.T), + unit=1.0 * result.units) if isinstance(other, Quantity): scale_l = self._unit.to_base_units() scale_r = other.to_base_units() From 6b2824ea25fd62f47ed9344a0bb5e3f30aace4ce Mon Sep 17 00:00:00 2001 From: Adnan-Ali-Ahmad Date: Thu, 17 Feb 2022 21:33:26 +0100 Subject: [PATCH 3/8] flake8 --- src/osyris/spatial/spatial.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/osyris/spatial/spatial.py b/src/osyris/spatial/spatial.py index fc1a6c9a..bc97b95d 100644 --- a/src/osyris/spatial/spatial.py +++ b/src/osyris/spatial/spatial.py @@ -36,7 +36,7 @@ def rotation_matrix_axis_angle(vec, angle): def spherical_to_cartesian(r, theta, phi): """ - Converts spherical components radius, colatitude and azimuth to x,y,z cartesian coordinates + Converts spherical components radius, colatitude and azimuth to x,y,z Returns an osyris array of shape (len(r), 3) containing a stacked xyz """ x, y, z = (r * np.cos(phi) * np.sin(theta), r * np.sin(phi) * np.sin(theta), @@ -94,9 +94,9 @@ def cartesian_to_cylindrical(position, origin=[0, 0, 0]): return radius, azimuth, elevation -def cylindrical_to_cartesian(r, azimuth, elevation): +def cylindrical_to_cartesian(radius, azimuth, elevation): """ - Converts cylindrical components radius, azimuth and elevation to x,y,z cartesian coordinates + Converts cylindrical components radius, azimuth and elevation to x,y,z Returns an osyris array of shape (len(r), 3) containing a stacked xyz """ x, y, z = (radius * np.cos(azimuth), radius * np.sin(azimuth), elevation) From 7888bc083389f58fe7370783c9e15298e4711443 Mon Sep 17 00:00:00 2001 From: Adnan-Ali-Ahmad Date: Thu, 17 Feb 2022 21:36:35 +0100 Subject: [PATCH 4/8] last typo fixed --- src/osyris/spatial/spatial.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/osyris/spatial/spatial.py b/src/osyris/spatial/spatial.py index bc97b95d..76267e4b 100644 --- a/src/osyris/spatial/spatial.py +++ b/src/osyris/spatial/spatial.py @@ -101,4 +101,4 @@ def cylindrical_to_cartesian(radius, azimuth, elevation): """ x, y, z = (radius * np.cos(azimuth), radius * np.sin(azimuth), elevation) xyz = np.vstack([x, y, z]) - return Array(values=np.transpose(xyz).values, unit=r.unit, name="position") + return Array(values=np.transpose(xyz).values, unit=radius.unit, name="position") From 41155c705cda7934f93859a88a229d995dfb8f15 Mon Sep 17 00:00:00 2001 From: Adnan-Ali-Ahmad Date: Thu, 17 Feb 2022 22:35:27 +0100 Subject: [PATCH 5/8] removed redundant function and normalized guiding vector in rotation_matrix() --- src/osyris/spatial/__init__.py | 2 +- src/osyris/spatial/spatial.py | 21 +++------------------ 2 files changed, 4 insertions(+), 19 deletions(-) diff --git a/src/osyris/spatial/__init__.py b/src/osyris/spatial/__init__.py index 7a6daeef..8db8a97b 100644 --- a/src/osyris/spatial/__init__.py +++ b/src/osyris/spatial/__init__.py @@ -3,4 +3,4 @@ # flake8: noqa -from .spatial import cartesian_to_spherical, spherical_to_cartesian, rotation_matrix, rotation_matrix_axis_angle +from .spatial import cartesian_to_spherical, spherical_to_cartesian, rotation_matrix diff --git a/src/osyris/spatial/spatial.py b/src/osyris/spatial/spatial.py index 76267e4b..e525c2a4 100644 --- a/src/osyris/spatial/spatial.py +++ b/src/osyris/spatial/spatial.py @@ -6,28 +6,13 @@ from ..config import units -def rotation_matrix(alpha, beta, gamma): - """ - Returns the 3D rotation matrix of angles 'alpha' (yaw) around x axis - 'beta' (pitch) around y axis - 'gamma' (roll) around z axis - """ - Rx = np.array([[1, 0, 0], [0, np.cos(alpha), -np.sin(alpha)], - [0, np.sin(alpha), np.cos(alpha)]]) - Ry = np.array([[np.cos(beta), 0, np.sin(beta)], [0, 1, 0], - [-np.sin(beta), 0, np.cos(beta)]]) - Rz = np.array([[np.cos(gamma), -np.sin(gamma), 0], - [np.sin(gamma), np.cos(gamma), 0], [0, 0, 1]]) - - return Rz @ Ry @ Rx - - -def rotation_matrix_axis_angle(vec, angle): +def rotation_matrix(vec, angle): """ Returns 3D rotation matrix of angle 'angle' around rotation vector 'vec'. """ if isinstance(vec, list): vec = np.array(vec) + vec = vec / np.linalg.norm(vec) R = np.cos(angle) * np.identity(3) + (np.sin(angle)) * np.cross( vec, np.identity(vec.shape[0]) * -1) + (1 - np.cos(angle)) * (np.outer(vec, vec)) @@ -97,7 +82,7 @@ def cartesian_to_cylindrical(position, origin=[0, 0, 0]): def cylindrical_to_cartesian(radius, azimuth, elevation): """ Converts cylindrical components radius, azimuth and elevation to x,y,z - Returns an osyris array of shape (len(r), 3) containing a stacked xyz + Returns an osyris array of shape (len(radius), 3) containing a stacked xyz """ x, y, z = (radius * np.cos(azimuth), radius * np.sin(azimuth), elevation) xyz = np.vstack([x, y, z]) From 3bdee6d53fa43c4aa69ae42b44d522e86cb32ce0 Mon Sep 17 00:00:00 2001 From: Adnan-Ali-Ahmad Date: Fri, 18 Feb 2022 22:10:47 +0100 Subject: [PATCH 6/8] slicing up PR --- src/osyris/core/array.py | 18 ------- src/osyris/spatial/__init__.py | 6 --- src/osyris/spatial/spatial.py | 89 ---------------------------------- 3 files changed, 113 deletions(-) delete mode 100644 src/osyris/spatial/__init__.py delete mode 100644 src/osyris/spatial/spatial.py diff --git a/src/osyris/core/array.py b/src/osyris/core/array.py index 23de304f..705b1ab5 100644 --- a/src/osyris/core/array.py +++ b/src/osyris/core/array.py @@ -256,24 +256,6 @@ def __mul__(self, other): unit=1.0 * result.units) return self.__class__(values=self._array * other, unit=self._unit) - def __matmul__(self, other): - if isinstance(other, self.__class__): - scale_l = self._unit.to_base_units() - scale_r = other._unit.to_base_units() - result = scale_l * scale_r - lhs = self._array - rhs = other._array * result.magnitude - lhs, rhs = self._broadcast(lhs, rhs) - return self.__class__(values=np.transpose(lhs @ rhs.T), - unit=1.0 * result.units) - if isinstance(other, Quantity): - scale_l = self._unit.to_base_units() - scale_r = other.to_base_units() - result = scale_l * scale_r - return self.__class__(values=self._array * result.magnitude, - unit=1.0 * result.units) - return self.__class__(values=self._array @ other, unit=self._unit) - def __imul__(self, other): if isinstance(other, self.__class__): scale_l = self._unit.to_base_units() diff --git a/src/osyris/spatial/__init__.py b/src/osyris/spatial/__init__.py deleted file mode 100644 index 8db8a97b..00000000 --- a/src/osyris/spatial/__init__.py +++ /dev/null @@ -1,6 +0,0 @@ -# SPDX-License-Identifier: BSD-3-Clause -# Copyright (c) 2022 Osyris contributors (https://github.com/nvaytet/osyris) - -# flake8: noqa - -from .spatial import cartesian_to_spherical, spherical_to_cartesian, rotation_matrix diff --git a/src/osyris/spatial/spatial.py b/src/osyris/spatial/spatial.py deleted file mode 100644 index e525c2a4..00000000 --- a/src/osyris/spatial/spatial.py +++ /dev/null @@ -1,89 +0,0 @@ -# SPDX-License-Identifier: BSD-3-Clause -# Copyright (c) 2022 Osyris contributors (https://github.com/nvaytet/osyris) - -import numpy as np -from ..core import Array -from ..config import units - - -def rotation_matrix(vec, angle): - """ - Returns 3D rotation matrix of angle 'angle' around rotation vector 'vec'. - """ - if isinstance(vec, list): - vec = np.array(vec) - vec = vec / np.linalg.norm(vec) - R = np.cos(angle) * np.identity(3) + (np.sin(angle)) * np.cross( - vec, - np.identity(vec.shape[0]) * -1) + (1 - np.cos(angle)) * (np.outer(vec, vec)) - return R - - -def spherical_to_cartesian(r, theta, phi): - """ - Converts spherical components radius, colatitude and azimuth to x,y,z - Returns an osyris array of shape (len(r), 3) containing a stacked xyz - """ - x, y, z = (r * np.cos(phi) * np.sin(theta), r * np.sin(phi) * np.sin(theta), - r * np.cos(theta)) - xyz = np.vstack([x, y, z]) - return Array(values=np.transpose(xyz).values, unit=r.unit, name="position") - - -def cartesian_to_spherical(position, origin=[0, 0, 0]): - """ - Converts cartesian components in 'position' to spherical components - Returns a list of osyris Arrays containing radius, colatitude and azimuth - """ - if isinstance(origin, Array): - centered_pos = position - origin - else: - centered_pos = position - Array(values=origin, unit=centered_pos.unit) - - X = centered_pos[:, 0].values - Y = centered_pos[:, 1].values - Z = centered_pos[:, 2].values - - radius = np.linalg.norm([X, Y, Z], axis=0) - colatitude = np.arctan2(np.linalg.norm([X, Y], axis=0), Z) - - azimuth = np.arctan2(Y, X) + np.pi - - radius = Array(values=radius, unit=position.unit, name='radius') - colatitude = Array(values=colatitude, unit=units('rad'), name='colatitude') - azimuth = Array(values=azimuth, unit=units('rad'), name='azimuth') - - return radius, colatitude, azimuth - - -def cartesian_to_cylindrical(position, origin=[0, 0, 0]): - """ - Converts cartesian components in 'position' to cylindrical components - Returns a list of osyris Arrays containing radius, azimuth and elevation - """ - if isinstance(origin, Array): - centered_pos = position - origin - else: - centered_pos = position - Array(values=origin, unit=centered_pos.unit) - - elevation = centered_pos[:, 2].values - - radius = np.linalg.norm([centered_pos[:, 0].values, centered_pos[:, 1].values], - axis=0) - azimuth = np.arctan2(centered_pos[:, 1].values, centered_pos[:, 0].values) - - radius = Array(values=radius, unit=centered_pos.unit, name='radius') - azimuth = Array(values=azimuth, unit=units('rad'), name='azimuth') - elevation = Array(values=elevation, unit=centered_pos.unit, name='elevation') - - return radius, azimuth, elevation - - -def cylindrical_to_cartesian(radius, azimuth, elevation): - """ - Converts cylindrical components radius, azimuth and elevation to x,y,z - Returns an osyris array of shape (len(radius), 3) containing a stacked xyz - """ - x, y, z = (radius * np.cos(azimuth), radius * np.sin(azimuth), elevation) - xyz = np.vstack([x, y, z]) - return Array(values=np.transpose(xyz).values, unit=radius.unit, name="position") From 90c6381b7da179dcc89bd9318fa3109f87ec33c5 Mon Sep 17 00:00:00 2001 From: Adnan-Ali-Ahmad Date: Fri, 18 Feb 2022 22:11:23 +0100 Subject: [PATCH 7/8] forgot to remove spatial import --- src/osyris/__init__.py | 1 - 1 file changed, 1 deletion(-) diff --git a/src/osyris/__init__.py b/src/osyris/__init__.py index 57e66ecc..0e9f23fd 100644 --- a/src/osyris/__init__.py +++ b/src/osyris/__init__.py @@ -6,4 +6,3 @@ from .config import config, units from .plot import histogram1d, histogram2d, plane, scatter, map, plot from .core import Array, Datagroup, Dataset, Plot -from .spatial import spatial From 2e938a12847d3ab2a0f8792d498c46d4476dd16d Mon Sep 17 00:00:00 2001 From: Adnan-Ali-Ahmad <60321972+Adnan-Ali-Ahmad@users.noreply.github.com> Date: Sun, 20 Feb 2022 16:07:19 +0100 Subject: [PATCH 8/8] Update src/osyris/core/array.py Co-authored-by: Neil Vaytet <39047984+nvaytet@users.noreply.github.com> --- src/osyris/core/array.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/osyris/core/array.py b/src/osyris/core/array.py index 705b1ab5..1ea20cef 100644 --- a/src/osyris/core/array.py +++ b/src/osyris/core/array.py @@ -49,8 +49,7 @@ def __getitem__(self, slice_): name=self._name) if self.shape[0] == 1: return out - if (isinstance(slice_, int) - or isinstance(slice_, np.integer)) and self.ndim > 1: + if isinstance(slice_, (int, np.integer)) and self.ndim > 1: return out.reshape(1, len(out)) return out