Note

Click here to download the full example code

Tutorial 02: Targets and options

Many useful local averages (see Tutorial 05: Kernels and 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.

An example: variants of the Vicsek model

In its most abstract form, the Vicsek model reads:

\[\mathrm{d}X^i_t = c_0 V^i_t \mathrm{d} t\]
\[\mathrm{d}V^i_t = \sigma\mathsf{P}(V^i_t)\circ ( J^i_t \mathrm{d}t + \mathrm{d} B^i_t),\]

where \(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.

First, some standard imports…

import time
import torch
import pprint
from matplotlib import pyplot as plt
from sisyphe.models import Vicsek
from sisyphe.display import display_kinetic_particles

use_cuda = torch.cuda.is_available()
dtype = torch.cuda.FloatTensor if use_cuda else torch.FloatTensor

The classical non-normalised Vicsek model

The target is:

\[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|}.\]

The parameters of the model…

N = 100000
L = 100

pos = L*torch.rand((N,2)).type(dtype)
vel = torch.randn(N,2).type(dtype)
vel = vel/torch.norm(vel,dim=1).reshape((N,1))

R = 3.
c = 3.
nu = 5.
sigma = 1.

dt = .01

The choice of the target is implemented in the keyword argument variant. For the classical normalised target, it is given by the following dictionary.

variant = {"name" : "normalised", "parameters" : {}}

simu = Vicsek(
    pos = pos.detach().clone(),
    vel = vel.detach().clone(),
    v = c,
    sigma = sigma,
    nu = nu,
    interaction_radius = R,
    box_size = L,
    dt = dt,
    variant = variant,
    block_sparse_reduction = True,
    number_of_cells = 40**2)

Finally run the simulation over 300 units of time….

frames = [0., 5., 10., 30., 42., 71., 100, 124, 161, 206, 257, 300]

s = time.time()
it, op = display_kinetic_particles(simu, frames, order=True)
e = time.time()
  • Vicsek(normalised) Parameters: N=100000 ; R=3.0 ; nu=5.0 ; sigma=1.0 ; v=3.0 Time=0.0
  • Vicsek(normalised) Parameters: N=100000 ; R=3.0 ; nu=5.0 ; sigma=1.0 ; v=3.0 Time=5.0
  • Vicsek(normalised) Parameters: N=100000 ; R=3.0 ; nu=5.0 ; sigma=1.0 ; v=3.0 Time=10.0
  • Vicsek(normalised) Parameters: N=100000 ; R=3.0 ; nu=5.0 ; sigma=1.0 ; v=3.0 Time=30.0
  • Vicsek(normalised) Parameters: N=100000 ; R=3.0 ; nu=5.0 ; sigma=1.0 ; v=3.0 Time=42.0
  • Vicsek(normalised) Parameters: N=100000 ; R=3.0 ; nu=5.0 ; sigma=1.0 ; v=3.0 Time=71.0
  • Vicsek(normalised) Parameters: N=100000 ; R=3.0 ; nu=5.0 ; sigma=1.0 ; v=3.0 Time=100.0
  • Vicsek(normalised) Parameters: N=100000 ; R=3.0 ; nu=5.0 ; sigma=1.0 ; v=3.0 Time=124.0
  • Vicsek(normalised) Parameters: N=100000 ; R=3.0 ; nu=5.0 ; sigma=1.0 ; v=3.0 Time=161.0
  • Vicsek(normalised) Parameters: N=100000 ; R=3.0 ; nu=5.0 ; sigma=1.0 ; v=3.0 Time=206.0
  • Vicsek(normalised) Parameters: N=100000 ; R=3.0 ; nu=5.0 ; sigma=1.0 ; v=3.0 Time=257.0
  • Vicsek(normalised) Parameters: N=100000 ; R=3.0 ; nu=5.0 ; sigma=1.0 ; v=3.0 Time=300.0
  • plot b targetsoptions

Out:

Progress:0%
Progress:1%
Progress:2%
Progress:3%
Progress:4%
Progress:5%
Progress:6%
Progress:7%
Progress:8%
Progress:9%
Progress:10%
Progress:11%
Progress:12%
Progress:13%
Progress:14%
Progress:15%
Progress:16%
Progress:17%
Progress:18%
Progress:19%
Progress:20%
Progress:21%
Progress:22%
Progress:23%
Progress:24%
Progress:25%
Progress:26%
Progress:27%
Progress:28%
Progress:29%
Progress:30%
Progress:31%
Progress:32%
Progress:33%
Progress:34%
Progress:35%
Progress:36%
Progress:37%
Progress:38%
Progress:39%
Progress:40%
Progress:41%
Progress:42%
Progress:43%
Progress:44%
Progress:45%
Progress:46%
Progress:47%
Progress:48%
Progress:49%
Progress:50%
Progress:51%
Progress:52%
Progress:53%
Progress:54%
Progress:55%
Progress:56%
Progress:57%
Progress:58%
Progress:59%
Progress:60%
Progress:61%
Progress:62%
Progress:63%
Progress:64%
Progress:65%
Progress:66%
Progress:67%
Progress:68%
Progress:69%
Progress:70%
Progress:71%
Progress:72%
Progress:73%
Progress:74%
Progress:75%
Progress:76%
Progress:77%
Progress:78%
Progress:79%
Progress:80%
Progress:81%
Progress:82%
Progress:83%
Progress:84%
Progress:85%
Progress:86%
Progress:87%
Progress:88%
Progress:89%
Progress:90%
Progress:91%
Progress:92%
Progress:93%
Progress:94%
Progress:95%
Progress:96%
Progress:97%
Progress:98%
Progress:99%

Print the total simulation time and the average time per iteration.

print('Total time: '+str(e-s)+' seconds')
print('Average time per iteration: '+str((e-s)/simu.iteration)+' seconds')

Out:

Total time: 145.67891550064087 seconds
Average time per iteration: 0.004855963850021362 seconds

Plot the histogram of the angles of the directions of motion.

angle = torch.atan2(simu.vel[:,1],simu.vel[:,0])
angle = angle.cpu().numpy()
h = plt.hist(angle, bins=1000)
plt.xlabel("angle")
plt.show()
plot b targetsoptions

After an initial clustering phase, the system self-organizes into a uniform flock.

Non-normalised Vicsek model

The target is:

\[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|}.\]

