
.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "examples\models\plot_beyeler2019_axonmap.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_examples_models_plot_beyeler2019_axonmap.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_examples_models_plot_beyeler2019_axonmap.py:


============================================================================
Beyeler et al. (2019): Axonal streaks with the axon map model
============================================================================

This example shows how to apply the
:py:class:`~pulse2percept.models.AxonMapModel` to an
:py:class:`~pulse2percept.implants.ArgusII` implant.

The axon map model assumes that electrical stimulation leads to percepts that
are elongated along the direction of the underlying nerve fiber bundle
trajectory. Because the layout of nerve fiber bundles in the human retina is
highly stereotyped [Jansonius2009]_, percept shape is predictable based on
(but also highly variable depending on) the location of the stimulating
electrode.

An axon's sensitivity to electrical stimulation is assumed to decay
exponentially:

*  with distance from the soma :math:`(x_{soma}, y_{soma})`, with spatial decay
   constant :math:`\lambda`,
*  with distance from the stimulated retinal location
   :math:`(x_{stim}, y_{stim})`, with spatial decay constant :math:`\rho`:

.. math::

    I_{axon}(x,y; \rho, \lambda) =& \exp \Big(
    -\frac{(x-x_{stim})^2 + (y-y_{stim})^2}{2 \rho^2} \Big) \\
                                    & \exp \Big(
    -\frac{(x-x_{soma})^2 + (y-y_{soma})^2}{2 \lambda^2} \Big).

The axon map model can be instantiated and run in three steps.

Creating the model
------------------

The first step is to instantiate the
:py:class:`~pulse2percept.models.AxonMapModel` class by calling its
constructor method.
The two most important parameters to set are ``rho`` and ``lam`` from
the equation above (here set to 150 micrometers and 500 micrometers,
respectively):

.. GENERATED FROM PYTHON SOURCE LINES 46-52

.. code-block:: Python


    import numpy as np
    from pulse2percept.implants import ArgusII
    from pulse2percept.models import AxonMapModel
    model = AxonMapModel(rho=150, lam=500)








.. GENERATED FROM PYTHON SOURCE LINES 54-56

Parameters you don't specify will take on default values. You can inspect
all current model parameters as follows:

.. GENERATED FROM PYTHON SOURCE LINES 56-59

.. code-block:: Python


    print(model)





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    AxonMapModel(ax_segments_range=(0, 50), 
                 axon_pickle='axons.pickle', 
                 axons_range=(-180, 180), eye='RE', 
                 grid_type='rectangular', ignore_pickle=False, 
                 lam=500, loc_od=(15.5, 1.5), 
                 min_ax_sensitivity=0.001, 
                 min_current_spread=1e-08, n_ax_segments=500, 
                 n_axons=1000, n_gray=None, n_jobs=14, 
                 n_threads=14, ndim=[2], noise=None, rho=150, 
                 spatial=AxonMapSpatial, step=0.25, 
                 temporal=None, thresh_percept=0, verbose=True, 
                 vfmap=Watson2014Map(ndim=2), xrange=(-15, 15), 
                 yrange=(-15, 15))




.. GENERATED FROM PYTHON SOURCE LINES 60-94

This reveals a number of other parameters to set, such as:

* ``xrange``, ``yrange``: the extent of the visual field to be simulated,
  specified as a range of x and y coordinates (in degrees of visual angle,
  or dva). For example, we are currently sampling x values between -20 dva
  and +20dva, and y values between -15 dva and +15 dva.
* ``step``: The resolution (in dva) at which to sample the visual field.
  For example, we are currently sampling at 0.25 dva in both x and y
  direction.
* ``loc_od_x``, ``loc_od_y``: the location of the center of the optic disc
  (in dva)
* ``thresh_percept``: You can also define a brightness threshold, below which
  the predicted output brightness will be zero. It is currently set to
  ``1/sqrt(e)``, because that will make the radius of the predicted percept
  equal to ``rho``.

A number of parameters control the amount of detail used when generating the
axon map:

