notebooks/ref_data_compare_from_paths.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-08-05T06:20:46.411543Z",
     "start_time": "2022-08-05T06:20:45.480943Z"
    }
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "import shutil\n",
    "import tempfile\n",
    "import subprocess\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from bokeh.plotting import figure, show, ColumnDataSource\n",
    "from bokeh.layouts import gridplot\n",
    "from bokeh.models.tools import HoverTool\n",
    "from bokeh.io import output_notebook\n",
    "\n",
    "output_notebook()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make sure you are in the root of `tardis-refdata` after running the following cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-08-05T06:20:46.422636Z",
     "start_time": "2022-08-05T06:20:46.419904Z"
    }
   },
   "outputs": [],
   "source": [
    "cd .."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define classes and functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-08-05T06:20:47.387178Z",
     "start_time": "2022-08-05T06:20:47.385246Z"
    }
   },
   "outputs": [],
   "source": [
    "# data pickled with protocol 5 can't be opened with `python<3.8.3`, use the backport\n",
    "\n",
    "if sys.version_info < (3, 8, 3):\n",
    "    import pickle5\n",
    "\n",
    "    sys.modules[\"pickle\"] = pickle5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-08-05T06:20:47.934867Z",
     "start_time": "2022-08-05T06:20:47.925248Z"
    }
   },
   "outputs": [],
   "source": [
    "def highlight_missing(val):\n",
    "    if val == True:\n",
    "        return 'background-color: #BCF5A9'\n",
    "    else:\n",
    "        return 'background-color: #F5A9A9'\n",
    "    \n",
    "def highlight_relative_difference(val):\n",
    "    ret = 'background-color: #BCF5A9'\n",
    "    if val is None:\n",
    "        ret = 'background-color: #BCF5A9'\n",
    "    elif val > 1e-2:\n",
    "        ret = 'background-color: #F2F5A9'\n",
    "    elif val > 1e-1:\n",
    "        ret = 'background-color: #F5D0A9'\n",
    "    elif val > 1:\n",
    "        ret = 'background-color: #F5A9A9'\n",
    "    return ret"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-08-05T06:20:48.567357Z",
     "start_time": "2022-08-05T06:20:48.533806Z"
    }
   },
   "outputs": [],
   "source": [
    "class ReferenceComparerFromPaths(object):\n",
    "    def __init__(self, ref1_path, ref2_path, compare_path='unit_test_data.h5'):\n",
    "        self.ref1_fname = ref1_path\n",
    "        self.ref2_fname = ref2_path\n",
    "        self.compare_path = compare_path\n",
    "        \n",
    "    def generate_test_table(self):\n",
    "        rd1_hdfs = pd.HDFStore(self.ref1_fname, mode='r')\n",
    "        rd2_hdfs = pd.HDFStore(self.ref2_fname, mode='r')\n",
    "        rd1_keys = rd1_hdfs.keys()\n",
    "        rd2_keys = rd2_hdfs.keys()\n",
    "        rd1_hdfs.close()\n",
    "        rd2_hdfs.close()\n",
    "        rd1_df = pd.DataFrame(index=rd1_keys, columns=['exists'])\n",
    "        rd2_df = pd.DataFrame(index=rd2_keys, columns=['exists'])\n",
    "        rd1_df['exists'] = True\n",
    "        rd2_df['exists'] = True\n",
    "        joined_df = rd1_df.join(rd2_df, how='outer', lsuffix='_1', rsuffix='_2')\n",
    "        joined_df = joined_df.fillna(False)\n",
    "        return joined_df\n",
    "    \n",
    "    def compare_refdata(self, test_table):\n",
    "        test_table['match'] = None\n",
    "        test_table['abs_diff_mean'] = None\n",
    "        test_table['abs_diff_max'] = None\n",
    "        test_table['rel_diff_mean'] = None\n",
    "        test_table['rel_diff_max'] = None\n",
    "        for row_id, row in test_table.iterrows():\n",
    "            \n",
    "            if row[['exists_1', 'exists_2']].all():\n",
    "                ref1_df = pd.read_hdf(self.ref1_fname, row_id)\n",
    "                ref2_df = pd.read_hdf(self.ref2_fname, row_id)\n",
    "                \n",
    "                if isinstance(ref1_df, pd.Series):\n",
    "                    try:\n",
    "                        pd.testing.assert_series_equal(ref1_df, ref2_df)\n",
    "                    except AssertionError:\n",
    "                        test_table.loc[row_id, 'match'] = False\n",
    "                        abs_diff = np.fabs(ref1_df - ref2_df)\n",
    "                        rel_diff = (abs_diff / np.fabs(ref1_df))[ref1_df != 0]\n",
    "                        test_table.loc[row_id, 'abs_diff_mean'] = abs_diff.mean()\n",
    "                        test_table.loc[row_id, 'abs_diff_max'] = abs_diff.max()\n",
    "                        test_table.loc[row_id, 'rel_diff_mean'] = rel_diff.mean()\n",
    "                        test_table.loc[row_id, 'rel_diff_max'] = rel_diff.max()\n",
    "                    else:\n",
    "                        test_table.loc[row_id, 'match'] = True\n",
    "\n",
    "                elif isinstance(ref1_df, pd.DataFrame):\n",
    "                    try:\n",
    "                        pd.testing.assert_frame_equal(ref1_df, ref2_df)\n",
    "                    except AssertionError:\n",
    "                        test_table.loc[row_id, 'match'] = False\n",
    "                        abs_diff = np.fabs(ref1_df - ref2_df)\n",
    "                        rel_diff = (abs_diff / np.fabs(ref1_df))[ref1_df != 0]\n",
    "                        test_table.loc[row_id, 'abs_diff_mean'] = abs_diff.mean(skipna=True).mean()\n",
    "                        test_table.loc[row_id, 'abs_diff_max'] = abs_diff.max(skipna=True).max()\n",
    "                        test_table.loc[row_id, 'rel_diff_mean'] = rel_diff.mean(skipna=True).mean()\n",
    "                        test_table.loc[row_id, 'rel_diff_max'] = rel_diff.max(skipna=True).max()\n",
    "                    else:\n",
    "                        test_table.loc[row_id, 'match'] = True\n",
    "\n",
    "                else:\n",
    "                    raise ValueError('Needs to be a Series or DataFrame but is' + str(type(ref1_df)))\n",
    "        return test_table\n",
    "        \n",
    "    def generate_test_table_from_direct_paths(self):\n",
    "        tt = self.generate_test_table()\n",
    "        tt = self.compare_refdata(tt)\n",
    "        return tt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-08-05T06:20:51.367171Z",
     "start_time": "2022-08-05T06:20:51.330209Z"
    }
   },
   "outputs": [],
   "source": [
    "REF1_PATH = os.environ['REF1_PATH']\n",
    "\n",
    "if not REF1_PATH:\n",
    "    raise ValueError\n",
    "\n",
    "REF2_PATH = os.environ['REF2_PATH']\n",
    "\n",
    "if not REF2_PATH:\n",
    "    raise ValueError"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-08-05T06:21:40.146638Z",
     "start_time": "2022-08-05T06:21:40.121845Z"
    }
   },
   "outputs": [],
   "source": [
    "REF1_PATH, REF2_PATH"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Test table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-08-08T11:35:59.118098Z",
     "start_time": "2022-08-08T11:35:46.084963Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "comparer = ReferenceComparerFromPaths(\n",
    "    ref1_path=REF1_PATH, \n",
    "    ref2_path=REF2_PATH, \n",
    ")\n",
    "tt = comparer.generate_test_table_from_direct_paths()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-08-08T11:38:38.696220Z",
     "start_time": "2022-08-08T11:38:38.679872Z"
    }
   },
   "outputs": [],
   "source": [
    "tt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-08-05T06:22:19.228480Z",
     "start_time": "2022-08-05T06:22:18.305631Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "tt[[\"exists_1\", \"exists_2\", 'rel_diff_mean', 'rel_diff_max', 'match']].style.applymap(\n",
    "    highlight_missing, subset=['exists_1', 'exists_2', 'match']).applymap(\n",
    "    highlight_relative_difference, subset=['rel_diff_mean', 'rel_diff_max'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Detailed inspection of the reference data\n",
    "\n",
    "If parts of the reference data show differences between revisions, you should invest some time examining these differences in detail. Often, visualizing the relevant data blocks already helps. \n",
    "\n",
    "You can use the following plotting routines as a blueprint and adjust and extend them to your needs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-08-08T11:40:40.990061Z",
     "start_time": "2022-08-08T11:40:40.975288Z"
    }
   },
   "outputs": [],
   "source": [
    "def compare_output_nu(df1, df2, mpl_backend=False):\n",
    "    nu_min = np.min([df1.min(), df2.min()])\n",
    "    nu_max = np.max([df1.max(), df2.max()])\n",
    "    \n",
    "    if mpl_backend:\n",
    "        plt.figure(figsize=(14, 6))\n",
    "        plt.subplot(121)\n",
    "        plt.plot(df1, df2, ',')\n",
    "        plt.xlabel(\"output_nu, ref 1\")\n",
    "        plt.ylabel(\"output_nu, ref 2\")\n",
    "        plt.subplot(122)\n",
    "        plt.hist(df1, bins=np.linspace(nu_min, nu_max, 100), histtype=\"step\", label=\"ref 1\")\n",
    "        plt.hist(df2, bins=np.linspace(nu_min, nu_max, 100), histtype=\"step\", label=\"ref 2\")\n",
    "        plt.xlabel(\"output_nu\")\n",
    "        plt.legend(frameon=False)\n",
    "        \n",
    "        return\n",
    "\n",
    "    \n",
    "    TOOLTIPS = [(\"(x,y)\", \"(@x, @y)\")]\n",
    "    hover = HoverTool(tooltips=TOOLTIPS)\n",
    "    \n",
    "    p = figure()\n",
    "    output_nu = ColumnDataSource(pd.DataFrame.from_records({'x': df1.values, \n",
    "                                                            'y': df2.values}))\n",
    "    p.circle('x', 'y', size=1, source=output_nu)\n",
    "    p.xaxis.axis_label = \"output_nu, ref 1\"\n",
    "    p.yaxis.axis_label = \"output_nu, ref 2\"\n",
    "    p.xaxis.formatter.precision = 1\n",
    "    p.yaxis.formatter.precision = 1\n",
    "    p.add_tools(hover)\n",
    "\n",
    "    # Step lines are hacky way to make histograms with Bokeh\n",
    "    arr_hist1, edges1 = np.histogram(df1.values, \n",
    "                                     bins = 100, \n",
    "                                     range = [nu_min, nu_max])\n",
    "    arr_hist2, edges2 = np.histogram(df1.values, \n",
    "                                     bins = 100, \n",
    "                                     range = [nu_min, nu_max])\n",
    "    \n",
    "    hist1 = ColumnDataSource(pd.DataFrame.from_records({'x': np.linspace(nu_min, nu_max, 100),\n",
    "                                                        'y': arr_hist1}))\n",
    "    hist2 = ColumnDataSource(pd.DataFrame.from_records({'x': np.linspace(nu_min, nu_max, 100),\n",
    "                                                        'y': arr_hist2}))\n",
    "    q = figure()\n",
    "    q.step('x', 'y', source=hist1, legend_label='ref 1')\n",
    "    q.step('x', 'y', source=hist2, legend_label='ref 2', color='#ff7f0e')\n",
    "    q.xaxis.axis_label = \"output_nu\"\n",
    "    q.xaxis.formatter.precision = 1\n",
    "    q.legend.click_policy=\"hide\"\n",
    "    \n",
    "    # Currently HoverTool does not work for step line glyph. See: https://github.com/bokeh/bokeh/issues/7419\n",
    "    q.add_tools(hover)\n",
    "    \n",
    "    plot = gridplot([p, q], ncols=2, plot_width=420, plot_height=360)\n",
    "    \n",
    "    show(plot)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-08-08T11:40:41.574233Z",
     "start_time": "2022-08-08T11:40:41.554823Z"
    }
   },
   "outputs": [],
   "source": [
    "def compare_spectrum(ref1_nu, ref1_L, ref2_nu, ref2_L, mpl_backend=False):\n",
    "    \n",
    "    if mpl_backend:\n",
    "        plt.figure(figsize=(14, 6))\n",
    "        plt.subplot(121)\n",
    "        plt.plot(ref1_nu, ref1_L, label=\"ref 1\")\n",
    "        plt.plot(ref2_nu, ref2_L, label=\"ref 2\")\n",
    "        plt.xlabel(\"nu\")\n",
    "        plt.ylabel(\"L\")\n",
    "        plt.legend(frameon=False)\n",
    "        plt.subplot(122)\n",
    "        plt.plot(ref1_nu, ref1_L / ref2_L)\n",
    "        plt.xlabel(\"nu\")\n",
    "        plt.ylabel(\"L ref 1 / L ref 2\")\n",
    "        \n",
    "        return\n",
    "    \n",
    "    \n",
    "    TOOLTIPS = [(\"(x,y)\", \"(@x, @y)\")]\n",
    "    hover = HoverTool(tooltips=TOOLTIPS)\n",
    "    \n",
    "    p = figure()\n",
    "    spectrum1 = ColumnDataSource(pd.DataFrame.from_records({'x': ref1_nu.values, \n",
    "                                                            'y': ref1_L}))\n",
    "    spectrum2 = ColumnDataSource(pd.DataFrame.from_records({'x': ref2_nu.values, \n",
    "                                                            'y': ref2_L}))\n",
    "    p.line('x', 'y', source=spectrum1, legend_label='ref 1')\n",
    "    p.line('x', 'y', source=spectrum2, legend_label='ref 2', color='#ff7f0e')\n",
    "    p.xaxis.axis_label = \"L\"\n",
    "    p.yaxis.axis_label = \"nu\"\n",
    "    p.xaxis.formatter.precision = 1\n",
    "    p.yaxis.formatter.precision = 1\n",
    "    p.legend.click_policy=\"hide\"\n",
    "    p.add_tools(hover)\n",
    "    \n",
    "    q = figure()\n",
    "    lum_ratio = ColumnDataSource(pd.DataFrame.from_records({'x': ref1_nu.values, \n",
    "                                                            'y': ref1_L.values/ref2_L.values}))\n",
    "    q.circle('x', 'y', size=1, source=lum_ratio)\n",
    "    q.xaxis.axis_label = \"nu\"\n",
    "    q.yaxis.axis_label = \"L ref 1 / L ref 2\"\n",
    "    q.xaxis.formatter.precision = 1\n",
    "    q.yaxis.formatter.precision = 1\n",
    "    q.add_tools(hover)\n",
    "    \n",
    "    \n",
    "    plot = gridplot([p, q], ncols=2, plot_width=420, plot_height=360)\n",
    "    \n",
    "    show(plot)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get the data and find all the entries for which differences exist:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-08-08T11:40:42.917613Z",
     "start_time": "2022-08-08T11:40:42.905070Z"
    }
   },
   "outputs": [],
   "source": [
    "tmp1 = pd.HDFStore(comparer.ref1_fname, \"r\")\n",
    "tmp2 = pd.HDFStore(comparer.ref2_fname, \"r\")\n",
    "\n",
    "diff_entries = tt.loc[(tt[\"match\"] == False) & (tt[\"exists_1\"] == True) & (tt[\"exists_2\"] == True)].index"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-08-08T11:40:44.316699Z",
     "start_time": "2022-08-08T11:40:44.156168Z"
    }
   },
   "outputs": [],
   "source": [
    "compare_output_nu(tmp1['/test_simulation/output_nu'], tmp2['/test_simulation/output_nu'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "compare_spectrum(tmp1['/test_runner_simple/spectrum/_frequency'][:-1], \n",
    "                 tmp1['/test_runner_simple/spectrum/luminosity'],\n",
    "                 tmp2['/test_runner_simple/spectrum/_frequency'][:-1], \n",
    "                 tmp2['/test_runner_simple/spectrum/luminosity'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "compare_spectrum(tmp1['/test_runner_simple_integral_macroatom_interp/spectrum/_frequency'][:-1], \n",
    "                 tmp1['/test_runner_simple_integral_macroatom_interp/spectrum_integrated/luminosity'],\n",
    "                 tmp2['/test_runner_simple_integral_macroatom_interp/spectrum/_frequency'][:-1], \n",
    "                 tmp2['/test_runner_simple_integral_macroatom_interp/spectrum_integrated/luminosity'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "compare_spectrum(tmp1['/test_runner_simple_integral_macroatom/spectrum/_frequency'][:-1], \n",
    "                 tmp1['/test_runner_simple_integral_macroatom/spectrum_integrated/luminosity'],\n",
    "                 tmp2['/test_runner_simple_integral_macroatom/spectrum/_frequency'][:-1], \n",
    "                 tmp2['/test_runner_simple_integral_macroatom/spectrum_integrated/luminosity'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "compare_spectrum(tmp1['/test_runner_simple_integral_downbranch/spectrum/_frequency'][:-1], \n",
    "                 tmp1['/test_runner_simple_integral_downbranch/spectrum_integrated/luminosity'],\n",
    "                 tmp2['/test_runner_simple_integral_downbranch/spectrum/_frequency'][:-1], \n",
    "                 tmp2['/test_runner_simple_integral_downbranch/spectrum_integrated/luminosity'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "compare_spectrum(tmp1['/test_runner_simple/spectrum_virtual/_frequency'][:-1], \n",
    "                 tmp1['/test_runner_simple/spectrum_virtual/luminosity'],\n",
    "                 tmp2['/test_runner_simple/spectrum_virtual/_frequency'][:-1], \n",
    "                 tmp2['/test_runner_simple/spectrum_virtual/luminosity'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if False in tt['match']:\n",
    "    os.environ[\"reference_comparision\"] = \"failed\""
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.13"
  },
  "toc": {
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": "block",
   "toc_window_display": false
  },
  "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": false
  },
  "vscode": {
   "interpreter": {
    "hash": "4a6df93559b6ae16f50884f690d4cbca3f2c8dfd9034cf12f41f359d9acc97fc"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}