<commit_before># -*- coding: utf-8 -*-
"""
================================================================
GLM for Spike Train Prediction in Primate Retinal Ganglion Cells
================================================================

* Original tutorial adapted from <NAME>, Princeton University
* Dataset provided by <NAME>, Stanford University
* The dataset is granted by the original authors for educational use only
* Please contact ``<EMAIL>`` if using beyond its purposes

The original MATLAB and Python tutorial can be found from
https://github.com/pillowlab/GLMspiketraintutorial.

These data were collected by <NAME> in the lab of
<NAME> at the Salk Institute. For full information see
Uzzell et al. [1]_, or Pillow et al. [2]_.

In this tutorial, we will demonstrate how to fit linear GLM and
Poisson GLM to predict the spike counts
recorded from primate retinal ganglion cells.
The dataset contains spike responses from 2 ON and 2 OFF parasol
retinal ganglion cells (RGCs) in primate retina, stimulated with
full-field `binary white noise`. Two experiment performed consisted of a
long (20-minute) binary stochastic (non-repeating) stimulus
which can be used for computing the spike-triggered average
(or characterizing some other model of the response).

References
----------
.. [1] <NAME>., and <NAME>. "Precision of
   spike trains in primate retinal ganglion cells."
   Journal of Neurophysiology 92.2 (2004)
.. [2] Pillow, <NAME>., et al. "Prediction and decoding of
   retinal ganglion cell responses with a probabilistic
   spiking model." Journal of Neuroscience 25.47 (2005)

"""

# Authors: <NAME> <<EMAIL>>
#          <NAME> <<EMAIL>>
# License: MIT

########################################################
#
# Import all the relevance libraries.

import os.path as op
import json

import numpy as np
from numpy.linalg import pinv
from scipy.linalg import hankel

from pyglmnet import GLM

import matplotlib.pyplot as plt
from tempfile import TemporaryDirectory

########################################################
#
# Fetch the dataset. The JSON file contains the
# following keys:  ``stim`` (binary stochastic stimulation),
# ``stim_times`` (time of the stimulation), and
# ``spike_times`` (recorded time of the spikes)

from pyglmnet.datasets import fetch_rgc_data

with TemporaryDirectory(prefix="tmp_glm-tools") as temp_dir:
    dpath = fetch_rgc_data(dpath=temp_dir, accept_rgcs_license=True)
    with open(op.join(dpath, 'data_RGCs.json'), 'r') as f:
        rgcs_dataset = json.loads(f.read())

stim = np.array(rgcs_dataset['stim'])
stim_times = np.array(rgcs_dataset['stim_times'])

n_cells = len(rgcs_dataset['spike_times'])
dt = stim_times[1] - stim_times[0]  # time between the stimulation
n_times = len(stim)  # total number of the stimulation
sfreq = 1. / dt  # frequency of the stimulation

########################################################
#
# You can pick a cell to work with and visualize the spikes for one second.
# In this case, we will pick cell number 2 (ON cell).

cell_idx = 2  # pick cell number 2 (ON cell)
spike_time = rgcs_dataset['spike_times']['cell_%d' % cell_idx]
spike_time = np.array(spike_time)
n_spikes = len(spike_time)  # number of spikes

# bin the spikes
t_bins = np.arange(n_times + 1) * dt
spikes_binned, _ = np.histogram(spike_time, t_bins)

print('Loaded RGC data: cell {}'.format(cell_idx))
print('Number of stim frames: {:d} ({:.1f} minutes)'.
      format(n_times, n_times * dt / 60))
print('Time bin size: {:.1f} ms'.format(dt * 1000))
print('Number of spikes: {} (mean rate = {:.1f} Hz)\n'.
      format(n_spikes, n_spikes / n_times * 60))

########################################################
#
# Visualize first 120 samples
sample_idx = np.arange(120)
t_sample = dt * sample_idx
tmax = t_sample.max()

fig, axes = plt.subplots(3, 1, sharex=True)
axes[0].step(t_sample, stim[sample_idx])
axes[0].set_xlim([0., tmax])
axes[0].set_title('Raw Stimulus, spikes, and'
                  ' spike counts')
axes[0].set_ylabel('Stimulation Intensity')

tspplot = spike_time[(spike_time >= t_sample.min()) &
                     (spike_time < t_sample.max())]
