<jupyter_start><jupyter_text>PyStorm Walkthrough<jupyter_code># basic imports for working with matrices and plotting
import numpy as np
import matplotlib.pyplot as plt
import scipy.interpolate
import nengo<jupyter_output><empty_output><jupyter_text>This guide will introduce you to PyStorm, a Python package that exposes abstractions of Braindrop's hardware constructs as network objects. Basic familiarity with PyStorm's network objects is assumed (see first 2 pages of __braindrop_pystorm_documentation.pdf__). The outline for this walkthrough is as follows:1. Creating a 1D pool and its tap-point matrix2. Sending inputs and measuring spikes3. Decoding functions4. Visualizing tuning curves *** NOTE ***Below we've provided a cell for resetting the chip's connection to the server. In most cases, this code will be unnecessary. However, sometimes when the initial handshake occurs (during the `HAL()` call in the tutorial code below), the USB connection fails and the Python call quickly raises the error messeage `Comm failed to init`. If this happens, run this cell below to reset the connection and try again. This may take a few tries for it to reconnect.<jupyter_code>!sudo /home/ee207student/usbreset.sh # reset the USB connection, which sometimes dies.
# if Comm below is repeatedly & quickly failing to init, you can run this cell.<jupyter_output>0000:08:00.0
tee: /sys/bus/pci/drivers/xhci_hcd/bind: No such device
Resetting devices from /sys/bus/pci/drivers/uhci_hcd...
/home/ee207student/usbreset.sh: line 20: echo: write error: No such device
/home/ee207student/usbreset.sh: line 21: echo: write error: No such device
Resetting devices from /sys/bus/pci/drivers/xhci_hcd...<jupyter_text>*** END NOTE ***  PyStorm Imports<jupyter_code>############################################################################
# Importing essential PyStorm libraries

# HAL, or Hardware Abstraction Layer, must be imported and initialized
from pystorm.hal import HAL
HAL = HAL()

# data_utils contains a useful function for binning spikes and estimating spike rates
from pystorm.hal import data_utils

# RunControl objects are responsible for executing command sequences on the chip, such as initiating
# tuning curve measurements
from pystorm.hal.run_control import RunControl

# Graph objects represent properties of your network
# so that HAL can allocate resources for them on the chip
from pystorm.hal.neuromorph import graph<jupyter_output><empty_output><jupyter_text>1 Creating a 1D Pool and Tap Point MatrixLet's review basic requirements for building and simulating an ensemble with PyStorm (for syntax definitions, see __braindrop_pystorm_documentation.pdf__).  1.1 Defining a pool's parameters and its tap-point matrixDefine your neuron array's width $W$ and height $H$, such that the ensemble's size is $N = WH$, and the dimensionality $D$ of its input space. Then, construct a $N$-by-$D$ tap-point matrix `tap_matrix`, which specifies an _anchor encoder_ for each tap-point used. For instance, in Problem Set 3, you defined tap points on a 1-D diffuser with plus and minus signs ($+1$ or $-1$). In $D$ dimensions, these anchor encoders are $D$-dimensional vectors in $\{-1, 0, +1\}^D$. When we use a 2-D diffuser to create two-dimensional encoders, we draw the tap-points' anchor encoders from $\{-1, 0, +1\}^2$.Although a $N$-by-$D$ matrix can specify all $N$ anchor encoders, a $N$-neuron pool has only $N/4$ tap-points—since there is 1 synaptic filter per $2\times2$ cluster of neurons. If you specify more than $N/4$ anchor encoders, HAL will throw an error when it tries to create your network. The code below implements an algorithm that assigns 1-D anchor encoders to tap-points. It's quite simple:- Divide the set of tap-points in half along the width axis - Assign one side $[+1]$ anchor encoders and the other half $[-1]$ anchor encoders. `plot_tap_matrix()` visualizes the results within the $W$-by-$H$ neuron array. Free parameter $\rho$ controls tap-point density by skipping synapses. Try changing this parameter to create different tap-point matrices.<jupyter_code>def tap_mat_1d(width, height, rho=0)
    '''
    Creates a structured 1D tap-point matrix
    
    Parameters
    ----------
    width: int
        Width of neuron array
    height: int
        Height of neuron array
    rho: int
        1/2**rho is the fraction of total tap-points the 
        algorithm should use
    
    Returns
    -------
    tap_matrix: numpy matrix
        width*height x 2 tap-point matrix
    '''
    # creating an empty tap-point matrix
    tap_matrix = np.zeros((N, Din))

    # one synaptic filter per 4 neurons
    for x in range(2**rho - 1, width, 2 * 2**rho):
        for y in range(2**rho - 1, height, 2 * 2**rho):
            n = y * width + x
            if x < width // 2:
                tap_matrix[n, 0] = 1
            else:
                tap_matrix[n, 0] = -1
    return tap_matrix

