path: root/slides/latex/examples/glimpse-at-scipy.tex

\documentclass{article}

\usepackage{amsmath} 
\usepackage{graphicx} 
\usepackage{color}
\usepackage{listings} 
\usepackage{url}
%\definecolor{darkgreen}{rgb}{0,0.5,0}

\lstset{language=Python,
    basicstyle=\ttfamily\bfseries,
    commentstyle=\color{red}\itshape,
  stringstyle=\color{green},
  showstringspaces=false,
  keywordstyle=\color{blue}\bfseries}

\title{A Glimpse at Scipy}
\author{FOSSEE}
\date{June 2010}

\begin{document}
\maketitle

\begin{abstract}
This document shows a glimpse of the features of Scipy that will be
explored during this course.
\end{abstract}

\section{Introduction}

SciPy is open-source software for mathematics, science, and
engineering.

SciPy (pronounced ``Sigh Pie'') is a collection of mathematical
algorithms and convenience functions built on the Numpy extension for
Python. It adds significant power to the interactive Python session by
exposing the user to high-level commands and classes for the
manipulation and visualization of data. With SciPy, an interactive
Python session becomes a data-processing and system-prototyping
environment rivaling sytems such as \emph{Matlab, IDL, Octave, R-Lab,
  and Scilab}. \cite{scipy}

%% \begin{quote}
%%         	In 1998, ... I came across Python and its numerical extension
%%       	(Numeric) while I was looking for ways to analyze large data sets
%%       	... using a high-level language. I quickly fell in love with
%%       	Python programming which is a remarkable statement to make about a
%%       	programming language. If I had not seen others with the same view,
%%       	I might have seriously doubted my sanity.  -- Travis Oliphant, Creator of Numpy    
%% \end{quote}

\subsection{Sub-packages of Scipy}

SciPy is organized into subpackages covering different scientific
computing domains. These are summarized in the \underline{table
  \ref{subpkg}}.

\begin{table}
  \caption{Sub-packages available in Scipy}
  \label{subpkg}
\begin{tabular}{|l|l|}
\hline

\textbf{Subpackage} & \textbf{Description}\\

\hline

\texttt{cluster} & Clustering algorithms\\

\hline

\texttt{constants} & Physical and mathematical constants\\

\hline

\texttt{fftpack} & Fast Fourier Transform routines\\

\hline

\texttt{integrate} & Integration and ordinary differential equation
solvers\\

\hline

\texttt{interpolate} & Interpolation and smoothing splines\\

\hline

\texttt{io} & Input and Output\\

\hline

\texttt{linalg} & Linear algebra\\

\hline

\texttt{maxentropy} & Maximum entropy methods\\

\hline

\texttt{ndimage} & N-dimensional image processing\\

\hline

\texttt{odr} & Orthogonal distance regression\\

\hline

\texttt{optimize} & Optimization and root-finding routines\\

\hline

\texttt{signal} & Signal processing\\

\hline

\texttt{sparse} & Sparse matrices and associated routines\\

\hline

\texttt{spatial} & Spatial data structures and algorithms\\

\hline

\texttt{special} & Special functions\\

\hline

\texttt{stats} & Statistical distributions and functions\\

\hline

\texttt{weave} & C/C++ integration\\

\hline
\end{tabular}
\end{table}

\subsection{Use of Scipy in this course}
Following is a partial list of tasks we shall perform using Scipy, in
this course.

\begin{enumerate}
  \item Plotting \footnote{using \texttt{pylab} - see Appendix
    \ref{mpl}}
  \item Matrix Operations
  \begin{itemize}
    \item Inverse
    \item Determinant
  \end{itemize}
  \item Solving Equations
  \begin{itemize}
    \item System of Linear equations
    \item Polynomials
    \item Non-linear equations
  \end{itemize}
  \item Integration 
  \begin{itemize}
    \item Quadrature
    \item ODEs
  \end{itemize}
\end{enumerate}
\section{A Glimpse of Scipy functions}

This section gives a brief overview of the tasks that are going to be
performed using Scipy, in future classes of this course.

\subsection{Matrix Operations}

Let $\mathbf{A}$ be the matrix 
\(
\begin{bmatrix}
1 &3 &5\\
2 &5 &1\\
2 &3 &8
\end{bmatrix}
\) 

To input $\mathbf{A}$ matrix into python, we do the following in
\texttt{ipython}\footnote{\texttt{ipython} must be started with
  \texttt{-pylab} flag}\\

\begin{lstlisting}
  In []: A = array([[1,3,5],[2,5,1],[2,3,8]])
\end{lstlisting}

\subsubsection{Inverse}