markerline, _, _ = axes[1].stem(spike_time[sample_idx],
                                [1] * len(spike_time[sample_idx]))
markerline.set_markerfacecolor('none')
axes[1].set_xlim([0., tmax])
axes[1].set_ylim([0, 2])
axes[1].set_ylabel('Spikes')

markerline, _, _ = axes[2].stem(t_sample, spikes_binned[sample_idx])
markerline.set_markerfacecolor('none')
axes[2].set_xlim([0., tmax])
axes[2].set_xlabel('Time (s)')
axes[2].set_ylabel('Binned Spike counts')
plt.show()

########################################################
#
# We can use ``scipy``'s function ``hankel`` to make our design matrix.
# Design matrix :math:`X` can be created using the stimulation and its history.
# Later in the tutorial, we will also incorporate spikes history into our
# design matrix.

n_t_filt = 25  # tweak this to see different results
stim_padded = np.pad(stim, (n_t_filt - 1, 0))
Xdsgn = hankel(stim_padded[0: -n_t_filt + 1], stim[-n_t_filt:])

plt.figure()
plt.imshow(Xdsgn[:50, :],
           cmap='binary',
           aspect='auto',
           interpolation='nearest')
plt.xlabel('lags before spike time')
plt.ylabel('time bin of response')
plt.title('Sample first 50 rows of design'
          ' matrix created using Hankel')
plt.show()

########################################################
# **Fitting and predicting with a linear-Gaussian GLM**
#
# For a general linear model, an observed spikes can be
# thought of an underlying parameter
# :math:`\beta_0, \beta` that control the spiking:
# :math:`y = (\beta_0 + X \beta) \ + \epsilon` and
# :math:`y \sim \text{Poiss}(\beta_0 + X \beta)` where
# :math:`X` is the stimulation and history of stimulation
#
# We can add an offset or a "constant" to our design matrix.
# This can be done by concatenating a column of 1 to
# our previously created design matrix.

Xdsgn_offset = np.hstack((np.ones((n_times, 1)), Xdsgn))

# compute whitened spike-triggered average (STA)
wsta_offset = pinv(Xdsgn_offset) @ spikes_binned

# predict spike counts
spikes_pred_lgGLM_offset = Xdsgn_offset @ wsta_offset

########################################################
# **Fitting and predicting with a Gaussian GLM**
#
# We can also assume that their is a non-linear function governing
# the underlying the firing patterns. Concreately, we can write down as
# :math:`y = q(\beta_0 + X \beta) + \epsilon`, and
# :math:`y \sim \text{Poiss}(\beta_0 + X \beta)`.
# Here we use :math:`exp()` as the inverse link function :math:`q`.
# We can use the :class:`GLM` class to predict the parameters.

# create possion GLM instance
glm_poisson = GLM(distr='poisson',
                  alpha=0.05,
                  learning_rate=1.0,
                  score_metric='pseudo_R2',
                  reg_lambda=1e-7)
glm_poisson.fit(Xdsgn, spikes_binned)
spikes_pred_poissonGLM = glm_poisson.predict(Xdsgn)

#############################################################################
# **Adding spikes history for predicting spike counts**
#
# We can even further predict the spikes by concatenating the spikes history.
#
# .. note::
#    The spike-history portion of the design
#    matrix had better be shifted so that we aren't allowed to
#    use the spike count on this time bin to predict itself!

n_t_filt = 25  # same as before, stimulation history
n_t_hist = 20  # spikes history

# using both stimulation history and spikes history
spikes_padded = np.pad(spikes_binned, (n_t_hist, 0))

Xstim = hankel(stim_padded[:-n_t_filt + 1], stim[-n_t_filt:])
Xspikes = hankel(spikes_padded[:-n_t_hist], stim[-n_t_hist:])
Xdsgn_hist = np.hstack((Xstim, Xspikes))  # design matrix with spikes history

# Now, we are ready to fit Poisson GLM with spikes history.
# create possion GLM instance
glm_poisson_hist = GLM(distr='poisson', alpha=0.05,
                       learning_rate=1.0,
                       score_metric='pseudo_R2',
                       reg_lambda=1e-7)
