def abs_angle_diff(v_i, v_j):
"""
Returns the absolute value of the angle between two 3D vectors.
Parameters
----------
v_i : :obj:`numpy.ndarray`
the first 3D array
v_j : :obj:`numpy.ndarray`
the second 3D array
"""
# compute angle distance
dot_prod = min(max(v_i.dot(v_j), -1), 1)
angle_diff = np.arccos(dot_prod)
return np.abs(angle_diff)
# dictionary of distance functions
python类arccos()的实例源码
def test_branch_cuts(self):
# check branch cuts and continuity on them
yield _check_branch_cut, np.log, -0.5, 1j, 1, -1, True
yield _check_branch_cut, np.log2, -0.5, 1j, 1, -1, True
yield _check_branch_cut, np.log10, -0.5, 1j, 1, -1, True
yield _check_branch_cut, np.log1p, -1.5, 1j, 1, -1, True
yield _check_branch_cut, np.sqrt, -0.5, 1j, 1, -1, True
yield _check_branch_cut, np.arcsin, [ -2, 2], [1j, 1j], 1, -1, True
yield _check_branch_cut, np.arccos, [ -2, 2], [1j, 1j], 1, -1, True
yield _check_branch_cut, np.arctan, [0-2j, 2j], [1, 1], -1, 1, True
yield _check_branch_cut, np.arcsinh, [0-2j, 2j], [1, 1], -1, 1, True
yield _check_branch_cut, np.arccosh, [ -1, 0.5], [1j, 1j], 1, -1, True
yield _check_branch_cut, np.arctanh, [ -2, 2], [1j, 1j], 1, -1, True
# check against bogus branch cuts: assert continuity between quadrants
yield _check_branch_cut, np.arcsin, [0-2j, 2j], [ 1, 1], 1, 1
yield _check_branch_cut, np.arccos, [0-2j, 2j], [ 1, 1], 1, 1
yield _check_branch_cut, np.arctan, [ -2, 2], [1j, 1j], 1, 1
yield _check_branch_cut, np.arcsinh, [ -2, 2, 0], [1j, 1j, 1], 1, 1
yield _check_branch_cut, np.arccosh, [0-2j, 2j, 2], [1, 1, 1j], 1, 1
yield _check_branch_cut, np.arctanh, [0-2j, 2j, 0], [1, 1, 1j], 1, 1
def test_branch_cuts_complex64(self):
# check branch cuts and continuity on them
yield _check_branch_cut, np.log, -0.5, 1j, 1, -1, True, np.complex64
yield _check_branch_cut, np.log2, -0.5, 1j, 1, -1, True, np.complex64
yield _check_branch_cut, np.log10, -0.5, 1j, 1, -1, True, np.complex64
yield _check_branch_cut, np.log1p, -1.5, 1j, 1, -1, True, np.complex64
yield _check_branch_cut, np.sqrt, -0.5, 1j, 1, -1, True, np.complex64
yield _check_branch_cut, np.arcsin, [ -2, 2], [1j, 1j], 1, -1, True, np.complex64
yield _check_branch_cut, np.arccos, [ -2, 2], [1j, 1j], 1, -1, True, np.complex64
yield _check_branch_cut, np.arctan, [0-2j, 2j], [1, 1], -1, 1, True, np.complex64
yield _check_branch_cut, np.arcsinh, [0-2j, 2j], [1, 1], -1, 1, True, np.complex64
yield _check_branch_cut, np.arccosh, [ -1, 0.5], [1j, 1j], 1, -1, True, np.complex64
yield _check_branch_cut, np.arctanh, [ -2, 2], [1j, 1j], 1, -1, True, np.complex64
# check against bogus branch cuts: assert continuity between quadrants
yield _check_branch_cut, np.arcsin, [0-2j, 2j], [ 1, 1], 1, 1, False, np.complex64
yield _check_branch_cut, np.arccos, [0-2j, 2j], [ 1, 1], 1, 1, False, np.complex64
yield _check_branch_cut, np.arctan, [ -2, 2], [1j, 1j], 1, 1, False, np.complex64
yield _check_branch_cut, np.arcsinh, [ -2, 2, 0], [1j, 1j, 1], 1, 1, False, np.complex64
yield _check_branch_cut, np.arccosh, [0-2j, 2j, 2], [1, 1, 1j], 1, 1, False, np.complex64
yield _check_branch_cut, np.arctanh, [0-2j, 2j, 0], [1, 1, 1j], 1, 1, False, np.complex64
def test_against_cmath(self):
import cmath
points = [-1-1j, -1+1j, +1-1j, +1+1j]
name_map = {'arcsin': 'asin', 'arccos': 'acos', 'arctan': 'atan',
'arcsinh': 'asinh', 'arccosh': 'acosh', 'arctanh': 'atanh'}
atol = 4*np.finfo(np.complex).eps
for func in self.funcs:
fname = func.__name__.split('.')[-1]
cname = name_map.get(fname, fname)
try:
cfunc = getattr(cmath, cname)
except AttributeError:
continue
for p in points:
a = complex(func(np.complex_(p)))
b = cfunc(p)
assert_(abs(a - b) < atol, "%s %s: %s; cmath: %s" % (fname, p, a, b))
def interp_z(z0, z1, ratio, interp='linear'):
if interp == 'linear':
z_t = (1 - ratio) * z0 + ratio * z1
if interp == 'slerp':
N = len(z0)
z_t = []
for i in range(N):
z0_i = z0[i]
z1_i = z1[i]
z0_n = z0_i / np.linalg.norm(z0_i)
z1_n = z1_i / np.linalg.norm(z1_i)
omega = np.arccos(np.dot(z0_n, z1_n))
sin_omega = np.sin(omega)
if sin_omega == 0:
z_i = interp_z(z0_i, z1_i, ratio, 'linear')
else:
z_i = np.sin((1 - ratio) * omega) / sin_omega * z0_i + np.sin(ratio * omega) / sin_omega * z1_i
z_t.append(z_i[np.newaxis,...])
z_t = np.concatenate(z_t, axis=0)
return z_t
def adjust(self, rh, rv, frequency, eps_1, mu1):
# in place modification of rh and rv for the rough soil reflectivity model of Wegmüller & Mätzler (1999)
# Calculate ksigma = wavenumber*soilp%sigma(standard deviation of surface height)
ksigma = 2*np.pi*frequency*np.sqrt((1/2.9979e8)**2*eps_1) * self.roughness_rms
ksigma = ksigma.real
# Calculation of rh with ksigma
rh *= np.exp(-ksigma**(np.sqrt(0.1 * mu1))) # H pola
# calculation of rv with rh (the model is valid for angle between 0-70°
mask = mu1 < np.cos(60*np.pi/180)
rv[~mask] = rh[~mask] * mu1[~mask]**0.655 # <-- * ou ** ??
rv[mask] = rh[mask] * (0.635-0.0014*(np.arccos(mu1[mask])*180/np.pi-60))
def angles(self):
"""
Returns the angles (alpha, beta, gamma) of the lattice.
"""
if self._angles is None:
# Angles
angles = np.zeros(3)
for i in range(3):
j = (i + 1) % 3
k = (i + 2) % 3
angles[i] = np.dot(
self._matrix[j],
self._matrix[k]) / (self.lengths[j] * self.lengths[k])
angles = np.clip(angles, -1.0, 1.0)
self._angles = np.arccos(angles) * 180. / np.pi
return self._angles
def _plot_mpl(scheme):
# pylint: disable=relative-import, unused-variable
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.set_aspect('equal')
flt = numpy.vectorize(float)
pts = flt(scheme.points)
wgs = flt(scheme.weights)
for p, w in zip(pts, wgs):
# <https://en.wikipedia.org/wiki/Spherical_cap>
w *= 4 * numpy.pi
theta = numpy.arccos(1.0 - abs(w) / (2*numpy.pi))
color = '#1f77b4' if w >= 0 else '#d62728'
_plot_spherical_cap_mpl(ax, p, theta, color)
ax.set_axis_off()
return
def _draw_frame(self, frame):
for sph in frame.sphlist:
glPushMatrix()
glColor3f(*sph.color)
glTranslatef(sph.pos[0], sph.pos[1], sph.pos[2])
glutSolidSphere(sph.radius, self.sphere_slices, self.sphere_slices)
glPopMatrix()
for cyl in frame.cyllist:
p1, p2 = cyl.p1, cyl.p2
vec = p2 - p1
vlen = float(np.sqrt(np.square(vec).sum()))
angle = np.arccos(vec[2] / vlen) * 180.0 / np.pi
if vec[2] < 0:
angle = -angle
glPushMatrix()
glColor3f(*cyl.color)
glTranslatef(p1[0], p1[1], p1[2])
glRotatef(angle, -vec[1]*vec[2], vec[0]*vec[2], 0.0)
gluCylinder(self._quadric, 3.0, 3.0, vlen, 10, 10)
glPopMatrix()
def sample_sphere3d(radius=1., n_samples=1):
"""
Sample points from 3D sphere.
@param radius: radius of the sphere
@type radius: float
@param n_samples: number of samples to return
@type n_samples: int
@return: n_samples times random cartesian coordinates inside the sphere
@rtype: numpy array
"""
from numpy.random import random
from numpy import arccos, transpose, cos, sin, pi, power
r = radius * power(random(n_samples), 1 / 3.)
theta = arccos(2. * (random(n_samples) - 0.5))
phi = 2 * pi * random(n_samples)
x = cos(phi) * sin(theta) * r
y = sin(phi) * sin(theta) * r
z = cos(theta) * r
return transpose([x, y, z])
def polar3d(x):
"""
Polar coordinate representation of a three-dimensional vector.
@param x: vector (i.e. rank one array)
@return: polar coordinates (radius and polar angles)
"""
if x.shape != (3,):
raise ValueError(x)
r = norm(x)
theta = numpy.arccos(x[2] / r)
phi = numpy.arctan2(x[1], x[0])
return numpy.array([r, theta, phi])
def hav_dist(locs1, locs2):
"""
Return a distance matrix between two set of coordinates.
Use geometric distance (default) or haversine distance (if longlat=True).
Parameters
----------
locs1 : numpy.array
The first set of coordinates as [(long, lat), (long, lat)].
locs2 : numpy.array
The second set of coordinates as [(long, lat), (long, lat)].
Returns
-------
mat_dist : numpy.array
The distance matrix between locs1 and locs2
"""
locs1 = np.radians(locs1)
locs2 = np.radians(locs2)
cos_lat1 = np.cos(locs1[..., 0])
cos_lat2 = np.cos(locs2[..., 0])
cos_lat_d = np.cos(locs1[..., 0] - locs2[..., 0])
cos_lon_d = np.cos(locs1[..., 1] - locs2[..., 1])
return 6367000 * np.arccos(
cos_lat_d - cos_lat1 * cos_lat2 * (1 - cos_lon_d))
def compute_document_similarity(X):
'''
From a matrix of unit distances, computes the cosine similarity
then changes to the angular distance (for a proper metric).
'''
S = cdist(X, X, metric='cosine')
S -= 1
S *= -1
S[S > 1] = 1.0
S[S < 0] = 0.0
# Set nan values to zero
S[np.isnan(S)] = 0
# Convert to angular distance (a proper metric)
S = 1 - (np.arccos(S) / np.pi)
assert(not np.isnan(S).any())
assert(not np.isinf(S).any())
return S
def value(self, xyz):
xyz = xyz.reshape(-1,3)
a = self.a
b = self.b
c = self.c
# vector from first atom to central atom
vector1 = xyz[a] - xyz[b]
# vector from last atom to central atom
vector2 = xyz[c] - xyz[b]
# norm of the two vectors
norm1 = np.sqrt(np.sum(vector1**2))
norm2 = np.sqrt(np.sum(vector2**2))
dot = np.dot(vector1, vector2)
# Catch the edge case that very rarely this number is -1.
if dot / (norm1 * norm2) <= -1.0:
if (np.abs(dot / (norm1 * norm2)) + 1.0) < -1e-6:
raise RuntimeError('Encountered invalid value in angle')
return np.pi
if dot / (norm1 * norm2) >= 1.0:
if (np.abs(dot / (norm1 * norm2)) - 1.0) > 1e-6:
raise RuntimeError('Encountered invalid value in angle')
return 0.0
return np.arccos(dot / (norm1 * norm2))
def value(self, xyz):
xyz = xyz.reshape(-1,3)
a = np.array(self.a)
b = self.b
c = np.array(self.c)
xyza = np.mean(xyz[a], axis=0)
xyzc = np.mean(xyz[c], axis=0)
# vector from first atom to central atom
vector1 = xyza - xyz[b]
# vector from last atom to central atom
vector2 = xyzc - xyz[b]
# norm of the two vectors
norm1 = np.sqrt(np.sum(vector1**2))
norm2 = np.sqrt(np.sum(vector2**2))
dot = np.dot(vector1, vector2)
# Catch the edge case that very rarely this number is -1.
if dot / (norm1 * norm2) <= -1.0:
if (np.abs(dot / (norm1 * norm2)) + 1.0) < -1e-6:
raise RuntimeError('Encountered invalid value in angle')
return np.pi
return np.arccos(dot / (norm1 * norm2))
def calc_fac_dfac(q0):
"""
Calculate the prefactor mapping the quaternion to the exponential map
and also its derivative. Takes the first element of the quaternion only
"""
# Ill-defined around q0=1.0
qm1 = q0-1.0
# if np.abs(q0) == 1.0:
# fac = 2
# dfac = -2/3
if np.abs(qm1) < 1e-8:
fac = 2 - 2*qm1/3
dfac = -2/3
else:
fac = 2*np.arccos(q0)/np.sqrt(1-q0**2)
dfac = -2/(1-q0**2)
dfac += 2*q0*np.arccos(q0)/(1-q0**2)**1.5
return fac, dfac
def getProjectedAngleInXYPlane(self, z=0, ref_axis=[0,1], centre=[0,0], inDeg=True):
'''
Project the OA vector to z=z, calculate the XY position, construct a
2D vector from [centre] to this XY and measure the angle subtended by
this vector from [ref_axis] (clockwise).
'''
ref_axis = np.array(ref_axis)
centre = np.array(centre)
point_vector_from_fit_centre = np.array(self.getXY(z=z)) - centre
dotP = np.dot(ref_axis, point_vector_from_fit_centre)
crossP = np.cross(ref_axis, point_vector_from_fit_centre)
angle = np.arccos(dotP/(np.linalg.norm(ref_axis)*np.linalg.norm(point_vector_from_fit_centre)))
if np.sign(crossP) > 0:
angle = (np.pi-angle) + np.pi
if inDeg:
dir_v = self._eval_direction_vector()
return np.degrees(angle)
else:
return angle
def _box_vectors_to_lengths_angles(box_vectors):
unitcell_lengths = []
for basis in box_vectors:
unitcell_lengths.append(np.array([np.linalg.norm(frame_v) for frame_v in basis]))
unitcell_angles = []
for vs in box_vectors:
angles = np.array([np.degrees(
np.arccos(np.dot(vs[i], vs[j])/
(np.linalg.norm(vs[i]) * np.linalg.norm(vs[j]))))
for i, j in [(0,1), (1,2), (2,0)]])
unitcell_angles.append(angles)
unitcell_angles = np.array(unitcell_angles)
return unitcell_lengths, unitcell_angles
def angular(embedding, other_embedding=None):
"""Compute angular distance
Parameters
----------
embedding : (n_samples, n_dimension) numpy array
other_embedding : (n_other_samples, n_dimension, ) numpy array, optional
L2-normalized embeddings
Returns
-------
distance : (n_samples, n_other_samples) numpy array or float
Angular distance
"""
n_samples, _ = embedding.shape
if other_embedding is None:
other_embedding = embedding
return arccos(ag_np.stack(
ag_np.sum(embedding[i] * other_embedding, axis=1)
for i in range(n_samples)))
def test_fit_from_points(self):
# Set up a collection of points in the X-Y plane.
np.random.seed(0)
points = np.hstack([
np.random.random((100, 2)),
np.zeros(100).reshape(-1, 1)
])
plane = Plane.fit_from_points(points)
# The normal vector should be closely aligned with the Z-axis.
z_axis = np.array([0., 0., 1.])
angle = np.arccos(
np.dot(plane.normal, z_axis) /
np.linalg.norm(plane.normal)
)
self.assertTrue(angle % np.pi < 1e-6)
def _gather_stats(mesh):
# The cosines of the angles are the negative dot products of
# the normalized edges adjacent to the angle.
norms = numpy.sqrt(mesh.ei_dot_ei)
normalized_ei_dot_ej = numpy.array([
mesh.ei_dot_ej[0] / norms[1] / norms[2],
mesh.ei_dot_ej[1] / norms[2] / norms[0],
mesh.ei_dot_ej[2] / norms[0] / norms[1],
])
# pylint: disable=invalid-unary-operand-type
angles = numpy.arccos(-normalized_ei_dot_ej) / (2*numpy.pi) * 360.0
hist, bin_edges = numpy.histogram(
angles,
bins=numpy.linspace(0.0, 180.0, num=19, endpoint=True)
)
return hist, bin_edges
def gen_single_ps(self):
"""
Generate single point source, and return its data as a list.
"""
# Redshift and luminosity
self.z, self.lumo = self.get_lumo_redshift()
# angular diameter distance
self.param = PixelParams(self.z)
self.dA = self.param.dA
# W/Hz/Sr to Jy
dA = self.dA * 3.0856775814671917E+22 # Mpc to meter
self.lumo = self.lumo / dA**2 / (10.0**-24) # [Jy]
# Position
x = np.random.uniform(0, 1)
self.lat = (np.arccos(2 * x - 1) / np.pi * 180 - 90) # [deg]
self.lon = np.random.uniform(0, np.pi * 2) / np.pi * 180 # [deg]
# Radius
self.radius = self.param.get_angle(self.get_radius()) # [rad]
# Area
self.area = np.pi * self.radius**2 # [sr]
ps_list = [self.z, self.dA, self.lumo, self.lat, self.lon,
self.area, self.radius]
return ps_list
def gen_single_ps(self):
"""
Generate single point source, and return its data as a list.
"""
# Redshift and luminosity
self.z, self.lumo = self.get_lumo_redshift()
# angular diameter distance
self.param = PixelParams(self.z)
self.dA = self.param.dA
# W/Hz/Sr to Jy
dA = self.dA * 3.0856775814671917E+22 # Mpc to meter
self.lumo = self.lumo / dA**2 / (10.0**-24) # [Jy]
# Position
x = np.random.uniform(0, 1)
self.lat = (np.arccos(2 * x - 1) / np.pi * 180 - 90) # [deg]
self.lon = np.random.uniform(0, np.pi * 2) / np.pi * 180 # [deg]
# lobe
lobe = self.gen_lobe()
# Area
self.area = np.pi * self.lobe_maj * self.lobe_min
ps_list = [self.z, self.dA, self.lumo, self.lat, self.lon, self.area]
ps_list.extend(lobe)
return ps_list
def gen_single_ps(self):
"""
Generate single point source, and return its data as a list.
"""
# Redshift and luminosity
self.z, self.lumo = self.get_lumo_redshift()
# angular diameter distance
self.param = PixelParams(self.z)
self.dA = self.param.dA
# W/Hz/Sr to Jy
dA = self.dA * 3.0856775814671917E+22 # Mpc to meter
self.lumo = self.lumo / dA**2 / (10.0**-24) # [Jy]
# Position
x = np.random.uniform(0, 1)
self.lat = (np.arccos(2 * x - 1) / np.pi * 180 - 90) # [deg]
self.lon = np.random.uniform(0, np.pi * 2) / np.pi * 180 # [deg]
# Area
npix = hp.nside2npix(self.nside)
self.area = 4 * np.pi / npix # [sr]
ps_list = [self.z, self.dA, self.lumo, self.lat, self.lon, self.area]
return ps_list
def gen_single_ps(self):
"""
Generate single point source, and return its data as a list.
"""
# Redshift and luminosity
self.z, self.lumo = self.get_lumo_redshift()
# angular diameter distance
self.param = PixelParams(self.z)
self.dA = self.param.dA
# Radius
self.radius = self.param.get_angle(self.get_radius()) # [rad]
# W/Hz/Sr to Jy
dA = self.dA * 3.0856775814671917E+22 # Mpc to meter
self.lumo = self.lumo / dA**2 / (10.0**-24) # [Jy]
# Position
x = np.random.uniform(0, 1)
self.lat = (np.arccos(2 * x - 1) / np.pi * 180 - 90) # [deg]
self.lon = np.random.uniform(0, np.pi * 2) / np.pi * 180 # [deg]
# Area
self.area = np.pi * self.radius**2 # [sr] ?
ps_list = [self.z, self.dA, self.lumo, self.lat, self.lon,
self.area, self.radius]
return ps_list
calculateDistanceInt.py 文件源码
项目：RecursiveHierarchicalClustering
作者: xychang
项目源码
文件源码
阅读 28
收藏 0
点赞 0
评论 0
def run(self):
# pr = cProfile.Profile()
print '[LOG]: start new thread '+str(self.threadID)
curTime = time.time()
distM = self.matrix[self.sfrom].dot(
self.matrix[self.sto].T).todense()
distM = np.maximum(
np.arccos(np.minimum(distM, np.ones(distM.shape))) /
(PI_VALUE/200)-0.01,
np.zeros(distM.shape)).astype(np.int8)
# np.savetxt(self.fo, distM, fmt = '%d')
np.save(self.fo + '.npy', distM)
print('[LOG]: thread %d finished after %d' %
(self.threadID, time.time() - curTime))
# self.pr.disable()
# # sortby = 'cumulative'
# # pstats.Stats(pr).strip_dirs().sort_stats(sortby).print_stats()
# self.pr.print_stats()
def process_coords_for_computations(self, coords_for_computations, t):
"""
"""
if self._teffext:
return coords_for_computations
x, y, z, r = coords_for_computations[:,0], coords_for_computations[:,1], coords_for_computations[:,2], np.sqrt((coords_for_computations**2).sum(axis=1))
theta = np.arccos(z/r)
phi = np.arctan2(y, x)
xi_r = self._radamp * Y(self._m, self._l, theta, phi) * np.exp(-1j*2*np.pi*self._freq*t)
xi_t = self._tanamp * self.dYdtheta(self._m, self._l, theta, phi) * np.exp(-1j*2*np.pi*self._freq*t)
xi_p = self._tanamp/np.sin(theta) * self.dYdphi(self._m, self._l, theta, phi) * np.exp(-1j*2*np.pi*self._freq*t)
new_coords = np.zeros(coords_for_computations.shape)
new_coords[:,0] = coords_for_computations[:,0] + xi_r * np.sin(theta) * np.cos(phi)
new_coords[:,1] = coords_for_computations[:,1] + xi_r * np.sin(theta) * np.sin(phi)
new_coords[:,2] = coords_for_computations[:,2] + xi_r * np.cos(theta)
return new_coords
def test_branch_cuts(self):
# check branch cuts and continuity on them
yield _check_branch_cut, np.log, -0.5, 1j, 1, -1, True
yield _check_branch_cut, np.log2, -0.5, 1j, 1, -1, True
yield _check_branch_cut, np.log10, -0.5, 1j, 1, -1, True
yield _check_branch_cut, np.log1p, -1.5, 1j, 1, -1, True
yield _check_branch_cut, np.sqrt, -0.5, 1j, 1, -1, True
yield _check_branch_cut, np.arcsin, [ -2, 2], [1j, 1j], 1, -1, True
yield _check_branch_cut, np.arccos, [ -2, 2], [1j, 1j], 1, -1, True
yield _check_branch_cut, np.arctan, [0-2j, 2j], [1, 1], -1, 1, True
yield _check_branch_cut, np.arcsinh, [0-2j, 2j], [1, 1], -1, 1, True
yield _check_branch_cut, np.arccosh, [ -1, 0.5], [1j, 1j], 1, -1, True
yield _check_branch_cut, np.arctanh, [ -2, 2], [1j, 1j], 1, -1, True
# check against bogus branch cuts: assert continuity between quadrants
yield _check_branch_cut, np.arcsin, [0-2j, 2j], [ 1, 1], 1, 1
yield _check_branch_cut, np.arccos, [0-2j, 2j], [ 1, 1], 1, 1
yield _check_branch_cut, np.arctan, [ -2, 2], [1j, 1j], 1, 1
yield _check_branch_cut, np.arcsinh, [ -2, 2, 0], [1j, 1j, 1], 1, 1
yield _check_branch_cut, np.arccosh, [0-2j, 2j, 2], [1, 1, 1j], 1, 1
yield _check_branch_cut, np.arctanh, [0-2j, 2j, 0], [1, 1, 1j], 1, 1
def test_branch_cuts_complex64(self):
# check branch cuts and continuity on them
yield _check_branch_cut, np.log, -0.5, 1j, 1, -1, True, np.complex64
yield _check_branch_cut, np.log2, -0.5, 1j, 1, -1, True, np.complex64
yield _check_branch_cut, np.log10, -0.5, 1j, 1, -1, True, np.complex64
yield _check_branch_cut, np.log1p, -1.5, 1j, 1, -1, True, np.complex64
yield _check_branch_cut, np.sqrt, -0.5, 1j, 1, -1, True, np.complex64
yield _check_branch_cut, np.arcsin, [ -2, 2], [1j, 1j], 1, -1, True, np.complex64
yield _check_branch_cut, np.arccos, [ -2, 2], [1j, 1j], 1, -1, True, np.complex64
yield _check_branch_cut, np.arctan, [0-2j, 2j], [1, 1], -1, 1, True, np.complex64
yield _check_branch_cut, np.arcsinh, [0-2j, 2j], [1, 1], -1, 1, True, np.complex64
yield _check_branch_cut, np.arccosh, [ -1, 0.5], [1j, 1j], 1, -1, True, np.complex64
yield _check_branch_cut, np.arctanh, [ -2, 2], [1j, 1j], 1, -1, True, np.complex64
# check against bogus branch cuts: assert continuity between quadrants
yield _check_branch_cut, np.arcsin, [0-2j, 2j], [ 1, 1], 1, 1, False, np.complex64
yield _check_branch_cut, np.arccos, [0-2j, 2j], [ 1, 1], 1, 1, False, np.complex64
yield _check_branch_cut, np.arctan, [ -2, 2], [1j, 1j], 1, 1, False, np.complex64
yield _check_branch_cut, np.arcsinh, [ -2, 2, 0], [1j, 1j, 1], 1, 1, False, np.complex64
yield _check_branch_cut, np.arccosh, [0-2j, 2j, 2], [1, 1, 1j], 1, 1, False, np.complex64
yield _check_branch_cut, np.arctanh, [0-2j, 2j, 0], [1, 1, 1j], 1, 1, False, np.complex64
def test_against_cmath(self):
import cmath
points = [-1-1j, -1+1j, +1-1j, +1+1j]
name_map = {'arcsin': 'asin', 'arccos': 'acos', 'arctan': 'atan',
'arcsinh': 'asinh', 'arccosh': 'acosh', 'arctanh': 'atanh'}
atol = 4*np.finfo(np.complex).eps
for func in self.funcs:
fname = func.__name__.split('.')[-1]
cname = name_map.get(fname, fname)
try:
cfunc = getattr(cmath, cname)
except AttributeError:
continue
for p in points:
a = complex(func(np.complex_(p)))
b = cfunc(p)
assert_(abs(a - b) < atol, "%s %s: %s; cmath: %s" % (fname, p, a, b))
