Skip to content

Commit 274fe8a

Browse files
feat(vectors): add support for jagged awkward arrays in vector operations
This commit introduces new functionality to handle jagged awkward arrays in both 2D and 3D vector operations. It includes the addition of new test cases to verify the behavior of vector operations with jagged arrays and an example script demonstrating practical usage in particle physics data analysis.
1 parent 58dcb9e commit 274fe8a

File tree

3 files changed

+183
-3
lines changed

3 files changed

+183
-3
lines changed

examples/07_jagged_arrays.py

+129
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,129 @@
1+
# Example demonstrating vector operations with jagged Awkward arrays
2+
import numpy as np
3+
try:
4+
import awkward as ak
5+
HAS_AWKWARD = True
6+
except ImportError:
7+
HAS_AWKWARD = False
8+
print("This example requires the awkward package. Please install it with 'pip install awkward'")
9+
exit(1)
10+
11+
from lvec import Vector2D, Vector3D
12+
13+
print("=== Jagged Arrays with Vector2D ===\n")
14+
15+
# Create jagged arrays for 2D vectors
16+
x_jagged = ak.Array([[1.0, 2.0], [3.0, 4.0, 5.0], [6.0]])
17+
y_jagged = ak.Array([[2.0, 3.0], [4.0, 5.0, 6.0], [7.0]])
18+
19+
# Create 2D vectors with jagged arrays
20+
vec2d_jagged = Vector2D(x_jagged, y_jagged)
21+
22+
print("Vector components:")
23+
print(f"x: {vec2d_jagged.x}")
24+
print(f"y: {vec2d_jagged.y}")
25+
26+
# Calculate properties
27+
print("\nVector properties:")
28+
print(f"Magnitude (r): {vec2d_jagged.r}")
29+
print(f"Azimuthal angle (phi): {vec2d_jagged.phi}")
30+
31+
# Vector operations with jagged arrays
32+
v2d_doubled = Vector2D(2 * x_jagged, 2 * y_jagged)
33+
34+
print("\nVector operations:")
35+
print(f"Original vector: ({vec2d_jagged.x}, {vec2d_jagged.y})")
36+
print(f"Doubled vector: ({v2d_doubled.x}, {v2d_doubled.y})")
37+
print(f"Sum: ({(vec2d_jagged + v2d_doubled).x}, {(vec2d_jagged + v2d_doubled).y})")
38+
39+
# Dot product with jagged arrays
40+
dot_product = vec2d_jagged.dot(v2d_doubled)
41+
print(f"Dot product: {dot_product}")
42+
43+
# Indexing into jagged arrays
44+
print("\nIndexing into first element of each subarray:")
45+
first_elements = vec2d_jagged[0]
46+
print(f"First elements: ({first_elements.x}, {first_elements.y})")
47+
48+
print("\n=== Jagged Arrays with Vector3D ===\n")
49+
50+
# Create jagged arrays for 3D vectors
51+
x_jagged = ak.Array([[1.0, 2.0], [3.0, 4.0, 5.0], [6.0]])
52+
y_jagged = ak.Array([[2.0, 3.0], [4.0, 5.0, 6.0], [7.0]])
53+
z_jagged = ak.Array([[3.0, 4.0], [5.0, 6.0, 7.0], [8.0]])
54+
55+
# Create 3D vectors with jagged arrays
56+
vec3d_jagged = Vector3D(x_jagged, y_jagged, z_jagged)
57+
58+
print("Vector components:")
59+
print(f"x: {vec3d_jagged.x}")
60+
print(f"y: {vec3d_jagged.y}")
61+
print(f"z: {vec3d_jagged.z}")
62+
63+
# Calculate properties
64+
print("\nVector properties:")
65+
print(f"Magnitude (r): {vec3d_jagged.r}")
66+
print(f"Cylindrical radius (rho): {vec3d_jagged.rho}")
67+
print(f"Azimuthal angle (phi): {vec3d_jagged.phi}")
68+
print(f"Polar angle (theta): {vec3d_jagged.theta}")
69+
70+
# Vector operations with jagged arrays
71+
v3d_doubled = Vector3D(2 * x_jagged, 2 * y_jagged, 2 * z_jagged)
72+
73+
print("\nVector operations:")
74+
print(f"Original vector: ({vec3d_jagged.x}, {vec3d_jagged.y}, {vec3d_jagged.z})")
75+
print(f"Doubled vector: ({v3d_doubled.x}, {v3d_doubled.y}, {v3d_doubled.z})")
76+
print(f"Sum: ({(vec3d_jagged + v3d_doubled).x}, {(vec3d_jagged + v3d_doubled).y}, {(vec3d_jagged + v3d_doubled).z})")
77+
78+
# Dot and cross products with jagged arrays
79+
dot_product = vec3d_jagged.dot(v3d_doubled)
80+
cross_product = vec3d_jagged.cross(v3d_doubled)
81+
82+
print(f"Dot product: {dot_product}")
83+
print(f"Cross product: ({cross_product.x}, {cross_product.y}, {cross_product.z})")
84+
85+
# Indexing into jagged arrays
86+
print("\nIndexing into first element of each subarray:")
87+
first_elements = vec3d_jagged[0]
88+
print(f"First elements: ({first_elements.x}, {first_elements.y}, {first_elements.z})")
89+
90+
print("\n=== Practical Example: Particle Physics Data ===\n")
91+
92+
# Simulate an event with varying number of particles per event
93+
nevents = 3
94+
95+
# Create jagged arrays for particle momenta in multiple events
96+
px = ak.Array([[10.0, 20.0, 30.0], [40.0, 50.0], [60.0, 70.0, 80.0, 90.0]])
97+
py = ak.Array([[15.0, 25.0, 35.0], [45.0, 55.0], [65.0, 75.0, 85.0, 95.0]])
98+
pz = ak.Array([[5.0, 10.0, 15.0], [20.0, 25.0], [30.0, 35.0, 40.0, 45.0]])
99+
100+
# Create 3D momentum vectors
101+
momenta = Vector3D(px, py, pz)
102+
103+
print(f"Number of events: {len(momenta.x)}")
104+
print(f"Particles per event: {ak.num(momenta.x)}")
105+
106+
# Calculate transverse momentum (pt) for all particles
107+
pt = momenta.rho
108+
print(f"\nTransverse momentum (pt): {pt}")
109+
110+
# Calculate total momentum per event
111+
total_px = ak.sum(momenta.x, axis=1)
112+
total_py = ak.sum(momenta.y, axis=1)
113+
total_pz = ak.sum(momenta.z, axis=1)
114+
115+
print("\nTotal momentum per event:")
116+
print(f"px: {total_px}")
117+
print(f"py: {total_py}")
118+
print(f"pz: {total_pz}")
119+
120+
# Create total momentum vector
121+
total_momentum = Vector3D(total_px, total_py, total_pz)
122+
print(f"\nTotal momentum magnitude: {total_momentum.r}")
123+
124+
# Find particles with pt > 50
125+
high_pt_mask = pt > 50
126+
high_pt_particles = momenta[high_pt_mask]
127+
128+
print(f"\nHigh-pt particles (pt > 50): {ak.num(high_pt_particles.x, axis=0)}")
129+
print(f"High-pt px values: {high_pt_particles.x}")

