session_07_industry_talk_ML.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "78457e09",
   "metadata": {
    "toc": true
   },
   "source": [
    "<h1>Table of Contents<span class=\"tocSkip\"></span></h1>\n",
    "<div class=\"toc\"><ul class=\"toc-item\"><li><span><a href=\"#Introduction\" data-toc-modified-id=\"Introduction-1\"><span class=\"toc-item-num\">1&nbsp;&nbsp;</span>Introduction</a></span><ul class=\"toc-item\"><li><span><a href=\"#Load-packages\" data-toc-modified-id=\"Load-packages-1.1\"><span class=\"toc-item-num\">1.1&nbsp;&nbsp;</span>Load packages</a></span></li><li><span><a href=\"#Load-dataset:-Groundwater-chemistry-South-Australia\" data-toc-modified-id=\"Load-dataset:-Groundwater-chemistry-South-Australia-1.2\"><span class=\"toc-item-num\">1.2&nbsp;&nbsp;</span>Load dataset: Groundwater chemistry South Australia</a></span></li></ul></li><li><span><a href=\"#Clustering\" data-toc-modified-id=\"Clustering-2\"><span class=\"toc-item-num\">2&nbsp;&nbsp;</span>Clustering</a></span><ul class=\"toc-item\"><li><span><a href=\"#Preprocess-the-data\" data-toc-modified-id=\"Preprocess-the-data-2.1\"><span class=\"toc-item-num\">2.1&nbsp;&nbsp;</span>Preprocess the data</a></span></li><li><span><a href=\"#Clustering-with-DBSCAN\" data-toc-modified-id=\"Clustering-with-DBSCAN-2.2\"><span class=\"toc-item-num\">2.2&nbsp;&nbsp;</span>Clustering with DBSCAN</a></span></li></ul></li><li><span><a href=\"#Classification\" data-toc-modified-id=\"Classification-3\"><span class=\"toc-item-num\">3&nbsp;&nbsp;</span>Classification</a></span><ul class=\"toc-item\"><li><span><a href=\"#Prepare-dataset\" data-toc-modified-id=\"Prepare-dataset-3.1\"><span class=\"toc-item-num\">3.1&nbsp;&nbsp;</span>Prepare dataset</a></span></li><li><span><a href=\"#Classify-with-Random-Forests\" data-toc-modified-id=\"Classify-with-Random-Forests-3.2\"><span class=\"toc-item-num\">3.2&nbsp;&nbsp;</span>Classify with Random Forests</a></span></li><li><span><a href=\"#Nitrate-class-prediction-based-on-TDS-and-major-ions\" data-toc-modified-id=\"Nitrate-class-prediction-based-on-TDS-and-major-ions-3.3\"><span class=\"toc-item-num\">3.3&nbsp;&nbsp;</span>Nitrate class prediction based on TDS and major ions</a></span></li></ul></li><li><span><a href=\"#Regression\" data-toc-modified-id=\"Regression-4\"><span class=\"toc-item-num\">4&nbsp;&nbsp;</span>Regression</a></span></li><li><span><a href=\"#Dimensionality-reduction\" data-toc-modified-id=\"Dimensionality-reduction-5\"><span class=\"toc-item-num\">5&nbsp;&nbsp;</span>Dimensionality reduction</a></span></li></ul></div>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "8ea317cd",
   "metadata": {},
   "source": [
    "# Introduction\n",
    "## Load packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "2e84aa0b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# preamble\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "# import functions we'll need later\n",
    "# allow interactive figures in notebook\n",
    "# %matplotlib notebook"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "ce4ad473",
   "metadata": {},
   "source": [
    "## Load dataset: Groundwater chemistry South Australia\n",
    "Pandas can also read data directly from excel files, using [pd.read_excel](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.read_excel.html).\n",
    "The dataset we are using is a dataset of groundwater chemistry of South Australia\n",
    "\n",
    "[Gray, David J. and Bardwell, Nicole (2016) Hydrogeochemistry of South Australia: Data Release: Accompanying Notes. CSIRO, Australia. EP156116 34p](https://data.csiro.au/collections/collection/CIcsiro:17862v1)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "ff778cc8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Drillhole</th>\n",
       "      <th>Date</th>\n",
       "      <th>Unit_No</th>\n",
       "      <th>Sample_No</th>\n",
       "      <th>CSIRO_No</th>\n",
       "      <th>Long</th>\n",
       "      <th>Lat</th>\n",
       "      <th>TDSc_mgL</th>\n",
       "      <th>pH</th>\n",
       "      <th>Eh_mV</th>\n",
       "      <th>...</th>\n",
       "      <th>18O_SMOW</th>\n",
       "      <th>2H_SMOW</th>\n",
       "      <th>Ra_BqL</th>\n",
       "      <th>13C_pmPDB</th>\n",
       "      <th>87Sr_rat</th>\n",
       "      <th>Alpha_BqL</th>\n",
       "      <th>Beta_BqL</th>\n",
       "      <th>14C_pmc</th>\n",
       "      <th>14C_age</th>\n",
       "      <th>222Rn_BqL</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>143</td>\n",
       "      <td>1941-07-18</td>\n",
       "      <td>473400001</td>\n",
       "      <td>97764</td>\n",
       "      <td>C00001</td>\n",
       "      <td>129.110034</td>\n",
       "      <td>-31.584249</td>\n",
       "      <td>159.870772</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>144</td>\n",
       "      <td>1941-07-18</td>\n",
       "      <td>473400002</td>\n",
       "      <td>1831820</td>\n",
       "      <td>C00002</td>\n",
       "      <td>129.346375</td>\n",
       "      <td>-31.521174</td>\n",
       "      <td>186.619224</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>146</td>\n",
       "      <td>1973-07-06</td>\n",
       "      <td>473400004</td>\n",
       "      <td>399452</td>\n",
       "      <td>C00003</td>\n",
       "      <td>129.385055</td>\n",
       "      <td>-31.640240</td>\n",
       "      <td>30817.910000</td>\n",
       "      <td>7.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>147</td>\n",
       "      <td>1973-06-30</td>\n",
       "      <td>473400005</td>\n",
       "      <td>1831822</td>\n",
       "      <td>C00004</td>\n",
       "      <td>129.157832</td>\n",
       "      <td>-31.665696</td>\n",
       "      <td>19035.280000</td>\n",
       "      <td>7.5</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>147</td>\n",
       "      <td>1973-07-21</td>\n",
       "      <td>473400005</td>\n",
       "      <td>399459</td>\n",
       "      <td>C00005</td>\n",
       "      <td>129.157832</td>\n",
       "      <td>-31.665696</td>\n",
       "      <td>18970.360000</td>\n",
       "      <td>7.1</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>32292</th>\n",
       "      <td>276087</td>\n",
       "      <td>2013-02-26</td>\n",
       "      <td>653102002</td>\n",
       "      <td>2026137</td>\n",
       "      <td>C32338</td>\n",
       "      <td>138.008903</td>\n",
       "      <td>-33.159667</td>\n",
       "      <td>11192.837466</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>32293</th>\n",
       "      <td>276089</td>\n",
       "      <td>2013-02-26</td>\n",
       "      <td>653102004</td>\n",
       "      <td>2026138</td>\n",
       "      <td>C32339</td>\n",
       "      <td>138.007531</td>\n",
       "      <td>-33.159979</td>\n",
       "      <td>40829.201102</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>32294</th>\n",
       "      <td>276091</td>\n",
       "      <td>2013-02-26</td>\n",
       "      <td>653102006</td>\n",
       "      <td>2026135</td>\n",
       "      <td>C32340</td>\n",
       "      <td>138.011364</td>\n",
       "      <td>-33.160587</td>\n",
       "      <td>21575.757576</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>32295</th>\n",
       "      <td>276092</td>\n",
       "      <td>2013-02-27</td>\n",
       "      <td>653102007</td>\n",
       "      <td>2026136</td>\n",
       "      <td>C32341</td>\n",
       "      <td>138.010962</td>\n",
       "      <td>-33.161361</td>\n",
       "      <td>26115.702479</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>32296</th>\n",
       "      <td>276093</td>\n",
       "      <td>2012-11-24</td>\n",
       "      <td>653102008</td>\n",
       "      <td>2025794</td>\n",
       "      <td>C32342</td>\n",
       "      <td>138.010177</td>\n",
       "      <td>-33.159359</td>\n",
       "      <td>41685.196752</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>32297 rows × 77 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "       Drillhole       Date    Unit_No  Sample_No CSIRO_No        Long  \\\n",
       "0            143 1941-07-18  473400001      97764   C00001  129.110034   \n",
       "1            144 1941-07-18  473400002    1831820   C00002  129.346375   \n",
       "2            146 1973-07-06  473400004     399452   C00003  129.385055   \n",
       "3            147 1973-06-30  473400005    1831822   C00004  129.157832   \n",
       "4            147 1973-07-21  473400005     399459   C00005  129.157832   \n",
       "...          ...        ...        ...        ...      ...         ...   \n",
       "32292     276087 2013-02-26  653102002    2026137   C32338  138.008903   \n",
       "32293     276089 2013-02-26  653102004    2026138   C32339  138.007531   \n",
       "32294     276091 2013-02-26  653102006    2026135   C32340  138.011364   \n",
       "32295     276092 2013-02-27  653102007    2026136   C32341  138.010962   \n",
       "32296     276093 2012-11-24  653102008    2025794   C32342  138.010177   \n",
       "\n",
       "             Lat      TDSc_mgL   pH  Eh_mV  ...  18O_SMOW  2H_SMOW  Ra_BqL  \\\n",
       "0     -31.584249    159.870772  NaN    NaN  ...       NaN      NaN     NaN   \n",
       "1     -31.521174    186.619224  NaN    NaN  ...       NaN      NaN     NaN   \n",
       "2     -31.640240  30817.910000  7.0    NaN  ...       NaN      NaN     NaN   \n",
       "3     -31.665696  19035.280000  7.5    NaN  ...       NaN      NaN     NaN   \n",
       "4     -31.665696  18970.360000  7.1    NaN  ...       NaN      NaN     NaN   \n",
       "...          ...           ...  ...    ...  ...       ...      ...     ...   \n",
       "32292 -33.159667  11192.837466  NaN    NaN  ...       NaN      NaN     NaN   \n",
       "32293 -33.159979  40829.201102  NaN    NaN  ...       NaN      NaN     NaN   \n",
       "32294 -33.160587  21575.757576  NaN    NaN  ...       NaN      NaN     NaN   \n",
       "32295 -33.161361  26115.702479  NaN    NaN  ...       NaN      NaN     NaN   \n",
       "32296 -33.159359  41685.196752  NaN    NaN  ...       NaN      NaN     NaN   \n",
       "\n",
       "       13C_pmPDB  87Sr_rat  Alpha_BqL  Beta_BqL  14C_pmc  14C_age  222Rn_BqL  \n",
       "0            NaN       NaN        NaN       NaN      NaN      NaN        NaN  \n",
       "1            NaN       NaN        NaN       NaN      NaN      NaN        NaN  \n",
       "2            NaN       NaN        NaN       NaN      NaN      NaN        NaN  \n",
       "3            NaN       NaN        NaN       NaN      NaN      NaN        NaN  \n",
       "4            NaN       NaN        NaN       NaN      NaN      NaN        NaN  \n",
       "...          ...       ...        ...       ...      ...      ...        ...  \n",
       "32292        NaN       NaN        NaN       NaN      NaN      NaN        NaN  \n",
       "32293        NaN       NaN        NaN       NaN      NaN      NaN        NaN  \n",
       "32294        NaN       NaN        NaN       NaN      NaN      NaN        NaN  \n",
       "32295        NaN       NaN        NaN       NaN      NaN      NaN        NaN  \n",
       "32296        NaN       NaN        NaN       NaN      NaN      NaN        NaN  \n",
       "\n",
       "[32297 rows x 77 columns]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fname = 'CSH_SA.xlsx'\n",
    "sheet = 'Data'\n",
    "chem = pd.read_excel(fname,sheet)\n",
    "chem"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "11b76bfb",
   "metadata": {},
   "source": [
    "# Clustering\n",
    "Cluster the data based on Alkalinity and TDS\n",
    "## Preprocess the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "cd0342db",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/javascript": "/* Put everything inside the global mpl namespace */\n/* global mpl */\nwindow.mpl = {};\n\nmpl.get_websocket_type = function () {\n    if (typeof WebSocket !== 'undefined') {\n        return WebSocket;\n    } else if (typeof MozWebSocket !== 'undefined') {\n        return MozWebSocket;\n    } else {\n        alert(\n            'Your browser does not have WebSocket support. ' +\n                'Please try Chrome, Safari or Firefox ≥ 6. ' +\n                'Firefox 4 and 5 are also supported but you ' +\n                'have to enable WebSockets in about:config.'\n        );\n    }\n};\n\nmpl.figure = function (figure_id, websocket, ondownload, parent_element) {\n    this.id = figure_id;\n\n    this.ws = websocket;\n\n    this.supports_binary = this.ws.binaryType !== undefined;\n\n    if (!this.supports_binary) {\n        var warnings = document.getElementById('mpl-warnings');\n        if (warnings) {\n            warnings.style.display = 'block';\n            warnings.textContent =\n                'This browser does not support binary websocket messages. ' +\n                'Performance may be slow.';\n        }\n    }\n\n    this.imageObj = new Image();\n\n    this.context = undefined;\n    this.message = undefined;\n    this.canvas = undefined;\n    this.rubberband_canvas = undefined;\n    this.rubberband_context = undefined;\n    this.format_dropdown = undefined;\n\n    this.image_mode = 'full';\n\n    this.root = document.createElement('div');\n    this.root.setAttribute('style', 'display: inline-block');\n    this._root_extra_style(this.root);\n\n    parent_element.appendChild(this.root);\n\n    this._init_header(this);\n    this._init_canvas(this);\n    this._init_toolbar(this);\n\n    var fig = this;\n\n    this.waiting = false;\n\n    this.ws.onopen = function () {\n        fig.send_message('supports_binary', { value: fig.supports_binary });\n        fig.send_message('send_image_mode', {});\n        if (fig.ratio !== 1) {\n            fig.send_message('set_device_pixel_ratio', {\n                device_pixel_ratio: fig.ratio,\n            });\n        }\n        fig.send_message('refresh', {});\n    };\n\n    this.imageObj.onload = function () {\n        if (fig.image_mode === 'full') {\n            // Full images could contain transparency (where diff images\n            // almost always do), so we need to clear the canvas so that\n            // there is no ghosting.\n            fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n        }\n        fig.context.drawImage(fig.imageObj, 0, 0);\n    };\n\n    this.imageObj.onunload = function () {\n        fig.ws.close();\n    };\n\n    this.ws.onmessage = this._make_on_message_function(this);\n\n    this.ondownload = ondownload;\n};\n\nmpl.figure.prototype._init_header = function () {\n    var titlebar = document.createElement('div');\n    titlebar.classList =\n        'ui-dialog-titlebar ui-widget-header ui-corner-all ui-helper-clearfix';\n    var titletext = document.createElement('div');\n    titletext.classList = 'ui-dialog-title';\n    titletext.setAttribute(\n        'style',\n        'width: 100%; text-align: center; padding: 3px;'\n    );\n    titlebar.appendChild(titletext);\n    this.root.appendChild(titlebar);\n    this.header = titletext;\n};\n\nmpl.figure.prototype._canvas_extra_style = function (_canvas_div) {};\n\nmpl.figure.prototype._root_extra_style = function (_canvas_div) {};\n\nmpl.figure.prototype._init_canvas = function () {\n    var fig = this;\n\n    var canvas_div = (this.canvas_div = document.createElement('div'));\n    canvas_div.setAttribute('tabindex', '0');\n    canvas_div.setAttribute(\n        'style',\n        'border: 1px solid #ddd;' +\n            'box-sizing: content-box;' +\n            'clear: both;' +\n            'min-height: 1px;' +\n            'min-width: 1px;' +\n            'outline: 0;' +\n            'overflow: hidden;' +\n            'position: relative;' +\n            'resize: both;' +\n            'z-index: 2;'\n    );\n\n    function on_keyboard_event_closure(name) {\n        return function (event) {\n            return fig.key_event(event, name);\n        };\n    }\n\n    canvas_div.addEventListener(\n        'keydown',\n        on_keyboard_event_closure('key_press')\n    );\n    canvas_div.addEventListener(\n        'keyup',\n        on_keyboard_event_closure('key_release')\n    );\n\n    this._canvas_extra_style(canvas_div);\n    this.root.appendChild(canvas_div);\n\n    var canvas = (this.canvas = document.createElement('canvas'));\n    canvas.classList.add('mpl-canvas');\n    canvas.setAttribute(\n        'style',\n        'box-sizing: content-box;' +\n            'pointer-events: none;' +\n            'position: relative;' +\n            'z-index: 0;'\n    );\n\n    this.context = canvas.getContext('2d');\n\n    var backingStore =\n        this.context.backingStorePixelRatio ||\n        this.context.webkitBackingStorePixelRatio ||\n        this.context.mozBackingStorePixelRatio ||\n        this.context.msBackingStorePixelRatio ||\n        this.context.oBackingStorePixelRatio ||\n        this.context.backingStorePixelRatio ||\n        1;\n\n    this.ratio = (window.devicePixelRatio || 1) / backingStore;\n\n    var rubberband_canvas = (this.rubberband_canvas = document.createElement(\n        'canvas'\n    ));\n    rubberband_canvas.setAttribute(\n        'style',\n        'box-sizing: content-box;' +\n            'left: 0;' +\n            'pointer-events: none;' +\n            'position: absolute;' +\n            'top: 0;' +\n            'z-index: 1;'\n    );\n\n    // Apply a ponyfill if ResizeObserver is not implemented by browser.\n    if (this.ResizeObserver === undefined) {\n        if (window.ResizeObserver !== undefined) {\n            this.ResizeObserver = window.ResizeObserver;\n        } else {\n            var obs = _JSXTOOLS_RESIZE_OBSERVER({});\n            this.ResizeObserver = obs.ResizeObserver;\n        }\n    }\n\n    this.resizeObserverInstance = new this.ResizeObserver(function (entries) {\n        var nentries = entries.length;\n        for (var i = 0; i < nentries; i++) {\n            var entry = entries[i];\n            var width, height;\n            if (entry.contentBoxSize) {\n                if (entry.contentBoxSize instanceof Array) {\n                    // Chrome 84 implements new version of spec.\n                    width = entry.contentBoxSize[0].inlineSize;\n                    height = entry.contentBoxSize[0].blockSize;\n                } else {\n                    // Firefox implements old version of spec.\n                    width = entry.contentBoxSize.inlineSize;\n                    height = entry.contentBoxSize.blockSize;\n                }\n            } else {\n                // Chrome <84 implements even older version of spec.\n                width = entry.contentRect.width;\n                height = entry.contentRect.height;\n            }\n\n            // Keep the size of the canvas and rubber band canvas in sync with\n            // the canvas container.\n            if (entry.devicePixelContentBoxSize) {\n                // Chrome 84 implements new version of spec.\n                canvas.setAttribute(\n                    'width',\n                    entry.devicePixelContentBoxSize[0].inlineSize\n                );\n                canvas.setAttribute(\n                    'height',\n                    entry.devicePixelContentBoxSize[0].blockSize\n                );\n            } else {\n                canvas.setAttribute('width', width * fig.ratio);\n                canvas.setAttribute('height', height * fig.ratio);\n            }\n            /* This rescales the canvas back to display pixels, so that it\n             * appears correct on HiDPI screens. */\n            canvas.style.width = width + 'px';\n            canvas.style.height = height + 'px';\n\n            rubberband_canvas.setAttribute('width', width);\n            rubberband_canvas.setAttribute('height', height);\n\n            // And update the size in Python. We ignore the initial 0/0 size\n            // that occurs as the element is placed into the DOM, which should\n            // otherwise not happen due to the minimum size styling.\n            if (fig.ws.readyState == 1 && width != 0 && height != 0) {\n                fig.request_resize(width, height);\n            }\n        }\n    });\n    this.resizeObserverInstance.observe(canvas_div);\n\n    function on_mouse_event_closure(name) {\n        /* User Agent sniffing is bad, but WebKit is busted:\n         * https://bugs.webkit.org/show_bug.cgi?id=144526\n         * https://bugs.webkit.org/show_bug.cgi?id=181818\n         * The worst that happens here is that they get an extra browser\n         * selection when dragging, if this check fails to catch them.\n         */\n        var UA = navigator.userAgent;\n        var isWebKit = /AppleWebKit/.test(UA) && !/Chrome/.test(UA);\n        if(isWebKit) {\n            return function (event) {\n                /* This prevents the web browser from automatically changing to\n                 * the text insertion cursor when the button is pressed. We\n                 * want to control all of the cursor setting manually through\n                 * the 'cursor' event from matplotlib */\n                event.preventDefault()\n                return fig.mouse_event(event, name);\n            };\n        } else {\n            return function (event) {\n                return fig.mouse_event(event, name);\n            };\n        }\n    }\n\n    canvas_div.addEventListener(\n        'mousedown',\n        on_mouse_event_closure('button_press')\n    );\n    canvas_div.addEventListener(\n        'mouseup',\n        on_mouse_event_closure('button_release')\n    );\n    canvas_div.addEventListener(\n        'dblclick',\n        on_mouse_event_closure('dblclick')\n    );\n    // Throttle sequential mouse events to 1 every 20ms.\n    canvas_div.addEventListener(\n        'mousemove',\n        on_mouse_event_closure('motion_notify')\n    );\n\n    canvas_div.addEventListener(\n        'mouseenter',\n        on_mouse_event_closure('figure_enter')\n    );\n    canvas_div.addEventListener(\n        'mouseleave',\n        on_mouse_event_closure('figure_leave')\n    );\n\n    canvas_div.addEventListener('wheel', function (event) {\n        if (event.deltaY < 0) {\n            event.step = 1;\n        } else {\n            event.step = -1;\n        }\n        on_mouse_event_closure('scroll')(event);\n    });\n\n    canvas_div.appendChild(canvas);\n    canvas_div.appendChild(rubberband_canvas);\n\n    this.rubberband_context = rubberband_canvas.getContext('2d');\n    this.rubberband_context.strokeStyle = '#000000';\n\n    this._resize_canvas = function (width, height, forward) {\n        if (forward) {\n            canvas_div.style.width = width + 'px';\n            canvas_div.style.height = height + 'px';\n        }\n    };\n\n    // Disable right mouse context menu.\n    canvas_div.addEventListener('contextmenu', function (_e) {\n        event.preventDefault();\n        return false;\n    });\n\n    function set_focus() {\n        canvas.focus();\n        canvas_div.focus();\n    }\n\n    window.setTimeout(set_focus, 100);\n};\n\nmpl.figure.prototype._init_toolbar = function () {\n    var fig = this;\n\n    var toolbar = document.createElement('div');\n    toolbar.classList = 'mpl-toolbar';\n    this.root.appendChild(toolbar);\n\n    function on_click_closure(name) {\n        return function (_event) {\n            return fig.toolbar_button_onclick(name);\n        };\n    }\n\n    function on_mouseover_closure(tooltip) {\n        return function (event) {\n            if (!event.currentTarget.disabled) {\n                return fig.toolbar_button_onmouseover(tooltip);\n            }\n        };\n    }\n\n    fig.buttons = {};\n    var buttonGroup = document.createElement('div');\n    buttonGroup.classList = 'mpl-button-group';\n    for (var toolbar_ind in mpl.toolbar_items) {\n        var name = mpl.toolbar_items[toolbar_ind][0];\n        var tooltip = mpl.toolbar_items[toolbar_ind][1];\n        var image = mpl.toolbar_items[toolbar_ind][2];\n        var method_name = mpl.toolbar_items[toolbar_ind][3];\n\n        if (!name) {\n            /* Instead of a spacer, we start a new button group. */\n            if (buttonGroup.hasChildNodes()) {\n                toolbar.appendChild(buttonGroup);\n            }\n            buttonGroup = document.createElement('div');\n            buttonGroup.classList = 'mpl-button-group';\n            continue;\n        }\n\n        var button = (fig.buttons[name] = document.createElement('button'));\n        button.classList = 'mpl-widget';\n        button.setAttribute('role', 'button');\n        button.setAttribute('aria-disabled', 'false');\n        button.addEventListener('click', on_click_closure(method_name));\n        button.addEventListener('mouseover', on_mouseover_closure(tooltip));\n\n        var icon_img = document.createElement('img');\n        icon_img.src = '_images/' + image + '.png';\n        icon_img.srcset = '_images/' + image + '_large.png 2x';\n        icon_img.alt = tooltip;\n        button.appendChild(icon_img);\n\n        buttonGroup.appendChild(button);\n    }\n\n    if (buttonGroup.hasChildNodes()) {\n        toolbar.appendChild(buttonGroup);\n    }\n\n    var fmt_picker = document.createElement('select');\n    fmt_picker.classList = 'mpl-widget';\n    toolbar.appendChild(fmt_picker);\n    this.format_dropdown = fmt_picker;\n\n    for (var ind in mpl.extensions) {\n        var fmt = mpl.extensions[ind];\n        var option = document.createElement('option');\n        option.selected = fmt === mpl.default_extension;\n        option.innerHTML = fmt;\n        fmt_picker.appendChild(option);\n    }\n\n    var status_bar = document.createElement('span');\n    status_bar.classList = 'mpl-message';\n    toolbar.appendChild(status_bar);\n    this.message = status_bar;\n};\n\nmpl.figure.prototype.request_resize = function (x_pixels, y_pixels) {\n    // Request matplotlib to resize the figure. Matplotlib will then trigger a resize in the client,\n    // which will in turn request a refresh of the image.\n    this.send_message('resize', { width: x_pixels, height: y_pixels });\n};\n\nmpl.figure.prototype.send_message = function (type, properties) {\n    properties['type'] = type;\n    properties['figure_id'] = this.id;\n    this.ws.send(JSON.stringify(properties));\n};\n\nmpl.figure.prototype.send_draw_message = function () {\n    if (!this.waiting) {\n        this.waiting = true;\n        this.ws.send(JSON.stringify({ type: 'draw', figure_id: this.id }));\n    }\n};\n\nmpl.figure.prototype.handle_save = function (fig, _msg) {\n    var format_dropdown = fig.format_dropdown;\n    var format = format_dropdown.options[format_dropdown.selectedIndex].value;\n    fig.ondownload(fig, format);\n};\n\nmpl.figure.prototype.handle_resize = function (fig, msg) {\n    var size = msg['size'];\n    if (size[0] !== fig.canvas.width || size[1] !== fig.canvas.height) {\n        fig._resize_canvas(size[0], size[1], msg['forward']);\n        fig.send_message('refresh', {});\n    }\n};\n\nmpl.figure.prototype.handle_rubberband = function (fig, msg) {\n    var x0 = msg['x0'] / fig.ratio;\n    var y0 = (fig.canvas.height - msg['y0']) / fig.ratio;\n    var x1 = msg['x1'] / fig.ratio;\n    var y1 = (fig.canvas.height - msg['y1']) / fig.ratio;\n    x0 = Math.floor(x0) + 0.5;\n    y0 = Math.floor(y0) + 0.5;\n    x1 = Math.floor(x1) + 0.5;\n    y1 = Math.floor(y1) + 0.5;\n    var min_x = Math.min(x0, x1);\n    var min_y = Math.min(y0, y1);\n    var width = Math.abs(x1 - x0);\n    var height = Math.abs(y1 - y0);\n\n    fig.rubberband_context.clearRect(\n        0,\n        0,\n        fig.canvas.width / fig.ratio,\n        fig.canvas.height / fig.ratio\n    );\n\n    fig.rubberband_context.strokeRect(min_x, min_y, width, height);\n};\n\nmpl.figure.prototype.handle_figure_label = function (fig, msg) {\n    // Updates the figure title.\n    fig.header.textContent = msg['label'];\n};\n\nmpl.figure.prototype.handle_cursor = function (fig, msg) {\n    fig.canvas_div.style.cursor = msg['cursor'];\n};\n\nmpl.figure.prototype.handle_message = function (fig, msg) {\n    fig.message.textContent = msg['message'];\n};\n\nmpl.figure.prototype.handle_draw = function (fig, _msg) {\n    // Request the server to send over a new figure.\n    fig.send_draw_message();\n};\n\nmpl.figure.prototype.handle_image_mode = function (fig, msg) {\n    fig.image_mode = msg['mode'];\n};\n\nmpl.figure.prototype.handle_history_buttons = function (fig, msg) {\n    for (var key in msg) {\n        if (!(key in fig.buttons)) {\n            continue;\n        }\n        fig.buttons[key].disabled = !msg[key];\n        fig.buttons[key].setAttribute('aria-disabled', !msg[key]);\n    }\n};\n\nmpl.figure.prototype.handle_navigate_mode = function (fig, msg) {\n    if (msg['mode'] === 'PAN') {\n        fig.buttons['Pan'].classList.add('active');\n        fig.buttons['Zoom'].classList.remove('active');\n    } else if (msg['mode'] === 'ZOOM') {\n        fig.buttons['Pan'].classList.remove('active');\n        fig.buttons['Zoom'].classList.add('active');\n    } else {\n        fig.buttons['Pan'].classList.remove('active');\n        fig.buttons['Zoom'].classList.remove('active');\n    }\n};\n\nmpl.figure.prototype.updated_canvas_event = function () {\n    // Called whenever the canvas gets updated.\n    this.send_message('ack', {});\n};\n\n// A function to construct a web socket function for onmessage handling.\n// Called in the figure constructor.\nmpl.figure.prototype._make_on_message_function = function (fig) {\n    return function socket_on_message(evt) {\n        if (evt.data instanceof Blob) {\n            var img = evt.data;\n            if (img.type !== 'image/png') {\n                /* FIXME: We get \"Resource interpreted as Image but\n                 * transferred with MIME type text/plain:\" errors on\n                 * Chrome.  But how to set the MIME type?  It doesn't seem\n                 * to be part of the websocket stream */\n                img.type = 'image/png';\n            }\n\n            /* Free the memory for the previous frames */\n            if (fig.imageObj.src) {\n                (window.URL || window.webkitURL).revokeObjectURL(\n                    fig.imageObj.src\n                );\n            }\n\n            fig.imageObj.src = (window.URL || window.webkitURL).createObjectURL(\n                img\n            );\n            fig.updated_canvas_event();\n            fig.waiting = false;\n            return;\n        } else if (\n            typeof evt.data === 'string' &&\n            evt.data.slice(0, 21) === 'data:image/png;base64'\n        ) {\n            fig.imageObj.src = evt.data;\n            fig.updated_canvas_event();\n            fig.waiting = false;\n            return;\n        }\n\n        var msg = JSON.parse(evt.data);\n        var msg_type = msg['type'];\n\n        // Call the  \"handle_{type}\" callback, which takes\n        // the figure and JSON message as its only arguments.\n        try {\n            var callback = fig['handle_' + msg_type];\n        } catch (e) {\n            console.log(\n                \"No handler for the '\" + msg_type + \"' message type: \",\n                msg\n            );\n            return;\n        }\n\n        if (callback) {\n            try {\n                // console.log(\"Handling '\" + msg_type + \"' message: \", msg);\n                callback(fig, msg);\n            } catch (e) {\n                console.log(\n                    \"Exception inside the 'handler_\" + msg_type + \"' callback:\",\n                    e,\n                    e.stack,\n                    msg\n                );\n            }\n        }\n    };\n};\n\nfunction getModifiers(event) {\n    var mods = [];\n    if (event.ctrlKey) {\n        mods.push('ctrl');\n    }\n    if (event.altKey) {\n        mods.push('alt');\n    }\n    if (event.shiftKey) {\n        mods.push('shift');\n    }\n    if (event.metaKey) {\n        mods.push('meta');\n    }\n    return mods;\n}\n\n/*\n * return a copy of an object with only non-object keys\n * we need this to avoid circular references\n * https://stackoverflow.com/a/24161582/3208463\n */\nfunction simpleKeys(original) {\n    return Object.keys(original).reduce(function (obj, key) {\n        if (typeof original[key] !== 'object') {\n            obj[key] = original[key];\n        }\n        return obj;\n    }, {});\n}\n\nmpl.figure.prototype.mouse_event = function (event, name) {\n    if (name === 'button_press') {\n        this.canvas.focus();\n        this.canvas_div.focus();\n    }\n\n    // from https://stackoverflow.com/q/1114465\n    var boundingRect = this.canvas.getBoundingClientRect();\n    var x = (event.clientX - boundingRect.left) * this.ratio;\n    var y = (event.clientY - boundingRect.top) * this.ratio;\n\n    this.send_message(name, {\n        x: x,\n        y: y,\n        button: event.button,\n        step: event.step,\n        modifiers: getModifiers(event),\n        guiEvent: simpleKeys(event),\n    });\n\n    return false;\n};\n\nmpl.figure.prototype._key_event_extra = function (_event, _name) {\n    // Handle any extra behaviour associated with a key event\n};\n\nmpl.figure.prototype.key_event = function (event, name) {\n    // Prevent repeat events\n    if (name === 'key_press') {\n        if (event.key === this._key) {\n            return;\n        } else {\n            this._key = event.key;\n        }\n    }\n    if (name === 'key_release') {\n        this._key = null;\n    }\n\n    var value = '';\n    if (event.ctrlKey && event.key !== 'Control') {\n        value += 'ctrl+';\n    }\n    else if (event.altKey && event.key !== 'Alt') {\n        value += 'alt+';\n    }\n    else if (event.shiftKey && event.key !== 'Shift') {\n        value += 'shift+';\n    }\n\n    value += 'k' + event.key;\n\n    this._key_event_extra(event, name);\n\n    this.send_message(name, { key: value, guiEvent: simpleKeys(event) });\n    return false;\n};\n\nmpl.figure.prototype.toolbar_button_onclick = function (name) {\n    if (name === 'download') {\n        this.handle_save(this, null);\n    } else {\n        this.send_message('toolbar_button', { name: name });\n    }\n};\n\nmpl.figure.prototype.toolbar_button_onmouseover = function (tooltip) {\n    this.message.textContent = tooltip;\n};\n\n///////////////// REMAINING CONTENT GENERATED BY embed_js.py /////////////////\n// prettier-ignore\nvar _JSXTOOLS_RESIZE_OBSERVER=function(A){var t,i=new WeakMap,n=new WeakMap,a=new WeakMap,r=new WeakMap,o=new Set;function s(e){if(!(this instanceof s))throw new TypeError(\"Constructor requires 'new' operator\");i.set(this,e)}function h(){throw new TypeError(\"Function is not a constructor\")}function c(e,t,i,n){e=0 in arguments?Number(arguments[0]):0,t=1 in arguments?Number(arguments[1]):0,i=2 in arguments?Number(arguments[2]):0,n=3 in arguments?Number(arguments[3]):0,this.right=(this.x=this.left=e)+(this.width=i),this.bottom=(this.y=this.top=t)+(this.height=n),Object.freeze(this)}function d(){t=requestAnimationFrame(d);var s=new WeakMap,p=new Set;o.forEach((function(t){r.get(t).forEach((function(i){var r=t instanceof window.SVGElement,o=a.get(t),d=r?0:parseFloat(o.paddingTop),f=r?0:parseFloat(o.paddingRight),l=r?0:parseFloat(o.paddingBottom),u=r?0:parseFloat(o.paddingLeft),g=r?0:parseFloat(o.borderTopWidth),m=r?0:parseFloat(o.borderRightWidth),w=r?0:parseFloat(o.borderBottomWidth),b=u+f,F=d+l,v=(r?0:parseFloat(o.borderLeftWidth))+m,W=g+w,y=r?0:t.offsetHeight-W-t.clientHeight,E=r?0:t.offsetWidth-v-t.clientWidth,R=b+v,z=F+W,M=r?t.width:parseFloat(o.width)-R-E,O=r?t.height:parseFloat(o.height)-z-y;if(n.has(t)){var k=n.get(t);if(k[0]===M&&k[1]===O)return}n.set(t,[M,O]);var S=Object.create(h.prototype);S.target=t,S.contentRect=new c(u,d,M,O),s.has(i)||(s.set(i,[]),p.add(i)),s.get(i).push(S)}))})),p.forEach((function(e){i.get(e).call(e,s.get(e),e)}))}return s.prototype.observe=function(i){if(i instanceof window.Element){r.has(i)||(r.set(i,new Set),o.add(i),a.set(i,window.getComputedStyle(i)));var n=r.get(i);n.has(this)||n.add(this),cancelAnimationFrame(t),t=requestAnimationFrame(d)}},s.prototype.unobserve=function(i){if(i instanceof window.Element&&r.has(i)){var n=r.get(i);n.has(this)&&(n.delete(this),n.size||(r.delete(i),o.delete(i))),n.size||r.delete(i),o.size||cancelAnimationFrame(t)}},A.DOMRectReadOnly=c,A.ResizeObserver=s,A.ResizeObserverEntry=h,A}; // eslint-disable-line\nmpl.toolbar_items = [[\"Home\", \"Reset original view\", \"fa fa-home\", \"home\"], [\"Back\", \"Back to previous view\", \"fa fa-arrow-left\", \"back\"], [\"Forward\", \"Forward to next view\", \"fa fa-arrow-right\", \"forward\"], [\"\", \"\", \"\", \"\"], [\"Pan\", \"Left button pans, Right button zooms\\nx/y fixes axis, CTRL fixes aspect\", \"fa fa-arrows\", \"pan\"], [\"Zoom\", \"Zoom to rectangle\\nx/y fixes axis\", \"fa fa-square-o\", \"zoom\"], [\"\", \"\", \"\", \"\"], [\"Download\", \"Download plot\", \"fa fa-floppy-o\", \"download\"]];\n\nmpl.extensions = [\"eps\", \"jpeg\", \"pgf\", \"pdf\", \"png\", \"ps\", \"raw\", \"svg\", \"tif\", \"webp\"];\n\nmpl.default_extension = \"png\";/* global mpl */\n\nvar comm_websocket_adapter = function (comm) {\n    // Create a \"websocket\"-like object which calls the given IPython comm\n    // object with the appropriate methods. Currently this is a non binary\n    // socket, so there is still some room for performance tuning.\n    var ws = {};\n\n    ws.binaryType = comm.kernel.ws.binaryType;\n    ws.readyState = comm.kernel.ws.readyState;\n    function updateReadyState(_event) {\n        if (comm.kernel.ws) {\n            ws.readyState = comm.kernel.ws.readyState;\n        } else {\n            ws.readyState = 3; // Closed state.\n        }\n    }\n    comm.kernel.ws.addEventListener('open', updateReadyState);\n    comm.kernel.ws.addEventListener('close', updateReadyState);\n    comm.kernel.ws.addEventListener('error', updateReadyState);\n\n    ws.close = function () {\n        comm.close();\n    };\n    ws.send = function (m) {\n        //console.log('sending', m);\n        comm.send(m);\n    };\n    // Register the callback with on_msg.\n    comm.on_msg(function (msg) {\n        //console.log('receiving', msg['content']['data'], msg);\n        var data = msg['content']['data'];\n        if (data['blob'] !== undefined) {\n            data = {\n                data: new Blob(msg['buffers'], { type: data['blob'] }),\n            };\n        }\n        // Pass the mpl event to the overridden (by mpl) onmessage function.\n        ws.onmessage(data);\n    });\n    return ws;\n};\n\nmpl.mpl_figure_comm = function (comm, msg) {\n    // This is the function which gets called when the mpl process\n    // starts-up an IPython Comm through the \"matplotlib\" channel.\n\n    var id = msg.content.data.id;\n    // Get hold of the div created by the display call when the Comm\n    // socket was opened in Python.\n    var element = document.getElementById(id);\n    var ws_proxy = comm_websocket_adapter(comm);\n\n    function ondownload(figure, _format) {\n        window.open(figure.canvas.toDataURL());\n    }\n\n    var fig = new mpl.figure(id, ws_proxy, ondownload, element);\n\n    // Call onopen now - mpl needs it, as it is assuming we've passed it a real\n    // web socket which is closed, not our websocket->open comm proxy.\n    ws_proxy.onopen();\n\n    fig.parent_element = element;\n    fig.cell_info = mpl.find_output_cell(\"<div id='\" + id + \"'></div>\");\n    if (!fig.cell_info) {\n        console.error('Failed to find cell for figure', id, fig);\n        return;\n    }\n    fig.cell_info[0].output_area.element.on(\n        'cleared',\n        { fig: fig },\n        fig._remove_fig_handler\n    );\n};\n\nmpl.figure.prototype.handle_close = function (fig, msg) {\n    var width = fig.canvas.width / fig.ratio;\n    fig.cell_info[0].output_area.element.off(\n        'cleared',\n        fig._remove_fig_handler\n    );\n    fig.resizeObserverInstance.unobserve(fig.canvas_div);\n\n    // Update the output cell to use the data from the current canvas.\n    fig.push_to_output();\n    var dataURL = fig.canvas.toDataURL();\n    // Re-enable the keyboard manager in IPython - without this line, in FF,\n    // the notebook keyboard shortcuts fail.\n    IPython.keyboard_manager.enable();\n    fig.parent_element.innerHTML =\n        '<img src=\"' + dataURL + '\" width=\"' + width + '\">';\n    fig.close_ws(fig, msg);\n};\n\nmpl.figure.prototype.close_ws = function (fig, msg) {\n    fig.send_message('closing', msg);\n    // fig.ws.close()\n};\n\nmpl.figure.prototype.push_to_output = function (_remove_interactive) {\n    // Turn the data on the canvas into data in the output cell.\n    var width = this.canvas.width / this.ratio;\n    var dataURL = this.canvas.toDataURL();\n    this.cell_info[1]['text/html'] =\n        '<img src=\"' + dataURL + '\" width=\"' + width + '\">';\n};\n\nmpl.figure.prototype.updated_canvas_event = function () {\n    // Tell IPython that the notebook contents must change.\n    IPython.notebook.set_dirty(true);\n    this.send_message('ack', {});\n    var fig = this;\n    // Wait a second, then push the new image to the DOM so\n    // that it is saved nicely (might be nice to debounce this).\n    setTimeout(function () {\n        fig.push_to_output();\n    }, 1000);\n};\n\nmpl.figure.prototype._init_toolbar = function () {\n    var fig = this;\n\n    var toolbar = document.createElement('div');\n    toolbar.classList = 'btn-toolbar';\n    this.root.appendChild(toolbar);\n\n    function on_click_closure(name) {\n        return function (_event) {\n            return fig.toolbar_button_onclick(name);\n        };\n    }\n\n    function on_mouseover_closure(tooltip) {\n        return function (event) {\n            if (!event.currentTarget.disabled) {\n                return fig.toolbar_button_onmouseover(tooltip);\n            }\n        };\n    }\n\n    fig.buttons = {};\n    var buttonGroup = document.createElement('div');\n    buttonGroup.classList = 'btn-group';\n    var button;\n    for (var toolbar_ind in mpl.toolbar_items) {\n        var name = mpl.toolbar_items[toolbar_ind][0];\n        var tooltip = mpl.toolbar_items[toolbar_ind][1];\n        var image = mpl.toolbar_items[toolbar_ind][2];\n        var method_name = mpl.toolbar_items[toolbar_ind][3];\n\n        if (!name) {\n            /* Instead of a spacer, we start a new button group. */\n            if (buttonGroup.hasChildNodes()) {\n                toolbar.appendChild(buttonGroup);\n            }\n            buttonGroup = document.createElement('div');\n            buttonGroup.classList = 'btn-group';\n            continue;\n        }\n\n        button = fig.buttons[name] = document.createElement('button');\n        button.classList = 'btn btn-default';\n        button.href = '#';\n        button.title = name;\n        button.innerHTML = '<i class=\"fa ' + image + ' fa-lg\"></i>';\n        button.addEventListener('click', on_click_closure(method_name));\n        button.addEventListener('mouseover', on_mouseover_closure(tooltip));\n        buttonGroup.appendChild(button);\n    }\n\n    if (buttonGroup.hasChildNodes()) {\n        toolbar.appendChild(buttonGroup);\n    }\n\n    // Add the status bar.\n    var status_bar = document.createElement('span');\n    status_bar.classList = 'mpl-message pull-right';\n    toolbar.appendChild(status_bar);\n    this.message = status_bar;\n\n    // Add the close button to the window.\n    var buttongrp = document.createElement('div');\n    buttongrp.classList = 'btn-group inline pull-right';\n    button = document.createElement('button');\n    button.classList = 'btn btn-mini btn-primary';\n    button.href = '#';\n    button.title = 'Stop Interaction';\n    button.innerHTML = '<i class=\"fa fa-power-off icon-remove icon-large\"></i>';\n    button.addEventListener('click', function (_evt) {\n        fig.handle_close(fig, {});\n    });\n    button.addEventListener(\n        'mouseover',\n        on_mouseover_closure('Stop Interaction')\n    );\n    buttongrp.appendChild(button);\n    var titlebar = this.root.querySelector('.ui-dialog-titlebar');\n    titlebar.insertBefore(buttongrp, titlebar.firstChild);\n};\n\nmpl.figure.prototype._remove_fig_handler = function (event) {\n    var fig = event.data.fig;\n    if (event.target !== this) {\n        // Ignore bubbled events from children.\n        return;\n    }\n    fig.close_ws(fig, {});\n};\n\nmpl.figure.prototype._root_extra_style = function (el) {\n    el.style.boxSizing = 'content-box'; // override notebook setting of border-box.\n};\n\nmpl.figure.prototype._canvas_extra_style = function (el) {\n    // this is important to make the div 'focusable\n    el.setAttribute('tabindex', 0);\n    // reach out to IPython and tell the keyboard manager to turn it's self\n    // off when our div gets focus\n\n    // location in version 3\n    if (IPython.notebook.keyboard_manager) {\n        IPython.notebook.keyboard_manager.register_events(el);\n    } else {\n        // location in version 2\n        IPython.keyboard_manager.register_events(el);\n    }\n};\n\nmpl.figure.prototype._key_event_extra = function (event, _name) {\n    // Check for shift+enter\n    if (event.shiftKey && event.which === 13) {\n        this.canvas_div.blur();\n        // select the cell after this one\n        var index = IPython.notebook.find_cell_index(this.cell_info[0]);\n        IPython.notebook.select(index + 1);\n    }\n};\n\nmpl.figure.prototype.handle_save = function (fig, _msg) {\n    fig.ondownload(fig, null);\n};\n\nmpl.find_output_cell = function (html_output) {\n    // Return the cell and output element which can be found *uniquely* in the notebook.\n    // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n    // IPython event is triggered only after the cells have been serialised, which for\n    // our purposes (turning an active figure into a static one), is too late.\n    var cells = IPython.notebook.get_cells();\n    var ncells = cells.length;\n    for (var i = 0; i < ncells; i++) {\n        var cell = cells[i];\n        if (cell.cell_type === 'code') {\n            for (var j = 0; j < cell.output_area.outputs.length; j++) {\n                var data = cell.output_area.outputs[j];\n                if (data.data) {\n                    // IPython >= 3 moved mimebundle to data attribute of output\n                    data = data.data;\n                }\n                if (data['text/html'] === html_output) {\n                    return [cell, data, j];\n                }\n            }\n        }\n    }\n};\n\n// Register the function which deals with the matplotlib target/channel.\n// The kernel may be null if the page has been refreshed.\nif (IPython.notebook.kernel !== null) {\n    IPython.notebook.kernel.comm_manager.register_target(\n        'matplotlib',\n        mpl.mpl_figure_comm\n    );\n}\n",
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div id='185b080a-cfb4-4c5f-b3d3-3d83f53dce9c'></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\Users\\vince\\anaconda3\\envs\\geopandas_env\\lib\\site-packages\\pandas\\core\\internals\\blocks.py:351: RuntimeWarning: divide by zero encountered in log10\n",
      "  result = func(self.values, **kwargs)\n"
     ]
    }
   ],
   "source": [
    "TDSHCO3 = chem[['TDSc_mgL','HCO3_mgL']].dropna() #only select TDS & HCO3, drop any rows with missing values\n",
    "plt.figure()\n",
    "plt.loglog(TDSHCO3['TDSc_mgL'],TDSHCO3['HCO3_mgL'],'.k',alpha=0.2)\n",
    "plt.xlabel('TDS (mg/L)')\n",
    "plt.ylabel('HCO3 (mg/L)')\n",
    "TDSHCO3_l = np.log10(TDSHCO3) #log10 transform"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "09d53017",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>TDSc_mgL</th>\n",
       "      <th>HCO3_mgL</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>2.203769</td>\n",
       "      <td>2.210811</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2.270956</td>\n",
       "      <td>2.180581</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>4.488803</td>\n",
       "      <td>2.201397</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>4.279559</td>\n",
       "      <td>2.434569</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>4.278076</td>\n",
       "      <td>2.421604</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>32284</th>\n",
       "      <td>3.433370</td>\n",
       "      <td>2.176091</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>32285</th>\n",
       "      <td>3.558409</td>\n",
       "      <td>2.397940</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>32286</th>\n",
       "      <td>3.456079</td>\n",
       "      <td>2.712429</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>32287</th>\n",
       "      <td>3.684425</td>\n",
       "      <td>2.621383</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>32296</th>\n",
       "      <td>4.619982</td>\n",
       "      <td>1.983716</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>26003 rows × 2 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "       TDSc_mgL  HCO3_mgL\n",
       "0      2.203769  2.210811\n",
       "1      2.270956  2.180581\n",
       "2      4.488803  2.201397\n",
       "3      4.279559  2.434569\n",
       "4      4.278076  2.421604\n",
       "...         ...       ...\n",
       "32284  3.433370  2.176091\n",
       "32285  3.558409  2.397940\n",
       "32286  3.456079  2.712429\n",
       "32287  3.684425  2.621383\n",
       "32296  4.619982  1.983716\n",
       "\n",
       "[26003 rows x 2 columns]"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "TDSHCO3_l"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "f6314391",
   "metadata": {},
   "source": [
    "## Clustering with DBSCAN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f70b6a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.cluster import DBSCAN\n",
    "X = np.array(TDSHCO3_l) #convert to numpy array\n",
    "X = X[~np.isinf(X).any(axis=1)] # remove -inf values\n",
    "db = DBSCAN(eps=0.1, min_samples=10).fit(X)\n",
    "labels = db.labels_\n",
    "# Number of clusters in labels, ignoring noise if present.\n",
    "n_clusters_ = len(set(labels)) - (1 if -1 in labels else 0)\n",
    "n_noise_ = list(labels).count(-1)\n",
    "\n",
    "print(\"Estimated number of clusters: %d\" % n_clusters_)\n",
    "print(\"Estimated number of noise points: %d\" % n_noise_)\n",
    "\n",
    "unique_labels = set(labels)\n",
    "core_samples_mask = np.zeros_like(labels, dtype=bool)\n",
    "core_samples_mask[db.core_sample_indices_] = True\n",
    "\n",
    "colors = [plt.cm.Spectral(each) for each in np.linspace(0, 1, len(unique_labels))]\n",
    "plt.figure()\n",
    "for k, col in zip(unique_labels, colors):\n",
    "    if k == -1:\n",
    "        # Black used for noise.\n",
    "        col = [0, 0, 0, 1]\n",
    "\n",
    "    class_member_mask = labels == k\n",
    "\n",
    "    xy = X[class_member_mask & core_samples_mask]\n",
    "    plt.plot(\n",
    "        xy[:, 0],\n",
    "        xy[:, 1],\n",
    "        \"o\",\n",
    "        markerfacecolor=tuple(col),\n",
    "        markeredgecolor=\"k\",\n",
    "        markersize=4,\n",
    "    )\n",
    "\n",
    "    xy = X[class_member_mask & ~core_samples_mask]\n",
    "    plt.plot(\n",
    "        xy[:, 0],\n",
    "        xy[:, 1],\n",
    "        \"o\",\n",
    "        markerfacecolor=tuple(col),\n",
    "        markeredgecolor=\"k\",\n",
    "        markersize=1,\n",
    "    )\n",
    "\n",
    "plt.title(f\"Estimated number of clusters: {n_clusters_}\\nEstimated number of noise points: {n_noise_}\")\n",
    "plt.xlabel('TDS (mg/L)')\n",
    "plt.ylabel('HCO3 (mg/L)')\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "95a95b58",
   "metadata": {},
   "source": [
    "# Classification\n",
    "Classify samples according to salinity category. Use major ions to estimate class.\n",
    "## Prepare dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "86749e39",
   "metadata": {},
   "outputs": [],
   "source": [
    "# select variables\n",
    "sal = chem[['TDSc_mgL','HCO3_mgL','Na_mgL','K_mgL','Mg_mgL','Ca_mgL','Cl_mgL','SO4_mgL','NO3N_mgL']].dropna()\n",
    "# create labels based on TDS\n",
    "def classifier(row):\n",
    "    if row[\"TDSc_mgL\"] <= 1000:\n",
    "        return 1\n",
    "    elif (row[\"TDSc_mgL\"] > 1000) and (row[\"TDSc_mgL\"] <= 3000):\n",
    "        return 2\n",
    "    elif (row[\"TDSc_mgL\"] > 3000) and (row[\"TDSc_mgL\"] <= 10000):\n",
    "        return 3\n",
    "    else:\n",
    "        return 4\n",
    "sal['class'] = sal.apply(classifier, axis=1)\n",
    "# check labels\n",
    "plt.figure()\n",
    "plt.semilogx(sal['TDSc_mgL'],sal['class'],'.k')\n",
    "plt.xlabel('TDS (mg/L)')\n",
    "plt.ylabel('Label')\n",
    "plt.yticks(np.arange(1,5),['Fresh','Slightly saline','Moderately saline', 'Highly saline'])\n",
    "plt.tight_layout()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "d3dfe113",
   "metadata": {},
   "source": [
    "It is not trivial to predict the salinity class based on major ion composition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8215585c",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(9,9))\n",
    "for i,ion in enumerate(['HCO3_mgL','Na_mgL','K_mgL','Mg_mgL','Ca_mgL','Cl_mgL','SO4_mgL','NO3N_mgL']):\n",
    "    plt.subplot(3,3,i+1)\n",
    "    plt.semilogx(sal[ion],sal['class'],'.k')\n",
    "    plt.title(ion)\n",
    "plt.tight_layout()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "6c5a14e8",
   "metadata": {},
   "source": [
    "## Classify with Random Forests\n",
    "To train a random forest classifier, we need to split the dataset into a training dataset and a validation dataset. Here I chose to use the first 4500 samples to train the classifier. The remainder of the dataset (~1000 samples) is used for validation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "588f2faf",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.ensemble import RandomForestClassifier\n",
    "X = np.array(sal[['HCO3_mgL','Na_mgL','K_mgL','Mg_mgL','Ca_mgL','Cl_mgL','SO4_mgL','NO3N_mgL']])\n",
    "Y = np.array(sal['class'])\n",
    "X_tr = X[0:4500,:]\n",
    "Y_tr = Y[0:4500]\n",
    "X_val = X[4500::,:]\n",
    "Y_val = Y[4500::]\n",
    "clf = RandomForestClassifier(n_estimators=10)\n",
    "clf = clf.fit(X_tr, Y_tr)\n",
    "val = clf.predict(X_val)\n",
    "# compute a summary table with the count of correctly and incorrectly classified samples\n",
    "conf = np.zeros((4,4))\n",
    "for i in np.arange(1,5):\n",
    "    for j in np.arange(1,5):\n",
    "        conf[i-1,j-1]=np.sum((Y_val==i)&(val==j))\n",
    "conff = pd.DataFrame(conf,columns=['RF fresh','RF slightly saline','RF moderately saline', 'RF highly saline'],\n",
    "                     index=['Fresh','Slightly saline','Moderately saline', 'Highly saline'])\n",
    "conff"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "0bad4ba7",
   "metadata": {},
   "source": [
    "## Nitrate class prediction based on TDS and major ions\n",
    "Let's try to do a similar exercise, but now label the dataset based on the nitrate concentration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb3a2ef1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create labels based on Nitrate\n",
    "def classifier(row):\n",
    "    if row[\"NO3N_mgL\"] <= 1:\n",
    "        return 1\n",
    "    elif (row[\"NO3N_mgL\"] > 1) and (row[\"NO3N_mgL\"] <= 25):\n",
    "        return 2\n",
    "    else:\n",
    "        return 3\n",
    "sal['classNO3'] = sal.apply(classifier, axis=1)\n",
    "# check labels\n",
    "plt.figure()\n",
    "plt.semilogx(sal['NO3N_mgL'],sal['classNO3'],'.k')\n",
    "plt.xlabel('NO3 (mg/L)')\n",
    "plt.ylabel('Label')\n",
    "plt.yticks(np.arange(1,4),['No NO3','Low NO3','High NO3'])\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "59e51452",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(9,9))\n",
    "for i,ion in enumerate(['TDSc_mgL','HCO3_mgL','Na_mgL','K_mgL','Mg_mgL','Ca_mgL','Cl_mgL','SO4_mgL']):\n",
    "    plt.subplot(3,3,i+1)\n",
    "    plt.semilogx(sal[ion],sal['classNO3'],'.k')\n",
    "    plt.title(ion)\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b0873824",
   "metadata": {},
   "outputs": [],
   "source": [
    "X = np.array(sal[['TDSc_mgL','HCO3_mgL','Na_mgL','K_mgL','Mg_mgL','Ca_mgL','Cl_mgL','SO4_mgL']])\n",
    "Y = np.array(sal['classNO3'])\n",
    "X_tr = X[0:4500,:]\n",
    "Y_tr = Y[0:4500]\n",
    "X_val = X[4500::,:]\n",
    "Y_val = Y[4500::]\n",
    "clf = RandomForestClassifier(n_estimators=10)\n",
    "clf = clf.fit(X_tr, Y_tr)\n",
    "val = clf.predict(X_val)\n",
    "# compute a summary table with the count of correctly and incorrectly classified samples\n",
    "conf = np.zeros((3,3))\n",
    "for i in np.arange(1,4):\n",
    "    for j in np.arange(1,4):\n",
    "        conf[i-1,j-1]=np.sum((Y_val==i)&(val==j))\n",
    "conff = pd.DataFrame(conf,columns=['RF no NO3','RF low NO3','RF high NO3'],\n",
    "                     index=['no NO3','low NO3','high NO3'])\n",
    "conff"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "3e87b34f",
   "metadata": {},
   "source": [
    "# Regression\n",
    "Instead of predicting classes of salinity, let's try to predict salinity based on major ions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96ae53c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.ensemble import RandomForestRegressor\n",
    "from sklearn.metrics import r2_score\n",
    "\n",
    "X = np.array(sal[['HCO3_mgL','Na_mgL','K_mgL','Mg_mgL','Ca_mgL','Cl_mgL','SO4_mgL','NO3N_mgL']])\n",
    "Y = np.log10(np.array(sal['TDSc_mgL']))\n",
    "X_tr = X[0:4500,:]\n",
    "Y_tr = Y[0:4500]\n",
    "X_val = X[4500::,:]\n",
    "Y_val = Y[4500::]\n",
    "\n",
    "clf = RandomForestRegressor(n_estimators=50)\n",
    "clf = clf.fit(X_tr, Y_tr)\n",
    "val = clf.predict(X_val)\n",
    "r2 = r2_score(Y_val,val)\n",
    "plt.figure()\n",
    "plt.loglog(10**Y_val,10**val,'.k')\n",
    "plt.title(f'Predict TDS from major ions \\n R2 = {r2:4.3f}')\n",
    "plt.xlabel('Measured TDS (mg/L)')\n",
    "plt.ylabel('Modelled TDS (mg/L)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "edf6256d",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.ensemble import RandomForestRegressor\n",
    "\n",
    "X = np.array(sal[['TDSc_mgL','HCO3_mgL','Na_mgL','K_mgL','Mg_mgL','Ca_mgL','Cl_mgL','SO4_mgL']])\n",
    "Y = np.array(sal['NO3N_mgL'])\n",
    "X_tr = X[0:4500,:]\n",
    "Y_tr = Y[0:4500]\n",
    "X_val = X[4500::,:]\n",
    "Y_val = Y[4500::]\n",
    "\n",
    "clf = RandomForestRegressor(n_estimators=50)\n",
    "clf = clf.fit(X_tr, Y_tr)\n",
    "val = clf.predict(X_val)\n",
    "r2 = r2_score(Y_val,val)\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(Y_val,val,'.k')\n",
    "plt.title(f'Predict NO3 from major ions \\n R2 = {r2:4.3f}')\n",
    "plt.xlabel('Measured NO3 (mg/L)')\n",
    "plt.ylabel('Modelled NO3 (mg/L)')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "a9998501",
   "metadata": {},
   "source": [
    "# Dimensionality reduction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74f244ec",
   "metadata": {},
   "outputs": [],
   "source": [
    "cols = ['Lat','Long','TDSc_mgL', 'pH', 'HCO3_mgL', 'Na_mgL', 'K_mgL', 'Mg_mgL', 'Ca_mgL', 'Cl_mgL', 'SO4_mgL','NO3N_mgL']\n",
    "dat = chem[cols].dropna()\n",
    "dat"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "8e32ef07",
   "metadata": {},
   "source": [
    "As we've seen before, the data are very skewed. Before doing any multivariate data analysis it is therefore recommended to normalise data. This is often done by rescaling all variables so their range falls between 0 and 1 or by transforming the variables by subtracting the mean and dividing by the standard deviation. We used log transform earlier for skewed data. This has some drawbacks, especially if there are 0s in the dataset.\n",
    "\n",
    "Another way to normalise data is to calculate the rank, i.e. the position if you were to rank them from smallest to largest. This can be easily done wiith the `rankdata` function from [scipy.stats](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.rankdata.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31ae56e1",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import rankdata\n",
    "# calculate rank of data\n",
    "dr = rankdata(np.array(dat[cols[2::]]),axis=0)\n",
    "dr"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "0330e2d0",
   "metadata": {},
   "source": [
    "Here we'll look into dimensionality reduction and manifold learning. The goal is to find a representation of the data in 2D such that samples that are similar are plotted close to each other and samples that are very different are plotted far apart. The [manifold learning page](https://scikit-learn.org/stable/modules/manifold.html#manifold) gives an overview of methods you can use. The method we'll be using is [t-distributed Stochastic Neighbor Embedding](https://scikit-learn.org/stable/modules/generated/sklearn.manifold.TSNE.html#sklearn.manifold.TSNE).\n",
    "\n",
    "We need to import the function from scikit learn, specify the parameters and then train the algorithm with our dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a16141d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.manifold import TSNE\n",
    "tsne = TSNE(n_components=2, init='pca')\n",
    "X = tsne.fit_transform(dr)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "2702c052",
   "metadata": {},
   "source": [
    "The result is a 2D numpy array with an x and y coordinate for each sample. We can visualise this with `plt.scatter` and color the plot with the rank of each of the variables:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "af84e901",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(8,8))\n",
    "for i,col in enumerate(cols[2::]):\n",
    "    ax = plt.subplot(4,3,i+1, aspect=1)\n",
    "    ax.scatter(X[:,0],X[:,1],0.5,dr[:,i],cmap='viridis')\n",
    "    ax.set_title(col)\n",
    "    ax.set_axis_off()\n",
    "plt.tight_layout()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "a9555e64",
   "metadata": {},
   "source": [
    "This is a spatial dataset, so we want to know show this information on a map. I've developed a 2D perceptually colormap that can be used to assign a unique color to each sample, based on the coordinates of the TSNE projection:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "91acc2bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "def percuniform_rgb(x,y):\n",
    "    '''\n",
    "    Create RGB values for x,y positions from perceptually uniform colour scheme\n",
    "    IN:\n",
    "        x: [nx1] array of x values\n",
    "        y: [nx1] array of y values\n",
    "    OUT:\n",
    "        rgb: [nx3] array of rgb values\n",
    "    '''\n",
    "    # rescale cartesian coordinates into range [-1,1]\n",
    "    # normalise based on max(range(x),range(y))\n",
    "    # multiply by 2 and subtract 1 to have data \n",
    "    # - centered on [0,0] \n",
    "    # - x and y each in range [-1,1]\n",
    "    range_x = x.max()-x.min()\n",
    "    range_y = y.max()-y.min()\n",
    "    range_m = max(range_x,range_y)\n",
    "    x_s = 2*((x-x.min())/range_m)-1\n",
    "    y_s = 2*((y-y.min())/range_m)-1\n",
    "    # load spline interpolant of colour scheme\n",
    "    rgb_interp = np.load('BivariateColourScheme.npy', allow_pickle=True, encoding='latin1').item()\n",
    "    # interpolate rgb values\n",
    "    rgb = np.zeros((len(x),3))\n",
    "    for i,col in enumerate(['R','G','B']):\n",
    "        rgb[:,i] = np.clip(rgb_interp[col].ev(x_s,y_s),0,1)\n",
    "    return(rgb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e25c9d93",
   "metadata": {},
   "outputs": [],
   "source": [
    "tsnergb = percuniform_rgb(X[:,0],X[:,1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37aeb8f3",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,ax = plt.subplots()\n",
    "ax.scatter(X[:,0],X[:,1],5,tsnergb)\n",
    "ax.set_aspect('equal')\n",
    "ax.set_title('Color based on position')\n",
    "ax.set_axis_off()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "861c95ce",
   "metadata": {},
   "source": [
    "We can now make a map of the samples, where each sample is colored based on its location in the TSNE plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ffc90d35",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "fig,ax = plt.subplots()\n",
    "# modify the default figure size\n",
    "fig.set_size_inches(10,10)\n",
    "ax.set_facecolor(\"whitesmoke\")\n",
    "ax.set_title('Color based on TSNE projection')\n",
    "# zoom in to SA data\n",
    "ax.set_xlim(chem['Long'].min(),chem['Long'].max())\n",
    "ax.set_ylim(chem['Lat'].min(),chem['Lat'].max())\n",
    "# use scatter to plot point with color based on recharge\n",
    "s = ax.scatter(dat['Long'],dat['Lat'],15,tsnergb)\n",
    "ax.set_xlabel('Longitude')\n",
    "ax.set_ylabel('Latitude')"
   ]
  }
 ],
 "metadata": {
  "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.16"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": true,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}