【Spikingjelly】SNN框架教程的代码解读_3

时间驱动：神经元

LIF神经元

clock_driven.examples.lif_fc_mnist.py在该系列文章第一篇中做过分析，里面介绍了LIF神经元的阈下动态方程和代码实现，参考链接：
【Spikingjelly】SNN框架教程的代码解读.

这一节比较简单，以时间驱动的方式逐步给与神经元输入，并查看它的膜电位和输出脉冲。

导入相关模块：

import torch
import torch.nn as nn
import numpy as np
from spikingjelly.clock_driven import neuron
from spikingjelly import visualizing
from matplotlib import pyplot as plt

a. 单个神经元

建立LIF神经元模型并且进行逐步仿真

    lif = neuron.LIFNode(tau=100.)

    lif.reset()
    x = torch.as_tensor([2.])
    T = 150
    s_list = []
    v_list = []
    for t in range(T):
        s_list.append(lif(x)) #保存每次生成的脉冲
        v_list.append(lif.v) #保存每次积分的膜电势

    visualizing.plot_one_neuron_v_s(np.asarray(v_list), np.asarray(s_list), v_threshold=lif.v_threshold,
                                    v_reset=lif.v_reset,dpi=200)
    plt.show()

其中visualizing.plot_one_neuron_v_s()用于观察此处单个神经元膜电位与输出脉冲

def plot_one_neuron_v_s(v: np.ndarray, s: np.ndarray, v_threshold=1.0, v_reset=0.0,
                        title='$V_{t}$ and $S_{t}$ of the neuron', dpi=200):
    '''
    :param v: shape=[T], 存放神经元不同时刻的电压
    :param s: shape=[T], 存放神经元不同时刻释放的脉冲
    :param v_threshold: 神经元的阈值电压
    :param v_reset: 神经元的重置电压。也可以为 ``None``
    :param title: 图的标题
    :param dpi: 绘图的dpi
    :return: 一个figure
    '''
    fig = plt.figure(dpi=dpi)
    ax0 = plt.subplot2grid((3, 1), (0, 0), rowspan=2)
    ax0.set_title(title)
    T = s.shape[0]
    t = np.arange(0, T)
    ax0.plot(t, v)
    ax0.set_xlim(-0.5, T - 0.5)
    ax0.set_ylabel('voltage')
    ax0.axhline(v_threshold, label='$V_{threshold}$', linestyle='-.', c='r')
    if v_reset is not None:
        ax0.axhline(v_reset, label='$V_{reset}$', linestyle='-.', c='g')
    ax0.legend()
    t_spike = s * t
    mask = (s == 1)  # eventplot中的数值是时间发生的时刻，因此需要用mask筛选出
    ax1 = plt.subplot2grid((3, 1), (2, 0))
    ax1.eventplot(t_spike[mask], lineoffsets=0, colors='r')
    ax1.set_xlim(-0.5, T - 0.5)

    ax1.set_xlabel('simulating step')
    ax1.set_ylabel('spike')
    ax1.set_yticks([])

    ax1.xaxis.set_major_locator(matplotlib.ticker.MaxNLocator(integer=True))
    return fig, ax0, ax1

这里$V_{reset}$为0.0，重置方式是hard reset。膜电势超过阈值，恢复到重置电压0.0V，结果如下：

b. 多个神经元

单层，32个神经元的膜电位和输出脉冲：

    lif.reset()
    x = torch.rand(size=[32]) * 4 #随机数，乘谁都可以
    T = 50
    s_list = []
    v_list = []
    for t in range(T):
        s_list.append(lif(x).unsqueeze(0))
        v_list.append(lif.v.unsqueeze(0))

    s_list = torch.cat(s_list)
    v_list = torch.cat(v_list)

    visualizing.plot_2d_heatmap(array=np.asarray(v_list), title='Membrane Potentials', xlabel='Simulating Step',
                                ylabel='Neuron Index', int_x_ticks=True, x_max=T, dpi=200)
    visualizing.plot_1d_spikes(spikes=np.asarray(s_list), title='Membrane Potentials', xlabel='Simulating Step',
                               ylabel='Neuron Index', dpi=200)
    plt.show()

这里膜电位的输入就是x，结果如下：

$$FiringRate=\frac{N}{T}$$

时间驱动：编码器

为什么要用编码器？
因为脉冲神经网络中的神经元信息传递都依靠脉冲，所以，需要将输入网络的实数值也编码成脉冲序列。
现有的脉冲编码方式主要有基于速率的编码方式（有时候也叫频率编码），如：泊松编码；还有基于时间的编码，如：延时编码

泊松编码器

同样该系列文章第一篇中做过分析event_driven/examples/tempotron_mnist.py，里面给出了泊松编码的实现方式：

class PoissonEncoder(StatelessEncoder):
    def __init__(self):
        """
        无状态的泊松编码器。输出脉冲的发放概率与输入 ``x`` 相同。
        """
        super().__init__()

    def forward(self, x: torch.Tensor):
        out_spike = torch.rand_like(x).le(x).to(x)
        return out_spike

你可能会好奇，forward中用产生随机数的方式，按照脉冲发放的概率生成脉冲，没有见到泊松编码中的泊松分布的概率密度函数：
$$P(X=k)=\frac{{\lambda }^k}{k!}{e}^{-\lambda },k=0,1,\cdots$$
首先理解泊松分布，可参考以下链接：泊松分布的现实意义是什么，为什么现实生活多数服从于泊松分布？