<commit_msg>Using Gaussian GLM instead<commit_after># -*- coding: utf-8 -*-
"""
================================================================
GLM for Spike Train Prediction in Primate Retinal Ganglion Cells
================================================================

* Original tutorial adapted from <NAME>, Princeton University
* Dataset provided by <NAME>, Stanford University
* The dataset is granted by the original authors for educational use only
* Please contact ``<EMAIL>`` if using beyond its purposes

The original MATLAB and Python tutorial can be found from
https://github.com/pillowlab/GLMspiketraintutorial.

These data were collected by <NAME> in the lab of
E.J. Chichilnisky at the Salk Institute. For full information see
Uzzell et al. [1]_, or Pillow et al. [2]_.

In this tutorial, we will demonstrate how to fit linear GLM and
Poisson GLM to predict the spike counts
recorded from primate retinal ganglion cells.
The dataset contains spike responses from 2 ON and 2 OFF parasol
retinal ganglion cells (RGCs) in primate retina, stimulated with
full-field `binary white noise`. Two experiment performed consisted of a
long (20-minute) binary stochastic (non-repeating) stimulus
which can be used for computing the spike-triggered average
(or characterizing some other model of the response).

References
----------
.. [1] <NAME>., and <NAME>. "Precision of
   spike trains in primate retinal ganglion cells."
   Journal of Neurophysiology 92.2 (2004)
.. [2] Pillow, <NAME>., et al. "Prediction and decoding of
   retinal ganglion cell responses with a probabilistic
   spiking model." Journal of Neuroscience 25.47 (2005)

"""

# Authors: <NAME> <<EMAIL>>
#          <NAME> <<EMAIL>>
# License: MIT

########################################################
#
# Import all the relevance libraries.

import os.path as op
import json

import numpy as np
from scipy.linalg import hankel

from pyglmnet import GLM

import matplotlib.pyplot as plt
from tempfile import TemporaryDirectory

########################################################
#
# Fetch the dataset. The JSON file contains the
# following keys:  ``stim`` (binary stochastic stimulation),
# ``stim_times`` (time of the stimulation), and
# ``spike_times`` (recorded time of the spikes)

from pyglmnet.datasets import fetch_rgc_data

with TemporaryDirectory(prefix="tmp_glm-tools") as temp_dir:
    dpath = fetch_rgc_data(dpath=temp_dir, accept_rgcs_license=True)
    with open(op.join(dpath, 'data_RGCs.json'), 'r') as f:
        rgcs_dataset = json.loads(f.read())

stim = np.array(rgcs_dataset['stim'])
stim_times = np.array(rgcs_dataset['stim_times'])

n_cells = len(rgcs_dataset['spike_times'])
dt = stim_times[1] - stim_times[0]  # time between the stimulation
n_times = len(stim)  # total number of the stimulation
sfreq = 1. / dt  # frequency of the stimulation

########################################################
#
# You can pick a cell to work with and visualize the spikes for one second.
# In this case, we will pick cell number 2 (ON cell).

cell_idx = 2  # pick cell number 2 (ON cell)
spike_time = rgcs_dataset['spike_times']['cell_%d' % cell_idx]
spike_time = np.array(spike_time)
n_spikes = len(spike_time)  # number of spikes

# bin the spikes
t_bins = np.arange(n_times + 1) * dt
spikes_binned, _ = np.histogram(spike_time, t_bins)

print('Loaded RGC data: cell {}'.format(cell_idx))
print('Number of stim frames: {:d} ({:.1f} minutes)'.
      format(n_times, n_times * dt / 60))
print('Time bin size: {:.1f} ms'.format(dt * 1000))
print('Number of spikes: {} (mean rate = {:.1f} Hz)\n'.
      format(n_spikes, n_spikes / n_times * 60))

########################################################
#
# Visualize first 120 samples
sample_idx = np.arange(120)
t_sample = dt * sample_idx
tmax = t_sample.max()

fig, axes = plt.subplots(3, 1, sharex=True)
axes[0].step(t_sample, stim[sample_idx])
axes[0].set_xlim([0., tmax])
axes[0].set_title('Raw Stimulus, spikes, and'
                  ' spike counts')
axes[0].set_ylabel('Stimulation Intensity')

tspplot = spike_time[(spike_time >= t_sample.min()) &
                     (spike_time < t_sample.max())]
markerline, _, _ = axes[1].stem(spike_time[sample_idx],
                                [1] * len(spike_time[sample_idx]))
