<jupyter_start><jupyter_text>Exercise to "The Neural Basis of Episodic Memory" from the GRK Lecture Series<NAME><jupyter_code>import numpy as np
from scipy.io import loadmat
from scipy.signal import hilbert
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
%matplotlib inline
from pylab import rcParams
rcParams['figure.figsize'] = 8, 6
rcParams.update({'font.size': 15})

def nearest_neighbors_sorted(x, y):
    """ Locate timestamp pairs """
    x, y = map(np.asarray, (x, y))
    y_idx = np.argsort(y)
    y = y[y_idx]
    nearest_neighbor = np.empty((len(x),), dtype=np.intp)
    for j, xj in enumerate(x) :
        idx = np.searchsorted(y, xj)
        if idx == len(y) or idx != 0 and y[idx] - xj > xj - y[idx-1] :
            idx -= 1
        nearest_neighbor[j] = y_idx[idx]
        y = np.delete(y, idx)
        y_idx = np.delete(y_idx, idx)
    return nearest_neighbor<jupyter_output><empty_output><jupyter_text>1. Phase PrecessionLoad data and initialize variables<jupyter_code>fpath = '/home/kashefy/phd/grk_lecture_series/08_kempter/data.mat'

dat = loadmat(fpath)
print dat.keys()
lfp = dat['lfp']
lfp_t = dat['lfp_t']
pos_x = dat['pos_x']
pos_y = dat['pos_y']
pos_t = dat['pos_t']
spk_t = dat['spk_t']<jupyter_output>['pos_x', 'pos_y', 'lfp_t', '__header__', 'pos_t', '__globals__', 'lfp', '__version__', 'spk_t']<jupyter_text>a) Plot Trajectory<jupyter_code>plt.figure()
plt.plot(pos_x, pos_y, c='g', label=r'$location_{animal}$',
         linewidth=0.4, zorder=1)
spk_t_idx = nearest_neighbors_sorted(spk_t.ravel(), pos_t.ravel())
plt.scatter(pos_x[spk_t_idx], pos_y[spk_t_idx], label=r'$location_{spike}$',
            s=80, c='r', marker='.', zorder=2)
plt.title("animal's trajectory")
plt.xlabel(r'$position_{x}$')
plt.ylabel(r'$position_{y}$')
plt.legend(bbox_to_anchor=(1.15, 1.15))
plt.grid(True)

plt.figure()
# from spike counts to firing rates
sampling_rate = 1./np.mean(np.diff(pos_t, axis=0)) # sampling rate is 50 hz
print('sampling rate: %.3f Hz' % (sampling_rate,))
counts, bins = np.histogram(pos_x[spk_t_idx])
rates = []
for s, e, c in zip(bins[:-1], bins[1:], counts):
    in_bin = np.logical_and(pos_x >= s, pos_x < e)
    rates.append(c*sampling_rate/np.count_nonzero(in_bin)) 
plt.bar(bins[:-1], rates, width=30, align='center')
plt.title(r'firing rate | $position_{x}$')
plt.ylabel('firing rate [Hz]')
plt.xlabel(r'$position_{x}$')
plt.grid(True)

plt.figure()
dx = np.diff(pos_x, axis=0)
dx = np.vstack((dx, dx[-1]))
plt.scatter(pos_x[dx>=0], pos_y[dx>=0], c='b', marker='.', label='left-to-right',
         linewidth=0.5, alpha=0.4, edgecolor='none', zorder=1)
plt.scatter(pos_x[dx<0], pos_y[dx<0], c='m', marker='.', label='right-to-left',
         linewidth=0.5, alpha=0.4, edgecolor='none', zorder=1)
spk_t_idx = nearest_neighbors_sorted(spk_t.ravel(), pos_t.ravel())
plt.scatter(pos_x[spk_t_idx], pos_y[spk_t_idx], label=r'$location_{spike}$',
            s=80, c='r', marker='.', zorder=2)
plt.title("animal's trajectory")
plt.xlabel(r'$position_{x}$')
plt.ylabel(r'$position_{y}$')
plt.legend(bbox_to_anchor=(1.15, 1.15))
plt.grid(True)<jupyter_output>sampling rate: 50.000 Hz<jupyter_text>b) estimate spike phases<jupyter_code>dx_lr_idx, _ = np.where(np.diff(pos_x, axis=0) > 0) # positive speeds, from left to right
pos_t_lr = pos_t[dx_lr_idx]
spk_t_lr_idx = nearest_neighbors_sorted(spk_t.ravel(), pos_t_lr.ravel())

