search

borntofrappe / the-nature-of-code

Recreate some of the examples from the book [The Nature of Code](https://natureofcode.com/book/) with Lua and Love2D.
4 stars 0 forks source link

readme

The Nature of Code introduces many concepts to simulate real-world phenomena with code. Here, I follow the book to learn about these concepts.

The notes which follow are my own. The demos are written in Lua and Love2D.

the-nature-of-code repository

Useful Links

00 - Randomness

Randomness provides a first, rudimentary way to simulate real phenomena.

Random

A random function returns a number in a range with the same probability as any number in the same range. The output isn"t truly random, but pseudo-random, whereby the function creates a series of numbers and returns one of them. The sequence repeats itself, but over a long stretch of time.

To move an entity at random, there exist several strategies:

Please note:

Probability

The probability of a single event is given by the number of outcomes representing the event divided by the number of possible outcomes. The probability of multiple events occurring in sequence is obtained by multiplying the single events. The concept is useful to describe the random functions, but also the distribution of other functions.

With the random function, you can obtain a certain probability with at least two methods:

Please note:

Normal distribution

A normal, or Gaussian, distribution returns a number starting from two values: a mean and a standard variation. The idea of the distribution is that numbers are more likely to approach the mean than deviate from it. How often the numbers distance themselves from the mean is described by the standard deviation. With a small deviation, the numbers gather around the mean. Opposedly, and with a large deviation, the numbers scatter away from the central value.

For the normal distribution, it is useful to remember the following:

Please note:

Custom distribution

To fit the needs of a simulation, you can customize a distribution in several ways. The Probability folder already introduces two possibilities, with the methods picking a number from a given set, or using the cumulative probability. Here, a custom distribution is built with the following algorithm:

This approach means that values with a greater probability are more likely to be accepted. The custom nature of the distribution comes from the way the probability is then computed.

It helps to consider the probability as the y value for a function in an (x, y) plane. When the probability is equal to the random value, as in the demo, the relationship is linear. When the probability is instead equal to the value squared, the relationship is squared.

Please note:

Perlin noise

The Perlin noise function allows to create a sequence of numbers connected to each other, with the goal of providing smooth random values. You pick numbers from the sequence, and the distance between the numbers dictates the difference between the two. The greater the distance, the more likely the numbers will differ. The smaller the offset, the more likely the numbers will resemble one another.

While it is possible to create a function implementing the logic of Perlin noise, Love2D provides a similar functionality in love.math.noise.

The function returns a sequence of numbers in the (0,1) range and accepts multiple arguments, to create noise in multiple dimensions. In one dimension, each number is related to the one coming before and after it. In two, it is connected to the numbers representing the possible neighbors in the (x,y) plane. In three, the neighbors considering a third dimension (z) as well.

Please note:

01 - Vectors

Vectors as introduced in the book are euclidean vectors, entities with a magnitude and a direction. They are introduced in the context of a plane with two dimensions, x and y, but fundamentally, they work in the same manner with additional dimensions.

Think of a vector with two components as an arrow. The length of the arrow describe its magnitude, while the angle relative to an axis its direction. A vector describing the position of a particle details where to position the object from the origin. A vector describing the velocity dictates where to move the same particle.

Please note:

Vector math

Vectors follow specific rules to compute mathematical operations.

Please note:

Velocity

With two vectors describing the position and velocity, it is possible to move an object at a constant rate.

Please note:

Acceleration

With three vectors, it is possible to move an object with an accelerated rate, speeding up or slowing down. The idea is to update the position with the velocity, and the velocity with the acceleration.

self.position:add(self.velocity)
self.velocity:add(self.acceleration)

With this setup, the goal of the simulation is to then set a particular acceleration. This value can be set arbitrarily, as with a constant or random value, or following actual physics, as with of gravity or wind strength.

Please note:

02 - Forces

The goal is to adapt the concept of forces and Newtonian physics to the simplified environment introduced with vectors.

Newton"s laws

In the simplified environment, a force is described as a vector which causes an object with mass to accelerate.

Newton"s first law, arguing that an object at rest stays at rest and an object in motion stays in motion, is adapted by saying that the vector describing the velocity stays constant. The only way for an object to stop, to reach an equilibrium is for the velocity to be affected by a series of forces which cancel its magnitude.

Newton"s third law, detailing an action/reaction pair for any force, is partly incorporated in the environment by occasionally including a force contrasting the original one.

Newton"s second law, providing a formula to compute a force on the basis of mass and acceleration, is finally essential to the simulation. This law states that the force is equal to mass times acceleration.

->      ->
F = m * a

In a first approximation, and assuming a mass equal to 1, an immediate way to apply the force is to therefore set the acceleration to the force itself

function Mover:applyForce(force)
  self.acceleration = force
end

Please note:

Force accumulation

The previous solution works, but only when a single force is applied. With multiple forces, only the last one is incorporated in the acceleration vector. The implementation is therefore modified to consider the effect of a force as cumulative (force accumulation)

function Mover:applyForce(force)
  self.acceleration:add(force)
end

Be careful that it is necessary to reset the acceleration vector. In this manner, the object considers the forces available in the specific frame, and not every force accumulated in the simulation.

self.acceleration:multiply(0)

Please note:

Mass

In a slightly more elaborated construct, a force is weighed by the object mass.

function Mover:applyForce(force)
  self.acceleration:add(LVector:divide(force, self.mass))
end

It is important to note, however, that forces like gravity already incorporate the mass in their value. For these forces, it is necessary to remove the mass"s influence.

local gravity = LVector:new(0, GRAVITY)
mover:applyForce(LVector:multiply(gravity, mover.mass))

Please note:

Creating forces

In a simulation, create a force with an arbitrary value, or following the guidance of an actual formula. In this last instance, the simulation needs to compute the magnitude and direction of the force vector from a given set of values.

Friction

Friction is applied on a moving object to progressively reduce its velocity. The actual formula computes the vector by considering the unit vector describing the velocity, a coefficient of friction (mu), and the magnitude of the normal force.

->                       ^
F = -1 * mu * ||N|| * velocity

In a first approximation, it is possible to simplify the force by considering its direction and magnitude. In terms of direction, friction has a direction opposite to the velocity vector. Notice how the unit vector is multiplied by -1.

->          ^
F = -1 * velocity

In terms of magnitude, the force is scaled with a constant value.

->          ^
F = -1 * velocity * c

A more elaborate simulation would try to compute the normal vector, and its eventual magnitude, would incorporate the coefficient according to the surface creating the friction, or again the normal force and its magnitude, but in the approximation, it is enough to scale the vector with a constant. By changing the constant, the simulation is able to then describe a surface with higher/smaller friction.

Please note:

Drag

A force of drag considers the density of the material through which the object is moving, rho, the magnitude of the object"s velocity, ||v||, the surface area subject to resistance A, a coefficient of drag C and the velocity"s unit vector ^v

->                                      ^
F = -1 / 2 * rho * ||v||^2 * A * C * velocity

Similarly to the Friction demo, however, the force can be simplified by considering direction and magnitude. For the direction, the force is again the opposite of the velocity vector.

->              ^
F = -1 / 2 * velocity

For the magnitude, the force is scaled according to the magnitude of the velocity, and a value summarising the other constants.

->          ^
F = -1 * velocity * ||v||^2 * c

Please note:

Gravitational attraction

The force of gravity depends on the mass of the objects involved, m1 and m2, the distance between said objects d, as well as a constant describing the gravitational force G. In terms of direction, the force finally depends on the direction of the vector connecting the two objects, ^r.

->                            ^
F = G * ((m1 * m2) / d ^ 2) * r

The unit vector connecting the objects describes the direction of ths force.

->                            ^
F = G * ((m1 * m2) / d ^ 2) * r

The constant, mass values and then distance finally influence the magnitude of the force.

F = G / (d ^ 2)

Please note:

03 - Oscillation

To discuss oscillating motion, it is first necessary to introduce angles, polar coordinates and trigonometry. Trigonometry relates to the study of the angles and sides of right triangles, and is useful to model angles, angular velocity and angular acceleration.

Angles

To rotate an entity, modify the coordinate system with translations and rotations.

love.graphics.translate(x, y)
love.graphics.rotate(angle)
--

In the snippet, the drawing operations following the two expressions will be rotated by amount described by angle.

Love2D, like the Processing library used in the course, works with an angle in radians, not degrees. Radians describe the angle in terms of the ratio between the length of the arc of a circle and its radius, with 1 radian being the angle at which the arc has the same length as the radius. To convert between the two, use the following formula:

radians = 2 * math.pi * (degrees / 360)

Alternatively, Lua provides the math.rad to convert to radians, math.deg to convert to degrees.

It is perhaps useful to note that pi is the ratio between a circle"s circumference and its diameter. It is roughly 3.14159 and is provided by Lua in the math library and math.pi.

Please note:

Angular motion

With the same logic described in the update function of the Mover entity, modify the angle with a variable describing its velocity and a variable describing its acceleration. The only difference is that the angle is represented by a single variable, and not a vector.

self.angle = self.angle + self.angularVelocity
self.angularVelocity = self.angularVelocity + self.angularAcceleration

To set the acceleration, similarly to the previous examples, use a hard-coded measure or consider the surrounding environment, the forces involved with actual formulas.

Please note:

Trigonometry

As prefaced at the top of the section, trigonometry relates to the study of the angles and sides of right triangles. In this light, the mnemonic device sohcahtoa is useful to remember the following formulae:

math.sin(theta) = opposite / hypothenuse
math.cos(theta) = adjacent / hypothenuse
math.tan(theta) = opposite / adjacent

To find the angle, the sine, cosine and tangent are used in their inverse form. For instance and knowing the sides of the right triangle, the angle can be computed as:

theta = math.atan(opposite / adjacent)

Please note:

Polar coordinats

While Love2D, similarly to the Processing library, renders elements in an (x, y) plane, with cartesian coordinates, it is useful to model the simulation in polar coordinates, considering a distance and angle (r, theta).

The trigonometry introduced in the previous section is useful to move from one set to the other.

Please note:

Amplitude and period

Oscillation, as the periodic movement between two points, can be defined in terms of amplitude and period.

Considering a sine or cosine function, the period is math.pi * 2, and the amplitude is 1.

Starting from the amplitude and period, it is possible to describe simple harmonic motion by updating a variable as follows.

x = amplitude * math.cos(math.pi * 2 * (frameCount / period))

The value returned by math.cos doesn"t exceed the (-1, 1) range, which means the variable is assigned a value in the (-amplitude, amplitude) interval. Inside of the parenthesis, dividing the incrementing variable frameCount by the period and multiplying the value by math.pi * 2 means the function completes a cycle as the count reaches the value of the period.

One other feature which defines oscillation is its frequency. This value describes the number of cycles per time unit, and is the inverse of the period. If an oscillation has a period of math.pi * 2, it completes a cycle in math.pi * 2, and has a frequency of math.pi * 2 / frame; it covers a certain distance in the span of a single frame.

Oscillation with angular velocity

Instead of mapping a variable according to amplitude and period, the idea is to consider an incrementing variable in the trigonometric functions for the sine or cosine distribution.

x = amplitude * math.cos(angle)

In this manner the simulation re-introduces the concept of angular velocity, and ultimately angular acceleration.

angle = angle + angularVelocity

It is still possible to define the period, as the amount of time it takes for angle to reach math.pi * 2.

Please note:

Applying trigonometry

With the knowledge accumulated in the chapter, and the chapters before it, the goal is to here apply the concepts in practical simulations.

Waves

In order to create a wave, all that is necessary is to increment the angle for every entity, assigning the horizontal coordinate to a fraction of the total width and the vertical coordinate to the sine or cosine of the angle. The amplitude remains relevant, in describing the height of the line.

Please note:

Pendulum

In the simulation, a pendulum is composed of a pivot, an arm and a bob.

 x  <-- pivot
  \
   \  <-- arm
    \
     o  <-- bob

The idea is to have the pivot function as the point around which the bob rotates, at a distance given by the inflexible arm.

With this structure, the pendulum is subject to a force of gravity, pulling the bob downwards. The fact that the bob is attached to the immovable pivot, however, means that the force of gravity doesn"t affect the round shape only in its y dimension.

 x
  \
   \
    \
     o
  ->/ |
  F/  |
  /90°|
  \   |
   \  |
    \a|
     \| gravity

Considering a right triangle whose hypothenuse describes the force of gravity, it is possible to decompose the vertical force in two different segments, and find the force F according to the angle a and the trigonometric functions introduced in the previous sections.

     x
     |a\
     |  \
     |   \
          o
opposite/ |
       /  |
      /   |
      \   | hypothenuse
       \  |
        \a|
         \|

The angle is updated with the same logic introduced in the force chapter.

self.angle = self.angle + self.angularVelocity * dt * UPDATE_SPEED
self.angularVelocity = self.angularVelocity + self.angularAcceleration * dt * UPDATE_SPEED

Instead of a vector, here is the angle to be updated by the velocity, and the velocity by the acceleration.

Velocity and acceleration are initialized to 0. The acceleration value is then set with the trigonometric functions prefaced in this very chapter. The relevant function is the sine function, knowing the angle and the hypothenuse, and realizing the opposite segment describes the desired force (remember the soh in sohcahtoa).

self.angularAcceleration = math.sin(self.angle) * GRAVITY * dt * -1

The value is multiplied by -1 since in Love2D, the coordinate system works left to right, top to bottom.

To create a more realistic simulation, the acceleration is also scaled according to the length of the arm.

self.angularAcceleration = math.sin(self.angle) * GRAVITY / ARM_LENGTH * dt * -1

Finally, and to have the pendulum slowly reduce its oscillation, the angular velocity is multiplied by a value slightly smaller than 1.

self.angularVelocity = self.angularVelocity * 0.995

Please note:

Spring

Instead of an inflexible arm, like the one introduced with the simulation of the pendulum, the idea is to have an elastic arm. Instead of using trigonometric functions to then describe the position of the bob and the length of the arm, the idea is to consider the force applied by the spring. In this manner, the influence of the spring can be paired with other forces, like gravity or wind.

Starting with Hooke"s law, the force applied by a spring is directly proportional to the extension of the spring. The farther the bob stretches the spring from its rest length, the greater the force.

->       ->
F = -k * x

k describes a constant to describe how elastic is the spring. It scales the vector describing the displacement to have a stronger or weaker force.

x describes the displacement from the spring"s rest length.

A spring tends to a state of equilibrium, described by the rest length. As the bob stretches the arm, the displacement causes an opposite force toward the original length.

    x
    |
    |
    | restLength
    |
  ->|
  F | currentLength
    o

Coming back to the formula, and with the logic introduced in the force chapter, it is necessary to evaluate the force"s magnitude and direction. In terms of magnitude, it is possible to directly use the constant k. In terms of direction, the vector x is obtained by comparing the current length of the arm against the rest length.

In code, the idea is to ultimately apply the force similarly to the force of gravity introduced in earlier demos.

bob:applyForce(spring)

The vector describing the force is then calculated as follows:

04 - Particle System

Starting from a single particle, the idea is to manage multiple entities, in concert. A system is useful to simulate complex phenomena, like fire, smoke, flocks of birds.

Particle

Picking up from the logic introduced with the Mover entity in the forces chapter, a single particle is built with three vectors: position, velocity and acceleration. The idea is to update the position with the velocity, and the velocity with the acceleration.

this = {
  ["position"] = position,
  ["velocity"] = velocity,
  ["acceleration"] = acceleration,
}

On top of these vectors, each particle is attributed a lifespan. This value is useful to determine when a particle dies off. In the context of a particle system, there is usually an emitter, producing particles or a stream of particles; as new and new particles are created, it is necessary to remove existing ones.

this = {
  ["lifespan"] = 1
}

In a first implementation, the lifespan is connected to the alpha channel of the particle, and decreased with every iteration.

function Particle:update(dt)
  self.lifespan = self.lifespan - dt * UPDATE_SPEED
end

function Particle:render()
  love.graphics.setColor(0.11, 0.11, 0.11, self.lifespan)
  -- draw circle
end

An additional method on the Particle entity finally describes whether the particle is dead or not, by returning true or false on the basis of the lifespan value.

function Particle:isDead()
  return self.lifespan == 0
end```

Please note:

Particles

From a single particle, the idea is to produce a new entity with every iteration of the love.update function, and remove particles when they eventually die; this last part is implemented with the :isDead method.

When removing a particle inside of a loop considering the collection, it is important to note that the collection is updated by translating the items to the left. Iterating through the table in order, the risk is to therefore skip a particle. One immediate way to fix this is to loop through the collection backwards.

for i = #particles, 1, -1 do
  particles[i]:update(dt)
  if particles[i]:isDead() then
    table.remove(particles, i)
  end
end

Please note:

Particle system

The collection of particles is stored and managed in a ParticleSystem entity. The system is initialized with an x and y coordinate, to spawn particles from a specific point of origin.

Please note::

Particle systems

Building on top of a particle system as a collection of particles, it is possible to build a collection of collections, of particle systems. This is beyond the scope of the chapter, but useful to describe a simulation in multiple layers of complexity.

Please note:

Inheritance and polymorphism

In the scope of object oriented programming, inheritance is useful to create a system in which entities pick up and expand the logic introduced by other entities. Case in point, SquareParticle and CircleParticle can inherit the behavior of the Particle entity, and render the particle with a different shape, respectively, a square and a circle.

Polymorphism, as introduced in the book, refers to how it is possible to have different classes like Dog and Cat are able to retain multiple types, multiple forms. A Dog class inheriting from an Animal class is therefore both a dog and an animal. This is useful to have a collection like an array or array list with a singular type, storing different types of animals and call the same function on every instance. What the function does, then, depends on the implementation in the different classes.

Please note:

Forces

The idea is to re-introduce the concepts of earlier chapters, but in the context of a particle system. The section is also useful to describe the structure of the simulation, and how certain methods are defined on the system, as a whole, and other methods on the particles, individually.

Apply force

One immediate way to influence the particle system is by applying a force to each and every particle. The logic is structured in two steps:

The function is defined on the Particle entity, similar to Mover in the dedicated chapter.

function Particle:applyForce(force)
  self.acceleration:add(LVector:divide(force, self.mass))
end

Please note:

Repel

Another way to affect the particles in the particle system is to create a separate structure, like a Repeller entity.

In love.update(dt), the system receives the repeller as an argument.

function love.update(dt)
  particleSystem:applyRepeller(repeller)
end

The dedicated function, then, loops through the set of particles to create the appropriate force of each entity.

function ParticleSystem:applyRepeller(repeller)
  for i, particle in ipairs(self.particles) do
    local force = repeller:repel(particle)
    particle:applyForce(force)
  end
end

Please note:

Blend mode

Particles are drawn with regular shapes, like circles or rectangles as done so far, or with images.

function Particle:render()
  love.graphics.setColor(1, 1, 1, self.lifespan)
  love.graphics.draw(img, self.position.x, self.position.y, 0, 1, 1, img:getWidth() / 2, img:getHeight() / 2)
end

Images are also useful to introduce blend modes, modifying the default pixel value of overlapping entities. With the add option, the rgb components are added to make the image brighter and eventually white. This is useful for instance to simulate fire particles.

function ParticleSystem:render()
  love.graphics.setBlendMode("add")
  -- draw particles
end

Please note:

05 - Physics Libraries

With vectors, forces, trigonometry it is possible to simulate a first environment with rudimentary physics. It is possible to refine the simulation considering more complex natural phenomena, but one alternative comes in the form of physics libraries. Here you find code by other developers already considering the issue of simulating life, simulating nature. A physics engine provides a level of complexity only grasped in the previous sections. The price is that you need to learn about the library, its requirements and also limitations.

Please note:

Box2D

Fundamentally, a simulation with Box2D works in two steps: set up and update. In the setup phase, you initialize the world, and populate the environment with however many entities are necessary. In the update phase, Box2D considers the underlying physics to update the world as necessary; there is no need to consider the position, velocity, acceleration and forces of the individual entities.

Box2D considers all the underlying physics, but it is however necessary to set up the world with the procedure and syntax prescribed by the library.

Core elements

A Box2D simulation starts with a world. This is where the simulation defines the features of the environment, like its gravity.

GRAVITY = 20
function love.load()
  world = love.physics.newWorld(0, GRAVITY)
end

Box2D works with meters, kilograms and real-world units. Since Love2D works with pixels instead, it is useful to adapt the measures with the setMeter function.

GRAVITY = 20
GRAVITY_METER = 9.81
function love.load()
  love.physics.setMeter(METER)
  world = love.physics.newWorld(0, GRAVITY * GRAVITY_METER)
end

Once initialized, the world is updated with the update() function, considering every single entity included in the simulation.

function love.update(dt)
  world:update(dt)
end

To populate the world, each entity needs three parts: a body, a shape and a fixture.

This is enough to have the world consider and update an object. To provide a visual then, Love2D provides different methods to retrieve the defining features of the bodies. In the context of a circle, body:getX(). body:getY() and shape:getRadius() allow to find the measures for the circle"s position and radius.

Particles

The demo works to create two entities dedicated to different shapes, a circle and a rectangle, and to populate a table with multiple copies of each. Adding multiple objects is also important to stress the importance of removing entities when they are no longer necessary. Removing items involves two steps:

Fixed

A previous section introduced how bodies have different types. By default an object is static, but it"s possible to modify this value already in the declaration of the body. This is what the previous demo achieved in love.physics.newBody().

love.physics.newBody(world, x, y, "dynamic")

Shortly, a body can be static, fixed in the world and not subject to its forces, dynamic, reacting to the world"s gravity, forces, and collisions, kinematic, not subject to forces, but manually moved through its velocity. Consider for instance a platform (fixed), a ball (dynamic) or a character directly controlled by the player (kinematic).

Please note:

Curvy boundary

A curved surface is introduced with a ChainShape. This particular shape accepts as argument a series of points, which are then connected to make up the object. Using a particular distribution or a trigonometric function, the effect is that the points produce the desired visual.

Please note:

Complex shapes

There are at least two different approaches to building complex shapes:

  1. use a PolygonShape, detailing the vertices of the desired outline;

  2. fix multiple shapes to the same body.

The second approach is the topic of the demo, and the reason the ComplexShape entity actually introduces two rectangle shapes.

Joint

A joint creates a connection between multiple bodies. There are different types, each with its own usefulness and defining features. For instance:

Please note:

Forces

Re-introducing the topic from previous chapters, it is possible to affect a Box2D simulation applying a force directly on a body.

Body:applyForce(fx, fy)

Please note:

Collision events

Often, it is helpful to react to a collision between bodies. Box2D provides an interface to execute some code in the lifecycle of a collision, by referring to a callcback function on the world.

world:setCallbacks(beginContact)

:setCallbacks actually accepts up to four arguments, to consider when a contact begins, ends, and when a collision is about to be resolved or has just been resolved.

world:setCallbacks(beginContact, endContact, preSolve, postSolve)

The functions must be defined in the code, and receive a series of arguments describing the collision and the objects involved. More accurately, beginContact receives the fixtures of the objects involved, and a table describing the collision.

function beginContact(f1, f2, collision)

end

The collision occurs between two objects coming into contact with each other. It is here extremely useful to know which objects are however involved. To this end, a fixture can describe a label with a userdata field.

fixture:setUserData("attractor")

The label is then evaluated in the body of the beginContact function.

Please note:

06 - Autonomous Agents

With the chapter the idea is to include entities capable of moving on their own, on the basis of a desire. These entities share three defining features:

The goal is to run a simulation without a pre-ordained structure, and analyze the interaction between the independent entities.

The demos are inspired by the book which itself cites as inspiration the work of Craig Reynolds on algorithmic steering behaviors.

Agency

In the folder you find several examples to illustrate the concept of desire. An entity might desire to move toward a target, or away from an obstacle; this wanting is materialized in a function applying a force on the basis of the surrounding environment.

Steering

In the steering example, a Vehicle entity is initilized with a structure eerily similar to the Mover or Particle entities introduced in previous sections. Every entity has a position, velocity and acceleration. Motion is however expressed in three layers:

  1. action selection: select an action on the basis of a goal or set of goals; for instance, compute the desired velocity as the difference between the position of the vehicle and target

    local desiredVelocity = LVector:subtract(target.position, self.position)

  2. steering: formulate a force to materialize the action; for instance, generate a force considering the desired velocity against the current velocity

    local steeringForce = LVector:subtract(desiredVelocity, self.velocity)

  3. locomotion: actually move the vehicle; for instance, apply a force to modify the entity's acceleration and velocity

    self:applyForce(steeringForce)

In the demo, the logic steering the vehicle toward the target is described in the steer() method. In the body of the function, the entity computes the desired velocity and steering force, but refines the movement with two additional variables:

Please note:

Arriving

The idea is to have the Vehicle entity steer toward the target, but then slow down as it gets closer and closer to the target's position.

The function evaluates the desired velocity, but also the distance of the vector.

local desiredVelocity = LVector:subtract(target.position, self.position)
local distance = desiredVelocity:getMagnitude()

The velocity is then multiplied by a value proportional to the actual distance. The smaller the distance, the slower the force. This however, only when the vehicle is in the range of the target.

if distance < RADIUS_SLOWDOWN then
  -- from [0, RADIUS_SLOWDOWN] to [0, maxSpeed]
  local speed = map(distance, 0, RADIUS_SLOWDOWN, 0, self.maxSpeed)
  desiredVelocity:multiply(speed)
end

Please note:

Pursuing

The vehicle moves toward the target considering both its position and velocity. In order to achieve this effect, the target is initialized with a velocity vector, and this vector is modified to have the entity move toward the mouse. The same vector is then incorporated in the pursue function of the vehicle in order to modify the desired velocity.

Please note:

Bouncing

The vehicle moves in the window at a constant speed, and changes this behavior in order to respect arbitrary boundaries. This is achieved by applying an force opposite to the desired velocity if the desired velocity would move the entity outside of the boundaries.

Please note:

Jittering

Following a suggestion included in the book, the vehicle moves toward a point in the vicinity of the target. By choosing a point around the target, signalled by a small circle the vehicle moves unpredictably, and yet pursuing the target.

Please note:

Flow field

With a flow field the window is divided in a certain number of columns and rows. In this grid, the cells describe a velocity, which is then picked up by the vehicle as it navigates the environment.

Each cell is attributed an angle, and there are multiple demos which differ in how this angle is computed, as well as a force vector. To compute this vector, the cosine and sine functions identify where the segment should start and end.

local x1 = math.cos(angle + math.pi)
local y1 = math.sin(angle + math.pi)
local x2 = math.cos(angle)
local y2 = math.sin(angle)

Incrementing the angle by math.pi allows to find the origin of the segment, half a rotation from the destination. This structure is useful to draw a line without the translate and rotate functions.

From this setup, all the vehicle needs is to find a cell, and apply a force matching the vector.

As mentioned, there are multiple demos which change how the angle is computed:

Please note:

Math

The Math folder introduces a few concepts essential for future demos.

Dot product

The dot product allows to compute the angle between two vectors. It is helpful to build a scalar projection, and ultimately important to introduce a simulation in which vehicles follow a given path. Defined as the multiplication of two vectors, the product is computed as follows:

By computing the product with the first formula, it is possible to find the angle solving the second formula for theta.

Please note:

Scalar projection

The goal is to find a point on a vector according to another vector, a projection.

->
a /|
 / |
/  |
---o---
   ->
   b

Cast a line from the vector a to the vector b so that the line creates a 90 degrees angle. With this line then, the dot product allows to compute the point making up the projection through the angle.

h /|
 / |
/t |
---o
 ->
 x

In terms of math, there are a couple of steps involved:

This means that ultimately, the projection is computed considering the dot product of the normalized b vector.

Please note:

Path following

The dot product and scalar projection are useful as building block for path following, another behavior studied by Craig Reynolds. The idea is to have a vehicle move in the window and follow the trajectory described by a path.

In the folder there are two demos: Straight and Segments, to show how a vehicle first follows a single straight line and then a series of connected segments.

Straight

A Path entity includes two vectors for each line, describing where the line should start and end. It also describes a radius to have the vehicle move toward the line with some margin.

In the Vehicle entity then, the follow function receives the path and modifies the velocity of the vehicle with the following logic:

Please note:

Segments

The idea is to have a path described not by two vectors, by a series of segments, each with a start and end point.

local segments = {}
for i = 1, POINTS_PATH do
-- create segment
end

To have the segments connected to one another, Path:new keeps a reference to the y coordinate where each segment should end.

local previousY = math.random(HEIGHT_MIN, HEIGHT_MAX)
  for i = 1, POINTS_PATH do
    local x1 = xStart + xIncrement * (i - 1)
    local x2 = xStart + xIncrement * i
    local y1 = previousY
    local y2 = math.random(HEIGHT_MIN, HEIGHT_MAX)
    previousY = y2
  end
end

From this setup, Path is equipped with a segments field, storing the desired x and y coordinates. The Vehicle entity then needs to loop through the collection to evaluate the projection on each and every line.

function Vehicle:follow(path)
  -- knowing desiredLocation
  for i, segment in ipairs(path.segments) do
    -- compute projection
  end
end

The idea is to consider here the closest projection which belongs to the actual path. Both conditions are necessary to avoid moving the entity towards the wrong segment.

To check if the projection belongs to the path, it is enough to check if the x coordinate falls between the beginning and the end of the segment.

if projection.x > segment.start.x and projection.x < segment.finish.x then
  -- consider distance
end

In order to consider the closest projection, then, the distance is evaluated against a variable intialized with a large value.

local recordDistance = math.huge
local recordProjection = nil

for i, segment in ipairs(path.segments) do
  -- evaluate distance
  if distance > RADIUS_PATH and distance < recordDistance then
    recordDistance = distance
    recordProjection = projection
  end
end

The projection is stored in yet another variable, in the same conditional, and is ultimately used to change the velocity of the vehicle.

if recordProjection then
  self:seek(recordProjection)
end

Please note:

Group behaviors

With a group behavior the goal is to model the movement of Vehicle entities relative to other entities.

Alignment

Align the movement of a vehicle relative to other the entities. In the align method, each entity receives the collection describing the instances.

function Vehicle:align(vehicles)
end

From this collection, the idea is to tally every vector describing the velocity.

local velocity = LVector:new(0, 0)
for i, vehicle in ipairs(vehicles) do
  velocity:add(vehicle.velocity)
end

Every vector except the one describing the entity itself. To this end, the new method is modified to have each entity distinguished with an id.

local this = {
  -- previous attributes
  ["id"] = math.random()
}

The velocity is then added only if the id do not match.

if vehicle.id ~= self.id then
  velocity:add(vehicle.velocity)
end

This has the desired effect of aligning all the entities. However, as explained in the book, it is useful to limit the area of alignment to those vehicles falling in an arbitrary radius. Here, it is helpful to consider the distance between the vectors.

if distance < DISTANCE_VEHICLE and vehicle.id ~= self.id then
  velocity:add(vehicle.velocity)
end

By increasing/decreasing VEHICLE_DISTANCE, the demo shows how a vehicle aligns itself with more/less neighbors, and form a larger/smaller cluster.

Please note:

Separation

Instead of aligning vehicles together, the entities are separated by applying a force away from the surrounding neighbors.

The code is similar to the alignment demo, in that the function loops through the collection of vehicles and considers only those vehicles closer than a given range. However, once a vehicle is within range, the vector is computed by subtracting the entities in reverse order.

local force = LVector:subtract(self.position, vehicle.position)

The idea is to describe a force away from the neighbor. Once computed, the force is added to a vector considering the cumulative vector for every neighbor.

force:divide(distance)
separationForce:add(force)

The force is also weighed by the distance, so that the closer a neighbor is, the more influence it will have on the cumulative vector.

This is enough to consider the influence of the neighbors. However, it is finally useful to scale the vector to avoid excessive values. One way to achieve this is to normalize the vector and scale the result by an arbitrary amount.

separationForce:normalize()
separationForce:multiply(MAX_SPEED)
self:applyForce(separationForce)

Please note:

Combining behaviors

The folder works to show how multiple behaviors can be paired to produce more complex simulations. The forces are computed in dedicated functions and applied in the applyBehaviors method.

function Vehicle:applyBehaviors(vehicles, target)
  local steeringForce = self:steer(target)
  local separationForce = self:separate(vehicles)

  self:applyForce(steeringForce)
  self:applyForce(separationForce)
end

By weighing the forces with different factors, the idea is to have a particular force take precedence. For instance, and in the exercise pairing the steering and separation forces, the vector pushing the vehicles away is multiplied by three. This has the net result of creating a swarm of vehicles moving toward the intended target, but giving precedence to avoid any overlap.

Cohesion

The exercises include forces already developed in previous demos, in order to steer vehicles toward a target, away from each other or to follow a path. There is however a new force in the Separation, alignment and cohesion folder, to express cohesion between the entities.

The idea of the cohesive force is to consider the surrounding neighbors and move the vehicle to the center of the group. This is achieved by considering the sum of the vectors describing the position, and computing the average.

The force is then described by the distance between the vehicle's current position and the newfound vector.

local cohesionForce = LVector:subtract(position, self.position)

07 - Cellular Automata

With a cellular automaton the book introduces a system of rules. Such a system has three foundational ingredients:

There are two particular examples in Wolfram's elementary cellular automata and Conway's game of life, but to get started, the folder includes a rudimentary system in Cellular automaton.

Cellular automaton

The example creates a one-dimensional grid, where a cell has up to two neighbors, described by the previous or following unit.

A cell has a boolean value, initialized at random.

local isAlive = math.random() > 0.5

At each iteration, the idea is to modify the state so that a cell is alive only if one of its neighbors is alive.

local isAlive = aliveNeighbors == 1

Wolfram

Wolfram systems are described by one-dimensional cellular automata, where a ruleset dictates the state of the cell as a function of the state of the same cell and its surrounding neighbors.

Each cell is initialized with a binary value, so that the generation is described by a sequence of 0s and 1s. At each iteration then, the idea is to consider the neighboring cell to form a string, like 001 or 010, and use the string to describe the state from a given set. For instance, a ruleset might be initialized as follows:

local ruleset = {
  ["111"] = 1,
  ["110"] = 0,
  ["101"] = 1,
  ["100"] = 0,
  ["011"] = 1,
  ["010"] = 0,
  ["001"] = 1,
  ["000"] = 1,
}

As a sequence of neighbors describes 101, the new cell receives the matching state of 1.

The specific sequence describes the rule, in decimal representation. For instance, ruleset 10101011 describes rule 171.

The outcome can be described in one of four categories:

Please note: the Wolfram folder contains a series of demos, each with its own purpose:

Game of life

The game of life provides two dimensional cellular automata. More than just an exercise, its works as a starting point for simulations inspired by nature, and specifically inspired by a system of biological reproduction. With a set of limited rules, the idea is to produce pattern and unpredictability similarly to the systems introduced in the Wolfram section. The systems ultimately settle to create uniformity or repeating the same pattern; oscillating as it were.

The two dimensional system is described by a two dimensional collection. Each iteration creates not a new generation, but a new frame of the system.

Each cell is initialized with a boolean at random.

grid[column][row] = math.random() > 0.5 and 1 or 0

With each iteration, then, the idea is to consider the available neighbors and change the state not on the basis of a ruleset, but considering the features of the neighborhood, and speficially the number of neighboring cells which are alive. The rules modify the grid as follows:

Please note: similarly to Wolfram, the folder dedicated to the game of life includes a series of demos, each with its own purpose:

08 - Fractals

Fractals are defined as geometric shapes that can be split into parts, each of which resembles the whole. These shapes share a few common features:

Recursive draw

To draw a fractal it is first necessary to introduce the concept of recursion, the repeated application of a rule to successive results. In practice, this concept sees a function calling itself, but with different arguments. The exemplary use case is that of a function computing the factorial.

function factorial(n)
  if n == 1 then
    return 1
  end
  return n * factorial(n - 1)
end

Trivially, the factorial of 5 is 5! and is computed as 5*4*3*2*1. It can be re-written as 5*4!, and 5*4*3! until 1. It is important to stress this last value because it describes the exit condition, where the recursion should stop. Without such a condition, the function would call itself indefinitely resulting in a stack overflow.

Acknowledging recursion, the demo shows how a draw function calls itself to draw circles at different locations and with different radii.

function draw(x, y, r)
  love.graphics.circle("line", x, y, r)
  if r > RADIUS_MIN then
    draw(x, y, r / 2)
    draw(x + r, y, r / 2)
    draw(x - r, y, r / 2)
    draw(x, y + r, r / 2)
  end
end

The function is first called with a value describing the largest circle.

draw(WINDOW_WIDTH / 2, WINDOW_HEIGHT / 2, RADIUS_MAX)

At each iteration, then, it draws more circles assuming the radius is more than an arbitrary threshold. As the draw function calls itself with different maeasures, every circle has a smaller circles to the right, to the left and to the bottom, resulting in an intriguing pattern.

Sierpinski triangle

With a function drawing a shape and recursively calling itself it is possible draw Sieprinski triangle, a particular type of fractal in which a triangle is subdivided into three, smaller triangle. The structure is based on an equilateral triangle, where the sides are all the same; to find the coordinates of this shape, it is necessary to know the height of the triangle.

local height = (side * 3 ^ 0.5) / 2

Once an equilateral triangle is drawn, the idea is to position the three, smaller triangles to subdivide the shape evenly.

draw(x - side / 4, y + height / 4, side / 2)
draw(x + side / 4, y + height / 4, side / 2)
draw(x, y - height / 4, side / 2)

Notice that the code retains an exit condition, once more to have the function eventually terminate its recursive pattern.

if side > SIDE_MIN then
  -- draw triangles
end

Cantor rule

A Cantor rule is expressed by a line recursively drawn by dividing the line in three parts, and erasing the middle third. The concept is visualized in the demo showing the different steps.

Koch line

Building on top of the example introduced with the Cantor rule, a Koch line divides a line in three parts, erases the middle portion and connects the remaining thirds with the sides of an equilateral triangle. It is here necessary to have a reference to the line(s), so that the new segments are created from known coordinates, and to this end, the code is modified to have the line(s) expressed with a dedicated entity and two vectors, detailing where the lines should start and end.

At first, the demo shows a single line spanning the width of the window. Following a mouse click, then, the idea is to have the line segmented following the rules mentioned above. Consider each step as a generation, much similarly to the CA described in the previous chapter.

Please note:

Koch snowflake

Following a suggestion from the book, the demo reiterates the concept introduced in the Koch line to produce a regular polygon where the sides are actually Koch lines. The Snoflake entity is initialized with a certain number of sides and generations, with default values describing a triangular shape after 5 iterations. Following a mouse press, the shape is re-initialized with a random number of sides in the [3, 10] range.

The demo includes a bit of math, mostly to find the coordinates of the line segments, but also to have the size of the polygon change depending on the number of sides. The idea is to have a variable RADIUS describe the radius of the circle wrapping around the polygon, and compute the length of the side and of the apothem in order to keep the shape inside the edges of the window and exactly in the center.

Tree

Trees are useful to describe fractals without perfect self-similarity. By modifying the shape with random values and probabilities the idea is to avoid a stylized, symmetric shape to find another way to emulate nature.

The folder includes multiple demos to explore the topic:

L-system

An L-system is a grammar-based sytem, a way to write strings that is adapted in the demos to map characters to particular drawing instructions. It is ultimately useful to describe production rules through strings, commanding an entity with a series of instructions.

Such a system is characterized by three defining features:

Generation after generation, the rule makes it possible to build a sentence with a variety of characters.

The folder includes a few demos to illustrate the point:

Please note:

09 - The Evolution of Code

The idea of the chapter is to have information passed from object to object in order to have the simulation evolve over time. The principles of biological evolution, and specifically darwinian evolution, are taken as inspiration to create a series of genetic algorithms and solve specific problems.

The first type of algorithm is useful to solve problems where the solution space is so vast to make a brute-force approach feasible. For instance, having a computer produce a certain sequence of letters matching an input string. The problem is trivial, but helps to demonstrate how a genetic algorithm works.

Inspired by actual biological evolution and specifically darwinian natural selection, genetic algorithms consider three key principles:

The principles are put into practice in a series of steps, and are detailed in the context of a specific problem.

Shakesperian monkey

The demo works to illustrate the brute force approach of finding a match for a four letter word. While ineffective, it also works to introduce the building blocks for the genetic algorithm:

With a four letter word and twenty-six possible characters, the odds of finding a match are already a measly 1 in 456 976, motivating a different approach.

Simplified genetic algorithm

Please note: the demo differs from the logic discussed in the book in that the program works through functions instead of dedicated entities, like Population or DNA classes. This is on purpose to have the project comparable with Shakesperian monkey.

The steps of the traditional genetic algorithm are included in the love.load and love.update functions.

And that is essentially it for the genetic algorithm. With each iteration, population includes fitter and fitter values, until a children is able to reproduce the input sentence.

The rest of the logic described in love.update is useful to:

  1. identify the word in the population with the best fit

  2. stop the iterative process when the input is reproduced

  3. store the best fit in a words collection, used to show the result in the window

Traditional genetic algorithm

The idea is to refactor the previous demo to implement the algorithm with DNA and Population entities. This allows the script to have a more general structure, one in which love.load initializes a population and love.update modifies the population with a series of functions.

population:select()
population:reproduce()

How the population is initialized, how the mating pool, children, parents, fitness value are calculated is then a matter delegated to the specific entities.

The exercise doesn't introduce concepts, but there are a couple of notable differences with respect to the previous demo:

Pool selection

The project updates the traditional genetic algorithm in the way it selects two parents elements. Instead of populating a selection table with a number of copies proportional to the fitness ratio, the script picks a dna at random and accepts or rejects the element according to a probability. The probability is mapped to the maximum fitness ratio, so that the greater the ratio, the more likely it will be for the element to be picked.

local parent1 = self:select(maxFitnessRatio)
local parent2 = self:select(maxFitnessRatio)

The select function continues to pick a dna until a suitable parent is found.

while true do
  local dna = self.population[math.random(#self.population)]
  local fitnessRatio = dna:getFitnessRatio(self.target)
  local probability = math.random() * maxFitnessRatio
  if probability < fitnessRatio then
    return dna
  end
end

The fitness ratio is also squared to preserve the intention of the previous demo, that is increase the importance of higher values.

local fitnessRatio = dna:getFitnessRatio(self.target)
fitnessRatio = fitnessRatio ^ 2

Evolutionary flow field

With this demo I intend to apply the logic of the traditional genetic algorithm to the demo created in chapter 6 to moves vehicles according to a random flow field (the-nature-of-code\06 Autonomous Agents\Flow field\Random).

The demo is modified in increments and as follows:

The project is useful to stress the difference between genotype (genetic material, DNA) and phenotype (visual representation, Vehicle) as introduced in the beginning of the chapter. In future demos the distinction is reiterated showing the genetic material with faces or again particles.

Interactive selection

In the exercise the genetic algorithm is tailored to have a population evolve as a result of user input. The demo shows a collection of faces, changing the size, color and position of a few prominent features. When the mouse hovers upon a specific face, the fitness value is increased by an arbitrary amount. When a new generation is created, the population consider the elements with higher fitness values so that the face becomes a mixture of the preferred patterns.

Please note:

Ecosystem simulation

The project sets up a genetic algorithm to simulate an ecosystem evolving over time. Instead of having a clear stopping point, where the population evolves abruptly from generation to generation, the goal is to have a series of entities live on their own.

Each entity, or bloop, has a lifespan, counting down to 0, at which point the particle dies. The lifespan is mapped to the opacity of the figure, so that the members of the population slowly fade out, but the value is increased when an entity collides with a pellet of food. Randomly and while still alive, each bloop has also a possibility to spawn a new particle, to which it passes its own genetic material in terms of size.

Variety is included in form of the size of the bloop. The size is inversely proportional to the speed, so that larger particles move slower than smaller ones. Eventually, the simulation shows a preference for a particular type of particle, preferring a smaller, faster entity than a larger, slower one.

Please note:

Smart rocket

The goal is to have a population of rockets navigate the environment toward an arbitrary target. The movement is influenced by a series of vectors describing the acceleration of the rocket at different frames. As the generation comes to an end, the fitness is computed on the basis of the distance between rocket and target, so that eventually, the population learns to follow a similar, efficient trajectory.

From this starting point, the demo is updated to have a more complex simulation, with one or more obstacles interrupting the movement of the entities. By default, the simulation includes a single obstacle, in the right half of the window, but following mouse input, it is possible to complicate the environment with additional platforms.

The fitness value takes account of the distance, but also the obstacles and target. In this light, hitting an obstacle, or the edges of the window, results in a greatly diminished fitness.

if self.hasCollided then
  return fitness * 0.1
end

Hitting the target, on the other hand, increases the fitness considering the number of frames it takes the rocket to reach the target.

if self.hasReached then
  return fitness * FRAMES / self.frames
end

self.frames is a counter variable for each rocket, incrementing each time the rocket moves. Dividing by this value means the program has a preference for entities reaching the target in less frames. Multiplying for this value, however, the preference is set for slower rockets. This could be preferable in a scenario where the goal is to get closer to the target, without crashing in it.

if self.hasReached then
  return fitness * self.frames / FRAMES
end

Please note:

10 - Neural Networks

A neural network is introduced as a complex, adaptive system; a system which learns over time to solve a specific problem.

Perceptron

A first example of neural network is a perceptron, a system with one or more inputs, a processor and a single output.

input 1
       -> processor -> output
input 2
input n

The perceptron is a type of network described as feed forward, meaning that learning happens by feeding the system a series of input values and measuring the output against a known result. As the result in known in advance, the network is also an example of supervised learning, meaning the system is trained to provide an acceptable answer on the basis of training data.

The network associates a weight to each input.

input1 x weight1
                  -> processor -> output
input2 x weight2

These weights are then adjusted as the network encounters a mistake, an error between guessed and known output.

In the specific demo, the structure is used to solve a classification problem, where the window is divided in two areas and the task is to categorize an input point to a specific region.

1 1/
1 / -1
 / -1
/

The perceptron receives as input an x and y coordinate and returns 1 or -1 according to whether the point is above or below a diagonal. The problem could be solved with linear algebra, but with a network, it is approached with a series of points describing the training data. The perceptron guesses the output starting from a random set of weights, and tunes these parameters to respect the training data.

A guess works considering the inputs, weights, and an activation function. The goal of an activation function is to produce the output in a desired form (in this instance 1 or -1). In the specific demo, the function considers the sign of the weighted sum for all the input values (sum of the values multiplied by the respective weights).

function Perceptron:activate(n)
  return n >= 0 and 1 or -1
end

It is important to note that evaluating the sign is but one of the possible activation functions. It is useful here since it provides a binary, clear way to return 1 or -1.

Training is performed with a series of points which are already assigned the correct output value.

function Point:new()
  local this = {
    ["label"] = x > y and 1 or -1
  }
end

As it receives a point, the perceptron evaluates the error as the difference between guess and answer.

local guess = self:guess(inputs)
local error = target - guess

The combinations are described by this short table. From this error, the weights are adjusted by multiplying the error for the input value.

for i = 1, #self.weights do
  self.weights[i] = self.weights[i] + error * inputs[i]
end

This is similar to the example describing the steering behavior for autonomous agents. The difference between desired and actual velocity provides the error, which is then used to steer the agent toward the desired target.

The correction is multiplied by a learning rate, to have the network assimilate the error, learn from its mistake, by moving the weights toward the desired values.

self.weights[i] = self.weights[i] + error * inputs[i] * LEARNING_RATE

This is again similar to the steering example. How much the error modifies the velocity describes how rapidly the vehicle changes its direction, and if it approaches the desired trajectory or overshoots the same value.

In the demo, the perceptron is highlighted with a series of points, colored with a fill or stroke according to the label 1 and -1. The perceptron makes a series of guesses and highlights them with a green or red outline. At first the guesses are wrong, but by pressing the mouse cursor, the perceptron is trained with more and more points, improving the decision with each additional datum.

Perceptron line

The first demo implementing the perceptron is modified to have the network consider a more general classification problem, where a point is assigned -1 or 1 on the basis of the equation of a line

y = slope * x + coefficient

The demo is updated to also consider a cartesian coordinate system, where the inputs are in the [-1, 1] range and the values are mapped to the pixel coordinates with a map function.

Given the new coordinate system, it is finally necessary to introduce a third input in a bias

function Point:new()
  local this = {
    ["inputs"] = {x, y, 1},
  }
end

The value is necessary to consider a valid guess for the (0, 0) origin.

Matrix

In order to develop a multilayer neural network, it is helpful to have a small library to manage data in matrix format, that is data in rows and columns. Matrix.lua implements matrix and matrix math through a series of functions.

The functions can be split in two main categories: those which modify the instance matrix, and those which return a new instance altogether. In the first category, it is possible to modify the instance with:

In the second category, it is possible to add and multiply matrices. What changes is that the input matrixes are not modified.

In terms of object oriented programming, it is necessary to decipher if the input of the add or multiply function is a matrix, since the logic is different from when the input is a scalar. To achieve this with lua lua, it would be enough to check the metatable of the input.

if getmetatable(input) == self

end

The isInstance function slightly refines this approach by returning a boolean, and considering the metatables of the input table in a loop.

function isInstance(instance, class)
  while instance do
    instance = getmetatable(instance)
    if instance == getmetatable(class) then
      return true
    end
  end
  return false
end

Since getmetatable returns the metatable of the instance, and nil when a metatable is no longer available, it is possible to cycle through all the "parent metatables", and to check multiple layers.