Category: Algorithms

Minimal Vicsek model in Python

The Vicsek model is one of the simplest models for active matter. It displays interesting features, such as swarming.

Large scale simulations are often needed in order to provide firm statements on the statistical properties of this kind of models. However, for a pedagogical and illustrative purpose it may be useful to have an elementary code to play with. For this purpose, I have written a relatively simple Python code which implements the model with a few clever tricks to make simulations of a few thousands of agents possible on a standard laptop.

We follow Gregoire and Chaté in the formalism: point-wise particles move synchronously at constant speed v in discretised time of steps Δt=1. The particles have an orientation described by an angle θ which evolves taking into account all particles k within a given radius of interaction

\theta_{j}^{t+1}=\arg \left[\sum_{k \sim j} e^{i \theta_{k}^{t}}\right]+\eta \xi_{j}^{t}

For the neighbourhood calculations, cell-lists would be ideal, but they are too complex for the kind of elementary code that we want to write. What we are going to use is the kd-tree quick nearest neighbour lookup algorithm as implemented in Scipy and some clever sparse matrix manipulation. For visualisation, we employ matplotlib, so that the resulting code is just 60 lines with only very popular libraries.

import numpy as np
import scipy as sp
from scipy import sparse
from scipy.spatial import cKDTree
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation


L = 32.0
rho = 3.0
N =	int(rho*L**2)
print(" N",N)

r0 = 1.0
deltat = 1.0
factor =0.5
v0 = r0/deltat*factor
iterations = 10000
eta = 0.15


pos = np.random.uniform(0,L,size=(N,2))
orient = np.random.uniform(-np.pi, np.pi,size=N)

fig, ax= plt.subplots(figsize=(6,6))

qv = ax.quiver(pos[:,0], pos[:,1], np.cos(orient[0]), np.sin(orient), orient, clim=[-np.pi, np.pi])

def animate(i):
	print(i)

	global orient
	tree = cKDTree(pos,boxsize=[L,L])
	dist = tree.sparse_distance_matrix(tree, max_distance=r0,output_type='coo_matrix')

	#important 3 lines: we evaluate a quantity for every column j
	data = np.exp(orient[dist.col]*1j)
	# construct  a new sparse marix with entries in the same places ij of the dist matrix
	neigh = sparse.coo_matrix((data,(dist.row,dist.col)), shape=dist.get_shape())
	# and sum along the columns (sum over j)
	S = np.squeeze(np.asarray(neigh.tocsr().sum(axis=1)))
	
	
	orient = np.angle(S)+eta*np.random.uniform(-np.pi, np.pi, size=N)


	cos, sin= np.cos(orient), np.sin(orient)
	pos[:,0] += cos*v0
	pos[:,1] += sin*v0

	pos[pos>L] -= L
	pos[pos<0] += L

	qv.set_offsets(pos)
	qv.set_UVC(cos, sin,orient)
	return qv,

anim = FuncAnimation(fig,animate,np.arange(1, 200),interval=1, blit=True)
plt.show()

The result is the following animation (the colour indicates the orientation):

Box Counting in Numpy

The fractal dimension of an object is a single scalar number that allows us to quantify how compact an object is , i.e. how wiggly a curve is, how wrinkled a surface is, how porous a complex volume is.

However, estimating the fractal dimension of an object is not an easy task, and many methods exist. The simplest method is box counting: the idea is to fully cover the object with many boxes of a given size, count how many boxes are needed to cover the object and repeat the process for many box sizes. The scaling of the number of boxes covering the object with the size of the boxes gives an estimate for the fractal dimension of the object.

Sierpinski

The algorithm has many limitations, but in its simplest form it can be easily implemented in Python. The idea is to simply bin the object in a histogram of variable bin sizes. This can be easily generalised to any dimensions, thanks to Numpy’s histrogramdd  function.

In the following code, I test the idea with a known fractal, Sierpinski’s triangle, which has an exact (Hausdorff) fractal dimension of log(3)/log(2).

import numpy as np
import pylab as pl

def rgb2gray(rgb):
	r, g, b = rgb[:,:,0], rgb[:,:,1], rgb[:,:,2]
	gray = 0.2989 * r + 0.5870 * g + 0.1140 * b
	return gray

image=rgb2gray(pl.imread("Sierpinski.png"))

# finding all the non-zero pixels
pixels=[]
for i in range(image.shape[0]):
	for j in range(image.shape[1]):
		if image[i,j]>0:
			pixels.append((i,j))

Lx=image.shape[1]
Ly=image.shape[0]
print (Lx, Ly)
pixels=pl.array(pixels)
print (pixels.shape)

# computing the fractal dimension
#considering only scales in a logarithmic list
scales=np.logspace(0.01, 1, num=10, endpoint=False, base=2)
Ns=[]
# looping over several scales
for scale in scales:
	print ("======= Scale :",scale)
	# computing the histogram
	H, edges=np.histogramdd(pixels, bins=(np.arange(0,Lx,scale),np.arange(0,Ly,scale)))
	Ns.append(np.sum(H>0))

# linear fit, polynomial of degree 1
coeffs=np.polyfit(np.log(scales), np.log(Ns), 1)

pl.plot(np.log(scales),np.log(Ns), 'o', mfc='none')
pl.plot(np.log(scales), np.polyval(coeffs,np.log(scales)))
pl.xlabel('log $\epsilon$')
pl.ylabel('log N')
pl.savefig('sierpinski_dimension.pdf')

print ("The Hausdorff dimension is", -coeffs[0]) #the fractal dimension is the OPPOSITE of the fitting coefficient
np.savetxt("scaling.txt", list(zip(scales,Ns)))