Source code for galaximview.basefuncs
import numpy as np
from numpy import linalg as la
[docs]def identity_function(x):
"""Returns the argument.
Parameters
----------
x : any
Argument.
Returns
-------
any
Argument.
"""
return x
[docs]def inverse_function(function, x):
"""Returns the value of the inverse function of 'function' at x for a few functions.
Parameters
----------
function : ufunc
Function whose inverse function will be computed at x.
x : float (depends on function)
Argument.
Returns
-------
float (depends on function)
Argument.
"""
if function == identity_function:
return identity_function(x)
elif function == np.log10:
return 10 ** x
elif function == np.log:
return np.exp(x)
elif function == np.sqrt:
return x ** 2
else:
raise ValueError("Function not inversable yet.")
[docs]def calc_com(pos, mass):
"""Computes centre of mass of particles.
Parameters
----------
pos : ndarray
(N,3) positions.
mass : ndarray
(N,) masses.
Returns
-------
ndarray
(3,) Centre of mass.
"""
CheckArrays.check_poslike_array(pos)
CheckArrays.check_masslike_array(mass)
if not np.any(mass):
raise ValueError("Mass array should not be zero everywhere.")
if any(mass < 0):
raise ValueError("Mass array should be positive everywhere.")
if np.shape(pos)[0] != np.shape(mass)[0]:
raise ValueError("Should be given pos array of shape (3,N) and mass array of shape (N,).")
com = np.zeros((3), dtype='float64')
for k in range(0, 3):
com[k] = np.sum(pos[:, k] * mass, dtype="float64") / np.sum(mass, dtype="float64")
return com
[docs]def get_basis_los(azim, elev):
"""Gets the basis on which to project to obtain 2D projection when the line-of-sight has an azimuth (angle between
unit vector of x-axis and line-of-sight) azim and an elevation (angle between (xy) plane and line-of-sight) elev.
When navigating "bove" or "below" z-axis in 3D plot, the elevation still increases (resp. diminishes).
Parameters
----------
azim : float
In [-180,180], azimuth of camera in matplotlib 3D plot.
elev : float
In [-180,180], elevation of camera in matplotlib 3D plot.
Returns
-------
ndarray
(3,3) Basis on which to project.
Notes
-----
For example:
- An x vs y projection plot (with ascending x to right and y to top), in original x and y coordinates ('face-on')
corresponds to a viewing-angle with azim = -90 and elev = 90.
- A z vs x plot corresonds to azim = -90 and elev = 0
"""
azim_in_radians = azim * np.pi / 180.
elev_in_radians = elev * np.pi / 180.
if (elev_in_radians > np.pi / 2.):
elev = np.pi - elev_in_radians
phi = azim_in_radians + np.pi
elif (elev_in_radians < -np.pi / 2.):
elev = - np.pi - elev_in_radians
phi = azim_in_radians + np.pi
else:
elev = elev_in_radians
phi = azim_in_radians
theta = np.pi / 2. - elev
elos = np.array([np.sin(theta) * np.cos(phi), np.sin(theta) * np.sin(phi), np.cos(theta)])
elos = BaseTools.normalise_vector_to_unity(elos)
ex, ey, ez = BaseTools.cart_ex_3d, BaseTools.cart_ey_3d, BaseTools.cart_ez_3d
if (np.allclose(elos, ez)):
elos = ez
exlos = BaseTools.rotate_around_vector(ex, elos, phi + np.pi / 2.)
eylos = BaseTools.rotate_around_vector(ey, elos, phi + np.pi / 2.)
elif (np.allclose(elos, -ez)):
elos = -ez
exlos = BaseTools.rotate_around_vector(ex, elos, np.pi / 2. - phi)
eylos = - BaseTools.rotate_around_vector(ey, elos, np.pi / 2. - phi)
else:
exlos = np.cross(ez, elos)
eylos = np.cross(elos, exlos)
exlos = BaseTools.normalise_vector_to_unity(exlos)
eylos = BaseTools.normalise_vector_to_unity(eylos)
if ((elev_in_radians > np.pi / 2.) | (elev_in_radians < -np.pi / 2.)):
if (np.allclose(elos, -ez) is False):
# switch view orientation
exlos = -exlos
eylos = -eylos
basis_los = np.vstack((exlos, eylos, elos))
return basis_los
[docs]class BaseTools:
"""Class grouping basic array manipulation tools."""
def __init__(self):
pass
cart_ex_3d = np.array([1., 0, 0])
cart_ey_3d = np.array([0, 1., 0])
cart_ez_3d = np.array([0, 0, 1.])
cart_ex_2d = np.array([1., 0])
cart_ey_2d = np.array([0, 1.])
[docs] @classmethod
def sph_kernel(cls, r, h):
"""General function for Gadget-2 type SPH kernel. Can check that quad(lambda x: BaseTools.sph_kernel(x,
h)*4*pi*x**2, 0,h) = 1 for chosen h.
Parameters
----------
r : float
Distance.
h : float
Smoothing length.
Returns
-------
float
SPH kernel.
"""
if ((r >= 0) & (r <= 0.5 * h)):
w = 1 - 6 * (r / h) ** 2 + 6 * (r / h) ** 3
elif ((r > 0.5 * h) & (r <= h)):
w = 2 * (1 - r / h) ** 3
else:
w = 0
w *= 8 / (np.pi * h ** 3)
return w
[docs] @classmethod
def average64_1d(cls, arr, weights=None):
"""Computes average of a 1D array weighted by weights, in float64 precision.
Parameters
----------
arr : ndarray
(N,) 1D array.
weights : ndarray, optional
(N,) weights.
Returns
-------
float
(Weighted) average of array.
Raises
------
ValueError
If Array and weights are not of same shape.
"""
if (weights is None):
weights = np.ones_like(arr)
CheckArrays.check_masslike_array(arr)
CheckArrays.check_masslike_array(weights)
if np.shape(arr)[0] != np.shape(weights)[0]:
raise ValueError("Should be given arrays of same shape (N,).")
average = np.sum(arr * weights, dtype='float64') / np.sum(weights, dtype='float64')
return average
[docs] @classmethod
def dispersion64_1d(cls, arr, weights=None):
"""Computes square root of variance of a 1D array weighted by weights, in float64 precision.
Parameters
----------
arr : ndarray
(N,) 1D array.
weights : ndarray, optional
(N,) weights.
Returns
-------
float64
(Weighted) dispersion of array.
Raises
------
ValueError
If Array and weights are not of same shape.
"""
CheckArrays.check_masslike_array(arr)
if (weights is None):
weights = np.ones_like(arr)
CheckArrays.check_masslike_array(weights)
if np.shape(arr)[0] != np.shape(weights)[0]:
raise ValueError("Should be given arrays of same shape (N,).")
disp = np.sum(arr ** 2 * weights, dtype='float64') / np.sum(weights, dtype='float64') - (
np.sum(arr * weights, dtype='float64') / np.sum(weights, dtype='float64')) ** 2
disp = np.sqrt(disp)
return disp
[docs] @classmethod
def subset_of_indexes(cls, indexes, fraction, seed=None):
"""Extracts random fraction of indexes with possibility to set the seed for identical reproduction.
Parameters
----------
indexes : ndarray
(N,) indexes array.
fraction : float
Fraction to select.
seed : int, optional
Seed for random numbers generator.
Returns
-------
ndarray
(Nselected,) randomly selected indexes such that Nselected=int(frac * N) with N the number if lembers of
1D input array.
"""
CheckArrays.check_indslike_array(indexes)
if ((fraction < 0) | (fraction > 1)):
raise ValueError("Fraction should be between 0 and 1.")
if (seed is not None):
np.random.seed(seed)
if (fraction != 1):
npartnew = int(fraction * np.size(indexes))
subsetindexes = np.random.randint(0, np.size(indexes), size=npartnew)
subsetindexes = indexes[subsetindexes]
else:
subsetindexes = indexes
return subsetindexes
[docs] @classmethod
def normalise_vector_to_unity(cls, vect):
"""Normalises array.
Parameters
----------
vect : ndarray
Array to normalise.
Returns
-------
ndarray
Array normalised to unity.
"""
if not np.any(vect):
raise ValueError("Can't normalise null array to unity.")
elif np.any(np.isnan(vect)):
raise ValueError("Can't normalise array containing nans to unity.")
else:
vect = vect.astype('float64')
norm = la.norm(vect)
if (norm == np.inf):
raise ValueError('Norm is infinite.')
elif (norm < np.finfo(float).eps):
print("Norm=", norm)
raise ValueError('Norm is too small.')
else:
vect /= la.norm(vect)
return vect
[docs] @classmethod
def cart_3d_basis(cls):
"""Returns np.vstack((BaseDefs.cart_ex_3d, BaseDefs.cart_ey_3d, BaseDefs.cart_ez_3d))
Returns
-------
ndarray
(3,3) cartesian 3D basis.
"""
return np.vstack((cls.cart_ex_3d, cls.cart_ey_3d, cls.cart_ez_3d))
[docs] @classmethod
def cart_2d_basis(cls):
"""Returns np.vstack((BaseDefs.cart_ex_2d, BaseDefs.cart_ey_2d))
Returns
-------
ndarray
(2,2) cartesian 2D basis.
"""
return np.vstack((cls.cart_ex_2d, cls.cart_ey_2d))
[docs] @classmethod
def invert_basis(cls, basis):
"""Returns the inverses basis (or any inverse of a (N,N) array).
Parameters
----------
basis : ndarray
(N,N) basis to invert.
Returns
-------
ndarray
(N,N) inverse basis.
"""
return la.inv(basis)
[docs] @classmethod
def check_ortho_normed_direct_basis(cls, vects_base):
"""Checks that (N,N) array with N = 2 or 3 (exclusively) contains N vstacked vectors each normed to unity,
orthonormal to one another and directly ordered.
Parameters
----------
vects_base : ndarray
(N_dims, N_dims) vertically stacked basis arrays. N_dims must be 2 or 3.
Returns
-------
bool
Returns True if basis passed the check.
Raises
------
ValueError
If requirements for shapes or values are not satisfied.
"""
if ((len(np.shape(vects_base)) != 2) | (np.shape(vects_base)[0] != np.shape(vects_base)[1])):
raise ValueError("Must be given vstacked basis vectors, vects_basis must be of shape (N_dim, N_dim).")
for k in range(0, np.shape(vects_base)[0]):
normedvect = cls.normalise_vector_to_unity(vects_base[k])
diff = abs(normedvect - vects_base[k, :])
if np.any(np.round(diff, 7) != 0):
raise ValueError("The basis vectors are not normed to unity.")
if np.shape(vects_base)[0] == 2:
dotprod = np.dot(vects_base[0], vects_base[1])
diff = abs(dotprod)
if np.any(np.round(diff, 7) != 0):
raise ValueError("The 2D basis vectors are not orthogonal.")
crossprod = np.cross(vects_base[0, :], vects_base[1, :])
if crossprod <= 0:
raise ValueError("The 2D basis vectors are not directly orthogonal.")
elif np.shape(vects_base)[0] == 3:
eztheor = np.cross(vects_base[0, :], vects_base[1, :])
diff = abs(vects_base[2, :] - eztheor)
if np.any(np.round(diff, 7) != 0):
raise ValueError("The basis is not orthogonal.")
else:
raise ValueError("Must be given vstacked basis vectors, vects_base must be of shape (N_dim, N_dim), "
"with N_dim=2 or N_dim=3.")
return True
[docs] @classmethod
def project_on_orthonormed_direct_basis(cls, vect_to_project, vects_base):
"""Projects pos array of shape (N_elements, N_dims) on vstacked base of shape (N_dim, N_dim).
Parameters
----------
vect_to_project : ndarray
(N_elements,N_dims) vertically stacked (N_dims,) arrays to be projecteded.
vects_base : ndarray
(N_dims, N_dims) vertically stacked basis arrays.
Returns
-------
ndarray
(N_elements, N_dims) projected array.
"""
cls.check_ortho_normed_direct_basis(vects_base)
vects_base = vects_base.astype('float64')
vect_to_project = vect_to_project.astype('float64')
if (np.shape(vect_to_project)[-1] != np.shape(vects_base)[0]):
raise ValueError("Must be given (N_elements,N_dim) array to be projected and a N_dim dimensions base.")
ndim = np.shape(vects_base)[0]
posproj = np.zeros_like(vect_to_project)
for k in range(0, ndim):
if (len(np.shape(posproj)) == 1):
posproj[k] = np.dot(vect_to_project, vects_base[k, :])
else:
posproj[:, k] = np.dot(vect_to_project, vects_base[k, :])
return posproj
[docs] @classmethod
def rotate_around_vector(cls, pos, elos, alpharot):
"""Performs rotation of every vector of shape (3,) of the pos array of shape(N,3) of angle alpharot around
vector elos of shape (3,).
Parameters
----------
pos : ndarray
(N,3) vertically stacked (3,) arrays to be rotated.
elos : ndarray
(3,) vector around whichrotation is performed.
alpharot : float
Angle of rotation.
Returns
-------
ndarray
(N,3) rotated array.
"""
if ((np.shape(elos) != (3,))):
raise ValueError("only works in 3D")
if (np.shape(pos)[-1] != 3):
raise ValueError("must be given positions of shape (N,3) and a 3D vector to rotate around")
elos = cls.normalise_vector_to_unity(elos)
idmatrix = np.identity(3)
qmatrix = np.array([[0, -elos[2], elos[1]], [elos[2], 0, -elos[0]], [-elos[1], elos[0], 0]])
rotmat = idmatrix + np.sin(alpharot) * qmatrix + (1 - np.cos(alpharot)) * np.dot(qmatrix, qmatrix)
return np.transpose(np.dot(rotmat, np.transpose(pos)))
[docs]class CheckArrays:
"""Class grouping methods that check that arrays given to some methods have the right shapes, signs."""
def __init__(self):
pass
[docs] @classmethod
def check_poslike_array(cls, arr, allow_nan=False, allow_inf=False):
"""Checks that array is of shape (N,3), with no inf nor nan values except if allow_inf or
allow_nan are set to True.
Parameters
----------
arr : ndarray
(N,3) input array.
allow_nan : bool, default = `False`
`True` to allow nan values.
allow_inf : bool, default = `False`
`True` to allow inf values.
Returns
-------
bool
Returns True if input array passed the check.
Raises
------
TypeError
If requirements of type and shape are not satisfied.
ValueError
If requirements for values (limits, nan, inf) are not satisfied.
"""
if type(arr) != np.ndarray:
raise TypeError("Should be given a numpy array.")
else:
if (allow_nan is False):
if np.any(np.isnan(arr)):
raise ValueError("Array contains nan values.")
if (allow_inf is False):
if np.any(np.isinf(arr)):
raise ValueError("Array contains inf values.")
if ((len(np.shape(arr)) != 2) | (np.shape(arr)[1] != 3)):
raise TypeError("Should be given positions type array of shape (N,3).")
return True
[docs] @classmethod
def check_cylposlike_array(cls, arr, allow_nan=False, allow_inf=False):
"""Checks that array is of shape (N,3), that the first column is positive (cylindrical radii) and the
second between -pi and pi (phi angle). (Last column is z.) With no inf nor nan values except if allow_inf or
allow_nan are set to True.
Parameters
----------
arr : ndarray
(N,3) input array.
allow_nan : bool, default = `False`
`True` to allow nan values.
allow_inf : bool, default = `False`
`True` to allow inf values.
Returns
-------
bool
Returns True if input array passed the check.
Raises
------
TypeError
If requirements of type and shape are not satisfied.
ValueError
If requirements for values (limits, nan, inf) are not satisfied.
"""
cls.check_poslike_array(arr, allow_nan=allow_nan, allow_inf=allow_inf)
if np.any(arr[:, 0] < 0):
raise ValueError("Radii contain negative values.")
if (np.any(arr[:, 1] < -np.pi)) | (np.any(arr[:, 1] > np.pi)):
raise ValueError("Phi angles contain value < -pi or > +pi.")
return True
[docs] @classmethod
def check_spherposlike_array(cls, arr, allow_nan=False, allow_inf=False):
"""Checks that array is of shape (N,3), that the first column is positive (spherical radii),
the second between 0 and pi (theta angle) and the third betwwen -pi and pi (phi angle). With no inf nor nan
values except if allow_inf or allow_nan are set to True.
Parameters
----------
arr : ndarray
(N,3) input array.
allow_nan : bool, default = `False`
`True` to allow nan values.
allow_inf : bool, default = `False`
`True` to allow inf values.
Returns
-------
bool
Returns True if input array passed the check.
Raises
------
TypeError
If requirements of type and shape are not satisfied.
ValueError
If requirements for values (limits, nan, inf) are not satisfied.
"""
cls.check_poslike_array(arr, allow_nan=allow_nan, allow_inf=allow_inf)
if np.any(arr[:, 0] < 0):
raise ValueError("Radii contain negative values.")
if (np.any(arr[:, 1] < 0)) | (np.any(arr[:, 1] > np.pi)):
raise ValueError("Theta angles contains value < 0 or > +pi.")
if (np.any(arr[:, 2] < -np.pi)) | (np.any(arr[:, 2] > np.pi)):
raise ValueError("Phi angles contains value < -pi or > +pi.")
return True
[docs] @classmethod
def check_masslike_array(cls, arr, allow_nan=False, allow_inf=False):
"""Checks that array is a 1D array with no inf or nan values except if allow_inf or allow_nan are set to True.
Parameters
----------
arr : ndarray
(N,) input 1D array.
allow_nan : bool, default = `False`
`True` to allow nan values.
allow_inf : bool, default = `False`
`True` to allow inf values.
Returns
-------
bool
Returns True if input array passed the check.
Raises
------
TypeError
If requirements of type and shape are not satisfied.
ValueError
If requirements for values (nan, inf) are not satisfied.
"""
if (type(arr) != np.ndarray):
raise TypeError("Should be given a numpy array.")
if (allow_nan is False):
if np.any(np.isnan(arr)):
raise ValueError("Array contains nan values.")
if (allow_inf is False):
if np.any(np.isinf(arr)):
raise ValueError("Array contains inf values.")
if len(np.shape(arr)) != 1:
raise TypeError("Should be given array of shape (N,).")
return True
[docs] @classmethod
def check_indslike_array(cls, arr, allow_nan=False, allow_inf=False, allow_neg=False):
"""Checks that array is an array of natural integers with no inf or nan values except if allow_neg,
allow_inf or allow_nan are set to True.
Parameters
----------
arr : ndarray
(N,) input array.
allow_nan : bool, default = `False`
`True` to allow nan values.
allow_inf : bool, default = `False`
`True` to allow inf values.
allow_neg : bool, default = `False`
`True` to allow negatvie values.
Returns
-------
bool
Returns True if input array passed the check.
Raises
------
TypeError
If requirements of type are not satisfied.
ValueError
If requirements for values (nan, inf, negative) are not satisfied.
"""
if type(arr) != np.ndarray:
raise TypeError("Should be given a numpy array.")
else:
if (allow_nan is False):
if np.any(np.isnan(arr)):
raise ValueError("Array contains nan values.")
if (allow_inf is False):
if np.any(np.isinf(arr)):
raise ValueError("Array contains inf values.")
dtype = arr.dtype
if (np.issubdtype(dtype, np.integer) is False):
raise TypeError("Should be given array of integers.")
if (allow_neg is False):
if np.any(arr < 0):
raise ValueError("Array elements should be positive.")
return True