Gain insights into utilizing NumPy for data manipulation and analytics. Learn about implementing concepts in both Python and NumPy through coding challenges and quizzes.

python.tar.gz

python3.6

Live3.6

If you're looking to grow your career in machine learning or data science in this day and age, adding a powerful library to your skill set is an important place to start. In that vein, Python has become one of the most widely used tools in the industry for serious data analytics, and NumPy is probably the most widely used data analytics library. With NumPy, you can manipulate data involving multi-dimensional arrays and matrices (think linear algebra).

Join us as we venture into the vast world of NumPy in this comprehensive course. Each lesson dive into the actual implementation of concepts in both pure Python and then NumPy, exploring how NumPy vectorization compares to traditional Python that uses a procedural and object-oriented approach.

Practice and test yourself along the way with in-browser coding challenges, quizzes, and more.

This course is intended for users who are already familiar with intermediate level Python.

From Python to Numpy

# Boid Class Implementation (Python)
Since each boid is an autonomous entity with several properties such as position and velocity, it seems natural to start by writing a `Boid` class:

import math
import random
from vec2 import vec2

class Boid:
    def __init__(self, x=0, y=0):
        self.position = vec2(x, y)
        angle = random.uniform(0, 2*math.pi)
        self.velocity = vec2(math.cos(angle), math.sin(angle))
        self.acceleration = vec2(0, 0)

The `vec2` object is a very simple class that handles all common vector operations with 2 components. It will save us some writing in the main `Boid` class. Note that there are some vector packages in the Python Package Index, but that would be an overkill for such a simple example.

Boid is a difficult case for regular Python because a boid has interaction with local neighbors. Since the boids are moving so, at every step, we need to calculate the distance of one boid with every other boid that comes within the interaction radius and then sort those distances. The prototypical way of writing the three rules is thus something like:

def separation(self, boids):
  count = 0
  for other in boids:
    d = (self.position - other.position).length()
    if 0 < d < desired_separation:
      count += 1
      # ...
  if count > 0:
    # ...

def alignment(self, boids): # ...
def cohesion(self, boids): # ...

To complete the picture, we can also create a `Flock` object:

class Flock:
  def __init__(self, count=150):
    self.boids = []
    for i in range(count):
      boid = Boid()
      self.boids.append(boid)

  def run(self):
    for boid in self.boids:
      boid.run(self.boids)

# Complete Solution

Merging all the logic from above, given below is the complete implementation of Boids. 

import math
from collections import namedtuple


def struct(name, members):
    cls = namedtuple(name, members)
    cls.__repr__ = lambda self: "%s(%s)" % (name, ','.join(str(s) for s in self))
    return cls


class vec2(struct('vec2', ('x', 'y'))):

    def __add__(self, other):
        if isinstance(other, vec2):
            return vec2(self.x+other.x, self.y+other.y)
        return vec2(self.x+other, self.y+other)

    def __sub__(self, other):
        if isinstance(other, vec2):
            return vec2(self.x-other.x, self.y-other.y)
        return vec2(self.x-other, self.y-other)

    def __mul__(self, other):
        if isinstance(other, vec2):
            return vec2(self.x*other.x, self.y*other.y)
        return vec2(self.x*other, self.y*other)

    def __truediv__(self, other):
        if isinstance(other, vec2):
            return vec2(self.x/other.x, self.y/other.y)
        return vec2(self.x/other, self.y/other)

    def length(self):
        return math.hypot(self.x, self.y)

    def normalized(self):
        length = self.length()
        if not length:
            length = 1.0
        return vec2(self.x/length, self.y/length)

    def limited(self, maxlength=1.0):
        length = self.length()
        if length > maxlength:
            return vec2(maxlength*self.x/length, maxlength*self.y/length)
        return self


# -----------------------------------------------------------------------------
# From Numpy to Python
# Copyright (2017) Nicolas P. Rougier - BSD license
# More information at https://github.com/rougier/numpy-book
# -----------------------------------------------------------------------------
import math
import random
from vec2 import vec2


