Simulates particles bonking into one another in a hilariously inefficient O(n^2) way. Dependencies: numpy, scipy, pygame

particles.py

import pygame
import math
import numpy
import random
import scipy
import scipy.spatial

SCREEN_SIZE = 300
clock = pygame.time.Clock()
screen = pygame.display.set_mode((SCREEN_SIZE, SCREEN_SIZE))#, pygame.FULLSCREEN)
running = True
physics_framerate = 50
physics_diff = 2
tick = 0

DEFAULT_FLIP = 4.0
DEFAULT_PUSHY = 0.2

# Atoms
ps_raw = []
ds_raw = []

# Springs
# a_1_index, a_2_index, pushy/pully flipover distance, pushy force per distance, pully force per distance, snap distance
springs = []

# Let's just make some soup
for x in range(100):
    ps_raw += [(random.random() * SCREEN_SIZE, random.random() * SCREEN_SIZE)]
    ds_raw += [((random.random() - 0.5) * 0.2, (random.random() - 0.5) * 0.2)]
    
def calc_offsets(xy1, xy2):
    d0 = numpy.subtract.outer(xy1[:,0], xy2[:,0])
    d1 = numpy.subtract.outer(xy1[:,1], xy2[:,1])
    return numpy.vstack(([d0.T], [d1.T])).T

# Numpyfication
ps = numpy.array(ps_raw)
ds = numpy.array(ds_raw)

flip_raw = [[DEFAULT_FLIP] * len(ps_raw)] * len(ps_raw)
push_raw = [[DEFAULT_PUSHY] * len(ps_raw)] * len(ps_raw)

flip = numpy.array(flip_raw)
push = numpy.array(push_raw)

numpy.seterr(all='ignore')

while running:
    tick += 1
    clock.tick(physics_framerate)
    event = pygame.event.poll()
    if event.type == pygame.QUIT:
        running = False
    
    # Physics
    # Movement
    ps += ds
    
    ps %= SCREEN_SIZE
    
    # Springs
    distances = scipy.spatial.distance.cdist(ps, ps)
    offsets = calc_offsets(ps, ps)
    # Scale offsets
    offsets /= numpy.vstack(([distances.T], [distances.T])).T
    # Get rid of nans
    offsets = numpy.nan_to_num(offsets)
    # Figure out which atoms actually push on one another
    distances -= flip
    push_bitmap = numpy.less(distances, 0)
    distances *= push_bitmap
    # Scale the distances by the push amount
    distances *= push
    # Calculate the push vectors for each relation
    offsets *= numpy.vstack(([distances.T], [distances.T])).T
    # Flatten push vectors into a list to apply to the delta-vs
    ds += numpy.add.reduce(offsets)
    
    # Draw
    if tick % physics_diff == 0:
        screen.fill((0, 0, 0))
        for p in ps:
            screen.set_at(p, (255, 255, 255))
        pygame.display.flip()