Skip to content

Instantly share code, notes, and snippets.

@thquinn

thquinn/sagemath_juggling.py

Last active September 23, 2025 07:59
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save thquinn/39c140d572dbc8d67eed702f3a8af836 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save thquinn/39c140d572dbc8d67eed702f3a8af836 to your computer and use it in GitHub Desktop.
Download ZIP
A not-very-fast symbolic physics simulation using SymPy.
Raw
sagemath_juggling.py
import os
import sage.all as sage
from itertools import combinations
from PIL import Image, ImageDraw
from typing import cast, Dict
from enum import Enum
os.system("clear")
class SimCircle:
center: sage.vector
radius: sage.Expression
velocity: sage.vector
def __init__(self, center: sage.vector, radius: sage.Expression, velocity: sage.vector = sage.vector([0, 0])):
self.center = center
self.radius = radius
self.velocity = velocity
self.first = False
def to_tuple(self):
return sage.vector([self.center[0], self.center[1], self.velocity[0], self.velocity[1]])
def __str__(self):
return f'SimCircle center: {self.center}, radius: {self.radius}, velocity: {self.velocity}'
class SimSegment:
p1: sage.vector
p2: sage.vector
def __init__(self, p1: sage.vector, p2: sage.vector):
self.p1 = p1
self.p2 = p2
def length(self):
return sage.sqrt((self.p2[0] - self.p1[0]) ** 2 + (self.p2[1] - self.p1[1]) ** 2)
def contains_box(self, p):
xs = sorted([self.p1[0], self.p2[0]])
ys = sorted([self.p1[1], self.p2[1]])
return bool(p[0] >= xs[0]) and bool(p[0] <= xs[1]) and bool(p[1] >= ys[0]) and bool(p[1] <= ys[1])
def projection(self, p):
AB = self.p2 - self.p1
AC = p - self.p1
AD = AB * (AB.dot_product(AC)) / AB.dot_product(AB)
return self.p1 + AD
def __str__(self):
return f'SimSegment p1: {self.p1}, p2: {self.p2}'
class SimState(Enum):
INCOMPLETE = 1
COMPLETE_WITH_PERIOD = 2
INVALID_MULTIPLE_COLLISIONS = 3
INVALID_SPAWN_OVERLAP = 4
INVALID_TIMEOUT = 5
class Sim:
SPAWN_RATE: sage.Expression = sage.Integer(1)
TIMEOUT: sage.Expression = SPAWN_RATE * sage.Integer(50)
SPAWN_HEIGHT: sage.Expression = sage.Integer(3)
DESPAWN_HEIGHT: sage.Expression = sage.Integer(-10)
CIRCLE_RADIUS: sage.Expression = sage.Rational('1/3')
GRAVITY: sage.Expression = sage.Integer(-3)
SEGMENT_LENGTH: sage.Expression = sage.Rational('3/2')
SEGMENT_SEPARATION: sage.Expression = sage.Integer(3)
COEFFICIENT_OF_RESTITUTION: sage.Expression = sage.Rational('1/1')
FRAME_TIME = sage.Rational('1/10')
t: sage.Expression
circles: list[SimCircle]
segments: list[SimSegment]
next_spawn_t: sage.Expression
previous_configurations: Dict[list, sage.Expression]
state: SimState
def __init__(self, theta1:sage.Expression, theta2:sage.Expression):
self.t = sage.Integer(0)
self.circles = []
left_segment_center = sage.vector([0, 0])
right_segment_center = sage.vector([Sim.SEGMENT_SEPARATION, 0])
self.segments = [
SimSegment(left_segment_center - sage.vector([sage.cos(theta1), -sage.sin(theta1)]) * Sim.SEGMENT_LENGTH / 2, left_segment_center + sage.vector([sage.cos(theta1), -sage.sin(theta1)]) * Sim.SEGMENT_LENGTH / 2),
SimSegment(right_segment_center - sage.vector([sage.cos(theta2), sage.sin(theta2)]) * Sim.SEGMENT_LENGTH / 2, right_segment_center + sage.vector([sage.cos(theta2), sage.sin(theta2)]) * Sim.SEGMENT_LENGTH / 2)
]
self.next_spawn_t = sage.Integer(0)
self.previous_configurations = dict()
self.state = SimState.INCOMPLETE
def iterate(self):
collisions = []
for circle_pair in combinations(self.circles, 2):
next_collision = self.t + circle_collision_next_time(*circle_pair)
collisions.append([next_collision, *circle_pair])
for segment in self.segments:
for circle in self.circles:
next_collision = self.t + circle_segment_collision_next_time(circle, segment, Sim.GRAVITY)
collisions.append([next_collision, circle, segment])
min_time_collisions = []
min_time: sage.Expression = sage.infinity
for collision in collisions:
if collision[0] < min_time:
min_time_collisions = []
min_time = collision[0]
if collision[0] <= min_time:
min_time_collisions.append(collision)
min_time = min_time.full_simplify() if min_time != sage.oo else sage.oo
print(f'min collision time: {min_time.n()}, next spawn: {self.next_spawn_t.n()}')
if min_time >= self.next_spawn_t:
self.advance_t_to(self.next_spawn_t)
self.next_spawn_t += Sim.SPAWN_RATE
# TODO: check for circles in spawn area
self.circles.append(SimCircle(sage.vector([0, Sim.SPAWN_HEIGHT]), Sim.CIRCLE_RADIUS))
if (len(self.circles) == 1):
self.circles[0].first = True
return -1
# TODO: check for circles involved in multiple collisions
self.advance_t_to(min_time)
for collision in min_time_collisions:
c1 = cast(SimCircle, collision[1])
if isinstance(collision[2], SimSegment):
seg = cast(SimSegment, collision[2])
perform_circle_point_collision(c1, seg.projection(c1.center), Sim.COEFFICIENT_OF_RESTITUTION)
else:
c2 = cast(SimCircle, collision[2])
perform_circle_collision(c1, c2, Sim.COEFFICIENT_OF_RESTITUTION)
self.circles = [c for c in self.circles if c.center[1] > Sim.DESPAWN_HEIGHT]
print('circles in the sim:', len(self.circles))
configuration = [c.to_tuple() for c in self.circles]
configuration.sort()
configuration = str(configuration)
if configuration in self.previous_configurations:
prev_config_time = self.previous_configurations[configuration]
config_delta = self.t - prev_config_time
self.state = SimState.COMPLETE_WITH_PERIOD
return config_delta / Sim.SPAWN_RATE
self.previous_configurations[configuration] = self.t
if self.t > Sim.TIMEOUT:
self.state = SimState.INVALID_TIMEOUT
return -1
def advance_t_to(self, next_t: sage.Expression):
if (self.t == 0):
self.save_img();
print(f'advancing to t={next_t.n()}')
next_frame = sage.ceil(self.t / Sim.FRAME_TIME) * Sim.FRAME_TIME
while (next_frame < next_t):
self.advance_t_to_impl(next_frame)
next_frame += Sim.FRAME_TIME
self.save_img();
self.advance_t_to_impl(next_t)
def advance_t_to_impl(self, next_t: sage.Expression):
delta: sage.Expression = next_t - self.t
if delta == 0:
return
for circle in self.circles:
circle.center = sage.vector([
circle.center[0] + circle.velocity[0] * delta,
circle.center[1] + circle.velocity[1] * delta + Sim.GRAVITY * delta ** 2 / 2
])
circle.velocity = sage.vector([
circle.velocity[0],
circle.velocity[1] + delta * Sim.GRAVITY
])
self.t = next_t
def save_img(self):
pixels_per_unit = 100
bounds_min = (-3, -3)
bounds_max = (6, 5)
def coor_convert(p: sage.vector):
px = (p[0].n() - bounds_min[0]) * pixels_per_unit
py = (bounds_max[1] - p[1].n()) * pixels_per_unit
return (px, py)
img = Image.new("RGBA", ((bounds_max[0] - bounds_min[0]) * pixels_per_unit, (bounds_max[1] - bounds_min[1]) * pixels_per_unit))
draw = ImageDraw.Draw(img)
draw.rectangle([(0,0), img.size], fill = (0, 0, 0) )
for segment in self.segments:
draw.line([*coor_convert(segment.p1), *coor_convert(segment.p2)], width=5, fill=(255, 255, 255))
for circle in self.circles:
draw.ellipse(
[
*coor_convert(circle.center - sage.vector([circle.radius, -circle.radius])),
*coor_convert(circle.center + sage.vector([circle.radius, -circle.radius]))
],
fill=(255, 255, 255)
)
img.save(f'img{self.t.n()}.png', 'PNG')
def circle_collision_next_time(c1: SimCircle, c2: SimCircle):
sage.forget()
print('circle collision next time in')
t = sage.var('t')
sage.assume(t > 0)
sage.assume(t, 'real')
eq = sage.sqrt( (c1.velocity[0] * t + c1.center[0] - c2.velocity[0] * t - c2.center[0]) ** 2 + (c1.velocity[1] * t + c1.center[1] - c2.velocity[1] * t - c2.center[1]) ** 2) == c1.radius + c2.radius
collision_times = eq.solve(t, to_poly_solve=True)
print('circle collision next time out')
return min([c.rhs() for c in collision_times]) if len(collision_times) > 0 else sage.oo
def circle_segment_collision_next_time(c: SimCircle, s: SimSegment, g: sage.Expression):
sage.forget()
if c.first:
print('circle-segment collision next time in')
print(c, s, g)
t = sage.var('t')
# sage.assume(t > 0)
# sage.assume(t, 'real')
p = sage.vector([c.center[0] + c.velocity[0] * t, c.center[1] + c.velocity[1] * t + g * t ** 2 / 2])
# Find collisions with the line defined by the segment.
eq = sage.abs_symbolic((s.p2[0] - s.p1[0]) * (p[1] - s.p1[1]) - (p[0] - s.p1[0]) * (s.p2[1] - s.p1[1])) / s.length() == c.radius
if c.first:
print('line equation', eq)
collision_times = eq.solve(t, to_poly_solve=True)
#####
# ^ TODO -- BUG! ^
#####
# to_poly_solve is missing solutions. e.g. the following:
# t = var('t')
# f = 2/3*abs(-1/245760000*sqrt(2)*((4*(5*sqrt(-32*sqrt(2) + 288) - 72)^2 - sqrt(2)*(sqrt(2)*((5*sqrt(-32*sqrt(2) + 288) - 72)^2 + 900*sqrt(2) + 720*sqrt(-32*sqrt(2) + 288) - 12384) + 1800) + 3600*sqrt(2) + 2880*sqrt(-32*sqrt(2) + 288) - 49536)*(sqrt(2)*(sqrt(2)*((5*sqrt(-32*sqrt(2) + 288) - 72)^2 + 900*sqrt(2) + 720*sqrt(-32*sqrt(2) + 288) - 12384) + 1800) - 3600*sqrt(2))*t*sqrt(-32*sqrt(2) + 288)/(abs(-1/2400*(5*sqrt(-32*sqrt(2) + 288) - 72)^2 + 1/9600*sqrt(2)*(sqrt(2)*((5*sqrt(-32*sqrt(2) + 288) - 72)^2 + 900*sqrt(2) + 720*sqrt(-32*sqrt(2) + 288) - 12384) + 1800) - 3/8*sqrt(2) - 3/10*sqrt(-32*sqrt(2) + 288) + 129/25)^2 + abs(-1/9600*sqrt(2)*(sqrt(2)*((5*sqrt(-32*sqrt(2) + 288) - 72)^2 + 900*sqrt(2) + 720*sqrt(-32*sqrt(2) + 288) - 12384) + 1800) + 3/8*sqrt(2))^2) + 69120000*sqrt(2) - 552960000) + 1/245760000*sqrt(2)*(((4*(5*sqrt(-32*sqrt(2) + 288) - 72)^2 - sqrt(2)*(sqrt(2)*((5*sqrt(-32*sqrt(2) + 288) - 72)^2 + 900*sqrt(2) + 720*sqrt(-32*sqrt(2) + 288) - 12384) + 1800) + 3600*sqrt(2) + 2880*sqrt(-32*sqrt(2) + 288) - 49536)^2*sqrt(-32*sqrt(2) + 288)/(abs(-1/2400*(5*sqrt(-32*sqrt(2) + 288) - 72)^2 + 1/9600*sqrt(2)*(sqrt(2)*((5*sqrt(-32*sqrt(2) + 288) - 72)^2 + 900*sqrt(2) + 720*sqrt(-32*sqrt(2) + 288) - 12384) + 1800) - 3/8*sqrt(2) - 3/10*sqrt(-32*sqrt(2) + 288) + 129/25)^2 + abs(-1/9600*sqrt(2)*(sqrt(2)*((5*sqrt(-32*sqrt(2) + 288) - 72)^2 + 900*sqrt(2) + 720*sqrt(-32*sqrt(2) + 288) - 12384) + 1800) + 3/8*sqrt(2))^2) - 46080000*sqrt(-32*sqrt(2) + 288))*t - 276480000*t^2 - 76800*(5*sqrt(-32*sqrt(2) + 288) - 72)^2 + 69120000*sqrt(2) - 55296000*sqrt(-32*sqrt(2) + 288) + 951091200)) == (1/3)
# solve(f, t, to_poly_solve=True)
# test on https://sagecell.sagemath.org/
# with to_poly_solve, returns [] — without, returns a simplified equation:
# 1/16*abs(18*sqrt(2)*t^2 + 3*sqrt(2)*t*sqrt(-32*sqrt(2) + 288) - 36*sqrt(2) - 8) == (1/2)
# that Wolfram can find real solutions to. adding assumptions doesn't help
#####
# maybe even worse, if you simplify the equation and square it on both sides you get:
# t = var('t')
# f = 1/256*(18*sqrt(2)*t^2 + 3*sqrt(2)*t*sqrt(-32*sqrt(2) + 288) - 36*sqrt(2) - 8)**2 == (1/4)
# s = solve(f, t, to_poly_solve=True)
# [s.rhs().n() for s in s]
# which returns four complex solutions. but if you plug the same equation into Wolfram, you get back the same 4 real solutions as the non-squared version, the smallest positive one is the one we want.
# does Sage sometimes return WRONG solutions, and exclude solutions even when to_poly_solve is False?
#####
# Filter out collisions where the collision point is not on the segment.
collision_times = [t_val for t_val in collision_times if s.contains_box(s.projection(p.subs(t = t_val.rhs())))]
if c.first:
print('collision times', collision_times)
# Add collisions with the endpoints themselves. This may include times where the circle contains part/all of the segment.
# We don't really care about that, though: as long as we always process the earliest event, we shouldn't be letting
# circles pass through segments to begin with.
for endpoint in [s.p1, s.p2]:
# TODO: can we square both sides of this kind of equation?
eq = sage.sqrt((endpoint[0] - p[0]) ** 2 + (endpoint[1] - p[1]) ** 2) == c.radius
if c.first:
print('endpoint equation', eq)
endpoint_solutions = eq.solve(t, to_poly_solve=True)
collision_times.extend(endpoint_solutions)
if c.first:
print('collision times w/ endpoints', collision_times)
collision_times = [c for c in collision_times if c.rhs() in sage.RR and c.rhs() > sage.Integer(0)]
if c.first:
print('collision times after filter', collision_times)
print('circle-segment collision next time out')
return min([c.rhs() for c in collision_times]) if len(collision_times) > 0 else sage.oo
def perform_circle_collision(c1: SimCircle, c2: SimCircle, e: sage.Expression):
dv = c2.velocity - c1.velocity
normal = (c2.center - c1.center).normalized()
vn = dv[0] * normal[0] + dv[1] * normal[1]
j = (1 + e) * vn / (1 / c1.radius + 1 / c2.radius)
c1.velocity += j * normal / c1.radius
c2.velocity -= j * normal / c2.radius
def perform_circle_point_collision(c: SimCircle, p: sage.vector, e: sage.Expression):
normal = (c.center - p).normalized()
vn = c.velocity[0] * normal[0] + c.velocity[1] * normal[1]
c.velocity -= normal * vn * (1 + e)
# c1 = SimCircle(sage.vector([0, 0]), 1, sage.vector([0, 0]))
# c2 = SimCircle(sage.vector([5, 0]), 1, sage.vector([-1, 0]))
# print(circle_collision_next_time(c1, c2))
# s = SimSegment(sage.vector([-2, -5]), sage.vector([20, -6]))
# print(circle_segment_collision_next_time(c1, s, -1))
sim = Sim(sage.pi / 4, sage.pi / 4)
while sim.state == SimState.INCOMPLETE:
print(f'iterating: {sim.t}')
period = sim.iterate()
if (sim.state == SimState.COMPLETE_WITH_PERIOD):
print(f'we actually found a period: {period}')
print(sim.state)
Raw
sagemath_juggling_v2.py
import os
from itertools import combinations
from PIL import Image, ImageDraw
from typing import cast, Dict
from enum import Enum
from operator import add
import time
import math
import glob
import contextlib
import sys
import sage.all as sage
type Value = sage.QQbar
type Vector = tuple[Value, Value]
os.system("clear")
class SimCircle:
center: Vector
radius: Value
velocity: Vector
def __init__(self, center: Vector, radius: Value, velocity: Vector = (sage.QQbar(0), sage.QQbar(0))):
self.center = center
self.radius = radius
self.velocity = velocity
self.first = False
def get_position_delta(self, t, g):
position_with_t = list(map(add, self.center, tuple(i * t for i in self.velocity)))
position_with_t[1] += g * t ** 2 / 2
return position_with_t
def to_tuple(self):
tuple = (self.center[0], self.center[1], self.velocity[0], self.velocity[1])
return tuple
def __str__(self):
return f'SimCircle center: {self.center}, radius: {self.radius}, velocity: {self.velocity}'
class SimSegment:
p1: Vector
p2: Vector
def __init__(self, p1: Vector, p2: Vector):
self.p1 = p1
self.p2 = p2
def project_point_onto_line(self, p):
# project point onto line (thanks MBo on StackOverflow: https://stackoverflow.com/a/64330724/999658)
abx = self.p2[0] - self.p1[0]
aby = self.p2[1] - self.p1[1]
acx = p[0] - self.p1[0]
acy = p[1] - self.p1[1]
coeff = (abx*acx + aby*acy) / (abx*abx+aby*aby)
dx = self.p1[0] + abx * coeff
dy = self.p1[1] + aby * coeff
return (dx, dy)
def __str__(self):
return f'SimSegment p1: {self.p1}, p2: {self.p2}'
def vector_add(v1: Vector, v2: Vector):
return (v1[0] + v2[0], v1[1] + v2[1])
def vector_subtract(v1: Vector, v2: Vector):
return (v1[0] - v2[0], v1[1] - v2[1])
def vector_normalize(v: Vector):
length = sage.sqrt(v[0] ** 2 + v[1] ** 2)
return (v[0] / length, v[1] / length)
def vector_dot(v1: Vector, v2: Vector):
return (v1[0] * v2[0] + v1[1] * v2[1])
def vector_scale(v: Vector, e: Value):
return (v[0] * e, v[1] * e)
def vector_simplify(v: Vector):
v[0].simplify()
v[1].simplify()
return v
class SimState(Enum):
INCOMPLETE = 1
COMPLETE_WITH_PERIOD = 2
INVALID_MULTIPLE_COLLISIONS = 3 # this doesn't really happen in practice; we're not bothering to detect it
INVALID_SPAWN_OVERLAP = 4
INVALID_TIMEOUT = 5
INVALID_OVER_COLLISION_BUDGET = 6
class Sim:
SPAWN_RATE: Value = sage.QQbar(1)
TIMEOUT: Value = SPAWN_RATE * 20 # the max findable period will be ~3 less than this value
COLLISION_BUDGET = 200
SPAWN_HEIGHT: Value = sage.QQbar(3)
DESPAWN_HEIGHT: Value = sage.QQbar(-3)
CIRCLE_RADIUS: Value = sage.QQbar(1) / 3
GRAVITY: Value = sage.QQbar(-3)
SEGMENT_LENGTH: Value = sage.QQbar(1) / 3
SEGMENT_SEPARATION: Value = sage.QQbar(5) / 2
COEFFICIENT_OF_RESTITUTION: Value = sage.QQbar(4) / 5
FRAME_TIME = sage.QQbar(1)
time_cc = 0
time_cs = 0
t: Value
collisions: int
circles: list[SimCircle]
segments: list[SimSegment]
next_spawn_t: Value
previous_configurations: Dict[list, Value]
state: SimState
gif: bool
def __init__(self, theta1:Value, theta2:Value, gif=False):
self.t = sage.QQbar(0)
self.collisions = 0
self.circles = []
left_segment_center = (sage.QQbar(0), sage.QQbar(0))
right_segment_center = (Sim.SEGMENT_SEPARATION, sage.QQbar(0))
left_offset = sage.QQbar(sage.cos(theta1) / 2), sage.QQbar(sage.sin(theta1) / -2)
right_offset = sage.QQbar(sage.cos(theta2) / -2), sage.QQbar(sage.sin(theta2) / -2)
self.segments = [
SimSegment(vector_subtract(left_segment_center, left_offset), vector_add(left_segment_center, left_offset)),
SimSegment(vector_add(right_segment_center, right_offset), vector_subtract(right_segment_center, right_offset)),
]
self.next_spawn_t = sage.QQbar(0)
self.previous_configurations = []
self.state = SimState.INCOMPLETE
self.gif = gif
def iterate(self):
collisions = []
for circle_pair in combinations(self.circles, 2):
start = time.time()
next_collision = self.t + circle_collision_next_time(*circle_pair)
Sim.time_cc += time.time() - start
collisions.append([next_collision, *circle_pair])
for segment in self.segments:
for circle in self.circles:
start = time.time()
next_collision = self.t + circle_segment_collision_next_time(circle, segment, Sim.GRAVITY)
Sim.time_cs += time.time() - start
collisions.append([next_collision, circle, segment])
min_time_collisions = []
min_time = 9999
for collision in collisions:
if collision[0] < min_time:
min_time_collisions = []
min_time = collision[0]
if collision[0] <= min_time:
min_time_collisions.append(collision)
print('found min_time', min_time, type(min_time))
if min_time >= self.next_spawn_t:
print('spawning at t=', self.next_spawn_t)
self.advance_t_to(self.next_spawn_t)
if not self.has_room_to_spawn():
self.state = SimState.INVALID_SPAWN_OVERLAP
return -1
self.next_spawn_t = self.next_spawn_t + Sim.SPAWN_RATE
self.circles.append(SimCircle((sage.QQbar(0), Sim.SPAWN_HEIGHT), Sim.CIRCLE_RADIUS))
if len(self.circles) == 1:
self.circles[0].first = True
return self.check_configuration()
min_time.simplify()
self.advance_t_to(min_time)
# sys.exit(0)
self.collisions += 1
if self.collisions > Sim.COLLISION_BUDGET:
self.state = SimState.INVALID_OVER_COLLISION_BUDGET
return -1
for collision in min_time_collisions:
c1 = cast(SimCircle, collision[1])
if isinstance(collision[2], SimSegment):
seg = cast(SimSegment, collision[2])
projection = seg.project_point_onto_line(c1.center)
projection = vector_simplify(projection)
perform_circle_point_collision(c1, projection, Sim.COEFFICIENT_OF_RESTITUTION)
else:
c2 = cast(SimCircle, collision[2])
perform_circle_collision(c1, c2, Sim.COEFFICIENT_OF_RESTITUTION)
c1.velocity = vector_simplify(c1.velocity)
c2.velocity = vector_simplify(c2.velocity)
self.circles = [c for c in self.circles if c.center[1] > Sim.DESPAWN_HEIGHT]
return -1
def has_room_to_spawn(self):
for circle in self.circles:
print(circle)
distance_squared = circle.center[0] ** 2 + (circle.center[1] - Sim.SPAWN_HEIGHT) ** 2
diameter_squared = (Sim.CIRCLE_RADIUS * 2) ** 2
if distance_squared <= diameter_squared:
return False
return True
def check_configuration(self):
configuration = [c.to_tuple() for c in self.circles]
configuration.sort()
for prev_con in self.previous_configurations:
if configuration == prev_con[1]:
prev_config_time = prev_con[0]
config_delta = self.t - prev_config_time
self.state = SimState.COMPLETE_WITH_PERIOD
return config_delta / Sim.SPAWN_RATE
self.previous_configurations.append([self.t, configuration])
if self.t > Sim.TIMEOUT:
self.state = SimState.INVALID_TIMEOUT
return -1
def advance_t_to(self, next_t: Value):
if self.t == 0:
self.save_img()
next_frame = math.ceil(self.t / Sim.FRAME_TIME) * Sim.FRAME_TIME
while self.gif and next_frame < next_t:
self.advance_t_to_impl(next_frame)
next_frame = next_frame + Sim.FRAME_TIME
self.save_img()
self.advance_t_to_impl(next_t) # simplifying here seems to take more time than it saves?
def advance_t_to_impl(self, next_t: Value):
delta: Value = next_t - self.t
if delta == 0:
return
for circle in self.circles:
circle.center = vector_simplify(circle.get_position_delta(delta, Sim.GRAVITY))
circle.velocity = vector_simplify((circle.velocity[0], circle.velocity[1] + delta * Sim.GRAVITY))
self.t = next_t
def save_img(self):
if not self.gif:
return
pixels_per_unit = 100
bounds_min = (-3, -3)
bounds_max = (6, 5)
def coor_convert(p: Vector):
px = (p[0] - bounds_min[0]) * pixels_per_unit
py = (bounds_max[1] - p[1]) * pixels_per_unit
return (px, py)
img = Image.new("RGBA", ((bounds_max[0] - bounds_min[0]) * pixels_per_unit, (bounds_max[1] - bounds_min[1]) * pixels_per_unit))
draw = ImageDraw.Draw(img)
draw.rectangle([(0,0), img.size], fill = (0, 0, 0) )
for segment in self.segments:
p1 = [p.n() for p in segment.p1]
p2 = [p.n() for p in segment.p2]
draw.line([*coor_convert(p1), *coor_convert(p2)], width=5, fill=(255, 255, 255))
radius_n = float(Sim.CIRCLE_RADIUS.real())
for circle in self.circles:
center = circle.center
coords = [
*coor_convert(list( map(add, center, (-radius_n, radius_n)) )),
*coor_convert(list( map(add, center, (radius_n, -radius_n)) ))
]
draw.ellipse(
coords,
fill=(255, 255, 255)
)
img.save(f'img{float(self.t):05.2f}.png', 'PNG')
manual_checks = True
def circle_collision_next_time(c1: SimCircle, c2: SimCircle):
if manual_checks:
dx_is_negative = (c2.center[0] - c1.center[0]) < 0
ddx_is_negative = (c2.velocity[0] - c1.velocity[0]) < 0
dy_is_negative = (c2.center[1] - c1.center[1]) < 0
ddy_is_negative = (c2.velocity[1] - c1.velocity[1]) < 0
if dx_is_negative == ddx_is_negative and dy_is_negative == ddy_is_negative:
return 9999
solution = circle_circle_quinn_solve(c1, c2)
if solution == -1:
return 9999
return solution
def circle_segment_collision_next_time(c: SimCircle, s: SimSegment, g: Value):
if manual_checks:
if c.center[0] < s.p1[0] - Sim.CIRCLE_RADIUS and c.velocity[0] <= 0:
return 9999
if c.center[0] > s.p2[0] + Sim.CIRCLE_RADIUS and c.velocity[0] >= 0:
return 9999
if c.center[1] < min(s.p1[1], s.p2[1]) - Sim.CIRCLE_RADIUS:
return 9999
solution = circle_segment_quinn_solve(c, s, g)
if solution == -1:
return 9999
return solution
def circle_circle_quinn_solve(c1: SimCircle, c2: SimCircle):
print(f'qc-c for {c1} {c2}')
dxi = c1.center[0] - c2.center[0]
dyi = c1.center[1] - c2.center[1]
dvx = c1.velocity[0] - c2.velocity[0]
dvy = c1.velocity[1] - c2.velocity[1]
A = dvx**2 + dvy**2
B = dxi * dvx + dyi * dvy
C = sage.sqrt(4 * A * Sim.CIRCLE_RADIUS**2 - (dvy * dxi - dvx * dyi) ** 2)
roots = [-(B + C) / A, (C - B) / A]
# get smallest positive real
roots = [r for r in roots if r.imag().is_zero() and r > 0]
if len(roots) == 0:
return -1
ret = min(roots)
print(f'qc-c done')
return ret
def circle_segment_quinn_solve(c: SimCircle, s: SimSegment, g: Value):
print(f'qc-s for {c} {s}')
x = c.center[0]
y = c.center[1]
vx = c.velocity[0]
vy = c.velocity[1]
x1 = s.p1[0]
x2 = s.p2[0]
y1 = s.p1[1]
y2 = s.p2[1]
r = Sim.CIRCLE_RADIUS
g = g
DX = (x1 - x2)
DY = (y1 - y2)
DX2 = DX ** 2
DY2 = DY ** 2
A = sage.sqrt(r ** 2 / (DX2 + DY2))
roots = []
D = vy * x1 - vy * x2 - vx * y1 + vx * y2
E = g * DX
B = sage.sqrt(DY2 * vx ** 2 + DX2 * vy ** 2 + 2 * DX2 * g * (A * DX - y) + 2 * DX * (vx * vy * (y2 - y1) + g * (A * DY2 - x2 * y1 + x * DY + x1 * y2)))
roots.append((B + D) / -E)
roots.append((B - D) / E)
C = sage.sqrt(DX2 * vy ** 2 + DY2 * (vx ** 2 + 2 * A * g * (x2 - x1)) - 2 * DX2 * g * (A * DX + y) + 2 * DX * (vx * vy * (y2 - y1) + g * (x * y1 - x2 * y1 - x * y2 + x1 * y2)))
roots.append((C + D) / -E)
roots.append((C - D) / E)
# filter out hits not on the segment
def inside_segment_box(p, x1, x2, y1, y2):
# project point onto line (thanks MBo on StackOverflow: https://stackoverflow.com/a/64330724/999658)
(dx, dy) = s.project_point_onto_line(p)
# loose "on line" check
min_x = min(x1, x2)
max_x = max(x1, x2)
min_y = min(y1, y2)
max_y = max(y1, y2)
return dx >= min_x and dx <= max_x and dy >= min_y and dy <= max_y
roots = [r for r in roots if r > 0]
roots = [r for r in roots if inside_segment_box(c.get_position_delta(r, g), x1, x2, y1, y2)]
# add point collisions
def point_roots(x1, y1):
print('point_roots input types', type(x1), type(y1))
print(f'qc-sp for {x1} {y1}')
GSQ = g**2
AA = vx * (x - x1) + vy * (y - y1)
DD = vx ** 2 + vy ** 2 + g * (y - y1)
EE = r ** 2 - x ** 2 + 2 * x * x1 - x1 ** 2 - y ** 2 + 2 * y * y1 - y1 ** 2
FF = 16 * DD ** 2 - 48 * EE * GSQ - 96 * AA * g * vy
GG = 128 * DD ** 3 + 1728 * AA ** 2 * GSQ + 1152 * DD * EE * GSQ - 1152 * AA * DD * g * vy - 1728 * EE * GSQ * vy ** 2
HH = sage.sqrt(-4 * FF ** 3 + GG ** 2)
def cube_root(d):
return d ** (sage.QQ(1) / 3)
II = cube_root(GG + HH)
cbrt2 = cube_root(sage.QQbar(2))
JJ = sage.sqrt(-((8 * DD)/(3 * GSQ)) + ((cbrt2 * FF)/(3 * GSQ * II)) + (II/(3 * cbrt2 * GSQ)) + ((4 * vy ** 2)/GSQ))
RADICAL_FIRST_HALF = -((16 * DD)/(3 * GSQ)) - (cbrt2 * FF)/(3 * GSQ * II) - II/(3 * cbrt2 * GSQ) + (8 * vy ** 2)/GSQ
RADICAL_SECOND_HALF = (-((64 * AA)/GSQ) + (64 * DD * vy)/g**3 - ((64 * vy**3)/g**3))/(4 * JJ)
roots = []
MINUS_RADICAL = sage.sqrt(RADICAL_FIRST_HALF - RADICAL_SECOND_HALF)
roots.append(-JJ/2 - vy/g - MINUS_RADICAL / 2)
roots.append(-JJ/2 - vy/g + MINUS_RADICAL / 2)
PLUS_RADICAL = sage.sqrt(RADICAL_FIRST_HALF + RADICAL_SECOND_HALF)
roots.append(JJ/2 - vy/g - PLUS_RADICAL / 2)
roots.append(JJ/2 - vy/g + PLUS_RADICAL / 2)
print(f'qc-sp done')
return roots
roots.extend(point_roots(x1, y1))
roots.extend(point_roots(x2, y2))
# get smallest positive real
for r in roots:
print('roottype', type(r))
roots = [r for r in roots if r.imag().is_zero() and r > 0]
if len(roots) == 0:
return -1
ret = min(roots)
print(f'qc-s done')
return ret
def perform_circle_collision(c1: SimCircle, c2: SimCircle, e: Value):
e = (e + 1) / 2
normal = vector_normalize(vector_subtract(c2.center, c1.center))
impulse = vector_dot(normal, vector_subtract(c2.velocity, c1.velocity))
impulse = vector_scale(impulse, e)
dv = vector_scale(impulse, normal)
c1.velocity = vector_simplify(vector_add(c1.velocity, dv))
c2.velocity = vector_simplify(vector_subtract(c2.velocity, dv))
def perform_circle_point_collision(c: SimCircle, p: Vector, e: Value):
normal = vector_normalize(vector_subtract(c.center, p))
dot = vector_dot(c.velocity, normal)
print(f'before c-p collision: {c}')
print(f'normal: {normal}, dot: {dot}')
c.velocity = vector_simplify(vector_subtract(c.velocity, vector_scale(normal, (1 + e) * dot)))
print(f'after c-p collision: {c}')
DEGREES_X = sage.QQ(45)
DEGREES_Y = sage.QQ(45)
DEGREE_RANGE = sage.QQ(88) # the entire span of angles: 45, 45, 2 -> [44 degrees, 46 degrees]
csv_filename = 'out.csv'
denominator = sage.QQ(1)
x_numerator = sage.QQ(0)
y_numerator = sage.QQ(0)
fresh_file = False
if os.path.isfile(csv_filename):
with open(csv_filename) as f:
lines = f.readlines()
if len(lines) > 1:
tokens = lines[-1].split(',')
denominator = sage.QQ(int(tokens[0]))
x_numerator = sage.QQ(int(tokens[1]))
y_numerator = sage.QQ(int(tokens[2]))
else:
fresh_file = True
with open(csv_filename, 'x') as f:
f.write('denominator,x_numerator,y_numerator,x_approx_degrees,y_approx_degrees,sim_state,period')
while True:
# increment numerators and denominator
if fresh_file:
fresh_file = False
else:
if y_numerator % 2 == 0 and denominator > 1:
x_numerator += 2
else:
x_numerator += 1
if x_numerator > denominator:
y_numerator += 1
if y_numerator > denominator:
denominator *= 2
y_numerator = 0
x_numerator = 1 - y_numerator % 2
# run simulation
x_degrees = 45#DEGREES_X - DEGREE_RANGE / 2 + DEGREE_RANGE * x_numerator / denominator
y_degrees = 45#DEGREES_Y - DEGREE_RANGE / 2 + DEGREE_RANGE * y_numerator / denominator
start = time.time()
sim = Sim(sage.pi / 180 * x_degrees, sage.pi / 180 * y_degrees, True)
while sim.state == SimState.INCOMPLETE:
print(f'iterating: {sim.t}')
period = sim.iterate()
if (sim.state == SimState.COMPLETE_WITH_PERIOD):
print(f'we actually found a period: {period}')
print(f'final sim state: {sim.state}')
print(f'time: {time.time() - start}s')
print(f'cctime: {Sim.time_cc}s, cstime: {Sim.time_cs}s')
if sim.gif:
with contextlib.ExitStack() as stack:
imgs = (stack.enter_context(Image.open(f))
for f in sorted(glob.glob('img*.png'))[int(period) * -10:])
img = next(imgs)
img.save(fp='img.gif', format='GIF', append_images=imgs, save_all=True, duration=50, loop=0)
for f in glob.glob("img*.png"):
os.remove(f)
break
else:
with open(csv_filename, "a") as myfile:
myfile.write(f"\n{denominator},{x_numerator},{y_numerator},{x_degrees},{y_degrees},{sim.state},{int(period)}")
Raw
sympy_juggling.py
# rewrite in SymEngine? https://github.com/symengine/symengine
# seems to not have an analogue to solve()
# or: https://www.sagemath.org/
# we might be able to specify a range for t in sp.solve to reduce solve times?
# algo:
# - start at t=0 with two line segments and an empty list of circles, and repeat the following
# - find collision times between all pairs of circles and all circle/segment pairs
# - filter out all collisions that don't take place at the min time
# - if the min time is less than or equal to the next spawn time:
# - if any one circle is involved in multiple collisions, the result is INVALID: MULTIPLE COLLISIONS
# - process all collisions
# - advance time to the min collision time
# - if the min time is greater than the next spawn time:
# - advance time to the next spawn time
# - remove all circles with y < some negative number
# - if we are at a spawn time, spawn a new circle
# - if the new circle overlaps an existing circle, the result is INVALID: SPAWN OVERLAP
# - check if we've seen this configuration of circles before
# - if we have, return (elapsed time since we last saw this configuration) / (spawn period)
# - otherwise, record the time and configuration
# - if the time is over some maximum, the result is INVALID: TIMEOUT
import sympy as sp
from itertools import combinations
from PIL import Image, ImageDraw
from typing import cast, Dict
from enum import Enum
import os
os.system("clear")
class SimCircle:
circle: sp.Circle
velocity: sp.Point
def __init__(self, circle: sp.Circle):
self.circle = circle
self.velocity = sp.Point2D(0, 0)
def to_tuple(self):
return [self.circle.center.x, self.circle.center.y, self.velocity.x, self.velocity.y]
def __str__(self) -> str:
return f'SimCircle {self.circle}, {self.velocity}'
class SimState(Enum):
INCOMPLETE = 1
COMPLETE_WITH_PERIOD = 2
INVALID_MULTIPLE_COLLISIONS = 3
INVALID_SPAWN_OVERLAP = 4
INVALID_TIMEOUT = 5
class Sim:
SPAWN_RATE: sp.Number = 1
TIMEOUT: sp.Number = SPAWN_RATE * 50
SPAWN_HEIGHT: sp.Number = 3
DESPAWN_HEIGHT: sp.Number = -10
CIRCLE_RADIUS: sp.Number = sp.Rational(1, 3)
GRAVITY: sp.Number = -3
SEGMENT_LENGTH: sp.Number = sp.Rational(3, 2)
SEGMENT_SEPARATION: sp.Number = 3
COEFFICIENT_OF_RESTITUTION: sp.Number = sp.Rational(9, 10)
FRAME_TIME = sp.Rational(1, 10)
t: sp.Number
circles: list[SimCircle]
segments: list[sp.Segment2D]
next_spawn_t: sp.Number
previous_configurations: Dict[list, sp.Number]
state: SimState
def __init__(self, theta1:sp.Number, theta2:sp.Number):
self.t = sp.S.Zero
self.circles = []
left_segment_center = sp.Point2D(0, 0)
right_segment_center = sp.Point2D(Sim.SEGMENT_SEPARATION, 0)
self.segments = [
sp.Segment2D(left_segment_center - sp.Point2D(sp.cos(theta1), -sp.sin(theta1)) * Sim.SEGMENT_LENGTH / 2, left_segment_center + sp.Point2D(sp.cos(theta1), -sp.sin(theta1)) * Sim.SEGMENT_LENGTH / 2),
sp.Segment2D(right_segment_center - sp.Point2D(sp.cos(theta2), sp.sin(theta2)) * Sim.SEGMENT_LENGTH / 2, right_segment_center + sp.Point2D(sp.cos(theta2), sp.sin(theta2)) * Sim.SEGMENT_LENGTH / 2)
]
self.next_spawn_t = sp.S(0)
self.previous_configurations = dict()
self.state = SimState.INCOMPLETE
def iterate(self):
collisions = []
for circle_pair in combinations(self.circles, 2):
next_collision = self.t + sp.Min(*circle_collision_times(*circle_pair))
collisions.append([next_collision, *circle_pair])
for segment in self.segments:
for circle in self.circles:
next_collision = self.t + sp.Min(*circle_line_segment_collision_times(circle, segment, Sim.GRAVITY))
collisions.append([next_collision, circle, segment])
min_time_collisions = []
min_time: sp.Number = sp.oo
for collision in collisions:
if collision[0] < min_time:
min_time_collisions = []
min_time = collision[0]
if collision[0] <= min_time:
min_time_collisions.append(collision)
print(f'min collision time: {min_time.evalf()}, next spawn: {self.next_spawn_t.evalf()}')
if min_time >= self.next_spawn_t:
self.advance_t_to(self.next_spawn_t)
self.next_spawn_t += Sim.SPAWN_RATE
# TODO: check for circles in spawn area
self.circles.append(SimCircle(sp.Circle(sp.Point2D(0, Sim.SPAWN_HEIGHT), Sim.CIRCLE_RADIUS)))
return -1
# TODO: check for circles involved in multiple collisions
self.advance_t_to(min_time)
for collision in min_time_collisions:
c1 = cast(SimCircle, collision[1])
if isinstance(collision[2], sp.Segment2D):
seg = cast(sp.Segment2D, collision[2])
perform_circle_point_collision(c1, seg.projection(c1.circle.center), Sim.COEFFICIENT_OF_RESTITUTION)
else:
c2 = cast(SimCircle, collision[2])
perform_circle_collision(c1, c2, Sim.COEFFICIENT_OF_RESTITUTION)
self.circles = [c for c in self.circles if c.circle.center.y > Sim.DESPAWN_HEIGHT]
configuration = [c.to_tuple() for c in self.circles]
configuration.sort()
configuration = str(configuration)
if configuration in self.previous_configurations:
prev_config_time = self.previous_configurations[configuration]
config_delta = prev_config_time - self.t
self.state = SimState.COMPLETE_WITH_PERIOD
return config_delta / Sim.SPAWN_RATE
self.previous_configurations[configuration] = self.t
if self.t > Sim.TIMEOUT:
self.state = SimState.INVALID_TIMEOUT
return -1
def advance_t_to(self, next_t: sp.Number):
if (self.t == 0):
self.save_img();
print(f'advancing to t={next_t.evalf()}')
next_frame = sp.ceiling(self.t / Sim.FRAME_TIME) * Sim.FRAME_TIME
while (next_frame < next_t):
self.advance_t_to_impl(next_frame)
next_frame += Sim.FRAME_TIME
self.save_img();
self.advance_t_to_impl(next_t)
def advance_t_to_impl(self, next_t: sp.Number):
next_t = sp.simplify(next_t)
delta: sp.Number = next_t - self.t
if delta.is_zero:
return
for circle in self.circles:
circle.circle = sp.Circle(
sp.Point2D(sp.simplify(circle.circle.center.x + circle.velocity.x * delta),
sp.simplify(circle.circle.center.y + circle.velocity.y * delta + Sim.GRAVITY * delta ** 2 / 2)),
circle.circle.radius
)
circle.velocity = sp.Point2D(
circle.velocity.x,
sp.simplify(circle.velocity.y + delta * Sim.GRAVITY),
)
self.t = next_t
def save_img(self):
pixels_per_unit = 100
bounds_min = (-3, -3)
bounds_max = (6, 5)
def coor_convert(p: sp.Point2D):
px = (p.x.evalf() - bounds_min[0]) * pixels_per_unit
py = (bounds_max[1] - p.y.evalf()) * pixels_per_unit
return (px, py)
img = Image.new("RGBA", ((bounds_max[0] - bounds_min[0]) * pixels_per_unit, (bounds_max[1] - bounds_min[1]) * pixels_per_unit))
draw = ImageDraw.Draw(img)
draw.rectangle([(0,0), img.size], fill = (0, 0, 0) )
for segment in self.segments:
draw.line([*coor_convert(segment.p1), *coor_convert(segment.p2)], width=5, fill=(255, 255, 255))
for circle in self.circles:
draw.ellipse([*coor_convert(circle.circle.center - sp.Point2D(circle.circle.radius, -circle.circle.radius)),
*coor_convert(circle.circle.center + sp.Point2D(circle.circle.radius, -circle.circle.radius))],
fill=(255, 255, 255))
img.save(f'img{self.t.evalf()}.png', 'PNG')
def circle_collision_times(c1: SimCircle, c2: SimCircle):
t = sp.symbols('t', real=True, positive=True)
eq = sp.Eq(sp.sqrt( (c1.velocity.x * t + c1.circle.center.x - c2.velocity.x * t - c2.circle.center.x) ** 2 + (c1.velocity.y * t + c1.circle.center.y - c2.velocity.y * t - c2.circle.center.y) ** 2), c1.circle.radius + c2.circle.radius)
collision_times = sp.solve(eq, t)
collision_times = [time for time in collision_times if time.is_real and time > 0]
return collision_times
def circle_line_segment_collision_times(c: SimCircle, s: sp.Segment2D, g: sp.Number):
print('=== circle_line_segment_collision_times ===')
print(c)
print(s)
t = sp.symbols('t', real=True, positive=True)
p = sp.Point2D(c.circle.center.x + c.velocity.x * t, c.circle.center.y + c.velocity.y * t + g * t ** 2 / 2)
# Find collisions with the line defined by the segment.
line = sp.Line2D(s)
eq = sp.Eq(line.distance(p), c.circle.radius)
print('colliding with line')
collision_times = list(sp.solve(eq, t))
# Filter out collisions where the collision point is not on the segment.
collision_times = [t_val for t_val in collision_times if s.contains(line.projection(p.subs(t, t_val)))]
# Add collisions with the endpoints themselves. This may include times where the circle contains part/all of the segment.
# We don't really care about that, though: as long as we always process the earliest event, we shouldn't be letting
# circles pass through segments to begin with.
for endpoint in [s.p1, s.p2]:
eq = sp.Eq(sp.simplify(endpoint.distance(p)), c.circle.radius)
print(f'colliding with endpoint: {eq}')
endpoint_solutions = sp.solve(eq, t, manual=True)
print(f'got back: {endpoint_solutions}')
collision_times.extend(endpoint_solutions)
collision_times = [time for time in collision_times if time.is_real and time > 0]
print('=== done ===')
return collision_times
def perform_circle_collision(c1: SimCircle, c2: SimCircle, e: sp.Number):
dv = c2.velocity - c1.velocity
normal = (c2.circle.center - c1.circle.center).unit
vn = dv.x * normal.x + dv.y * normal.y
j = (1 + e) * vn / (1 / c1.circle.radius + 1 / c2.circle.radius)
c1.velocity += j * normal / c1.circle.radius
c2.velocity -= j * normal / c2.circle.radius
def perform_circle_point_collision(c: SimCircle, p: sp.Point2D, e: sp.Number):
normal = (c.circle.center - p).unit
vn = c.velocity.x * normal.x + c.velocity.y * normal.y
c.velocity -= normal * vn * (1 + e)
sim = Sim(sp.pi / 4, sp.pi / 4)
while sim.state == SimState.INCOMPLETE:
print(f'iterating: {sim.t}')
period = sim.iterate()
if (period > 0):
print(f'we actually found a period: {period}')
print(sim.state)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment