{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "8af9a4bc",
   "metadata": {},
   "source": [
    "# Example: Preprocessing of a full dataset with multiple subjects\n",
    "\n",
    "This is an example of a real dataset of a study from our cognitive neuroscience lab at the University of Tromsø.\n",
    "\n",
    "The notebook cannot be run without the underlying data (which can be found at <https://osf.io/7wcej/>) but it can be used as an example to adapt for your own studies. To load and share the data of a study more easily, consider using a [configuration file](../user_guide/studies.rst)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "326f82b4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import glob\n",
    "import pypillometry as pp\n",
    "from pypillometry import u\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "\n",
    "plt.rcParams.update({\"figure.max_open_warning\": 0})  # suppress a warning"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "af456208",
   "metadata": {},
   "outputs": [],
   "source": [
    "# load all raw datasets\n",
    "exclude = [\"20_pav\", \"8_pav\", \"1_pav_pilot\"]  # these subjects did not have usable data\n",
    "\n",
    "# here, data is already stored in pypillometry's own format as .pd files\n",
    "pdfiles = glob.glob(\"data/eyedata/*.pd\")\n",
    "datasets = [pp.EyeData.from_file(fname) for fname in pdfiles]\n",
    "\n",
    "# datasets are stored in a dict with the subject number as key\n",
    "datasets = {d.name.split(\"_\")[0]: d for d in datasets if d.name not in exclude}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7d6d84bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# preprocessing parameters (durations in ms)\n",
    "default_param = {\"min_duration\": 10,        # min duration of a blink\n",
    "                 \"vel_onset_min_duration\": 6,   # duration of onset/offset for blink detection\n",
    "                 \"vel_offset_min_duration\": 4,\n",
    "                 \"vel_onset\": -5,           # velocity thresholds for onset and offset\n",
    "                 \"vel_offset\": 5,\n",
    "                 \"strategies\": [\"zero\", \"velocity\"],  # strategies for blink detection\n",
    "                 \"distance\": 100,           # minimum distance between two blinks\n",
    "                 \"margin\": (50, 150),       # margins for Mahot algorithm\n",
    "                 \"cutoff\": 5}               # lowpass-filter cutoff (Hz)\n",
    "\n",
    "# all subjects get the same set of default parameters initially\n",
    "params = {subj: default_param.copy() for subj in datasets.keys()}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e0609837",
   "metadata": {},
   "outputs": [],
   "source": [
    "# fine-tuning of the parameters per subject (this is done iteratively by\n",
    "# inspecting the per-subject PDF files generated below and changing parameters)\n",
    "params[\"4\"].update(vel_onset=-10, min_duration=5)\n",
    "params[\"6\"].update(strategies=[\"zero\"])\n",
    "params[\"7\"].update(vel_onset=-10)\n",
    "params[\"9\"].update(vel_onset=-10)\n",
    "params[\"10\"].update(vel_onset=-10)\n",
    "params[\"12\"].update(margin=(50, 200), vel_onset=-10)\n",
    "params[\"13\"].update(vel_onset=-10)\n",
    "params[\"14\"].update(margin=(50, 220), vel_onset=-10)\n",
    "params[\"15\"].update(margin=(50, 350), vel_onset=-10)\n",
    "params[\"16\"].update(margin=(50, 300), vel_onset=-10, min_duration=5)\n",
    "params[\"17\"].update(vel_onset=-10)\n",
    "params[\"21\"].update(vel_onset=-100, vel_offset=100,\n",
    "                    vel_onset_min_duration=4, vel_offset_min_duration=4)\n",
    "params[\"23\"].update(vel_onset=-10, vel_onset_min_duration=4, min_duration=5)\n",
    "params[\"24\"].update(margin=(50, 320), vel_onset=-10, min_duration=1)\n",
    "params[\"25\"].update(margin=(50, 350), vel_onset=-10)\n",
    "params[\"29\"].update(margin=(100, 400), vel_onset=-10)\n",
    "params[\"30\"].update(margin=(100, 300), vel_onset=-20)\n",
    "params[\"31\"].update(margin=(100, 300), vel_onset=-10)\n",
    "params[\"32\"].update(vel_onset=-10)\n",
    "params[\"34\"].update(margin=(70, 350), vel_onset=-10)\n",
    "params[\"35\"].update(margin=(50, 250), vel_onset=-10, min_duration=5)\n",
    "params[\"36\"].update(margin=(100, 300), vel_onset=-10, min_duration=5)\n",
    "params[\"37\"].update(margin=(80, 250), vel_onset=-10, min_duration=1)\n",
    "params[\"39\"].update(vel_onset=-10, min_duration=5)\n",
    "params[\"40\"].update(margin=(60, 200), vel_onset=-10)\n",
    "params[\"41\"].update(margin=(60, 230), vel_onset=-20, min_duration=3)\n",
    "params[\"42\"].update(vel_onset=-10, min_duration=0)\n",
    "params[\"43\"].update(vel_onset=-10, min_duration=0)\n",
    "params[\"46\"].update(margin=(80, 220), vel_onset=-20, min_duration=0)\n",
    "params[\"47\"].update(vel_onset=-10)\n",
    "params[\"49\"].update(vel_onset=-10)\n",
    "params[\"50\"].update(vel_onset=-10, min_duration=0)\n",
    "params[\"51\"].update(margin=(50, 220), vel_onset=-10)\n",
    "params[\"53\"].update(vel_onset=-10, min_duration=0)\n",
    "params[\"54\"].update(margin=(50, 220), vel_onset=-10)\n",
    "params[\"55\"].update(margin=(50, 250), vel_onset=-10)\n",
    "params[\"56\"].update(vel_onset=-10, min_duration=0)\n",
    "params[\"57\"].update(margin=(80, 150), vel_onset=-10, min_duration=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "775382e2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# the preprocessing pipeline, with the parameters of one subject\n",
    "def preproc(pars):\n",
    "    return (pp.pupil.blinks_detect(min_duration=pars[\"min_duration\"],\n",
    "                                   strategies=pars[\"strategies\"],\n",
    "                                   vel_onset=pars[\"vel_onset\"], vel_offset=pars[\"vel_offset\"],\n",
    "                                   vel_onset_min_duration=pars[\"vel_onset_min_duration\"],\n",
    "                                   vel_offset_min_duration=pars[\"vel_offset_min_duration\"])\n",
    "            | pp.pupil.blinks_merge_close(distance=pars[\"distance\"])\n",
    "            | pp.pupil.blinks_interpolate(margin=pars[\"margin\"],\n",
    "                                          vel_onset=pars[\"vel_onset\"], vel_offset=pars[\"vel_offset\"])\n",
    "            | pp.pupil.lowpass_filter(cutoff=pars[\"cutoff\"]))\n",
    "\n",
    "# run the preprocessing pipeline on all subjects and produce PDFs for inspection\n",
    "preprocs = {}\n",
    "for subj, d in datasets.items():\n",
    "    print(d.name)\n",
    "    dp = d | preproc(params[subj])\n",
    "    dp.plot.plot_pupil_segments(pdf_file=\"pics/%s.pdf\" % d.name, interv=1 * u.min)\n",
    "    preprocs[subj] = dp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "295e26e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# write down notes about each subject when going through the PDFs\n",
    "notes = {\n",
    "    \"1\": \"last 5 mins missing; else quite ok, a few semi-blinks not properly detected\",\n",
    "    \"2\": \"many blinks after min 17; but reconstruction works very well\",\n",
    "    \"3\": \"fine\",\n",
    "    \"4\": \"fine (after tuning); a couple of short blinks that do not go all the way down to zero\",\n",
    "    \"5\": \"an awful lot of blinks; reconstruction seems to work but subj should probably be excluded\",\n",
    "    \"6\": \"fine (after tuning); very nice data with few blinks\",\n",
    "    \"7\": \"ok after tuning; weird stuff min 32-33\",\n",
    "    \"8\": \"missing data\",\n",
    "    \"9\": \"fine (after tuning); worse data in second half but ok\",\n",
    "    \"10\": \"many blinks, probably exclude\",\n",
    "    \"11\": \"usable but quite a few blinks (especially in the last part from min 33)\",\n",
    "    \"12\": \"fine after tuning but many blinks in second half\",\n",
    "    \"13\": \"usable but quite a few blinks; some segments are bad but should be ok due to exclusion crit\",\n",
    "    \"14\": \"fine\",\n",
    "    \"15\": \"exclude, many blinks and long opening-period after blinks\",\n",
    "    \"16\": \"fine; very nice data with few blinks\",\n",
    "    \"17\": \"usable but quite a few blinks\",\n",
    "    \"18\": \"missing\",\n",
    "    \"19\": \"usable (borderline), many blinks but nice reconstruction\",\n",
    "    \"20\": \"missing\",\n",
    "    \"21\": \"fine; except from min 37; incredibly fast transients!\",\n",
    "    \"22\": \"exclude; way to many blinks to do anything with this\",\n",
    "    \"23\": \"fine\",\n",
    "    \"24\": \"fine; very slow opening-transient\",\n",
    "    \"25\": \"borderline; a lot of blinks (more towards end) and a few weird artifacts in the second half\",\n",
    "    \"26\": \"fine; great, clean data!\",\n",
    "    \"27\": \"exclude\",\n",
    "    \"28\": \"missing\",\n",
    "    \"29\": \"usable; slow opening-transient\",\n",
    "    \"30\": \"fine\",\n",
    "    \"31\": \"fine\",\n",
    "    \"32\": \"fine\",\n",
    "    \"33\": \"usable; many blinks in second half\",\n",
    "    \"34\": \"fine\",\n",
    "    \"35\": \"borderline; lots of blinks during later segments\",\n",
    "    \"36\": \"borderline; esp near the end\",\n",
    "    \"37\": \"borderline; blinks consist of many consecutive, short-blinks\",\n",
    "    \"38\": \"missing\",\n",
    "    \"39\": \"fine (after tuning)\",\n",
    "    \"40\": \"borderline; many blinks throughout\",\n",
    "    \"41\": \"fine\",\n",
    "    \"42\": \"fine; has some of the weird triple-blinks\",\n",
    "    \"43\": \"fine; missing from min 34 (last block?)\",\n",
    "    \"44\": \"fine\",\n",
    "    \"45\": \"exclude; many blinks throughout\",\n",
    "    \"46\": \"usable\",\n",
    "    \"47\": \"fine\",\n",
    "    \"48\": \"missing\",\n",
    "    \"49\": \"fine\",\n",
    "    \"50\": \"fine\",\n",
    "    \"51\": \"fine\",\n",
    "    \"52\": \"fine\",\n",
    "    \"53\": \"fine\",\n",
    "    \"54\": \"fine\",\n",
    "    \"55\": \"usable; borderline towards the end\",\n",
    "    \"56\": \"usable, worse towards end\",\n",
    "    \"57\": \"exclude\"\n",
    "}\n",
    "exclude_preproc = [\"5\", \"10\", \"15\", \"22\", \"27\", \"45\", \"57\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8dda0875",
   "metadata": {},
   "outputs": [],
   "source": [
    "# save table with preprocessing parameters per subject\n",
    "df = pd.DataFrame(params).T\n",
    "df.index.name = \"subj\"\n",
    "df = df.reset_index()\n",
    "df[\"subj\"] = pd.to_numeric(df.subj)\n",
    "df.sort_values(\"subj\").to_csv(\"preproc_params.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "148fc468",
   "metadata": {},
   "outputs": [],
   "source": [
    "# save table with notes\n",
    "pd.DataFrame(notes.items(), columns=[\"subj\", \"notes\"]).to_csv(\"preproc_notes.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc647b13",
   "metadata": {},
   "outputs": [],
   "source": [
    "# write final preprocessed files for further analysis\n",
    "preprocs = {subj: d for subj, d in preprocs.items() if subj not in exclude_preproc}\n",
    "pp.write_file(preprocs, \"data/export/eye_preproc.pd\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "868534a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# save table with a final summary of the preprocessed data\n",
    "def count_events(d, prefix):\n",
    "    return sum(lab.startswith(prefix) for lab in d.events.labels)\n",
    "\n",
    "tab = pd.DataFrame([{\"subj\": int(subj),\n",
    "                     \"duration_minutes\": d.get_duration(units=\"min\"),\n",
    "                     \"nevents\": d.nevents(),\n",
    "                     \"ncues\": count_events(d, \"C_\"),\n",
    "                     \"nfb\": count_events(d, \"F_\"),\n",
    "                     \"nresp\": count_events(d, \"R_\"),\n",
    "                     \"perc_interpolated\": 100 * d.data.get_layer(\"left_pupil\", \"interpolated\").mean(),\n",
    "                     \"blinks_per_min\": d.blink_stats()[\"left_pupil\"][\"n\"] / d.get_duration(units=\"min\"),\n",
    "                     \"note\": notes[subj]}\n",
    "                    for subj, d in preprocs.items()])\n",
    "tab.sort_values(\"subj\").to_csv(\"preproc_summary.csv\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python"
  },
  "nbsphinx": {
   "execute": "never"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}