{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Boids models of swarm intelligence\n",
    "\n",
    "\n",
    "<img src='./assets/Auklet_flock_Shumagins_1986.jpg' width='50%'>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above photograph of Auklet birds in flight show the classic flocking pattern observed when large groups of birds or fish travel together. Flocking is often cited as an example of [swarm intelligence](https://en.wikipedia.org/wiki/Swarm_intelligence) and the Boids models created by [Craig Reynolds](http://www.red3d.com/cwr/boids/) (1986) is one of the most well-known computational model of such behavior. In this model, each bird is represented by an agent in the simulation (called a boid) which follows a set of local rules. By defining how each boid responds to its neighbors, large groups of boids exhibit complex, emergent behaviors.\n",
    "\n",
    "In this notebook, we will set up a boid simulation and visualize and interact with it using HoloViews. The code used here is a highly condensed version of the boids code in the 'From Python to Numpy' book by Nicolas Rougier that you can find [here](https://www.labri.fr/perso/nrougier/from-python-to-numpy/#boids). This is an excellent resource for learning how to build simulations with NumPy and for learning how exactly the boids code used in this notebook works.\n",
    "\n",
    "We start by importing HoloViews and NumPy and loading the extension:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import holoviews as hv\n",
    "from holoviews import opts\n",
    "import numpy as np\n",
    "\n",
    "hv.extension('bokeh')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Defining the simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our boids are simply points with an associated velocity that live in a 2D toroidal universe where the edges of the world wrap around. Our world has a width and a height, and our boids have a velocity ``vel`` and a position ``pos``.\n",
    "\n",
    "The following class defines the initial state of our boids simulation where we have defined a simple random initial state, giving our boids and initial randomized position and velocity:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def radarray(N):\n",
    "    \"Draw N random samples between 0 and 2pi radians\"\n",
    "    return np.random.uniform(0, 2*np.pi, N)\n",
    "\n",
    "class BoidState(object):\n",
    "    def __init__(self, N=500, width=400, height=400):\n",
    "        self.width, self.height, self.iteration = width, height, 0\n",
    "        self.vel = np.vstack([np.cos(radarray(N)),  np.sin(radarray(N))]).T\n",
    "        r = min(width, height)/2*np.random.uniform(0, 1, N)\n",
    "        self.pos = np.vstack([width/2 +  np.cos(radarray(N))*r,  \n",
    "                              height/2 + np.sin(radarray(N))*r]).T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To keep our simulation code as short as possible, we define two helper functions that we will be reusing shortly:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def count(mask, n): \n",
    "    return np.maximum(mask.sum(axis=1), 1).reshape(n, 1)\n",
    "\n",
    "def limit_acceleration(steer, n, maxacc=0.03):\n",
    "    norm = np.sqrt((steer*steer).sum(axis=1)).reshape(n, 1)\n",
    "    np.multiply(steer, maxacc/norm, out=steer, where=norm > maxacc)\n",
    "    return norm, steer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can define a highly condensed ``flock`` method on the ``Boids`` class which runs a single step of our boids flocking simulation. This code applies the following three local rules to all of our boid agents:\n",
    "\n",
    "* separation: Each boid steers to avoid crowding in its neighborhood\n",
    "* alignment: Each boid steers towards the average heading of its localized neighbors\n",
    "* cohesion: Each boid steers toward the average position (center of mass) of its localized neighbors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Boids(BoidState):\n",
    "    \n",
    "    def flock(self, min_vel=0.5, max_vel=2.0):\n",
    "        n = len(self.pos)\n",
    "        dx = np.subtract.outer(self.pos[:,0], self.pos[:,0])\n",
    "        dy = np.subtract.outer(self.pos[:,1], self.pos[:,1])\n",
    "        dist = np.hypot(dx, dy)\n",
    "        mask_1, mask_2 = (dist > 0) * (dist < 25), (dist > 0) * (dist < 50)\n",
    "        target = np.dstack((dx, dy))\n",
    "        target = np.divide(target, dist.reshape(n,n,1)**2, out=target, where=dist.reshape(n,n,1) != 0)\n",
    "        steer = (target*mask_1.reshape(n, n, 1)).sum(axis=1) / count(mask_1, n)\n",
    "        norm = np.sqrt((steer*steer).sum(axis=1)).reshape(n, 1)\n",
    "        steer = max_vel*np.divide(steer, norm, out=steer, where=norm != 0) - self.vel\n",
    "        norm, separation = limit_acceleration(steer, n)\n",
    "        target = np.dot(mask_2, self.vel)/count(mask_2, n)\n",
    "        norm = np.sqrt((target*target).sum(axis=1)).reshape(n, 1)\n",
    "        target = max_vel * np.divide(target, norm, out=target, where=norm != 0)\n",
    "        steer = target - self.vel\n",
    "        norm, alignment = limit_acceleration(steer, n)\n",
    "        target = np.dot(mask_2, self.pos)/ count(mask_2, n)\n",
    "        desired = target - self.pos\n",
    "        norm = np.sqrt((desired*desired).sum(axis=1)).reshape(n, 1)\n",
    "        desired *= max_vel / norm\n",
    "        steer = desired - self.vel\n",
    "        norm, cohesion = limit_acceleration(steer, n)\n",
    "        self.vel += 1.5 * separation + alignment + cohesion\n",
    "        norm = np.sqrt((self.vel*self.vel).sum(axis=1)).reshape(n, 1)\n",
    "        np.multiply(self.vel, max_vel/norm, out=self.vel, where=norm > max_vel)\n",
    "        np.multiply(self.vel, min_vel/norm, out=self.vel, where=norm < min_vel)\n",
    "        self.pos += self.vel + (self.width, self.height)\n",
    "        self.pos %= (self.width, self.height)\n",
    "        self.iteration += 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualizing the boid simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our simulation consists of points (boids) in 2D space that have a heading. The natural HoloViews element to visualize this data is the ``VectorField``. We start by setting some plot and style options for ``VectorField`` elements in this notebook:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "opts.defaults(opts.VectorField(xaxis=None, yaxis=None, scale=0.08, \n",
    "                               normalize_lengths=False, rescale_lengths=False))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's initialize the simulation with 500 boids:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "boids = Boids(500)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can write a simple function that takes our ``boids`` simulation and returns a ``VectorField``, labelling it with the simulation iteration number. We can now use this to visualize the randomized initial state:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def boids_vectorfield(boids, iteration=1):\n",
    "    angle = (np.arctan2(boids.vel[:, 1], boids.vel[:, 0]))\n",
    "    return hv.VectorField([boids.pos[:,0], boids.pos[:,1], \n",
    "                           angle, np.ones(boids.pos[:,0].shape)], extents=(0,0,400,400), \n",
    "                          label='Iteration: %s' % boids.iteration)\n",
    "\n",
    "boids_vectorfield(boids)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have our ``boids_vectorfield`` function, we can easily define a ``flock`` function that steps the flock simulation and returns the resulting ``VectorField``. This can be used in a ``DynamicMap`` together with the streams system as described in the Responding to Events and Custom Interactivity user guides:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from holoviews.streams import Stream\n",
    "\n",
    "def flock():\n",
    "    boids.flock()\n",
    "    return boids_vectorfield(boids)\n",
    "\n",
    "dmap = hv.DynamicMap(flock, streams=[Stream.define('Next')()])\n",
    "dmap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initially, the output above shows the result of the simulation at iteration zero. By updating the stream (which has no parameters), we can now drive our simulation forwards using the ``event`` method on our ``DynamicMap``:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dmap.periodic(0.01, timeout=60, block=True) # Run the simulation for 60 seconds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python",
   "pygments_lexer": "ipython3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}