|
5 | 5 | # system to decide whether to invest in loop peeling etc. Here we write |
6 | 6 | # the logical structure, and allow optimizing for size as well. |
7 | 7 |
|
8 | | -import math |
9 | | - |
10 | 8 | import numpy as np |
11 | | -from numba import njit |
12 | | - |
13 | | -# Define the np.ndarray type as an numpy array of 3 floats |
14 | | -# vec3d = np.empty(3, dtype=float) |
| 9 | +from numba import float64, int32, njit |
15 | 10 |
|
16 | | -# and 2 floats |
17 | | -# = np.empty(2, dtype=float) |
18 | 11 |
|
19 | | -@njit |
20 | | -def norm(x: float, y: float, z: float) -> float: |
21 | | - """Return the norm of a 3D vector given by 3 float components.""" |
22 | | - return vec_norm(vec_set(x, y, z)) |
| 12 | +@njit(float64(float64[:]),fastmath=True, cache=True, nogil=True) |
| 13 | +def vec_norm(vec: np.ndarray) -> float: |
| 14 | + """vec_norm() gives the norm of a vector.""" |
| 15 | + return np.sqrt(vec[0]**2 + vec[1]**2 + vec[2]**2) |
23 | 16 |
|
24 | | -@njit |
| 17 | +@njit(float64[:](float64,float64,float64),fastmath=True, cache=True, nogil=True) |
25 | 18 | def vec_set(x: float, y: float, z: float) -> np.ndarray: |
26 | 19 | """Set the components of a 3D vector from separate doubles.""" |
27 | 20 | return np.array([x, y, z]) |
28 | 21 |
|
| 22 | +@njit(float64(float64,float64,float64),fastmath=True, cache=True, nogil=True) |
| 23 | +def norm(x: float, y: float, z: float) -> float: |
| 24 | + """Return the norm of a 3D vector given by 3 float components.""" |
| 25 | + return vec_norm(vec_set(x, y, z)) |
29 | 26 |
|
| 27 | + |
| 28 | +@njit(float64[:](float64[:]),fastmath=True, cache=True, nogil=True) |
30 | 29 | def vec_copy(src: np.ndarray) -> np.ndarray: |
31 | 30 | """Copy one 3D vector into another.""" |
32 | 31 | return src.copy() |
33 | 32 |
|
34 | | -@njit |
| 33 | +@njit(float64[:](float64[:],float64[:]),fastmath=True, cache=True, nogil=True) |
35 | 34 | def vec_subt(from_: np.ndarray, sub: np.ndarray) -> np.ndarray: |
36 | 35 | """Subtract two 3D vectors.""" |
37 | 36 | return from_ - sub |
38 | 37 |
|
39 | | -@njit |
| 38 | +@njit(float64[:](float64[:],float64[:]),fastmath=True, cache=True, nogil=True) |
40 | 39 | def vec_add(vec1: np.ndarray, vec2: np.ndarray) -> np.ndarray: |
41 | 40 | """Add two 3D vectors.""" |
42 | 41 | return vec1 + vec2 |
43 | 42 |
|
44 | | -@njit |
| 43 | +@njit(float64[:](float64[:],float64),fastmath=True, cache=True, nogil=True) |
45 | 44 | def vec_scalar_mul(vec: np.ndarray, scalar: float) -> np.ndarray: |
46 | 45 | """vec_scalar_mul(np.ndarray, scalar) multiplies a vector by a scalar.""" |
47 | 46 | return vec * scalar |
48 | 47 |
|
49 | | -@njit |
| 48 | +@njit(float64(float64[:],float64[:]),fastmath=True, cache=True, nogil=True) |
50 | 49 | def vec_diff_norm(vec1: np.ndarray, vec2: np.ndarray) -> float: |
51 | 50 | """vec_diff_norm() gives the norm of the difference between two vectors.""" |
52 | | - # return np.linalg.norm(vec1 - vec2) |
53 | 51 | vec = vec1 - vec2 |
54 | | - return math.sqrt(vec[0]**2 + vec[1]**2 + vec[2]**2) |
| 52 | + return np.sqrt(vec[0]**2 + vec[1]**2 + vec[2]**2) |
55 | 53 |
|
56 | | -@njit |
57 | | -def vec_norm(vec: np.ndarray) -> float: |
58 | | - """vec_norm() gives the norm of a vector.""" |
59 | | - return math.sqrt(vec[0]**2 + vec[1]**2 + vec[2]**2) |
60 | 54 |
|
61 | | -@njit |
| 55 | +@njit(float64(float64[:],float64[:]),fastmath=True, cache=True, nogil=True) |
62 | 56 | def vec_dot(vec1: np.ndarray, vec2: np.ndarray) -> float: |
63 | 57 | """vec_dot() gives the dot product of two vectors as lists of floats.""" |
64 | 58 | # return np.dot(vec1, vec2) |
65 | 59 | return float(vec1[0]*vec2[0] + vec1[1]*vec2[1] + vec1[2]*vec2[2]) |
66 | 60 |
|
67 | | -@njit |
| 61 | +@njit(float64[:](float64[:],float64[:]),fastmath=True, cache=True, nogil=True) |
68 | 62 | def vec_cross(vec1: np.ndarray, vec2: np.ndarray) -> np.ndarray: |
69 | 63 | """Cross product of two vectors.""" |
70 | 64 | # return np.cross(vec1, vec2) |
71 | 65 | return np.array([vec1[1]*vec2[2] - vec1[2]*vec2[1], |
72 | 66 | vec1[2]*vec2[0] - vec1[0]*vec2[2], |
73 | 67 | vec1[0]*vec2[1] - vec1[1]*vec2[0]]) |
74 | 68 |
|
75 | | -@njit |
| 69 | +@njit("boolean(float64[:],float64[:],float64)",fastmath=True, cache=True, nogil=True) |
76 | 70 | def vec_cmp(vec1: np.ndarray, vec2: np.ndarray, tol: float = 1e-6) -> bool: |
77 | 71 | """vec_cmp() checks whether two vectors are equal within a tolerance.""" |
78 | 72 | return np.allclose(vec1, vec2, atol=tol) |
79 | 73 |
|
80 | | -@njit |
| 74 | +@njit(float64[:](float64[:]),fastmath=True, cache=True, nogil=True) |
81 | 75 | def unit_vector(vec: np.ndarray) -> np.ndarray: |
82 | 76 | """Normalize a vector to a unit vector.""" |
83 | 77 | magnitude = vec_norm(vec) |
84 | 78 | if magnitude == 0: |
85 | 79 | return vec # Avoid division by zero for zero vectors |
86 | 80 | return vec / magnitude |
87 | 81 |
|
88 | | -@njit |
89 | | -def vec_init(length=3) -> np.ndarray: |
| 82 | +@njit(float64[:](int32),fastmath=True, cache=True, nogil=True) |
| 83 | +def vec_init(length: int=3) -> np.ndarray: |
90 | 84 | """Initialize a vector to zero.""" |
91 | 85 | return np.zeros(length, dtype=float) |
0 commit comments