lvec/backends.py

+5-1
Original file line numberDiff line numberDiff line change
@@ -82,4 +82,8 @@ def backend_atan2(y, x, lib):
8282

8383
def backend_log(x, lib):
8484
"""Compute natural logarithm using appropriate backend."""
85-
return ak.log(x) if lib == 'ak' else np.log(x)
85+
return np.log(x)
86+
87+
def backend_exp(x, lib):
88+
"""Compute exponential using appropriate backend."""
89+
return np.exp(x)

lvec/tests/test_vectors.py

+49-2
Original file line numberDiff line numberDiff line change
@@ -22,13 +22,33 @@ def test_vec2d_numpy():
2222
assert np.all(np.isclose(v.r, np.sqrt(x**2 + y**2)))
2323

2424
def test_vec2d_awkward():
25-
"""Test 2D vector with awkward arrays"""
25+
"""Test 2D vector with regular awkward arrays"""
2626
x = ak.Array([[1.0, 2.0], [3.0, 4.0]])
2727
y = ak.Array([[2.0, 3.0], [4.0, 5.0]])
2828
v = Vec2D(x, y)
2929
assert ak.all(v.x == x)
3030
assert ak.all(v.y == y)
3131

32+
def test_vec2d_jagged_awkward():
33+
"""Test 2D vector with jagged awkward arrays"""
34+
x = ak.Array([[1.0, 2.0], [3.0, 4.0, 5.0]])
35+
y = ak.Array([[2.0, 3.0], [4.0, 5.0, 6.0]])
36+
v = Vec2D(x, y)
37+
assert ak.all(v.x == x)
38+
assert ak.all(v.y == y)
39+
40+
# Test vector operations
41+
v2 = Vec2D(2 * x, 2 * y)
42+
v_sum = v + v2
43+
assert ak.all(v_sum.x == 3 * x)
44+
assert ak.all(v_sum.y == 3 * y)
45+
46+
# Test magnitude and angle
47+
r = v.r
48+
phi = v.phi
49+
assert ak.all(np.sqrt(x**2 + y**2) == r)
50+
assert ak.all(np.arctan2(y, x) == phi)
51+
3252
def test_vec2d_operations():
3353
"""Test 2D vector operations"""
3454
v1 = Vec2D(1.0, 2.0)
@@ -67,7 +87,7 @@ def test_vec3d_numpy():
6787
assert np.all(np.isclose(v.r, np.sqrt(x**2 + y**2 + z**2)))
6888

