Genetic algorithm for finding bounce paths

Prototype: Bounce GA
https://github.com/yDgunz/BounceGA

Overview

Bounce juggling against a set of arbitrary surfaces seemed like a cool feature to add to the gunswap animator. However, given the scope of this addition, it seemed logical to first create a prototype that can calculate a single bounce path given a set of surfaces, start and end position, and a number of other inputs discussed below.

Simulating the Bounce Path

In order to simulate the bounce path we first need to define some parameters, primarily the surfaces and the start and end positions of the ball. Assuming the surfaces are flat, we can describe them by a position vector, a normal vector, and a scale. For simplicity we will also assume that the surfaces have an edge parallel to the x-z axis. The figure below (TODO) shows the start position (p0), final position (pT), and a surface described by position vector (Sp), normal vector (Sn) and scale (Ss).

The bounce path begins with the ball leaving p0 with a velocity of v0. Moving forward in time in increments of dt we can calculate a new position p(t+dt) = p(t)+v(t)*dt and a new velocity v(t+dt) = v(t) + g*dt. This applies until we reach the first bounce.

The first bounce occurs when the ball is colliding with a surface .We can determine if a ball and surface are colliding if the distance from position of the ball to any point on the surface is less than or equal to the radius of the ball. We can describe this condition as (p(t)-Sp) dot Sn <= R. The figure below shows this (TODO).

Assuming an elastic collision, the new velocity will be v(t) - Sn*C*2*v(t) dot Sn, where C is the coefficient of restitution (link). From here we can continue to calculate the new position and velocity as described above until any other collisions occur.

Given the correct velocity v0 we will ultimately end up at pT at time t = T. So this whole problem boils down to finding the correct velocity v0.

We could find v0 by generating random velocities and running the simulation described above until we've found a velocity where the final position in the simulation falls within some minimum success threshold. However, considering a three dimensional bounce path with a small success threshold, this method may take too long to find a solution. It seemed perfect for application of a genetic algorithm.

Genetic Algorithm

A genetic algorithm is an optimization technique that mimics the process of natural selection. The approach described below may be more aptly described as an evolutionary algorithm or differential evolution. I'm not really familiar with the terminology and ultimately whatever you call this heuristic, it got the job done.

Let's first define a population of candidate solutions called G0. Each member of G0 is an initial velocity vector for which we run the simulation described above. The members are assigned a fitness which is the distance from pT the ball ends up if it starts with the member's initial velocity.

If none of the members of G0 have a fitness that falls within the success threshold, we need to generate a new population G1. We do this by selecting members from G0 with a bias towards members with a higher fitness. We then randomly mutate these members and calculate their new fitness values. If none of the members of this new generation, G1, meet the success criteria then we repeat this process until a solution is found, or some limit of attempts is met.

Input Syntax

There are two categories of inputs necessary to run this genetic algorithm. We need the parameters for the bounce path simulation (the fitness function) and the parameters for the genetic algorithm such as population size, max number of generations, etc.

The prototype uses the following input structure:


Fitness Parameters

p0.x, p0.y, p0.z, pT.x, pT.y, pT.z, minT, maxT
R, C
numBounces[, bounceOrder]
tossSign, catchSign
normal.x, normal.y, normal.z, position.x, position.y, position.z, scale - (surface 0)
...
normal.x, normal.y, normal.z, position.x, position.y, position.z, scale - (surface N)

• p0 - starting position
• pT - ending position
• minT - minimum amount of time for bounce path
• maxT - maximum amount of time for bounce path
• R - ball radius
• C - ball coefficient of restitution
• numBounces - expected number of bounces
• bounceOrder (optional) - comma separated list of surface indices representing the bounce order. So 1,0 means bounce off surface 1 first, then surface 0
• tossSign - 1 means the toss is up, -1 means the toss is down, u means it doesn't matter
• catchSign - 1 means the catch is made while the prop is traveling up, -1 means the catch is made while the prop is traveling down, u means it doesn't matter
• surfaceNormal - vector for the surface normal
• surfacePosition - vector for the position of the center of the surface
• surfaceScale - all surfaces are squares with an edge parallel to the x-z plane, the scale provides the half-width of the square

GA Parameters

maxGenerations, populationSize, mutationChance, mutationScale, initialScale, fitnessThreshold, noGA

• maxGenerations - max number of generations to run GA for
• populationSize - total size of each generation
• mutationChance - percentage chance that a member will mutate
• mutationScale - maximum magnitude of mutation
• initialScale - maximum magnitude of newly generated members
• fitnessThreshold - maximum distance from pT to be considered a success
• noGA (optional) - set to 1 to turn off selection by fitness and mutation (ie. iterating through random solutions)

Conclusion

The genetic algorithm has proven to be an adequate solution to the arbitrary bounce path problem. Next steps are to integrate this into the gunswap animator and to figure out the best way of inputting all the necessary parameters. Thus far most of the gunswap inputs have been handled through form elements, but there's a lot here and I may consider reworking some of the existing inputs in favor of a more text based approach like in the bounce path prototype.