* ``n_axons``: the number of axons to generate
* ``axons_range``: the range of angles (in degrees) to use at which axon
  trajectories emanate from the center of the optic disc
* ``n_ax_segments``: the number of segments each generated axon should have
* ``n_ax_segments_range``: the range of distances (in dva) to use, measured
  from the center of the optic disc, at which axon segments should be placed
* ``axons_pickle``: path to a pickle file where previously generated axon
  maps are stored

To change parameter values, either pass them directly to the constructor
above or set them by hand.

Then build the model. This is a necessary step before you can actually use
the model to predict a percept, as it performs a number of expensive setup
computations (e.g., building the axon map, calculating electric potentials):

.. GENERATED FROM PYTHON SOURCE LINES 94-97

.. code-block:: Python


    model.build()





.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    AxonMapModel(ax_segments_range=(0, 50), 
                 axon_pickle='axons.pickle', 
                 axons_range=(-180, 180), eye='RE', 
                 grid_type='rectangular', ignore_pickle=False, 
                 lam=500, loc_od=(np.float64(15.5), 1.5), 
                 min_ax_sensitivity=0.001, 
                 min_current_spread=1e-08, n_ax_segments=500, 
                 n_axons=1000, n_gray=None, n_jobs=14, 
                 n_threads=14, ndim=[2], noise=None, rho=150, 
                 spatial=AxonMapSpatial, step=0.25, 
                 temporal=None, thresh_percept=0, verbose=True, 
                 vfmap=Watson2014Map(ndim=2), xrange=(-15, 15), 
                 yrange=(-15, 15))



.. GENERATED FROM PYTHON SOURCE LINES 98-117

.. important ::

    You need to build a model only once. After that, you can apply any number
    of stimuli -- or even apply the model to different implants -- without
    having to rebuild (which takes time).

    However, if you change important model parameters outside the constructor
    (e.g., by directly setting ``model.lam = 100``), you will have to
    call ``model.build()`` again for your changes to take effect.

Assigning a stimulus
--------------------
The second step is to specify a visual prosthesis from the
:py:mod:`~pulse2percept.implants` module.

In the following, we will create an
:py:class:`~pulse2percept.implants.ArgusII` implant. By default, the implant
will be centered over the fovea (at x=0, y=0) and aligned with the horizontal
meridian (rot=0):

.. GENERATED FROM PYTHON SOURCE LINES 117-120

.. code-block:: Python


    implant = ArgusII()








.. GENERATED FROM PYTHON SOURCE LINES 121-123

You can inspect the location of the implant with respect to the underlying
nerve fiber bundles using the built-in plot methods:

.. GENERATED FROM PYTHON SOURCE LINES 123-128

.. code-block:: Python


    model.plot()
    implant.plot()





.. image-sg:: /examples/models/images/sphx_glr_plot_beyeler2019_axonmap_001.png
   :alt: plot beyeler2019 axonmap
   :srcset: /examples/models/images/sphx_glr_plot_beyeler2019_axonmap_001.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    <Axes: xlabel='x (microns)', ylabel='y (microns)'>



.. GENERATED FROM PYTHON SOURCE LINES 129-138

By default, the plots will be added to the current Axes object.
Alternatively, you can pass ``ax=`` to specify in which Axes to plot.

The easiest way to assign a stimulus to the implant is to pass a NumPy array
that specifies the current amplitude to be applied to every electrode in the
implant.

For example, the following sends 1 microamp to all 60 electrodes of the
implant:

.. GENERATED FROM PYTHON SOURCE LINES 138-141

.. code-block:: Python


    implant.stim = np.ones(60)








.. GENERATED FROM PYTHON SOURCE LINES 142-146

Predicting the percept
----------------------
The third step is to apply the model to predict the percept resulting from
the specified stimulus. Note that this may take some time on your machine:

.. GENERATED FROM PYTHON SOURCE LINES 146-149

.. code-block:: Python


    percept = model.predict_percept(implant)








.. GENERATED FROM PYTHON SOURCE LINES 150-157