class Boid:
    def __init__(self, x, y):
        self.acceleration = vec2(0, 0)
        angle = random.uniform(0, 2*math.pi)
        self.velocity = vec2(math.cos(angle), math.sin(angle))
        self.position = vec2(x, y)
        self.r = 2.0
        self.max_velocity = 2
        self.max_acceleration = 0.03

    def seek(self, target):
        desired = target - self.position
        desired = desired.normalized()
        desired *= self.max_velocity
        steer = desired - self.velocity
        steer = steer.limited(self.max_acceleration)
        return steer

    # Wraparound
    def borders(self):
        x, y = self.position
        x = (x+self.width) % self.width
        y = (y+self.height) % self.height
        self.position = vec2(x,y)

    # Separation
    # Method checks for nearby boids and steers away
    def separate(self, boids):
        desired_separation = 25.0
        steer = vec2(0, 0)
        count = 0

        # For every boid in the system, check if it's too close
        for other in boids:
            d = (self.position - other.position).length()
            # If the distance is greater than 0 and less than an arbitrary
            # amount (0 when you are yourself)
            if 0 < d < desired_separation:
                # Calculate vector pointing away from neighbor
                diff = self.position - other.position
                diff = diff.normalized()
                steer += diff/d  # Weight by distance
                count += 1       # Keep track of how many

        # Average - divide by how many
        if count > 0:
            steer /= count

        # As long as the vector is greater than 0
        if steer.length() > 0:
            # Implement Reynolds: Steering = Desired - Velocity
            steer = steer.normalized()
            steer *= self.max_velocity
            steer -= self.velocity
            steer = steer.limited(self.max_acceleration)

        return steer

    # Alignment
    # For every nearby boid in the system, calculate the average velocity
    def align(self, boids):
        neighbor_dist = 50
        sum = vec2(0, 0)
        count = 0
        for other in boids:
            d = (self.position - other.position).length()
            if 0 < d < neighbor_dist:
                sum += other.velocity
                count += 1

        if count > 0:
            sum /= count
            # Implement Reynolds: Steering = Desired - Velocity
            sum = sum.normalized()
            sum *= self.max_velocity
            steer = sum - self.velocity
            steer = steer.limited(self.max_acceleration)
            return steer
        else:
            return vec2(0, 0)

    # Cohesion
    # For the average position (i.e. center) of all nearby boids, calculate
    # steering vector towards that position
    def cohesion(self, boids):
        neighbor_dist = 50
        sum = vec2(0, 0)  # Start with empty vector to accumulate all positions
        count = 0
        for other in boids:
            d = (self.position - other.position).length()
            if 0 < d < neighbor_dist:
                sum += other.position  # Add position
                count += 1
        if count > 0:
            sum /= count
            return self.seek(sum)
        else:
            return vec2(0, 0)

    def flock(self, boids):
        sep = self.separate(boids)  # Separation
        ali = self.align(boids)  # Alignment
        coh = self.cohesion(boids)  # Cohesion

        # Arbitrarily weight these forces
        sep *= 1.5
        ali *= 1.0
        coh *= 1.0

        # Add the force vectors to acceleration
        self.acceleration += sep
        self.acceleration += ali
        self.acceleration += coh

    def update(self):
        # Update velocity
        self.velocity += self.acceleration
        # Limit speed
        self.velocity = self.velocity.limited(self.max_velocity)
        self.position += self.velocity
        # Reset acceleration to 0 each cycle
        self.acceleration = vec2(0, 0)

    def run(self, boids):
        self.flock(boids)
        self.update()
        self.borders()


class Flock:
    def __init__(self, count=150, width=640, height=360):
        self.width = width
        self.height = height
        self.boids = []
        for i in range(count):
            boid = Boid(width/2, height/2)
            boid.width = width
            boid.height = height
            self.boids.append(boid)

    def run(self):
        for boid in self.boids:
            # Passing the entire list of boids to each boid individually
            boid.run(self.boids)

    def cohesion(self, boids):
        P = np.zeros((len(boids),2))
        for i, boid in enumerate(self.boids):
            P[i] = boid.cohesion(self.boids)
        return P

        


import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

n=50
flock = Flock(n)
P = np.zeros((n,2))

def update(*args):
    flock.run()
    for i,boid in enumerate(flock.boids):
        P[i] = boid.position
    scatter.set_offsets(P)

fig = plt.figure(figsize=(8, 5))
ax = fig.add_axes([0.0, 0.0, 1.0, 1.0], frameon=True)
scatter = ax.scatter(P[:,0], P[:,1],
                     s=30, facecolor="red", edgecolor="None", alpha=0.5)

ax.set_xlim(0,640)
ax.set_ylim(0,360)
ax.set_xticks([])
ax.set_yticks([])
Writer = animation.writers['ffmpeg']
writer = Writer(fps=15, metadata=dict(artist='Me'), bitrate=1800)
anim = animation.FuncAnimation(fig, update, interval=10)
anim.save('output/output.mp4', writer=writer)
plt.show()

# Drawbacks 
Using this approach, we can have up to 50 boids until the computation time becomes too slow for a smooth animation. As you may have guessed, we can do much better using NumPy, but let me first point out the main problem with this Python implementation.


- If you look at the code, you will certainly notice there is a lot of redundancy. More precisely, we do not exploit the fact that the *Euclidean* distance is reflexive, that is, $|x-y| = |y-x|$.

- In this naive Python implementation, each rule (function) computes $n^2$ distances while $\frac{n^2}{2}$ would be sufficient if properly cached.
- Furthermore, each rule re-computes every distance without caching the result for the other functions. In the end, we are computing $3n^2$ distances instead of $\frac{n^2}{2}$.

In this lesson we will try to implement the Boids class using the traditional Pythonic approach and analyze it in terms of efficiency. 

Introduction

Anatomy of an Array

Code Vectorization

Problem Vectorization

Custom Vectorization

Beyond NumPy

Conclusion

Coding Example: Implement the behavior of Boids (Python approach)

Boid Class Implementation (Python)