Merge pull request #3 from EPFL-ENAC/refactor-scalars-use-math

Refactor scalars use math

pNeuma_simulator/contact_distance/contact_distance.py

@@ -1,3 +1,5 @@
+from math import atan2, cos, sin, sqrt
+
 import numpy as np
 from numba import jit
 
@@ -35,7 +37,7 @@ def calc_dtc(
     """
 
     # distance from i to j
-    s_i_j = np.sqrt((x_j - x_i) ** 2 + (y_j - y_i) ** 2)
+    s_i_j = sqrt((x_j - x_i) ** 2 + (y_j - y_i) ** 2)
     # distance of closest approach between i and j
     min_d = ellipses(l_j, w_j, l_i, w_i, x_j, y_j, x_i, y_i, theta_j, theta_i)
     dtc = s_i_j - min_d
@@ -78,43 +80,43 @@ def ellipses(
     Returns:
         float: distance between the centers when two ellipses are externally tangent
     """
-    theta3 = np.atan2(y2 - y1, x2 - x1)
-    cs1, sn1 = np.cos(theta1), np.sin(theta1)
-    cs3, sn3 = np.cos(theta3), np.sin(theta3)
-    k1d = np.cos(theta3 - theta1)
-    k2d = np.cos(theta3 - theta2)
-    k1k2 = np.cos(theta2 - theta1)
+    theta3 = atan2(y2 - y1, x2 - x1)
+    cs1, sn1 = cos(theta1), sin(theta1)
+    cs3, sn3 = cos(theta3), sin(theta3)
+    k1d = cos(theta3 - theta1)
+    k2d = cos(theta3 - theta2)
+    k1k2 = cos(theta2 - theta1)
     # eccentricity of ellipses
     e1 = 1 - (b1**2 / a1**2)
     e2 = 1 - (b2**2 / a2**2)
     # component of A'
     eta = a1 / b1 - 1
     a11 = b1**2 / b2**2 * (1 + 0.5 * (1 + k1k2) * (eta * (2 + eta) - e2 * (1 + eta * k1k2) ** 2))
-    a12 = b1**2 / b2**2 * 0.5 * np.sqrt(1 - k1k2**2) * (eta * (2 + eta) + e2 * (1 - eta**2 * k1k2**2))
+    a12 = b1**2 / b2**2 * 0.5 * sqrt(1 - k1k2**2) * (eta * (2 + eta) + e2 * (1 - eta**2 * k1k2**2))
     a22 = b1**2 / b2**2 * (1 + 0.5 * (1 - k1k2) * (eta * (2 + eta) - e2 * (1 - eta * k1k2) ** 2))
     # eigenvalues of A'
-    lambda1 = 0.5 * (a11 + a22) + 0.5 * np.sqrt((a11 - a22) ** 2 + 4 * a12**2)
-    lambda2 = 0.5 * (a11 + a22) - 0.5 * np.sqrt((a11 - a22) ** 2 + 4 * a12**2)
+    lambda1 = 0.5 * (a11 + a22) + 0.5 * sqrt((a11 - a22) ** 2 + 4 * a12**2)
+    lambda2 = 0.5 * (a11 + a22) - 0.5 * sqrt((a11 - a22) ** 2 + 4 * a12**2)
     # major and minor axes of transformed ellipse
-    b2p = 1 / np.sqrt(lambda1)
-    a2p = 1 / np.sqrt(lambda2)
+    b2p = 1 / sqrt(lambda1)
+    a2p = 1 / sqrt(lambda2)
     deltap = a2p**2 / b2p**2 - 1
     if abs(k1k2) == 1:
         if a11 > a22:
-            kpmp = 1 / np.sqrt(1 - e1 * k1d**2) * b1 / a1 * k1d
+            kpmp = 1 / sqrt(1 - e1 * k1d**2) * b1 / a1 * k1d
         else:
-            kpmp = (sn3 * cs1 - cs3 * sn1) / np.sqrt(1 - e1 * k1d**2)
+            kpmp = (sn3 * cs1 - cs3 * sn1) / sqrt(1 - e1 * k1d**2)
     elif deltap != 0:
         kpmp = (
-            a12 / np.sqrt(1 + k1k2) * (b1 / a1 * k1d + k2d + (b1 / a1 - 1) * k1d * k1k2)
-            + (lambda1 - a11) / np.sqrt(1 - k1k2) * (b1 / a1 * k1d - k2d - (b1 / a1 - 1) * k1d * k1k2)
-        ) / np.sqrt(2 * (a12**2 + (lambda1 - a11) ** 2) * (1 - e1 * k1d**2))
+            a12 / sqrt(1 + k1k2) * (b1 / a1 * k1d + k2d + (b1 / a1 - 1) * k1d * k1k2)
+            + (lambda1 - a11) / sqrt(1 - k1k2) * (b1 / a1 * k1d - k2d - (b1 / a1 - 1) * k1d * k1k2)
+        ) / sqrt(2 * (a12**2 + (lambda1 - a11) ** 2) * (1 - e1 * k1d**2))
     else:
         kpmp = 0
     if kpmp == 0 or deltap == 0:
         Rc = a2p + 1
         # The distance of closest approach