def plot_tap_matrix(matrix, width, height):
    '''
    Plots single dimension of a tap-point matrix with a given width and height
    '''
    plt.imshow(matrix.reshape((height, width)), cmap='bwr', vmin=-1, vmax=1, interpolation='none', aspect='equal')
    plt.colorbar(fraction=0.046, pad=0.04)
    ax = plt.gca()
    ax.set_yticks(np.arange(0, height, 1))
    ax.set_xticks(np.arange(0, width, 1))
    ax.set_yticklabels(np.arange(1, height+1, 1));
    ax.set_xticklabels(np.arange(1, width+1, 1));
    ax.set_yticks(np.arange(-.5, height, 1), minor=True);
    ax.set_xticks(np.arange(-.5, width, 1), minor=True);

    plt.grid(which='minor', color='k', linestyle='--')

############################################################################
# Creating pool parameters and a tap-point matrix

# width, height, and size of a 1D pool
width = 16
height = 16
width_height = (width, height)
N = width * height
Din = 1
rho = 0
tap_matrix = tap_mat_1d(width, height, rho=rho)

# plotting the matrix
plt.figure(figsize=(6,6))
plot_tap_matrix(tap_matrix, width, height)<jupyter_output><empty_output><jupyter_text>1.2 Creating the pool and its connectionsNow that we've selected pool's width (16), height (16), dimensionality (1), and tap-point matrix, we may create an actual network. We do this by defining input and output nodes, for incoming and outgoing signals, a pool of neurons,  a bucket to accumulate its neurons' spikes, and then connect these network objects to create a graph object.<jupyter_code>############################################################################
# Creating a pool and its connections

# create the network as a graph object, named "net" here, which tracks connections and configuration parameters
# for the network's objects.
net = graph.Network("net")

# decoders are initially zero, we remap them later (without touching the rest of the network)
# using HAL.remap_weights()
Dout = 1
decoders = np.zeros((Dout, N))

# create and name the input and specify its dimensionality
i1 = net.create_input("i1", Din)

# create and name the neuron pool and supply the computed tap-point matrix.
# we also explicitly set the biases here; you will optimize these in
# the next problem. All neurons have their bias set to -3.
p1 = net.create_pool("p1", tap_matrix, biases=-3)

# create and name the bucket and specify its dimensionality
b1 = net.create_bucket("b1", Dout)

# create and name the output and specify its dimensionality
o1 = net.create_output("o1", Dout)

# connect the input to the pool
# "None" is for the weights field; only connections to buckets can have weights
net.create_connection("c_i1_to_p1", i1, p1, None)

# connect the pool to the bucket and specify initial weights
decoder_conn = net.create_connection("c_p1_to_b1", p1, b1, decoders)

# connect the bucket to the output
net.create_connection("c_b1_to_o1", b1, o1, None)<jupyter_output><empty_output><jupyter_text>1.3 Setting up HAL to map the defined networkHAL instantiates our network on Braindrop and allocates resources to it. To do this, we have to first configure Braindrop's and the host machine's time-bins (i.e., time-resolution). Then, we `map(net)`, the Graph object (defined in 2.1.3). We must redo this call any time we change `net`'s properties.<jupyter_code>############################################################################
# Setting up the HAL interface and mapping the network

# misc driver parameters: downstream refers to Braindrop and upstream refers to the host
downstream_time_res = 10000 # ns, or 10us
upstream_time_res = 1000000 # ns, or 1ms

HAL.set_time_resolution(downstream_time_res, upstream_time_res)

