03 - Numpy Essential - instructor-with answer.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Numpy and Matplotlib essential"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As suggested in one of the previous challenges, the *numpy* library is providing an object called *array* which is very similar object to the list. The main difference is the set of operations that can be performed on them. Numpy arrays are oriented towards computation. \n",
    "\n",
    "In this lesson we are going to learn a little more about how to use this fundamental library to do any numerical analysis in Python.\n",
    "\n",
    "Numpy extends the Python language by providing new types (array, matrix, masked_array...), functions and methods to realise efficient numerical calculation using Python.\n",
    "\n",
    "As we saw in the Introduction course, the most basic numpy type is called an *array*. The most basic array is a multi-dimensional object which contains numerical data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Numpy](http://www.numpy.org/) is **the** numerical library for Python. This library is too big to be taught in a day. That will be extremely boring and not useful at all. \n",
    "\n",
    "This library is at the base of all the other python libraries used in science and data science: \n",
    "\n",
    "- [scipy](https://www.scipy.org/) Fundamental library for scientific computing (interation, optimisation...)\n",
    "- [pandas](http://pandas.pydata.org/) data structure and data analysis tools\n",
    "- [matplotlib](https://matplotlib.org/) Python 2D plotting\n",
    "\n",
    "And more specialised one like:\n",
    "\n",
    "- [astropy](http://www.astropy.org/) for the astronomy\n",
    "- [h5py](https://www.h5py.org/) to interact with the HDF5 format\n",
    "- [scikit-learn](http://scikit-learn.org/) for Machine Learning\n",
    "- [TensorFlow](https://www.tensorflow.org) for Deep Learning\n",
    "\n",
    "We are going to learn some of the basic commands not seen in the first course. We will introduce some of these libraries but keep in mind that we are just covering the basics so you can understand how to start using this library. If you find that your interest is piqued and/or they are applicable to your problem, you can use the documentation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To start, we are going to import the two libraries *numpy* and *matplotlib* that will be used in this episode."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pylab as plt\n",
    "\n",
    "# Only on jupyter notebook\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "odds_np = np.arange(1,10,2)\n",
    "\n",
    "print('A numpy array:', odds_np)\n",
    "print('type of the object:', type(odds_np))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**note**\n",
    "\n",
    "'numpy.arange' is very similar to python's built in 'range' function but it returns a numpy array instead of an iterator."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style='background:#B1E0A8; padding:10px 10px 10px 10px;'>\n",
    "<H2> Challenges </H2>\n",
    "\n",
    " <ol>\n",
    " <li> \n",
    " What is the result of adding two lists together?\n",
    " </li>\n",
    " <li>\n",
    " What happens when you add two numpy arrays?\n",
    " </li>\n",
    " <br>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "odds = [ i for i in range(10) if (i % 2) ]\n",
    "print(\"list of odds:\", odds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "odds+odds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "odds_np + odds_np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Work with matrix like object\n",
    "\n",
    "*Numpy* arrays are not exactly matrices. They are multidimensional objects to store data and they are not necessarily numerical. In contrast, a *matrix* object has to be numerical. In practice, a 2D numerical numpy array is very similar to a numpy matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style='background:#B1E0A8; padding:10px 10px 10px 10px;'>\n",
    "<H2> Challenges </H2>\n",
    "\n",
    " <ol>\n",
    " <li> \n",
    " Create a numpy array which contains odd numbers greater than 1 and less than 20. \n",
    " </li>\n",
    "    <li>\n",
    "        Then create an array containing only the first 9 elements from the first array.\n",
    "    </li>\n",
    " <br>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "odds_np = np.arange(1,20,2)\n",
    "\n",
    "arr = odds_np[:9]\n",
    "print(arr)\n",
    "print(arr.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reshaping arrays\n",
    "\n",
    "An array shape can be modified using the method *reshape*. Here we are creating a square array 3x3 from the previous array *arr*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "arr = arr.reshape((3,3)) # notice the usage of a tuple (3,3).\n",
    "print('The square array:\\n', arr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we have used a new object called a tuple. A tuple is very similar to a list but it is designed to be *immutable*. (If you remember from the previous course this means that individual items cannot be modified.)\n",
    "\n",
    "We can verify the type of the array created above and verify that it is, indeed, a numpy array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "type(arr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can check the type of the element which composed that array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(arr.dtype)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is possible to modify the type of the element by using the method *astype*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "arr = arr.astype(float)\n",
    "print(arr.dtype)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An array can be converted into a matrix object. A numpy matrix does have the characteristic of the mathematical matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mat = np.matrix(arr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "type(mat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A multiplication of matrix with the operator *\\** will provide the expected mathematical answer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mat * mat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style='background:#B1E0A8; padding:10px 10px 10px 10px;'>\n",
    "<H2> Challenges </H2>\n",
    "\n",
    " <ol>\n",
    " <li> \n",
    " What happens when you multiply two arrays together? How is it different to multiplying matrices?\n",
    " </li>\n",
    " \n",
    " <br>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "arr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "arr * arr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Matrix multiplication is also possible for arrays -- use the *matmul* function from numpy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.matmul(arr, arr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.matmul(mat, mat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Masked array"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One very useful tool that numpy provides to analyse experimental data is the masked array. When you are taking data from an experiment you always have some data which are not present, or with a bad signal to noise ratio, or that you cannot use for some other reason.\n",
    "\n",
    "Masked arrays associate a numpy array with another array composed only of boolean values (True or False). These tell numpy whether to use (or not) the respective element.\n",
    "\n",
    "To demonstrate this we are going to create a Gaussian function and use it to generate an example dataset and generate a plot. We will then add some noise to it and use a masked array to filter out the noisy data.\n",
    "\n",
    "Reminder: the Gaussian function is define by:\n",
    "\n",
    "\\begin{equation*}\n",
    "g(x) = \\frac{1}{\\sqrt{2\\pi\\sigma^{2}}} e^{-\\frac{(x-\\mu)^2}{2\\sigma^2}}\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style='background:#B1E0A8; padding:10px 10px 10px 10px;'>\n",
    "<H2> Challenges </H2>\n",
    "\n",
    " <ol>\n",
    " <li> \n",
    " Create a function called `gauss` which will take three arguments (inputs): *x*, *mu*, and *sigma*, as defined above. \n",
    "     \n",
    "(x is an array, mu is the position of the centre of the curve/peak and sigma is the width of the bell)\n",
    " </li>\n",
    " <li>\n",
    " Create a numpy array using the 'numpy' function 'linspace' which will contain 1000 points equally spaced between x=-100 and x=100.<br>\n",
    " Hint: You can print the help documentation of a function with 'help(name_of_the_function)'\n",
    " </li>\n",
    " <li>\n",
    " Using the above gauss function and the array, create a list  which contains the value of the gauss from x=-100 to x=100.\n",
    " </li>\n",
    " <li>\n",
    " Use the 'matplotlib' library to plot the curve with mu=0 and sigma=10. \n",
    " </li>\n",
    " <br>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Challenge 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gauss(x, mu=0, sigma=1):\n",
    "    return (1./(sigma * np.sqrt(2 * np.pi)) * \n",
    "            np.exp( - (x - mu)**2 / (2 * sigma**2)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Challenge 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(-100, 100, 1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Challenge 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu = 0   # the position of the center of the peak \n",
    "sigma = 10  # # The width of the 'bell'\n",
    "g = gauss(x, mu, sigma)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Challenge 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(x, g)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Noisy signal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we are going to add some random noise to that curve. \n",
    "To do it we can use the numpy function *normal* from the module *random* provided by numpy library:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = 100*g   # To have something visible we are multiplying the function by 100.\n",
    "\n",
    "# Creation of the noise\n",
    "noisy = np.random.normal(g)\n",
    "plt.plot(x, g+noisy)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A way to calculate the 'Signal to Noise' ratio of the previous data set is by dividing the noisy data by the standard deviation of the noisy data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get standard dev of noisy data\n",
    "rms = np.std(noisy)\n",
    "# Divide noisy data by std dev\n",
    "SN = noisy / rms\n",
    "# Plot\n",
    "plt.plot(SN)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will next create a mask that masks the data where the SN (signal to noise) is less than 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask = SN < 1\n",
    "print('Signal to noise:',SN[:10])\n",
    "print('Mask:',mask[:10])\n",
    "print('Mask shape:',mask.shape)\n",
    "noisy_ma = np.ma.array(noisy, mask=mask)\n",
    "plt.plot(noisy_ma)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Working with images"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are going to learn some commands that deal with images. \n",
    "Since most scientific domains use their own file format, we obviously can not learn all of them.\n",
    "We will use a typical astronomical image format: the 'fits' file. \n",
    "\n",
    "In the *data* directory you should find a file called *502nmos.fits*. \n",
    "\n",
    "We can verify that the file is indeed here:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cd data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ls *.fits"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To be able to manipulate data in this file, we need to import a library which will be able to open the file and put the data in a numpy array. \n",
    "\n",
    "This is a good occasion to install a new library. To realise that, open a terminal. On Microsoft Windows, can you start the **Anaconda prompt** software: \n",
    "\n",
    "![MS Windows terminal](images/anaconda-prompt.png \"Anaconda Prompt terminal\")\n",
    "\n",
    "On Unix system, you can start a terminal or having one started through the Jupyter notebook:\n",
    "\n",
    "![Starting a Terminal with Jupyter](images/jupyter_terminal.png \"Starting a terminal with Jupyter Notebook\")\n",
    "\n",
    "![Terminal using Jupyter](images/jupyter_terminal2.png \"Terminal with Jupyter Notebook\")\n",
    "\n",
    "We are going to use the command **pip** which will allows you to install any package available on the python software repository: [Pypi](https://pypi.python.org/pypi). Here we will install the library [pyfits](http://www.astropy.org/) which will provide us the tools to open the fits images files. \n",
    "\n",
    "**Note:**\n",
    "\n",
    "This is only an example to teach you how to install a library not yet available and how to open a specific image format. I am not expecting you to use this specific library in the furutre. Each domain, has their own file format (most of the time): \n",
    "\n",
    "- [NetCDF](https://www.unidata.ucar.edu/netcdf/)\n",
    "- [HDF5](https://www.hdfgroup.org/HDF5/)\n",
    "- MS Excel\n",
    "- SQL\n",
    "- ...\n",
    "\n",
    "All of these formats can be open using python but you will have to install an additional library to do it.\n",
    "\n",
    "In the previously open terminal:\n",
    "\n",
    "```bash\n",
    "\n",
    "pip install -U pyfits --user\n",
    "```\n",
    "\n",
    "\n",
    "The *pip* command will look at the *pypi* repository if the library *pyfits* is available, the *-U* option will look if there are an upgrade available (if *pyfits* is already present on your system) and the *--user* option will install the library in the user account not on the system wide (you do not have to be administrator to install a new python library)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "<div style='background:#B1E0A8; padding:10px 10px 10px 10px;'>\n",
    "<H2> Challenges </H2>\n",
    "\n",
    " <ol>\n",
    " <li>\n",
    " What is the main format for your data?\n",
    " </li>\n",
    " <li>\n",
    " Find a python library which will allows you to open this format and convert it in a numpy array.\n",
    " </li>\n",
    " </div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some available libraries:\n",
    "\n",
    "- FITS, VOTable...: [astropy](http://docs.astropy.org/en/stable/index.html#files-i-o-and-communication)\n",
    "- CSV, HDF5, MS Excel, SQL...: [pandas](http://pandas.pydata.org/pandas-docs/stable/io.html)\n",
    "- NetCDF: [netcdf4](http://unidata.github.io/netcdf4-python/)\n",
    "- Matlab mat: [scipy](https://docs.scipy.org/doc/scipy/reference/tutorial/io.html)\n",
    "- Tiff: [Pillow](https://pillow.readthedocs.io/en/latest/), [pylibtiff](https://github.com/pearu/pylibtiff),[matplotlib]()..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are going to the data directory:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "#Windows\n",
    "#os.chdir('C:\\Work\\Python-advanced\\data')\n",
    "#os.chdir('data')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.getcwd()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Working with binary file\n",
    "\n",
    "We will import the new library and we can now open the fits file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#import astropy.io.fits as pyfits\n",
    "\n",
    "try: \n",
    "    import pyfits\n",
    "except ImportError:\n",
    "    import astropy.io.fits as pyfits\n",
    "    \n",
    "im1 = pyfits.open('data/502nmos.fits')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note:**\n",
    "\n",
    "By default pyfits opens a file with the option *memmap=True*. \n",
    "This option opens the fits file without copying the data into memory and allows us to open very large files which will not fit in the physical memory."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fits file are composed of a list of HDUs (Header and data units). We can list the information with the method *info*.\n",
    "\n",
    "Here we will be only interested on the primary HDU which is an image and is called *PRIMARY*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "im1.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can check what is the type of the image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "type(im1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note**\n",
    "HDUList is a new object that you probably are not yet familiar with."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can access the first HDU (the one of interest) which contains the image as an element in a list using the index:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "HDU_copy1 = im1[0]     # As mentioned an HDU is a list and we want to have access to the first element."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is also possible to access it by using it's name (or key) as with a dictionary: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "HDU_copy2 = im1['PRIMARY']  # HDU can also be accessed through their name like dictionary"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can verified that both copy are the same."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "HDU_copy1 == HDU_copy2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fits data, as with most scientific data images containers, are comprised of a header and a data part. The header contains metadata relevant to the observations and to the data itself. You can print them on screen by using the attribute *header*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "im1[0].header"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To access the data itself (here an image of a nebulae), we access the *data* attribute."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "imdata = im1[0].data\n",
    "print(type(imdata))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use matplotlib to show what the image look like"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(imdata, origin='lower', cmap='gray')\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, it is difficult to see anything on this image. This happens very often in astronomical images where very bright objects are saturating the CCD and a linear output will show a limited number of pixels due to the contrast. \n",
    "\n",
    "We are going to improve the visible output by doing some simple analysis on the image which will help to side step the contrast problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('mean value im1:', imdata.mean())\n",
    "print('median value im1:', np.median(imdata))  # Note: median is not provided by a method!\n",
    "print('max value im1:', imdata.max())\n",
    "print('min value im1:', imdata.min())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The previous results give us some useful information. We can see that the range between the minimum and maximum value is really big and also that the maximum value is very far from the median or median value in a pixel. That probably means that a very small number of pixels have very high values.\n",
    "\n",
    "We can check it by plotting a histogram which will give us the number of pixels per number of photons.\n",
    "\n",
    "**Note**\n",
    "The ravel() function flattens our 2x2 image data (where x and y were pixel co-ordinates and the values are photon counts) \n",
    "...and turns it into a flat 1 dimensional array of photon counts (again: one value for each pixel but we drop the pixel co-ordinate information).  \n",
    "See imdata.ravel?, help(imdata.ravel) or numpy.ravel for more information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hist = plt.hist(imdata.ravel(), bins=100)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The histogram confirms our suspicion: many pixels have very low photons.\n",
    "\n",
    "When we have such contrasting data, there are several ways we can improve the visualisation. \n",
    "\n",
    "For example:\n",
    "\n",
    "1) Plotting the y-axis (counts) in a logarithm.\n",
    "\n",
    "**Note**\n",
    "Plotting a scale using a logorithm is a well known way of reducing contrast in the data. For example scientists using earth quake data, sound level data and natural science dat often use this technique.\n",
    "See: https://en.wikipedia.org/wiki/Logarithmic_scale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hist = plt.hist(imdata.ravel(), bins=100)\n",
    "plt.yscale('log')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2) We can also manually modify the upper limit in number of pixels with photons (y axis) using the method ylim:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hist = plt.hist(imdata.ravel(), bins=100)\n",
    "plt.ylim(0,1e5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3) It is also possible to limit the range in the number of photons in a pixel (x axis) using the *range* argument in the *hist* function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hist = plt.hist(imdata.ravel(), bins=100, range=(1,30))\n",
    "plt.ylim(0,1e5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4) Using the previous graphic we can limit the range of photons when plotting the previous image and improve the contrast using the *vmax* argument: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(imdata, origin='lower', cmap='gray', vmax=25)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next plot does not provide a lot of useful information in our particular analysis but we include it as it is a classical plot that you can obtain with visualisation tools. \n",
    "\n",
    "Axis 0: is the sums of the pixels in the rows (axis=0) of the images, divided by the number of rows ('imdata.shape[0]').\n",
    "Axis 1: is a sum the columns of the images and divides them by the number of columns. \n",
    "\n",
    "With this plot you can distinguish the 2 bright objects (stars) located in the lower left part of the image and the brightest, in the more central location."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(imdata.sum(axis=0)/imdata.shape[0], label='axis0')\n",
    "plt.plot(imdata.sum(axis=1)/imdata.shape[1], label='axis1')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "5) Another option to improve the contrast visually is to plot the image after a logarithm conversion (that will flatten the image). The problem is that the image can have negative value. \n",
    "\n",
    "**NOTE** \n",
    "The  *imdata - imdata.min() + 1* assures us that every value in the image will be strictly greater than 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "logim = np.log(imdata - imdata.min() + 1)\n",
    "plt.imshow(logim, origin='lower', cmap='gray')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Final method: Applying a masked array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a dataset that has values greater than 25 masked out.\n",
    "immasked = np.ma.masked_greater(imdata, 25)\n",
    "\n",
    "# We could also do it this way, where we create the mask explicitly\n",
    "# and then apply it.\n",
    "#mask = imdata > 25\n",
    "#masked1 = np.ma.array(imdata, mask=mask)\n",
    "\n",
    "print(immasked)\n",
    "print(immasked.mask)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Matplotlib is aware of this numpy object and will plot the image by looking **only** at the pixels with the value **True** in the associated mask:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(immasked, cmap='gray', origin='lower')\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To access the mask created above, it is possible using the attribute **mask**\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(immasked.mask, cmap='gray', origin='lower')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('original average:', imdata.mean())\n",
    "print('Masked average:', immasked.mean())\n",
    "print()\n",
    "print('original max:', imdata.max())\n",
    "print('Masked max:', immasked.max())\n",
    "print()\n",
    "print('original min:', imdata.min())\n",
    "print('Masked min:', immasked.min())\n",
    "print()\n",
    "print('original median:', np.mean(imdata))\n",
    "print('Masked median:', np.ma.median(immasked))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style='background:#B1E0A8; padding:10px 10px 10px 10px;'>\n",
    "<H2> Challenges </H2>\n",
    "\n",
    " <ol>\n",
    " <li>\n",
    " Create a new masked array 'immasked2' where there is a lower limit set to zero\n",
    "     (Similar to immasked = np.ma.masked_greater(imdata, 25) above, but call the 'masked_less_equal' function.\n",
    " </li>\n",
    "  <li>\n",
    " Update the immasked2.mask to apply both masks (ignore values greater than 25 OR (use the '|' operator for 'or') less than or equal to zero). Show this image.\n",
    " </li>\n",
    " </div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "immasked2 = np.ma.masked_less_equal(imdata, 0)\n",
    "\n",
    "immasked2.mask = immasked2.mask | immasked.mask\n",
    "\n",
    "plt.imshow(immasked2, cmap='gray', origin='lower')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Dtype"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "We have already come across the 'dtype' which stands for data type. \n",
    "This not only provides useful information but is also a very powerful tool, especially when loading a data file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cd data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using this tool we can load a complex file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%more informations.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('informations.txt') as f:\n",
    "    data = f.readlines()\n",
    "\n",
    "for line in data:\n",
    "    print(line)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "data = np.loadtxt('informations.txt', \n",
    "                  skiprows=1, \n",
    "                  dtype={'names':('firstname', 'surname', 'age', 'phone'), \n",
    "                         'formats':('S8','S8','i','i')})\n",
    "\n",
    "\n",
    "print(data)\n",
    "print(type(data))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What is the shape of this numpy array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can access to each row using the indexes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And you can access to th different columns by using the keyword as in a dictionary:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data['age']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data['age'][0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(data['firstname'][0], data['phone'][0])"
   ]
  }
 ],
 "metadata": {
  "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}