v8 and v9

Released v8 a while back, but never made the release notes. Just pushed a v9 tonight for some stuff that was sitting around locally. For these two releases the big updates were:

- Changed prop toss/spin orientation syntax. The siteswap now just encodes the number of spins, so you can say "5S3" for 5 club triples, but you can't say "5S3XY" to define the toss/spin orientation. The hand movement string is now used for the prop orientation. You can now add {x,y,z,th} to the end of each hand movement point to define the axis of rotation (x,y,z) and the rotation (th). There's a lot going on with this, and it's kind of half finished, so I'll follow-up in a future blog post.
- Changed the default prop toss/spin orientation. It now looks much more realistic for clubs. Even if I don't continue on with the ability to change the orientation, I think the default now looks pretty nice.
- Added a bunch of advanced controls for the dwell path parameters. You can now modify the impact the toss/catch velocities have on the dwell path, including whether or not to match velocities on toss and catch. These parameters are really meaningless to the average user who would just be interested in a good looking animation, so I'm going to see if there are any magic numbers that always look nice.
- Added more examples

v7

The big feature in v7 is moving elbows. See previous blog post for more details about this. Also included some other minor changes, but nothing particularly noteworthy.

Finding Elbow Positions

Intro

One of the glaring issues with the gunswap animator was the lack of elbow movement. It's fairly easy to calculate the hand positions of a juggler throughout a pattern; they are well defined by the dwell path of the pattern, and when a prop is in flight the hand is just moving towards the next catch. The elbows are a bit trickier though. My goal was to define a function that will return the elbow position given the hand position.

This problem is called "inverse kinematics" and it is fairly well studied. However, a google search on the topic yields papers like this (http://www.ro.feri.uni-mb.si/predmeti/robotizacija/knjiga/inverzna.pdf). As someone who got a C in linear algebra over 7 years ago, I felt a bit in over my head.

Since Juggling Lab has already solved this problem, I decided to sift through their source code. I've used Juggling Lab as a reference for a lot of parts of this project - the bulk of my siteswap syntax is based on the Juggling Lab syntax. After a little searching I found the section of code I was looking for: http://sourceforge.net/p/jugglinglab/code/HEAD/tree/trunk/source/jugglinglab/renderer/Juggler.java#l128.

Unfortunately I was having a hard time following the Juggling Lab code, and didn't want to just start copy/pasting things into gunswap without understanding them. So I sought help from the main developer of Juggling Lab, Jack Boyce. Here's what he had to say:

"If you think about the connection between the shoulder and the hand as a linkage of two stiff pieces (upper arm and lower arm), then the position of the elbow becomes constrained once the hand and shoulder locations are specified. There is only one degree of freedom left, the rotation of the entire arm around the shoulder-hand axis. (Imagine a chicken flapping its wings to show the motion I'm talking about.)"

Up to this point I had been trying to work out complicated solutions considering the individual rotation angles of the shoulder and elbow. However, thinking about the the problem with one degree of freedom - the "chicken wing angle" - was the key to unlocking the solution.

Solution

Below is how I came up with a function getElbowPosition(S,H,l,w,hand) that returns elbow coordinates given inputs for the shoulder position (S), hand position (H), forearm/bicep length (l), "chicken wing" angle (w) and hand. The hand input is just 0/1 (left/right) indicating which direction the chicken wing angle should rotate. A higher chicken wing angle means your elbows are lifted higher. A chicken wing angle of 0 means the plane created by your shoulder/elbow/hand is parallel to the y axis. The juggler in Fig 1. has a high chicken wing angle.

Fig 1.

Fig 2.

Fig 3.

Fig 4.

Below is a step-by-step explanation of the equations in Fig 4.

1) The first step is to transform H to a coordinate system where S is at the origin. This gives a new hand position H'.

2) Another transformation. H'' is the hand position in a coordinate system rotated along the y-axis (the dashed lines in Fig. 2). H'' is a simpler vector to work with because it has a 0 z component.

3) Theta is used to go back to H' from H''. This will be useful later on. Note that in the code I actually use Math.atan2 since Math.atan only goes from -pi/2 to pi/2 and Math.atan2 determines the angle based on the quadrant the input coordinates are in.

4) h is the distance from the shoulder/hand axis to the elbow. This is the height of an isosceles triangle with sides of length l. http://mathworld.wolfram.com/IsoscelesTriangle.html

5) Now we find unit vectors u1 and u2 that are both orthogonal to each other and the shoulder/hand axis. These are simple to find because H'' has a 0 z component.

6) Here is where the chicken wing angle comes into play. This is the position of the elbow as defined by an equation for a circle around the shoulder/hand axis. http://math.stackexchange.com/questions/73237/parametric-equation-of-a-circle-in-3d-space

7) Transforming E'' to E' using theta which we solved for in (3). 

8) Transforming E' to E which is our final result.

Explanation of gunswap code

Here's an explanation of the gunswap code as of v6.0.

Siteswap.js

Siteswap.js is the core siteswap parsing library. It exports 1 function:

