Create_optimized_threshold_dict.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Optimized threshold for each output"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This script creates an optimized threshold dictionary per output. \n",
    "\n",
    "It explores SHERPA studies specifically performed to find the optimal threshold for each output. The ***optimal threshold*** *is the most stringent threshold that allows the CausalNN achieving the minim validation loss (i.e., using two-sig-figures of the mse).*\n",
    "\n",
    "In addition it plots the NNs' loss (performance) vs complexity (num. of trainable parameters)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Packages, Modules & Definitions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys; sys.path.append(\"..\")\n",
    "\n",
    "import os, pickle, glob, yaml\n",
    "import numpy    as np\n",
    "import numpy.ma as ma\n",
    "import matplotlib.pyplot as plt\n",
    "from pathlib import Path\n",
    "import csv\n",
    "from math import log10, floor\n",
    "from sklearn    import preprocessing\n",
    "from collections import Counter\n",
    "from os.path import exists\n",
    "from matplotlib.ticker import (MultipleLocator, FormatStrFormatter,\n",
    "                               AutoMinorLocator)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def count_inputs(input_list_file):\n",
    "    c = Counter()\n",
    "    with open(input_list_file) as f:\n",
    "        for line in f:\n",
    "            input = line.strip()\n",
    "            c[input] += 1\n",
    "    return c['1']\n",
    "\n",
    "def calc_parameters(inputs, outputs, nl, nodes):\n",
    "    l_hidden = [nodes]*nl\n",
    "    layers   = l_hidden + [outputs]\n",
    "    parameters = 0\n",
    "    for iL in range(len(layers)):\n",
    "        if iL == 0: parameters += ( inputs * layers[iL] + layers[iL] )\n",
    "        else: parameters += ( layers[iL-1] * layers[iL] + layers[iL] )\n",
    "    return parameters\n",
    "\n",
    "def round_sig(x, sig=2):\n",
    "    return round(x, sig-int(floor(log10(abs(x))))-1)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Options"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SHE_alg    = 'GridSearch'\n",
    "method     = 'sig' # 'sig'; 'aic'\n",
    "sig        = 6\n",
    "pdf        = True\n",
    "\n",
    "## Explore num. of hidden layers and num. of nodes (per layer) with a fixed quantile-thrs (.59)\n",
    "# SHE_path   = f'./hyperparameter_tuning_{SHE_alg}_SHERPA/'\n",
    "# SHE_date   = '20220327_'\n",
    "# # SHE_date   = '20220327_f'\n",
    "# # SHE_date   = '20220327_phq-'\n",
    "# # SHE_date   = '20220327_tphystnd-'\n",
    "\n",
    "## Explore optimal quantile-thrs (.59)\n",
    "SHE_path   = f'./hyperparameter_tuning_{SHE_alg}_threshold_SHERPA/'\n",
    "SHE_date   = '202204*_'\n",
    "# SHE_date   = '20220416_prect'\n",
    "# SHE_date   = '20220417_f*'\n",
    "# SHE_date   = '2022041[6,7]_[p,f][r,l,s]*'\n",
    "# SHE_date   = '20220417_tphystnd-*'\n",
    "# SHE_date   = '20220417_tphystnd-??.*'\n",
    "# SHE_date   = '20220417_phq-??.*'\n",
    "# SHE_date   = '20220417_tphystnd-1??.*'\n",
    "# SHE_date   = '20220417_phq-1??.*'\n",
    "# SHE_date   = '20220417_tphystnd-[2,3,4,5,6,7,8,9]??.*'\n",
    "# SHE_date   = '20220417_phq-[2,3,4,5,6,7,8,9]??.*'\n",
    "# SHE_date   = '20220417_phq-??.*'\n",
    "# SHE_date   = '20220417_phq-???.*'\n",
    "# SHE_date   = '20220417_phq-609.*'\n",
    "# SHE_date   = '20220417_tphystnd-[1,2,3,4,5,6,7]??.*'\n",
    "# SHE_date   = '20220417_tphystnd-[1,2,3,5]?.*'\n",
    "# SHE_date   = '20220417_tphystnd-3.*'\n",
    "# SHE_date   = '2022041[6,7]_*'\n",
    "# SHE_date   = '20220417_tphystnd-524.*'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "thrsType = 'quantile' if pdf else 'spatially'\n",
    "outPath = \"./nn_config/\"+thrsType+\"_dicts/\"\n",
    "Path(outPath).mkdir(parents=True, exist_ok=True)\n",
    "# outFile = outPath+f\"001_{thrsType}_thresholds_sig-{['aic',sig][method=='sig']}.yml\"\n",
    "outFile = outPath+f\"001_{thrsType}_thresholds_mse-{['aic',sig][method=='sig']}.yml\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SHERPA Studies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "thresholds_dict = {}; thresholds_dict['thresholds_dict'] = {}\n",
    "loss_complexity_dict = {}\n",
    "loss       = {'absolute':[],'normalized':[]}\n",
    "thresholds = []\n",
    "\n",
    "SHE_studies = glob.glob(SHE_path+SHE_date+'*')\n",
    "\n",
    "for i, study in enumerate(SHE_studies):\n",
    "    print(f\"Process study: {study}\")\n",
    "    output   = study.split('/')[-1].split('_')[-1]\n",
    "    \n",
    "    loss_complexity_dict[output] = {}\n",
    "    loss_output = []; loss_sig_output = []\n",
    "    threholds_output = []\n",
    "    aic_output = []\n",
    "    num_inputs_output = []\n",
    "    \n",
    "    file = open(f\"{study}/results.csv\")\n",
    "    csvreader = csv.reader(file)\n",
    "    \n",
    "    header = next(csvreader)\n",
    "#     print(header)\n",
    "    for row in csvreader:\n",
    "        if 'COMPLETED' == row[1]:\n",
    "            \n",
    "            threshold = row[5]\n",
    "            input_list_file = glob.glob(study+f\"/{row[0]}/*/*/*/*input_list.txt\")[0]\n",
    "            num_inputs      = count_inputs(input_list_file)\n",
    "#             print(f\"{threshold}-thr: {num_inputs} inputs\")\n",
    "            \n",
    "            cpx = calc_parameters(num_inputs, 1, int(row[3]), int(row[4]))\n",
    "            aic = 2*num_inputs + 1950*np.log(float(row[-1]))\n",
    "            loss_complexity_dict[output][row[0]] = {\n",
    "                 'threshold':threshold,\n",
    "                'num_inputs':num_inputs,\n",
    "                'num_layers':int(row[3]),\n",
    "                 'num_nodes':int(row[4]),\n",
    "                      'loss':float(row[-1]),\n",
    "            'loss_round_sig':round_sig(float(row[-1]),sig),\n",
    "                'complexity':cpx,\n",
    "                       'aic':aic,\n",
    "            }\n",
    "            loss_output.append(float(row[-1]))\n",
    "            loss_sig_output.append(round_sig(float(row[-1]),sig))\n",
    "            threholds_output.append(threshold)\n",
    "            aic_output.append(aic)\n",
    "            num_inputs_output.append(num_inputs)\n",
    "\n",
    "    # All output metrics\n",
    "    loss['absolute'].append(loss_output)\n",
    "    loss_norm       = preprocessing.minmax_scale(loss_output, feature_range=(0, 1), axis=0, copy=True)\n",
    "    loss['normalized'].append(loss_norm)\n",
    "    thresholds.append(threholds_output)\n",
    "\n",
    "    # Best case\n",
    "    if method == 'sig':\n",
    "        loss_output_min = np.min(loss_sig_output)\n",
    "        best_trials     = [\n",
    "            (id,i) for i,id in enumerate(loss_complexity_dict[output]) \n",
    "            if loss_complexity_dict[output][id]['loss_round_sig'] <= loss_output_min\n",
    "        ]\n",
    "        best_thrs  = [(loss_complexity_dict[output][i[0]]['threshold'],i[-1]) for i in best_trials]\n",
    "        best_trial = best_trials[best_thrs.index(max(best_thrs))]\n",
    "    elif method == 'aic':\n",
    "        aic_output      = np.array(aic_output)\n",
    "        aic_output_min  = np.min(aic_output)\n",
    "        best_trials     = [\n",
    "            (id,i) for i,id in enumerate(loss_complexity_dict[output]) \n",
    "            if loss_complexity_dict[output][id]['aic'] == aic_output_min\n",
    "        ]\n",
    "        best_thrs  = [(loss_complexity_dict[output][i[0]]['threshold'],i[-1]) for i in best_trials]\n",
    "        best_trial = best_trials[best_thrs.index(min(best_thrs))]\n",
    "    \n",
    "    # Best output case\n",
    "    loss_complexity_dict[output]['best'] = {\n",
    "              'trial':best_trial[0],\n",
    "          'threshold':loss_complexity_dict[output][best_trial[0]]['threshold'],\n",
    "         'num_inputs':loss_complexity_dict[output][best_trial[0]]['num_inputs'],\n",
    "         'num_layers':loss_complexity_dict[output][best_trial[0]]['num_layers'],\n",
    "          'num_nodes':loss_complexity_dict[output][best_trial[0]]['num_nodes'],\n",
    "               'loss':loss_complexity_dict[output][best_trial[0]]['loss'],\n",
    "     'loss_round_sig':loss_complexity_dict[output][best_trial[0]]['loss_round_sig'],\n",
    "          'loss_norm':loss_norm[best_trial[-1]],\n",
    "    }\n",
    "    \n",
    "    # Optmized threshold dictionary\n",
    "    thresholds_dict['thresholds_dict'][output] = \\\n",
    "    float(loss_complexity_dict[output][best_trial[0]]['threshold'])\n",
    "    \n",
    "    file.close()\n",
    "\n",
    "# Adding manually phq: 3.64; 7.59 # Note those levels are constant and phq (q) is negligible\n",
    "thresholds_dict['thresholds_dict']['phq-3.64'] = .95\n",
    "thresholds_dict['thresholds_dict']['phq-7.59'] = .95\n",
    "\n",
    "loss['absolute']   = np.array(loss['absolute']).flatten()\n",
    "loss['normalized'] = np.array(loss['normalized']).flatten()\n",
    "thresholds         = np.array(thresholds).flatten()\n",
    "\n",
    "\n",
    "# Optmized threshold dictionary: Save?\n",
    "if not exists(outFile):\n",
    "    print(f\"Create {outFile} dictionary.\")\n",
    "    with open(outFile, 'w') as outfile:\n",
    "        yaml.dump(thresholds_dict, outfile, default_flow_style=False)\n",
    "else:\n",
    "    print(f\"{outFile} dictionary exists.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Scatter plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def loss_complexity_scatter_plot(\n",
    "    loss_complexity_dict,\n",
    "    loss,\n",
    "    thresholds,\n",
    "    best=True,\n",
    "    normalized=True,\n",
    "    title='',\n",
    "    lims=False,\n",
    "    save=False,\n",
    "    savenm=False,\n",
    "    **kwargs\n",
    "):\n",
    "      \n",
    "    fig, ax = plt.subplots(1, 1, figsize=(5, 5))\n",
    "\n",
    "    colors_dict = {\n",
    "        'flnt':'orange',\n",
    "        'flns':'orange',\n",
    "        'fsnt':'chocolate',\n",
    "        'fsns':'chocolate',\n",
    "#          'phq':'steelblue',\n",
    "         'phq':'indigo',\n",
    "       'prect':'b',\n",
    "    'tphystnd':'r',\n",
    "    }\n",
    "    markers_dict = {        \n",
    "        1:'1',2:'2',3:'3',4:'>',5:'<',6:'^',7:'v',8:'+',9:'*',10:'.',\n",
    "    }\n",
    "\n",
    "    lims = lims if lims else 1\n",
    "    \n",
    "    # Normalization [0,1]? Note different outputs will have different loss.\n",
    "    if normalized:\n",
    "        loss = loss['normalized']\n",
    "    else:\n",
    "#         loss = np.log(loss['absolute'])\n",
    "        loss = loss['absolute']\n",
    "    \n",
    "    # Scatter Plot\n",
    "    ax.scatter(\n",
    "        thresholds,\n",
    "        loss,\n",
    "        color='grey',\n",
    "    )\n",
    "\n",
    "    # Best outputs case\n",
    "    if best:\n",
    "        for output in sorted(loss_complexity_dict):\n",
    "            if normalized:\n",
    "                loss       = loss_complexity_dict[output]['best']['loss_norm']\n",
    "            else:\n",
    "#                 loss       = loss_complexity_dict[output]['best']['loss_round_sig']\n",
    "                loss       = loss_complexity_dict[output]['best']['loss']\n",
    "#                 loss = np.log(loss_complexity_dict[output]['best']['loss'])\n",
    "            threshold  = loss_complexity_dict[output]['best']['threshold']\n",
    "            num_inputs = loss_complexity_dict[output]['best']['num_inputs']\n",
    "            color  = colors_dict[output.split('-')[0]]\n",
    "            marker = markers_dict[loss_complexity_dict[output]['best']['num_layers']]\n",
    "            ax.scatter(\n",
    "                threshold,\n",
    "                loss,\n",
    "                color=color,\n",
    "                marker=marker,\n",
    "                label=output+f\" ({threshold}: {num_inputs})\",\n",
    "                alpha=.8,\n",
    "            )\n",
    "    \n",
    "    # Attributes\n",
    "    # Make a plot with major ticks that are multiples of 20 and minor ticks that\n",
    "    # are multiples of 5.  Label major ticks with '%d' formatting but don't label\n",
    "    # minor ticks.\n",
    "    ax.xaxis.set_major_locator(MultipleLocator(2))\n",
    "\n",
    "#     if lims != False:\n",
    "#         ax.set_xlim(-.01, lims); ax.set_ylim(-.01, lims)\n",
    "#         ax.axline((-.01,lims) ,(lims, -.01), color='k', linestyle='--')\n",
    "#     else:\n",
    "#         ax.axline((0,lims) ,(lims, 0), color='k', linestyle='--')\n",
    "\n",
    "    if normalized:\n",
    "        ax.set_ylabel(\"Normalized loss (mse)\"); ax.set_xlabel(\"Thresholds\")\n",
    "    else:\n",
    "        ax.set_ylabel(\"Loss (mse)\"); ax.set_xlabel(\"Thresholds\")\n",
    "        ax.set_yscale('log')\n",
    "#     ax.legend(ncol=3,fontsize='medium',bbox_to_anchor=(.5, -.5), loc='lower center')\n",
    "\n",
    "    # Save fig\n",
    "    fig.suptitle(title)\n",
    "    if save is not False:\n",
    "        Path(save).mkdir(parents=True, exist_ok=True)\n",
    "        savenm = f'{statsnm}_single_level-{idx_lab}_nTime-{nTime}.png' if not savenm else savenm\n",
    "        fig.savefig(f\"{save}/{savenm}\",dpi=1000.)\n",
    "        print(f\"{save}/{savenm}\")\n",
    "\n",
    "#         # Legend: colors & markers\n",
    "#         fig_aux = plt.figure(figsize=(5, 5))\n",
    "#         if best:\n",
    "#             for output in sorted(loss_complexity_dict):\n",
    "#                 if normalized:\n",
    "#                     loss       = loss_complexity_dict[output]['best']['loss_norm']\n",
    "#                 else:\n",
    "#                     loss       = loss_complexity_dict[output]['best']['loss']\n",
    "#                 threshold  = loss_complexity_dict[output]['best']['threshold']\n",
    "#                 num_inputs = loss_complexity_dict[output]['best']['num_inputs']\n",
    "#                 color  = colors_dict[output.split('-')[0]]\n",
    "#                 marker = markers_dict[loss_complexity_dict[output]['best']['num_layers']]\n",
    "#                 ax.scatter(\n",
    "#                     threshold,\n",
    "#                     loss,\n",
    "#                     color=color,\n",
    "#                     marker=marker,\n",
    "#                     label=output+f\" ({threshold}: {num_inputs})\",\n",
    "#                     alpha=.8,\n",
    "#                 )\n",
    "        \n",
    "#         for j, jKey in enumerate({**colors_dict, **markers_dict}):\n",
    "#             color = colors_dict[jKey] if jKey > 10 else 'k'\n",
    "#             marker=markers_dict[jKey] if jKey < 10 else 'o'\n",
    "#             lab   = 'nodes-'+str(jKey) if jKey > 10 else 'hl-'+str(jKey)\n",
    "#             plt.scatter(\n",
    "#                 -100.,-100.,\n",
    "#                 label=lab,\n",
    "#                 alpha=.8,\n",
    "#                 color=color,\n",
    "#                 marker=marker,\n",
    "#                 s=50,\n",
    "#             )\n",
    "\n",
    "#         plt.xlim(0, 10)\n",
    "#         plt.ylim(0, 10)\n",
    "#         plt.legend(ncol=3,fontsize='medium',loc='center',frameon=False)\n",
    "#         fig_aux.savefig(f\"{save}/legend_aux.png\",dpi=1000.)\n",
    "            \n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "normalized=False # False; True\n",
    "loss_complexity_scatter_plot(\n",
    "    loss_complexity_dict,\n",
    "    loss,\n",
    "    thresholds,\n",
    "    best=True,\n",
    "    normalized=normalized,\n",
    "#     lims=.5,\n",
    "#     save='./plots', savenm=f\"{SHE_path.split('/')[1]}_all-outputs_sig-{['aic',sig][method=='sig']}_norm-{normalized}.png\",\n",
    "#     save='./plots', savenm=f\"{SHE_path.split('/')[1]}_prect_sig-{['aic',sig][method=='sig']}_norm-{normalized}.png\",\n",
    "#     save='./plots', savenm=f\"{SHE_path.split('/')[1]}_rad_sig-{['aic',sig][method=='sig']}_norm-{normalized}.png\",\n",
    "#     save='./plots', savenm=f\"{SHE_path.split('/')[1]}_dTdt-3hPa-87hPa_sig-{['aic',sig][method=='sig']}_norm-{normalized}.png\",\n",
    "#     save='./plots', savenm=f\"{SHE_path.split('/')[1]}_dQdt-3hPa-87hPa_sig-{['aic',sig][method=='sig']}_norm-{normalized}.png\",\n",
    "#     save='./plots', savenm=f\"{SHE_path.split('/')[1]}_dTdt-100hPa-197hPa_sig-{['aic',sig][method=='sig']}_norm-{normalized}.png\",\n",
    "#     save='./plots', savenm=f\"{SHE_path.split('/')[1]}_dQdt-100hPa-197hPa_sig-{['aic',sig][method=='sig']}_norm-{normalized}.png\",\n",
    "#     save='./plots', savenm=f\"{SHE_path.split('/')[1]}_dTdt-200hPa-992hPa_sig-{['aic',sig][method=='sig']}_norm-{normalized}.png\",\n",
    "#     save='./plots', savenm=f\"{SHE_path.split('/')[1]}_dQdt-200hPa-992hPa_sig-{['aic',sig][method=='sig']}_norm-{normalized}.png\",\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "causalnncam",
   "language": "python",
   "name": "causalnncam"
  },
  "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}