The inverse of a matrix $\mathbf{A}$ is the matrix $\mathbf{B}$ such
that $\mathbf{A}\mathbf{B} = \mathbf{I}$ where $\mathbf{I}$ is the
identity matrix consisting of ones down the main diagonal. Usually
$\mathbf{B}$ is denoted $\mathbf{B} = \mathbf{A}^{-1}$ . In SciPy, the
matrix inverse of matrix $\mathbf{A}$ is obtained using

\lstinline+inv(A)+.
\begin{lstlisting}
  In []: inv(A)
  Out[]: 
  array([[-1.48,  0.36,  0.88],
         [ 0.56,  0.08, -0.36],
         [ 0.16, -0.12,  0.04]])
\end{lstlisting}

\subsubsection{Determinant}

The determinant of a square matrix $\mathbf{A}$ is denoted
$\left|\mathbf{A}\right|$. Suppose $a_{ij}$ are the elements of the
matrix $\mathbf{A}$ and let
$\mathbf{M}_{ij}=\left|\mathbf{A}_{ij}\right|$ be the determinant of
the matrix left by removing the $i^{th}$ row and $j^{th}$ column from
$\mathbf{A}$. Then for any row $i$

    \[ \left|\mathbf{A}\right|=\sum_{j}\left(-1\right)^{i+j}a_{ij}\mathbf{M}_{ij} \]

This is a recursive way to define the determinant where the base case
is defined by accepting that the determinant of a $1\times1$ matrix is
the only matrix element. In SciPy the determinant can be calculated
with $det$ . For example, the determinant of

    \[ \mathbf{A}=\begin{bmatrix} 1 & 3 & 5\\ 2 & 5 & 1\\ 2 & 3 & 8\end{bmatrix}\]

is

    \begin{eqnarray*} 
      |\mathbf{A}| & = & 1\begin{vmatrix} 5 & 1\\ 3 & 8\end{vmatrix}
                        -3\begin{vmatrix} 2 & 1\\ 2 & 8\end{vmatrix}
                        +5\begin{vmatrix}2 & 5\\ 2 & 3\end{vmatrix}\\  
                   & = & 1(5\cdot8-3\cdot1)-3(2\cdot8-2\cdot1)+5(2\cdot3-2\cdot5)=-25
    \end{eqnarray*}

In SciPy, this is computed as shown below

\begin{lstlisting}
  In []: A = array([[1, 3, 5], [2, 5, 1], [2, 3, 8]])
  In []: det(A)
  Out[]: -25.0
\end{lstlisting}

\subsection{Solving Equations}

\subsubsection{Linear Equations}

Solving linear systems of equations is straightforward using the scipy
command \lstinline+solve+. This command expects an input matrix and a
right-hand-side vector. The solution vector is then computed. An
option for entering a symmetrix matrix is offered which can speed up
the processing when applicable.  As an example, suppose it is desired
to solve the following simultaneous equations:

    \begin{eqnarray} x+3y+5z & = & 10\\ 2x+5y+z & = & 8\\ 2x+3y+8z & = & 3\end{eqnarray}

We could find the solution vector using a matrix inverse:

    \[ \left[\begin{array}{c} x\\ y\\ z\end{array}\right]=\left[\begin{array}{ccc} 1 & 3 & 5\\ 2 & 5 & 1\\ 2 & 3 & 8\end{array}\right]^{-1}\left[\begin{array}{c} 10\\ 8\\ 3\end{array}\right]=\frac{1}{25}\left[\begin{array}{c} -232\\ 129\\ 19\end{array}\right]=\left[\begin{array}{c} -9.28\\ 5.16\\ 0.76\end{array}\right] \]

However, it is better to use the solve command which can be faster and
more numerically stable. In this case it however gives the same
answer.

\begin{lstlisting}
  In []: A = array([[1, 3, 5], [2, 5, 1], [2, 3, 8]])
  In []: b = array([[10], [8], [3]])
  In []: dot(inv(A), b)
  Out[]: 
  array([[-9.28],
         [ 5.16],
         [ 0.76]])

  In []: solve(A,b)
  Out[]: 
  array([[-9.28],
         [ 5.16],
         [ 0.76]])
\end{lstlisting}

\subsubsection{Polynomials}

Solving a polynomial is straightforward in scipy using the
\lstinline+roots+ command. It expects the coefficients of the
polynomial in their decreasing order. For example, let's find the
roots of $x^3 - 2x^2 - \frac{1}{2}x + 1$ are $2$, $\sqrt{2}$ and
$-\sqrt{2}$. This is easy to see.

\begin{align*}
x^3 - 2x^2 - \frac{1}{2}x + 1 = 0\\
x^2(x-2) - \frac{1}{2}(x-2) = 0\\
(x-2)(x^2 - \frac{1}{2}) = 0\\
(x-2)(x - \frac{1}{\sqrt{2}})(x + \frac{1}{\sqrt{2}}) = 0
\end{align*}

