04_OSEMN.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Data science is OSEMN\n",
    "\n",
    "According to a popular model, the elements of data science are\n",
    "\n",
    "* Obtaining data\n",
    "* Scrubbing data\n",
    "* Exploring data\n",
    "* Modeling data\n",
    "* iNterpreting data\n",
    "\n",
    "and hence the acronym OSEMN, pronounced as “Awesome”.\n",
    "\n",
    "We will start with the **O**, moving towards the rest later, but first let's have a quick look at what it all boils down to"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "data = np.loadtxt('populations.txt')\n",
    "year, hares, lynxes, carrots = data.T # trick: columns to variables\n",
    "\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "plt.figure(figsize=(12,6))\n",
    "plt.axes([0.2, 0.1, 0.5, 0.8]) \n",
    "plt.plot(year, hares, year, lynxes, year, carrots) \n",
    "plt.legend(('Hare', 'Lynx', 'Carrot'), loc=(1.05, 0.5)) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import Image\n",
    "Image(\"cases_vs_searches.jpeg\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By plotting the data a clear (and reasonable) correlations between pray and predator becomes evident. How can it be quantified? Is that statistical significant? What about the correlation between carrots and hares? Is that evident? Is that significant?\n",
    "\n",
    "Finding correlations in data is the main goal of data science, though that is not the end of the story: as this precious [site](http://tylervigen.com/spurious-correlations) demonstrates, **correlations is not causation**. \n",
    "\n",
    "(I've been invited to a school of Philosophy of Science to talk about the role of ML in Physics and they even asked me to write a summary of that. You find it [here](https://www.dropbox.com/s/lcjpgsrsoi0iyq3/Fisica_ML.docx), in Italian.. (humanists like that better than English))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Exercise*: write an algorithm that determins and quantifies a correlation between two time series. Use as an example the hare-lynx-carrot dataset."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**N.B:** \n",
    "If in the cells below you import a packpage not yet installed, you can either install it the usual way, or run a cell like the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# uncomment and set  NAME_OF_THE_PACKAGE to what you need\n",
    "'''\n",
    "import subprocess\n",
    "subprocess.call(['pip', 'install', 'NAME_OF_THE_PACKAGE'])\n",
    "'''"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Obtaining and processing (remote) data\n",
    "\n",
    "Accessing data is a really serious business. Data can sit on public or on remote machines. In the case of the former, things may be straightforward, whereas in the latter case you need to worry about a few things.\n",
    "\n",
    "In both cases, depending on the size of the dataset, the managment of the dataset can become extremely complicated. We won't deal here with large datasets, which require a whole course per se.., but still care should be put. In particular, it is not wise to keep (and even worse commit) data into a git repository!\n",
    "\n",
    "The suggestion is then to create a directory somewhere and copy the example datasets there. From a terminal:\n",
    "\n",
    "```bash\n",
    "\n",
    "# create a data directory in your home directory\n",
    "mkdir ~/data/\n",
    "\n",
    "# check the content (it's empty now of course)\n",
    "ls -ltr ~/data/\n",
    "\n",
    "# in the case you need to move there:\n",
    "cd ~/data/\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Download data from a server\n",
    "\n",
    "A nice set of interesting datasets can be found on this [server](https://archive.ics.uci.edu/ml/datasets.php) that collects training/test data for machine learning developments. Several of those pertein physical sciences, it is worth browsing through those.\n",
    "\n",
    "You can download any of those, in the following we will consider a dataset from the MAGIC experiment. For that we will the `wget` command"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get the dataset and its description on the proper data directory\n",
    "!wget https://archive.ics.uci.edu/ml/machine-learning-databases/magic/magic04.data -P ~/data/\n",
    "!wget https://archive.ics.uci.edu/ml/machine-learning-databases/magic/magic04.names -P ~/data/    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# print the description. This can (and better) be done from a terminal\n",
    "!cat ~/data/magic04.names"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is possible to download and load remote files via their url's directly from within python (and thus on a jupyter session). This is a rather powerful tool as it allows http communications, IO streaming and so on.\n",
    "\n",
    "Care should be put as the dataset is stored in memory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import urllib.request\n",
    "url ='https://archive.ics.uci.edu/ml/machine-learning-databases/magic/magic04.names'\n",
    "with urllib.request.urlopen(url) as data_file:\n",
    "    #print (data_file.read(300))\n",
    "    for line in data_file:\n",
    "        print (line)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Copy data from a remote machine\n",
    "\n",
    "Often datasets are not available on websites but rather they are sitting on some remote machine. Several tools are there that can allow you to get hold off remote data, even from within python (e.g. [paramiko](https://www.paramiko.org/)), but best in this case is to get a local copy. E.g. from a terminal:\n",
    "\n",
    "```bash\n",
    "scp guest104@gatep.fisica.unipd.it:~/data/data_000637.* ~/data/\n",
    "```\n",
    "\n",
    "by issuing that command you are immediately exposed to the most relevant problem in obtaining the data: permissions/authorization.\n",
    "\n",
    "Indeed that will not work (as you don't have an account on that machine and I'd be put into jail if I gave you the password), still you'll need that file later, so \"wget\" it:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!wget https://www.dropbox.com/s/69xe1d5f19nvtw3/data_000637.dat -P ~/data/\n",
    "\n",
    "# copy the interpreted version as well\n",
    "!wget https://www.dropbox.com/s/xvjzaxzz3ysphme/data_000637.txt -P ~/data/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Secondily (essentially a further consequence of the same issue), the remote machine itself may have accessibility restrictions, e.g. being behind a firewall. In that case you may need to use a tunnel:\n",
    "\n",
    "``` bash\n",
    "ssh -L 1234:<address of R known to G>:22 <user at G>@<address of G> \n",
    "\n",
    "scp -P 1234 <user at R>@127.0.0.1:/path/to/file file-name-to-be-copied\n",
    "```\n",
    "\n",
    "In summary, just getting the data is a complicated business."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data Formats\n",
    "\n",
    "datasets can be stored in a gazillion different ways, often they have formats which are application dependent, even though more and more standards are being established. Python have \"readers\" for most of the formats, another reason for being the optimal programming language for data analysis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Text files \n",
    "\n",
    "Plain text files are commonly used for \"readibility\", at the price of a very poor storing efficiency due to their low entropy. [UTF-8](https://en.wikipedia.org/wiki/UTF-8) is the most common encoding.\n",
    "\n",
    "Reading (and writing) text files in python is straightforward:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "file_name = \"/Users/mzanetti/data/magic04.data\"\n",
    "\n",
    "# mode can be specified for writing, reading or both\n",
    "with open(file_name, mode='r') as f:\n",
    "    # print-out the whole file\n",
    "    # print (f.read()) \n",
    "    for line in f:\n",
    "        ## print line by line\n",
    "        print (line)\n",
    "        ## each line is a string, you need to split it yourself\n",
    "        # for c in line.split(): print(c) # check the functionalities of the split() method "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### CSV files\n",
    "\n",
    "If you are lucky text files are already framed into a defined structured, in a \"table-like\" manner. These files are colled \"comma separated values\" (csv), even though the separator may well not be the \",\" symbol.\n",
    "Python have package to deal with that:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import csv\n",
    "\n",
    "with open('/Users/mzanetti/data/magic04.data') as data_file:\n",
    "    for line in csv.reader(data_file, delimiter=','): # the delimiter is often guessed by the reader\n",
    "        # again note that elements of each line are treated as strings\n",
    "        # if you need to convert them into numbers, you need to to that yourself\n",
    "        fLength,fWidth,fSize,\\\n",
    "        fConc,fConc1,fAsym,\\\n",
    "        fM3Long,fM3Trans,fAlpha,fDist = map(float,line[:-1])\n",
    "        category = line[-1]\n",
    "        print (fLength,fWidth,fSize,fConc,fConc1,fAsym,fM3Long,fM3Trans,fAlpha,fDist)\n",
    "        print (category)\n",
    "        break"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "More often than not, csv files have comments (e.g. starting with '#'), which cannot be interpreted by the reader. Tricks like:\n",
    "\n",
    "```python\n",
    "csv.reader(row for row in f if not row.startswith('#'))\n",
    "```\n",
    "\n",
    "may be useful"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Binary (hexadecimal) files\n",
    "\n",
    "The output of sensors often is stored as hexadecimal files. Information is packed in a well defined format (similarly to how floating point numbers are formatted).\n",
    "To read and process hexadecimal files in python you need to use the \"b\" option of `open` and progress along the file at step of defined lenght (depending on the size of the words information is packed into)\n",
    "\n",
    "There are several tool to display and edit hex/bin files, e.g. this [one](https://hexed.it/)\n",
    "\n",
    "The following is an example from data collected from an FPGA implementing a TDC. Relevant infomation are the coordinates of the TDC channels and their time measurements."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import struct, time\n",
    "\n",
    "with  open('/Users/mzanetti/data/data_000637.dat','rb') as file:\n",
    "    file_content=file.read()\n",
    "    word_counter=0\n",
    "    word_size = 8 # size of the word in bytes\n",
    "    for i in range(0, len(file_content), word_size):\n",
    "        word_counter+=1\n",
    "        if word_counter>100: break\n",
    "        time.sleep(0.1)\n",
    "        thisInt = struct.unpack('<q', file_content[i:i+word_size])[0]\n",
    "        head = (thisInt >> 62) & 0x3\n",
    "        if head == 1:\n",
    "            fpga     = (thisInt >> 58) & 0xF\n",
    "            tdc_chan = (thisInt >> 49) & 0x1FF\n",
    "            orb_cnt  = (thisInt >> 17) & 0xFFFFFFFF\n",
    "            bx       = (thisInt >> 5 ) & 0xFFF\n",
    "            tdc_meas = (thisInt >> 0 ) & 0x1F\n",
    "            if i==0 : print ('{0},{1},{2},{3},{4},{5}'.format('HEAD', 'FPGA', 'TDC_CHANNEL', 'ORB_CNT', 'BX', 'TDC_MEAS'))\n",
    "            print ('{0},{1},{2},{3},{4},{5}'.format(head, fpga, tdc_chan, orb_cnt, bx, tdc_meas))\n",
    "        else:\n",
    "            print ('ERROR! head =', head)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### JSON files\n",
    "\n",
    "JSON is JavaScript Object Notation - a format used widely for web-based resource sharing. It is very similar in structure to a Python nested dictionary. Here is an example from http://json.org/example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%file example.json\n",
    "{\n",
    "    \"glossary\": {\n",
    "        \"title\": \"example glossary\",\n",
    "            \"GlossDiv\": {\n",
    "            \"title\": \"S\",\n",
    "                    \"GlossList\": {\n",
    "                \"GlossEntry\": {\n",
    "                    \"ID\": \"SGML\",\n",
    "                                    \"SortAs\": \"SGML\",\n",
    "                                    \"GlossTerm\": \"Standard Generalized Markup Language\",\n",
    "                                    \"Acronym\": \"SGML\",\n",
    "                                    \"Abbrev\": \"ISO 8879:1986\",\n",
    "                                    \"GlossDef\": {\n",
    "                        \"para\": \"A meta-markup language, used to create markup languages such as DocBook.\",\n",
    "                                            \"GlossSeeAlso\": [\"GML\", \"XML\"]\n",
    "                    },\n",
    "                                    \"GlossSee\": \"markup\"\n",
    "                }\n",
    "            }\n",
    "        }\n",
    "    }\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!cat example.json"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import json\n",
    "data = json.load(open('example.json'))\n",
    "print (data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# and can be parsed using standard key lookups\n",
    "data['glossary']['GlossDiv']['GlossList']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### HDF5\n",
    "\n",
    "HDF5 is a hierarchical format often used to store complex scientific data. For instance, Matlab now saves its data to HDF5. It is particularly useful to store complex hierarchical data sets with associated metadata, for example, the results of a computer simulation experiment.\n",
    "\n",
    "The main concepts associated with HDF5 are\n",
    "\n",
    "* file: container for hierachical data - serves as ‘root’ for tree\n",
    "* group: a node for a tree\n",
    "* dataset: array for numeric data - can be huge\n",
    "* attribute: small pieces of metadata that provide additional context\n",
    "\n",
    "Now let's create a dummy dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import h5py \n",
    "\n",
    "\n",
    "#Now mock up some simple dummy data to save to our file.\n",
    "d1 = np.random.random(size = (1000,20))\n",
    "d2 = np.random.random(size = (1000,200))\n",
    "\n",
    "print (d1.shape, d2.shape)\n",
    "\n",
    "\n",
    "hf = h5py.File('data.h5', 'w')\n",
    "hf.create_dataset('dataset_1', data=d1)\n",
    "hf.create_dataset('dataset_2', data=d2)\n",
    "hf.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and read it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hf = h5py.File('data.h5', 'r')\n",
    "\n",
    "print (hf.keys())\n",
    "\n",
    "n1 = hf.get('dataset_1')\n",
    "print (n1)\n",
    "\n",
    "n1 = np.array(n1)\n",
    "print (n1.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "here is another example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import h5py\n",
    "import numpy as np\n",
    "import os\n",
    "\n",
    "# creating a HDF5 file\n",
    "import datetime\n",
    "\n",
    "if not os.path.exists('example.hdf5'):\n",
    "\n",
    "    with h5py.File('example.hdf5') as f:\n",
    "        project = f.create_group('project')\n",
    "        expt1 = project.create_group('expt1')\n",
    "        expt2 = project.create_group('expt2')\n",
    "        expt1.create_dataset('counts', (100,), dtype='i')\n",
    "        expt2.create_dataset('values', (1000,), dtype='f')\n",
    "\n",
    "        expt1['counts'][:] = range(100)\n",
    "        expt2['values'][:] = np.random.random(1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with h5py.File('example.hdf5') as f:\n",
    "    project = f['project']\n",
    "    print (project['expt1']['counts'][:10])\n",
    "    print (project['expt2']['values'][:10])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pandas\n",
    "\n",
    "the most convenient tool to read and process formatted dataset is however Pandas. In the following a couple of examples. Pandas will be the main subject of the next classes\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "file_name=\"/Users/mzanetti/data/data_000637.txt\"\n",
    "data=pd.read_csv(file_name,nrows=10,skiprows=range(1,1))\n",
    "data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!cat /Users/mzanetti/data/data_000637.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "file_name=\"/Users/mzanetti/data/magic04.data\"\n",
    "data=pd.read_csv(file_name,nrows=1000)\n",
    "data.columns=['fLength','fWidth','fSize',\n",
    "        'fConc','fConc1','fAsym',\n",
    "        'fM3Long','fM3Trans','fAlpha','fDist','category']\n",
    "data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "data.plot.scatter(\"fLength\",\"fWidth\",)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.hist(\"fAlpha\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ROOT\n",
    "\n",
    "[ROOT](https://root.cern.ch/) needs a special mention. It is still nowadays, and by far, the most convenient tool to store and manage complex datasets pertaining physics experiments where \"events\" are recorded, in particular High Energy, Nuclear, Astro physics. It allows a nested structure, with complex data objects (classes) and references between them.\n",
    "\n",
    "ROOT per se is obnoxious as it has been developed in the years as way-too-many-purposes package, but it's I/O is formidable.\n",
    "\n",
    "Installing ROOT is a pain. \n",
    "\n",
    "(bare) ROOT files can be opened with non-ROOT library, [uproot](https://uproot.readthedocs.io/en/latest/index.html) (check its git [repo](https://github.com/scikit-hep/uproot))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import uproot\n",
    "\n",
    "\n",
    "events = uproot.open(\"https://scikit-hep.org/uproot3/examples/Zmumu.root\")[\"events\"]\n",
    "events"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "events.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "array = events[\"E1\"].array(library=\"np\")\n",
    "array\n",
    "#plt.hist(array)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For more complicated examples, you can take a look at the CMS experiment [open data](http://opendata.cern.ch/docs/about-cms); sure you can find the Higgs boson in there, extracting the signal is way easier than obtaining, storing, and interpreting the data.."
   ]
  }
 ],
 "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}