{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Threshold data from Horsager et al. (2009)\n\nThis example shows how to use the Horsager et al. (2009) dataset.\n\n[Horsager2009]_ used a set of psychophysical detection tasks to determine\nperceptual thresholds in two Argus I users.\n\n.. important ::\n\n        You will need to install [Pandas](https://pandas.pydata.org)\n        (``pip install pandas``) for this dataset.\n\nA typical task involved turning on a single electrode in the array and asking\nthe subject if they saw something. If the subject did not detect a visual\npercept (i.e., a phosphene), the stimulus amplitude was increased.\nUsing a staircase procedure, the researchers were able to determine the current\namplitude at which subjects were able to detect a phosphene 50% of the time.\nThis amplitude is called the threshold current, and it is stored in the\n\"stim_amp\" column of the dataset.\n\n.. note ::\n\n    Mean threshold currents were extracted from the main publication and its\n    supplementary materials using Webplot Digitizer.\n    Therefore, these data are provided without any guarantee of correctness or\n    completeness.\n\n## Loading the dataset\n\nThe dataset can be loaded as a Pandas ``DataFrame``:\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from pulse2percept.datasets import load_horsager2009\ndata = load_horsager2009()\nprint(data)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Inspecting the DataFrame tells us that there are 552 threshold measurements\n(the rows) each with 21 different attributes (the columns).\n\nThese attributes include specifiers such as \"subject\", \"electrode\", and\n\"stim_type\". We can print all column names using:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "data.columns"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. note ::\n\n    The meaning of all column names is explained in the docstring of\n    the :py:func:`~pulse2percept.datasets.load_horsager2009` function.\n\nFor example, \"stim_type\" corresponds to the different stimulus types that\nwere used in the paper:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "data.stim_type.unique()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The first entry, \"single_pulse\", corresponds to the single biphasic pulse\nused to produce Figure 3. We can verify which figure a data point came from\nby inspecting the \"source\" column.\n\nTo select all the \"single_pulse\" rows, we can index into the DataFrame as\nfollows:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print(data[data.stim_type == 'single_pulse'])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "This leaves us with 80 rows, some of which come from Figure 3, others from\nFigure S3.1 in the supplementary material.\n\nAn alternative to indexing into the DataFrame is to load only a subset of\nthe data:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print(load_horsager2009(stim_types='single_pulse'))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Plotting the data\n\nArguably the most important column is \"stim_amp\". This is the current\namplitude of the different stimuli (single pulse, pulse trains, etc.) used\nat threshold.\n\nWe might be interested in seeing how threshold amplitude varies as a function\nof pulse duration. We could either use Matplotlib to generate a scatter plot\nor use pulse2percept's own visualization function:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from pulse2percept.viz import scatter_correlation\nscatter_correlation(data.pulse_dur, data.stim_amp)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ":py:func:`~pulse2percept.viz.scatter_correlation` above generates a scatter\nplot of the stimulus amplitude as a function of pulse duration, and performs\nlinear regression to calculate a correlation $r$ and a $p$ value.\nAs expected from the literature, now it becomes evident that stimulus\namplitude is negatively correlated with pulse duration (no matter the exact\nstimulus used).\n\n## Recreating the stimuli\n\nTo recreate the stimulus used to obtain a specific data point, we need to use\nthe values specified in the different columns of a particular row.\n\nFor example, the first row of the dataset specifies the stimulus used to\nobtain the threshold current for subject S05 on Electrode C3 using a single\nbiphasic pulse:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "row = data.loc[0, :]\nprint(row)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Using :py:class:`~pulse2percept.stimuli.BiphasicPulse`, we can create the\npulse to specifications:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from pulse2percept.stimuli import BiphasicPulse\nstim = BiphasicPulse(row['stim_amp'], row['pulse_dur'],\n                     interphase_dur=row['interphase_dur'],\n                     stim_dur=row['stim_dur'], electrode=row['electrode'])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can zoom in on the first 10 ms to see the stimulus waveform:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nstim.plot(time=np.arange(0, 10, 0.01))"
      ]
    }
  ],
  "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.7.9"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}