spk_t_lr = pos_t_lr[spk_t_lr_idx]
lfp_t_idx = nearest_neighbors_sorted(spk_t_lr.ravel(), lfp_t.ravel())

pos_x_lr = pos_x[dx_lr_idx]
phase = np.angle(hilbert(lfp.ravel()))
plt.plot(pos_x_lr[spk_t_lr_idx], phase[lfp_t_idx], '.')
plt.title(r'Spike phase vs. x position for left-right runs')
plt.xlabel('x position for left-right runs')
plt.ylabel('spike phase [rad]')<jupyter_output><empty_output><jupyter_text>Grid-cell Pattern Formation  a) Distribute input centers<jupyter_code>N = 900
L = 4.0 # [m]
a = 1.0
A_pos = 10.1
A_neg = 10.0
sigma_pos = 0.295 # [m]
sigma_neg = 0.3 # [m]
tau = 2.0 # [sec]
B = 0.3
dt = 0.050 # [sec]
T = 200 # [sec]

x, y = np.mgrid[-L/2:L/2:L/np.sqrt(N),-L/2:L/2:L/np.sqrt(N)]<jupyter_output><empty_output><jupyter_text>b) Compute Correlation Matrix<jupyter_code>def pdist(x1, x2, y1, y2, L):
    """ Compute periodic distance"""
    dx = np.abs(x1-x2)
    x_p = min(dx, L-dx)
    dy = np.abs(y1-y2)
    y_p = min(dy, L-dy)
    return np.sqrt(x_p*x_p + y_p*y_p)

c = np.empty([N, N])
for i in xrange(N):
    for j in xrange(N):
        d_ij = pdist(x.ravel()[i], x.ravel()[j], y.ravel()[i], y.ravel()[j], L)
        d_ij_sq = d_ij*d_ij
        c[i, j] = A_pos * np.exp(-d_ij_sq/(2. * sigma_pos * sigma_pos)) - A_neg * np.exp(-d_ij_sq / (2. * sigma_neg * sigma_neg))<jupyter_output><empty_output><jupyter_text>c) Simulate Weight Dynamic<jupyter_code>w_avg = B/(a-np.sum(c)/N)
w = np.random.rand(N) * 2. * w_avg
time_ = np.arange(0, T, dt)

t_plot = time_[np.linspace(0, len(time_), num=16, endpoint=False, dtype='int')] # time steps at which to plot weight vector
gs = gridspec.GridSpec(int(np.sqrt(len(t_plot))), int(np.sqrt(len(t_plot))))
gs.update(hspace=0.5, wspace=0.15)
plt.figure(figsize=(16, 12))

w2 = np.copy(w)

t_plot_idx = 0
for t in time_:    
    if t_plot_idx < len(t_plot) and t == t_plot[t_plot_idx]:
        plt.subplot(gs[t_plot_idx])
        plt.imshow(np.reshape(np.copy(w), (int(np.sqrt(N)), int(np.sqrt(N)))))
        plt.title('t=%.2fs' % (t,))
        plt.colorbar()
        t_plot_idx += 1
        
    w_t = np.copy(w)
    for i in xrange(N):
        w[i] = w_t[i] + dt/tau * (-a*w_t[i] + np.sum(c[i] * w_t + B))
    w = np.abs(w)
plt.suptitle('Evolution of weight vector at different time steps t')
plt.imshow(np.reshape(w, (int(np.sqrt(N)), int(np.sqrt(N)))))
plt.title('Final weight vector')
plt.xlabel("input center j")
plt.ylabel("input center i")
plt.colorbar()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.18554051738906974, 'fraction_numerical': 0.02158643138598595, 'mean_word_length': 4.5759312320916905, 'pattern_counts': {'":': 0, '<': 32, '<?xml version=': 0, '>': 33, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '589202', 'license': 'BSD-2-Clause', 'n_tokens_mistral': 2687, 'n_tokens_neox': 2538, 'n_words': 470, 'path': 'grk_kempter_episodic_mem.ipynb', 'repo_name': 'kashefy/grk_kempter_episodic_mem', 'size': '5845'}