{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Nanduri et al. (2012): Frequency vs. amplitude modulation\n\nThis example shows how to use the\n:py:class:`~pulse2percept.models.Nanduri2012Model`.\n\nThe model introduced in [Nanduri2012]_ assumes that electrical stimulation\nleads to percepts that quickly increase in brightness (over the time course\nof ~100ms) and then slowly fade away (over the time course of seconds).\nThe model also assumes that amplitude and frequency modulation have different\neffects on the perceived brightness and size of a phosphene.\n\n## Generating a pulse train\n\nThe first step is to build a pulse train using the\n:py:class:`~pulse2percept.stimuli.PulseTrain` class.\nWe want to generate a 20Hz pulse train (0.45ms pulse duration, cathodic-first)\nat 30uA that lasts for a second:\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from pulse2percept.stimuli import BiphasicPulseTrain, Stimulus\ntsample = 0.005  # sampling time step (ms)\nphase_dur = 0.45  # duration of the cathodic/anodic phase (ms)\nstim_dur = 1000  # stimulus duration (ms)\namp_th = 30  # threshold current (uA)\nstim = BiphasicPulseTrain(20, amp_th, phase_dur, interphase_dur=phase_dur,\n                          stim_dur=stim_dur)\n\n# Configure Matplotlib:\nimport matplotlib.pyplot as plt\nplt.style.use('ggplot')\nfrom matplotlib import rc\nrc('font', size=12)\n\n# Plot the stimulus in the range t=[0, 60] ms:\nstim.plot(time=(0, 60))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Creating an implant\n\nBefore we can run the Nanduri model, we need to create a retinal implant to\nwhich we can assign the above pulse train.\n\nFor the purpose of this exercise, we will create an\n:py:class:`~pulse2percept.implants.ElectrodeArray` consisting of a single\n:py:class:`~pulse2percept.implants.DiskElectrode` with radius=260um centered\nat (x,y) = (0,0); i.e., centered over the fovea:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from pulse2percept.implants import DiskElectrode, ElectrodeArray\nearray = ElectrodeArray(DiskElectrode(0, 0, 0, 260))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Usually we would use a predefined retinal implant such as\n:py:class:`~pulse2percept.implants.ArgusII` or\n:py:class:`~pulse2percept.implants.AlphaIMS`. Alternatively, we can wrap the\nelectrode array created above with a\n:py:class:`~pulse2percept.implants.ProsthesisSystem` to create our own\nretinal implant. We will also assign the above created stimulus to it:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from pulse2percept.implants import ProsthesisSystem\nimplant = ProsthesisSystem(earray, stim=stim)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Running the model\n\nInteracting with a model always involves three steps:\n\n1.  **Initalize** the model by passing the desired parameters.\n2.  **Build** the model to perform (sometimes expensive) one-time setup\n    computations.\n3.  **Predict** the percept for a given stimulus.\n\nIn the following, we will run the Nanduri model on a single pixel, (0, 0):\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from pulse2percept.models import Nanduri2012Model\nmodel = Nanduri2012Model(xrange=(0, 0), yrange=(0, 0))\nmodel.build()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "<div class=\"alert alert-info\"><h4>Note</h4><p>You can also directly instantiate\n    :py:class:`~pulse2percept.models.Nanduri2012Temporal` and pass a stimulus\n    to it. However, please note that there might be subtle differences; e.g.,\n    :py:class:`~pulse2percept.models.Nanduri2012Model` will pass the stimulus\n    through the spatial model first.</p></div>\n\nAfter building the model, we are ready to predict the percept.\n\nThe input to the model is an implant containing a stimulus (usually a pulse\ntrain), which is processed in time by a number of linear filtering steps as\nwell as a stationary nonlinearity (a sigmoid).\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\npercept = model.predict_percept(implant, t_percept=np.arange(1000))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The output of the model is a :py:class:`~pulse2percept.percepts.Percept`\nobject that contains a time series with the predicted brightness of the\nvisual percept at every time step.\n\n\"Brightness of a stimulus\" is defined as the maximum brightness value\nencountered over time:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "bright_th = percept.data.max()\nbright_th"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plotting the percept next to the applied stimulus reveals that the model\npredicts the perceived brightness to increase rapidly (within ~100ms) and\nthen drop off slowly (over the time course of seconds).\nThis is consistent with behavioral reports from Argus II users.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, ax = plt.subplots(figsize=(12, 5))\nax.plot(implant.stim.time,\n        -0.02 + 0.01 * implant.stim.data[0, :] / implant.stim.data.max(),\n        linewidth=3, label='pulse')\nax.plot(percept.time, percept.data[0, 0, :], linewidth=3, label='percept')\nax.plot([0, stim_dur], [bright_th, bright_th], 'k--', label='max brightness')\nax.plot([0, stim_dur], [0, 0], 'k')\n\nax.set_xlabel('time (s)')\nax.set_ylabel('predicted brightness (a.u.)')\nax.set_yticks(np.arange(0, 0.14, 0.02))\nax.set_xlim(0, stim_dur)\nfig.legend(loc='center')\nfig.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "<div class=\"alert alert-info\"><h4>Note</h4><p>In psychophysical experiments, brightness is usually expressed in\n    relative terms; i.e., compared to a reference stimulus. In the model,\n    brightness has arbitrary units.</p></div>\n\n## Brightness as a function of frequency/amplitude\n\nThe Nanduri paper reports that amplitude and frequency have different effects\non the perceived brightness of a phosphene.\n\nTo study these effects, we will apply the model to a number of amplitudes and\nfrequencies:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Use the following pulse duration (ms):\npdur = 0.45\n# Generate values in the range [0, 50] uA with a step size of 5 uA or smaller:\namps = np.linspace(0, 50, 11)\n# Initialize an empty list that will contain the predicted brightness values:\nbright_amp = []\nfor amp in amps:\n    # For each value in the `amps` vector, now stored as `amp`, do the\n    # following:\n    # 1. Generate a pulse train with amplitude `amp`, 20 Hz frequency, 0.5 s\n    #    duration, pulse duration `pdur`, and interphase gap `pdur`:\n    implant.stim = BiphasicPulseTrain(20, amp, pdur, interphase_dur=pdur,\n                                      stim_dur=stim_dur)\n    # 2. Run the temporal model:\n    percept = model.predict_percept(implant)\n    # 3. Find the largest value in percept, this will be the predicted\n    # brightness:\n    bright_pred = percept.data.max()\n    # 4. Append this value to `bright_amp`:\n    bright_amp.append(bright_pred)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We then repeat the procedure for a whole range of frequencies:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Generate values in the range [0, 100] Hz with a step size of 10 Hz or\n# smaller:\nfreqs = np.linspace(0, 100, 11)\n# Initialize an empty list that will contain the predicted brightness values:\nbright_freq = []\nfor freq in freqs:\n    # For each value in the `amps` vector, now stored as `amp`, do the\n    # following:\n    # 1. Generate a pulse train with amplitude `amp`, 20 Hz frequency, 0.5 s\n    #    duration, pulse duration `pdur`, and interphase gap `pdur`:\n    implant.stim = BiphasicPulseTrain(freq, 20, pdur, interphase_dur=pdur,\n                                      stim_dur=stim_dur)\n    # 2. Run the temporal model\n    percept = model.predict_percept(implant)\n    # 3. Find the largest value in percept, this will be the predicted\n    # brightness:\n    bright_pred = percept.data.max()\n    # 4. Append this value to `bright_amp`:\n    bright_freq.append(bright_pred)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plotting the two curves side-by-side reveals that the model predicts\nbrightness to saturate quickly with increasing amplitude, but to scale\nlinearly with stimulus frequency:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, ax = plt.subplots(ncols=2, sharey=True, figsize=(12, 5))\n\nax[0].plot(amps, bright_amp, 'o-', linewidth=4)\nax[0].set_xlabel('amplitude (uA)')\nax[0].set_ylabel('predicted brightness (a.u.)')\n\nax[1].plot(freqs, bright_freq, 'o-', linewidth=4)\nax[1].set_xlabel('frequency (Hz)')\n\nfig.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Phosphene size as a function of amplitude/frequency\n\nThe paper also reports that phosphene size is affected differently by\namplitude vs. frequency modulation.\n\nTo introduce space into the model, we need to re-instantiate the model and\nthis time provide a range of (x,y) values to simulate. These values are\nspecified in degrees of visual angle (dva). They are sampled at ``xydva``\ndva:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "model = Nanduri2012Model(xystep=0.5, xrange=(-4, 4), yrange=(-4, 4))\nmodel.build()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We will again apply the model to a whole range of amplitude and frequency\nvalues taken directly from the paper:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Use the amplitude values from the paper:\namp_factors = [1, 1.25, 1.5, 2, 4, 6]\n\n# Output brightness in 1ms time steps:\nt_percept = np.arange(0, stim_dur, 1)\n\n# Initialize an empty list that will contain the brightest frames:\nframes_amp = []\n\nfor amp_f in amp_factors:\n    # For each value in the `amp_factors` vector, now stored as `amp_f`, do:\n    # 1. Generate a pulse train with amplitude `amp_f` * `amp_th`, frequency\n    #    20Hz, 0.5s duration, pulse duration `pdur`, and interphase gap `pdur`:\n    implant.stim = BiphasicPulseTrain(20, amp_f * amp_th, pdur,\n                                      interphase_dur=pdur, stim_dur=stim_dur)\n    # 2. Run the temporal model:\n    percept = model.predict_percept(implant, t_percept=t_percept)\n    # 3. Save the brightest frame:\n    frames_amp.append(percept.max(axis='frames'))\n\n# Use the amplitude values from the paper:\nfreqs = [40.0 / 3, 20, 2.0 * 40 / 3, 40, 80, 120]\n\n# Initialize an empty list that will contain the brightest frames:\nframes_freq = []\n\nfor freq in freqs:\n    # For each value in the `freqs` vector, now stored as `freq`, do:\n    # 1. Generate a pulse train with amplitude 1.25 * `amp_th`, frequency\n    #    `freq`, 0.5s duration, pulse duration `pdur`, and interphase gap\n    #    `pdur`:\n    implant.stim = BiphasicPulseTrain(freq, 1.25 * amp_th, pdur,\n                                      interphase_dur=pdur, stim_dur=stim_dur)\n    # 2. Run the temporal model:\n    percept = model.predict_percept(implant, t_percept=t_percept)\n    # 3. Save the brightest frame:\n    frames_freq.append(percept.max(axis='frames'))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "This allows us to reproduce Fig. 7 of [Nanduri2012]_:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, axes = plt.subplots(nrows=2, ncols=len(amp_factors), figsize=(16, 6))\n\nfor ax, amp, frame in zip(axes[0], amp_factors, frames_amp):\n    ax.imshow(frame, vmin=0, vmax=0.3, cmap='gray')\n    ax.set_title(f'{amp:.2g} xTh / 20 Hz', fontsize=16)\n    ax.set_xticks([])\n    ax.set_yticks([])\naxes[0][0].set_ylabel('amplitude\\nmodulation')\n\nfor ax, freq, frame in zip(axes[1], freqs, frames_freq):\n    ax.imshow(frame, vmin=0, vmax=0.3, cmap='gray')\n    ax.set_title(f'1.25xTh / {freq} Hz', fontsize=16)\n    ax.set_xticks([])\n    ax.set_yticks([])\naxes[1][0].set_ylabel('frequency\\nmodulation')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Phosphene size as a function of brightness\n\nLastly, the above data can also be visualized as a function of brightness to\nhighlight the difference between frequency and amplitude modulation (see\nFig.8 of [Nanduri2012]_):\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt.plot([np.max(frame) for frame in frames_amp],\n         [np.sum(frame >= bright_th) for frame in frames_amp],\n         'o-', label='amplitude modulation')\nplt.plot([np.max(frame) for frame in frames_freq],\n         [np.sum(frame >= bright_th) for frame in frames_freq],\n         'o-', label='frequency modulation')\nplt.xlabel('brightness (a.u.)')\nplt.ylabel('area (# pixels)')\nplt.legend()"
      ]
    }
  ],
  "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
}