# invoke HAL.map(), make tons of neuromorph/driver calls under the hood
# to map the network
print("calling map")
HAL.map(net)<jupyter_output>calling map<jupyter_text>2 Sending Inputs and Measuring Spike RatesIn order to drive Braindrop with a stimulus $x(t)$, PyStorm requires you to specifiy its values at particular times they occur and a value for the parameter $f_\mathrm{max}$. The stimulus $-1 < x(t) < +1$ is streamed to Braindrop as a periodic train of impulses with instantaneous rate $-f_\mathrm{max} < f_\mathrm{max}x(t) < +f_\mathrm{max}$, where the minus sign denotes an impulse with area -1 instead of +1.  2.1 Example input signalsThe functions `make_cons_stimulus()` and `make_sine_stimulus()` demonstrate two different inputs that are useful for measuring neurons' firing properties.`make_cons_stimulus()` creates a D-dimensional stimulus that steps through all points in a uniformly-gridded space $[-f_\mathrm{max},+f_\mathrm{max}]^D$. This signal's value equals each of the $Q^D$ points for a constant amount of time. It's used to sweep a $D$-dimensional space and measure the neurons' responses. That yields estimates for their firing-rates using bin-widths equal to `hold_time`. Such box-car filtering isn't what the synapse does, however. Thus it yields a wrong estimate for firing-rates. Training decoders using these firing-rate estimates therefore does not produce decoders that operate with the synapses' firing-rate estimates.`make_sine_stimulus()` creates a multidimensional sinusoidal stimulus with amplitude $f_\mathrm{max}$, period `base_period`, and duration given in `cycles`. Note that it's constant over a step size `input_dt`. Each dimension's ($D_i$) sinusoid oscillates at frequency $\phi^{-D_i}/T$, where $\phi$ is the golden ratio $\frac{1+\sqrt{5}}{2}$ and $T$ is the first-dimension's period, which `base_period` signifies. This scaling ensures no pair of dimensions ever synchronize. As a result, the full domain $[-1,1]^D$ is swept. When used to train decoders with `base_period` equal to $2\pi\tau_s$, where $\tau_s$ is the synapses' time-constant, it yields reasonable firing-rate noise estimates (after lowpass filtering time-constant $\tau_s$). Training decoders will be demonstrated in section 3.<jupyter_code>def make_cons_stimulus(Din, hold_time, points_per_dim, fmax=1000):
    '''
    Creates multidimensional staircase-like signal
    
    Parameters
    ----------
    Din: int
        Input dimensionality
    hold_time: float (seconds)
        Duration that each input is applied
    points_per_dim: int
        Number of points per dimension
    fmax: float
        Units converting x = +/-1 to firing rate
    
    Returns
    -------
    inp_times: numpy array
        Length-S array of times a new input starts (in seconds) 
    inputs: numpy matrix
        S-by-Din matrix of inputs (in units of firing rate)
    '''
    total_points = points_per_dim ** Din
    stim_rates_1d = np.linspace(-fmax, fmax, points_per_dim).astype(int)
    if Din > 1:
        stim_rates_mesh = np.meshgrid(*([stim_rates_1d]*Din))
    else:
        stim_rates_mesh = [stim_rates_1d]
    stim_rates = [r.flatten() for r in stim_rates_mesh]
    inputs = np.array(stim_rates).T
    
    inp_times = np.arange(0, total_points) * hold_time
    return inp_times, inputs

def make_sine_stimulus(Din, base_period, cycles, fmax=1000, input_dt=1e-3):
    '''
    Creates multidimensional sinusoidal input
    
    Parameters
    ----------
    Din: int
        Input dimensionality
    base_period: float (seconds)
        Period of oscillation for first dimension
    cycles: float
        Number of cycles the signal is generated for
    fmax: float
        Units converting x = +/-1 to firing rate
    input_dt: float
        Step-size between each consecutive input (in seconds)
        
    Returns
    -------
    inp_times: numpy array
        Length-S array of times a new input starts (in seconds) 
    inputs: numpy matrix
        S-by-Din matrix of inputs (in units of firing rate)
    '''
    phi = (1 + np.sqrt(5))/2.0
    freq = 1/base_period
    f = ((1/phi)**np.arange(Din)) * freq
    stim_func = lambda t: np.sin(2*np.pi*f*t)
    
    sim_time = cycles * base_period

    inp_times = np.arange(int(sim_time/input_dt))*input_dt
    inputs = np.zeros((len(inp_times), Din))

    for index, t in enumerate(inp_times):
        v = stim_func(t)

        for i, vv in enumerate(v):
            ff = int(vv*fmax)
            if ff == 0:
                ff = 1

            inputs[index, i] = float(ff)

    return inp_times, inputs

############################################################################
# Plotting example stimuli

Dstim = 2
hold_time = 0.1
Q = 10
times, inp = make_cons_stimulus(Dstim, hold_time, Q)
plt.figure()
plt.plot(inp[:,0], inp[:,1], '--o')
plt.xlabel(r'$x_1$')
plt.ylabel(r'$x_2$')
plt.title('Constant Sweep, Total Length = %.1f sec' % times[-1])

