From 66111e506eeee3bd2ca8782e43c9df623761e858 Mon Sep 17 00:00:00 2001 From: Evan Kiefl Date: Mon, 3 Nov 2025 21:32:49 -0800 Subject: [PATCH 01/22] selected_event -> event --- pooltool/evolution/event_based/introspection.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pooltool/evolution/event_based/introspection.py b/pooltool/evolution/event_based/introspection.py index 61b50df8..e6b9c1d7 100644 --- a/pooltool/evolution/event_based/introspection.py +++ b/pooltool/evolution/event_based/introspection.py @@ -112,7 +112,7 @@ def _get_collision_events_from_cache( class SimulationSnapshot: step_number: int system: System - selected_event: Event + event: Event collision_cache: CollisionCache transition_cache: TransitionCache engine: PhysicsEngine @@ -230,7 +230,7 @@ def simulate_with_snapshots( snapshot = SimulationSnapshot( step_number=step, system=system_pre_evolve, - selected_event=event, + event=event, collision_cache=collision_cache_snapshot, transition_cache=transition_cache_snapshot, engine=engine, From 7207e42ed8ec5dc8921abfbea498dd5850922d39 Mon Sep 17 00:00:00 2001 From: Evan Kiefl Date: Mon, 3 Nov 2025 21:46:18 -0800 Subject: [PATCH 02/22] Improve make_kiss fn - Tests not passing --- pooltool/constants.py | 1 - pooltool/physics/resolve/ball_ball/core.py | 54 ++++++++++++++----- pooltool/physics/resolve/ball_cushion/core.py | 15 +++--- .../physics/resolve/transition/__init__.py | 10 ++-- .../event_based/test_introspection.py | 10 ++-- 5 files changed, 61 insertions(+), 29 deletions(-) diff --git a/pooltool/constants.py b/pooltool/constants.py index f180a69e..1ce73b30 100644 --- a/pooltool/constants.py +++ b/pooltool/constants.py @@ -12,7 +12,6 @@ np.set_printoptions(precision=16, suppress=True) EPS = np.finfo(float).eps * 100 -EPS_SPACE = 1e-9 # Ball states stationary: int = 0 diff --git a/pooltool/physics/resolve/ball_ball/core.py b/pooltool/physics/resolve/ball_ball/core.py index c0d49ea7..e825b805 100644 --- a/pooltool/physics/resolve/ball_ball/core.py +++ b/pooltool/physics/resolve/ball_ball/core.py @@ -1,7 +1,6 @@ from abc import ABC, abstractmethod from typing import Protocol -import pooltool.constants as const import pooltool.ptmath as ptmath from pooltool.objects.ball.datatypes import Ball @@ -28,18 +27,49 @@ class CoreBallBallCollision(ABC): def make_kiss(self, ball1: Ball, ball2: Ball) -> tuple[Ball, Ball]: """Translate the balls so they are (almost) touching - This makes a correction such that if the balls are not _exactly_ 2*R apart, they - are moved equally along their line of centers such that they are. Then, to avoid - downstream float precision round-off errors, a small epsilon of additional - distance (constants.EPS_SPACE) is put between them, ensuring the balls are - non-intersecting. - """ - r1, r2 = ball1.state.rvw[0], ball2.state.rvw[0] - n = ptmath.unit_vector(r2 - r1) + Uses binary search to find a time offset that positions the balls at the target + separation distance. The balls are moved equally along their line of centers + (traced forward/backward in time along their velocity vectors) until they are + separated by 2*R + spacer, where spacer is a small epsilon to avoid float + precision errors. + + Args: + ball1: First ball in the collision + ball2: Second ball in the collision - correction = 2 * ball1.params.R - ptmath.norm3d(r2 - r1) + const.EPS_SPACE - ball2.state.rvw[0] += correction / 2 * n - ball1.state.rvw[0] -= correction / 2 * n + Returns: + Modified ball1 and ball2 with adjusted positions + """ + r1 = ball1.state.rvw[0] + r2 = ball2.state.rvw[0] + v1 = ball1.state.rvw[1] + v2 = ball2.state.rvw[1] + + tmag = ball1.params.R / max(ptmath.norm3d(v1), ptmath.norm3d(v2)) + tmin = -tmag + tmax = tmag + + spacer = 1e-12 + distance_target = 1e-15 + + max_iter = 100 + for _ in range(max_iter): + t = (tmin + tmax) / 2 + r1 = r1 - t * v1 + r2 = r2 - t * v2 + distance = ptmath.norm3d(r2 - r1) + error = distance - (2 * ball1.params.R + spacer) + + if abs(error) < distance_target: + break + + if error > 0: + tmax = t + else: + tmin = t + + ball1.state.rvw[0] = r1 + ball2.state.rvw[0] = r2 return ball1, ball2 diff --git a/pooltool/physics/resolve/ball_cushion/core.py b/pooltool/physics/resolve/ball_cushion/core.py index af06bb59..d67c5158 100644 --- a/pooltool/physics/resolve/ball_cushion/core.py +++ b/pooltool/physics/resolve/ball_cushion/core.py @@ -3,7 +3,6 @@ import numpy as np -import pooltool.constants as const import pooltool.ptmath as ptmath from pooltool.objects.ball.datatypes import Ball from pooltool.objects.table.components import ( @@ -57,7 +56,7 @@ def make_kiss(self, ball: Ball, cushion: LinearCushionSegment) -> Ball: This makes a correction such that if the ball is not a distance R from the cushion, the ball is moved along the normal such that it is. To avoid downstream float precision round-off error, a small epsilon of additional distance - (constants.EPS_SPACE) is put between them, ensuring the cushion and ball are + (`spacer`) is put between them, ensuring the cushion and ball are separated post-resolution. """ normal = cushion.get_normal_xy(ball.state.rvw) @@ -72,10 +71,10 @@ def make_kiss(self, ball: Ball, cushion: LinearCushionSegment) -> Ball: ) c[2] = ball.state.rvw[0, 2] + spacer = 1e-9 + # Move the ball to exactly meet the cushion - correction = ( - ball.params.R - ptmath.norm3d(ball.state.rvw[0] - c) + const.EPS_SPACE - ) + correction = ball.params.R - ptmath.norm3d(ball.state.rvw[0] - c) + spacer ball.state.rvw[0] -= correction * normal return ball @@ -107,7 +106,7 @@ def make_kiss(self, ball: Ball, cushion: CircularCushionSegment) -> Ball: This makes a correction such that if the ball is not a distance R from the cushion, the ball is moved along the normal such that it is. To avoid downstream float precision round-off error, a small epsilon of additional distance - (constants.EPS_SPACE) is put between them, ensuring the cushion and ball are + (`spacer`) is put between them, ensuring the cushion and ball are separated post-resolution. """ normal = cushion.get_normal_xy(ball.state.rvw) @@ -115,12 +114,14 @@ def make_kiss(self, ball: Ball, cushion: CircularCushionSegment) -> Ball: # orient the normal so it points away from playing surface normal = normal if np.dot(normal, ball.state.rvw[1]) > 0 else -normal + spacer = 1e-9 + c = np.array([cushion.center[0], cushion.center[1], ball.state.rvw[0, 2]]) correction = ( ball.params.R + cushion.radius - ptmath.norm3d(ball.state.rvw[0] - c) - - const.EPS_SPACE + - spacer ) ball.state.rvw[0] += correction * normal diff --git a/pooltool/physics/resolve/transition/__init__.py b/pooltool/physics/resolve/transition/__init__.py index eebb63c5..90845cae 100644 --- a/pooltool/physics/resolve/transition/__init__.py +++ b/pooltool/physics/resolve/transition/__init__.py @@ -15,6 +15,8 @@ from pooltool.objects.ball.datatypes import Ball from pooltool.physics.resolve.models import BallTransitionModel +TOLERANCE = 1e-12 + class BallTransitionStrategy(Protocol): """Ball transition models must satisfy this protocol""" @@ -47,8 +49,8 @@ def resolve(self, ball: Ball, transition: EventType, inplace: bool = False) -> B # angular velocity components are nearly 0. Then set them to exactly 0. v = ball.state.rvw[1] w = ball.state.rvw[2] - assert (np.abs(v) < const.EPS_SPACE).all() - assert (np.abs(w[:2]) < const.EPS_SPACE).all() + assert (np.abs(v) < TOLERANCE).all() + assert (np.abs(w[:2]) < TOLERANCE).all() ball.state.rvw[1, :] = [0.0, 0.0, 0.0] ball.state.rvw[2, :2] = [0.0, 0.0] @@ -58,8 +60,8 @@ def resolve(self, ball: Ball, transition: EventType, inplace: bool = False) -> B # set them to exactly 0. v = ball.state.rvw[1] w = ball.state.rvw[2] - assert (np.abs(v) < const.EPS_SPACE).all() - assert (np.abs(w) < const.EPS_SPACE).all() + assert (np.abs(v) < TOLERANCE).all() + assert (np.abs(w) < TOLERANCE).all() ball.state.rvw[1, :] = [0.0, 0.0, 0.0] ball.state.rvw[2, :] = [0.0, 0.0, 0.0] diff --git a/tests/evolution/event_based/test_introspection.py b/tests/evolution/event_based/test_introspection.py index 216924e0..f751a65d 100644 --- a/tests/evolution/event_based/test_introspection.py +++ b/tests/evolution/event_based/test_introspection.py @@ -35,8 +35,8 @@ def test_selected_event_in_all_possible_events(): all_events = snapshot.get_prospective_events() first_event = all_events[0] - assert snapshot.selected_event.event_type == first_event.event_type - assert snapshot.selected_event.time == first_event.time + assert snapshot.event.event_type == first_event.event_type + assert snapshot.event.time == first_event.time def test_pre_evolve_equals_snapshot_system(): @@ -59,7 +59,7 @@ def test_post_evolve_advances_time(): for step in range(len(seq)): snapshot = seq[step] - event = snapshot.selected_event + event = snapshot.event post_evolve = snapshot.post_evolve_system(event) assert post_evolve.t == event.time @@ -75,7 +75,7 @@ def test_post_resolve_of_n_equals_pre_evolve_of_n_plus_1(): next_snapshot = seq[step + 1] post_resolve = current_snapshot.post_resolve_system( - current_snapshot.selected_event + current_snapshot.event ) pre_evolve_next = next_snapshot.pre_evolve_system() @@ -89,7 +89,7 @@ def test_system_state_progression(): for step in range(len(seq)): snapshot = seq[step] - event = snapshot.selected_event + event = snapshot.event pre_evolve = snapshot.pre_evolve_system() post_evolve = snapshot.post_evolve_system(event) From 935cfe9a33a0b153431e85e6e7cd9e7811a977ef Mon Sep 17 00:00:00 2001 From: Evan Kiefl Date: Wed, 5 Nov 2025 20:56:00 -0800 Subject: [PATCH 03/22] Intersecting balls are now events - This is made possible by improvements to the make_kiss function - In brief, a spacer between balls is enforced by translating colliding balls along their trajectories until the spacer is achieved. --- pooltool/evolution/event_based/simulate.py | 4 +- pooltool/physics/resolve/ball_ball/core.py | 197 +++++++++++++++--- pooltool/ptmath/roots/quadratic.py | 49 +++++ .../event_based/test_introspection.py | 4 +- tests/evolution/event_based/test_simulate.py | 24 +-- 5 files changed, 234 insertions(+), 44 deletions(-) diff --git a/pooltool/evolution/event_based/simulate.py b/pooltool/evolution/event_based/simulate.py index 5b013d18..dfc321fd 100755 --- a/pooltool/evolution/event_based/simulate.py +++ b/pooltool/evolution/event_based/simulate.py @@ -100,6 +100,7 @@ def step(self) -> Event: if self.max_events > 0 and self.num_events > self.max_events: self.shot.stop_balls() + self.shot._update_history(null_event(time=self.shot.t)) self.done = True self.num_events += 1 @@ -367,8 +368,7 @@ def get_next_ball_ball_collision( ptmath.norm3d(ball1_state.rvw[0] - ball2_state.rvw[0]) < ball1_params.R + ball2_params.R ): - # If balls are intersecting, avoid internal collisions - cache[ball_pair] = np.inf + cache[ball_pair] = shot.t else: ball_pairs.append(ball_pair) collision_coeffs.append( diff --git a/pooltool/physics/resolve/ball_ball/core.py b/pooltool/physics/resolve/ball_ball/core.py index e825b805..f41c3143 100644 --- a/pooltool/physics/resolve/ball_ball/core.py +++ b/pooltool/physics/resolve/ball_ball/core.py @@ -1,6 +1,8 @@ from abc import ABC, abstractmethod from typing import Protocol +import numpy as np + import pooltool.ptmath as ptmath from pooltool.objects.ball.datatypes import Ball @@ -25,51 +27,190 @@ class CoreBallBallCollision(ABC): """Operations used by every ball-ball collision resolver""" def make_kiss(self, ball1: Ball, ball2: Ball) -> tuple[Ball, Ball]: - """Translate the balls so they are (almost) touching + """Position balls at precise target separation before collision resolution. + + This method adjusts ball positions so they are separated by exactly 2*R + + spacer, where R is the ball radius and ``spacer`` is a small epsilon to prevent + ball intersection that occurs due to floating-point precision if an explicit + spacer is not added. + + The primary method solves a quadratic equation to find the time offset that + positions the balls at the target separation. Balls are moved along their + trajectories (position + velocity * time) to this configuration. Quadratic terms + are assumed negligible. + + If the required movement exceeds 10x the spacer, which can occur if balls are + moving with nearly the same velocity, a naive fallback strategy is used that + moves the balls uniformly along the line of centers until they're separated by + an amount ``spacer``. + + Algorithm: + 1. Calculate quadratic coefficients for separation equation + 2. Solve for time offset that achieves target separation + 3. Move balls to corrected positions: r_new = r + t * v + 4. If movement > 10x spacer (fallback): + - Identify which ball is "chasing" (higher line-of-centers velocity) and + which is being "chased" (lower line-of-centers velocity) + - Apply full displacement to chased ball, other stays fixed - Uses binary search to find a time offset that positions the balls at the target - separation distance. The balls are moved equally along their line of centers - (traced forward/backward in time along their velocity vectors) until they are - separated by 2*R + spacer, where spacer is a small epsilon to avoid float - precision errors. + Returns: + tuple[Ball, Ball]: + ``ball1`` and ``ball2`` modified in place with adjusted positions. + """ + r1 = ball1.state.rvw[0] + r2 = ball2.state.rvw[0] + v1 = ball1.state.rvw[1] + v2 = ball2.state.rvw[1] + + spacer = 1e-6 + initial_distance = ptmath.norm3d(r2 - r1) + target_distance = 2 * ball1.params.R + spacer + + Bx = v2[0] - v1[0] + By = v2[1] - v1[1] + Bz = v2[2] - v1[2] + Cx = r2[0] - r1[0] + Cy = r2[1] - r1[1] + Cz = r2[2] - r1[2] + alpha = Bx**2 + By**2 + Bz**2 + beta = 2 * Bx * Cx + 2 * By * Cy + 2 * Bz * Cz + gamma = Cx**2 + Cy**2 + Cz**2 - (2 * ball1.params.R + spacer) ** 2 + roots = ptmath.roots.quadratic.solve_complex(alpha, beta, gamma) + + imag_mag = np.abs(roots.imag) + real_mag = np.abs(roots.real) + keep = (imag_mag / real_mag) < 1e-3 + roots = roots[keep].real + t = roots[np.abs(roots).argmin()] + + r1_corrected = r1 + t * v1 + r2_corrected = r2 + t * v2 + + # This fallback exists for cases when the balls are moving nearly in unison + # (i.e. nearly same velocity). In these cases, the amount of time to create a + # distance `spacer` between them can be high enough to significantly displace + # both balls. + for r, r_corrected in zip([r1, r2], [r1_corrected, r2_corrected]): + if ptmath.norm3d(r - r_corrected) > 10 * spacer: + line_of_centers = (r2 - r1) / initial_distance + total_displacement = target_distance - initial_distance + + # Determine which ball is chasing (has higher radial velocity). + v1_radial = np.dot(v1, line_of_centers) + v2_radial = np.dot(v2, line_of_centers) + + if v1_radial > v2_radial: + # Ball 1 is chasing ball 2, push ball 2 forward. + r1_corrected = r1 + r2_corrected = r2 + total_displacement * line_of_centers + else: + # Ball 2 is chasing ball 1, push ball 1 backward. + r1_corrected = r1 - total_displacement * line_of_centers + r2_corrected = r2 + + break + + ball1.state.rvw[0] = r1_corrected + ball2.state.rvw[0] = r2_corrected + + return ball1, ball2 + + def resolve_continually_touching( + self, ball1: Ball, ball2: Ball + ) -> tuple[Ball, Ball]: + """Prevent repeated collision detection for nearly-touching balls moving in unison. + + This method is called to handle rare cases where balls are moving with very + similar velocities. This can happen in some edge cases when frozen balls in a + perfect line are given energy along their line, (e.g. Newton's cradle). Without + intervention, the balls repeatedly trigger events microseconds apart that stall + progression of the simulation, sometimes indefinitely, via an explosion of + events. + + This is an unfortunate consequence of modeling non-instantaneous multibody + collisions using instantaneous pairwise collisions, and resolving it requires + some phenomonelogical intervention that hopefully appears to be realistic, + despite it not being grounded in theory. + + The solution is applied in ths method, and uses a momentum transfer mechanism: + the "chased" ball (slower in the line-of-centers direction) steals a fraction of + the "chaser's" radial momentum. This creates gradual separation over time, + preventing the balls from triggering repeated collision events while maintaining + physically plausible behavior. + + Algorithm: + 1. Projects velocities onto line of centers to get radial components + 2. If radial relative velocity is below threshold (< 1mm/s): + - Identifies which ball is "chasing" (higher radial velocity) + - Chased ball steals 10% of chaser's radial momentum + - Chaser loses this momentum, chased gains it + 3. Tangential velocity components remain unchanged Args: ball1: First ball in the collision ball2: Second ball in the collision Returns: - Modified ball1 and ball2 with adjusted positions + Modified ball1 and ball2 with adjusted velocities + + Notes: + - Practically speaking, this is a no-op method for all but the most + contrived simulation conditions. For one such condition, see + `sandbox/newtons_cradle.py` """ r1 = ball1.state.rvw[0] r2 = ball2.state.rvw[0] v1 = ball1.state.rvw[1] v2 = ball2.state.rvw[1] - tmag = ball1.params.R / max(ptmath.norm3d(v1), ptmath.norm3d(v2)) - tmin = -tmag - tmax = tmag + v1_speed = ptmath.norm3d(v1) + v2_speed = ptmath.norm3d(v2) + both_moving = v1_speed > 0 and v2_speed > 0 - spacer = 1e-12 - distance_target = 1e-15 + if not both_moving: + return ball1, ball2 - max_iter = 100 - for _ in range(max_iter): - t = (tmin + tmax) / 2 - r1 = r1 - t * v1 - r2 = r2 - t * v2 - distance = ptmath.norm3d(r2 - r1) - error = distance - (2 * ball1.params.R + spacer) + theft_fraction = 0.1 + velocity_similarity_threshold = 0.9 - if abs(error) < distance_target: - break + line_of_centers = ptmath.unit_vector(r2 - r1) + + # Velocities projected onto the line of centers (loc). + v1_loc = np.dot(v1, line_of_centers) + v2_loc = np.dot(v2, line_of_centers) + + cosine_similarity = np.dot(v1, v2) / (v1_speed * v2_speed) + velocities_aligned = cosine_similarity > velocity_similarity_threshold - if error > 0: - tmax = t + if abs(v2_loc - v1_loc) < 0.01 and velocities_aligned: + if v1_loc > v2_loc: + chaser_loc_vel = v1_loc + ball1_is_chaser = True else: - tmin = t + chaser_loc_vel = v2_loc + ball1_is_chaser = False - ball1.state.rvw[0] = r1 - ball2.state.rvw[0] = r2 + # Chased ball steals fraction of chaser's line of centers momentum + # FIXME: We assume equal mass, so transfer velocity directly + theft_fraction = theft_fraction + stolen_loc_velocity = chaser_loc_vel * theft_fraction + + if ball1_is_chaser: + v1_loc_new = v1_loc - stolen_loc_velocity + v2_loc_new = v2_loc + stolen_loc_velocity + else: + v1_loc_new = v1_loc + stolen_loc_velocity + v2_loc_new = v2_loc - stolen_loc_velocity + + v1_corrected = v1 - v1_loc * line_of_centers + v1_loc_new * line_of_centers + v2_corrected = v2 - v2_loc * line_of_centers + v2_loc_new * line_of_centers + + momentum_before = v1 + v2 + momentum_after = v1_corrected + v2_corrected + assert np.allclose(momentum_before, momentum_after, rtol=1e-10) + + ball1.state.rvw[1] = v1_corrected + ball2.state.rvw[1] = v2_corrected return ball1, ball2 @@ -81,8 +222,10 @@ def resolve( ball2 = ball2.copy() ball1, ball2 = self.make_kiss(ball1, ball2) + ball1, ball2 = self.solve(ball1, ball2) + ball1, ball2 = self.resolve_continually_touching(ball1, ball2) - return self.solve(ball1, ball2) + return ball1, ball2 @abstractmethod def solve(self, ball1: Ball, ball2: Ball) -> tuple[Ball, Ball]: diff --git a/pooltool/ptmath/roots/quadratic.py b/pooltool/ptmath/roots/quadratic.py index 2f8c885b..0c770832 100644 --- a/pooltool/ptmath/roots/quadratic.py +++ b/pooltool/ptmath/roots/quadratic.py @@ -1,7 +1,9 @@ +import cmath import math import numpy as np from numba import jit +from numpy.typing import NDArray import pooltool.constants as const @@ -19,3 +21,50 @@ def solve(a: float, b: float, c: float) -> tuple[float, float]: u1 = (-bp - delta**0.5) / a u2 = -u1 - b / a return u1, u2 + + +# TODO: In the branch `3d`, which will eventually be merged into main, this function has +# replaced `solve`. +@jit(nopython=True, cache=const.use_numba_cache) +def solve_complex(a: float, b: float, c: float) -> NDArray[np.complex128]: + _a = complex(a) + _b = complex(b) + _c = complex(c) + + roots = np.full(2, np.nan, dtype=np.complex128) + + if abs(_a) != 0: + # Quadratic case + d = _b * _b - 4 * _a * _c + sqrt_d = cmath.sqrt(d) + + # Sign trick to reduce catastrophic cancellation + sign_b = 1.0 if _b.real >= 0 else -1.0 + + r1_num = -_b - sign_b * sqrt_d + r1_den = 2 * _a + + # Fallback if numerator is tiny + if abs(r1_num) < 1e-14 * abs(r1_den): + r1_num = -_b + sign_b * sqrt_d + + r1 = r1_num / r1_den + + # Use product identity for x2 + if abs(r1) < 1e-14: + r2 = (-_b + sqrt_d) / (2 * _a) + else: + r2 = (_c / _a) / r1 + + roots[0] = r1 + roots[1] = r2 + return roots + + if abs(_b) != 0: + # Linear case + r1 = -_c / _b + roots[0] = r1 + return roots + + # Equation is just c=0. Either zero or infinite solutions. Returns nans + return roots diff --git a/tests/evolution/event_based/test_introspection.py b/tests/evolution/event_based/test_introspection.py index f751a65d..875b6476 100644 --- a/tests/evolution/event_based/test_introspection.py +++ b/tests/evolution/event_based/test_introspection.py @@ -74,9 +74,7 @@ def test_post_resolve_of_n_equals_pre_evolve_of_n_plus_1(): current_snapshot = seq[step] next_snapshot = seq[step + 1] - post_resolve = current_snapshot.post_resolve_system( - current_snapshot.event - ) + post_resolve = current_snapshot.post_resolve_system(current_snapshot.event) pre_evolve_next = next_snapshot.pre_evolve_system() assert post_resolve == pre_evolve_next diff --git a/tests/evolution/event_based/test_simulate.py b/tests/evolution/event_based/test_simulate.py index 9cba0f99..0d4879e3 100644 --- a/tests/evolution/event_based/test_simulate.py +++ b/tests/evolution/event_based/test_simulate.py @@ -434,16 +434,15 @@ def true_time_to_collision(eps, V0, mu_r, g): @pytest.mark.parametrize( "solver", [quartic.QuarticSolver.NUMERIC, quartic.QuarticSolver.HYBRID] ) -def test_no_ball_ball_collisions_for_intersecting_balls(solver: quartic.QuarticSolver): - """Two already intersecting balls don't collide +def test_ball_ball_collisions_for_intersecting_balls(solver: quartic.QuarticSolver): + """Two already intersecting balls collide. - In this instance, no further collision is detected because the balls are already - intersecting. Otherwise perpetual internal collisions occur, keeping the two balls - locked. + Previously, intersecting balls were prevented from colliding to avoid perpetual + internal collisions. Now, with the improved make_kiss implementation, intersecting + balls are properly separated and collide normally. - This test doesn't make sure that balls don't intersect, it tests the safeguard that - prevents already intersecting balls from colliding with their internal walls, which - keeps them intersected like links in a chain. + This test verifies that intersecting balls are detected as a collision at time == + shot.t , - ~ , , - ~ , , ' ' ,, ' ' , @@ -478,12 +477,13 @@ def test_no_ball_ball_collisions_for_intersecting_balls(solver: quartic.QuarticS _assert_rolling(system.balls["cue"].state.rvw, system.balls["cue"].params.R) assert ( - get_next_event(system, quartic_solver=solver).event_type != EventType.BALL_BALL + get_next_event(system, quartic_solver=solver).event_type == EventType.BALL_BALL ) - assert ( - get_next_ball_ball_collision(system, CollisionCache(), solver=solver).time - == np.inf + collision_event = get_next_ball_ball_collision( + system, CollisionCache(), solver=solver ) + assert collision_event.time != np.inf + assert collision_event.time == 0 def test_ball_history_immutability(): From 94591f63e4e2a67249c1c7797073b27969906f4a Mon Sep 17 00:00:00 2001 From: Evan Kiefl Date: Wed, 5 Nov 2025 20:59:43 -0800 Subject: [PATCH 04/22] Add newton's cradle - Breaks when range(5) --- sandbox/newtons_cradle.py | 32 ++++++++++++++++++++++++++++++++ 1 file changed, 32 insertions(+) create mode 100644 sandbox/newtons_cradle.py diff --git a/sandbox/newtons_cradle.py b/sandbox/newtons_cradle.py new file mode 100644 index 00000000..c2d041be --- /dev/null +++ b/sandbox/newtons_cradle.py @@ -0,0 +1,32 @@ +import pooltool as pt + +table = pt.Table.default() +ball_radius = pt.BallParams.default().R + +balls = [ + pt.Ball.create( + str(i + 1), + xy=(0.5 * table.w, 0.4 * table.l + 2 * i * ball_radius), + ballset=pt.objects.BallSet("pooltool_pocket"), + ) + for i in range(4) +] + +balls.append( + pt.Ball.create( + "cue", + xy=(0.5 * table.w, 0.1 * table.l), + ballset=pt.objects.BallSet("pooltool_pocket"), + ) +) + +system = pt.System( + balls=balls, + table=table, + cue=pt.Cue.default(), +) + +engine = pt.physics.PhysicsEngine() +system.strike(V0=2, phi=pt.aim.at_ball(system, "1", cut=0)) +pt.simulate(system, inplace=True, engine=engine) +pt.show(system) From 5ab06d3fc1e170266dd3f128de16b928a1ebb378 Mon Sep 17 00:00:00 2001 From: Evan Kiefl Date: Fri, 14 Nov 2025 08:27:06 -0800 Subject: [PATCH 05/22] Remove ** operators in numba fns --- pooltool/evolution/event_based/solve.py | 28 +++++++++++++------------ 1 file changed, 15 insertions(+), 13 deletions(-) diff --git a/pooltool/evolution/event_based/solve.py b/pooltool/evolution/event_based/solve.py index 5ac90fea..c1a80ed1 100644 --- a/pooltool/evolution/event_based/solve.py +++ b/pooltool/evolution/event_based/solve.py @@ -154,11 +154,11 @@ def ball_ball_collision_coeffs( Bx, By = b2x - b1x, b2y - b1y Cx, Cy = c2x - c1x, c2y - c1y - a = Ax**2 + Ay**2 + a = Ax * Ax + Ay * Ay b = 2 * Ax * Bx + 2 * Ay * By - c = Bx**2 + 2 * Ax * Cx + 2 * Ay * Cy + By**2 + c = Bx * Bx + 2 * Ax * Cx + 2 * Ay * Cy + By * By d = 2 * Bx * Cx + 2 * By * Cy - e = Cx**2 + Cy**2 - 4 * R**2 + e = Cx * Cx + Cy * Cy - 4 * R * R return a, b, c, d, e @@ -237,16 +237,16 @@ def ball_linear_cushion_collision_time( B = lx * bx + ly * by if direction == 0: - C = l0 + lx * cx + ly * cy + R * np.sqrt(lx**2 + ly**2) + C = l0 + lx * cx + ly * cy + R * np.sqrt(lx * lx + ly * ly) root1, root2 = ptmath.roots.quadratic.solve(A, B, C) roots = [root1, root2] elif direction == 1: - C = l0 + lx * cx + ly * cy - R * np.sqrt(lx**2 + ly**2) + C = l0 + lx * cx + ly * cy - R * np.sqrt(lx * lx + ly * ly) root1, root2 = ptmath.roots.quadratic.solve(A, B, C) roots = [root1, root2] else: - C1 = l0 + lx * cx + ly * cy + R * np.sqrt(lx**2 + ly**2) - C2 = l0 + lx * cx + ly * cy - R * np.sqrt(lx**2 + ly**2) + C1 = l0 + lx * cx + ly * cy + R * np.sqrt(lx * lx + ly * ly) + C2 = l0 + lx * cx + ly * cy - R * np.sqrt(lx * lx + ly * ly) root1, root2 = ptmath.roots.quadratic.solve(A, B, C1) root3, root4 = ptmath.roots.quadratic.solve(A, B, C2) roots = [root1, root2, root3, root4] @@ -308,11 +308,13 @@ def ball_circular_cushion_collision_coeffs( bx, by = v * cos_phi, v * sin_phi cx, cy = rvw[0, 0], rvw[0, 1] - A = 0.5 * (ax**2 + ay**2) + A = 0.5 * (ax * ax + ay * ay) B = ax * bx + ay * by - C = ax * (cx - a) + ay * (cy - b) + 0.5 * (bx**2 + by**2) + C = ax * (cx - a) + ay * (cy - b) + 0.5 * (bx * bx + by * by) D = bx * (cx - a) + by * (cy - b) - E = 0.5 * (a**2 + b**2 + cx**2 + cy**2 - (r + R) ** 2) - (cx * a + cy * b) + E = 0.5 * (a * a + b * b + cx * cx + cy * cy - (r + R) * (r + R)) - ( + cx * a + cy * b + ) return A, B, C, D, E @@ -381,11 +383,11 @@ def ball_pocket_collision_coeffs( bx, by = v * cos_phi, v * sin_phi cx, cy = rvw[0, 0], rvw[0, 1] - A = 0.5 * (ax**2 + ay**2) + A = 0.5 * (ax * ax + ay * ay) B = ax * bx + ay * by - C = ax * (cx - a) + ay * (cy - b) + 0.5 * (bx**2 + by**2) + C = ax * (cx - a) + ay * (cy - b) + 0.5 * (bx * bx + by * by) D = bx * (cx - a) + by * (cy - b) - E = 0.5 * (a**2 + b**2 + cx**2 + cy**2 - r**2) - (cx * a + cy * b) + E = 0.5 * (a * a + b * b + cx * cx + cy * cy - r * r) - (cx * a + cy * b) return A, B, C, D, E From 8a3002a03367eab042e7ccfbf7c6cfd82c24cdb5 Mon Sep 17 00:00:00 2001 From: Evan Kiefl Date: Mon, 17 Nov 2025 07:42:24 -0800 Subject: [PATCH 06/22] Use is_overlapping --- pooltool/evolution/event_based/simulate.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/pooltool/evolution/event_based/simulate.py b/pooltool/evolution/event_based/simulate.py index d1356c94..7b956064 100755 --- a/pooltool/evolution/event_based/simulate.py +++ b/pooltool/evolution/event_based/simulate.py @@ -353,9 +353,8 @@ def get_next_ball_ball_collision( and ball2_state.s in const.nontranslating ): cache[ball_pair] = np.inf - elif ( - ptmath.norm3d(ball1_state.rvw[0] - ball2_state.rvw[0]) - < ball1_params.R + ball2_params.R + elif ptmath.is_overlapping( + ball1_state.rvw, ball2_state.rvw, ball1_params.R, ball2_params.R, ): cache[ball_pair] = shot.t else: From ce71146ae4dda909e69bb850967ef24d0a2fef8d Mon Sep 17 00:00:00 2001 From: Evan Kiefl Date: Mon, 17 Nov 2025 07:43:05 -0800 Subject: [PATCH 07/22] Avoid use of ** operator --- pooltool/ptmath/utils.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pooltool/ptmath/utils.py b/pooltool/ptmath/utils.py index ec3b0181..7c50d003 100644 --- a/pooltool/ptmath/utils.py +++ b/pooltool/ptmath/utils.py @@ -266,7 +266,7 @@ def point_on_line_closest_to_point( @jit(nopython=True, cache=const.use_numba_cache) def squared_norm3d(vec: NDArray[np.float64]) -> float: """Calculate the squared norm of a 3D vector""" - return vec[0] ** 2 + vec[1] ** 2 + vec[2] ** 2 + return vec[0] * vec[0] + vec[1] * vec[1] + vec[2] * vec[2] @jit(nopython=True, cache=const.use_numba_cache) @@ -290,7 +290,7 @@ def norm3d(vec: NDArray[np.float64]) -> float: @jit(nopython=True, cache=const.use_numba_cache) def squared_norm2d(vec: NDArray[np.float64]) -> float: """Calculate the squared norm of a 2D vector""" - return vec[0] ** 2 + vec[1] ** 2 + return vec[0] * vec[0] + vec[1] * vec[1] @jit(nopython=True, cache=const.use_numba_cache) From a3f28562a0ae34d6b8ddfa691cd716c2fc826c65 Mon Sep 17 00:00:00 2001 From: Evan Kiefl Date: Mon, 17 Nov 2025 08:57:56 -0800 Subject: [PATCH 08/22] Add speed dependent spacer -- cradle tests pass --- pooltool/physics/resolve/ball_ball/core.py | 26 ++-- tests/evolution/event_based/test_simulate.py | 141 +++++++++++++++++++ 2 files changed, 159 insertions(+), 8 deletions(-) diff --git a/pooltool/physics/resolve/ball_ball/core.py b/pooltool/physics/resolve/ball_ball/core.py index f41c3143..7e5e4466 100644 --- a/pooltool/physics/resolve/ball_ball/core.py +++ b/pooltool/physics/resolve/ball_ball/core.py @@ -62,7 +62,12 @@ def make_kiss(self, ball1: Ball, ball2: Ball) -> tuple[Ball, Ball]: v1 = ball1.state.rvw[1] v2 = ball2.state.rvw[1] - spacer = 1e-6 + approach_speed = np.dot(v1 - v2, ptmath.unit_vector(r2 - r1)) + + spacer_min = 0.5e-6 + spacer_max = 1e-6 + speed_factor = max(0, min(approach_speed, 10)) + spacer = spacer_min + speed_factor * (spacer_max - spacer_min) initial_distance = ptmath.norm3d(r2 - r1) target_distance = 2 * ball1.params.R + spacer @@ -72,15 +77,20 @@ def make_kiss(self, ball1: Ball, ball2: Ball) -> tuple[Ball, Ball]: Cx = r2[0] - r1[0] Cy = r2[1] - r1[1] Cz = r2[2] - r1[2] - alpha = Bx**2 + By**2 + Bz**2 + alpha = Bx * Bx + By * By + Bz * Bz beta = 2 * Bx * Cx + 2 * By * Cy + 2 * Bz * Cz - gamma = Cx**2 + Cy**2 + Cz**2 - (2 * ball1.params.R + spacer) ** 2 - roots = ptmath.roots.quadratic.solve_complex(alpha, beta, gamma) - - imag_mag = np.abs(roots.imag) - real_mag = np.abs(roots.real) + gamma = ( + Cx * Cx + + Cy * Cy + + Cz * Cz + - (2 * ball1.params.R + spacer) * (2 * ball1.params.R + spacer) + ) + roots_complex = ptmath.roots.quadratic.solve_complex(alpha, beta, gamma) + + imag_mag = np.abs(roots_complex.imag) + real_mag = np.abs(roots_complex.real) keep = (imag_mag / real_mag) < 1e-3 - roots = roots[keep].real + roots = roots_complex[keep].real t = roots[np.abs(roots).argmin()] r1_corrected = r1 + t * v1 diff --git a/tests/evolution/event_based/test_simulate.py b/tests/evolution/event_based/test_simulate.py index b5000c8b..5ed1a679 100644 --- a/tests/evolution/event_based/test_simulate.py +++ b/tests/evolution/event_based/test_simulate.py @@ -4,6 +4,8 @@ import pooltool.constants as const import pooltool.ptmath as ptmath +from pooltool.evolution.event_based.introspection import simulate_with_snapshots +from pooltool import aim, events from pooltool.events import EventType, ball_ball_collision, ball_pocket_collision from pooltool.evolution.event_based.cache import CollisionCache from pooltool.evolution.event_based.simulate import ( @@ -14,6 +16,8 @@ ) from pooltool.evolution.event_based.solve import ball_ball_collision_time from pooltool.objects import Ball, BilliardTableSpecs, Cue, Table +from pooltool.objects.ball.params import BallParams +from pooltool.objects.ball.sets import BallSet from pooltool.ptmath.roots import quadratic from pooltool.system import System from tests.evolution.event_based.test_data import TEST_DIR @@ -540,3 +544,140 @@ def test_stick_ball_event_detection(): assert simulated.events[1].event_type == EventType.STICK_BALL assert simulated.events[1].time == 0 assert simulated.events[2].time > 0 + + +def test_newtons_cradle_backspin(): + """Test Newtons's cradle when incoming ball has backspin. + + This is a easier simulation scenario than test_newtons_cradle_rolling_spin because + the incoming ball's spin, after the first collision, pulls it away from the line of + centers, avoiding a chain of followup collisions that ultimately push the whole line + of balls forwards. + """ + # Create a newton's cradle system, where balls 1, 2, and 3 are in a line. the cue + # ball is placed on the same line of centers some distance away and has momentum + # towards the 1 ball. + table = Table.default() + ball_radius = BallParams.default().R + + balls = [ + Ball.create( + str(i + 1), + xy=(0.5 * table.w, 0.4 * table.l + 2 * i * ball_radius), + ballset=BallSet("pooltool_pocket"), + ) + for i in range(3) + ] + + balls.append( + Ball.create( + "cue", + xy=(0.5 * table.w, 0.1 * table.l), + ballset=BallSet("pooltool_pocket"), + ) + ) + + system = System( + balls=balls, + table=table, + cue=Cue.default(), + ) + + system.strike(V0=2, b=-0.5, phi=aim.at_ball(system, "1", cut=0)) + system = simulate(system, max_events=10) + + # Define a 1-millisecond window that starts from the first collision. + first_collision = events.filter_type(system.events, EventType.BALL_BALL)[0] + start_time = first_collision.time - 1e-9 + end_time = first_collision.time + 1e-3 + + collision_chain = events.filter_events( + system.events, + events.by_time(start_time, after=True), + events.by_time(end_time, after=False), + events.by_type(EventType.BALL_BALL), + ) + + # The last ball should be involved in exactly one collision, since a near-elastic + # collision should effectively halt the second-last ball while sending the last ball + # flying away, precluding followup collisions. + assert len(events.filter_ball(collision_chain, "3")) == 1 + + # We don't assert how the second wave collision waves play out, but we know the + # collision sequence up until the last ball is hit: a chain of events where momentum + # is transfered sequentially through the line of balls. + expected_collision_order = [{"cue", "1"}, {"1", "2"}, {"2", "3"}] + + assert len(collision_chain) >= len(expected_collision_order) + + for idx, expected in enumerate(expected_collision_order): + actual = set(collision_chain[idx].ids) + assert expected == actual + + +def test_newtons_cradle_rolling_spin(): + """Test Newtons's cradle when incoming ball has rolling spin. + + This is a harder simulation scenario than test_newtons_cradle_backspin because the + incoming ball's spin, after the first collision, triggers many followup collisions + that ultimately push the who line of balls forwards. + """ + # Create a newton's cradle system, where balls 1, 2, and 3 are in a line. the cue + # ball is placed on the same line of centers some distance away and has momentum + # towards the 1 ball. + table = Table.default() + ball_radius = BallParams.default().R + + balls = [ + Ball.create( + str(i + 1), + xy=(0.5 * table.w, 0.4 * table.l + 2 * i * ball_radius), + ballset=BallSet("pooltool_pocket"), + ) + for i in range(3) + ] + + balls.append( + Ball.create( + "cue", + xy=(0.5 * table.w, 0.1 * table.l), + ballset=BallSet("pooltool_pocket"), + ) + ) + + system = System( + balls=balls, + table=table, + cue=Cue.default(), + ) + + system.strike(V0=2, phi=aim.at_ball(system, "1", cut=0)) + simulate(system, inplace=True, max_events=10) + + # Define a 1-millisecond window that starts from the first collision. + first_collision = events.filter_type(system.events, EventType.BALL_BALL)[0] + start_time = first_collision.time - 1e-9 + end_time = first_collision.time + 1e-3 + + collision_chain = events.filter_events( + system.events, + events.by_time(start_time, after=True), + events.by_time(end_time, after=False), + events.by_type(EventType.BALL_BALL), + ) + + # The last ball should be involved in exactly one collision, since a near-elastic + # collision should effectively halt the second-last ball while sending the last ball + # flying away, precluding followup collisions. + assert len(events.filter_ball(collision_chain, "3")) == 1 + + # We don't assert how the second wave collision waves play out, but we know the + # collision sequence up until the last ball is hit: a chain of events where momentum + # is transfered sequentially through the line of balls. + expected_collision_order = [{"cue", "1"}, {"1", "2"}, {"2", "3"}] + + assert len(collision_chain) >= len(expected_collision_order) + + for idx, expected in enumerate(expected_collision_order): + actual = set(collision_chain[idx].ids) + assert expected == actual From edb8aac670080c145e80324e05d61314844126b0 Mon Sep 17 00:00:00 2001 From: Evan Kiefl Date: Sun, 14 Dec 2025 17:34:00 -0800 Subject: [PATCH 09/22] Checkpoint --- .../evolution/event_based/introspection.py | 15 +------ pooltool/physics/resolve/ball_ball/core.py | 34 ++++++++------- sandbox/newtons_cradle.py | 6 +-- .../event_based/test_introspection.py | 10 ++--- .../resolve/ball_ball/test_ball_ball.py | 41 ++++++++++++++++++- 5 files changed, 69 insertions(+), 37 deletions(-) diff --git a/pooltool/evolution/event_based/introspection.py b/pooltool/evolution/event_based/introspection.py index 6b6d21af..16d0d204 100644 --- a/pooltool/evolution/event_based/introspection.py +++ b/pooltool/evolution/event_based/introspection.py @@ -111,7 +111,7 @@ def _get_collision_events_from_cache( class SimulationSnapshot: step_number: int system: System - event: Event + next_event: Event collision_cache: CollisionCache transition_cache: TransitionCache engine: PhysicsEngine @@ -227,7 +227,7 @@ def simulate_with_snapshots( snapshot = SimulationSnapshot( step_number=step, system=system_pre_evolve, - event=event, + next_event=event, collision_cache=collision_cache_snapshot, transition_cache=transition_cache_snapshot, engine=engine, @@ -240,14 +240,3 @@ def simulate_with_snapshots( step += 1 return sim.shot, snapshot_sequence - - -if __name__ == "__main__": - from pathlib import Path - - import pooltool as pt - - output = Path("test.json") - simulate_with_snapshots(pt.System.example(), output) - seq = SimulationSnapshotSequence.load(output) - system = pt.simulate(pt.System.example()) diff --git a/pooltool/physics/resolve/ball_ball/core.py b/pooltool/physics/resolve/ball_ball/core.py index 7e5e4466..8290d645 100644 --- a/pooltool/physics/resolve/ball_ball/core.py +++ b/pooltool/physics/resolve/ball_ball/core.py @@ -36,8 +36,8 @@ def make_kiss(self, ball1: Ball, ball2: Ball) -> tuple[Ball, Ball]: The primary method solves a quadratic equation to find the time offset that positions the balls at the target separation. Balls are moved along their - trajectories (position + velocity * time) to this configuration. Quadratic terms - are assumed negligible. + trajectories (position + velocity * time) to this configuration. Acceleration + terms are assumed negligible. If the required movement exceeds 10x the spacer, which can occur if balls are moving with nearly the same velocity, a naive fallback strategy is used that @@ -101,22 +101,26 @@ def make_kiss(self, ball1: Ball, ball2: Ball) -> tuple[Ball, Ball]: # distance `spacer` between them can be high enough to significantly displace # both balls. for r, r_corrected in zip([r1, r2], [r1_corrected, r2_corrected]): - if ptmath.norm3d(r - r_corrected) > 10 * spacer: + if ptmath.norm3d(r - r_corrected) > 2 * spacer: line_of_centers = (r2 - r1) / initial_distance total_displacement = target_distance - initial_distance - # Determine which ball is chasing (has higher radial velocity). - v1_radial = np.dot(v1, line_of_centers) - v2_radial = np.dot(v2, line_of_centers) + ## Determine which ball is chasing (has higher radial velocity). + # v1_radial = np.dot(v1, line_of_centers) + # v2_radial = np.dot(v2, line_of_centers) - if v1_radial > v2_radial: - # Ball 1 is chasing ball 2, push ball 2 forward. - r1_corrected = r1 - r2_corrected = r2 + total_displacement * line_of_centers - else: - # Ball 2 is chasing ball 1, push ball 1 backward. - r1_corrected = r1 - total_displacement * line_of_centers - r2_corrected = r2 + # if v1_radial > v2_radial: + # # Ball 2 is chasing ball 1, pull 2 backwards. + # r1_corrected = r1 + # r2_corrected = r2 + total_displacement * line_of_centers + # else: + # # Ball 1 is chasing ball 2, pull 1 backwards. + # r1_corrected = r1 - total_displacement * line_of_centers + # r2_corrected = r2 + + correction = 2 * ball1.params.R - ptmath.norm3d(r2 - r1) + spacer + r1_corrected = r1 - correction / 2 * ptmath.unit_vector(r2 - r1) + r2_corrected = r2 + correction / 2 * ptmath.unit_vector(r2 - r1) break @@ -233,7 +237,7 @@ def resolve( ball1, ball2 = self.make_kiss(ball1, ball2) ball1, ball2 = self.solve(ball1, ball2) - ball1, ball2 = self.resolve_continually_touching(ball1, ball2) + # ball1, ball2 = self.resolve_continually_touching(ball1, ball2) return ball1, ball2 diff --git a/sandbox/newtons_cradle.py b/sandbox/newtons_cradle.py index c2d041be..719b853e 100644 --- a/sandbox/newtons_cradle.py +++ b/sandbox/newtons_cradle.py @@ -1,4 +1,5 @@ import pooltool as pt +from pooltool.evolution.event_based.introspection import simulate_with_snapshots table = pt.Table.default() ball_radius = pt.BallParams.default().R @@ -9,7 +10,7 @@ xy=(0.5 * table.w, 0.4 * table.l + 2 * i * ball_radius), ballset=pt.objects.BallSet("pooltool_pocket"), ) - for i in range(4) + for i in range(3) ] balls.append( @@ -28,5 +29,4 @@ engine = pt.physics.PhysicsEngine() system.strike(V0=2, phi=pt.aim.at_ball(system, "1", cut=0)) -pt.simulate(system, inplace=True, engine=engine) -pt.show(system) +system, snapshots = simulate_with_snapshots(system) diff --git a/tests/evolution/event_based/test_introspection.py b/tests/evolution/event_based/test_introspection.py index 875b6476..186aa3b0 100644 --- a/tests/evolution/event_based/test_introspection.py +++ b/tests/evolution/event_based/test_introspection.py @@ -35,8 +35,8 @@ def test_selected_event_in_all_possible_events(): all_events = snapshot.get_prospective_events() first_event = all_events[0] - assert snapshot.event.event_type == first_event.event_type - assert snapshot.event.time == first_event.time + assert snapshot.next_event.event_type == first_event.event_type + assert snapshot.next_event.time == first_event.time def test_pre_evolve_equals_snapshot_system(): @@ -59,7 +59,7 @@ def test_post_evolve_advances_time(): for step in range(len(seq)): snapshot = seq[step] - event = snapshot.event + event = snapshot.next_event post_evolve = snapshot.post_evolve_system(event) assert post_evolve.t == event.time @@ -74,7 +74,7 @@ def test_post_resolve_of_n_equals_pre_evolve_of_n_plus_1(): current_snapshot = seq[step] next_snapshot = seq[step + 1] - post_resolve = current_snapshot.post_resolve_system(current_snapshot.event) + post_resolve = current_snapshot.post_resolve_system(current_snapshot.next_event) pre_evolve_next = next_snapshot.pre_evolve_system() assert post_resolve == pre_evolve_next @@ -87,7 +87,7 @@ def test_system_state_progression(): for step in range(len(seq)): snapshot = seq[step] - event = snapshot.event + event = snapshot.next_event pre_evolve = snapshot.pre_evolve_system() post_evolve = snapshot.post_evolve_system(event) diff --git a/tests/physics/resolve/ball_ball/test_ball_ball.py b/tests/physics/resolve/ball_ball/test_ball_ball.py index 6d8815df..c5778c97 100644 --- a/tests/physics/resolve/ball_ball/test_ball_ball.py +++ b/tests/physics/resolve/ball_ball/test_ball_ball.py @@ -1,4 +1,5 @@ import math +from pathlib import Path import attrs import numpy as np @@ -10,6 +11,10 @@ from pooltool.physics.resolve.ball_ball.frictional_inelastic import FrictionalInelastic from pooltool.physics.resolve.ball_ball.frictional_mathavan import FrictionalMathavan from pooltool.physics.resolve.ball_ball.frictionless_elastic import FrictionlessElastic +from pooltool.ptmath.utils import norm3d +from pooltool.serialize import conversion + +MAKE_KISS_TEST_DIR = Path(__file__).parent / "make_kiss_data" def vector_from_magnitude_and_direction(magnitude: float, angle_radians: float): @@ -197,7 +202,7 @@ def test_gearing_z_spin( assert np.allclose( np.cross(ob_f.vel, unit_normal), np.zeros_like(unit_normal), atol=1e-3 ), "Gearing english shouldn't cause throw" - assert abs(ob_f.avel[2]) < 5e-3, "Gearing english shouldn't cause induced side-spin" + assert abs(ob_f.avel[2]) < 1e-2, "Gearing english shouldn't cause induced side-spin" @pytest.mark.parametrize("model", [FrictionalInelastic()]) @@ -255,3 +260,37 @@ def test_low_relative_surface_velocity( assert ptmath.norm3d(cb_v_c_f - ob_v_c_f) < 1e-3, ( "Final relative contact velocity should be zero" ) + + +@pytest.mark.parametrize( + "model", + [ + FrictionalInelastic(), + FrictionalMathavan(), + FrictionlessElastic(), + ], +) +@pytest.mark.parametrize( + "case", + [ + "case_01", + ], +) +def test_make_kiss(model: BallBallCollisionStrategy, case: str): + ball1 = conversion.structure_from( + MAKE_KISS_TEST_DIR / f"{case}_ball_1.msgpack", Ball + ) + ball2 = conversion.structure_from( + MAKE_KISS_TEST_DIR / f"{case}_ball_2.msgpack", Ball + ) + + spacer = norm3d(ball1.xyz - ball2.xyz) - ball1.params.R * 2 + print(spacer) + print(ball1.xyz) + print(ball2.xyz) + + model.make_kiss(ball1, ball2) + spacer = norm3d(ball1.xyz - ball2.xyz) - ball1.params.R * 2 + print(spacer) + print(ball1.xyz) + print(ball2.xyz) From d7101367cc15cbb6d0fc3466c9cff5603f95d99b Mon Sep 17 00:00:00 2001 From: Evan Kiefl Date: Mon, 26 Jan 2026 16:07:33 -0700 Subject: [PATCH 10/22] Remove half-baked tests --- .../resolve/ball_ball/test_ball_ball.py | 39 ------------------- 1 file changed, 39 deletions(-) diff --git a/tests/physics/resolve/ball_ball/test_ball_ball.py b/tests/physics/resolve/ball_ball/test_ball_ball.py index c5778c97..a024eadb 100644 --- a/tests/physics/resolve/ball_ball/test_ball_ball.py +++ b/tests/physics/resolve/ball_ball/test_ball_ball.py @@ -1,5 +1,4 @@ import math -from pathlib import Path import attrs import numpy as np @@ -11,10 +10,6 @@ from pooltool.physics.resolve.ball_ball.frictional_inelastic import FrictionalInelastic from pooltool.physics.resolve.ball_ball.frictional_mathavan import FrictionalMathavan from pooltool.physics.resolve.ball_ball.frictionless_elastic import FrictionlessElastic -from pooltool.ptmath.utils import norm3d -from pooltool.serialize import conversion - -MAKE_KISS_TEST_DIR = Path(__file__).parent / "make_kiss_data" def vector_from_magnitude_and_direction(magnitude: float, angle_radians: float): @@ -260,37 +255,3 @@ def test_low_relative_surface_velocity( assert ptmath.norm3d(cb_v_c_f - ob_v_c_f) < 1e-3, ( "Final relative contact velocity should be zero" ) - - -@pytest.mark.parametrize( - "model", - [ - FrictionalInelastic(), - FrictionalMathavan(), - FrictionlessElastic(), - ], -) -@pytest.mark.parametrize( - "case", - [ - "case_01", - ], -) -def test_make_kiss(model: BallBallCollisionStrategy, case: str): - ball1 = conversion.structure_from( - MAKE_KISS_TEST_DIR / f"{case}_ball_1.msgpack", Ball - ) - ball2 = conversion.structure_from( - MAKE_KISS_TEST_DIR / f"{case}_ball_2.msgpack", Ball - ) - - spacer = norm3d(ball1.xyz - ball2.xyz) - ball1.params.R * 2 - print(spacer) - print(ball1.xyz) - print(ball2.xyz) - - model.make_kiss(ball1, ball2) - spacer = norm3d(ball1.xyz - ball2.xyz) - ball1.params.R * 2 - print(spacer) - print(ball1.xyz) - print(ball2.xyz) From ef866c451d94d94cfc16f9b889beaee72b0f1be5 Mon Sep 17 00:00:00 2001 From: Evan Kiefl Date: Mon, 26 Jan 2026 16:08:41 -0700 Subject: [PATCH 11/22] Linting + doc touchups --- pooltool/evolution/event_based/simulate.py | 5 ++++- pooltool/physics/resolve/ball_cushion/core.py | 4 ++-- tests/evolution/event_based/test_simulate.py | 5 ++--- 3 files changed, 8 insertions(+), 6 deletions(-) diff --git a/pooltool/evolution/event_based/simulate.py b/pooltool/evolution/event_based/simulate.py index 7b956064..a1f81c81 100755 --- a/pooltool/evolution/event_based/simulate.py +++ b/pooltool/evolution/event_based/simulate.py @@ -354,7 +354,10 @@ def get_next_ball_ball_collision( ): cache[ball_pair] = np.inf elif ptmath.is_overlapping( - ball1_state.rvw, ball2_state.rvw, ball1_params.R, ball2_params.R, + ball1_state.rvw, + ball2_state.rvw, + ball1_params.R, + ball2_params.R, ): cache[ball_pair] = shot.t else: diff --git a/pooltool/physics/resolve/ball_cushion/core.py b/pooltool/physics/resolve/ball_cushion/core.py index 6c456763..cd7ad45c 100644 --- a/pooltool/physics/resolve/ball_cushion/core.py +++ b/pooltool/physics/resolve/ball_cushion/core.py @@ -56,7 +56,7 @@ def make_kiss(self, ball: Ball, cushion: LinearCushionSegment) -> Ball: This makes a correction such that if the ball is not a distance R from the cushion, the ball is moved along the normal such that it is. To avoid downstream float precision round-off error, a small epsilon of additional distance - (`spacer`) is put between them, ensuring the cushion and ball are + (``spacer``) is put between them, ensuring the cushion and ball are separated post-resolution. """ normal = cushion.get_normal_xy(ball.xyz) @@ -106,7 +106,7 @@ def make_kiss(self, ball: Ball, cushion: CircularCushionSegment) -> Ball: This makes a correction such that if the ball is not a distance R from the cushion, the ball is moved along the normal such that it is. To avoid downstream float precision round-off error, a small epsilon of additional distance - (`spacer`) is put between them, ensuring the cushion and ball are + (``spacer``) is put between them, ensuring the cushion and ball are separated post-resolution. """ normal = cushion.get_normal_xy(ball.xyz) diff --git a/tests/evolution/event_based/test_simulate.py b/tests/evolution/event_based/test_simulate.py index 5ed1a679..32cd6776 100644 --- a/tests/evolution/event_based/test_simulate.py +++ b/tests/evolution/event_based/test_simulate.py @@ -4,7 +4,6 @@ import pooltool.constants as const import pooltool.ptmath as ptmath -from pooltool.evolution.event_based.introspection import simulate_with_snapshots from pooltool import aim, events from pooltool.events import EventType, ball_ball_collision, ball_pocket_collision from pooltool.evolution.event_based.cache import CollisionCache @@ -407,8 +406,8 @@ def true_time_to_collision(eps, V0, mu_r, g): assert diff < 10e-12 # Less than 10 femptosecond difference -def test_no_ball_ball_collisions_for_intersecting_balls(): - """Two already intersecting balls don't collide +def test_ball_ball_collision_for_intersecting_balls(): + """Two already intersecting balls don't collide. Previously, intersecting balls were prevented from colliding to avoid perpetual internal collisions. Now, with the improved make_kiss implementation, intersecting From 626315d989598e4c7f613b8990bae301720f5557 Mon Sep 17 00:00:00 2001 From: Evan Kiefl Date: Mon, 26 Jan 2026 16:54:31 -0700 Subject: [PATCH 12/22] Add spacer --- pooltool/ptmath/utils.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/pooltool/ptmath/utils.py b/pooltool/ptmath/utils.py index a402c6ce..253288d5 100644 --- a/pooltool/ptmath/utils.py +++ b/pooltool/ptmath/utils.py @@ -409,9 +409,13 @@ def get_ball_energy(rvw: NDArray[np.float64], R: float, m: float) -> float: def is_overlapping( - rvw1: NDArray[np.float64], rvw2: NDArray[np.float64], R1: float, R2: float + rvw1: NDArray[np.float64], + rvw2: NDArray[np.float64], + R1: float, + R2: float, + min_spacer: float = 0.0, ) -> bool: - return norm3d(rvw1[0] - rvw2[0]) < (R1 + R2) + return norm3d(rvw1[0] - rvw2[0]) < (R1 + R2 + min_spacer) @jit(nopython=True, cache=const.use_numba_cache) From fd62486e5588d9286dbd86634541ac8488df0599 Mon Sep 17 00:00:00 2001 From: Evan Kiefl Date: Mon, 26 Jan 2026 17:18:53 -0700 Subject: [PATCH 13/22] Refactor the fallback (midpoint shift) --- pooltool/constants.py | 4 ++ pooltool/physics/resolve/ball_ball/core.py | 53 ++++++---------------- 2 files changed, 18 insertions(+), 39 deletions(-) diff --git a/pooltool/constants.py b/pooltool/constants.py index 1ce73b30..bbc3d852 100644 --- a/pooltool/constants.py +++ b/pooltool/constants.py @@ -13,12 +13,16 @@ EPS = np.finfo(float).eps * 100 +MIN_DIST = 1e-6 +"""The minimum distance between balls.""" + # Ball states stationary: int = 0 """The stationary motion state label A ball with this motion state is both motionless and not in a pocket. """ + spinning: int = 1 """The spinning motion state label diff --git a/pooltool/physics/resolve/ball_ball/core.py b/pooltool/physics/resolve/ball_ball/core.py index 8290d645..fc891d53 100644 --- a/pooltool/physics/resolve/ball_ball/core.py +++ b/pooltool/physics/resolve/ball_ball/core.py @@ -3,6 +3,7 @@ import numpy as np +import pooltool.constants as const import pooltool.ptmath as ptmath from pooltool.objects.ball.datatypes import Ball @@ -39,19 +40,17 @@ def make_kiss(self, ball1: Ball, ball2: Ball) -> tuple[Ball, Ball]: trajectories (position + velocity * time) to this configuration. Acceleration terms are assumed negligible. - If the required movement exceeds 10x the spacer, which can occur if balls are - moving with nearly the same velocity, a naive fallback strategy is used that - moves the balls uniformly along the line of centers until they're separated by - an amount ``spacer``. + If the midpoint (collision point) shifts by more than 5x the spacer, which can + occur if balls are moving with nearly the same velocity, a naive fallback + strategy is used that moves the balls uniformly along the line of centers until + they're separated by an amount ``spacer``. Algorithm: 1. Calculate quadratic coefficients for separation equation 2. Solve for time offset that achieves target separation 3. Move balls to corrected positions: r_new = r + t * v - 4. If movement > 10x spacer (fallback): - - Identify which ball is "chasing" (higher line-of-centers velocity) and - which is being "chased" (lower line-of-centers velocity) - - Apply full displacement to chased ball, other stays fixed + 4. If midpoint shifts more than 5x spacer (fallback): + - Balls moved along line of centers until separated by amount ``spacer``. Returns: tuple[Ball, Ball]: @@ -62,14 +61,7 @@ def make_kiss(self, ball1: Ball, ball2: Ball) -> tuple[Ball, Ball]: v1 = ball1.state.rvw[1] v2 = ball2.state.rvw[1] - approach_speed = np.dot(v1 - v2, ptmath.unit_vector(r2 - r1)) - - spacer_min = 0.5e-6 - spacer_max = 1e-6 - speed_factor = max(0, min(approach_speed, 10)) - spacer = spacer_min + speed_factor * (spacer_max - spacer_min) - initial_distance = ptmath.norm3d(r2 - r1) - target_distance = 2 * ball1.params.R + spacer + spacer = const.MIN_DIST Bx = v2[0] - v1[0] By = v2[1] - v1[1] @@ -100,29 +92,12 @@ def make_kiss(self, ball1: Ball, ball2: Ball) -> tuple[Ball, Ball]: # (i.e. nearly same velocity). In these cases, the amount of time to create a # distance `spacer` between them can be high enough to significantly displace # both balls. - for r, r_corrected in zip([r1, r2], [r1_corrected, r2_corrected]): - if ptmath.norm3d(r - r_corrected) > 2 * spacer: - line_of_centers = (r2 - r1) / initial_distance - total_displacement = target_distance - initial_distance - - ## Determine which ball is chasing (has higher radial velocity). - # v1_radial = np.dot(v1, line_of_centers) - # v2_radial = np.dot(v2, line_of_centers) - - # if v1_radial > v2_radial: - # # Ball 2 is chasing ball 1, pull 2 backwards. - # r1_corrected = r1 - # r2_corrected = r2 + total_displacement * line_of_centers - # else: - # # Ball 1 is chasing ball 2, pull 1 backwards. - # r1_corrected = r1 - total_displacement * line_of_centers - # r2_corrected = r2 - - correction = 2 * ball1.params.R - ptmath.norm3d(r2 - r1) + spacer - r1_corrected = r1 - correction / 2 * ptmath.unit_vector(r2 - r1) - r2_corrected = r2 + correction / 2 * ptmath.unit_vector(r2 - r1) - - break + midpoint = (r1 + r2) / 2 + midpoint_corrected = (r1_corrected + r2_corrected) / 2 + if ptmath.norm3d(midpoint - midpoint_corrected) > 5 * spacer: + correction = 2 * ball1.params.R - ptmath.norm3d(r2 - r1) + spacer + r1_corrected = r1 - correction / 2 * ptmath.unit_vector(r2 - r1) + r2_corrected = r2 + correction / 2 * ptmath.unit_vector(r2 - r1) ball1.state.rvw[0] = r1_corrected ball2.state.rvw[0] = r2_corrected From 15a02334e9309d2065d94b94b4be95a7aa7b2a7f Mon Sep 17 00:00:00 2001 From: Evan Kiefl Date: Mon, 26 Jan 2026 17:20:47 -0700 Subject: [PATCH 14/22] Remove continually touching fn --- pooltool/physics/resolve/ball_ball/core.py | 100 --------------------- 1 file changed, 100 deletions(-) diff --git a/pooltool/physics/resolve/ball_ball/core.py b/pooltool/physics/resolve/ball_ball/core.py index fc891d53..d9a3d13f 100644 --- a/pooltool/physics/resolve/ball_ball/core.py +++ b/pooltool/physics/resolve/ball_ball/core.py @@ -104,105 +104,6 @@ def make_kiss(self, ball1: Ball, ball2: Ball) -> tuple[Ball, Ball]: return ball1, ball2 - def resolve_continually_touching( - self, ball1: Ball, ball2: Ball - ) -> tuple[Ball, Ball]: - """Prevent repeated collision detection for nearly-touching balls moving in unison. - - This method is called to handle rare cases where balls are moving with very - similar velocities. This can happen in some edge cases when frozen balls in a - perfect line are given energy along their line, (e.g. Newton's cradle). Without - intervention, the balls repeatedly trigger events microseconds apart that stall - progression of the simulation, sometimes indefinitely, via an explosion of - events. - - This is an unfortunate consequence of modeling non-instantaneous multibody - collisions using instantaneous pairwise collisions, and resolving it requires - some phenomonelogical intervention that hopefully appears to be realistic, - despite it not being grounded in theory. - - The solution is applied in ths method, and uses a momentum transfer mechanism: - the "chased" ball (slower in the line-of-centers direction) steals a fraction of - the "chaser's" radial momentum. This creates gradual separation over time, - preventing the balls from triggering repeated collision events while maintaining - physically plausible behavior. - - Algorithm: - 1. Projects velocities onto line of centers to get radial components - 2. If radial relative velocity is below threshold (< 1mm/s): - - Identifies which ball is "chasing" (higher radial velocity) - - Chased ball steals 10% of chaser's radial momentum - - Chaser loses this momentum, chased gains it - 3. Tangential velocity components remain unchanged - - Args: - ball1: First ball in the collision - ball2: Second ball in the collision - - Returns: - Modified ball1 and ball2 with adjusted velocities - - Notes: - - Practically speaking, this is a no-op method for all but the most - contrived simulation conditions. For one such condition, see - `sandbox/newtons_cradle.py` - """ - r1 = ball1.state.rvw[0] - r2 = ball2.state.rvw[0] - v1 = ball1.state.rvw[1] - v2 = ball2.state.rvw[1] - - v1_speed = ptmath.norm3d(v1) - v2_speed = ptmath.norm3d(v2) - both_moving = v1_speed > 0 and v2_speed > 0 - - if not both_moving: - return ball1, ball2 - - theft_fraction = 0.1 - velocity_similarity_threshold = 0.9 - - line_of_centers = ptmath.unit_vector(r2 - r1) - - # Velocities projected onto the line of centers (loc). - v1_loc = np.dot(v1, line_of_centers) - v2_loc = np.dot(v2, line_of_centers) - - cosine_similarity = np.dot(v1, v2) / (v1_speed * v2_speed) - velocities_aligned = cosine_similarity > velocity_similarity_threshold - - if abs(v2_loc - v1_loc) < 0.01 and velocities_aligned: - if v1_loc > v2_loc: - chaser_loc_vel = v1_loc - ball1_is_chaser = True - else: - chaser_loc_vel = v2_loc - ball1_is_chaser = False - - # Chased ball steals fraction of chaser's line of centers momentum - # FIXME: We assume equal mass, so transfer velocity directly - theft_fraction = theft_fraction - stolen_loc_velocity = chaser_loc_vel * theft_fraction - - if ball1_is_chaser: - v1_loc_new = v1_loc - stolen_loc_velocity - v2_loc_new = v2_loc + stolen_loc_velocity - else: - v1_loc_new = v1_loc + stolen_loc_velocity - v2_loc_new = v2_loc - stolen_loc_velocity - - v1_corrected = v1 - v1_loc * line_of_centers + v1_loc_new * line_of_centers - v2_corrected = v2 - v2_loc * line_of_centers + v2_loc_new * line_of_centers - - momentum_before = v1 + v2 - momentum_after = v1_corrected + v2_corrected - assert np.allclose(momentum_before, momentum_after, rtol=1e-10) - - ball1.state.rvw[1] = v1_corrected - ball2.state.rvw[1] = v2_corrected - - return ball1, ball2 - def resolve( self, ball1: Ball, ball2: Ball, inplace: bool = False ) -> tuple[Ball, Ball]: @@ -212,7 +113,6 @@ def resolve( ball1, ball2 = self.make_kiss(ball1, ball2) ball1, ball2 = self.solve(ball1, ball2) - # ball1, ball2 = self.resolve_continually_touching(ball1, ball2) return ball1, ball2 From ead73740569521d7a8596739aeb1df81a394ce0e Mon Sep 17 00:00:00 2001 From: Evan Kiefl Date: Mon, 26 Jan 2026 17:33:54 -0700 Subject: [PATCH 15/22] Prioritize event between overlapping balls -- breaks test --- pooltool/evolution/event_based/simulate.py | 10 +- pooltool/physics/resolve/ball_ball/core.py | 111 +++++++++++++-------- 2 files changed, 72 insertions(+), 49 deletions(-) diff --git a/pooltool/evolution/event_based/simulate.py b/pooltool/evolution/event_based/simulate.py index a1f81c81..2fdd1a8c 100755 --- a/pooltool/evolution/event_based/simulate.py +++ b/pooltool/evolution/event_based/simulate.py @@ -348,11 +348,6 @@ def get_next_ball_ball_collision( if ball1_state.s == const.pocketed or ball2_state.s == const.pocketed: cache[ball_pair] = np.inf - elif ( - ball1_state.s in const.nontranslating - and ball2_state.s in const.nontranslating - ): - cache[ball_pair] = np.inf elif ptmath.is_overlapping( ball1_state.rvw, ball2_state.rvw, @@ -360,6 +355,11 @@ def get_next_ball_ball_collision( ball2_params.R, ): cache[ball_pair] = shot.t + elif ( + ball1_state.s in const.nontranslating + and ball2_state.s in const.nontranslating + ): + cache[ball_pair] = np.inf else: dtau_E = solve.ball_ball_collision_time( rvw1=ball1_state.rvw, diff --git a/pooltool/physics/resolve/ball_ball/core.py b/pooltool/physics/resolve/ball_ball/core.py index d9a3d13f..60a6f560 100644 --- a/pooltool/physics/resolve/ball_ball/core.py +++ b/pooltool/physics/resolve/ball_ball/core.py @@ -27,6 +27,24 @@ def solve(self, ball1: Ball, ball2: Ball) -> tuple[Ball, Ball]: class CoreBallBallCollision(ABC): """Operations used by every ball-ball collision resolver""" + def _apply_fallback_positioning( + self, + ball1: Ball, + ball2: Ball, + r1: np.ndarray, + r2: np.ndarray, + spacer: float, + ) -> tuple[np.ndarray, np.ndarray]: + """Apply fallback positioning by moving balls along line of centers. + + This fallback strategy moves balls uniformly along the line of centers until + they're separated by the target distance (2*R + spacer). + """ + correction = 2 * ball1.params.R - ptmath.norm3d(r2 - r1) + spacer + r1_corrected = r1 - correction / 2 * ptmath.unit_vector(r2 - r1) + r2_corrected = r2 + correction / 2 * ptmath.unit_vector(r2 - r1) + return r1_corrected, r2_corrected + def make_kiss(self, ball1: Ball, ball2: Ball) -> tuple[Ball, Ball]: """Position balls at precise target separation before collision resolution. @@ -40,17 +58,18 @@ def make_kiss(self, ball1: Ball, ball2: Ball) -> tuple[Ball, Ball]: trajectories (position + velocity * time) to this configuration. Acceleration terms are assumed negligible. - If the midpoint (collision point) shifts by more than 5x the spacer, which can - occur if balls are moving with nearly the same velocity, a naive fallback - strategy is used that moves the balls uniformly along the line of centers until - they're separated by an amount ``spacer``. + If both balls are non-translating, or if the midpoint (collision point) shifts + by more than 5x the spacer (which can occur if balls are moving with nearly the + same velocity), a naive fallback strategy is used that moves the balls + uniformly along the line of centers until they're separated by an amount + ``spacer``. Algorithm: - 1. Calculate quadratic coefficients for separation equation - 2. Solve for time offset that achieves target separation - 3. Move balls to corrected positions: r_new = r + t * v - 4. If midpoint shifts more than 5x spacer (fallback): - - Balls moved along line of centers until separated by amount ``spacer``. + 1. If both balls are non-translating, apply fallback + 2. Otherwise, calculate quadratic coefficients for separation equation + 3. Solve for time offset that achieves target separation + 4. Move balls to corrected positions: r_new = r + t * v + 5. If midpoint shifts more than 5x spacer, apply fallback Returns: tuple[Ball, Ball]: @@ -63,41 +82,45 @@ def make_kiss(self, ball1: Ball, ball2: Ball) -> tuple[Ball, Ball]: spacer = const.MIN_DIST - Bx = v2[0] - v1[0] - By = v2[1] - v1[1] - Bz = v2[2] - v1[2] - Cx = r2[0] - r1[0] - Cy = r2[1] - r1[1] - Cz = r2[2] - r1[2] - alpha = Bx * Bx + By * By + Bz * Bz - beta = 2 * Bx * Cx + 2 * By * Cy + 2 * Bz * Cz - gamma = ( - Cx * Cx - + Cy * Cy - + Cz * Cz - - (2 * ball1.params.R + spacer) * (2 * ball1.params.R + spacer) - ) - roots_complex = ptmath.roots.quadratic.solve_complex(alpha, beta, gamma) - - imag_mag = np.abs(roots_complex.imag) - real_mag = np.abs(roots_complex.real) - keep = (imag_mag / real_mag) < 1e-3 - roots = roots_complex[keep].real - t = roots[np.abs(roots).argmin()] - - r1_corrected = r1 + t * v1 - r2_corrected = r2 + t * v2 - - # This fallback exists for cases when the balls are moving nearly in unison - # (i.e. nearly same velocity). In these cases, the amount of time to create a - # distance `spacer` between them can be high enough to significantly displace - # both balls. - midpoint = (r1 + r2) / 2 - midpoint_corrected = (r1_corrected + r2_corrected) / 2 - if ptmath.norm3d(midpoint - midpoint_corrected) > 5 * spacer: - correction = 2 * ball1.params.R - ptmath.norm3d(r2 - r1) + spacer - r1_corrected = r1 - correction / 2 * ptmath.unit_vector(r2 - r1) - r2_corrected = r2 + correction / 2 * ptmath.unit_vector(r2 - r1) + if ( + ball1.state.s in const.nontranslating + and ball2.state.s in const.nontranslating + ): + r1_corrected, r2_corrected = self._apply_fallback_positioning( + ball1, ball2, r1, r2, spacer + ) + else: + Bx = v2[0] - v1[0] + By = v2[1] - v1[1] + Bz = v2[2] - v1[2] + Cx = r2[0] - r1[0] + Cy = r2[1] - r1[1] + Cz = r2[2] - r1[2] + alpha = Bx * Bx + By * By + Bz * Bz + beta = 2 * Bx * Cx + 2 * By * Cy + 2 * Bz * Cz + gamma = ( + Cx * Cx + + Cy * Cy + + Cz * Cz + - (2 * ball1.params.R + spacer) * (2 * ball1.params.R + spacer) + ) + roots_complex = ptmath.roots.quadratic.solve_complex(alpha, beta, gamma) + + imag_mag = np.abs(roots_complex.imag) + real_mag = np.abs(roots_complex.real) + keep = (imag_mag / real_mag) < 1e-3 + roots = roots_complex[keep].real + t = roots[np.abs(roots).argmin()] + + r1_corrected = r1 + t * v1 + r2_corrected = r2 + t * v2 + + midpoint = (r1 + r2) / 2 + midpoint_corrected = (r1_corrected + r2_corrected) / 2 + if ptmath.norm3d(midpoint - midpoint_corrected) > 5 * spacer: + r1_corrected, r2_corrected = self._apply_fallback_positioning( + ball1, ball2, r1, r2, spacer + ) ball1.state.rvw[0] = r1_corrected ball2.state.rvw[0] = r2_corrected From 9c9234f8b3b5506bb7f8011980feabf07cbf473a Mon Sep 17 00:00:00 2001 From: Evan Kiefl Date: Mon, 26 Jan 2026 21:40:27 -0700 Subject: [PATCH 16/22] Resolve simultaneous events according to an event priority heuristic --- pooltool/evolution/event_based/simulate.py | 127 +++++++++++++++------ 1 file changed, 92 insertions(+), 35 deletions(-) diff --git a/pooltool/evolution/event_based/simulate.py b/pooltool/evolution/event_based/simulate.py index 2fdd1a8c..a20b9e63 100755 --- a/pooltool/evolution/event_based/simulate.py +++ b/pooltool/evolution/event_based/simulate.py @@ -50,6 +50,61 @@ def _system_has_energy(system: System) -> bool: ) +def get_event_priority(event: Event, shot: System) -> tuple[int, float]: + """Compute priority for an event to resolve ties among simultaneous events. + + Returns a tuple (tier, energy) where: + - Lower tier = higher priority + - Higher energy = higher priority within the same tier + + Priority tiers: + - Tier 1: STICK_BALL (always first) + - Tier 2: Transitions and BALL_POCKET (can resolve without affecting others) + - Tier 3: BALL_BALL and ball-cushion collisions + + Args: + event: The event to compute priority for. + shot: The system state at the time the event was detected. + + Returns: + A tuple of (tier, energy) for sorting. + """ + event_type = event.event_type + + if event_type == EventType.NONE: + return (99, 0.0) + + if event_type == EventType.STICK_BALL: + return (1, shot.cue.V0**2) + + if event_type == EventType.BALL_POCKET: + ball_id = event.ids[0] + ball = shot.balls[ball_id] + energy = ptmath.get_ball_energy(ball.state.rvw, ball.params.R, ball.params.m) + return (2, energy) + + if event_type.is_transition(): + ball_id = event.ids[0] + ball = shot.balls[ball_id] + energy = ptmath.get_ball_energy(ball.state.rvw, ball.params.R, ball.params.m) + return (2, energy) + + if event_type == EventType.BALL_BALL: + ball1_id, ball2_id = event.ids + v1 = shot.balls[ball1_id].state.rvw[1] + v2 = shot.balls[ball2_id].state.rvw[1] + energy = ptmath.squared_norm3d(v1 - v2) + return (3, energy) + + if event_type in (EventType.BALL_LINEAR_CUSHION, EventType.BALL_CIRCULAR_CUSHION): + ball_id = event.ids[0] + ball = shot.balls[ball_id] + energy = ptmath.get_ball_energy(ball.state.rvw, ball.params.R, ball.params.m) + return (3, energy) + + return (99, 0.0) + + @attrs.define class _SimulationState: shot: System @@ -253,8 +308,8 @@ def get_next_event( if collision_cache is None: collision_cache = CollisionCache.create() - # Start by assuming next event doesn't happen - event = null_event(time=np.inf) + # Collect all candidate events from each detection function. + candidates: list[Event] = [] # Stick-ball collisions only occur at t=0 (shot initiation), so we skip this # check after the first timestep as an optimization. Other collision types are @@ -262,41 +317,43 @@ def get_next_event( # at t=0, we still call the remaining detection functions to fully populate the # collision cache, which is needed by debug/introspection tools. if shot.t == 0: - stick_ball_event = get_next_stick_ball_collision( - shot, collision_cache=collision_cache + candidates.append( + get_next_stick_ball_collision(shot, collision_cache=collision_cache) ) - if stick_ball_event.time < event.time: - event = stick_ball_event - transition_event = transition_cache.get_next() - if transition_event.time < event.time: - event = transition_event - - ball_ball_event = get_next_ball_ball_collision( - shot, collision_cache=collision_cache + candidates.append(transition_cache.get_next()) + candidates.append( + get_next_ball_ball_collision(shot, collision_cache=collision_cache) ) - if ball_ball_event.time < event.time: - event = ball_ball_event - - ball_circular_cushion_event = get_next_ball_circular_cushion_event( - shot, collision_cache=collision_cache + candidates.append( + get_next_ball_circular_cushion_event(shot, collision_cache=collision_cache) ) - if ball_circular_cushion_event.time < event.time: - event = ball_circular_cushion_event - - ball_linear_cushion_event = get_next_ball_linear_cushion_collision( - shot, collision_cache=collision_cache + candidates.append( + get_next_ball_linear_cushion_collision(shot, collision_cache=collision_cache) ) - if ball_linear_cushion_event.time < event.time: - event = ball_linear_cushion_event - - ball_pocket_event = get_next_ball_pocket_collision( - shot, collision_cache=collision_cache + candidates.append( + get_next_ball_pocket_collision(shot, collision_cache=collision_cache) ) - if ball_pocket_event.time < event.time: - event = ball_pocket_event - return event + # Find the earliest time among all candidates. + min_time = min(event.time for event in candidates) + + if min_time == np.inf: + return null_event(time=np.inf) + + # Filter to only events occurring at the earliest time. + simultaneous = [e for e in candidates if e.time == min_time] + + if len(simultaneous) == 1: + return simultaneous[0] + + # When multiple events occur at the same time, select by priority tier, then by + # energy within the tier (higher energy first). + def sort_key(e: Event) -> tuple[int, float]: + tier, energy = get_event_priority(e, shot) + return (tier, -energy) + + return min(simultaneous, key=sort_key) def get_next_stick_ball_collision( @@ -348,6 +405,11 @@ def get_next_ball_ball_collision( if ball1_state.s == const.pocketed or ball2_state.s == const.pocketed: cache[ball_pair] = np.inf + elif ( + ball1_state.s in const.nontranslating + and ball2_state.s in const.nontranslating + ): + cache[ball_pair] = np.inf elif ptmath.is_overlapping( ball1_state.rvw, ball2_state.rvw, @@ -355,11 +417,6 @@ def get_next_ball_ball_collision( ball2_params.R, ): cache[ball_pair] = shot.t - elif ( - ball1_state.s in const.nontranslating - and ball2_state.s in const.nontranslating - ): - cache[ball_pair] = np.inf else: dtau_E = solve.ball_ball_collision_time( rvw1=ball1_state.rvw, From 2ae543bab9690fbd01f3e5c630631c567a5ea3fe Mon Sep 17 00:00:00 2001 From: Evan Kiefl Date: Tue, 27 Jan 2026 13:42:45 -0800 Subject: [PATCH 17/22] Add safety factor (saves newtons cradle, at least partially) --- pooltool/ptmath/roots/_quartic_numba.py | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/pooltool/ptmath/roots/_quartic_numba.py b/pooltool/ptmath/roots/_quartic_numba.py index 52c28fcb..3bf60d6c 100644 --- a/pooltool/ptmath/roots/_quartic_numba.py +++ b/pooltool/ptmath/roots/_quartic_numba.py @@ -1,8 +1,14 @@ -"""1:1 exact translation of the "1010" quartic root-finding algorithm. +"""Translation of the "1010" quartic root-finding algorithm with modifications. The original implementation is written in C, and this module was written by Claude Code. -The 1:1 correspondence has been tested to floating point precision on a test of 100,000 -difficult to determine quartics. + +Modifications from the original algorithm: + - The threshold for falling back to an alternative factorization method has been + loosened by a safety factor (d2_safety_factor). The original threshold was too + strict for certain edge cases where the discriminant d2 is very small but nonzero, + causing the algorithm to produce duplicate roots instead of four distinct roots. + This is particularly relevant for Newton's cradle-like simulations where balls + move with very similar velocities. Solve speed: @@ -58,6 +64,7 @@ cubic_rescal_fact = 3.488062113727083e102 quart_rescal_fact = 7.156344627944542e76 macheps = 2.2204460492503131e-16 +d2_safety_factor = 100 @jit(nopython=True, cache=const.use_numba_cache) @@ -630,7 +637,8 @@ def solve(a: float, b: float, c: float, d: float, e: float) -> NDArray[np.comple realcase_0 = -1 if realcase_0 == -1 or ( - abs(d2) <= macheps * max(abs(2.0 * b_p / 3.0), abs(phi0), l1 * l1) + abs(d2) + <= d2_safety_factor * macheps * max(abs(2.0 * b_p / 3.0), abs(phi0), l1 * l1) ): d3 = d_p - l3 * l3 if realcase_0 == 1: From bd821996b150740d433f7925925823e54fd1ab59 Mon Sep 17 00:00:00 2001 From: Evan Kiefl Date: Tue, 27 Jan 2026 13:46:28 -0800 Subject: [PATCH 18/22] Add continually touching back - Saves newtons cradle --- pooltool/physics/resolve/ball_ball/core.py | 100 +++++++++++++++++++++ 1 file changed, 100 insertions(+) diff --git a/pooltool/physics/resolve/ball_ball/core.py b/pooltool/physics/resolve/ball_ball/core.py index 60a6f560..f40d2145 100644 --- a/pooltool/physics/resolve/ball_ball/core.py +++ b/pooltool/physics/resolve/ball_ball/core.py @@ -127,6 +127,105 @@ def make_kiss(self, ball1: Ball, ball2: Ball) -> tuple[Ball, Ball]: return ball1, ball2 + def resolve_continually_touching( + self, ball1: Ball, ball2: Ball + ) -> tuple[Ball, Ball]: + """Prevent repeated collision detection for nearly-touching balls moving in unison. + + This method is called to handle rare cases where balls are moving with very + similar velocities. This can happen in some edge cases when frozen balls in a + perfect line are given energy along their line, (e.g. Newton's cradle). Without + intervention, the balls repeatedly trigger events microseconds apart that stall + progression of the simulation, sometimes indefinitely, via an explosion of + events. + + This is an unfortunate consequence of modeling non-instantaneous multibody + collisions using instantaneous pairwise collisions, and resolving it requires + some phenomonelogical intervention that hopefully appears to be realistic, + despite it not being grounded in theory. + + The solution is applied in this method, and uses a momentum transfer mechanism: + the "chased" ball (slower in the line-of-centers direction) steals a fraction of + the "chaser's" radial momentum. This creates gradual separation over time, + preventing the balls from triggering repeated collision events while maintaining + physically plausible behavior. + + Algorithm: + 1. Projects velocities onto line of centers to get radial components + 2. If radial relative velocity is below threshold (< 1mm/s): + - Identifies which ball is "chasing" (higher radial velocity) + - Chased ball steals 10% of chaser's radial momentum + - Chaser loses this momentum, chased gains it + 3. Tangential velocity components remain unchanged + + Args: + ball1: First ball in the collision + ball2: Second ball in the collision + + Returns: + Modified ball1 and ball2 with adjusted velocities + + Notes: + - Practically speaking, this is a no-op method for all but the most + contrived simulation conditions. For one such condition, see + `sandbox/newtons_cradle.py` + """ + r1 = ball1.state.rvw[0] + r2 = ball2.state.rvw[0] + v1 = ball1.state.rvw[1] + v2 = ball2.state.rvw[1] + + v1_speed = ptmath.norm3d(v1) + v2_speed = ptmath.norm3d(v2) + both_moving = v1_speed > 0 and v2_speed > 0 + + if not both_moving: + return ball1, ball2 + + theft_fraction = 0.10 + velocity_similarity_threshold = 0.9 + + line_of_centers = ptmath.unit_vector(r2 - r1) + + # Velocities projected onto the line of centers (loc). + v1_loc = np.dot(v1, line_of_centers) + v2_loc = np.dot(v2, line_of_centers) + + cosine_similarity = np.dot(v1, v2) / (v1_speed * v2_speed) + velocities_aligned = cosine_similarity > velocity_similarity_threshold + + if abs(v2_loc - v1_loc) < 0.01 and velocities_aligned: + if v1_loc > v2_loc: + chaser_loc_vel = v1_loc + ball1_is_chaser = True + else: + chaser_loc_vel = v2_loc + ball1_is_chaser = False + + # Chased ball steals fraction of chaser's line of centers momentum + # FIXME: We assume equal mass, so transfer velocity directly + theft_fraction = theft_fraction + stolen_loc_velocity = chaser_loc_vel * theft_fraction + + if ball1_is_chaser: + v1_loc_new = v1_loc - stolen_loc_velocity + v2_loc_new = v2_loc + stolen_loc_velocity + else: + v1_loc_new = v1_loc + stolen_loc_velocity + v2_loc_new = v2_loc - stolen_loc_velocity + + v1_corrected = v1 - v1_loc * line_of_centers + v1_loc_new * line_of_centers + v2_corrected = v2 - v2_loc * line_of_centers + v2_loc_new * line_of_centers + + momentum_before = v1 + v2 + momentum_after = v1_corrected + v2_corrected + assert np.allclose(momentum_before, momentum_after, rtol=1e-10) + + ball1.state.rvw[1] = v1_corrected + ball2.state.rvw[1] = v2_corrected + + return ball1, ball2 + def resolve( self, ball1: Ball, ball2: Ball, inplace: bool = False ) -> tuple[Ball, Ball]: @@ -136,6 +235,7 @@ def resolve( ball1, ball2 = self.make_kiss(ball1, ball2) ball1, ball2 = self.solve(ball1, ball2) + ball1, ball2 = self.resolve_continually_touching(ball1, ball2) return ball1, ball2 From e12471dcd6650c97cf2a050f9c69498d8aa8a4f7 Mon Sep 17 00:00:00 2001 From: Evan Kiefl Date: Tue, 27 Jan 2026 14:10:35 -0800 Subject: [PATCH 19/22] Update layouts to support continuous jumps --- pooltool/layouts.py | 48 +++++++++++++++++++++++++++++---------------- 1 file changed, 31 insertions(+), 17 deletions(-) diff --git a/pooltool/layouts.py b/pooltool/layouts.py index 55466deb..f08a0c70 100755 --- a/pooltool/layouts.py +++ b/pooltool/layouts.py @@ -52,48 +52,59 @@ def translation_map(cls) -> dict[Dir, tuple[float, float]]: } +Translation = Dir | float + + class Jump: @staticmethod - def LEFT(quantity: int = 1) -> list[Dir]: + def LEFT(quantity: int = 1) -> list[Translation]: return [Dir.LEFT] * quantity @staticmethod - def RIGHT(quantity: int = 1) -> list[Dir]: + def RIGHT(quantity: int = 1) -> list[Translation]: return [Dir.RIGHT] * quantity @staticmethod - def UP(quantity: int = 1) -> list[Dir]: + def UP(quantity: int = 1) -> list[Translation]: return [Dir.UP] * quantity @staticmethod - def DOWN(quantity: int = 1) -> list[Dir]: + def DOWN(quantity: int = 1) -> list[Translation]: return [Dir.DOWN] * quantity @staticmethod - def UPLEFT(quantity: int = 1) -> list[Dir]: + def UPLEFT(quantity: int = 1) -> list[Translation]: return [Dir.UPLEFT] * quantity @staticmethod - def UPRIGHT(quantity: int = 1) -> list[Dir]: + def UPRIGHT(quantity: int = 1) -> list[Translation]: return [Dir.UPRIGHT] * quantity @staticmethod - def DOWNRIGHT(quantity: int = 1) -> list[Dir]: + def DOWNRIGHT(quantity: int = 1) -> list[Translation]: return [Dir.DOWNRIGHT] * quantity @staticmethod - def DOWNLEFT(quantity: int = 1) -> list[Dir]: + def DOWNLEFT(quantity: int = 1) -> list[Translation]: return [Dir.DOWNLEFT] * quantity @staticmethod - def eval(translations: list[Dir], radius: float) -> tuple[float, float]: + def ANGLE(degrees: float, quantity: int = 1) -> list[Translation]: + radians = np.radians(degrees) + return [radians] * quantity + + @staticmethod + def eval(translations: list[Translation], radius: float) -> tuple[float, float]: mapping = Dir.translation_map assert isinstance(mapping, dict) - dx, dy = 0, 0 + dx, dy = 0.0, 0.0 - for direction in translations: - i, j = mapping[direction] + for translation in translations: + if isinstance(translation, Dir): + i, j = mapping[translation] + else: + i, j = 2 * np.cos(translation), 2 * np.sin(translation) dx += i * radius dy += j * radius @@ -106,7 +117,9 @@ class Pos: Attributes: loc: - A sequence of translations. + A sequence of translations. Each translation is either a Dir enum + for discrete directions, or a float representing an angle in radians + (0 = right, pi/2 = up). Use Jump.ANGLE(degrees) for convenience. relative_to: This defines what the translation is with respect to. This can either be another Pos, or a 2D coordinate, normalized by the table's @@ -114,7 +127,7 @@ class Pos: so (0.0, 0.0) is bottom-left and (1.0, 1.0) is top right. """ - loc: list[Dir] + loc: list[Translation] relative_to: Pos | tuple[float, float] @@ -130,7 +143,7 @@ class BallPos(Pos): ids: set[str] -JumpSequence = list[tuple[list[Dir], set[str]]] +JumpSequence = list[tuple[list[Translation], set[str]]] def ball_cluster_blueprint(seed: BallPos, jump_sequence: JumpSequence) -> list[BallPos]: @@ -153,10 +166,10 @@ def _get_ball_ids(positions: list[BallPos]) -> set[str]: return ids -def _get_anchor_translation(pos: Pos) -> tuple[tuple[float, float], list[Dir]]: +def _get_anchor_translation(pos: Pos) -> tuple[tuple[float, float], list[Translation]]: """Traverse the position's parent hierarchy until the anchor is found""" - translation_from_anchor: list[Dir] = [] + translation_from_anchor: list[Translation] = [] translation_from_anchor.extend(pos.loc) parent = pos.relative_to @@ -543,6 +556,7 @@ def get_rack( "Jump", "Pos", "BallPos", + "Translation", "ball_cluster_blueprint", "generate_layout", "get_rack", From 1a310a013d782ff378bc1ec9e21eddb4f4d9a17e Mon Sep 17 00:00:00 2001 From: Evan Kiefl Date: Tue, 27 Jan 2026 14:12:10 -0800 Subject: [PATCH 20/22] argparsify newtons_cradle.py --- sandbox/newtons_cradle.py | 78 +++++++++++++++++++++++++++------------ 1 file changed, 55 insertions(+), 23 deletions(-) diff --git a/sandbox/newtons_cradle.py b/sandbox/newtons_cradle.py index 719b853e..6f98cb6f 100644 --- a/sandbox/newtons_cradle.py +++ b/sandbox/newtons_cradle.py @@ -1,32 +1,64 @@ +import argparse + +import numpy as np + import pooltool as pt -from pooltool.evolution.event_based.introspection import simulate_with_snapshots +from pooltool.layouts import BallPos, Jump, ball_cluster_blueprint, generate_layout -table = pt.Table.default() -ball_radius = pt.BallParams.default().R -balls = [ - pt.Ball.create( - str(i + 1), - xy=(0.5 * table.w, 0.4 * table.l + 2 * i * ball_radius), - ballset=pt.objects.BallSet("pooltool_pocket"), +def main(): + parser = argparse.ArgumentParser( + description="Simulate an imperfect Newton's cradle. WARNING: if angle-variation " + "is low and n-balls is high, the event count will skyrocket.", + formatter_class=argparse.ArgumentDefaultsHelpFormatter, + ) + parser.add_argument("--n-balls", type=int, default=4, help="Number of balls") + parser.add_argument( + "--angle-variation", type=float, default=0.3, help="Max angle offset (degrees)" ) - for i in range(3) -] + parser.add_argument("--seed", type=int, default=42, help="Random seed") + args = parser.parse_args() + + rng = np.random.default_rng(args.seed) + + table = pt.Table.default() + ball_params = pt.BallParams.default() + + jump_sequence = [ + ( + Jump.ANGLE(90 + rng.uniform(-args.angle_variation, args.angle_variation)), + {str(i + 2)}, + ) + for i in range(args.n_balls - 1) + ] -balls.append( - pt.Ball.create( - "cue", - xy=(0.5 * table.w, 0.1 * table.l), + blueprint = ball_cluster_blueprint( + seed=BallPos([], (0.5, 0.4), {"1"}), + jump_sequence=jump_sequence, + ) + + cue = BallPos([], (0.5, 0.1), {"cue"}) + blueprint.append(cue) + + balls = generate_layout( + blueprint, + table, ballset=pt.objects.BallSet("pooltool_pocket"), + ball_params=ball_params, + spacing_factor=0, ) -) -system = pt.System( - balls=balls, - table=table, - cue=pt.Cue.default(), -) + system = pt.System( + balls=balls, + table=table, + cue=pt.Cue.default(), + ) + + system.strike(V0=2, phi=pt.aim.at_ball(system, "1", cut=0)) + pt.simulate(system, inplace=True) + print(len(system.events)) + pt.show(system) + -engine = pt.physics.PhysicsEngine() -system.strike(V0=2, phi=pt.aim.at_ball(system, "1", cut=0)) -system, snapshots = simulate_with_snapshots(system) +if __name__ == "__main__": + main() From ebe412cebcc2a0471e33165d5a0d208492726222 Mon Sep 17 00:00:00 2001 From: Evan Kiefl Date: Tue, 27 Jan 2026 14:32:39 -0800 Subject: [PATCH 21/22] Revert to original tolerance --- tests/physics/resolve/ball_ball/test_ball_ball.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/physics/resolve/ball_ball/test_ball_ball.py b/tests/physics/resolve/ball_ball/test_ball_ball.py index a024eadb..6d8815df 100644 --- a/tests/physics/resolve/ball_ball/test_ball_ball.py +++ b/tests/physics/resolve/ball_ball/test_ball_ball.py @@ -197,7 +197,7 @@ def test_gearing_z_spin( assert np.allclose( np.cross(ob_f.vel, unit_normal), np.zeros_like(unit_normal), atol=1e-3 ), "Gearing english shouldn't cause throw" - assert abs(ob_f.avel[2]) < 1e-2, "Gearing english shouldn't cause induced side-spin" + assert abs(ob_f.avel[2]) < 5e-3, "Gearing english shouldn't cause induced side-spin" @pytest.mark.parametrize("model", [FrictionalInelastic()]) From 9828421f5492eaa73d4cafae10dd6f3ecfb760cb Mon Sep 17 00:00:00 2001 From: Evan Kiefl Date: Tue, 27 Jan 2026 15:19:08 -0800 Subject: [PATCH 22/22] Address rabbitai review --- pooltool/physics/resolve/ball_ball/core.py | 1 - pooltool/physics/resolve/transition/__init__.py | 10 +++++----- tests/evolution/event_based/test_simulate.py | 2 +- 3 files changed, 6 insertions(+), 7 deletions(-) diff --git a/pooltool/physics/resolve/ball_ball/core.py b/pooltool/physics/resolve/ball_ball/core.py index f40d2145..0a90ebad 100644 --- a/pooltool/physics/resolve/ball_ball/core.py +++ b/pooltool/physics/resolve/ball_ball/core.py @@ -204,7 +204,6 @@ def resolve_continually_touching( # Chased ball steals fraction of chaser's line of centers momentum # FIXME: We assume equal mass, so transfer velocity directly - theft_fraction = theft_fraction stolen_loc_velocity = chaser_loc_vel * theft_fraction if ball1_is_chaser: diff --git a/pooltool/physics/resolve/transition/__init__.py b/pooltool/physics/resolve/transition/__init__.py index 90845cae..5c5eee19 100644 --- a/pooltool/physics/resolve/transition/__init__.py +++ b/pooltool/physics/resolve/transition/__init__.py @@ -15,7 +15,7 @@ from pooltool.objects.ball.datatypes import Ball from pooltool.physics.resolve.models import BallTransitionModel -TOLERANCE = 1e-12 +_TOLERANCE = 1e-12 class BallTransitionStrategy(Protocol): @@ -49,8 +49,8 @@ def resolve(self, ball: Ball, transition: EventType, inplace: bool = False) -> B # angular velocity components are nearly 0. Then set them to exactly 0. v = ball.state.rvw[1] w = ball.state.rvw[2] - assert (np.abs(v) < TOLERANCE).all() - assert (np.abs(w[:2]) < TOLERANCE).all() + assert (np.abs(v) < _TOLERANCE).all() + assert (np.abs(w[:2]) < _TOLERANCE).all() ball.state.rvw[1, :] = [0.0, 0.0, 0.0] ball.state.rvw[2, :2] = [0.0, 0.0] @@ -60,8 +60,8 @@ def resolve(self, ball: Ball, transition: EventType, inplace: bool = False) -> B # set them to exactly 0. v = ball.state.rvw[1] w = ball.state.rvw[2] - assert (np.abs(v) < TOLERANCE).all() - assert (np.abs(w) < TOLERANCE).all() + assert (np.abs(v) < _TOLERANCE).all() + assert (np.abs(w) < _TOLERANCE).all() ball.state.rvw[1, :] = [0.0, 0.0, 0.0] ball.state.rvw[2, :] = [0.0, 0.0, 0.0] diff --git a/tests/evolution/event_based/test_simulate.py b/tests/evolution/event_based/test_simulate.py index 32cd6776..16aa3ebb 100644 --- a/tests/evolution/event_based/test_simulate.py +++ b/tests/evolution/event_based/test_simulate.py @@ -407,7 +407,7 @@ def true_time_to_collision(eps, V0, mu_r, g): def test_ball_ball_collision_for_intersecting_balls(): - """Two already intersecting balls don't collide. + """Two already intersecting balls collide. Previously, intersecting balls were prevented from colliding to avoid perpetual internal collisions. Now, with the improved make_kiss implementation, intersecting