scipy/IntroductionSciPy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![NASA](http://www.nasa.gov/sites/all/themes/custom/nasatwo/images/nasa-logo.svg)\n",
    "\n",
    "<center>\n",
    "<h1><font size=\"+3\">GSFC Python Bootcamp</font></h1>\n",
    "</center>\n",
    "\n",
    "---\n",
    "\n",
    "<CENTER>\n",
    "<H1 style=\"color:red\">\n",
    "Introduction to SciPy\n",
    "</H1>\n",
    "</CENTER>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import print_function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font color='red'>What will be Covered?</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<OL>\n",
    "<LI> What is SciPy?\n",
    "<LI> SciPy Organization\n",
    "<LI> Sample SciPy Modules\n",
    "<UL>\n",
    "    <LI> interpolate\n",
    "    <LI> integrate\n",
    "    <LI> optimize\n",
    "    <LI> fftpack\n",
    "    <LI> Stats\n",
    "</UL>\n",
    "</OL>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font color='red'>Reference Documents </font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<OL>\n",
    "<LI> <A HREF=\"http://www.sam.math.ethz.ch/~raoulb/teaching/PythonTutorial/intro_scipy.html\">First Introduction to SciPy</A>\n",
    "<LI> <A HREF=\"http://www.scipy.org/\">Home of SciPy</A>\n",
    "<LI> <A HREF=\"http://docs.scipy.org/doc/scipy/reference/#reference\">SciPy Reference</A>\n",
    "</OL>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In case you prefer a video:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import YouTubeVideo\n",
    "YouTubeVideo(\"oYTs9HwFGbY\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font color='red'>What is SciPy?</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<UL>\n",
    "<LI> Collection of mathematical algorithms and convenience functions built on the Numpy extension for Python.\n",
    "<LI> Provides high-level commands and classes for manipulating and visualizing data.\n",
    "<LI> Contains subpackages covering different scientific domains such as Linear Algebra, Signal Processing, Statistics, etc. \n",
    "</UL>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### In this presentation, we only want to show the capabilities of SciPy."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font color='red'>SciPy Organization</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following packages are available in SciPy:\n",
    "\n",
    "<table>\n",
    "<tr> <td> <b> Subpackage </b></td> <td><b> Description </b></td></tr>\n",
    "<tr> <td>cluster</td> <td> \tClustering algorithms</td> </tr>\n",
    "<tr> <td>constants</td> <td> \tPhysical and mathematical constants</td> </tr>\n",
    "<tr> <td>fftpack</td> <td> \tFast Fourier Transform routines</td> </tr>\n",
    "<tr> <td>integrate</td> <td> \tIntegration and ordinary differential equation solvers</td> </tr>\n",
    "<tr> <td>interpolate</td> <td> \tInterpolation and smoothing splines</td> </tr>\n",
    "<tr> <td>io</td> <td> \tInput and Output</td> </tr>\n",
    "<tr> <td>linalg</td> <td> \tLinear algebra</td> </tr>\n",
    "<tr> <td>ndimage</td> <td> \tN-dimensional image processing</td> </tr>\n",
    "<tr> <td>odr</td> <td> \tOrthogonal distance regression</td> </tr>\n",
    "<tr> <td>optimize</td> <td> \tOptimization and root-finding routines</td> </tr>\n",
    "<tr> <td>signal</td> <td> \tSignal processing</td> </tr>\n",
    "<tr> <td>sparse</td> <td> \tSparse matrices and associated routines</td> </tr>\n",
    "<tr> <td>spatial</td> <td> \tSpatial data structures and algorithms</td> </tr>\n",
    "<tr> <td>special</td> <td> \tSpecial functions</td> </tr>\n",
    "<tr> <td>stats</td> <td> \tStatistical distributions and functions</td> </tr>\n",
    "<tr> <td>weave</td> <td> \tC/C++ integration</td> </tr>\n",
    "</table>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Scipy sub-packages need to be imported separately, for example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import scipy as sp\n",
    "from scipy import linalg, optimize, interpolate, special, integrate, fftpack"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font color='red'>Interpolation Problems</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### B-Spline Interpolation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# make ascending spiral in 3-space\n",
    "t=np.linspace(0,1.75*2*np.pi,100)\n",
    "\n",
    "x = np.sin(t)\n",
    "y = np.cos(t)\n",
    "z = t\n",
    "\n",
    "# add noise\n",
    "x+= np.random.normal(scale=0.1, size=x.shape)\n",
    "y+= np.random.normal(scale=0.1, size=y.shape)\n",
    "z+= np.random.normal(scale=0.1, size=z.shape)\n",
    "\n",
    "# spline parameters\n",
    "s = 3.0   # smoothness parameter\n",
    "k = 2     # spline order\n",
    "nest = -1 # estimate of number of knots needed (-1 = maximal)\n",
    "\n",
    "# Find the B-spline representation \n",
    "tckp,u = interpolate.splprep([x,y,z],s=s,k=k,nest=-1)\n",
    "\n",
    "# evaluate spline, including interpolated points\n",
    "xnew,ynew,znew = interpolate.splev(np.linspace(0,1,400),tckp)\n",
    "\n",
    "plt.subplot(2,2,1)\n",
    "data,=plt.plot(x,y,'bo-',label='data')\n",
    "fit,=plt.plot(xnew,ynew,'r-',label='fit')\n",
    "plt.legend()\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "\n",
    "plt.subplot(2,2,2)\n",
    "data,=plt.plot(x,z,'bo-',label='data')\n",
    "fit,=plt.plot(xnew,znew,'r-',label='fit')\n",
    "plt.legend()\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('z')\n",
    "\n",
    "plt.subplot(2,2,3)\n",
    "data,=plt.plot(y,z,'bo-',label='data')\n",
    "fit,=plt.plot(ynew,znew,'r-',label='fit')\n",
    "plt.legend()\n",
    "plt.xlabel('y')\n",
    "plt.ylabel('z')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Bivariate B-Spine Interpolation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x,y = np.mgrid[-1:1:20j,-1:1:20j]\n",
    "z = (x+y)*np.exp(-6.0*(x*x+y*y))\n",
    "plt.subplot(2,1,1)\n",
    "plt.pcolor(x,y,z)\n",
    "plt.colorbar()\n",
    "plt.title(\"Sparsely sampled function.\")\n",
    "\n",
    "# Interpolate function over new 70x70 grid\n",
    "\n",
    "xnew,ynew = np.mgrid[-1:1:70j,-1:1:70j]\n",
    "tck = interpolate.bisplrep(x,y,z,s=0)\n",
    "znew = interpolate.bisplev(xnew[:,0],ynew[0,:],tck)\n",
    "plt.subplot(2,1,2)\n",
    "plt.pcolor(xnew,ynew,znew)\n",
    "plt.colorbar()\n",
    "plt.title(\"Interpolated function.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font color='red'>Integration Problems</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Compute an integral"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%latex\n",
    "We want to compute the integral:\n",
    "    \n",
    "$$ I = \\int_{x=\\pi}^{2\\pi}\\int_{y=0}^{\\pi}{(y\\sin x + x\\cos y)dydx} $$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def integrand(y, x):\n",
    "    'y must be the first argument, and x the second.'\n",
    "    return y * np.sin(x) + x * np.cos(y)\n",
    "\n",
    "ans, err = integrate.dblquad(integrand, \n",
    "                             np.pi, 2*np.pi,  # x limits\n",
    "                   lambda x: 0,               # y limits\n",
    "                   lambda x: np.pi)\n",
    "print(ans)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Solving ODEs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%latex\n",
    "Consider the equation of motion (for the angle $\\theta$) of a damped pendulum:\n",
    "\n",
    "\\begin{eqnarray*}\n",
    "\\frac{d^2\\theta}{dt^2} = -\\frac{1}{Q} \\frac{d\\theta}{dt} + \\sin{\\theta} + d \\cos{\\Omega t}\n",
    "\\end{eqnarray*}\n",
    "\n",
    "It can be written as a system of ODEs:\n",
    "\\begin{eqnarray*}\n",
    "\\left\\{\n",
    "\\begin{array}{rcl}\n",
    "\\frac{d\\theta}{dt} & = & \\omega \\\\\n",
    "\\frac{d\\omega}{dt} & = & -\\frac{1}{Q} \\omega + \\sin{\\theta} + d \\cos{\\Omega t}\n",
    "\\end{array}\n",
    "\\right.\n",
    "\\end{eqnarray*}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use the integrate.odeint function to solve the above system:\n",
    "\n",
    "           integrate.odeint(func, y0, t, args=())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(y, t, params):\n",
    "    theta, omega = y      # unpack current values of y\n",
    "    Q, d, Omega = params  # unpack parameters\n",
    "    derivs = [omega,      # list of dy/dt=f functions\n",
    "             -omega/Q + np.sin(theta) + d*np.cos(Omega*t)]\n",
    "    return derivs\n",
    "\n",
    "# Parameters\n",
    "Q = 2.0          # quality factor (inverse damping)\n",
    "d = 1.5          # forcing amplitude\n",
    "Omega = 0.65     # drive frequency\n",
    "\n",
    "# Initial values\n",
    "theta0 = 0.0     # initial angular displacement\n",
    "omega0 = 0.0     # initial angular velocity\n",
    "\n",
    "# Bundle parameters for ODE solver\n",
    "params = [Q, d, Omega]\n",
    "\n",
    "# Bundle initial conditions for ODE solver\n",
    "y0 = [theta0, omega0]\n",
    "\n",
    "# Make time array for solution\n",
    "tStop = 200.\n",
    "tInc = 0.05\n",
    "t = np.arange(0., tStop, tInc)\n",
    "\n",
    "# Call the ODE solver\n",
    "psoln = integrate.odeint(f, y0, t, args=(params,))\n",
    "\n",
    "# Plot results\n",
    "fig = plt.figure(1, figsize=(8,8))\n",
    "\n",
    "# Plot theta as a function of time\n",
    "ax1 = fig.add_subplot(311)\n",
    "ax1.plot(t, psoln[:,0])\n",
    "ax1.set_xlabel('time')\n",
    "ax1.set_ylabel('theta')\n",
    "\n",
    "# Plot omega as a function of time\n",
    "ax2 = fig.add_subplot(312)\n",
    "ax2.plot(t, psoln[:,1])\n",
    "ax2.set_xlabel('time')\n",
    "ax2.set_ylabel('omega')\n",
    "\n",
    "# Plot omega vs theta\n",
    "ax3 = fig.add_subplot(313)\n",
    "twopi = 2.0*np.pi\n",
    "ax3.plot(psoln[:,0]%twopi, psoln[:,1], '.', ms=1)\n",
    "ax3.set_xlabel('theta')\n",
    "ax3.set_ylabel('omega')\n",
    "ax3.set_xlim(0., twopi)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font color='red'>Optimization Problems</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Find the local maxima of a function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want to compute the local maxima of the Bessel function J(x) given initial guess."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.arange(0,10,0.01)\n",
    "\n",
    "for k in np.arange(0.5,5.5):\n",
    "    y = special.jv(k,x)\n",
    "    plt.plot(x,y)\n",
    "\n",
    "plt.title('Different Bessel functions')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.arange(0,10,0.01)\n",
    "\n",
    "for k in np.arange(0.5,5.5):\n",
    "    y = special.jv(k,x)\n",
    "    plt.plot(x,y)\n",
    "    f = lambda x: -special.jv(k,x)\n",
    "    x_max = optimize.fminbound(f,0,6)\n",
    "    plt.plot([x_max], [special.jv(k,x_max)],'ro')\n",
    "\n",
    "plt.title('Different Bessel functions and their local maxima')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Find the fixed point of a function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def func(x, c1, c2):\n",
    "    return np.sqrt(c1/(x+c2))\n",
    "\n",
    "c1 = np.array([10])\n",
    "c2 = np.array([3])\n",
    "\n",
    "fPoint = optimize.fixed_point(func, [1.2], args=(c1,c2))\n",
    "print(fPoint, func(fPoint, c1, c2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Least Square Problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Parametric function: 'v' is the parameter vector, 'x' the independent varible\n",
    "fp = lambda v, x: v[0]/(x**v[1])*np.sin(v[2]*x)\n",
    "\n",
    "## Noisy function (used to generate data to fit)\n",
    "v_real = [1.5, 0.1, 2.]\n",
    "fn = lambda x: fp(v_real, x)\n",
    "\n",
    "## Error function\n",
    "e = lambda v, x, y: (fp(v,x)-y)\n",
    "\n",
    "## Generating noisy data to fit\n",
    "n = 30\n",
    "xmin = 0.1\n",
    "xmax = 5\n",
    "x = np.linspace(xmin,xmax,n)\n",
    "y = fn(x) + np.random.rand(len(x))*0.2*(fn(x).max()-fn(x).min())\n",
    "\n",
    "## Initial parameter value\n",
    "v0 = [3., 1, 4.]\n",
    "\n",
    "## Fitting\n",
    "v, success = optimize.leastsq(e, v0, args=(x,y), maxfev=10000)\n",
    "\n",
    "## Plot\n",
    "def plot_fit():\n",
    "    print('Estimater parameters: ', v)\n",
    "    print('Real parameters: ', v_real)\n",
    "    X = np.linspace(xmin,xmax,n*5)\n",
    "    plt.plot(x,y,'ro', X, fp(v,X))\n",
    "    plt.legend(['Orig','LeastSq'], loc=4)\n",
    "\n",
    "plot_fit()\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <font color='red'>FFT</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%latex\n",
    "\n",
    "The FFT $X[k]$ of length $N$ of the length-$N$ sequence $x[n]$ is defined as:\n",
    "    \n",
    "    $$ X[k]=\\sum_{n=o}^{N-1}{x[n]e^{-j\\frac{2\\pi}{N}nk}}$$\n",
    "\n",
    "    and the inverse transform is defined as:\n",
    "        \n",
    "    $$ x[n]=\\frac{1}{N}\\sum_{k=o}^{N-1}{X[k]e^{j\\frac{2\\pi}{N}kn}}$$    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# construct signal and plot in the time domain\n",
    "#---------------------------------------------\n",
    "plt.figure(figsize=(6,12))\n",
    "t = np.linspace(0, 1, 1001)\n",
    "y = np.sin(2*np.pi*t*6) + np.sin(2*np.pi*t*10) + np.sin(2*np.pi*t*13)\n",
    "plt.subplot(311)\n",
    "plt.plot(t, y, 'b-')\n",
    "plt.xlabel(\"TIME (sec)\")\n",
    "plt.ylabel(\"SIGNAL MAGNITUDE\")\n",
    "\n",
    "# compute FFT and plot the magnitude spectrum\n",
    "#--------------------------------------------\n",
    "F = fftpack.fft(y)\n",
    "N = len(t)                     # number of samples\n",
    "dt = 0.001                     # inter-sample time difference\n",
    "w = fftpack.fftfreq(N, dt)     # gives us a list of frequencies for the FFT\n",
    "ipos = np.where(w>0)\n",
    "freqs = w[ipos]                # only look at positive frequencies\n",
    "mags = np.abs(F[ipos])         # magnitude spectrum\n",
    "plt.subplot(312)\n",
    "plt.plot(freqs, mags, 'b-')\n",
    "plt.ylabel(\"POWER\")\n",
    "\n",
    "plt.subplot(313)\n",
    "plt.plot(freqs, mags, 'b-')\n",
    "plt.xlim([0, 50])          # replot but zoom in on freqs 0-50 Hz\n",
    "plt.ylabel(\"POWER\")\n",
    "plt.xlabel(\"FREQUENCY (Hz)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise: <font color='red'>Solving a Linear Regression Problem</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want to analyze the US population over the years. \n",
    "In particular, we want to study the percentage of the urban population \n",
    "since the first recorded census.\n",
    "<P>\n",
    "<B> Please connect to the website:</B> \n",
    "    <P>\n",
    "[http://www.u-s-history.com/pages/h980.html](http://www.u-s-history.com/pages/h980.html)\n",
    "<P>\n",
    "<P>\n",
    "We want to construct a linear regression model for the percentage of the US urban population.\n",
    "Here is what we need to do:\n",
    "<UL>\n",
    "<LI> Use the Numpy function <em>genfromtxt</em> to read the file <em>populationUS.txt</em> in order to get the years and urban percentages.\n",
    "<LI> Use the SciPy statistcal function <em>linregress</em> to perform the linear regression.\n",
    "<LI> Do a scatter plot and graph the best fitted line.\n",
    "</UL>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load linRegres_populationUrban.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Modify the file <em>linRegres_population.py</em> to:\n",
    "<UL>\n",
    "<LI> Also read the rural percentages\n",
    "<LI> Perform the linear regression on the rural percentage data\n",
    "<LI> Do the scatter plot and graph the best fitted line.\n",
    "<LI> Determine the year when the two lines intersect.\n",
    "</UL>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}