Define the corresponding dictionary…

kappa_0 = 15.

variant = {"name" : "max_kappa", "parameters" : {"kappa_max" : kappa_0}}

simu = Vicsek(
    pos = pos.detach().clone(),
    vel = vel.detach().clone(),
    v = c,
    sigma = sigma,
    nu = nu,
    interaction_radius = R,
    box_size = L,
    dt = dt,
    variant = variant,
    block_sparse_reduction = True,
    number_of_cells = 40**2)

Finally run the simulation over 300 units of time….

frames = [0., 5., 10., 30., 42., 71., 100, 124, 161, 206, 257, 300]

s = time.time()
it, op = display_kinetic_particles(simu, frames, order=True)
e = time.time()
  • Vicsek(max_kappa) Parameters: N=100000 ; R=3.0 ; nu=5.0 ; sigma=1.0 ; v=3.0 Time=0.0
  • Vicsek(max_kappa) Parameters: N=100000 ; R=3.0 ; nu=5.0 ; sigma=1.0 ; v=3.0 Time=5.0
  • Vicsek(max_kappa) Parameters: N=100000 ; R=3.0 ; nu=5.0 ; sigma=1.0 ; v=3.0 Time=10.0
  • Vicsek(max_kappa) Parameters: N=100000 ; R=3.0 ; nu=5.0 ; sigma=1.0 ; v=3.0 Time=30.0
  • Vicsek(max_kappa) Parameters: N=100000 ; R=3.0 ; nu=5.0 ; sigma=1.0 ; v=3.0 Time=42.0
  • Vicsek(max_kappa) Parameters: N=100000 ; R=3.0 ; nu=5.0 ; sigma=1.0 ; v=3.0 Time=71.0
  • Vicsek(max_kappa) Parameters: N=100000 ; R=3.0 ; nu=5.0 ; sigma=1.0 ; v=3.0 Time=100.0
  • Vicsek(max_kappa) Parameters: N=100000 ; R=3.0 ; nu=5.0 ; sigma=1.0 ; v=3.0 Time=124.0
  • Vicsek(max_kappa) Parameters: N=100000 ; R=3.0 ; nu=5.0 ; sigma=1.0 ; v=3.0 Time=161.0
  • Vicsek(max_kappa) Parameters: N=100000 ; R=3.0 ; nu=5.0 ; sigma=1.0 ; v=3.0 Time=206.0
  • Vicsek(max_kappa) Parameters: N=100000 ; R=3.0 ; nu=5.0 ; sigma=1.0 ; v=3.0 Time=257.0
  • Vicsek(max_kappa) Parameters: N=100000 ; R=3.0 ; nu=5.0 ; sigma=1.0 ; v=3.0 Time=300.0
  • plot b targetsoptions