We do it in scipy as shown below:
\begin{lstlisting}
  In []: coeff = array([1, -2, -2, 4])
  In []: roots(coeff)
\end{lstlisting}

\subsubsection{Non-linear Equations}

To find a root of a set of non-linear equations, the command
\lstinline+fsolve+ is needed. For example, the following example finds
the roots of the single-variable transcendental equation

    \[ x+2\cos\left(x\right)=0,\]

and the set of non-linear equations

    \begin{align} 
      x_{0}\cos\left(x_{1}\right) &= 4,\\ 
      x_{0}x_{1}-x_{1} &= 5
    \end{align}

The results are $x=-1.0299$ and $x_{0}=6.5041,\, x_{1}=0.9084$ .

\begin{lstlisting}
In []: def func(x):
   ...:     return x + 2*cos(x)

In []: def func2(x):
  ...:     out = [x[0]*cos(x[1]) - 4]
  ...:     out.append(x[1]*x[0] - x[1] - 5)
  ...:     return out

In []: from scipy.optimize import fsolve
In []: x0 = fsolve(func, 0.3)
In []: print x0
-1.02986652932

In []: x02 = fsolve(func2, [1, 1])
In []: print x02
[ 6.50409711  0.90841421]
\end{lstlisting}

\subsection{Integration}
% To be done in the lab. 
\subsubsection{Quadrature}

The function \texttt{quad} is provided to integrate a function of one
variable between two points. The points can be $\pm\infty$ ($\pm$
\texttt{inf}) to indicate infinite limits. For example, suppose you
wish to integrate the expression $e^{\sin(x)}$ in the interval
$[0,2\pi]$, i.e. $\int_0^{2\pi}e^{\sin(x)}dx$, it could be computed
using

\begin{lstlisting}
In []: def func(x):
   ...:     return exp(sin(x))

In []: from scipy.integrate import quad
In []: result = quad(func, 0, 2*pi)
In []: print result
(7.9549265210128457, 4.0521874164521979e-10)
\end{lstlisting}

\subsubsection{ODE}

We wish to solve an (a system of) Ordinary Differential Equation. For
this purpose, we shall use \lstinline{odeint}. As an illustration, let
us solve the ODE 

\begin{align} 
  \frac{dy}{dt} = ky(L-y)\\ 
  L = 25000,\,k = 0.00003,\,y(0) = 250 \nonumber 
\end{align}

We solve it in scipy as shown below. 

\begin{lstlisting}
In []: from scipy.integrate import odeint
In []: def f(y, t):
  ...:     k, L = 0.00003, 25000
  ...:     return k*y*(L-y)
  ...:
In []: t = linspace(0, 12, 60)
In []: y0 = 250
In []: y = odeint(f, y0, t)
\end{lstlisting}

Note: To solve a system of ODEs, we need to change the function to
return the right hand side of all the equations and the system and the
pass the required number of initial conditions to the
\lstinline{odeint} function.

\appendix

\section{Plotting using Pylab}\label{mpl}

The following piece of code, produces the plot in Figure \ref{fig:sin}
using \texttt{pylab}\cite{pylab} in \texttt{ipython}\footnote{start
  \texttt{ipython} with \texttt{-pylab} flag}\cite{ipy}

\begin{lstlisting}
In []: x = linspace(0, 2*pi, 50)
In []: plot(x, sin(x))
In []: title('Sine Curve between 0 and $\pi$')
In []: legend(['sin(x)'])
\end{lstlisting}

\begin{figure}[h!]
  \begin{center}
    \includegraphics[scale=0.4]{sine}    
  \end{center}
  \caption{Sine curve}
  \label{fig:sin}
\end{figure}



\begin{thebibliography}{9}
  \bibitem{scipy} 
    Eric Jones and Travis Oliphant and Pearu Peterson and others,
    \emph{SciPy: Open source scientific tools for Python}, 2001 -- , 
    \url{http://www.scipy.org/} 

  \bibitem{pylab}
   John D. Hunter, ``Matplotlib: A 2D Graphics Environment,''
   \emph{Computing in Science and Engineering}, vol. 9, no. 3,
   pp. 90-95, May/June 2007, doi:10.1109/MCSE.2007.55 

  \bibitem{ipy}
    Fernando Perez, Brian E. Granger, ``IPython: A System for
    Interactive Scientific Computing,'' \emph{Computing in Science and
    Engineering}, vol.~9, no.~3, pp.~21-29, May/June 2007,
    doi:10.1109/MCSE.2007.53. 

\end{thebibliography}
\end{document}