Blog home

$$ $$

The Zen of Python

In [1]:
import this

The Zen of Python, by Tim Peters

Beautiful is better than ugly.
Explicit is better than implicit.
Simple is better than complex.
Complex is better than complicated.
Flat is better than nested.
Sparse is better than dense.
Readability counts.
Special cases aren't special enough to break the rules.
Although practicality beats purity.
Errors should never pass silently.
Unless explicitly silenced.
In the face of ambiguity, refuse the temptation to guess.
There should be one-- and preferably only one --obvious way to do it.
Although that way may not be obvious at first unless you're Dutch.
Now is better than never.
Although never is often better than *right* now.
If the implementation is hard to explain, it's a bad idea.
If the implementation is easy to explain, it may be a good idea.
Namespaces are one honking great idea -- let's do more of those!
In [2]:
%matplotlib inline
from __future__ import division
import numpy as np
from numpy.random import rand
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.axes3d import Axes3D
from matplotlib import cm

Plotting using matplotlib

In [3]:
x = np.linspace(0, 4*np.pi, 64)
plt.plot(x, np.sin(x), '*-');

We can improve on the above plot in several ways. E.g

  • setting axis labels.
  • setting axis limits.
  • chossing color of our choice.
  • putting legends
  • ...

In the plot below we have included some of these features along with another plot in the same figure. This can be used to compare between two plots on the same figure.

In [4]:
x = np.linspace(0, 4*np.pi, 64)
plt.plot(x, np.sin(x), color="#348ABD", linewidth=2, linestyle="-", label='sin(x)');
plt.plot(x, np.cos(x), color="#A60628", linewidth= 3, linestyle="-", label='cos(x)');
plt.xlim([0, 4*np.pi]); plt.xlabel('x'); plt.legend(loc="lower left");

Plotting special functions using scipy

Here we plot the error function and the complementary error function. It can be seen from the plot that the two functions sum to unity.

In [5]:
from scipy.special import erf, erfc

f = plt.figure(figsize=(12, 5), dpi=80, facecolor='w', edgecolor='k')
x = np.linspace(0, 3, 64)
plt.plot(x, erf(x), color="#A60628", linewidth = 2, label='erf(x)')
plt.plot(x, erfc(x), color="#348ABD", linewidth = 2, label='erfc(x)')
plt.xlabel('x');plt.legend(loc="best");

The Lennard-Jones potential

$$ {\displaystyle V_{\text{LJ}}=4\varepsilon \left[\left({\frac {\sigma }{r}}\right)^{12}-\left({\frac {\sigma }{r}}\right)^{6}\right]=\varepsilon \left[\left({\frac {r_{\text{m}}}{r}}\right)^{12}-2\left({\frac {r_{\text{m}}}{r}}\right)^{6}\right],}$$

The above can be used as a purely repuslive potential as $$ \displaystyle V_{\text{WCA}}=\epsilon\left[\left(\frac{r_{\text m}}{r}\right)^{12}-2\left(\frac{r_{\text m}}{r}\right)^{6}\right]+\epsilon, $$ for $r$ less than $r_m$ and zero otherwise. We now plot these two potentials for different parameter values.

In [6]:
def plotLJ(sigma, epsilon):
    r = np.arange(0.1, 5, 0.01)
    rr = sigma/r
    V = 4*epsilon*(np.power(rr, 12) - np.power(rr, 6))
    plt.plot(r, V/epsilon, linewidth=2, label=' $\sigma$=%s, $\epsilon$%s'%(sigma, epsilon))
    plt.ylim(-1*epsilon, 5*epsilon); plt.xlabel('r'); plt.ylabel('V/$\epsilon$');   
    
    
plotLJ(sigma=1, epsilon=1)
plotLJ(sigma=2, epsilon=2)

plt.grid(); plt.legend(loc="upper right");
In [7]:
def plotWCA(rmin, epsilon):
    r = np.arange(0.1, rmin, 0.01)
    rr = rmin/r
    V = epsilon*(np.power(rr, 12) - 2*np.power(rr, 6)) + epsilon
    plt.plot(r, V/epsilon, linewidth=2, label=' $r_m$=%s, $\epsilon$%s'%(rmin, epsilon))
    plt.ylim(-1*epsilon, 5*epsilon); plt.xlabel('r'); plt.ylabel('V/$\epsilon$');   
    
    
plotWCA(rmin=2, epsilon=1)
plotWCA(rmin=3, epsilon=1)

plt.grid(); plt.legend(loc="lower left");

Plotting deformation of a two dimensional grid

Here we demonstrate that the bulk density only changes if the divergence of the strain tensor $u$ is non-zero. $$ $$ $$ \nabla \cdot u \neq 0$$ $$ $$ This is essentially because of the fact that the change in volume is related to the divergence of the strain.