Out:

Progress:0%
Progress:1%
Progress:2%
Progress:3%
Progress:4%
Progress:5%
Progress:6%
Progress:7%
Progress:8%
Progress:9%
Progress:10%
Progress:11%
Progress:12%
Progress:13%
Progress:14%
Progress:15%
Progress:16%
Progress:17%
Progress:18%
Progress:19%
Progress:20%
Progress:21%
Progress:22%
Progress:23%
Progress:24%
Progress:25%
Progress:26%
Progress:27%
Progress:28%
Progress:29%
Progress:30%
Progress:31%
Progress:32%
Progress:33%
Progress:34%
Progress:35%
Progress:36%
Progress:37%
Progress:38%
Progress:39%
Progress:40%
Progress:41%
Progress:42%
Progress:43%
Progress:44%
Progress:45%
Progress:46%
Progress:47%
Progress:48%
Progress:49%
Progress:50%
Progress:51%
Progress:52%
Progress:53%
Progress:54%
Progress:55%
Progress:56%
Progress:57%
Progress:58%
Progress:59%
Progress:60%
Progress:61%
Progress:62%
Progress:63%
Progress:64%
Progress:65%
Progress:66%
Progress:67%
Progress:68%
Progress:69%
Progress:70%
Progress:71%
Progress:72%
Progress:73%
Progress:74%
Progress:75%
Progress:76%
Progress:77%
Progress:78%
Progress:79%
Progress:80%
Progress:81%
Progress:82%
Progress:83%
Progress:84%
Progress:85%
Progress:86%
Progress:87%
Progress:88%
Progress:89%
Progress:90%
Progress:91%
Progress:92%
Progress:93%
Progress:94%
Progress:95%
Progress:96%
Progress:97%
Progress:98%
Progress:99%

Print the total simulation time and the average time per iteration.

print('Total time: '+str(e-s)+' seconds')
print('Average time per iteration: '+str((e-s)/simu.iteration)+' seconds')

Out:

Total time: 156.8030505180359 seconds
Average time per iteration: 0.005226768350601196 seconds

Plot the histogram of the angles of the directions of motion.

angle = torch.atan2(simu.vel[:,1],simu.vel[:,0])
angle = angle.cpu().numpy()
h = plt.hist(angle, bins=1000)
plt.xlabel("angle")
plt.show()
plot b targetsoptions

The system self-organizes into a strongly clustered flock with band-like structures.

Nematic Vicsek model

The target is:

\[J^i_t = \kappa (V^i_t\cdot \overline{\Omega}^i_t)\overline{\Omega}^i_t,\]

where \(\overline{\Omega}^i_t\) is any unit eigenvector associated to the maximal eigenvalue of the average Q-tensor:

\[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)}.\]

Define the corresponding dictionary…

variant = {"name" : "nematic", "parameters" : {}}

simu = Vicsek(
    pos = pos.detach().clone(),
    vel = vel.detach().clone(),
    v = c,
    sigma = sigma,
    nu = nu,
    interaction_radius = R,
    box_size = L,
    dt = dt,
    variant = variant,
    block_sparse_reduction = True,
    number_of_cells = 40**2)

Finally run the simulation over 100 units of time… The color code indicates the angle of the direction of motion between \(-\pi\) and \(\pi\).

frames = [0., 5., 10., 30., 42., 71., 100]

s = time.time()
it, op = display_kinetic_particles(simu, frames, order=True, color=True)
e = time.time()
  • Vicsek(nematic) Parameters: N=100000 ; R=3.0 ; nu=5.0 ; sigma=1.0 ; v=3.0 Time=0.0
  • Vicsek(nematic) Parameters: N=100000 ; R=3.0 ; nu=5.0 ; sigma=1.0 ; v=3.0 Time=5.0
  • Vicsek(nematic) Parameters: N=100000 ; R=3.0 ; nu=5.0 ; sigma=1.0 ; v=3.0 Time=10.0
  • Vicsek(nematic) Parameters: N=100000 ; R=3.0 ; nu=5.0 ; sigma=1.0 ; v=3.0 Time=30.0
  • Vicsek(nematic) Parameters: N=100000 ; R=3.0 ; nu=5.0 ; sigma=1.0 ; v=3.0 Time=42.0
  • Vicsek(nematic) Parameters: N=100000 ; R=3.0 ; nu=5.0 ; sigma=1.0 ; v=3.0 Time=71.0
  • Vicsek(nematic) Parameters: N=100000 ; R=3.0 ; nu=5.0 ; sigma=1.0 ; v=3.0 Time=100.0
  • plot b targetsoptions