Dstim = 2
base_period = 0.1
cycles = 30
times, inp = make_sine_stimulus(Dstim, base_period, cycles)
plt.figure()
plt.plot(inp[:,0], inp[:,1])
plt.xlabel(r'$x_1$')
plt.ylabel(r'$x_2$')
plt.title('Sinusoidal Sweep, Total Length = %.1f sec' % times[-1]);<jupyter_output><empty_output><jupyter_text>2.2 Sending inputs and collecting spikesTo send stimuli to Braindrop, we use the `RunControl` object, telling it what the stimulus is, where it should go, and for how long it should be measured. It returns spike or output data back to the host, binned with a resolution defined by the `upstream_time_res` parameter (set earlier).In this example, you will use `make_sine_stimulus()` as stimulus to your pool, with `base_period` equal to $2\pi\tau_s$, where $\tau_s$ is the average synaptic time-constant. In order to select the number of cycles, note that each cycle results in approximately 2 uncorrelated samples of the input space. If you want $Q$ samples per dimension, you should use $(Q/2)^D$ cycles to measure the firing rate. To be conservative, quadruple this number.<jupyter_code>############################################################################
# Sending inputs and collecting spikes

# define input parameters
Din = 1
fmax = 1000 # Hz
input_dt = 1e-3 # s
tau_syn = 0.02 # 20 ms
base_period = 2*np.pi*tau_syn # approximately 125 ms

# approximate samples per dimension
Q = 20
cycles = 4* ((Q / 2)**Din) # 4*10^D cycles, or 5 sec for D=1, 50 sec for D=2

# create the input
inp_times, inputs_sin = make_sine_stimulus(Din, base_period, cycles, fmax=fmax, input_dt=input_dt)

# package together the inputs and their timestamps
# times must be sent in ns (and thus are multiplied by 1e9)
FUDGE = 2
curr_time = HAL.get_time()
times = curr_time + FUDGE * 1e9 + inp_times * 1e9
times_w_end = np.hstack((times, times[-1] + input_dt * 1e9))
input_vals = {i1 : (times, inputs_sin)}

# create the RunControl object for our defined network
rc = RunControl(HAL, net)

# send inputs and collect spikes
print("getting spikes")
_, spikes_and_bin_times = rc.run_input_sweep(input_vals, get_raw_spikes=True, 
                                             end_time=times_w_end[-1], rel_time=False)
print("done")

# spike_bin_times is length-S, bin-size set by upstream_time_res
# spike_dict is a dictionary, spike_ct is a S-by-N matrix of spike counts within S bins
spike_dict, spike_bin_times = spikes_and_bin_times
spike_ct = spike_dict[p1]

# convert spike counts to spike trains
spikes = spike_ct / input_dt

# truncate inputs to length of measured outputs
inp_times_trunc = inp_times[:spikes.shape[0]]
inputs_sin_trunc = inputs_sin[:spikes.shape[0]]<jupyter_output>getting spikes
done<jupyter_text>3 Decoding FunctionsLet's analyze how to decode functions using the above architecture as an example. It depicts a `graph` object with two pools, `poolA` and `poolB`, each with their own set of tap-point matrices and synapses. The input node sends spikes at rate $f_\mathrm{max}x(t)$ to `poolA`'s synapses, specified by the pool's tap-point matrix. The synaptic current spreads through the diffuser and drives neurons in `poolA` to spike. These spikes are collected in `bucketAB` and thinned to decode a desired function $f(x)$. Once `bucketAB` emits its a decoded spike-train, it is filtered through `poolB`'s synapses (specified in its tap-point matrix) before spreading through the diffuser and driving `poolB`'s neurons. If we want to compute and output some other function $g(x)$ on `poolB`'s spike-train, we connect `poolB` to `bucketBOut`. This bucket thins the spike-train with its own set of decode weights, and the result is sent to the host through an output node.We first build this graph.input is u(t)input to pool A is $\tau_a\dot{x} + x = u(t)$Apply $f$ to $x(t)$, $f(x(t))$Then filter throught $\tau_b$, $\tau_b\dot{y} + y = f(x(t))$<jupyter_code>DinA = 1
DoutA = 1
widthA = 16
heightA = 16
rhoA = 0
Na = widthA * heightA
tap_mat_A = tap_mat_1d(widthA, heightA, rho=rhoA)
xy_loc_A = (0, 0) # select pool starting at top-left = (0, 0)

DinB = 1
DoutB = 1
widthB = 16
heightB = 16
rhoB = 0
Nb = widthB * heightB
tap_mat_B = tap_mat_1d(widthB, heightB, rho=rhoB)
xy_loc_B = (16, 0) # select pool starting at top-left = (16,0)

