{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Beyeler et al. (2019): Axonal streaks with the axon map model\n\nThis example shows how to apply the\n:py:class:`~pulse2percept.models.AxonMapModel` to an\n:py:class:`~pulse2percept.implants.ArgusII` implant.\n\nThe axon map model assumes that electrical stimulation leads to percepts that\nare elongated along the direction of the underlying nerve fiber bundle\ntrajectory. Because the layout of nerve fiber bundles in the human retina is\nhighly stereotyped [Jansonius2009]_, percept shape is predictable based on\n(but also highly variable depending on) the location of the stimulating\nelectrode.\n\nAn axon's sensitivity to electrical stimulation is assumed to decay\nexponentially:\n\n*  with distance from the soma $(x_{soma}, y_{soma})$, with spatial decay\n   constant $\\lambda$,\n*  with distance from the stimulated retinal location\n   $(x_{stim}, y_{stim})$, with spatial decay constant $\\rho$:\n\n\\begin{align}I_{axon}(x,y; \\rho, \\lambda) =& \\exp \\Big(\n    -\\frac{(x-x_{stim})^2 + (y-y_{stim})^2}{2 \\rho^2} \\Big) \\\\\n                                    & \\exp \\Big(\n    -\\frac{(x-x_{soma})^2 + (y-y_{soma})^2}{2 \\lambda^2} \\Big).\\end{align}\n\nThe axon map model can be instantiated and run in three steps.\n\n## Creating the model\n\nThe first step is to instantiate the\n:py:class:`~pulse2percept.models.AxonMapModel` class by calling its\nconstructor method.\nThe two most important parameters to set are ``rho`` and ``axlambda`` from\nthe equation above (here set to 150 micrometers and 500 micrometers,\nrespectively):\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nfrom pulse2percept.implants import ArgusII\nfrom pulse2percept.models import AxonMapModel\nmodel = AxonMapModel(rho=150, axlambda=500)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Parameters you don't specify will take on default values. You can inspect\nall current model parameters as follows:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print(model)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "This reveals a number of other parameters to set, such as:\n\n* ``xrange``, ``yrange``: the extent of the visual field to be simulated,\n  specified as a range of x and y coordinates (in degrees of visual angle,\n  or dva). For example, we are currently sampling x values between -20 dva\n  and +20dva, and y values between -15 dva and +15 dva.\n* ``xystep``: The resolution (in dva) at which to sample the visual field.\n  For example, we are currently sampling at 0.25 dva in both x and y\n  direction.\n* ``loc_od_x``, ``loc_od_y``: the location of the center of the optic disc\n  (in dva)\n* ``thresh_percept``: You can also define a brightness threshold, below which\n  the predicted output brightness will be zero. It is currently set to\n  ``1/sqrt(e)``, because that will make the radius of the predicted percept\n  equal to ``rho``.\n\nA number of parameters control the amount of detail used when generating the\naxon map:\n\n* ``n_axons``: the number of axons to generate\n* ``axons_range``: the range of angles (in degrees) to use at which axon\n  trajectories emanate from the center of the optic disc\n* ``n_ax_segments``: the number of segments each generated axon should have\n* ``n_ax_segments_range``: the range of distances (in dva) to use, measured\n  from the center of the optic disc, at which axon segments should be placed\n* ``axons_pickle``: path to a pickle file where previously generated axon\n  maps are stored\n\nIn addition, you can choose the parallelization back end used to speed up\nsimulations:\n\n* ``engine``:\n   * 'serial': single-core processing (no parallelization)\n   * 'joblib': parallelization using the `JobLib`_ library\n   * 'dask': parallelization using the `Dask`_ library\n\n* ``scheduler``:\n   * 'threading': a scheduler backed by a thread pool\n   * 'multiprocessing': a scheduler backed by a process pool\n\n\nTo change parameter values, either pass them directly to the constructor\nabove or set them by hand, like this:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "model.engine = 'serial'"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Then build the model. This is a necessary step before you can actually use\nthe model to predict a percept, as it performs a number of expensive setup\ncomputations (e.g., building the axon map, calculating electric potentials):\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "model.build()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. important ::\n\n    You need to build a model only once. After that, you can apply any number\n    of stimuli -- or even apply the model to different implants -- without\n    having to rebuild (which takes time).\n\n    However, if you change important model parameters outside the constructor\n    (e.g., by directly setting ``model.axlambda = 100``), you will have to\n    call ``model.build()`` again for your changes to take effect.\n\n## Assigning a stimulus\nThe second step is to specify a visual prosthesis from the\n:py:mod:`~pulse2percept.implants` module.\n\nIn the following, we will create an\n:py:class:`~pulse2percept.implants.ArgusII` implant. By default, the implant\nwill be centered over the fovea (at x=0, y=0) and aligned with the horizontal\nmeridian (rot=0):\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "implant = ArgusII()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "You can inspect the location of the implant with respect to the underlying\nnerve fiber bundles using the built-in plot methods:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "model.plot()\nimplant.plot()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "By default, the plots will be added to the current Axes object.\nAlternatively, you can pass ``ax=`` to specify in which Axes to plot.\n\nThe easiest way to assign a stimulus to the implant is to pass a NumPy array\nthat specifies the current amplitude to be applied to every electrode in the\nimplant.\n\nFor example, the following sends 1 microamp to all 60 electrodes of the\nimplant:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "implant.stim = np.ones(60)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Predicting the percept\nThe third step is to apply the model to predict the percept resulting from\nthe specified stimulus. Note that this may take some time on your machine:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "percept = model.predict_percept(implant)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The resulting percept is stored in a\n:py:class:`~pulse2percept.percepts.Percept` object, which is similar in\norganization to the :py:class:`~pulse2percept.stimuli.Stimulus` object:\nthe ``data`` container is a 3D NumPy array (Y, X, T) with labeled axes\n``xdva``, ``ydva``, and ``time``.\n\nThe percept can be plotted as follows:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ax = percept.plot()\nax.set_title('Predicted percept')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "A major prediction of the axon map model is that the percept changes\ndepending on the location of the implant. You can convince yourself of that\nby re-running the model on an implant shifted and rotated across the retina:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "implant = ArgusII(x=-50, y=50, rot=-45)\nmodel.plot()\nimplant.plot()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The resulting percepts should look very different from the previous example:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "implant.stim = np.ones(60)\npercept = model.predict_percept(implant)\nax = percept.plot()\nax.set_title('Predicted percept')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. important::\n\n    When specifying the rotation of the implant, positive angles will result\n    in counterclockwise rotations **on the retinal surface**.\n\n    However, because the superior (inferior) retina is mapped onto the lower\n    (upper) visual field, a counterclockwise orientation on the retina is\n    equivalent to a clockwise orientation of the percept in visual field\n    coordinates.\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "You can also use the axon map model to imitate\n:py:class:`~pulse2percept.models.ScoreboardModel` by setting lambda to a small\nvalue.\nHowever, you may have to increase the number of axons and number of segments\nper axon to get a smooth percept out:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "model = AxonMapModel(rho=200, axlambda=10, n_axons=3000, n_ax_segments=3000)\nmodel.build()\npercept = model.predict_percept(implant)\nax = percept.plot()\nax.set_title('Predicted percept')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "This is of course not very computationally efficient, because the model is\nstill performing all the axon map calculations.\nIn this case, you might be better off using\n:py:class:`~pulse2percept.models.ScoreboardModel`.\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
}