1. 事实上，泊松过程又被称为泊松流，当一个脉冲流满足独立增量性、增量平稳性和普通性时，这样的脉冲流就是一个泊松流。泊松编码器将输入数据 x 编码为发放次数分布符合泊松过程的脉冲序列。
2. 泊松分布的均值$\lambda$在脉冲序列中，意义为时间T时间步内平均发放的脉冲数N，如果N/T也就是我们常说的平均脉冲发放率
3. 上述公式，描述的是T时间内，发放k个脉冲的概率。而我们产生脉冲序列时，产生的是一个时间步，这个时间步是否发放脉冲。代码实现中forward产生的是一个时间步的脉冲，这个forward执行人为选的总时间步T次。不同时间步出现脉冲的个数是相互独立的，为此每个时间步内是否发放脉冲都有一定的概率，这个概率和图像的其与二维图像的像素值成正比
4. 总结：泊松编码将28*28的输入图像，编码为长度为28*28的0，1脉冲序列，每个时间步是否发脉冲都有自己的概率，相互独立，这个概率正比与输入图像的值，概率值在[0，1]。总体来看，T时间步内平均发放脉冲数为$\lambda$

泊松分布适合于描述单位时间内随机事件发生的次数，这正好时间T内发放几个脉冲对应，用于编码产生我们的脉冲序列。

a. 单独的时间步长

输入图像lena512.bmp，仿真20个时间步，得到20个脉冲矩阵。观察泊松编码的效果

    import torch
    import numpy as np
    import matplotlib
    import matplotlib.pyplot as plt
    from PIL import Image
    from spikingjelly.clock_driven import encoding
    from spikingjelly import visualizing

    # 读入lena图像
    lena_img = np.array(Image.open('lena512.bmp')) / 255
    x = torch.from_numpy(lena_img)

    pe = encoding.PoissonEncoder()
    w, h = x.shape
    encode_method = 'Possion' #不同编码器的演示教程
    if encode_method == 'Possion':
        single_timestep = False
        sum_timestep = True
        if single_timestep:
            # 仿真20个时间步长，将图像编码为脉冲矩阵并输出
            out_spike = torch.full((20, w, h), 0, dtype=torch.bool)
            T = 20
            for t in range(T):
                out_spike[t] = pe(x)

            plt.figure()
            plt.imshow(x, cmap='gray')
            plt.axis('off')

            visualizing.plot_2d_spiking_feature_map(out_spike.float().numpy(), 4, 5, 30, 'PoissonEncoder')
            plt.axis('off')
            plt.show()

b. 多个时间步长叠加

同样对lena灰度图进行编码，仿真512个时间步长，将每一步得到的脉冲矩阵叠加，得到第1、128、256、384、512步叠加得到的结果并画图：

        if sum_timestep:
            # 仿真512个时间步长，将编码的脉冲矩阵逐次叠加，得到第1、128、256、384、512次叠加的结果并输出
            superposition = torch.full((w, h), 0, dtype=torch.float)
            superposition_ = torch.full((5, w, h), 0, dtype=torch.float)
            T = 512
            for t in range(T):
                superposition += pe(x).float()
                if t == 0 or t == 127 or t == 255 or t == 387 or t == 511:
                    superposition_[int((t + 1) / 128)] = superposition

            # 归一化
            for i in range(5):
                min_ = superposition_[i].min()
                max_ = superposition_[i].max()
                superposition_[i] = (superposition_[i] - min_) / (max_ - min_)

            # 画图
            visualizing.plot_2d_spiking_feature_map(superposition_.numpy(), 1, 5, 30, 'PoissonEncoder')
            plt.axis('off')

            plt.show()

结果如下，可见当仿真足够的步长，泊松编码器得到的脉冲叠加后几乎可以重构出原始图像。

周期编码器

周期编码器是有状态编码器，首次forward时使用encode函数对 T 个时刻的输入序列x 进行编码得到 spike，在第 t 次调用 forward 时会输出 spike[t % T]。
有状态编码器实现如下：

class StatefulEncoder(base.MemoryModule):
    def __init__(self, T: int):
        """
        :param T: 编码周期。通常情况下，与SNN的仿真周期（总步长一致）
        :type T: int
        有状态编码器的基类。有状态编码器 ``encoder = StatefulEncoder(T)``，编码器会在首次调用 ``encoder(x)`` 时对 ``x` 进行编码。在
        第 ``t`` 次调用 ``encoder(x)`` 时会输出 ``spike[t % T]``
        """
        super().__init__()
        assert isinstance(T, int) and T >= 1
        self.T = T
        self.register_memory('spike', None)
        self.register_memory('t', 0)

    def forward(self, x: torch.Tensor = None):
        if self.spike is None:
            self.encode(x)

        t = self.t
        self.t += 1
        if self.t >= self.T:
            self.t = 0
        return self.spike[t]

    @abstractmethod
    def encode(self, x: torch.Tensor):
        raise NotImplementedError

    def extra_repr(self) -> str:
        return f'T={self.T}'

class PeriodicEncoder(StatefulEncoder):
    def __init__(self, spike: torch.Tensor):
        """
        周期性编码器，在第 ``t`` 次调用时输出 ``spike[t % T]``，其中 ``T = spike.shape[0]``
        """
        super().__init__(spike.shape[0])
        self.encode(spike)

    def encode(self, spike: torch.Tensor):
        self.spike = spike
        self.T = spike.shape[0]

延时编码器和带权相位编码器

原文教程中说的很详细，此处不再展开。

参考

原文教程：事件驱动

SNN系列｜编码篇(1)频率编码

泊松分布的现实意义是什么，为什么现实生活多数服从于泊松分布？

一种应用于脉冲神经网络的输入脉冲编码方法