# decoders are initially zero, we remap them later (without touching the rest of the network)
# using HAL.remap_weights()
decoders_AB = np.zeros((DoutA, Na))
decoders_BOut = np.zeros((DoutB, Nb))

net = graph.Network("net")
i1 = net.create_input("i1", DinA)
pA = net.create_pool("pA", tap_mat_A, biases=-3, user_xy_loc=xy_loc_A)
pB = net.create_pool("pB", tap_mat_B, biases=-3, user_xy_loc=xy_loc_B)

bAB = net.create_bucket("bAB", DoutA)
bBOut = net.create_bucket("bBOut", DoutB)

o1 = net.create_output("o1", Dout)

net.create_connection("c_i1_to_p1", i1, pA, None)
decoder_AB_conn = net.create_connection("c_pA_to_bAB", pA, bAB, decoders_AB)
net.create_connection("c_bAB_to_pB", bAB, pB)
decoder_BOut_conn = net.create_connnection("c_pB_to_bBOut", pB, bBOut, decoders_BOut)
net.create_connection("c_bBOut_to_o1", bBOut, o1, None)

HAL.map(net)<jupyter_output><empty_output><jupyter_text>3.1 Training decodersTo train decoders, first filter the input signal $x(t)$ by the synaptic time-constant (since the spikes were generated by filtering the input through the synapse as well), and then apply the transformation $f(x)$ to the resulting time-series. You could then take your spike-train matrix $\mathbf{S}$ and transformed and filtered input $\mathbf{f}(x_\mathrm{filt})$ find decoders $\mathbf{d}$ such that $\mathbf{S}\mathbf{d} = \mathbf{f}(x_\mathrm{filt})$. However, decoding directly on the spike train will result in extremely noisy function decoders. Instead, you should convolve both the spikes and the the target function with a longer time constant filter. This procedure will give a $S$-by-$N$ matrix of firing rates $\mathbf{A}$ and a length-$S$ vector of target firing-rates $\mathbf{f}_\mathrm{filt}(x_\mathrm{filt})$, allowing you to solve $\mathbf{A}\mathbf{d} = \mathbf{f}_\mathrm{filt}(x_\mathrm{filt})$ for decoders $\mathbf{d}$. After finding this decoder, you must check to ensure that all decoding weights are at most 1, since this is the largest weight that can be implemented by thinning. If this is not the case, rescale all the decoders by $\max(|\mathbf{d}|)$ and redefine $f_\mathrm{max}' = f_\mathrm{max}/\max(|\mathbf{d}|)$. This will be used to rescale our firing rate back to the units of the function space by dividing the decoded output by $f_\mathrm{max}'$ instead of $f_\mathrm{max}$.<jupyter_code>############################################################################
# Training decoders and evaluating training error

# create filter objects
f_in = nengo.synapses.Lowpass(tau_syn)
tau_out = 0.1
f_out = nengo.synapses.Lowpass(tau_out)

# here we demonstrate the target function f(x) = x
func = lambda x: x

# first filter the input signal
x_filt = f_in.filt(inputs_sin_trunc, y0=0, dt=input_dt)
# then apply the function
fx_filt = func(x_filt)

# create a firing-rate matrix by filtering with f_out
A = f_out.filt(spikes, y0=0, dt=input_dt)

# create a filtered target function
f_filt = f_out.filt(fx_filt, y0=0, dt=input_dt)

# use nengo's built-in solver
solver = nengo.solvers.LstsqL2(reg=0.0001)
# use solver to find decoders
dec, _ = solver(A, f_filt)

# make sure decoders are <= 1
max_d = np.max(np.abs(dec))
if max_d > 1:
    # if larger, rescale
    fmax_out = fmax / max_d
    dec /= max_d
    print("rescaling by", max_d)
else:
    # otherwise, we use the same scaling
    fmax_out = fmax

# apply the decode on the training data
fhat_filt = A @ dec

# rescale to output units
fhat_filt_out = fhat_filt / fmax_out
f_filt_out = f_filt / fmax

# compute training error
rmse = np.sqrt(np.mean((f_filt_out - fhat_filt_out)**2))

# plotting 4 periods of data
plt.figure()
plt.plot(inp_times_trunc[:int(4*base_period/input_dt)], 
         fhat_filt_out[:int(4*base_period/input_dt)], label='decoded')
plt.plot(inp_times_trunc[:int(4*base_period/input_dt)], 
         f_filt_out[:int(4*base_period/input_dt)], '--k', label='true')
