03-scipy-tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![SciPy](images/scipy_logo.png)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The SciPy framework builds on top NumPy and provides a large number of numerical algorithms for working with data. Some of the topics that SciPy covers are:\n",
    "\n",
    "* Special functions ([scipy.special](http://docs.scipy.org/doc/scipy/reference/special.html))\n",
    "* Integration/ODEs ([scipy.integrate](http://docs.scipy.org/doc/scipy/reference/integrate.html))\n",
    "* Optimization ([scipy.optimize](http://docs.scipy.org/doc/scipy/reference/optimize.html))\n",
    "* Interpolation ([scipy.interpolate](http://docs.scipy.org/doc/scipy/reference/interpolate.html))\n",
    "* Fourier Transforms ([scipy.fftpack](http://docs.scipy.org/doc/scipy/reference/fftpack.html))\n",
    "* Signal Processing ([scipy.signal](http://docs.scipy.org/doc/scipy/reference/signal.html))\n",
    "* Linear Algebra ([scipy.linalg](http://docs.scipy.org/doc/scipy/reference/linalg.html))\n",
    "* Sparse Eigenvalue Problems ([scipy.sparse](http://docs.scipy.org/doc/scipy/reference/sparse.html))\n",
    "* Statistics ([scipy.stats](http://docs.scipy.org/doc/scipy/reference/stats.html))\n",
    "* Multi-dimensional image processing ([scipy.ndimage](http://docs.scipy.org/doc/scipy/reference/ndimage.html))\n",
    "* File IO ([scipy.io](http://docs.scipy.org/doc/scipy/reference/io.html))\n",
    "\n",
    "\n",
    "In many cases, you will want to import specific names from `scipy` subpackages. However, as a start, it is helpful to do the following import:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy as sp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Approach"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One of the most important skills in data science is to be able to find Python functions and classes in a module and learn how to use them yourself. Here are some recommended steps on how to go about this:\n",
    "\n",
    "* Find the online documentation for the package you are using.\n",
    "* Try to find the subpackage or even the function that looks like will do the job.\n",
    "* Import the module, function or class and use tab completion and `?` to explore it.\n",
    "* Try using the function or class for an extremely simple case where you know the answer.\n",
    "* Then try using for your real problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data Fitting\n",
    "\n",
    "Let's randomize a exponential decay function and try to do a linear fit.\n",
    "\n",
    "We have\n",
    "\n",
    "$$D(t) = N e^{(-t/\\tau)}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Define the function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def D(t, N, tau):\n",
    "    return N * np.exp(-t/tau)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* randomize data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ndata = 50\n",
    "tdata = np.linspace(0., 5., ndata)\n",
    "data_pure = D(tdata, N=7., tau=3.)\n",
    "data_measure = data_pure + 0.3 * np.random.normal(size=ndata)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(tdata, data_pure)\n",
    "plt.plot(tdata, data_measure)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Look for fit scipy module and fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import curve_fit\n",
    "\n",
    "popt, pcov = curve_fit(D, tdata, data_measure)\n",
    "\n",
    "print (popt)\n",
    "N_fit, tau_fit = popt\n",
    "data_fit = D(tdata, N=N_fit, tau=tau_fit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* plot result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(tdata, data_pure)\n",
    "plt.plot(tdata, data_measure)\n",
    "plt.plot(tdata, data_fit, label=\"fit\")\n",
    "\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "##   <span style=\"color:blue\">Exercise: simple integration</span>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate \n",
    "\n",
    "$$I = \\int_0^\\infty e^{-x} dx$$\n",
    "\n",
    "hint: [GIYF](https://www.google.es/search?q=scipy%20integrate%20function&rct=j)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# %load files/sol_integ.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# integrate tab-completion / help example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Real life example: solar activity period"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Adapted from [Python for scientific use, Part II: Data analysis](http://linuxgazette.net/115/andreasen.html) \n",
    "by Anders Andreasen. See also a more modern approach to the subject using [Pandas](http://pandas.pydata.org/),  the Python Data Analysis Library, on [this Josh Hemann notebook](http://nbviewer.ipython.org/gist/jhemann/4569783)\n",
    "\n",
    "In the next example we will use the Fast Fourier Transform (FFT) in order to transform time-dependent data into the frequency domain. By doing so, it is possible to analyse if any predominant frequencies exists - i.e. if there is any periodicity in the data. \n",
    "\n",
    "The example data we will use is the sunspot activity measurements from year 1700 to year 2016 provided by the Sunspot index and Long-term Solar Observation. The data set is the yearly mean total sunspot observations available via HTTP.\n",
    "\n",
    "The datafile is also included in `files/SN_y_tot_V2.0.txt`. The sunspot data have been used for illustrating the power of FFT with respect to finding a periodicity in sunspot activity in various computer languages, e.g. Matlab and BASIC. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Retrieve data\n",
    "\n",
    "Solar spot data are available online from the [datafiles section](http://www.sidc.be/silso/datafiles) of the Sunspot index and Long-term Solar Observation webpage. We will use the yearly mean total sunspot number [plain text version.]( versionhttp://www.sidc.be/silso/DATA/SN_y_tot_V2.0.txt)\n",
    "\n",
    "A quick web search point us to the `urllib` python module and to a very simple retrival procedure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import urllib\n",
    "url = \"http://www.sidc.be/silso/DATA/SN_y_tot_V2.0.txt\"\n",
    "raw_data = urllib.request.urlopen(url)\n",
    "solarspots = np.loadtxt(raw_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "years = solarspots[:,0]\n",
    "spots = solarspots[:,1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10,6))\n",
    "\n",
    "plt.plot(years, spots)\n",
    "plt.plot(years, spots, \"ro\", label=\"spots\" )\n",
    "\n",
    "plt.xlabel(\"year\")\n",
    "plt.ylabel(\"yearly sunspots\")\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "\n",
    "#plt.savefig(\"yearly_sunsports.png\", dpi=300)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### FFT\n",
    "\n",
    "Let's take a look at the [SciPy Fast Fourier Transform module manual](http://docs.scipy.org/doc/scipy/reference/tutorial/fftpack.html), and straight from there...\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.fftpack import fft, fftfreq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spots_fft = fft(spots)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(spots_fft.real, spots_fft.imag)\n",
    "\n",
    "#plt.xlim(-5000,5000)\n",
    "plt.xlabel(\"real(FFT)\")\n",
    "plt.ylabel(\"imag(FFT)\")\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to construct a periodogram, i.e. a graph of power vs. frequency, we first compute the power of the FFT signal which is simply the **FFT signal squared**. We only need the part of the signal ranging from zero to a frequency equal to the Nyquist frequency, which is equal to half the maximum frequency, since frequencies above the Nyquist frequency correspond to negative frequencies. The frequency range is calculated from 0-N/2 as N/(2T) where N is the number of samples and T is the sampling time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = spots_fft.size\n",
    "timestep = 1 # year\n",
    "freq = fftfreq(N, d=timestep)   # 1/year\n",
    "power = np.abs(spots_fft)**2   # proportional to this\n",
    "\n",
    "freq = freq[1:N/2]\n",
    "power = power[1:N/2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(freq, power)\n",
    "\n",
    "plt.xlabel(\"Frequency [1/year]\")\n",
    "plt.ylabel(r\"$|FFT|^2$\")\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Which in years becomes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "period = 1./freq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(period, power)\n",
    "\n",
    "plt.xlabel(\"Period [years]\")\n",
    "plt.ylabel(r\"$|FFT|^2$\")\n",
    "plt.grid()\n",
    "\n",
    "#plt.xlim( (0,30) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solar spots period"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "max_power_index = np.argmax(power)\n",
    "max_power_period = period[max_power_index]\n",
    "print ('Maximum solar spots period is %1.2f years' %(max_power_period) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Scipy direct method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After doing the exercise, while looking for information, I found that scipy allows us to directly compute the periodgram with [`scipy.signal.periodgram`](http://scipy.github.io/devdocs/generated/scipy.signal.periodogram.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import signal\n",
    "frec=1.  # sampling frequency: 1 year\n",
    "f, Pxx_spectrum = signal.periodogram(spots, frec, scaling='spectrum')\n",
    "plt.plot(f, Pxx_spectrum)\n",
    "plt.xlim(.0,.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "max_power_index = np.argmax(Pxx_spectrum)\n",
    "max_power_period = 1./f[max_power_index]\n",
    "print ('Maximum solar spots period is %1.2f years' %(max_power_period) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Resources"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* [SciPy Website](http://www.scipy.org)\n",
    "* [SciPy Reference Documentation](http://docs.scipy.org/doc/scipy/reference/)\n",
    "* [Python Scientific Lecture Notes](http://scipy-lectures.github.io/index.html), Edited by Valentin Haenel,\n",
    "Emmanuelle Gouillart and Gaël Varoquaux.\n",
    "* [Lectures on Scientific Computing with Python](https://github.com/jrjohansson/scientific-python-lectures), J.R. Johansson.\n",
    "* [Introduction to Scientific Computing in Python](http://nbviewer.ipython.org/github/jakevdp/2014_fall_ASTR599/tree/master/), Jake Vanderplas."
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Raw Cell Format",
  "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}