{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# Tutorial 02: Targets and options\n\nMany useful **local averages** (see `tuto_averages`) are pre-defined and may be directly used in the simulation of new or classical models. This tutorial showcases the basic usage of the **target methods** to simulate the variants of the Vicsek model. \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## An example: variants of the Vicsek model\n\nIn its most abstract form, the Vicsek model reads: \n\n\\begin{align}\\mathrm{d}X^i_t = c_0 V^i_t \\mathrm{d} t\\end{align}\n\n\\begin{align}\\mathrm{d}V^i_t = \\sigma\\mathsf{P}(V^i_t)\\circ ( J^i_t \\mathrm{d}t + \\mathrm{d} B^i_t),\\end{align}\n\nwhere $c_0$ is the speed, $\\mathsf{P}(v)$ is the orthogonal projection on the orthogonal plane to $v$ and $\\sigma$ is the diffusion coefficient. The drift coefficient $J^i_t$ is called a **target**. In synchronous and asynchronous Vicsek models, the target is the center of the sampling distribution. A comparison of classical targets is shown below.\n\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First, some standard imports...\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import time \nimport torch\nimport pprint\nfrom matplotlib import pyplot as plt\nfrom sisyphe.models import Vicsek\nfrom sisyphe.display import display_kinetic_particles\n\nuse_cuda = torch.cuda.is_available()\ndtype = torch.cuda.FloatTensor if use_cuda else torch.FloatTensor" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The classical non-normalised Vicsek model\n\nThe target is: \n\n\\begin{align}J^i_t = \\kappa\\frac{\\sum_{j=1}^N K(|X^j_t-X^i_t|)V^j_t}{|\\sum_{j=1}^N K(|X^j_t-X^i_t|)V^j_t|}.\\end{align}\n\nThe parameters of the model... \n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "N = 100000\nL = 100 \n\npos = L*torch.rand((N,2)).type(dtype)\nvel = torch.randn(N,2).type(dtype)\nvel = vel/torch.norm(vel,dim=1).reshape((N,1))\n\nR = 3.\nc = 3.\nnu = 5.\nsigma = 1.\n\ndt = .01" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The choice of the target is implemented in the keyword argument ``variant``. For the classical normalised target, it is given by the following dictionary.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "variant = {\"name\" : \"normalised\", \"parameters\" : {}}\n\nsimu = Vicsek(\n pos = pos.detach().clone(),\n vel = vel.detach().clone(), \n v = c, \n sigma = sigma, \n nu = nu, \n interaction_radius = R,\n box_size = L,\n dt = dt,\n variant = variant,\n block_sparse_reduction = True,\n number_of_cells = 40**2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally run the simulation over 300 units of time.... \n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "frames = [0., 5., 10., 30., 42., 71., 100, 124, 161, 206, 257, 300]\n\ns = time.time()\nit, op = display_kinetic_particles(simu, frames, order=True)\ne = time.time()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Print the total simulation time and the average time per iteration. \n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "print('Total time: '+str(e-s)+' seconds')\nprint('Average time per iteration: '+str((e-s)/simu.iteration)+' seconds')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot the histogram of the angles of the directions of motion. \n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "angle = torch.atan2(simu.vel[:,1],simu.vel[:,0])\nangle = angle.cpu().numpy()\nh = plt.hist(angle, bins=1000)\nplt.xlabel(\"angle\")\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "After an initial clustering phase, the system self-organizes into a uniform flock. \n\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Non-normalised Vicsek model\n\nThe target is: \n\n\\begin{align}J^i_t = \\frac{\\frac{1}{N}\\sum_{j=1}^N K(|X^j_t-X^i_t|)V^j_t}{\\frac{1}{\\kappa}+\\frac{1}{\\kappa_0}|\\frac{1}{N}\\sum_{j=1}^N K(|X^j_t-X^i_t|)V^j_t|}.\\end{align}\n\nDefine the corresponding dictionary...\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "kappa_0 = 15.\n\nvariant = {\"name\" : \"max_kappa\", \"parameters\" : {\"kappa_max\" : kappa_0}}\n\nsimu = Vicsek(\n pos = pos.detach().clone(),\n vel = vel.detach().clone(), \n v = c, \n sigma = sigma, \n nu = nu, \n interaction_radius = R,\n box_size = L,\n dt = dt,\n variant = variant,\n block_sparse_reduction = True,\n number_of_cells = 40**2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally run the simulation over 300 units of time.... \n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "frames = [0., 5., 10., 30., 42., 71., 100, 124, 161, 206, 257, 300]\n\ns = time.time()\nit, op = display_kinetic_particles(simu, frames, order=True)\ne = time.time()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Print the total simulation time and the average time per iteration. \n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "print('Total time: '+str(e-s)+' seconds')\nprint('Average time per iteration: '+str((e-s)/simu.iteration)+' seconds')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot the histogram of the angles of the directions of motion. \n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "angle = torch.atan2(simu.vel[:,1],simu.vel[:,0])\nangle = angle.cpu().numpy()\nh = plt.hist(angle, bins=1000)\nplt.xlabel(\"angle\")\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The system self-organizes into a strongly clustered flock with band-like structures. \n\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Nematic Vicsek model \nThe target is: \n\n\\begin{align}J^i_t = \\kappa (V^i_t\\cdot \\overline{\\Omega}^i_t)\\overline{\\Omega}^i_t,\\end{align}\n\nwhere $\\overline{\\Omega}^i_t$ is any unit eigenvector associated to the maximal eigenvalue of the average Q-tensor:\n\n\\begin{align}Q^i_t = \\frac{1}{N}\\sum_{j=1}^N K(|X^j_t-X^i_t|){\\left(V^j_t\\otimes V^j_t - \\frac{1}{d} I_d\\right)}.\\end{align}\n\nDefine the corresponding dictionary...\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "variant = {\"name\" : \"nematic\", \"parameters\" : {}}\n\nsimu = Vicsek(\n pos = pos.detach().clone(),\n vel = vel.detach().clone(), \n v = c, \n sigma = sigma, \n nu = nu, \n interaction_radius = R,\n box_size = L,\n dt = dt,\n variant = variant,\n block_sparse_reduction = True, \n number_of_cells = 40**2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally run the simulation over 100 units of time... The color code indicates the angle of the direction of motion between $-\\pi$ and $\\pi$. \n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "frames = [0., 5., 10., 30., 42., 71., 100]\n\ns = time.time()\nit, op = display_kinetic_particles(simu, frames, order=True, color=True)\ne = time.time()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Print the total simulation time and the average time per iteration. \n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "print('Total time: '+str(e-s)+' seconds')\nprint('Average time per iteration: '+str((e-s)/simu.iteration)+' seconds')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot the histogram of the angles of the directions of motion. \n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# sphinx_gallery_thumbnail_number = -1\n\nangle = torch.atan2(simu.vel[:,1],simu.vel[:,0])\nangle = angle.cpu().numpy()\nh = plt.hist(angle, bins=1000)\nplt.xlabel(\"angle\")\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There are two modes separated by an angle $\\pi$ which indicates that two groups of equal size are moving in opposite direction. This a *nematic flock*. \n\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The target dictionary \n\nSeveral other targets are implemented. The complete list of available targets can be found in the **dictionary of targets** :attr:`target_method `.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "pprint.pprint(simu.target_method)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Custom targets can be added to the dictionary of targets using the method :meth:`add_target_method() `. Then a target can be readily used in a simulation using the method :meth:`compute_target() `. \n\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Options\n\nCustomized targets can also be defined by applying an **option** which modifies an existing target. See `examplemill`. \n\n" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.5" } }, "nbformat": 4, "nbformat_minor": 0 }