plt.xlabel(r'$t$ [s]')
plt.legend()
plt.title('RMSE = %.3f' % rmse)<jupyter_output><empty_output><jupyter_text>3.2 Changing the decoder and testing with new inputsIn order to update the network's decoder to the weights from 3.1, you must reassign the weights to the decoder connection (which connects the pool to the bucket) and call HAL to remap the network with the new weights.To test this model, we've chosen the sinusoidal signal generated by `make_sine_stimulus()`, this time with a base period of 1 second for a duration of 4 cycles. We package this input along with its timestamps as in section 2, and use the RunControl object to send these inputs to Braindrop. Instead of measuring spikes directly from the pool, we set `raw_spikes=False` to receive spikes from the output node `o1`. This will pass the neuron's spikes through the `bucket`, which decodes by thinning the pool's merged spike-train in Braindrop's accumulators. The thinned spike-train ultimately reaches the output node `o1` and is send back to the host, binned according to the resolution set by `upstream_time_res`.<jupyter_code>############################################################################
# Testing the trained decoders

# all we have to is change the decoders
decoder_conn.reassign_weights(dec.T)

# then call remap
print("remapping weights")
HAL.remap_weights()

# define input parameters
Din = 1
base_period = 1 # seconds
cycles = 4

# create the input
inp_times, inputs_sin = make_sine_stimulus(Din, base_period, cycles, fmax=fmax, input_dt=input_dt)

# package together inputs and their timestamps
# times must be sent in ns (and thus are multiplied by 1e9)
FUDGE = 2
curr_time = HAL.get_time()
times = curr_time + FUDGE * 1e9 + inp_times * 1e9
times_w_end = np.hstack((times, times[-1] + input_dt * 1e9))
input_vals = {i1 : (times, inputs_sin)}

# create the RunControl object for our defined network
rc = RunControl(HAL, net)

# send the input and collect outputs from the output node
# out_dict is a dictionary of output spikes for each output node,
# out_bin_times is a array of bin
outputs_and_bin_times, _ = rc.run_input_sweep(input_vals, get_raw_spikes=False, 
                                              end_time=times_w_end[-1], rel_time=False)

# out_bin_times is length-S, bin-size set by upstream_time_res
# out_dict is a dictionary, out_ct is a S-by-N matrix of output spike counts within S bins
out_dict, out_bin_times = outputs_and_bin_times
out_ct = out_dict[o1]

# convert spike counts to spike train
out_spikes = out_ct / input_dt

# truncate inputs to length of measured outputs
inp_times_trunc = inp_times[:spikes.shape[0]]
inputs_sin_trunc = inputs_sin[:spikes.shape[0]]<jupyter_output>remapping weights<jupyter_text>3.3 Plotting the decoded outputIn order to plot the decoded spike-train, we estimate the instantaneous spike-rate by filtering with `f_out`, which has a time constant set by `tau_out`. Then for comparison, we compute the desired function by first filtering by `f_in` to model the initial filtering by Braindrop's synapses, applying the transformation $f(x)$, then filtering by `f_out` to align with the filtered spike-train.When plotting, these are rescaled back the same range by dividing the desired function by `fmax` and the decoded function by `fmax_out`.<jupyter_code>############################################################################
# Visualizing the decoded spikes

# filter the decoded spikes to generate the output signal
fh_filt = f_out.filt(out_spikes, y0=0, dt=input_dt)

# create the desired function for comparison
x_filt = f_in.filt(inputs_sin_trunc, y0=0, dt=input_dt)
fx_filt = func(x_filt)
f_filt = f_out.filt(fx_filt, y0=0, dt=input_dt)

# rescale by fmax and fmax_out to return to consistent units
fh_filt_out = fh_filt / fmax_out
f_filt_out = f_filt / fmax

# compute rmse
rmse = np.sqrt(np.mean((f_filt_out - fh_filt_out)**2))

plt.figure()
plt.title('RMSE = %.3f' % rmse)
plt.plot(inp_times_trunc, fh_filt_out, label='decoded')
plt.plot(inp_times_trunc, f_filt_out, label='true')
plt.xlabel(r'$t$ [s]')<jupyter_output><empty_output><jupyter_text>4 Visualizing Tuning CurvesGenerating tuning curves from spike data can be done in one of two ways, depending on the way spike data is represented. We'll begin with the most straightforward method: driving our ensemble with $Q^D$ inputs uniformly covering $[-f_\mathrm{max}, +f_\mathrm{max}]$ for a fixed duration and estimating each neurons firing rate. Then we'll demonstrate how we can still extract this firing-rate information when the tuning curves are represented parametrically.   4.1 Measuring tuning curves from constant input signalsTo create this input signal, we use `make_cons_stimulus()`. This function naturally traverses a uniformly-gridded domain and holds inputs for constant durations. After generating the stimulus and creating timestamps for the input, we'll use the RunControl function `run_and_sweep()` to measure binned spike counts. Then, we'll use a utility function in `pystorm.hal.data_utils`, `bins_to_rates()`, to estimate each neuron's spike rate from the bin spike counts.<jupyter_code>############################################################################
# Measuring tuning curves and plotting