The resulting percept is stored in a
:py:class:`~pulse2percept.percepts.Percept` object, which is similar in
organization to the :py:class:`~pulse2percept.stimuli.Stimulus` object:
the ``data`` container is a 3D NumPy array (Y, X, T) with labeled axes
``xdva``, ``ydva``, and ``time``.

The percept can be plotted as follows:

.. GENERATED FROM PYTHON SOURCE LINES 157-161

.. code-block:: Python


    ax = percept.plot()
    ax.set_title('Predicted percept')




.. image-sg:: /examples/models/images/sphx_glr_plot_beyeler2019_axonmap_002.png
   :alt: Predicted percept
   :srcset: /examples/models/images/sphx_glr_plot_beyeler2019_axonmap_002.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    Text(0.5, 1.0, 'Predicted percept')



.. GENERATED FROM PYTHON SOURCE LINES 162-165

A major prediction of the axon map model is that the percept changes
depending on the location of the implant. You can convince yourself of that
by re-running the model on an implant shifted and rotated across the retina:

.. GENERATED FROM PYTHON SOURCE LINES 165-170

.. code-block:: Python


    implant = ArgusII(x=-50, y=50, rot=-45)
    model.plot()
    implant.plot()




.. image-sg:: /examples/models/images/sphx_glr_plot_beyeler2019_axonmap_003.png
   :alt: plot beyeler2019 axonmap
   :srcset: /examples/models/images/sphx_glr_plot_beyeler2019_axonmap_003.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    <Axes: xlabel='x (microns)', ylabel='y (microns)'>



.. GENERATED FROM PYTHON SOURCE LINES 171-172

The resulting percepts should look very different from the previous example:

.. GENERATED FROM PYTHON SOURCE LINES 172-178

.. code-block:: Python


    implant.stim = np.ones(60)
    percept = model.predict_percept(implant)
    ax = percept.plot()
    ax.set_title('Predicted percept')




.. image-sg:: /examples/models/images/sphx_glr_plot_beyeler2019_axonmap_004.png
   :alt: Predicted percept
   :srcset: /examples/models/images/sphx_glr_plot_beyeler2019_axonmap_004.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    Text(0.5, 1.0, 'Predicted percept')



.. GENERATED FROM PYTHON SOURCE LINES 179-188

.. important::

    When specifying the rotation of the implant, positive angles will result
    in counterclockwise rotations **on the retinal surface**.

    However, because the superior (inferior) retina is mapped onto the lower
    (upper) visual field, a counterclockwise orientation on the retina is
    equivalent to a clockwise orientation of the percept in visual field
    coordinates.

.. GENERATED FROM PYTHON SOURCE LINES 190-195

You can also use the axon map model to imitate
:py:class:`~pulse2percept.models.ScoreboardModel` by setting lambda to a small
value.
However, you may have to increase the number of axons and number of segments
per axon to get a smooth percept out:

.. GENERATED FROM PYTHON SOURCE LINES 195-202

.. code-block:: Python


    model = AxonMapModel(rho=200, lam=10, n_axons=3000, n_ax_segments=3000)
    model.build()
    percept = model.predict_percept(implant)
    ax = percept.plot()
    ax.set_title('Predicted percept')




.. image-sg:: /examples/models/images/sphx_glr_plot_beyeler2019_axonmap_005.png
   :alt: Predicted percept
   :srcset: /examples/models/images/sphx_glr_plot_beyeler2019_axonmap_005.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    Text(0.5, 1.0, 'Predicted percept')



.. GENERATED FROM PYTHON SOURCE LINES 203-207

This is of course not very computationally efficient, because the model is
still performing all the axon map calculations.
In this case, you might be better off using
:py:class:`~pulse2percept.models.ScoreboardModel`.


.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 6.609 seconds)


.. _sphx_glr_download_examples_models_plot_beyeler2019_axonmap.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: plot_beyeler2019_axonmap.ipynb <plot_beyeler2019_axonmap.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_beyeler2019_axonmap.py <plot_beyeler2019_axonmap.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: plot_beyeler2019_axonmap.zip <plot_beyeler2019_axonmap.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_