-        dist = Rc * b1 / np.sqrt(1 - e1 * k1d**2)
+        dist = Rc * b1 / sqrt(1 - e1 * k1d**2)
         return dist
     else:
         # coefficients of quartic for q
@@ -129,9 +131,9 @@ def ellipses(
         rts = np.real(np.roots(np.array([A, B, C, D, E], dtype=np.complex128)))
         qq = rts[rts > 0][0]
         # substitute for R'
-        Rc = np.sqrt(
+        Rc = sqrt(
             (qq**2 - 1) / deltap * (1 + b2p * (1 + deltap) / qq) ** 2 + (1 - (qq**2 - 1) / deltap) * (1 + b2p / qq) ** 2
         )
         # The distance of closest approach
-        dist = Rc * b1 / np.sqrt(1 - e1 * k1d**2)
+        dist = Rc * b1 / sqrt(1 - e1 * k1d**2)
         return dist

pNeuma_simulator/gang/navigation.py

@@ -1,3 +1,5 @@
+from math import exp
+
 import numpy as np
 from numba import jit
 from numpy import where
@@ -141,7 +143,7 @@ def decay(vel: float, theta: float) -> np.ndarray:
         alphas (numpy.ndarray): An array of decay angles.
     """
 
-    phi_max = int(np.exp(params.XM * norm(vel) * params.factor + params.CM))
+    phi_max = int(exp(params.XM * norm(vel) * params.factor + params.CM))
     # half degree resolution
     phi_range = np.linspace(phi_max, -phi_max, 4 * phi_max + 1)
     alphas = np.radians(phi_range) - theta

pNeuma_simulator/gang/particle.py

@@ -1,5 +1,7 @@
+from math import atan2, cos, degrees, sin
+
 from matplotlib.patches import Ellipse
-from numpy import array, atan2, cos, degrees, sin
+from numpy import array
 
 from pNeuma_simulator import params
 

pNeuma_simulator/initialization/equilibrium.py

@@ -1,4 +1,5 @@
-# from numpy import concatenate, repeat, cos, pi, isfinite, minimum, all
+from math import cos, pi
+
 import numpy as np
 from scipy.optimize import root_scalar
 from scipy.stats import distributions
@@ -118,7 +119,7 @@ def g(x):
         if total == lanes * L:
             f_eq = f(s_eq, lam, v0, d)
             # Adaptation time
-            prefactor = 1 + np.cos(2 * np.pi / n_veh)
+            prefactor = 1 + cos(2 * pi / n_veh)
             tau = 1 / (prefactor * f_eq)
             if all(np.isfinite(tau)):
                 tau = np.minimum(tau, 1 / lam)

pNeuma_simulator/initialization/poissondisc.py

@@ -1,4 +1,5 @@
 from copy import deepcopy
+from math import cos, pi, sin, sqrt
 from typing import List
 
 import numpy as np
@@ -250,10 +251,10 @@ def get_point(self, ref_pt: Particle):
         i = 0
         while i < self.k:
             # https://meyavuz.wordpress.com/2018/11/15/
-            rho = np.sqrt(self.rng.uniform(r**2, 4 * r**2))
-            theta = self.rng.uniform(0, 2 * np.pi)
-            x = ref_pt.x + rho * np.cos(theta)
-            y = ref_pt.y + rho / self.aspect * np.sin(theta)
+            rho = sqrt(self.rng.uniform(r**2, 4 * r**2))
+            theta = self.rng.uniform(0, 2 * pi)
+            x = ref_pt.x + rho * cos(theta)
+            y = ref_pt.y + rho / self.aspect * sin(theta)
             pt = Particle(x, y, 0.0, 0.0, "Moto")
             if not (
                 -self.width / 2 < pt.x < self.width / 2

pNeuma_simulator/simulate.py

@@ -1,6 +1,6 @@
 import warnings
 from copy import deepcopy
-from math import isinf, radians
+from math import cos, inf, isinf, pi, radians, sin
 from typing import Callable
 
 import numpy as np
@@ -97,8 +97,8 @@ def main(n_cars: int, n_moto: int, seed: int, parallel: Callable, COUNT: int = 5
         # Field of View analysis
         ##############################
         for image in agents + images:
-            cos_angle = np.cos(np.pi - image.theta)
-            sin_angle = np.sin(np.pi - image.theta)
+            cos_angle = cos(pi - image.theta)
+            sin_angle = sin(pi - image.theta)
             xc = params.xv - image.x
             yc = params.yv - image.y
             # https://stackoverflow.com/questions/37031356/
@@ -168,11 +168,11 @@ def main(n_cars: int, n_moto: int, seed: int, parallel: Callable, COUNT: int = 5
             # Distance from walls
             k_w = tangent_dist(theta_i, 0, l_i, w_i)
             if theta_i >= radians(0.5):
-                gap_w = (params.lane - y_i - k_w) / np.sin(theta_i)
+                gap_w = (params.lane - y_i - k_w) / sin(theta_i)
             elif theta_i <= -radians(0.5):
-                gap_w = (params.lane + y_i - k_w) / np.sin(-theta_i)
+                gap_w = (params.lane + y_i - k_w) / sin(-theta_i)
             else:
-                gap_w = np.inf
+                gap_w = inf
             interactions = agent.interactions
             if len(interactions) > 0:
                 gaps = []
@@ -212,9 +212,9 @@ def main(n_cars: int, n_moto: int, seed: int, parallel: Callable, COUNT: int = 5
                                 )
                             gap = s_i_j - min_d
                         else:
-                            gap = np.inf
+                            gap = inf
                     else:
-                        gap = np.inf
+                        gap = inf
                     gaps.append(gap)
                 if np.isfinite(gaps).sum() > 0:
                     leader = neighbors[np.argmin(gaps)]

pNeuma_simulator/utils/utils.py

@@ -1,3 +1,5 @@
+from math import cos, sin, sqrt
+
 import numpy as np
 from numba import jit
 
@@ -16,10 +18,10 @@ def direction(theta_i: float) -> tuple:
     """
 
     # Current direction vector
-    e_i = np.array([np.cos(theta_i), np.sin(theta_i)])
+    e_i = np.array([cos(theta_i), sin(theta_i)])
     # Normal vector to current direction
     # https://stackoverflow.com/questions/1243614/
-    e_i_n = np.array([np.sin(theta_i), -np.cos(theta_i)])
+    e_i_n = np.array([sin(theta_i), -cos(theta_i)])
     return e_i, e_i_n
 
 
@@ -55,10 +57,10 @@ def tangent_dist(theta_i: float, theta_j: float, a_i: float, b_i: float) -> floa
     Returns:
         k_prime (float): The tangent distance between theta_i and theta_j.
     """
-    if round(np.cos(theta_j - theta_i), 15) != 0:
-        c_prime = np.sin(theta_j - theta_i) / np.cos(theta_j - theta_i)
-        d_prime = np.sqrt(b_i**2 + a_i**2 * c_prime**2)
-        k_prime = d_prime / np.sqrt(c_prime**2 + 1)
+    if round(cos(theta_j - theta_i), 15) != 0:
+        c_prime = sin(theta_j - theta_i) / cos(theta_j - theta_i)
+        d_prime = sqrt(b_i**2 + a_i**2 * c_prime**2)
+        k_prime = d_prime / sqrt(c_prime**2 + 1)
     else:
         k_prime = a_i
     return k_prime