# the longer the hold time, the more accurately we estimate the spike rate
hold_time = 0.5

# stimulus parameters
Din = 1
Q = 10
fmax = 1000

# generate the input signal
inp_times, inputs_cons = make_cons_stimulus(Din, hold_time, Q, fmax=fmax)

# package together the stopwatch times and input stimuli
# times must be sent in ns (and thus are multiplied by 1e9)
FUDGE = 2
curr_time = HAL.get_time()
times = curr_time + FUDGE * 1e9 + inp_times * 1e9
times_w_end = np.hstack((times, times[-1] + hold_time * 1e9))
input_vals = {i1 : (times, inputs_cons)}

# send inputs and collect spikes
print("getting spikes")
_, spikes_and_bin_times = rc.run_input_sweep(input_vals, get_raw_spikes=True, 
                                             end_time=times_w_end[-1], rel_time=False)
print("done")

# parse spike data into bin times and binned spike counts
spike_dict, spike_bin_times = spikes_and_bin_times
spike_ct = spike_dict[p1]

# its arguments are spike counts, original bin times, bin edges, and a discard fraction
# it averages spike count over the holding time, discarding some fraction of initial spikes
# it returns a Q**Din-by-N matrix, exactly corresponding to the tuning curve matrix
A = data_utils.bins_to_rates(spike_ct, spike_bin_times, times_w_end, init_discard_frac=.2)

# plot the tuning curves
x = input_cons / fmax
plt.figure(figsize=(10,4))
plt.subplot(121)
plt.plot(x, A);
plt.xlabel(r'$x$')
plt.ylabel('Firing Rate [Hz]');

plt.subplot(122)
plt.plot(x, A)
plt.ylim([0, 2e3])
plt.xlabel(r'$x$')
plt.ylabel('Firing Rate [Hz]');<jupyter_output>getting spikes
done<jupyter_text>4.2 Estimating tuning curves from parametric dataWe'll now show you how you can estimate tuning curves when neurons' firing-rates are not explicit functions of the input $\mathbf{x}$. Instead, both inputs $\mathbf{x}$ and a $\mathbf{a}$ are functions of time $t$. So long as the input set $\mathcal{X} = \{\mathbf{x}(t) | t \leq T\}$ adequately samples the domain $[-f_\mathrm{max}, +f_\mathrm{max}]$, you can estimate the pool's firing-rate $\mathbf{a}(\mathbf{x}')$ for a new input $\mathbf{x}'$ by interpolating over $\{\mathbf{a}(\mathbf{x}_j), \mathbf{x}_j \in \mathcal{X} \cap\mathcal{N}(\mathbf{x}')\}$.In the code below, we demonstrate this principle using inputs $x(t)$ generated by `make_sine_inputs()`. To improve the accuracy of the firing-rate estimates, we slow the input oscillation with a longer `base_period` and additional cycles.<jupyter_code>def tuning_curve_from_parametric(x_times, x_vals, binned_spk_ct, curve_pts, 
                                 sigma=0.1, threshold=1e-3, default_value=None):
    """
    Estimate the tuning curve, given time-varying spike counts 
    and a time-varying signal
    
    Parameters
    ----------
    x_times: Python list
        list of times (in seconds) that we know the input for
    x_vals: numpy array
        the corresponding input for each time in x_times, in abstract units (not fmax)
    binned_spk_ct: numpy array
        spike counts for each bin
    curve_pts: numpy array
        the x values we want a frequency estimate for, in abstract units (not fmax)
    sigma: float
        std.dev of the gaussian weighting (in x-space)
    threshold: float
        if the total weight is less than this, we have no good estimate
        for this x point
    default_value: int or None
        the value to use if the total weight is less that the threshold 
        (usually 0 or None)

    Returns
    -------
    curve: numpy array
        tuning curve for the given neuron's binned spike counts
    """
    if len(curve_pts.shape) == 1:
        curve_pts = curve_pts.reshape(curve_pts.shape[0],1)
    if len(x_vals.shape) == 1:
        x_vals = x_vals.reshape(x_vals.shape[0],1)
    
    # convert binned spike counts to list of spike times
    spike_times = []
    for j, ct in enumerate(binned_spk_ct[:-1]):
        if ct > 0:
            # ct spikes within a dt are uniformly spaced within the dt
            spike_times.extend(list(np.linspace(x_times[j], x_times[j+1], ct)))
            
    # number of dimensions
    D = x_vals.shape[1]
    
    # throw away spikes that occur at the same time
    t = np.unique(spike_times)
    
    # compute frequency for each inter-spike interval
    freq = 1.0/np.diff(t)
    
    # the midpoints of each inter-spike interval
    t_mid = (t[:-1] + t[1:])/2
    
    # find the x values for each time point in t_mid
    x_mid = scipy.interpolate.griddata(x_times, x_vals, (t_mid,), method='nearest')
    
    # compute the L2 norm (squared) for each pair of points
    #  in x_mid and curve_pts.  That is, find the squared distance in 
    #  x-space between each of the points we're interested in (curve_pts)
    #  and each point we have (x_mid)
    v = np.sum([(x_mid[:,i][:,None]-curve_pts[:,i][None,:])**2 for i in range(D)], axis=0)
    
    # compute the weighting for each point (gaussian weighting)
    w = np.exp(-v/(2*sigma)**2)/(np.sqrt(2*np.pi*sigma**2))
    
    # for each point in curve_pts, determine the frequency using
    #  a weighted average (where the weight is the gaussian based on distance)
    total = np.sum(w, axis=0)
    safe_total = np.where(total<threshold, 1.0, total)
    curve = np.sum(w*freq[:,None], axis=0)/safe_total
    
    # for values where the weights are below some theshold, 
    #  set them to some default value (usually None or 0)
    curve[total<threshold] = default_value
    
    return curve
