path: root/day2/handout.tex

\documentclass[12pt]{article}
\usepackage{amsmath}
\title{Python Workshop\\Problems and Exercises}
\author{Asokan Pichai\\Prabhu Ramachandran}
\begin{document}
\maketitle

\section{Matrices and Arrays \& 2D Plotting}
\subsection{Matrices and Arrays}
\subsubsection{Basic Numpy}
\begin{verbatim}
# Simple array math example
>>> import numpy as np
>>> a = np.array([1,2,3,4])
>>> b = np.arange(2,6)
>>> b
array([2,3,4,5])
>>> a*2 + b + 1 # Basic math!
array([5, 8, 11, 14])

# Pi and e are defined.
>>> x = np.linspace(0.0, 10.0, 1000)
>>> x *= 2*np.pi/10 # inplace.
# apply functions to array.
>>> y = np.sin(x)
>>> z = np.exp(y)

>>> x = np.array([1., 2, 3, 4])
>>> np.size(x)
4
>>> x.dtype # What is a.dtype?
dtype('float64')
>>> x.shape
(4,)
>>> print np.rank(x), x.itemsize
1 8
>>> x[0] = 10
>>> print x[0], x[-1]
10.0 4.0

>>> a = np.array([[ 0, 1, 2, 3],
...            [10,11,12,13]])
>>> a.shape # (rows, columns)
(2, 4)
# Accessing and setting values
>>> a[1,3] 
13
>>> a[1,3] = -1
>>> a[1] # The second row
array([10,11,12,-1])

>>> b=np.array([[0,2,4,2],[1,2,3,4]])
>>> np.add(a,b,a)
>>> np.sum(x,axis=1)

>>> np.greater(a,4)
>>> np.sqrt(a)
\end{verbatim}

\subsubsection{Array Creation}
\begin{verbatim}
>>> np.array([2,3,4])  
array([2, 3, 4])

>>> np.linspace(0, 2, 4)   
array([0.,0.6666667,1.3333333,2.])

>>>np.ones([2,2])
array([[ 1.,  1.],
     [ 1.,  1.]])

>>>a = np.array([[1,2,3],[4,5,6]])
>>>np.ones_like(a)
array([[1, 1, 1],
       [1, 1, 1]])
\end{verbatim}
\subsubsection{Slicing, Striding Arrays}
\begin{verbatim}
>>> a = np.array([[1,2,3], [4,5,6], 
               [7,8,9]])
>>> a[0,1:3]
array([2, 3])
>>> a[1:,1:]
array([[5, 6],
       [8, 9]])
>>> a[:,2]
array([3, 6, 9])
>>> a[...,2]
array([3, 6, 9])

>>> a[0::2,0::2]
array([[1, 3],
       [7, 9]])
# Slices are references to the 
# same memory!
\end{verbatim}
\subsubsection{Random Numbers}
\begin{verbatim}
>>> np.random.rand(3,2)
array([[ 0.96276665,  0.77174861],
       [ 0.35138557,  0.61462271],
       [ 0.16789255,  0.43848811]])
>>> np.random.randint(1,100)
42
\end{verbatim}

\subsubsection{Problem Set}
\begin{verbatim}
    >>> from scipy import misc
    >>> A=misc.imread(name)
    >>> misc.imshow(A)
\end{verbatim}
\begin{enumerate}
  \item Convert an RGB image to Grayscale. $ Y = 0.5R + 0.25G + 0.25B $.
  \item Scale the image to 50\%.
  \item Introduce some random noise.
  \item Smooth the image using a mean filter.
  \\\small{Each element in the array is replaced by mean of all the neighbouring elements}
  \\\small{How fast does your code run?}
\end{enumerate}

\subsection{2D Plotting}
\subsubsection{Basic 2D Plotting}
\begin{verbatim}
$ ipython -pylab
>>> x = linspace(0, 2*pi, 1000)
>>> plot(x, sin(x)) 
>>> plot(x, sin(x), 'ro')
>>> xlabel(r'$\chi$', color='g')
# LaTeX markup!
>>> ylabel(r'sin($\chi$)', color='r')
>>> title('Simple figure', fontsize=20)
>>> savefig('/tmp/test.eps')
\end{verbatim}
\subsubsection{Tweaking plots}
\begin{verbatim}
# Set properties of objects:
>>> l, = plot(x, sin(x))
# Why "l,"?
>>> setp(l, linewidth=2.0, color='r')
>>> l.set_linewidth(2.0)
>>> draw() # Redraw.
>>> setp(l) # Print properties.
>>> clf() # Clear figure.
>>> close() # Close figure.
\end{verbatim}