6989
def test_vec3d_awkward():
70-
"""Test 3D vector with awkward arrays"""
90+
"""Test 3D vector with regular awkward arrays"""
7191
x = ak.Array([[1.0, 2.0], [3.0, 4.0]])
7292
y = ak.Array([[2.0, 3.0], [4.0, 5.0]])
7393
z = ak.Array([[3.0, 4.0], [5.0, 6.0]])
@@ -76,6 +96,33 @@ def test_vec3d_awkward():
7696
assert ak.all(v.y == y)
7797
assert ak.all(v.z == z)
7898

99+
def test_vec3d_jagged_awkward():
100+
"""Test 3D vector with jagged awkward arrays"""
101+
x = ak.Array([[1.0, 2.0], [3.0, 4.0, 5.0]])
102+
y = ak.Array([[2.0, 3.0], [4.0, 5.0, 6.0]])
103+
z = ak.Array([[3.0, 4.0], [5.0, 6.0, 7.0]])
104+
v = Vec3D(x, y, z)
105+
assert ak.all(v.x == x)
106+
assert ak.all(v.y == y)
107+
assert ak.all(v.z == z)
108+
109+
# Test vector operations
110+
v2 = Vec3D(2 * x, 2 * y, 2 * z)
111+
v_sum = v + v2
112+
assert ak.all(v_sum.x == 3 * x)
113+
assert ak.all(v_sum.y == 3 * y)
114+
assert ak.all(v_sum.z == 3 * z)
115+
116+
# Test magnitude and angles
117+
r = v.r
118+
rho = v.rho
119+
theta = v.theta
120+
phi = v.phi
121+
assert ak.all(np.sqrt(x**2 + y**2 + z**2) == r)
122+
assert ak.all(np.sqrt(x**2 + y**2) == rho)
123+
assert ak.all(np.arctan2(rho, z) == theta)
124+
assert ak.all(np.arctan2(y, x) == phi)
125+
79126
def test_vec3d_operations():
80127
"""Test 3D vector operations"""
81128
v1 = Vec3D(1.0, 2.0, 3.0)

0 commit comments

Comments
 (0)