############################################################################
# Estimating tuning curves and plotting

# the longer the period and higher the cycle count, the more densely
# we sample the domain
base_period = 1
cycles = 10

# stimulus parameters
Din = 1
fmax = 1000
input_dt = 1e-3

# generate the input signal
inp_times, inputs_sin = make_sine_stimulus(Din, base_period, cycles, fmax=fmax, input_dt=input_dt)

# package together the stopwatch times and input stimuli
# times must be sent in ns (and thus are multiplied by 1e9)
FUDGE = 2
curr_time = HAL.get_time()
times = curr_time + FUDGE * 1e9 + inp_times * 1e9
times_w_end = np.hstack((times, times[-1] + input_dt * 1e9))
input_vals = {i1 : (times, inputs_sin)}

# send inputs and collect spikes
print("getting spikes")
_, spikes_and_bin_times = rc.run_input_sweep(input_vals, get_raw_spikes=True, 
                                             end_time=times_w_end[-1], rel_time=False)
print("done")

# parse spike data into bin times and binned spike counts
spike_dict, spike_bin_times = spikes_and_bin_times
spike_ct = spike_dict[p1]

# we define the x-values that we will estimate firing rates for
Q = 20
x = np.linspace(-1, 1, Q)
A = np.zeros((Q, N))

# we go through each neuron's binned spike counts and extract its tuning curve
for n in range(N):
    # since interpolation radius is not defined in units of fmax, we convert inputs_sin
    # and x back to abstract units
    A[:,n] = tuning_curve_from_parametric(inp_times, inputs_sin / fmax, spike_ct[:,n], x)

# plot the tuning curves
plt.figure(figsize=(10,4))
plt.subplot(121)
plt.plot(x, A);
plt.xlabel(r'$x$')
plt.ylabel('Firing Rate [Hz]');

plt.subplot(122)
plt.plot(x, A)
plt.ylim([0, 2e3])
plt.xlim([-1, 1])
plt.xlabel(r'$x$')
plt.ylabel('Firing Rate [Hz]');<jupyter_output>getting spikes
done
{'fraction_non_alphanumeric': 0.11869418825940566, 'fraction_numerical': 0.01402430098082272, 'mean_word_length': 4.2580049261083746, 'pattern_counts': {'":': 0, '<': 59, '<?xml version=': 0, '>': 54, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '56744', 'license': 'MIT', 'n_tokens_mistral': 11324, 'n_tokens_neox': 10543, 'n_words': 3956, 'path': 'pystorm/examples/pystorm_walkthrough.ipynb', 'repo_name': 'Stanford-BIS/pystorm', 'size': '34155'}