{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fitting features in GITT data\n",
    "\n",
    "In this example, we show how to extract features from GITT data to be used in objective functions instead of directly fitting the data. We use the features defined in [\"Bayesian parameterization of continuum battery models from featurized electrochemical measurements considering noise.\", Kuhn et al., 2023.](https://chemistry-europe.onlinelibrary.wiley.com/doi/full/10.1002/batt.202200374)  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pybamm\n",
    "import pandas as pd\n",
    "import ionworkspipeline as iwp\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.patches import FancyArrowPatch\n",
    "from scipy.integrate import cumulative_trapezoid\n",
    "from functools import partial"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate data"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We begin by generating synthetic data for a half-cell. For this example, we simply generate a few cycles of a GITT experiment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = pybamm.lithium_ion.SPM({\"working electrode\": \"positive\"})\n",
    "\n",
    "parameter_values = iwp.ParameterValues(\"Xu2019\")\n",
    "\n",
    "\n",
    "def j0(c_e, c_s_surf, c_s_max, T):\n",
    "    j0_ref = pybamm.Parameter(\n",
    "        \"Positive electrode reference exchange-current density [A.m-2]\"\n",
    "    )\n",
    "    c_e_init = pybamm.Parameter(\"Initial concentration in electrolyte [mol.m-3]\")\n",
    "\n",
    "    return (\n",
    "        j0_ref\n",
    "        * (c_e / c_e_init) ** 0.5\n",
    "        * (c_s_surf / c_s_max) ** 0.5\n",
    "        * (1 - c_s_surf / c_s_max) ** 0.5\n",
    "    )\n",
    "\n",
    "\n",
    "parameter_values.update(\n",
    "    {\n",
    "        \"Positive electrode reference exchange-current density [A.m-2]\": 5,\n",
    "        \"Positive electrode exchange-current density [A.m-2]\": j0,\n",
    "        \"Positive particle diffusivity [m2.s-1]\": 1e-14,\n",
    "    },\n",
    "    check_already_exists=False,\n",
    ")\n",
    "\n",
    "C_rate = 1 / 15\n",
    "pause_h = 30 / 60 / 60\n",
    "step_h = 0.1\n",
    "rest_h = 0.2\n",
    "V_min = parameter_values[\"Lower voltage cut-off [V]\"]\n",
    "V_max = parameter_values[\"Upper voltage cut-off [V]\"]\n",
    "experiment = pybamm.Experiment(\n",
    "    [\n",
    "        (\n",
    "            f\"Rest for {pause_h} hour (1 second period)\",\n",
    "            f\"Discharge at {C_rate} C for {step_h} hour or until {V_min} V (1 second period)\",\n",
    "            f\"Rest for {rest_h} hour (1 second period)\",\n",
    "        ),\n",
    "    ]\n",
    "    * 3\n",
    ")\n",
    "\n",
    "sim = iwp.Simulation(model, parameter_values=parameter_values, experiment=experiment)\n",
    "_ = sim.solve()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We save the results to a dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "variables = [\"Time [s]\", \"Time [h]\", \"Current [A]\", \"Voltage [V]\"]\n",
    "df = pd.DataFrame(sim.solution.get_data_dict(variables))\n",
    "df = df.rename(columns={\"Cycle\": \"Cycle number\", \"Step\": \"Step number\"})\n",
    "t_h = df[\"Time [h]\"].values\n",
    "currents = df[\"Current [A]\"].values\n",
    "Q = cumulative_trapezoid(currents, t_h, initial=0)\n",
    "df[\"Capacity [A.h]\"] = Q"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extracting features"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we extract the features from the data and save them to a dataframe using the function `iwp.objectives.pulse.gitt_features_objective_variables`. When doing so, we pass in the step numbers which correspond to the pulse and rest steps. \n",
    "\n",
    "Let's look at the features of the first cycle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = df[df[\"Cycle number\"] == 0]\n",
    "data[\"Overpotential [mV]\"] = (data[\"Voltage [V]\"] - data[\"Voltage [V]\"].iloc[0]) * 1000\n",
    "\n",
    "pulse_step_num = 1\n",
    "gitt_features = iwp.objectives.pulse.gitt_features_objective_variables(\n",
    "    data,\n",
    "    step_nums=[pulse_step_num],\n",
    "    full_output=True,\n",
    ")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We then plot the features against the synthetic data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "# plot overpotential during the pulse and rest steps\n",
    "t = data[\"Time [s]\"].values\n",
    "eta = data[\"Overpotential [mV]\"].values\n",
    "ax.plot(t, eta)\n",
    "\n",
    "# plot the GITT features\n",
    "\n",
    "# ohmic drop\n",
    "eta_r = gitt_features[\"U0_gitt_step_1\"][0]\n",
    "x_loc = 200\n",
    "ax.add_artist(\n",
    "    FancyArrowPatch(\n",
    "        (x_loc, eta_r),\n",
    "        (x_loc, 0),\n",
    "        arrowstyle=\"<->\",\n",
    "        color=\"orange\",\n",
    "        mutation_scale=15,\n",
    "        linewidth=2,\n",
    "    )\n",
    ")\n",
    "ax.hlines(eta_r, xmin=t[0], xmax=x_loc, color=\"orange\", linestyle=\":\")\n",
    "ax.annotate(\"ohmic drop\", (x_loc + 25, eta_r / 2), color=\"orange\")\n",
    "\n",
    "# square-root slope (GITT)\n",
    "t_gitt = gitt_features[\"t_gitt_step_1\"]\n",
    "U_gitt = gitt_features[\"U_gitt_step_1\"]\n",
    "ax.plot(t_gitt, U_gitt, color=\"green\", linestyle=\"--\", linewidth=2)\n",
    "ax.annotate(\"square-root \\nslope (GITT)\", (70, -3.5), color=\"green\")\n",
    "\n",
    "# exponential relaxation\n",
    "t_exp = gitt_features[\"t_exp_step_1\"]\n",
    "U_exp = gitt_features[\"U_exp_step_1\"]\n",
    "ax.plot(t_exp, U_exp, color=\"red\", linestyle=\"--\", linewidth=2)\n",
    "ax.annotate(\"exponential \\nrelaxation\", (50, -8), color=\"red\")\n",
    "\n",
    "# square-root slope (ICI)\n",
    "t_ici = gitt_features[\"t_ici_step_1\"]\n",
    "U_ici = gitt_features[\"U_ici_step_1\"]\n",
    "ax.plot(t_ici, U_ici, color=\"purple\", linestyle=\"--\", linewidth=2)\n",
    "ax.annotate(\"square-root \\nslope (ICI)\", (340, -5.5), color=\"purple\")\n",
    "\n",
    "# concentration overpotential\n",
    "eta_c = gitt_features[\"U0_ici_step_1\"][0]\n",
    "x_loc = 600\n",
    "ax.add_artist(\n",
    "    FancyArrowPatch(\n",
    "        (x_loc, eta_c),\n",
    "        (x_loc, 0),\n",
    "        arrowstyle=\"<->\",\n",
    "        color=\"pink\",\n",
    "        mutation_scale=15,\n",
    "        linewidth=2,\n",
    "    )\n",
    ")\n",
    "ax.hlines(eta_c, xmin=t_ici[0], xmax=x_loc, color=\"pink\", linestyle=\":\")\n",
    "ax.annotate(\"concentration \\noverpotential\", (x_loc + 25, eta_c / 2), color=\"pink\")\n",
    "\n",
    "ax.set_xlim([0, 900])\n",
    "ax.set_ylim([-10, 1])\n",
    "ax.set_xlabel(\"Time [s]\")\n",
    "ax.set_ylabel(\"Overpotential [mV]\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fitting the features"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To use these features when fitting data, we can pass the function `gitt_features_objective_variable` to the `Pulse` objective class using the `objective variables` key in the `options` dictionary. Any keyword arguments can be passed to the function using the `objective_variables_kwargs` key in the `options` dictionary. Other useful objective variables functions include `overpotential_objective_variables`, which calculates the overpotentials at the specified times, and `ici_features_objective_variables`, which calculates the features for the Intermittent Current Interruption experiment (see [\"Implementing intermittent current interruption into Li-ion cell modelling for improved battery diagnostics.\", Yin et al., 2022](https://www.sciencedirect.com/science/article/pii/S0013468622010477) for details).\n",
    "\n",
    "Let's compare fitting the features with fitting the data directly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "parameters = {\n",
    "    \"Positive particle diffusivity [m2.s-1]\": iwp.Parameter(\n",
    "        \"D_s\", initial_value=1e-14, bounds=(5e-15, 5e-13)\n",
    "    ),\n",
    "    \"Positive electrode reference exchange-current density [A.m-2]\": iwp.Parameter(\n",
    "        \"j0\", initial_value=5, bounds=(0.1, 10)\n",
    "    ),\n",
    "}\n",
    "parameters_for_pipeline = {\n",
    "    k: v for k, v in parameter_values.items() if k not in parameters.keys()\n",
    "}\n",
    "\n",
    "objective_variables_fns = [\n",
    "    iwp.objectives.pulse.voltage_objective_variables,\n",
    "    partial(\n",
    "        iwp.objectives.pulse.gitt_features_objective_variables, step_nums=[1]\n",
    "    ),  # pulse step\n",
    "]\n",
    "params_fit_list = []\n",
    "for objective_variables_fn in objective_variables_fns:\n",
    "    objectives = iwp.objectives.pulse.get_pulse_objectives_by_cycle(\n",
    "        data,\n",
    "        options={\n",
    "            \"model\": model,\n",
    "            \"objective variables\": objective_variables_fn,\n",
    "        },\n",
    "    )\n",
    "    gitt_half_cell = iwp.DataFit(objectives, parameters=parameters)\n",
    "    results = gitt_half_cell.run(parameters_for_pipeline)\n",
    "    params_fit_list.append(results)\n",
    "    print(results)\n",
    "\n",
    "    gitt_half_cell.plot_fit_results()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For synthetic data, fitting the voltage directly gives better results than fitting the features. However, in practice, fitting the features can be more robust to noise and other model or experimental errors."
   ]
  }
 ],
 "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}