--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/latex/examples/glimpse-at-scipy.tex Tue Aug 31 18:45:36 2010 +0530
@@ -0,0 +1,426 @@
+\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}