{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Detecting and correcting artefacts in RR time series\n\nThis example describes artefacts correction in RR time series.\n\nThe function `correct_rr()` automatically detect artefacts using the method proposed\nby Lipponen & Tarvainen (2019) [#]_. Shorts, extra, long, missed \nand ectopic beats are corrected using either inserion/delection of RR intervals or \ninterpolation of the RR time series, and the\ndetection procedure is run again using cleaned intervals. Importantly, when using \nthis method the signal length can be altered after the interpolation, introducing \nmisalignement with eg. triggers from the experiment. For this reason, it is only \nrecommended to use it in the context of \"bloc design\" study, when heart rate\nvariability is measured for a given time interval (usually > 5 minutes).\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Author: Nicolas Legrand <nicolas.legrand@cfin.au.dk>\n# Licence: Apache-2.0"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\nimport numpy as np\nimport pandas as pd\nimport seaborn as sns\n\nfrom systole import import_dataset1\nfrom systole.correction import correct_rr\nfrom systole.detection import ecg_peaks\nfrom systole.hrv import frequency_domain\nfrom systole.plots import plot_frequency, plot_rr\nfrom systole.utils import input_conversion"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ecg_df = import_dataset1(modalities=['ECG', 'Stim'])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "signal, peaks = ecg_peaks(ecg_df.ecg, method='pan-tompkins', sfreq=1000)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "rr_ms = input_conversion(peaks, input_type=\"peaks\", output_type=\"rr_ms\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "are using Matplotlib as plotting backend.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plot_rr(rr_ms, input_type='rr_ms', figsize=(13, 5))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "peaks and extra peaks by manually increasing or decreasing the length of some RR \nintervals.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "np.random.seed(123)  # For result reproductibility\n\ncorrupted_rr = rr_ms.copy()  # Create a new RR intervals vector\n\n# Randomly select 50 intervals in the time series and multiply them by 2 (missed peaks)\ncorrupted_rr[np.random.choice(len(corrupted_rr), 50)] *= 2\n\n# Randomly select 50 intervals in the time series and divide them by 3 (extra peaks)\ncorrupted_rr[np.random.choice(len(corrupted_rr), 50)] /= 3"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "using `show_artefacts=True` so the artefacts detection runs automatically and shows\nin the plot.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plot_rr(\n    corrupted_rr, input_type='rr_ms', show_artefacts=True, line=False, figsize=(13, 5)\n    )\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "and short RR intervals and they are later correctly detected. We can now apply the RR\ntime series correction method. This function will automatically detect possible\nartefacts in the RR intervals and reconstruct the most probable value using time\nseries interpolation. The number of iteration is set to `2` by default, we add it \nhere for clarity.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "corrected_rr, _ = correct_rr(corrupted_rr)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plot_rr(corrected_rr, input_type='rr_ms', show_artefacts=True,\n        line=False, figsize=(13, 5))\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "This does not means that the new values match exactly the RR intervals, and the new \ncorrected time series will always slightly differs from the original one. However, we\ncan estimate how large this difference is by comparing the true, corrupted and \ncorrected time series a posteriori. Here, instead of comparing the time series side \nby side, we can have a look at some HRV metrics that are known to be affected by RR\nartefacts, like the high frequency HRV.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "_, axs = plt.subplots(1, 3, figsize=(13, 5), sharey=True)\nfor i, rr, lab in zip(range(3), \n                 [rr_ms, corrupted_rr, corrected_rr],\n                 [\"Original\", \"Corrupted\", \"Corrected\"]):\n    plot_frequency(rr, input_type=\"rr_ms\", ax=axs[i])\n    axs[i].set_title(lab)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "similar frequency dynamic. To get a more quantitative view of this resultm e can\nsimply repeat this process of RR corruption-correction many time and check the HF-HRV\nparameters estimated at each steps.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Clean the RR time series before simulation\ncorrected_rr, _ = correct_rr(rr_ms.copy())\n\nsimulation_df = pd.DataFrame([])\nfor i in range(20):\n    \n    # Measure HF-HRV for corrupted RR intervals time series\n    corrupted_rr = corrected_rr.copy()\n    corrupted_rr[np.random.choice(len(corrupted_rr), 50)] *= 2\n    corrupted_rr[np.random.choice(len(corrupted_rr), 50)] /= 3\n    corrupted_hrv = frequency_domain(corrupted_rr, input_type=\"rr_ms\")\n    corrupted_hf = corrupted_hrv[corrupted_hrv.Metric == \"hf_power_nu\"].Values.iloc[0]\n    \n    # Measure HF-HRV for corrected RR intervals time series\n    corrected, _ = correct_rr(corrupted_rr, verbose=False)\n    corrected_hrv = frequency_domain(corrected, input_type=\"rr_ms\")\n    corrected_hf = corrected_hrv[corrected_hrv.Metric == \"hf_power_nu\"].Values.iloc[0]\n\n    simulation_df = pd.concat(\n        [\n            simulation_df, \n            pd.DataFrame(\n                {\"HF-HRV (n.u.)\": [corrupted_hf, corrected_hf],\n                 \"Data Quality\": [\"Corrupted\", \"Corrected\"]\n                }\n                )\n        ]\n    )\n\ninitial_hrv = frequency_domain(corrected_rr, input_type=\"rr_ms\")\ninitial_hf = initial_hrv[initial_hrv.Metric == \"hf_power_nu\"].Values.iloc[0]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt.figure(figsize=(5, 8))\nsns.boxplot(data=simulation_df, x=\"Data Quality\", y=\"HF-HRV (n.u.)\", palette=\"vlag\")\nsns.stripplot(data=simulation_df, x=\"Data Quality\", y=\"HF-HRV (n.u.)\",\n              size=8, color=\".3\", linewidth=0)\nplt.axhline(y=initial_hf, linestyle=\"--\", color=\"gray\")\nplt.title(\"HF-HV recovery \\n after RR artefacts correction\")\nplt.annotate(\"True HF-HRV\", xy=(0, initial_hf), xytext=(-0.4, initial_hf - 0.05),\n             arrowprops = dict(facecolor ='grey', shrink = 0.05))\nsns.despine()\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "affected by the presence of simulated artefacts, the proposed correction methods\nallows to recover the true parameter with some precision.\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## References\n.. [#] Lipponen, J. A., & Tarvainen, M. P. (2019). A robust algorithm for\n  heart rate variability time series artefact correction using novel\n  beat classification. Journal of Medical Engineering & Technology,\n  43(3), 173\u2013181. https://doi.org/10.1080/03091902.2019.1640306\n\n"
      ]
    }
  ],
  "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.11.16"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}