From ac6d9d99aa37385a640234eedcde52599158483f Mon Sep 17 00:00:00 2001 From: Alan Pawlak Date: Fri, 20 Jan 2023 13:18:16 +0000 Subject: [PATCH 1/2] Ensure valid range for arccos function -1 to 1 --- src/sphstat/descriptives.py | 4 ++-- src/sphstat/distributions.py | 4 ++-- src/sphstat/twosample.py | 2 +- src/sphstat/utils.py | 4 ++-- 4 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/sphstat/descriptives.py b/src/sphstat/descriptives.py index 4db30eb..1b8fa2b 100644 --- a/src/sphstat/descriptives.py +++ b/src/sphstat/descriptives.py @@ -128,7 +128,7 @@ def rotationmatrix_withaxis(ua: np.array, ub: np.array) -> np.ndarray: ca /= np.linalg.norm(ca) B = ub.reshape((3, 1)) @ ca.reshape((1, 3)) B -= B.T - th = np.arccos(ab) + th = np.arccos(np.clip(ab, -1.0, 1.0)) A = np.eye(3) + np.sin(th) * B + (np.cos(th) - 1) * ((ub.reshape((3, 1)) @ ub.reshape((1, 3))) + ca.reshape((3, 1)) @ ca.reshape((1, 3))) return A @@ -381,7 +381,7 @@ def errmin(thph, sample): pts = sample['points'] err = 0 for pt in pts: - err += np.arccos(np.dot(pt, x)) + err += np.arccos(np.clip(np.dot(pt, x), -1.0, 1.0)) return err # We will use Nelder-Mead to minimise the arc lengths res = minimize(errmin, np.array([ms[0], ms[1]]), args = (sample), method='Nelder-Mead', tol=1e-6) diff --git a/src/sphstat/distributions.py b/src/sphstat/distributions.py index 904e01b..4dd8f37 100644 --- a/src/sphstat/distributions.py +++ b/src/sphstat/distributions.py @@ -289,7 +289,7 @@ def watson(numsamp: int, lamb: float, mu: float, nu: float, kappa: float) -> dic U = rng.uniform(0, 1, 1) V = rng.uniform(0, 1, 1) S = 1 / kappa * np.log(U / C + 1) - the = np.arccos(S) + the = np.arccos(np.clip(S, -1.0, 1.0)) phi = 2 * np.pi * rng.uniform(0, 1, 1) sample['tetas'].append(the) sample['phis'].append(phi) @@ -310,7 +310,7 @@ def watson(numsamp: int, lamb: float, mu: float, nu: float, kappa: float) -> dic U = rng.uniform(0, 1, 1) V = rng.uniform(0, 1, 1) S = (1 / c1) * np.tan(c2 * U) - the = np.arccos(S) + the = np.arccos(np.clip(S, -1.0, 1.0)) phi = 2 * np.pi * rng.uniform(0, 1, 1) sample['tetas'].append(the) sample['phis'].append(phi) diff --git a/src/sphstat/twosample.py b/src/sphstat/twosample.py index 2e92ea6..193aa38 100644 --- a/src/sphstat/twosample.py +++ b/src/sphstat/twosample.py @@ -898,5 +898,5 @@ def errors(samplecart: dict, srcpos: tuple) -> list: errs = [] srcvec = sph2cart(srcpos[0], srcpos[1]) for pt in samplecart['points']: - errs.append(np.arccos(np.dot(pt, srcvec))) + errs.append(np.arccos(np.clip(np.dot(pt, srcvec), -1.0. 1.0))) return errs diff --git a/src/sphstat/utils.py b/src/sphstat/utils.py index 98ffa89..2d25529 100644 --- a/src/sphstat/utils.py +++ b/src/sphstat/utils.py @@ -258,7 +258,7 @@ def angle(pt1: np.ndarray, pt2: np.ndarray) -> float: assert type(pt2) == np.ndarray assert np.isclose(np.linalg.norm(pt1), 1) assert np.isclose(np.linalg.norm(pt2), 1) - return np.arccos(np.dot(pt1, pt2)) + return np.arccos(np.clip(np.dot(pt1, pt2), -1.0, 1.0)) def cart2sph(pt: np.array) -> tuple: @@ -270,7 +270,7 @@ def cart2sph(pt: np.array) -> tuple: :rtype: tuple """ assert np.isclose(np.linalg.norm(pt), 1) - th = np.arccos(pt[2]) + th = np.arccos(np.clip(pt[2], -1.0, 1.0)) if pt[1] == 0: ph = 0 return th, ph From 034d292ec7b0793335afc510e5e851c7db5fe8a0 Mon Sep 17 00:00:00 2001 From: Alan Pawlak Date: Fri, 20 Jan 2023 13:39:03 +0000 Subject: [PATCH 2/2] Typo --- src/sphstat/twosample.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/sphstat/twosample.py b/src/sphstat/twosample.py index 193aa38..80e0542 100644 --- a/src/sphstat/twosample.py +++ b/src/sphstat/twosample.py @@ -898,5 +898,5 @@ def errors(samplecart: dict, srcpos: tuple) -> list: errs = [] srcvec = sph2cart(srcpos[0], srcpos[1]) for pt in samplecart['points']: - errs.append(np.arccos(np.clip(np.dot(pt, srcvec), -1.0. 1.0))) + errs.append(np.arccos(np.clip(np.dot(pt, srcvec), -1.0, 1.0))) return errs