Also, note that if the divergence of $u$ is zero then the body has only shear deformation but the density in the bulk remains constant. It is only in the case when the divergenc is non-zero we see the effect on this deformation in the bulk of the body.

In [8]:
def plotDeformation(N, fn):
    X, Y = np.meshgrid(np.arange(N), np.arange(N))
    if fn==1:
        u_x = np.sin(Y)
    elif fn==2:
        u_x = np.cos(X)
    u_y = 0
    X = X + u_x
    Y = Y + u_y
    T = np.arctan2(Y, X)
    ax.scatter(X, Y, s=75, c=T, alpha= 1, marker='o',  cmap=cm.coolwarm);
    plt.axis('off');

f = plt.figure(figsize=(18, 5), dpi=80, facecolor='w', edgecolor='k')
N = 20

ax = f.add_subplot(1, 2, 1)
plotDeformation(N,1)

ax = f.add_subplot(1, 2, 2)
plotDeformation(N,2)

Histograms

In [9]:
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 5))

N = 2**20
x = np.random.random(N)
ax1.hist(x, 100, normed=True, color='#A60628', alpha=0.4);

x = np.random.randn(N)
ax2.hist(x ,100, normed=True, color='#348ABD', alpha=0.4);

plt.suptitle("Histogram of random numbers drawn from uniform and Gaussian distribution",fontsize=20);

k-means and clustering

In [10]:
## first plot a random distribution which has clusters
N = 1024;
x1 = np.random.randn(N); y1 = np.random.randn(N)
x2=x1+5; y2=y1; x3=x1+5; y3=y1+5; x4=x1; y4=y1+5;
X = np.transpose(np.array([np.concatenate((x1, x2, x3, x4)), np.concatenate((y1, y2, y3, y4))]))
plt.scatter(X[:,0], X[:,1]);
In [11]:
# now using sklearn to plot these clusters
from sklearn.cluster import KMeans
km = KMeans(n_clusters=4)
km.fit(X)
labels = km.predict(X)
centrs =  km.cluster_centers_

plt.scatter(X[:,0], X[:,1], c=labels);
plt.scatter(centrs[:,0], centrs[:,1], marker='*', s=512);

3D plots

Here we give a demo of subplots in 3d using matplolib. The disturbances are are due to an initial Gaussian disturbance.

In [12]:
xx = np.linspace(-4, 4, 80)
yy = xx
x, y = np.meshgrid(xx, yy)

f = plt.figure(num=None, figsize=(15, 10), dpi=80, facecolor='w', edgecolor='k')

def plot_3dG(c_, n_):
    ax = f.add_subplot(2, 2, n_, projection='3d', )
    r = np.sqrt(x**2 + y**2)
    c = c_
    r1 = abs(r + c)
    r2 = abs(r - c)
    heights = 0.5*(np.exp(-r1*r1) + np.exp(-r2*r2))
    ax.plot_surface(x, y, heights, rstride=1, cstride=1, cmap=cm.coolwarm,
            linewidth=0, antialiased=False)

c, n = 0, 1           #plot 1
plot_3dG(c, n)

c, n = 2, 2           #plot 2
plot_3dG(c, n)

Contour plots

In the plot below we take the same equation as above but plot them using wireframe. Also plotted are the correseponding contours in all the direction.

In [13]:
xx = np.linspace(-2, 2, 80)
X, Y = np.meshgrid(xx, xx)
r = np.sqrt(X**2 + Y**2)
Z = 2*np.exp(-1.02*r*r)

f = plt.figure(num=None, figsize=(15, 10), dpi=80, facecolor='w', edgecolor='k')

ax = f.add_subplot(1,1,1, projection='3d')

ax.plot_wireframe(X, Y, Z, color="#348ABD", rstride=4, cstride=4, alpha=0.8)
cset = ax.contour(X, Y, Z, zdir='x', offset=-2,  alpha=0.6)
cset = ax.contour(X, Y, Z, zdir='y', offset=2,   alpha=0.6)
cset = ax.contour(X, Y, Z, zdir='z', offset=-2,  alpha=0.6)

ll = 2; ax.set_xlim3d(-ll, ll);    ax.set_ylim3d(-ll, ll);     ax.set_zlim3d(-ll, ll);

xkcd plots

In [14]:
with plt.xkcd():
    x = np.linspace(0, 10*np.pi, 256)
    plt.plot(x, np.sin(x)*np.exp(-0.1*x), label='sin(x)exp(-0.1x)')
    plt.plot(x, np.cos(x)*np.exp(-0.1*x), label='cos(x)exp(-0.1x)')
    plt.xlabel('x');   plt.xlim([0, 10*np.pi]); plt.ylabel(' f(x)');
    plt.legend(loc='lower right');  plt.title('xkcd plots :D');
plt.rcdefaults()