S2_TD3_python_notebook_en_student.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# <span style=\"color:blue\"> <center>Student / TD3 : 4TPU279U $-$ Bachelor 1st year $-$ spring 2023</center></span>\n",
    "# <center>Introduction to python programming</center>\n",
    "# <hr style=\"border:1px solid black\"><center>  Data Processing and Plot </center><hr style=\"border:1px solid black\">\n",
    "</br>\n",
    "\n",
    "<div style=\"text-align: right\"> Credits: R. Boisgard, L. Truflandier, Philippe Paillou, Julien Burgin, Sara Zein, Leo Delmarre, Simon Villain-Guillot </div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr style=\"border:1px solid grey\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following topics will be cover in this notebook:\n",
    "- Process data files\n",
    "- Basics in plotting\n",
    "- Survival guide for practical work"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <hr style=\"border:1px solid black\">  Process data files  <hr style=\"border:1px solid black\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using python intrinsic functions\n",
    "> Having a data file **located in the directory** `data_examples` located at the same place than the jupypter-notebook file the later can be *opened* using the built-in function [`open()`](https://docs.python.org/3/library/functions.html#open).\n",
    "> An example is given below with the data file `data_sin.csv` where the second argument in `open()` -namely `'r'`- stands for *reading* processing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "file = open('data_examples/data_sin.txt','r')\n",
    "print(type(file))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Remark: to visualize the content of `data_sin.txt` in a jupyter-notebook you can clic on `File/Open` and choose the file from the current directory*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> the content of the `data_sin.txt` can be imported using the `_io.TextIOWrapper`\n",
    "function `readlines()`. In that case data  is stored in `data` as a list of strings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = file.readlines()\n",
    "print(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> To get ride of the first line we can use `readline()` (without s!) first. This will extract the first line of the file. Then `readlines()` can be safely used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "file = open('data_examples/data_sin.txt','r')\n",
    "data = file.readline()\n",
    "data = file.readlines()\n",
    "print(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> The character string `\\n` indicates that the current line ends at that point and a new line starts right after it. This can be examplified by using `print()`. For instance:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Hello\\nworld!')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Having inspected `data_sin.txt` we know that we have 2 columns of data and that the first line refer to the type of data. Unfortunatly, `data` is a single list containing concatenated column elements."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> In order to separate column elements, we may recognize the space `' '` as separator. Note that `'.'` is also a seprator for identifying decimal numbers. To separate the columns we can re-process and operate *line by line* by looping over the rows of data and use the function `split`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = [ ] ; y = [ ]\n",
    "for line in data:\n",
    "    words = line.split(' ')\n",
    "    #print(words)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Since `words` is a list of 2 elements it becomes easy to store data into two **lists** $x$ and $y$ and convert strings to floats."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = [ ] ; y = [ ]\n",
    "for line in data:\n",
    "    words = line.split(' ')\n",
    "    x.append( float(words[0]) )\n",
    "    y.append( float(words[1]) )\n",
    "print(x)\n",
    "print()\n",
    "print(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> To summarize how to process a data file containing a header and two columns of data\n",
    "separated by a space, we have:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "file = open('data_examples/data_sin.txt','r')\n",
    "file.readline()\n",
    "\n",
    "x = [ ] ; y = [ ]\n",
    "for line in file:\n",
    "    words = line.split(' ')\n",
    "    x.append( float(words[0]) )\n",
    "    y.append( float(words[1]) )\n",
    "#print(x)\n",
    "#print()\n",
    "#print(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> For later use, note that if we want to operate on $x$ and $y$ using mathematical operators we will need to convert them into `numpy array`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import array\n",
    "x = array(x)\n",
    "y = array(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Note that even the command `open` can also be used for writing data but is not well adapted to handle list and arrays. We should prefer to use the function `savetxt` of the numpy libray. An example is given below*\n",
    "\n",
    "For instance if you want to create a file containing 2 columns $(\\theta,f)$ with $f(\\theta)=\\cos(x)$ with $\\theta\\in[0,\\pi]$ for 32 data points, you can use:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import pi, cos, linspace, savetxt, column_stack, sin\n",
    "\n",
    "theta = linspace(0.0,2*pi,64)\n",
    "f     = cos(theta)\n",
    "\n",
    "savetxt('data_examples/data_cos.txt', column_stack([theta,f]), fmt='%1.4f')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Note on CSV file\n",
    "A comma-separated values ([CSV](https://en.wikipedia.org/wiki/Comma-separated_valuesCSV)) file is a delimited text file that uses a comma to separate values. Each line of the file is a data record. A CSV file typically stores tabular data (numbers and text) in plain text, in which case each line will have the same number of fields. Alternative delimiter-separated files are often given a \".csv\" extension despite the use of a non-comma field separation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Exercice 3.0 :  English CSV format  <hr style=\"border:1px solid grey\">\n",
    "1. Open the data file `data_sin_en.csv` located in `data_examples` directory.\n",
    "2. Extract the 2 columns and store the results into 2 arrays. Note that the separator is a coma `','`.\n",
    "3. Print the results."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr style=\"border:1px solid grey\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Exercice 3.1 : French CSV format  <hr style=\"border:1px solid grey\">\n",
    "1. Open the data file `data_sin_fr.csv` located in `data_examples` directory.\n",
    "2. Extract the 2 columns and store the results into 2 arrays. \n",
    "3. Print the results.\n",
    "*Note that this is a French CSV data file. In this case the separator is the semi-column character \";\". The decimal separator is the coma \",\". To store the data into array you will need to replace the\n",
    "coma by a dot. This can be done with the function `replace()`.* An example is given below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test = '3,14116'\n",
    "test.replace(',','.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr style=\"border:1px solid grey\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using Pandas"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Pandas](https://en.wikipedia.org/wiki/Pandas_(software)) is a library written for Python programming language for data manipulation and analysis. In particular, it offers data structures and operations for manipulating numerical tables and time series which are refered as to *DataFrames*.\n",
    "\n",
    "[Pandas](https://pandas.pydata.org) is mainly used for data analysis and associated manipulation of tabular data in DataFrames. Pandas allows importing data from various file formats such as CSV, JSON, Parquet, SQL database tables or queries, and Microsoft Excel. The pandas library is built upon another library NumPy, which is oriented to efficiently working with arrays."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Read and process `CSV` file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pandas import read_csv\n",
    "\n",
    "data = read_csv('data_examples/harelynx.csv')\n",
    "print(type(data))\n",
    "print(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Read and process `EXCEL` file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pandas import read_excel\n",
    "\n",
    "data = read_excel('data_examples/harelynx.xls')\n",
    "print(type(data))\n",
    "print(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> From `panda DataFrame` to `numpy array`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_numpy_array = data.to_numpy()\n",
    "print(type(data_numpy_array))\n",
    "#print(data_array)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> From `numpy array` to `pandas DataFrame`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pandas import DataFrame\n",
    "\n",
    "data_pandas = DataFrame(data_numpy_array,columns=['t', 'lievre', 'lynx'])\n",
    "print(data_pandas)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Print `pandas DataFrame` to `CSV` file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# with the index and the header (default)\n",
    "data_pandas.to_csv('data_examples/lievrelynx_hi.csv') \n",
    "\n",
    "# without the index and the header\n",
    "data_pandas.to_csv('data_examples/lievrelynx.csv', header=False, index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Basic statistics with `pandas`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(data_pandas['lievre'].describe())\n",
    "print()\n",
    "print(data_pandas['lynx'].describe())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Exercice 3.2 :   <hr style=\"border:1px solid grey\">\n",
    "Given a set of data $\\{x\\}_{i=1}^{N}$, the mean ($\\mu$) and the standard deviation ($\\sigma$) are given by:\n",
    "$$\\mu = \\frac{1}{N}\\sum_{i=1}^N x_i$$ and $$ \\sigma=\\sqrt{\\frac{1}{N-1}\\sum_{i=1}^N(x_i-\\mu)^2}$$\n",
    "> Load the data file `lynxhare.csv` using `pandas` and convert to a numpy array.\n",
    "\n",
    "> Generate two 1D numpy arrays called `hare` and `lynx`.\n",
    "\n",
    "> Using the python built-in function (`min`, `max`, `sum`, `len`) compute the min/max, $\\mu$ and $\\sigma$ for the hare and lynx population. \n",
    "\n",
    "> Print the results with an appropriate format and compare to `pandas` values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy  import sqrt\n",
    "from pandas import read_csv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr style=\"border:1px solid grey\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <hr style=\"border:1px solid black\">  Basic in plotting  <hr style=\"border:1px solid black\">\n",
    "Below some examples how to do plots with python with an increasing complexity. A huge gallery of examples can be found here [https://matplotlib.org/stable/gallery/index.html](https://matplotlib.org/stable/gallery/index.html). 1D plots can be drawn using the (sub)-library \n",
    "[`pyplot`](https://matplotlib.org/stable/tutorials/introductory/pyplot.html) from [`matplotlib`](https://matplotlib.org)\n",
    "\n",
    "> Note that the alias `plt` shall be used througout this notebook for the `pyplot`.\n",
    "respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Simple 1D graphics\n",
    "\n",
    "> Let's plot $y=\\sin(2\\pi t)\\exp(-t)$ for $t\\in[0,5]$ with $n=40$ grid points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import pi, exp, cos, linspace\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "n    =  40\n",
    "tmin =   0\n",
    "tmax =   5\n",
    "\n",
    "t = linspace(tmin,tmax,n)\n",
    "y = cos(2*pi*t)*exp(-t)\n",
    "\n",
    "plt.plot(t,y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> By default, the function [`plot()`](https://matplotlib.org/2.1.1/api/_as_gen/matplotlib.pyplot.plot.html) interpolate (straight segments) between data point and the data are not represented by markers."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> In the examples below you can choose: [`linestyle`](https://matplotlib.org/stable/gallery/lines_bars_and_markers/linestyles.html), [`markers`](https://matplotlib.org/stable/api/markers_api.html), [`color`](https://matplotlib.org/stable/tutorials/colors/colors.html) and add a `label` for each plot. The [`legend()`](https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.legend.html) function should be called for the label to be added on the plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(t, y, marker='o', linestyle='', color='orange', label='no line')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(t, y, marker='o', linestyle='solid', color='orange', label='with line')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Let's play with line and marker options "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(t, y, marker='p', linestyle='solid', color='red', \\\n",
    "         linewidth=1.5, ms=10.0, mew=1.0, mec='black', mfc='orange')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Exercice 3.3 :  <hr style=\"border:1px solid grey\">\n",
    "Found the purpose of `linewidth`, `ms`, `mew`, `mec`and `mfc` in the example above "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr style=\"border:1px solid grey\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Exercice 3.4 :  <hr style=\"border:1px solid grey\">\n",
    "Using the functions of **Exercice 2.1** in **TD2** plot:\n",
    "$$f(x)=4x^2\\exp(-2x)$$ \n",
    "\n",
    "$$g(x)=\\frac{x^2}{2}\\left(1-\\frac{x}{2}\\right)^2\\exp(-x)$$\n",
    "\n",
    "$$h(x)= \\frac{4x^2}{243}\\left(3 - 2x + \\frac{2}{9}x^2\\right)^2 \\exp\\left(-\\frac{2x}{3}\\right)$$\n",
    "for $x\\in[0,20]$ with $n=100$ grid points, following the rules:\n",
    "1. For both curves, there is no markers\n",
    "2. The curves must be solid-orange, -red and -purple for $f(x)$, $g(x)$ and $h(x)$, respectively.\n",
    "3. The legend must be printed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import exp, linspace\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr style=\"border:1px solid grey\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Below an example to set ranges to $x\\in[0,12]$ and $y\\in[0,0.3]$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pylab as plt\n",
    "\n",
    "plt.plot(x,f,marker='',linestyle='solid', color='orange', linewidth=1.5, label='f(x)')\n",
    "plt.plot(x,g,marker='',linestyle='solid', color='red',    linewidth=1.5, label='g(x)')\n",
    "plt.plot(x,h,marker='',linestyle='solid', color='purple', linewidth=1.5, label='h(x)')\n",
    "\n",
    "plt.legend()\n",
    "\n",
    "plt.xlim(0,12)\n",
    "plt.ylim(0,0.3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> $x$- and $y$-axis labels are added using `xlabel()` and `ylabel()`. \n",
    "\n",
    "> The dollar symbols in the label stand for [$\\LaTeX$](https://en.wikibooks.org/wiki/LaTeX/Mathematics) formatting of equations.\n",
    "\n",
    "> $x$- and $y$-axis ticks can be set explicitely using `xticks()` and `yticks()` along with `arange()`\n",
    "\n",
    "> > The figure can be printed in save using the [`savefig()`](https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.savefig.html) function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pylab as plt\n",
    "from numpy import arange\n",
    "\n",
    "plt.plot(x,f,marker='',linestyle='solid', color='orange', linewidth=1.5, \n",
    "         label='$f(x)=4x^2\\exp(-2x)$')\n",
    "plt.plot(x,g,marker='',linestyle='solid', color='red',    linewidth=1.5,\n",
    "         label='$g(x)=\\\\frac{x^2}{2}\\\\left(1-\\\\frac{x}{2}\\\\right)^2\\\\exp(-x)$')\n",
    "plt.plot(x,h,marker='',linestyle='solid', color='purple', linewidth=1.5,\n",
    "         label='$h(x)=\\\\frac{4x^2}{243}\\\\left(3 - 2x + \\\\frac{2}{9}x^2\\\\right)^2 \\\\exp\\\\left(-\\\\frac{2x}{3}\\\\right)$')\n",
    "plt.legend()\n",
    "\n",
    "plt.xlabel('$x$')\n",
    "plt.ylabel('$function$')\n",
    "\n",
    "# x-ticks for x in [0,20] with step 1\n",
    "plt.xticks(arange(0.0, 20, 1 ))\n",
    "\n",
    "# y-ticks for y in [0,0.3] with step 0.02\n",
    "plt.yticks(arange(0.0, 0.3, 0.02))\n",
    "\n",
    "plt.xlim(0,15)\n",
    "plt.ylim(0,0.3)\n",
    "\n",
    "plt.savefig('figure_examples/figure_test.png', dpi=300, format='png', transparent=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Exercice 3.5 :   <hr style=\"border:1px solid grey\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We provide a data file `comete_C2022E3ZTF.txt` which gather observations of the C/2022 E3 (ZTF) comet from the 22th of March 2022 to the 25th of February 2023. These data are available from the [LESIA website](https://lesia.obspm.fr/comets/lib/all-obs-table.php?Code=CK22E030&y1=1908&m1=01). We want to plot the distance with respect to the sun ($r$) and to the earth ($d$) as a function of the [Julian day](https://fr.wikipedia.org/wiki/Jour_julien) ($j$). A visualisation of the trajectory is available [here](https://ssd.jpl.nasa.gov/tools/sbdb_lookup.html#/?sstr=C%2F2022%20E3&view=VOP).\n",
    "\n",
    "1. Go online and check which where are located $j$, $r$ and $d$ in the table.\n",
    "2. Check the data file located in the directory `data_examples`and identify the corresponding rows. \n",
    "3. Load the data and extract the 3 columns.\n",
    "4. On the same graph plot $r$ and $d$ as a function of $j$.\n",
    "5. Add axis labels and a legend.\n",
    "\n",
    "***Warning:*** *you may have trouble with some data which can't be converted to float. Find a way to circumvent\n",
    "the problem by using, for instance, the function `replace()` to replace the character by a float, say 0.0.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pylab as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr style=\"border:1px solid grey\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Exercice 3.6 :   <hr style=\"border:1px solid grey\">\n",
    "From the data file `harelynx.csv` located in `data_examples`, plot the variation of population of hares and lynx as a function of time.\n",
    "Follow the guidelines:\n",
    "1. the data from the file `lynxlievre.csv` should be loaded using pandas.\n",
    "2. the data shall be converted to a `numpy array`.\n",
    "3. from the `numpy array` 3 columns will be extracted $(t,x,y)$.\n",
    "4. on the same graphic plot the $x(t)$ for the hare and $y(t)$ for the lynx.\n",
    "5. compute the mean of the 2 populations using the function [`mean()`](https://numpy.org/doc/stable/reference/generated/numpy.mean.html) from numpy library.\n",
    "6. plot the mean as horizontal lines using [`axhline()`](https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.axhline.html) function from `matplotlib.pyplot`.\n",
    "6. add the legend and a title.\n",
    "7. adapts the ticks spacing along $x$ and $y$, the size of the figure, the color and type of markers.\n",
    "\n",
    "**An template is given below:**\n",
    "<img src=\"./figure_examples/figure_lynxhare_ref.png\" alt=\"Alternative text\" width=\"500\"/>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pandas import read_csv\n",
    "from numpy  import arange, mean\n",
    "import matplotlib.pylab as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr style=\"border:1px solid grey\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <hr style=\"border:1px solid black\">  Survival guide for practical work <hr style=\"border:1px solid black\">\n",
    "</br>\n",
    "\n",
    "Given an example, in this survival guide we will:\n",
    "\n",
    "1. Plot the data using python. Add axis labels, units, legend.\n",
    "2. Fit a linear model to the data. \n",
    "3. Superimpose this line on the plot and save the figure.\n",
    "4. Use the linear model equation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In water-resources engineering, the sizing of reservoirs depends on accurate estimates of\n",
    "water flow in the river that is being impounded. For some rivers, long-term historical\n",
    "records of such flow data are difficult to obtain. In contrast, meteorological data on\n",
    "precipitation are often available for many years past. Therefore, it is often useful to\n",
    "determine a relationship between flow and precipitation. This relationship can then be\n",
    "used to estimate flows for years when only precipitation measurements were made. The\n",
    "following data are available for a river: [<sup>1</sup>](#fn1)\n",
    "\n",
    "*Dans le domaine de l'ingénierie des ressources en eau, le dimensionnement des réservoirs dépend d'une estimation précise du débit de la rivière que l'on considère. Pour certains cours d'eau, il est difficile d'obtenir des enregistrements historiques à long terme du débit. En revanche, les données météorologiques sur les précipitations sont souvent disponibles pour de plus longues périodes. Par conséquent, il est souvent utile de déterminer une relation entre le débit et les précipitations. Cette relation peut alors être utilisée pour estimer les débits pour les années où seules les précipitations ont été mesurées. Les données suivantes sont disponibles pour une rivière :*\n",
    "\n",
    "Precip. (cm/yr)|  88.9 | 108.5 | 104.1 | 139.7 | 127.0 | 94.0 | 116.8 | 99.1 |\n",
    "---------------|-------|-------|-------|-------|-------|------|-------|------|\n",
    "Flow (m$^3$/s) |  14.6 |  16.7 |  15.3 |  23.2 |  19.5 | 16.1 |  18.1 | 16.6 |\n",
    "\n",
    "\n",
    "\n",
    "<span id=\"fn1\"> $^1$ Example taken from *Applied Numerical Methods with MATLAB* by S. C. Chapra.</span>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <hr style=\"border:1px solid black\"> Load the data <hr style=\"border:1px solid black\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Create numpy array with the raw data\n",
    "\n",
    "> Here we first create list using the data which will be converted to `numpy arrays`\n",
    "\n",
    "> Warning : we used English floating-point notation !"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import array\n",
    "\n",
    "prec = [88.9, 108.5, 104.1, 139.7, 127.0, 94.0, 116.8, 99.1]\n",
    "flow = [14.6,  16.7,  15.3,  23.2,  19.5, 16.1,  18.1, 16.6]\n",
    "\n",
    "x = array(prec)\n",
    "y = array(flow)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### ... or load the data from a file using python intrinsic function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "file = open('data_examples/prec_flow.csv','r')\n",
    "file.readline()\n",
    "\n",
    "x = [ ] ; y = [ ]\n",
    "for line in file:\n",
    "    words = line.split(',')\n",
    "    x.append( float(words[0]) )\n",
    "    y.append( float(words[1]) )\n",
    "print(prec)\n",
    "print()\n",
    "print(flow)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### ... or load the data from a file using pandas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pandas import read_csv\n",
    "\n",
    "data_pandas = read_csv('data_examples/prec_flow.csv')\n",
    "print(data_pandas)\n",
    "#data_numpy_array = data_pandas.to_numpy()\n",
    "#x = data_numpy_array[:,0]\n",
    "#y = data_numpy_array[:,1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <hr style=\"border:1px solid black\"> Plot the graphics <hr style=\"border:1px solid black\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Using Matplotlib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pylab as plt\n",
    "\n",
    "plt.plot(x,y,marker='o',color='black', linestyle='')\n",
    "plt.xlabel(\"Precipitation (cm/yr)\")\n",
    "plt.ylabel(\"Flow (m$^3$/s)\")\n",
    "\n",
    "plt.xticks(arange(min(x)+1.1,max(x)+1,5))\n",
    "plt.yticks(arange(min(y)+0.4,max(y),1))\n",
    "plt.title('Water flow as a function of precipitation')\n",
    "\n",
    "#plt.xlim()\n",
    "#plt.ylim()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Using Pandas"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Pandas can be used to performed plot with the function [`DataFrame.plot`](https://pandas.pydata.org/docs/reference/api/pandas.DataFrame.plot.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_pandas.plot(x='prec',y='flow',kind='scatter',color='black')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <hr style=\"border:1px solid black\"> Perform a linear fit <hr style=\"border:1px solid black\">\n",
    "\n",
    "> The fit is performed using `curve_fit` from the library `scipy.optimize`\n",
    "\n",
    "> Given a the data points $\\{x_i,y_i\\}^N_{i=1}$, the standard error $\\{\\sigma_i\\}^N_{i=1}$ on $\\{y_i\\}^N_{i=1}$ and the model function $f_\\textrm{model}(x;a_0,a_1,...a_m)$ where $\\{a_k\\}_{k=1}^m$ are the coefficient to optimized. Given $\\chi^2$:\n",
    "$$\\chi^2=\\sum_{i=1}^{N}\\frac{\\left[y_i-f_\\textrm{model}\\left(x_i;a_0,a_1,...a_m\\right)\\right]^2}{\\sigma^2_i}$$\n",
    "`curve_fit`perform the $\\chi^2$ minimization with respect to these coefficients."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we define the model function ; in this case a linear function: $f(x)=a_0 +a_1x$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fmodel(x,a0,a1):\n",
    "    return a0 + a1*x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we call the function `curve_fit` \n",
    "\n",
    "> with inputs $(f_\\textrm{model},x,y,\\sigma)$. \n",
    "\n",
    "> In output we have the optimal values for the coefficients (`popt`) and the (co)variance (`pcov`) of the coefficient estimates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import curve_fit\n",
    "from numpy import ones, sqrt\n",
    "\n",
    "# standard error is set to 1 in this example\n",
    "sig = ones((len(x)))\n",
    "\n",
    "# call curve_fit to obtain the optimized coefficients and covariance\n",
    "popt, pcov = curve_fit(f=fmodel, xdata=x, ydata=y, sigma=sig)\n",
    "\n",
    "# extract optimized coefficients\n",
    "[a0, a1] = popt\n",
    "\n",
    "# extract standard deviation errors\n",
    "sa0 = sqrt(pcov[0,0])\n",
    "sa1 = sqrt(pcov[1,1])\n",
    "\n",
    "print('a0 = %.4f +/- %.4f'%(a0,sa0))\n",
    "print('a1 = %.4f +/- %.4f'%(a1,sa1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <hr style=\"border:1px solid black\"> Plot the linear model <hr style=\"border:1px solid black\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given $x_\\textrm{model}\\in[\\min(x),\\max(x)]$ with $N=128$ compute $y_\\textrm{model}=f_\\textrm{model}(x_\\textrm{model};a_0,a_1)$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import linspace\n",
    "\n",
    "x_model = linspace(min(x),max(x),128)\n",
    "y_model = fmodel(x=x_model,a0=a0,a1=a1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the experimental data and the linear model function ; save the figure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pylab as plt\n",
    "\n",
    "plt.plot(x,      y,       marker='o', color='black', linestyle='', label='exp. data')\n",
    "plt.plot(x_model,y_model, marker= '', color='blue',  linestyle='-',\n",
    "         label='linear regression with $f(x)= %.3f + %.3f x$'%(a0,a1))\n",
    "\n",
    "plt.xlabel(\"Precipitation (cm/yr)\")\n",
    "plt.ylabel(\"Flow (m$^3$/s)\")\n",
    "\n",
    "plt.xticks(arange(min(x)+1.1,max(x)+1,5))\n",
    "plt.yticks(arange(min(y)+0.4,max(y),1))\n",
    "\n",
    "plt.legend()\n",
    "plt.title('Water flow as a function of precipitation')\n",
    "\n",
    "plt.savefig('figure_examples/figure_prec_flow.png', dpi=300, format='png', transparent=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <hr style=\"border:1px solid black\"> Use the linear model <hr style=\"border:1px solid black\">\n",
    "\n",
    "Evaluate the annual water flow if the precipitation is 120 cm/year."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prec = 120\n",
    "flow = fmodel(x=prec,a0=a0,a1=a1)\n",
    "print('flow = %.4f m^3/s for prec. = %.4f'%(prec,flow))"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Diaporama",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}