markerline.set_markerfacecolor('none')
axes[1].set_xlim([0., tmax])
axes[1].set_ylim([0, 2])
axes[1].set_ylabel('Spikes')

markerline, _, _ = axes[2].stem(t_sample, spikes_binned[sample_idx])
markerline.set_markerfacecolor('none')
axes[2].set_xlim([0., tmax])
axes[2].set_xlabel('Time (s)')
axes[2].set_ylabel('Binned Spike counts')
plt.show()

########################################################
#
# We can use ``scipy``'s function ``hankel`` to make our design matrix.
# Design matrix :math:`X` can be created using the stimulation and its history.
# Later in the tutorial, we will also incorporate spikes history into our
# design matrix.

n_t_filt = 25  # tweak this to see different results
stim_padded = np.pad(stim, (n_t_filt - 1, 0))
Xdsgn = hankel(stim_padded[0: -n_t_filt + 1], stim[-n_t_filt:])

plt.figure()
plt.imshow(Xdsgn[:50, :],
           cmap='binary',
           aspect='auto',
           interpolation='nearest')
plt.xlabel('lags before spike time')
plt.ylabel('time bin of response')
plt.title('Sample first 50 rows of design'
          ' matrix created using Hankel')
plt.show()

########################################################
# **Fitting and predicting with a linear-Gaussian GLM**
#
# For a general linear model, an observed spikes can be
# thought of an underlying parameter
# :math:`\beta_0, \beta` that control the spiking:
# :math:`y = (\beta_0 + X \beta) \ + \epsilon` and
# :math:`y \sim \text{Poiss}(\beta_0 + X \beta)` where
# :math:`X` is the stimulation and history of stimulation
#
# You can simply use linear Gaussian GLM to fit as follows

glm_lgGLM = GLM(distr='gaussian')
glm_lgGLM.fit(Xdsgn, spikes_binned)

# predict spike counts
spikes_pred_lgGLM_offset = glm_lgGLM.predict(Xdsgn)

########################################################
# **Fitting and predicting with a Gaussian GLM**
#
# We can also assume that their is a non-linear function governing
# the underlying the firing patterns. Concreately, we can write down as
# :math:`y = q(\beta_0 + X \beta) + \epsilon`, and
# :math:`y \sim \text{Poiss}(\beta_0 + X \beta)`.
# Here we use :math:`exp()` as the inverse link function :math:`q`.
# We can use the :class:`GLM` class to predict the parameters.

# create possion GLM instance
glm_poisson = GLM(distr='poisson',
                  alpha=0.05,
                  learning_rate=1.0,
                  score_metric='pseudo_R2',
                  reg_lambda=1e-7)
glm_poisson.fit(Xdsgn, spikes_binned)

# predict spike counts
spikes_pred_poissonGLM = glm_poisson.predict(Xdsgn)

#############################################################################
# **Adding spikes history for predicting spike counts**
#
# We can even further predict the spikes by concatenating the spikes history.
#
# .. note::
#    The spike-history portion of the design
#    matrix had better be shifted so that we aren't allowed to
#    use the spike count on this time bin to predict itself!

n_t_filt = 25  # same as before, stimulation history
n_t_hist = 20  # spikes history

# using both stimulation history and spikes history
spikes_padded = np.pad(spikes_binned, (n_t_hist, 0))

Xstim = hankel(stim_padded[:-n_t_filt + 1], stim[-n_t_filt:])
Xspikes = hankel(spikes_padded[:-n_t_hist], stim[-n_t_hist:])
Xdsgn_hist = np.hstack((Xstim, Xspikes))  # design matrix with spikes history

# Now, we are ready to fit Poisson GLM with spikes history.
# create possion GLM instance
glm_poisson_hist = GLM(distr='poisson', alpha=0.05,
                       learning_rate=1.0,
                       score_metric='pseudo_R2',
                       reg_lambda=1e-7)

{'fraction_non_alphanumeric': 0.18177825186672522, 'fraction_numerical': 0.013302378113823178, 'mean_word_length': 4.375379426644182, 'pattern_counts': {'":': 0, '<': 32, '<?xml version=': 0, '>': 32, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2120628', 'n_tokens_mistral': 5337, 'n_tokens_neox': 4907, 'n_words': 1607}