\subsubsection{Working with text}
\begin{verbatim}
>>> w = arange(-2,2,.1)
>>> plot(w,exp(-(w*w))*cos)
>>> ylabel('$f(\omega)$')
>>> xlabel('$\omega$')
>>> title(r"$f(\omega)=e^{-\omega^2}
            cos({\omega^2})$")
>>> annotate('maxima',xy=(0, 1), 
             xytext=(1, 0.8), 
             arrowprops=dict(
             facecolor='black', 
             shrink=0.05))
\end{verbatim}

\subsubsection{Legends}
\begin{verbatim}
>>> x = linspace(0, 2*np.pi, 1000)
>>> plot(x, cos(5*x), 'r--', 
         label='cosine')
>>> plot(x, sin(5*x), 'g--', 
         label='sine')
>>> legend() 
# Or use:
>>> legend(['cosine', 'sine'])
\end{verbatim}

\subsubsection{Multiple figures}
\begin{verbatim}
>>> figure(1)
>>> plot(x, sin(x))
>>> figure(2)
>>> plot(x, tanh(x))
>>> figure(1)
>>> title('Easy as 1,2,3')

\end{verbatim}

\subsubsection{Problem Set}
  \begin{enumerate}
    \item Write a function that plots any regular n-gon given n.
    \item Consider the logistic map, $f(x) = kx(1-x)$, plot it for
          $k=2.5, 3.5$ and $4$ in the same plot.

  \item Consider the iteration $x_{n+1} = f(x_n)$ where $f(x) =
        kx(1-x)$.  Plot the successive iterates of this process.
  \item Plot this using a cobweb plot as follows:
  \begin{enumerate}
    \item Start at $(x_0, 0)$
    \item Draw line to $(x_i, f(x_i))$; 
    \item Set $x_{i+1} = f(x_i)$
    \item Draw line to $(x_i, x_i)$
    \item Repeat from 2 for as long as you want 
  \end{enumerate}
\end{enumerate}

\section{Advanced Numpy}

\subsection{Broadcasting}
\begin{verbatim}
>>> a = np.arange(4)
>>> b = np.arange(5)
>>> a+b #Does this work?
>>> a+3
>>> c=np.array([3])
>>> a+c #Works!
>>> b+c #But how?
>>> a.shape, b.shape, c.shape

>>> a = np.arange(4)
>>> a+3
array([3, 4, 5, 6])

>>> x = np.ones((3, 5))
>>> y = np.ones(8)
>>> (x[..., None] + y).shape
(3, 5, 8)

\end{verbatim}

\subsection{Copies \& Views}
\begin{verbatim}
>>> a = np.array([[1,2,3],[4,5,6]])
>>> b = a
>>> b is a
>>> b[0,0]=0; print a
>>> c = a.view()
>>> c is a
>>> c.base is a
>>> c.flags.owndata
>>> d = a.copy()
>>> d.base is a
>>> d.flags.owndata

>>> a = np.arange(1,9)
>>> a.shape=3,3
>>> b = a[0,1:3]
>>> c = a[0::2,0::2]
>>> a.flags.owndata
>>> b.flags.owndata
>>> b.base
>>> c.base is a

>>> b = a[np.array([0,1,2])]
array([[1, 2, 3],
       [4, 5, 6],
       [7, 8, 9]])
>>> b.flags.owndata
>>> abool=np.greater(a,2)
>>> c = a[abool]
>>> c.flags.owndata

\end{verbatim}

\section{Scipy}
\subsection{Linear Algebra}
\begin{verbatim}
>>> import scipy as sp
>>> from scipy import linalg
>>> A=sp.mat(np.arange(1,10))
>>> A.shape=3,3
>>> linalg.inv(A)
>>> linalg.det(A)
>>> linalg.norm(A)
>>> linalg.expm(A) #logm
>>> linalg.sinm(A) #cosm, tanm, ...

>>> A = sp.mat(np.arange(1,10))
>>> A.shape=3,3
>>> linalg.lu(A)
>>> linalg.eig(A)
>>> linalg.eigvals(A)
\end{verbatim}

\subsection{Solving Linear Equations}

\begin{align*}
  3x + 2y - z  & = 1 \\
  2x - 2y + 4z  & = -2 \\
  -x + \frac{1}{2}y -z & = 0
\end{align*}
To Solve this, 
\begin{verbatim}
>>> A = sp.mat([[3,2,-1],[2,-2,4]
              ,[-1,1/2,-1]])
>>> B = sp.mat([[1],[-2],[0]])
>>> linalg.solve(A,B)
\end{verbatim}

\subsection{Integrate}
\subsubsection{Quadrature}
Calculate the area under $(sin(x) + x^2)$ in the range $(0,1)$
\begin{verbatim}
>>> def f(x):
        return np.sin(x)+x**2
>>> integrate.quad(f, 0, 1)
\end{verbatim}

\subsubsection{ODE Integration}
Numerically solve ODEs\\
\begin{align*}
\frac{dx}{dt} &=-e^{-t}x^2\\ 
         x(0) &=2    
\end{align*}
\begin{verbatim}
>>> def dx_dt(x,t):
...     return -np.exp(-t)*x**2
>>> x=integrate.odeint(dx_dt, 2, t)
>>> plt.plot(x,t)
\end{verbatim}

\subsection{Interpolation}
\subsubsection{1D Interpolation}
\begin{verbatim}
>>> from scipy import interpolate
>>> interpolate.interp1d?
>>> x = np.arange(0,2*np.pi,np.pi/4)
>>> y = np.sin(x)
>>> fl = interpolate.interp1d(
            x,y,kind='linear')
>>> fc = interpolate.interp1d(
             x,y,kind='cubic')
>>> fl(np.pi/3)
>>> fc(np.pi/3)
\end{verbatim}

\subsubsection{Splines}
Plot the Cubic Spline of $sin(x)$
\begin{verbatim}
>>> x = np.arange(0,2*np.pi,np.pi/4)
>>> y = np.sin(x)
>>> tck = interpolate.splrep(x,y)
>>> X = np.arange(0,2*np.pi,np.pi/50)
>>> Y = interpolate.splev(X,tck,der=0)
>>> plt.plot(x,y,'o',x,y,X,Y)
>>> plt.show()
\end{verbatim}

\subsection{Signal \& Image Processing}
Applying a simple median filter
\begin{verbatim}
>>> from scipy import signal, ndimage
>>> from scipy import lena
>>> A=lena().astype('float32')
>>> B=signal.medfilt2d(A)
>>> imshow(B)
\end{verbatim}
Zooming an array - uses spline interpolation
\begin{verbatim}
>>> b=ndimage.zoom(A,0.5)
>>> imshow(b)
\end{verbatim}

\section{3D Data Visualization}
\subsection{Using mlab}
\begin{verbatim}
>>> from enthought.mayavi import mlab

>>> mlab.test_<TAB>
>>> mlab.test_contour3d()
>>> mlab.test_contour3d??
\end{verbatim}

\subsubsection{Plotting Functions}
\begin{verbatim}
>>> from numpy import *
>>> t = linspace(0, 2*pi, 50)
>>> u = cos(t)*pi
>>> x, y, z = sin(u), cos(u), sin(t)

>>> x = mgrid[-3:3:100j,-3:3:100j]
>>> z = sin(x*x + y*y)
>>> mlab.surf(x, y, z)
\end{verbatim}

\subsubsection{Large 2D Data}
\begin{verbatim}
>>> mlab.mesh(x, y, z)

>>> phi, theta = numpy.mgrid[0:pi:20j, 
...                         0:2*pi:20j]
>>> x = sin(phi)*cos(theta)
>>> y = sin(phi)*sin(theta)
>>> z = cos(phi)
>>> mlab.mesh(x, y, z, 
...           representation=
...           'wireframe')
\end{verbatim}

\subsubsection{Large 3D Data}
\begin{verbatim}
>>> x, y, z = ogrid[-5:5:64j, 
...                -5:5:64j, 
...                -5:5:64j]
>>> mlab.contour3d(x*x*0.5 + y*y + 
                   z*z*2)

>>> mlab.test_quiver3d()
\end{verbatim}

\subsection{Motivational Problem}
Atmospheric data of temperature over the surface of the earth. Let temperature ($T$) vary linearly with height ($z$)\\
$T = 288.15 - 6.5z$

\begin{verbatim}
lat = linspace(-89, 89, 37)
lon = linspace(0, 360, 37)
z = linspace(0, 100, 11)
\end{verbatim}

\begin{verbatim}
x, y, z = mgrid[0:360:37j,-89:89:37j,
                0:100:11j]
t = 288.15 - 6.5*z
mlab.contour3d(x, y, z, t)
mlab.outline()
mlab.colorbar()
\end{verbatim}

\subsection{Lorenz equation}
\begin{eqnarray*}
\frac{d x}{dt} &=& s (y-x)\\
\frac{d y}{d t} &=& rx -y -xz\\
\frac{d z}{d t} &=& xy - bz\\
\end{eqnarray*}

Let $s=10,$
$r=28,$ 
$b=8./3.$

\begin{verbatim}
x, y, z = mgrid[-50:50:20j,-50:50:20j,
                -10:60:20j]

\end{verbatim}

\end{document}