Out:

Progress:0%
Progress:1%
Progress:2%
Progress:3%
Progress:4%
Progress:5%
Progress:6%
Progress:7%
Progress:8%
Progress:9%
Progress:10%
Progress:11%
Progress:12%
Progress:13%
Progress:14%
Progress:15%
Progress:16%
Progress:17%
Progress:18%
Progress:19%
Progress:20%
Progress:21%
Progress:22%
Progress:23%
Progress:24%
Progress:25%
Progress:26%
Progress:27%
Progress:28%
Progress:29%
Progress:30%
Progress:31%
Progress:32%
Progress:33%
Progress:34%
Progress:35%
Progress:36%
Progress:37%
Progress:38%
Progress:39%
Progress:40%
Progress:41%
Progress:42%
Progress:43%
Progress:44%
Progress:45%
Progress:46%
Progress:47%
Progress:48%
Progress:49%
Progress:50%
Progress:51%
Progress:52%
Progress:53%
Progress:54%
Progress:55%
Progress:56%
Progress:57%
Progress:58%
Progress:59%
Progress:60%
Progress:61%
Progress:62%
Progress:63%
Progress:64%
Progress:65%
Progress:66%
Progress:67%
Progress:68%
Progress:69%
Progress:70%
Progress:71%
Progress:72%
Progress:73%
Progress:74%
Progress:75%
Progress:76%
Progress:77%
Progress:78%
Progress:79%
Progress:80%
Progress:81%
Progress:82%
Progress:83%
Progress:84%
Progress:85%
Progress:86%
Progress:87%
Progress:88%
Progress:89%
Progress:90%
Progress:91%
Progress:92%
Progress:93%
Progress:94%
Progress:95%
Progress:96%
Progress:97%
Progress:98%
Progress:99%

Print the total simulation time and the average time per iteration.

print('Total time: '+str(e-s)+' seconds')
print('Average time per iteration: '+str((e-s)/simu.iteration)+' seconds')

Out:

Total time: 56.330119132995605 seconds
Average time per iteration: 0.00563301191329956 seconds

Plot the histogram of the angles of the directions of motion.

# sphinx_gallery_thumbnail_number = -1

angle = torch.atan2(simu.vel[:,1],simu.vel[:,0])
angle = angle.cpu().numpy()
h = plt.hist(angle, bins=1000)
plt.xlabel("angle")
plt.show()
plot b targetsoptions

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.

The target dictionary

Several other targets are implemented. The complete list of available targets can be found in the dictionary of targets target_method.

pprint.pprint(simu.target_method)

Out:

{'max_kappa': <bound method KineticParticles.max_kappa of <sisyphe.models.Vicsek object at 0x7efb4e14e970>>,
 'mean_field': <bound method KineticParticles.mean_field of <sisyphe.models.Vicsek object at 0x7efb4e14e970>>,
 'morse': <bound method Particles.morse_target of <sisyphe.models.Vicsek object at 0x7efb4e14e970>>,
 'motsch_tadmor': <bound method KineticParticles.motsch_tadmor of <sisyphe.models.Vicsek object at 0x7efb4e14e970>>,
 'nematic': <bound method KineticParticles.nematic of <sisyphe.models.Vicsek object at 0x7efb4e14e970>>,
 'normalised': <bound method KineticParticles.normalised of <sisyphe.models.Vicsek object at 0x7efb4e14e970>>,
 'overlapping_repulsion': <bound method Particles.overlapping_repulsion_target of <sisyphe.models.Vicsek object at 0x7efb4e14e970>>,
 'quadratic_potential': <bound method Particles.quadratic_potential_target of <sisyphe.models.Vicsek object at 0x7efb4e14e970>>}

Custom targets can be added to the dictionary of targets using the method add_target_method(). Then a target can be readily used in a simulation using the method compute_target().

Options

Customized targets can also be defined by applying an option which modifies an existing target. See Mills.

Total running time of the script: ( 6 minutes 25.941 seconds)

Gallery generated by Sphinx-Gallery