Skip to content

Instantly share code, notes, and snippets.

@ubuntor

ubuntor/approx_bounce.py

Created September 15, 2024 01:00
Show Gist options
  • Save ubuntor/c7de0572d0e68fb5dcf4e94d3efba7d4 to your computer and use it in GitHub Desktop.
Save ubuntor/c7de0572d0e68fb5dcf4e94d3efba7d4 to your computer and use it in GitHub Desktop.
Download ZIP
Raw
approx_bounce.py
import math
import numpy as np
import jax.numpy as jnp
from cr.sparse.pursuit import omp
NUM_SAMPLES = 512 # reduce this if it's using too much memory!
def resample(signal, samples):
resampled = []
for i in range(samples):
x = i/samples
for j in range(len(signal)):
if signal[j][0] > x:
x0, y0 = signal[j-1]
x1, y1 = signal[j]
y = y0 + ((x-x0)/(x1-x0)) * (y1-y0)
resampled.append(y)
break
return resampled
BOUNCE_SAMPLES = 2048
bounce_raw = []
assert BOUNCE_SAMPLES % 2 == 0
# generate first half of bounce and reflect it to get the second half
for i in range(BOUNCE_SAMPLES//2 + 1):
t = 2 * i/BOUNCE_SAMPLES
x = 0.5 * (0*((1-t)**3) + 0.8*3*t*((1-t)**2) + 1*3*(t**2)*(1-t) + 1*(t**3))
y = -0.25 * (1 - (0*((1-t)**3) + 0*3*t*((1-t)**2) + 1*3*(t**2)*(1-t) + 1*(t**3)))
bounce_raw.append((x,y))
bounce_raw += [(1-x,y) for x,y in bounce_raw[:-1][::-1]]
bounce = resample(bounce_raw, NUM_SAMPLES)
def harmonic(signal, n):
return np.array([signal[(i*n) % len(signal)] for i in range(len(signal))])
def normalize(signal):
mean = signal.mean()
signal -= mean
norm = np.linalg.norm(signal)
return signal / norm, norm, mean
def approximate(target, num_components):
target = jnp.array(np.array(resample(target, NUM_SAMPLES)))
atom_params = [] # frequency, phase
atoms = []
for harmonic_num in range(1, NUM_SAMPLES//2 + 1): # bounce is symmetric across x=0.5, so higher harmonics are redundant
h, norm, mean = normalize(harmonic(bounce, harmonic_num))
for rotation in range(NUM_SAMPLES // math.gcd(NUM_SAMPLES, harmonic_num)): # some phases are redundant
phase = rotation/NUM_SAMPLES
atoms.append(np.roll(h, -rotation))
atom_params.append((norm, mean, harmonic_num, phase))
print(f'{len(atoms)} atoms')
atoms = jnp.array(np.array(atoms).transpose())
solution = omp.matrix_solve_jit(atoms, target, num_components)
print(f'found solution: residual norm (squared): {solution.r_norm_sqr}')
components = []
total_offset = 0
for coeff, idx in zip(solution.x_I, solution.I):
norm, mean, frequency, phase = atom_params[idx]
coeff /= norm
total_offset += coeff * mean
components.append((coeff, frequency, phase))
return components
def make_1d(target, total_period, num_components):
components = approximate(target, num_components)
components.sort(key=lambda x: (x[0] < 0, -abs(x[0]))) # sort positive then negative, both by decreasing magnitude
is_negated = False
num_divs = len(components)
for i, (coeff, frequency, phase) in enumerate(components):
if coeff < 0 and not is_negated:
print('<div style="transform: rotate(180deg); pointer-events:none">')
is_negated = True
num_divs += 1
height = abs(coeff)
period = total_period * (1 / frequency)
phase *= total_period
delay = (phase % period) - period
color = 'hsl(210deg 100% 75%)' if is_negated else 'hsl(0deg 100% 75%)'
width = 100 + 5*(len(components)-i-1)
print(f'<div style="animation: {period}s bounce {delay}s infinite; height: {height}px; width: {width}px; border-radius: 25px; border: 2px solid {color}; display: flex; justify-content: center; align-items: center; pointer-events:none">')
if is_negated:
print('<div style="transform: rotate(180deg); font-size: 70px; pointer-events:none">:eggbug:</div>')
else:
print('<div style="font-size: 70px; pointer-events:none">:eggbug:</div>')
print('</div>'*num_divs)
def make_2d(target_x, target_y, total_period, num_components):
components_x = approximate(target_x, num_components)
components_x.sort(key=lambda x: (x[0] < 0, -abs(x[0])))
components_y = approximate(target_y, num_components)
components_y.sort(key=lambda x: (x[0] < 0, -abs(x[0])))
# this is garbage, but i'm speedrunning writing this, so...
is_negated = False
num_divs = len(components_y)
for i, (coeff, frequency, phase) in enumerate(components_y):
if coeff < 0 and not is_negated:
print('<div style="transform: rotate(180deg); pointer-events:none">')
is_negated = True
num_divs += 1
height = abs(coeff)
period = total_period * (1 / frequency)
phase *= total_period
delay = (phase % period) - period
color = 'hsl(210deg 100% 75%)' if is_negated else 'hsl(0deg 100% 75%)'
width = 100 + 5*(len(components_y)-i-1)
print(f'<div style="animation: {period}s bounce {delay}s infinite; height: {height}px; width: {width}px; border-radius: 25px; border: 2px solid {color}; display: flex; justify-content: center; align-items: center; pointer-events:none">')
print('<div style="transform: rotate(90deg); pointer-events:none">')
num_divs += 1
is_negated = False
num_divs += len(components_x)
for i, (coeff, frequency, phase) in enumerate(components_x):
if coeff < 0 and not is_negated:
print('<div style="transform: rotate(180deg); pointer-events:none">')
is_negated = True
num_divs += 1
height = abs(coeff)
period = total_period * (1 / frequency)
phase *= total_period
delay = (phase % period) - period
color = 'hsl(210deg 100% 75%)' if is_negated else 'hsl(0deg 100% 75%)'
width = 100 + 5*(len(components_x)-i-1)
print(f'<div style="animation: {period}s bounce {delay}s infinite; height: {height}px; width: {width}px; border-radius: 25px; border: 2px solid {color}; display: flex; justify-content: center; align-items: center; pointer-events:none">')
if is_negated:
print('<div style="transform: rotate(180deg); font-size: 70px; pointer-events:none">:eggbug:</div>')
else:
print('<div style="font-size: 70px; pointer-events:none">:eggbug:</div>')
print('</div>'*num_divs)
# target: pairs of (t, y) where 0 <= t <= 1, y is in pixels (positive y is up)
# target will be linearly interpolated
target = [(0, -100), (0.5, 100), (1, -100)] # example: triangle wave
#target = [(i/256, 100*math.sin(2*math.pi * i/256)) for i in range(257)] # example: sine wave
#target = [(0, -100), (0.5, -100), (0.5, 100), (1, 100)] # example: square wave
#target = [(0, -100), (1, 100)] # example: sawtooth wave
# example: square motion
#target_x = [(0, -100), (0.25, 100), (0.5, 100), (0.75, -100), (1, -100)]
#target_y = [(0, -100), (0.25, -100), (0.5, 100), (0.75, 100), (1, -100)]
# example: circle motion
#target_x = [(i/256, 100*math.cos(2*math.pi * i/256)) for i in range(257)]
#target_y = [(i/256, 100*math.sin(2*math.pi * i/256)) for i in range(257)]
# example: diamond motion
#target_x = [(0, -100), (0.5, 100), (1, -100)]
#target_y = [(0, 0), (0.25, -100), (0.5, 0), (0.75, 100), (1, 0)]
# example: dvd screensaver
#target_x = [(i/512, 250*(abs((i/256*5)%2-1)-0.5)) for i in range(513)]
#target_y = [(i/512, 200*(abs((i/256*4)%2-1)-0.5)) for i in range(513)]
'''
# example: eggbug
import svgpathtools
paths, attrs = svgpathtools.svg2paths('shadow_eggbug.svg')
target_x = []
target_y = []
for i in range(4097):
t = i/4096
p = paths[0].point(t)
target_x.append((t, 20*p.real))
target_y.append((t, 20*p.imag))
'''
make_1d(target, 4, 30)
#make_2d(target_x, target_y, 15, 30)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment