![[OLD FALL 2018] 15-104 • Introduction to Computing for Creative Practice](wp-content/uploads/2020/08/stop-banner.png)
These examples (provided by Golan, with only minor changes by Roger) illustrate the power of mutual interaction.
So far this week, we looked at particles that obey rules independently. Next, we considered springs where particles could be coupled and influenced by additional forces. Now, we consider interactions between every pair of particles!
Mutual Repulsion with Optional Gravity
Here, we create nParticles particles with mutual interaction. The idea here is that particles try to move away from each other. To implement this idea, we represent “try to move away” as a force. The closer particles are together, the greater the force that pushes them apart. In fact, we use the inverse square law and make the force proportional to 1/(distance*distance). As distance increases, this force approaches zero, but at small distances, e.g. distance = 1, it is significant. In fact, the force is infinite at zero distance, so to avoid huge forces and the resulting high velocities, we do not apply the force when distance is very small.
Implementation
We use a nested loop: for each particle, we want to compute the force applied by each other particle. All the forces add together, and in fact, the final net force may be zero due to particles pushing in opposite directions.
Looking at the inner for loop, you would expect j < myParticles.length, but instead, j runs from 0 to i, where i is the index of the particle. Why? If we want to compute the force from every particle, why does the inner loop only consider some of the particles? In fact, the inner loop could loop over every particle, but notice that when you know the force of particle A on particle B, you also know the force of particle B on particle A. This particular implementation optimizes the computation by computing forces only in one direction, from let’s say A to B. When force is applied to A, we simply apply the negative force to B, eliminating the need to (re)compute half of the forces.
Looking at the end of the inner loop, you see
|
1 2 |
ithParticle.addForce(repulsionForcex, repulsionForcey); // add in forces jthParticle.addForce(-repulsionForcex, -repulsionForcey); |
In the inner loop, j is always less than i, which eliminates half of the interactions, but we make up for this by updating two particles on each iteration of the inner loop. (Also, we eliminate the force of a particle on itself, i.e. where i === j.
// ========== PARTICLE METHODS =============
// Update the position based on force and velocity
function particleUpdate() {
if (this.bFixed == false) {
this.vx *= this.damping;
this.vy *= this.damping;
this.limitVelocities();
this.handleBoundaries();
this.px += this.vx;
this.py += this.vy;
}
}
// Prevent particle velocity from exceeding maxSpeed
function particleLimitVelocities() {
if (this.bLimitVelocities) {
var speed = sqrt(this.vx * this.vx + this.vy * this.vy);
var maxSpeed = 10;
if (speed > maxSpeed) {
this.vx *= maxSpeed / speed;
this.vy *= maxSpeed / speed;
}
}
}
// do boundary processing if enabled
function particleHandleBoundaries() {
if (this.bPeriodicBoundaries) {
if (this.px > width) this.px -= width;
if (this.px < 0) this.px += width;
if (this.py > height) this.py -= height;
if (this.py < 0) this.py += height;
} else if (this.bHardBoundaries) {
if (this.px >= width) {
this.vx = -abs(this.vx);
}
if (this.px <= 0) {
this.vx = abs(this.vx);
}
if (this.py >= height) {
this.vy = -abs(this.vy);
}
if (this.py <= 0) {
this.vy = abs(this.vy);
}
}
}
// draw the particle as a white circle
function particleDraw() {
fill(255);
ellipse(this.px, this.py, 9, 9);
}
// add a force to the particle using F = mA
function particleAddForce(fx, fy) {
var ax = fx / this.mass;
var ay = fy / this.mass;
this.vx += ax;
this.vy += ay;
}
// make a new particle
function makeParticle(x, y, dx, dy) {
var p = {px: x, py: y, vx: dx, vy: dy,
mass: 1.0, damping: 0.9,
bFixed: false,
bLimitVelocities: false,
bPeriodicBoundaries: false,
bHardBoundaries: false,
addForce: particleAddForce,
update: particleUpdate,
limitVelocities: particleLimitVelocities,
handleBoundaries: particleHandleBoundaries,
draw: particleDraw
}
return p;
}
// ========== THE MAIN PROGRAM ============
// Mutual repulsion, with optional gravity
var myParticles = [];
var nParticles = 400;
function setup() {
createCanvas(400, 400);
for (var i = 0; i < nParticles; i++) {
var rx = random(width);
var ry = random(height);
myParticles[i] = makeParticle(rx, ry, 0, 0);
}
}
function keyPressed() {
for (var i = 0; i < myParticles.length; i++) {
myParticles[i].px = random(width);
myParticles[i].py = random(height);
}
}
function draw() {
background(200);
var gravityForcex = 0;
var gravityForcey = 0.05;
var mutualRepulsionAmount = 2.5;
for (var i = 0; i < myParticles.length; i++) {
var ithParticle = myParticles[i];
var px = ithParticle.px;
var py = ithParticle.py;
if (mouseIsPressed) {
ithParticle.addForce(gravityForcex, gravityForcey);
}
for (var j = 0; j < i; j++) {
var jthParticle = myParticles[j];
var qx = jthParticle.px;
var qy = jthParticle.py;
var dx = px - qx;
var dy = py - qy;
var dh = sqrt(dx * dx + dy * dy);
if (dh > 1.0) {
var componentInX = dx / dh;
var componentInY = dy / dh;
var proportionToDistanceSquared = 1.0 / (dh * dh);
var repulsionForcex = mutualRepulsionAmount *
componentInX * proportionToDistanceSquared;
var repulsionForcey = mutualRepulsionAmount *
componentInY * proportionToDistanceSquared;
// add in forces
ithParticle.addForce( repulsionForcex, repulsionForcey);
jthParticle.addForce(-repulsionForcex, -repulsionForcey);
}
}
}
for (var i = 0; i < myParticles.length; i++) {
myParticles[i].bPeriodicBoundaries = false;
myParticles[i].bElasticBoundaries = true;
myParticles[i].update(); // update all locations
}
for (var i = 0; i < myParticles.length; i++) {
myParticles[i].draw(); // draw all particles
}
}
Efficiency
How long does it take to compute the draw function? Let’s think about what’s involved:
- There are fixed costs per draw: clearing the canvas, initializing variables, etc.
- There are per-particle costs: computing forces on the particle, drawing each particle.
In many cases, graphics operations dominate the cost of computation. Think about the cost of computing a few numbers like x, y, velocity compared to the cost of computing and drawing the outline of an ellipse pixel-by-pixel, and filling the ellipse, again pixel-by-pixel. Filling an area with pixels is at least some kind of nested loop, and the work is at least proportional to the number of pixels involved.
If there are not too many particles, you could correctly assume the cost of draw() is proportional to the number of particles. In computer science, we use N to represent the “size of the problem,” e.g. the number of particles (or N could be the number of characters of text to be processed by some text processing program, or N could be the number of data records to be searched in a database application, etc., but here, N is clearly the number of particles.)
One of the great achievements of computer science is the idea that we can think about how run-time depends upon N. We just said run-time is proportional to N, or in other words, run time is roughly N times some number, e.g. maybe it’s N * 0.0001s, meaning drawing takes about 0.0001 seconds per particle. (Yes, there’s also some fixed costs to clear the canvas, etc., but these matter less and less as N gets larger and larger, so the cost is dominated by N.)
But wait! What about the inner loop and force computation? How many force updates are there? If there are N particles, and each particle requires N/2 force computations, the total number of steps is N * N/2. Yikes. Removing the “/2” which is just a scale factor that does not change with N, we have N^2 (N-squared) operations. If N is small and computation is fast, this is no big deal compared to the huge overhead of drawing all those pixels for each particle, but what if N is larger? Suppose N = 10,000? How many force computations? 10,000 * 10,000 = 100,000,000. Double yikes! My computer takes about 2.5 seconds to do this much work. That’s nowhere close to animation rates. What can you do about it? Probably using another programming language (e.g. C or C++) would speed this up 10-fold. Doing the computation in parallel (more programming, more hardware) could get us up to animation speeds. But then, what about N = 100,000? Now, we’ve added 100 times more work and we’re in trouble again.
Conclusion: Algorithms matter! Computers are not infinitely fast, and you now know enough to write interesting computations that will bring your computer to its knees. Some simple analysis can tell you how your performance will depend on problem size (N). It’s the nature of the scaling (is it proportional to N? N squared? 2-to-the-N power?) that is most important when N gets large.
Flocking
This is a brief introduction to a very interesting phenomenon called “flocking.” A more complete introduction to flocking is available online, as is a great deal of literature. Our goal is to illustrate the main concepts. We might ask you about flocking on our written exam, but do not expect you to be able to implement a complete flocking algorithm without further study.
Flocking is inspired by flocks of birds that achieve mass coordination of movement even though each bird is making independent decisions. How does the flocking behavior emerge from local behavior? In the flocking model, often referred to as “boids” using the term coined by the inventor of the model, Craig Reynolds, there are 3 rules:
- Separation — avoid collisions with neighbors
- Alignment — steer in the same direction as your neighbors
- Cohesion — steer to the center of your neighbors to stay together
Here’s an example created with a flocking algorithm:
Let’s implement these rules in a particle system. The following is an example. There are comments in the code. This is based on code by Conrad Parker, but I found the separation routine caused jittery motion as particles moved in and out of the neighborhood, so I modified the force calculation to increase smoothly from 0 at the edge of the neighborhood. I also added an obstacle in the middle of the canvas that the particles avoid. Finally, the particle object implementation has some changes to allow forces to be summed before modifying velocity (requiring force to be explicitly zeroed at the beginning of each draw). For more info, see comments in the code.
In this program,
- “S” key randomizes particle locations,
- “G” key toggles grouping, and
- “F” key toggles following.
// ========== PARTICLE METHODS =============
// Update the position based on force and velocity
function particleUpdate() {
if (this.bFixed == false) {
this.vx *= this.damping;
this.vy *= this.damping;
this.limitVelocities();
this.handleBoundaries();
this.px += this.vx;
this.py += this.vy;
}
}
// Prevent particle velocity from exceeding maxSpeed
function particleLimitVelocities() {
if (this.bLimitVelocities) {
var speed = sqrt(this.vx * this.vx + this.vy * this.vy);
var maxSpeed = 10;
if (speed > maxSpeed) {
this.vx *= maxSpeed / speed;
this.vy *= maxSpeed / speed;
}
}
}
// do boundary processing if enabled
function particleHandleBoundaries() {
if (this.bPeriodicBoundaries) {
if (this.px > width) this.px -= width;
if (this.px < 0) this.px += width;
if (this.py > height) this.py -= height;
if (this.py < 0) this.py += height;
} else if (this.bHardBoundaries) {
if (this.px >= width) {
this.vx = -abs(this.vx);
}
if (this.px <= 0) {
this.vx = abs(this.vx);
}
if (this.py >= height) {
this.vy = -abs(this.vy);
}
if (this.py <= 0) {
this.vy = abs(this.vy);
}
}
}
// draw the particle as a white circle
function particleDraw() {
fill(255);
ellipse(this.px, this.py, 9, 9);
}
// add a force to the particle using F = mA
function particleAddForce(fx, fy) {
var ax = fx / this.mass;
var ay = fy / this.mass;
this.vx += ax;
this.vy += ay;
}
// make a new particle
function makeParticle(x, y, dx, dy) {
var p = {px: x, py: y, vx: dx, vy: dy,
mass: 1.0, damping: 0.9,
bFixed: false,
bLimitVelocities: false,
bPeriodicBoundaries: false,
bHardBoundaries: false,
addForce: particleAddForce,
update: particleUpdate,
limitVelocities: particleLimitVelocities,
handleBoundaries: particleHandleBoundaries,
draw: particleDraw
}
return p;
}
// ========== THE MAIN PROGRAM ============
// Mutual repulsion, with optional gravity
var myParticles = [];
var nParticles = 30;
var group = true;
var follow = true;
function setup() {
createCanvas(400, 400);
for (var i = 0; i < nParticles; i++) {
var rx = random(width);
var ry = random(height);
var p = makeParticle(rx, ry, 0, 0);
p.bElasticBoundaries = true;
p.damping = 0.99;
myParticles[i] = p;
}
frameRate(10);
}
function keyPressed() {
if (key === "S") {
for (var i = 0; i < myParticles.length; i++) {
myParticles[i].px = random(width);
myParticles[i].py = random(height);
}
} else if (key === "G") {
group = !group;
} else if (key === "F") {
follow = !follow;
}
}
// apply separation rule to particle p
var desiredSeparation = 50;
var separationFactor = 0.01;
function separate(p) {
// rule applies if separation is < desired separation,
// apply an opposing force to achieve greater separation
// opposing force grows as the separation becomes less
for (var i = 0; i < myParticles.length; i++) {
var boid = myParticles[i]; // get each other particle
var d = dist(p.px, p.py, boid.px, boid.py);
if (d > 1 & d < desiredSeparation) {
// divide by distance so that the total force is 1
var fx = (p.px - boid.px) / d;
var fy = (p.py - boid.py) / d;
// scale force by (desiredSeparation - d), which
// is 0 at desiredSeparation and grows as distance
// becomes less
var factor = (desiredSeparation - d) * separationFactor;
p.addForce(fx * factor, fy * factor);
}
}
}
// desired velocity is the average velocity of others
var alignmentFactor = 0.01;
function alignment(p) {
// similar in separate, get the average velocity of neighbors
// (note: we're getting velocity here, not position)
var sumx = 0; // sum of x and y velocities
var sumy = 0;
for (var i = 0; i < myParticles.length; i++) {
var boid = myParticles[i]; // get each other particle
sumx += boid.vx;
sumy += boid.vy;
}
p.addForce((sumx / myParticles.length - p.vx) * alignmentFactor,
(sumy / myParticles.length - p.vy) * alignmentFactor);
}
// steer toward center of flock
var cohesionFactor = 0.0001;
function cohesion(p) {
// get average position of all particles
var sumx = 0; // sum of x and y locations
var sumy = 0;
var count = 0; // how many neighbors
for (var i = 0; i < myParticles.length; i++) {
var boid = myParticles[i]; // get each other particle
sumx += boid.px;
sumy += boid.py;
}
var factor = cohesionFactor;
if (!group) {
factor = -cohesionFactor;
}
p.addForce((sumx / myParticles.length - p.px) * factor,
(sumy / myParticles.length - p.py) * factor);
}
var seekFactor = 0.0003;
function seek(p) {
var desiredX = mouseX - p.px;
var desiredY = mouseY - p.py;
p.addForce(desiredX * seekFactor, desiredY * seekFactor);
}
// avoid is similar to separate but it acts to move away from
// the given point at x, y, with desired separation size s
var avoidFactor = 0.05;
function avoid(p, x, y, s) {
var d = dist(p.px, p.py, x, y);
if (d > 1 & d < s) {
// divide by distance so that the total force is 1
var fx = (p.px - x) / d;
var fy = (p.py - y) / d;
// scale force by (s - d), which
// is 0 at s and grows as distance becomes less
var factor = (s - d) * avoidFactor;
p.addForce(fx * factor, fy * factor);
}
}
//==========================================================
function draw() {
background(200);
text("Type F or G", 10, 15);
if (follow) text("(F)ollow", 10, 30);
if (group) text("(G)roup", 10, 45);
ellipseMode(CENTER);
color(200, 0, 0);
ellipse(height / 2, width / 2, 30, 30);
color(0);
for (var i = 0; i < myParticles.length; i++) {
var ithParticle = myParticles[i];
ithParticle.fx = 0;
ithParticle.fy = 0;
separate(ithParticle);
alignment(ithParticle);
cohesion(ithParticle);
// this avoids the point in the middle of canvas
avoid(ithParticle, height / 2, width / 2, 40);
if (follow) seek(ithParticle);
}
for (var i = 0; i < myParticles.length; i++) {
var p = myParticles[i]
p.update(); // update all locations
p.draw(); // draw all particles
}
}
/*function draw() {
background(200);
var gravityForcex = 0;
var gravityForcey = 0.05;
var mutualRepulsionAmount = 2.5;
for (var i = 0; i < myParticles.length; i++) {
var ithParticle = myParticles[i];
var px = ithParticle.px;
var py = ithParticle.py;
if (mouseIsPressed) {
ithParticle.addForce(gravityForcex, gravityForcey);
}
for (var j = 0; j < i; j++) {
var jthParticle = myParticles[j];
var qx = jthParticle.px;
var qy = jthParticle.py;
var dx = px - qx;
var dy = py - qy;
var dh = sqrt(dx * dx + dy * dy);
if (dh > 1.0) {
var componentInX = dx / dh;
var componentInY = dy / dh;
var proportionToDistanceSquared = 1.0 / (dh * dh);
var repulsionForcex = mutualRepulsionAmount *
componentInX * proportionToDistanceSquared;
var repulsionForcey = mutualRepulsionAmount *
componentInY * proportionToDistanceSquared;
// add in forces
ithParticle.addForce( repulsionForcex, repulsionForcey);
jthParticle.addForce(-repulsionForcex, -repulsionForcey);
}
}
}
for (var i = 0; i < myParticles.length; i++) {
myParticles[i].bPeriodicBoundaries = false;
myParticles[i].bElasticBoundaries = true;
myParticles[i].update(); // update all locations
}
for (var i = 0; i < myParticles.length; i++) {
myParticles[i].draw(); // draw all particles
}
}
*/
/*
// Mutual repulsion, with optional gravity
function keyPressed() {
if (key === "S") {
for (var i = 0; i < myParticles.length; i++) {
myParticles[i].px = random(width);
myParticles[i].py = random(height);
}
} else if (key === "G") {
group = !group;
} else if (key === "F") {
follow = !follow;
}
}
function Particle(){
this.px = 0;
this.py = 0;
this.vx = 0;
this.vy = 0;
this.fx = 0;
this.fy = 0;
this.damping = 0.96;
this.mass = 1.0;
this.bLimitVelocities = true;
this.bPeriodicBoundaries = false;
this.bElasticBoundaries = false;
this.set = function(x,y) {
this.px = x;
this.py = y;
this.vx = 0;
this.vy = 0;
this.damping = 0.96;
this.mass = 1.0;
}
this.addForce = function(fx,fy) {
this.fx += fx;
this.fy += fy;
}
this.update = function() {
this.vx += this.fx/this.mass;
this.vy += this.fy/this.mass;
this.vx *= this.damping;
this.vy *= this.damping;
this.limitVelocities();
this.handleBoundaries();
this.px += this.vx;
this.py += this.vy;
}
this.limitVelocities = function() {
if (this.bLimitVelocities) {
var speed = sqrt(this.vx * this.vx + this.vy * this.vy);
var maxSpeed = 10;
if (speed>maxSpeed) {
this.vx *= maxSpeed/speed;
this.vy *= maxSpeed/speed;
}
}
}
this.inBounds = function() {
return(this.px < width & this.px > 0 &&
this.py < height && this.py > 0);
}
this.handleBoundaries = function() {
if (this.bPeriodicBoundaries) {
if (this.px>width) this.px -= width;
if (this.px<0) this.px += width;
if (this.py>height) this.py -= height;
if (this.py<0) this.py += height;
} else if (this.bElasticBoundaries) {
if (this.px>width) {
this.px = width*2-this.px;
this.vx = -this.vx;
}
if (this.px<0) {
this.px = -this.px;
this.vx = -this.vx;
}
if (this.py>height) {
this.py = height * 2 - this.py;
this.vy = -this.vy;
}
if (this.py<0) {
this.py = -this.py;
this.vy = -this.vy;
}
}
}
this.render = function() {
var d = max(0.1, sqrt(this.vx * this.vx + this.vy * this.vy)) / 10;
fill(0);
line(this.px, this.py, this.px - this.vx / d, this.py - this.vy / d);
}
}
*/