data_vis_and_exploration.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import reciprocal\n",
    "\n",
    "from eis.EISDataIO import ECM_from_raw_strings\n",
    "from utils_preprocessing import eis_dataframe_from_csv\n",
    "from utils_preprocessing import interpolate_to_freq_range, sort_circuits\n",
    "\n",
    "from utils import plot_freq_range, meas_points_print\n",
    "from utils import umap_plots, plot_all_nyquist\n",
    "\n",
    "# autoreload modules when code is run\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "# fix random seed for reproducibility\n",
    "seed = 42\n",
    "np.random.seed(seed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def print_information(df, df_test):\n",
    "    print(f\"Train data:{len(df)} spectra\")\n",
    "    print(f\"Test data:{len(df_test)} spectra\")\n",
    "    print(f\"Total data:{len(df) + len(df_test)} spectra\")\n",
    "\n",
    "    df_all = pd.concat([df, df_test])\n",
    "    for circuit in df_all[\"Circuit\"].unique():\n",
    "        print(f\"{circuit}: {len(df_all[df_all['Circuit'] == circuit])}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the csv into a data frame.\n",
    "filtered = 1\n",
    "save_figs = 0\n",
    "subplots_spectra = 0\n",
    "df_f = eis_dataframe_from_csv(\"data/train_data_filtered.csv\")\n",
    "df_test_f = eis_dataframe_from_csv(\"data/test_data_filtered.csv\")\n",
    "# Print the number of spectra per circuit in the data frame.\n",
    "print_information(df_f, df_test_f)\n",
    "\n",
    "df_uf = eis_dataframe_from_csv(\"data/train_data.csv\")\n",
    "df_test_uf = eis_dataframe_from_csv(\"data/test_data.csv\")\n",
    "# Print the number of spectra per circuit in the data frame.\n",
    "print_information(df_uf, df_test_uf)\n",
    "\n",
    "if filtered: \n",
    "    df = df_f.copy()\n",
    "    df_test = df_test_f.copy()\n",
    "else:\n",
    "    df = df_uf.copy()\n",
    "    df_test = df_test_uf.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The parameter values are drawn from independent reciprocal distributions.\n",
    "# They are a bit tricky to visualize with histograms, thus here's an example:\n",
    "def plot_stats(array, ax, cirucit, param_name, bins=20):\n",
    "    ax.set_title(f\"{cirucit} {param_name}\")\n",
    "    ax.set_xlabel(\"Parameter value\")\n",
    "    # Set x axis to log scale.\n",
    "    ax.set_xscale(\"log\")\n",
    "    ax.set_ylabel(\"Count\")\n",
    "    ax.hist(array, bins=bins)\n",
    "    return ax \n",
    "\n",
    "\n",
    "a, b = 10e-8, 10e-5\n",
    "mean, var, skew, kurt = reciprocal.stats(a, b, moments='mvsk')\n",
    "x = np.linspace(reciprocal.ppf(0.01, a, b), reciprocal.ppf(0.99, a, b), 100)\n",
    "\n",
    "fig, ax = plt.subplots(1, 1)\n",
    "ax.plot(x, reciprocal.pdf(x, a, b), 'r-', lw=1, alpha=0.6, label='reciprocal pdf')\n",
    "ax.set_title(\"Reciprocal distribution, PDF\")\n",
    "\n",
    "fig, ax = plt.subplots(1, 1)\n",
    "r = reciprocal.rvs(a, b, size=500)\n",
    "plot_stats(r, ax, \"\", \"Reciprocal Test\")\n",
    "\n",
    "fig, ax = plt.subplots(1, 1)\n",
    "r = reciprocal.rvs(1.59, 10.1, size=300)\n",
    "# round r to 3 decimals\n",
    "r = np.round(r, 3)\n",
    "plot_stats(r, ax, \"\", \"Reciprocal Test\")\n",
    "\n",
    "# Now draw from a uniform distribution. \n",
    "fig, ax = plt.subplots(1, 1)\n",
    "u = np.random.uniform(0.5, 1.01, size=1000)\n",
    "# round r to 3 decimals\n",
    "u = np.round(u, 3)\n",
    "plot_stats(u, ax, \"\", \"Uniform Test\")\n",
    "\n",
    "# Now draw from a uniform distribution. \n",
    "fig, ax = plt.subplots(1, 1)\n",
    "rt = reciprocal.rvs(0.5, 1.01, size=1000)\n",
    "# round r to 3 decimals\n",
    "rt = np.round(rt, 3)\n",
    "plot_stats(rt, ax, \"\", \"Reciprocal Test\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get all the parameter distributions of the circuits.\n",
    "\n",
    "# Merge the data frames.\n",
    "df_all = df.append(df_test, ignore_index=True)\n",
    "# Sort \n",
    "df_all = sort_circuits(df_all)\n",
    "\n",
    "# turn string into dict of str float pairs\n",
    "parameter_dict = {}\n",
    "max_nb_params = 0\n",
    "for circuit in df_all[\"Circuit\"].unique():\n",
    "    param_names = ECM_from_raw_strings(circuit, df_all.loc[df_all[\"Circuit\"]==circuit].iloc[0].Parameters).param_names\n",
    "    # make empty pandas df with param_names as columns\n",
    "    parameter_dict[circuit] = pd.DataFrame(columns=param_names)\n",
    "    if len(param_names) > max_nb_params:\n",
    "        max_nb_params = len(param_names)\n",
    "    # add rows with param_values\n",
    "    for i, sample in enumerate(df_all.loc[df_all[\"Circuit\"]==circuit].itertuples()):\n",
    "        parameter_dict[circuit].loc[i] = ECM_from_raw_strings(sample.Circuit, sample.Parameters).param_values\n",
    "\n",
    "# Check the covariance between the columns \n",
    "for circuit in df_all[\"Circuit\"].unique():  \n",
    "    print(parameter_dict[circuit].corr())\n",
    "    \n",
    "\n",
    "fig, ax = plt.subplots(9, max_nb_params, figsize=(45, 45))\n",
    "for i, circuit in enumerate(df_all[\"Circuit\"].unique()):\n",
    "    for j, param_name in enumerate(parameter_dict[circuit].columns):\n",
    "        ax[i, j] = plot_stats(parameter_dict[circuit][param_name].values.flatten(), ax[i, j], circuit, param_name)\n",
    "fig.tight_layout()\n",
    "\n",
    "# Conclusions:\n",
    "# 1. Parameters are independent\n",
    "# 2. Parameters are drawn from a reciprocal distribution, but are rounded to 3 significant digits.\n",
    "# 3. The time constants are drwan form a unifrom distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(np.min(parameter_dict[\"RC-RC-RCPE-RCPE\"][\"C1\"]))\n",
    "print(np.max(parameter_dict[\"RC-RC-RCPE-RCPE\"][\"C1\"]))\n",
    "# Comment: This was done in an initial preprocessing of the data because these parameters were missing. \n",
    "# Thus they do not represent the generating distribution of the parameters. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the minimum distance between the time constants of spectra with multiple time constants.\n",
    "\n",
    "# Circuit with multiple time constants: \n",
    "circuits_multiple_t = [\"L-R-RCPE-RCPE\", \"L-R-RCPE-RCPE-RCPE\", \"RC-RC-RCPE-RCPE\", \"RCPE-RCPE\", \"RCPE-RCPE-RCPE\", \"RCPE-RCPE-RCPE-RCPE\"]\n",
    "column_names = [[\"CPE1_t\", \"CPE2_t\"], [\"CPE1_t\", \"CPE2_t\", \"CPE3_t\"], [\"tau1\", \"tau2\", \"CPE1_t\", \"CPE2_t\"], [\"CPE1_t\", \"CPE2_t\"], [\"CPE1_t\", \"CPE2_t\", \"CPE3_t\"], [\"CPE1_t\", \"CPE2_t\", \"CPE3_t\", \"CPE4_t\"]]\n",
    "# count the letter C in the circuit string\n",
    "count_C = lambda x: x.count(\"C\")\n",
    "\n",
    "parameter_dict_time_constants = {}\n",
    "for i, circuit in enumerate(circuits_multiple_t):\n",
    "    t_ind = 0 \n",
    "    parameter_dict_time_constants[circuit] = pd.DataFrame(columns=column_names[i])\n",
    "    print(f\"{circuit}\")\n",
    "    for j, param_name in enumerate(parameter_dict[circuit].columns):\n",
    "        if circuit == \"RC-RC-RCPE-RCPE\":\n",
    "            if param_name==\"R1\":\n",
    "                tau = parameter_dict[circuit][param_name].values.flatten() * parameter_dict[circuit][\"C1\"].values.flatten()\n",
    "                parameter_dict_time_constants[circuit][column_names[i][t_ind]] = tau\n",
    "                t_ind += 1\n",
    "            if param_name==\"R2\":\n",
    "                tau = parameter_dict[circuit][param_name].values.flatten() * parameter_dict[circuit][\"C2\"].values.flatten()\n",
    "                parameter_dict_time_constants[circuit][column_names[i][t_ind]] = tau\n",
    "                t_ind += 1\n",
    "        if param_name[-1] == \"t\":\n",
    "            print(f\"{param_name} {parameter_dict[circuit][param_name].values.flatten()}\")\n",
    "            parameter_dict_time_constants[circuit][column_names[i][t_ind]] = parameter_dict[circuit][param_name].values.flatten()\n",
    "            t_ind += 1\n",
    "        # ax[i, j] = plot_stats(parameter_dict[circuit][param_name].values.flatten(), ax[i, j], circuit, param_name)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "parameter_dict_time_constants[circuit]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from itertools import combinations \n",
    "\n",
    "for key in parameter_dict_time_constants.keys():\n",
    "    tds = list(combinations(parameter_dict_time_constants[key].columns,2))\n",
    "    td_diff= pd.concat([np.abs(parameter_dict_time_constants[key][c[1]].sub(parameter_dict_time_constants[key][c[0]])) for c in tds], axis=1, keys=tds)\n",
    "    # appned td_diff columns to parameter_dict_time_constants[key]\n",
    "    td_diff.columns = td_diff.columns.map('-'.join)\n",
    "    # add \"diff\" to the column names\n",
    "    td_diff.columns = [f\"diff{c}\" for c in td_diff.columns]\n",
    "    parameter_dict_time_constants[key] = pd.concat([parameter_dict_time_constants[key], td_diff], axis=1)\n",
    "\n",
    "# draw 1000 values uniform at random \n",
    "t_example = np.random.uniform(0.5, 1.0, size=1000)\n",
    "t2_example = np.random.uniform(0.5, 1.0, size=1000)\n",
    "td_example = np.abs(t_example - t2_example)\n",
    "\n",
    "fig, ax = plt.subplots(1, 1)\n",
    "ax.hist(td_example, bins=50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "min_diff_taus = np.min(parameter_dict_time_constants[\"RC-RC-RCPE-RCPE\"][\"difftau1-tau2\"])\n",
    "print(f\"Min difference {min_diff_taus}\")\n",
    "diff_taus = np.sort(parameter_dict_time_constants[\"RC-RC-RCPE-RCPE\"][\"difftau1-tau2\"])\n",
    "print(f\"Sorted differences {diff_taus[:10]}\")\n",
    "\n",
    "median_diff_taus = np.median(parameter_dict_time_constants[\"RC-RC-RCPE-RCPE\"][\"difftau1-tau2\"])\n",
    "print(f\"Median difference {median_diff_taus}\")\n",
    "mean_diff_taus = np.mean(parameter_dict_time_constants[\"RC-RC-RCPE-RCPE\"][\"difftau1-tau2\"])\n",
    "print(f\"Mean difference {mean_diff_taus}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "min_diff_ts = np.min(parameter_dict_time_constants[\"L-R-RCPE-RCPE\"][\"diffCPE1_t-CPE2_t\"])\n",
    "print(f\"Min difference {min_diff_ts}\")\n",
    "diff_ts = parameter_dict_time_constants[\"L-R-RCPE-RCPE\"][\"diffCPE1_t-CPE2_t\"]\n",
    "diff_ts = np.sort(diff_ts)\n",
    "print(f\"Sorted differences {diff_ts[:10]}\")\n",
    "median_diff_ts = np.median(parameter_dict_time_constants[\"L-R-RCPE-RCPE\"][\"diffCPE1_t-CPE2_t\"])\n",
    "print(f\"Median difference {median_diff_ts}\")\n",
    "mean_diff_ts = np.mean(parameter_dict_time_constants[\"L-R-RCPE-RCPE\"][\"diffCPE1_t-CPE2_t\"])\n",
    "print(f\"Mean difference {mean_diff_ts}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(6, 6, figsize=(20, 20))\n",
    "for i, key in enumerate(parameter_dict_time_constants.keys()):\n",
    "    # Gett all the columns with \"diff\" in the name\n",
    "    diff_cols = [col for col in parameter_dict_time_constants[key].columns if \"diff\" in col]\n",
    "    print(diff_cols)\n",
    "    ax[i, 0].set_ylabel(f\"{key} \\n Count\")\n",
    "    for j, col in enumerate(diff_cols):\n",
    "        ax[i, j].set_title(f\"|{col[4:]}|\")\n",
    "        ax[i, j].hist(parameter_dict_time_constants[key][col].values.flatten(), bins=50)\n",
    "        if col == \"difftau1-tau2\":\n",
    "            # set the axis to log scale\n",
    "            ax[i, j].set_xscale(\"log\")\n",
    "        #ax[i, j] = plot_stats(parameter_dict_time_constants[key][col].values.flatten(), ax[i, j], \"\", col, bins=50)\n",
    "    if j < 6:\n",
    "        for k in range(j+1, 6):\n",
    "            ax[i, k].set_visible(False)\n",
    "fig.tight_layout()\n",
    "# ax[i, j].set_title(\"Distance between time constants\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# legend\n",
    "# Font size\n",
    "plot_freq_range(df, save=save_figs, verbose=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Include legend \n",
    "meas_points_print(df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Conclusion: There is no correspondance between the measurement points and the circuit \n",
    "# and there is no correspondance between the frequency ranges and the circuit.\n",
    "\n",
    "# Nevertheless, to turn the data into a matrix for machine learning, we need to interpolate to\n",
    "# a common frequency range (i.e., constant frequency values for all samples).\n",
    "interpolated_basis = np.geomspace(10, 1e5, num=30)\n",
    "\n",
    "df_sorted = sort_circuits(interpolate_to_freq_range(df, interpolated_basis))\n",
    "df_test_sorted = sort_circuits(interpolate_to_freq_range(df_test, interpolated_basis))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Umap random state fixed (plots differ slightly from publication, because the seed was previously not fixed).\n",
    "# Even though the seed is fixed, subple changes can occur (probabaly due to the stochastic nature of UMAP.)\n",
    "umap_plots(df_sorted, save=save_figs, random_state=42)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visualize data as images (used for the CNN)\n",
    "if subplots_spectra:\n",
    "    rows = 67 \n",
    "    cols = 111\n",
    "    fig, axs = plt.subplots(rows, cols, sharex=False, sharey=False, figsize=(4*cols, 4*rows), frameon=False)\n",
    "    for i in range(rows):\n",
    "        if i == 0:\n",
    "            axs[i, 0].set_ylabel(\"zreal\")\n",
    "        for j in range(cols):\n",
    "            axs[i, j].set_axis_off()\n",
    "            axs[i, j].plot(df['zreal'][i*cols+j], -df['zimag'][i*cols+j], linewidth=3, color='black')\n",
    "\n",
    "    axs[i, j].set_xlabel('Frequency (Hz)')\n",
    "    if save_figs:\n",
    "        fig.savefig('figures/eis_art.pdf')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Showcase outliers by plotting all spectra in a single plot\n",
    "ax = plot_all_nyquist(df, x_name='zreal', y_name='zimag', title='Spectra Labeled Data')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Normalize all spectra such that zreal and -zimag are in the range [0, 1]\n",
    "\n",
    "df['zreal_norm'] = df['zreal'].apply(lambda x: (x - x.min()) / (x.max() - x.min()))\n",
    "df['zimag_norm'] = df['zimag'].apply(lambda x: (x - x.min()) / (x.max() - x.min()))\n",
    "df['zreal_norm2'] = df['zreal'].apply(lambda x: (x) / (x.max()))\n",
    "df['zimag_norm2'] = df['zimag'].apply(lambda x: (x) / (- x.min()))\n",
    "\n",
    "ax = plot_all_nyquist(\n",
    "    df,\n",
    "    x_name='zreal_norm',\n",
    "    y_name='zimag_norm',\n",
    "    title='Spectra Labeled Data (Normalized)',\n",
    "    linewidth=0.022,\n",
    "    save_name='eis_art2_norm',\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lw = 0.1\n",
    "alpha = 0.2\n",
    "label_fontsize = 16\n",
    "# 5, 2 , 10, 23\n",
    "# Make a 5 row 2 column subplot\n",
    "rows = 3\n",
    "cols = 3\n",
    "fig, axs = plt.subplots(rows, cols, sharex=True, sharey=True, subplot_kw=dict(box_aspect=1), frameon=False)\n",
    "# Set the figure size\n",
    "fig.set_size_inches(15, 15)\n",
    "# unique circuits\n",
    "circuits = np.sort(df['Circuit'].unique())\n",
    "# Iterate over the axes and the circuits\n",
    "for i, ax in enumerate(axs.flat):\n",
    "    ax.set_aspect('equal')\n",
    "    if i <= len(circuits):\n",
    "        # get df for circuit\n",
    "        df_circuit = df[df['Circuit'] == circuits[i]].copy()\n",
    "        ax = plot_all_nyquist(\n",
    "            df_circuit,\n",
    "            ax,\n",
    "            'zreal_norm',\n",
    "            'zimag_norm',\n",
    "            title=circuits[i],\n",
    "            linewidth=lw,\n",
    "            alpha=alpha,\n",
    "            label_fontsize=label_fontsize\n",
    "        )\n",
    "        # if i even remove y labels\n",
    "        if i % cols != 0:\n",
    "            ax.set_ylabel('')\n",
    "            # Remove y ticks\n",
    "            ax.set_yticks([])\n",
    "        # if i not in last row remove x labels\n",
    "        if i not in [i for i in range(len(circuits)) if i+1 > len(circuits) - cols]:\n",
    "            ax.set_xlabel('')\n",
    "            # Remove x ticks\n",
    "            ax.set_xticks([])\n",
    "    else:\n",
    "        ax.set_axis_off()\n",
    "\n",
    "fig.set_tight_layout(True)\n",
    "fig.savefig('figures/spiderweb_eis_circuits.pdf')\n",
    "fig.savefig('figures/spiderweb_eis_circuits.eps')\n",
    "fig.savefig('figures/spiderweb_eis_circuits.png', dpi=300)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.7.13",
   "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.15"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "db3441f3efc18b286252faba8848b9af7d49db2dc5505bf705b457c48ff7d418"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}