{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Generating simple pulses and pulse trains\n\nThis example shows how to build and visualize basic types of stimuli such as\n:py:class:`~pulse2percept.stimuli.MonophasicPulse`,\n:py:class:`~pulse2percept.stimuli.BiphasicPulse` or a\n:py:class:`~pulse2percept.stimuli.PulseTrain` for a given implant.\n\nA monophasic pulse has a single phase and can be either anodic (by definition:\nhas a positive current amplitude) or cathodic (negative current amplitude).\n\nA biphasic pulse is generally charge-balanced for safety reasons (i.e., the\nnet current must sum to zero over time) and defined as either anodic-first\nor cathodic-first.\n\nMultiple pulses can form a pulse train.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Simplest stimulus\n---------------------\n:py:class:`~pulse2percept.stimuli.Stimulus` is the base class to generate\ndifferent types of stimuli. The simplest way to instantiate a Stimulus is\nto pass a scalar value which is interpreted as the current amplitude\nfor a single electrode.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Let's start by importing necessary modules\nfrom pulse2percept.stimuli import (MonophasicPulse, BiphasicPulse,\n                                   Stimulus, PulseTrain)\n\nimport numpy as np\n\nstim = Stimulus(10)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Parameters we don't specify will take on default values. We can inspect\nall current model parameters as follows:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print(stim)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "This command also reveals a number of other parameters to set, such as:\n\n* ``electrodes``: We can either specify the electrodes in the source\n  or within the stimulus. If none are specified it looks up the source\n  electrode.\n\n* ``metadata``: Optionally we can include metadata to the stimulus we\n  generate as a dictionary.\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": [
        "stim.metadata = {'name': 'A simple stimulus', 'date': '2020-01-01'}\nstim"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "A monophasic pulse\n--------------------\nWe can specify the arguments of the monophasic pulse as follows:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "pulse_type = 'anodic'  # anodic: positive amplitude, cathodic: negative\npulse_dur = 4.6 / 1000  # pulse duration in seconds\ndelay_dur = 10.0 / 1000  # pulse delivered after delay in seconds\nstim_dur = 0.5  # stimulus duration in seconds (pulse padded with zeros)\ntime_step = 0.1 / 1000  # temporal sampling step in seconds"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The sampling step ``time_step`` defines at which temporal resolution the\nstimulus is resolved. In the above example, the time step is 0.1 ms.\n\nBy calling Stimulus with a ``MonophasicPulse`` source, we can generate a\nsingle pulse:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "monophasic_stim = Stimulus(MonophasicPulse(ptype=pulse_type, pdur=pulse_dur,\n                                           delay_dur=delay_dur,\n                                           stim_dur=stim_dur,\n                                           tsample=time_step))\nprint(monophasic_stim)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Here, ``data`` is a 2D NumPy array where rows are electrodes and columns are\nthe points in time. Since we did not specify any electrode names in\n``MonophasicPulse``, the number of electrodes is inferred from the input\nsource type. There is only one row in the above example, denoting a single\nelectrode.\n\nBy default, the :py:class:`~pulse2percept.stimuli.MonophasicPulse` object\nautomatically assumes a current amplitude of 1 uA.\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can visualize the generated pulse using Matplotlib:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\nfig, ax = plt.subplots(figsize=(8, 5))\nax.plot(monophasic_stim.time, monophasic_stim.data[0, :])\nax.set_xlabel('Time (s)')\nax.set_ylabel('Amplitude ($\\mu$A)')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "A biphasic pulse\n------------------\nSimilarly, we can generate a biphasic pulse by changing the source of the\nstimulus to :py:class:`~pulse2percept.stimuli.BiphasicPulse`. This time\nparameter ``ptype`` can either be 'anodicfirst' or 'cathodicfirst'.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# set relevant parameters\npulse_type = 'cathodicfirst'\nbiphasic_stim = Stimulus(BiphasicPulse(ptype=pulse_type, pdur=pulse_dur,\n                                       tsample=time_step))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "If we visualize this stimulus, we can see the difference between a monophasic\nand biphasic pulse:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Create a figure with two subplots\nfig, axes = plt.subplots(1, 2, figsize=(20, 10))\n\n# First, plot monophasic pulse\naxes[0].plot(monophasic_stim.time, monophasic_stim.data[0])\nax.set_xlabel('Time (s)')\nax.set_ylabel('Amplitude ($\\mu$A)')\n\n# Second, plot biphasic pulse\naxes[1].plot(biphasic_stim.time, biphasic_stim.data[0])\nax.set_xlabel('Time (s)')\nax.set_ylabel('Amplitude ($\\mu$A)')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Changing pulse amplitude\n----------------------------------\nFor any given pulse, we can modify the amplitude by indexing into the ``data``\nrow that corresponds to the desired electrode. In the above example, we only\nhave one electrode (index 0).\nLet's say we want the amplitude of the monophasic pulse to be 10 micro amps.\nWe have two options: either change the values of the ``data`` array directly:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# get the data structure by indexing the electrode at 0\nmonophasic_stim.data[0] = 10*monophasic_stim.data[0]\nprint(monophasic_stim)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Or we can create a NumPy array and assign that to the data structure of the\nstimulus:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# recreate the same stimulus with an amplitude 1 microAmps.\nmonophasic_stim = Stimulus(MonophasicPulse(ptype='anodic', pdur=pulse_dur,\n                                           delay_dur=delay_dur,\n                                           stim_dur=stim_dur,\n                                           tsample=time_step))\nmonophasic_stim.data[0] = 10*np.ones_like(monophasic_stim.data[0])\nprint(monophasic_stim)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Similarly, let's say we want the cathodic part of the biphasic pulse to be -5\nmicro amps, and the anodic part to be +20 micro amps (note that this stimulus\nwouldn't be charge-balanced).\n\nWe first need to find the halfway point where the current switches from\ncathodic to anodic. To do that we first get the length of the pulse by\nindexing the single electrode at 0\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "length = len(biphasic_stim.data[0])\nprint(length)\n\n# Find the halfway where cathodic turns into anodic pulse\nhalf = int(len(biphasic_stim.data[0])/2)\nprint(\"Halfway index is\", half)\n\n# change the first half of the pulse to be 5 times larger\nbiphasic_stim.data[0][0:half] = 5*biphasic_stim.data[0][0:half]\n\n# change the second half to be 20 times larger\nbiphasic_stim.data[0][half:length] = 20*biphasic_stim.data[0][half:length]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's plot the monophasic and biphasic pulses again:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Create a figure with two subplots\nfig, axes = plt.subplots(nrows=2, figsize=(25, 15))\n\n# First, plot monophasic pulse\naxes[0].plot(monophasic_stim.time, monophasic_stim.data[0])\naxes[0].set_xlabel('Time (s)')\naxes[0].set_ylabel('Amplitude ($\\mu$A)')\n# Second, plot biphasic pulse\naxes[1].plot(biphasic_stim.time, biphasic_stim.data[0])\naxes[1].set_xlabel('Time (s)')\naxes[1].set_ylabel('Amplitude ($\\mu$A)')\nfig.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Generating standard pulse trains\n----------------------------------\nThe easiest way to generate a pulse train is to use the \n:py:class:`~pulse2percept.stimuli.PulseTrain` object, which allows for\nvarious stimulus attributes to be specified:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "time_step = 0.1 / 1000  # temporal sampling in seconds\nfreq = 20  # frequency in Hz\namp = 100  # maximum amplitude of the pulse train in microAmps\ndur = 0.2  # total duration of the pulse train in seconds\npulse_type = 'cathodicfirst'  # whether the first phase is positive or negative\npulse_order = 'gapfirst'  # whether the train starts with gap or a pulse.\n\n# Define the pulse train with given parameters\nptrain = PulseTrain(tsample=time_step,\n                    freq=freq,\n                    dur=dur,\n                    amp=amp,\n                    pulsetype=pulse_type,\n                    pulseorder=pulse_order)\n\n# Create a new stimulus where the pulse train is the source\nptrain_stim = Stimulus(ptrain)\n\n# Visualize:\nfig, ax = plt.subplots(figsize=(8, 5))\nax.plot(ptrain_stim.time, ptrain_stim.data[0, :])\nax.set_xlabel('Time (s)')\nax.set_ylabel('Amplitude ($\\mu$A)')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Alternatively, we are free to specify a discrete set of points in time and\nthe current amplitude we would like to apply at those times.\n\nIt is important to note that the :py:class:`~pulse2percept.stimuli.Stimulus`\nobject will linearly interpolate between specified time points.\nFor example, the following generates a simple sawtooth stimulus:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "stim = Stimulus([[0, -10, 10, -10, 10, -10, 0]],\n                time=[0, 0.1, 0.3, 0.5, 0.7, 0.9, 1.0])\nfig, ax = plt.subplots(figsize=(8, 5))\nax.plot(stim.time, stim.data[0, :])\nax.set_xlabel('Time (s)')\nax.set_ylabel('Amplitude ($\\mu$A)')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "For a biphasic pulse, we need to specify both the rising edge (low-to-high\ntransition) and falling edge (high-to-low transition) of the signal:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "stim = Stimulus([[0, 0, 10, 10,  0, 0]],\n                time=[0, 0.1, 0.1, 0.2, 0.2, 1.0])\nfig, ax = plt.subplots(figsize=(8, 5))\nax.plot(stim.time, stim.data[0, :])\nax.set_xlabel('Time (s)')\nax.set_ylabel('Amplitude ($\\mu$A)')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can thus generate arbitrarily complex stimuli:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "stim = Stimulus([[0, 0, 20, 20, -5, -5, 0, 0, 0, 20, 20, -5, -5, 0, 0]],\n                time=[0, 0.1, 0.1, 0.2, 0.2, 0.6, 0.6, 1.0, 1.1, 1.1, 1.2, 1.2, 1.6, 1.6, 2.0])\nfig, ax = plt.subplots(figsize=(8, 5))\nax.plot(stim.time, stim.data[0, :])\nax.set_xlabel('Time (s)')\nax.set_ylabel('Amplitude ($\\mu$A)')"
      ]
    }
  ],
  "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.3"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}