CreateSiteswap(siteswapStr [, options])

This function takes a set of inputs describing a siteswap and returns a JSON object with calculated properties of the siteswap.

Inputs

siteswapStr is the only required input and is a siteswap string that adheres to the gunswap notation.

options is a JSON object containing additional options.

Examples

var s1 = CreateSiteswap('531');
var s2 = CreateSiteswap('51',{validationOnly:true});
var s3 = CreateSiteswap('441',{beatDuration:.3,dwellRation:.5});

You don't need to specify all of the options. Here are the list of options and their defaults.

validationOnly: true 
Passing in false will prevent the generation of prop positions and rotations.

numSteps: 1000
Defines the size of the prop positions/rotations arrays. A value of 1000 will calculate the props' positions at 1000 equally spaced steps in the full cycle of the pattern.

beatDuration: .2
Defines the length of a single beat in seconds. A value of .2 will result in a 3 ball cascade with throws lasting .6 seconds.

dwellRatio: .5
Defines the amount of time (as a ratio of the beat duration) the a prop is held between throws. The dwell duration cannot exceed the beat duration, hence it's input as a ratio.

props: [{type: 'ball', radius: .05, C: .95}]
Array of props. Valid types are 'ball', 'club' and 'ring'. Obviously radius and C don't really make any sense for clubs/rings. C if the coefficient of restitution (ie. the bounciness) and is just used for determining the ball path in bounce patterns. If the input array length doesn't match the number of props for the pattern the last object in the array is continuously appended or removed to get a match.

dwellPath: 
This requires more explanation - perhaps another blog post soon.

Return object

The CreateSiteswap function returns a big JSON object with the following properties.

siteswap
The same siteswap string that was provided as an input.

validSyntax
validPattern
multiplex
sync
pass
All of the above are just booleans describing the pattern.

numJugglers
numProps
maxHeight
More fairly simple properties of the siteswap.

tosses
Array of tosses that occur at each beat. 

beats
Siteswap string broken down into its component beats.

states
State array used for pattern validation. See previous blog post.

propOrbits
For each prop an array of throws and what beat they occur at. 

propPositions
Array of x,y,z positions for each prop.

jugglerHandPositions
Array of x,y,z positions for each juggler's hands.

jugglers
Array of jugglers and their x,z positions (y is assumed to be 0).

validationOnly
numSteps
beatDuration
dwellDuration
props
dwellPath
The inputs.

errorMessage
If everything worked, this will be undefined. Otherwise this contains any error messages.

The inner workings of this function are an enigma. But the basic flow is:

1. Validate syntax (http://gunswap.blogspot.com/2014/02/using-regex-to-parse-siteswap-strings.html)
2. Validate pattern (http://gunswap.blogspot.com/2014/03/siteswap-validation.html)
3. Generate prop positions (no post yet)

SiteswapAnimator.js

SiteswapAnimator.js handles all the animation. It exports the SiteswapAnimator class.

Constructor

SiteswapAnimator(containerId)

containerId is the div that the canvas will be appended to. The constructor handles all the scene initialization.

Functions

go(siteswap[,options])

This function runs the animation. siteswap is supposed to be the return object from a call to CreateSiteswap. options is a JSON object containing any animation options. It currently just contains a boolean motionBlur.

zoomIn()
zoomOut()
resize(width,height)
updateAnimationSpeed(animationSpeed)
updateCameraMode(cameraMode)

All of the above function should be fairly obvious. animationSpeed is a value from 0 to 1 where 1 is full speed and 0 is paused. cameraMode is either 'sky' or 'juggler'.

Example

var animator = new SiteswapAnimator('myDiv');
animator.go(SiteswapJS.CreateSiteswap('531'));

index.js

This just instantiates a SiteswapAnimator, reads user inputs, and passes the results of a CreateSiteswap call to the animator.

util.js

This contains general helper code.

v6.0

Quick turnaround on this one, though not much to report. Another new layout for the UI along with the ability to change color/type of individual props. Also added motion blur.

v5.0

Version 5:

- Added support for siteswaps up to "z" (for passing use a capital "P")
- Changed dwell duration to dwell ratio
- Animation automatically starts on page load. If you pass in a siteswap in the URL query string then it'll use that siteswap. For example: ydgunz.github.io/gunswap?siteswap=531. Make sure to encode the URI component (just use encodeURIComponent), so to animate 24[54] you would put in ydgunz.github.io/gunswap?siteswap=%5B33%5D.
- Cubic bezier dwell paths! This is one of the ideas that I've been excited about for a while and it's finally coming together. This will need it's own post to explain everything, but the general idea is that using cubic bezier interpolation to define dwell paths will allow patterns like mills mess, factory, etc. I've added a bezier version of the cascade and the 1-count factory (I think some people call this "waterfall"?) to the hand movement select list. More to come on this.
- I spent some time trying to make more of an app out of this so that folks can save patterns. But I'm not really sure where that's going so I stopped. I think I'm going to spend some time documenting the existing code and figuring out how to clean it up some more. Oh the luxuries of a project with no deadline.