diff -r 9b5c25c0fec3 -r 313bebeb7862 latex/examples/glimpse-at-scipy.tex --- /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}