{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Phosphene drawings from Beyeler et al. (2019)\n\nThis example shows how to use the Beyeler et al. (2019) dataset.\n\n[Beyeler2019]_ asked Argus I/II users to draw what they see in response to\nsingle-electrode stimulation.\n\n.. important ::\n\n    For this dataset you will need to install both\n    [Pandas](https://pandas.pydata.org) (``pip install pandas``) and\n    [HDF4 for Python](https://www.h5py.org) (``pip install h5py``).\n\n## Loading the dataset\n\nDue to its size (66 MB), the dataset is not included with pulse2percept, but\ncan be downloaded from the Open Science Framework (OSF).\n\nBy default, the dataset will be stored in a local directory\n\u2018~/pulse2percept_data/\u2019 within your user directory (but a different path can be\nspecified).\nThis way, the datasets is only downloaded once, and future calls to the fetch\nfunction will load the dataset from your local copy.\n\nThe data itself will be provided as a Pandas ``DataFrame``:\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from pulse2percept.datasets import fetch_beyeler2019\n\ndata = fetch_beyeler2019()\nprint(data)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Inspecting the DataFrame tells us that there are 400 phosphene drawings\n(the rows) each with 16 different attributes (the columns).\n\nThese attributes include specifiers such as \"subject\", \"electrode\", and\n\"image\". 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.fetch_beyeler2019` function.\n\nFor example, \"subject\" contains the different subject IDs used in the study:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "data.subject.unique()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "To select all drawings from Subject 2, we can index into the DataFrame as\nfollows:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print(data[data.subject == 'S2'])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "This leaves us with 110 rows, each of which correspond to one phosphene\ndrawings from a number of different electrodes and trials.\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(fetch_beyeler2019(subjects='S2'))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Plotting the data\n\nArguably the most important column is \"image\". This is the phosphene drawing\nobtained during a particular trial.\n\nEach phosphene drawing is a 2D black-and-white NumPy array, so we can just\nplot it using Matplotlib like any other image:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\nplt.imshow(data.loc[0, 'image'], cmap='gray')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "However, we might be more interested in seeing how phosphene shape differs\nfor different electrodes.\nFor this we can use :py:func:`~pulse2percept.viz.plot_argus_phosphenes` from\nthe :py:mod:`~pulse2percept.viz` module.\nIn addition to the ``data`` matrix, the function will also want an\n:py:class:`~pulse2percept.implants.ArgusII` object implanted at the correct\nlocation.\n\nConsulting [Beyeler2019]_ tells us that the prosthesis was roughly implanted\nin the following location:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from pulse2percept.implants import ArgusII\nargus = ArgusII(x=-1331, y=-850, rot=-28.4, eye='RE')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "For now, let's focus on the data from Subject 2:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "data = fetch_beyeler2019(subjects='S2')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Passing both ``data`` and ``argus`` to\n:py:func:`~pulse2percept.viz.plot_argus_phosphenes` will then allow the\nfunction to overlay the phosphene drawings over a schematic of the implant.\nHere, phosphene drawings from different trials are averaged, and aligned with\nthe center of the electrode that was used to obtain the drawing:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from pulse2percept.viz import plot_argus_phosphenes\nplot_argus_phosphenes(data, argus)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Great! We have just reproduced a panel from Figure 2 in [Beyeler2019]_.\n\nAs [Beyeler2019]_ went on to show, the orientation of these phosphenes is\nwell aligned with the map of nerve fiber bundles (NFBs) in each subject's\neye.\n\nTo see how the phosphene drawings line up with the NFBs, we can also pass an\n:py:class:`~pulse2percept.models.AxonMapModel` to the function.\nOf course, we need to make sure that we use the correct dimensions. Subject\nS2 had their optic disc center located 16.2 deg nasally, 1.38 deg superior\nfrom the fovea:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from pulse2percept.models import AxonMapModel\nmodel = AxonMapModel(loc_od=(16.2, 1.38))\nplot_argus_phosphenes(data, argus, axon_map=model)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Predicting phosphene shape\n\nIn addition, the :py:class:`~pulse2percept.models.AxonMapModel` is well\nsuited to predict the shape of individual phosphenes. Using the values given\nin [Beyeler2019]_, we can tailor the axon map parameters to Subject 2:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nmodel = AxonMapModel(rho=315, axlambda=500, loc_od=(16.2, 1.38),\n                     xrange=(-30, 30), yrange=(-22.5, 22.5),\n                     thresh_percept=1 / np.sqrt(np.e))\nmodel.build()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now we need to activate one electrode at a time, and predict the resulting\npercept. We could build a :py:class:`~pulse2percept.stimuli.Stimulus` object\nwith a for loop that does just that, or we can use the following trick.\n\nThe stimulus' data container is a (electrodes, timepoints) shaped 2D NumPy\narray. Activating one electrode at a time is therefore the same as an\nidentity matrix whose size is equal to the number of electrodes. In code:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Find the names of all the electrodes in the dataset:\nelectrodes = data.electrode.unique()\n# Activate one electrode at a time:\nimport numpy as np\nfrom pulse2percept.stimuli import Stimulus\nargus.stim = Stimulus(np.eye(len(electrodes)), electrodes=electrodes)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Using the model's\n:py:func:`~pulse2percept.models.AxonMapModel.predict_percept`, we then get\na Percept object where each frame is the percept generated from activating\na single electrode:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "percepts = model.predict_percept(argus)\npercepts.play()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Finally, we can visualize the ground-truth and simulated phosphenes\nside-by-side:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from pulse2percept.viz import plot_argus_simulated_phosphenes\nfig, (ax_data, ax_sim) = plt.subplots(ncols=2, figsize=(15, 5))\nplot_argus_phosphenes(data, argus, scale=0.75, ax=ax_data)\nplot_argus_simulated_phosphenes(percepts, argus, scale=1.25, ax=ax_sim)\nax_data.set_title('Ground-truth phosphenes')\nax_sim.set_title('Simulated phosphenes')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Analyzing phosphene shape\n\nThe phosphene drawings also come annotated with different shape descriptors:\narea, orientation, and elongation.\nElongation is also called eccentricity in the computer vision literature,\nwhich is not to be confused with retinal eccentricity. It is simply a number\nbetween 0 and 1, where 0 corresponds to a circle and 1 corresponds to an\ninfinitesimally thin line (note that the Methods section of [Beyeler2019]_\ngot it wrong).\n\n[Beyeler2019]_ made the point that if each phosphene could be considered a\npixel (or essentially a blob), as is so often assumed in the literature, then\nmost phosphenes should have zero elongation.\n\nInstead, using Matplotlib's histogram function, we can convince ourselves\nthat most phosphenes are in fact elongated:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "data = fetch_beyeler2019()\ndata.eccentricity.plot(kind='hist')\nplt.xlabel('phosphene elongation')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Phosphenes are not pixels!\nAnd with that we have just reproduced Fig. 3C of [Beyeler2019]_.\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.7.9"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}