Pybind11类型转换

Pybind11帮助我们方便的实现C++和Python之间的调用，无论是相互的接口暴露还是数据传递都要求我们正确的在C++和Python类型之间进行转换（官方文档）。

三种场景

1. C++原生类型暴露给python
2. Python原生类型暴露给C++
3. C++和Python原生类型间相互转换

C++原生类型暴露给Python

使用py::class_将自定义的C++类型暴露给Python（参考这里）。C++类型传递给Python时会在原生C++类型外加一层wrapper，从Python取回时只用去掉wrapper即可。

Python原生类型暴露给C++

在C++中取用Python的原生类型（例如：tuple、list），例如：

1
2
3
4
5
6

// C++中使用Python对象
void print_list(py::list my_list) {
  for (auto item : my_list) {
    std::cout << item << " ";
  }
}

Python中传入Python对象给C++:

1
2

>>> print_list([1, 2, 3])
1 2 3

这个例子里Python中的list类型没有进行任何转换，只是被封装进了C++的py::list类中。目前Pybind11支持以下类型(参考这里)：handle, object, bool_, int_, float_, str, bytes, tuple, list, dict, slice, none, capsule, iterable, iterator, function, buffer, array, array_t。

C++和Python原生类型间相互转换

有些场景C++和Python都用的是各自的原生类型，Pybind11支持常见C++原生类型和Python原生类型间的相互转换，例如：

1
2
3
4
5

void print_vector(const std::vector<int> &v) {
  for (auto item : v) {
    std::cout << item << "\n";
  }
}

Python中调用：

1
2

>>> print_vector([1, 2, 3])
1 2 3

容易注意到Python中传入的为list类型，而C++中处理的类型为const std::vector<int> &。Pybind11这种默认类型转换会在原生类型间拷贝数据，例如上面Python中的调用会首先将list中的数据拷贝到一个std::vector<int>中，然后进行C++中的调用。

NOTE: 默认的拷贝可能开销很大，可以通过手写opaque types重载类型转换的wrapper来避免不必要的拷贝。opaque types还没有使用过，有机会单独挖坑写。

STL库

Pybind11默认（pybind11/pybind11.h）支持std::pair<>，std::tuple<>和list, set, dict之间的自动转换。引入pybind11/stl.h可以新增std::vector<>, std::deque<>, std::list<>, std::array<>, std::set<>, std::unordered_set<>, std::map<>和std::unordered_map<>的与Python的自动转换。

绑定STL容器

pybind11/stl_bind.h可以帮助我们将STL容器作为原生对象暴露出来。例如：

1
2
3
4
5
6
7
8
9
10
11
12

#include <pybind11/stl_bind.h>

PYBIND11_MAKE_OPAQUE(std::vector<int>);
PYBIND11_MAKE_OPAQUE(std::map<std::string, double>);

// ...

// 绑定
py::bind_vector<std::vector<int>>(m, "VectorInt");
py::bind_map<std::map<std::string, double>>(m, "MapStringDouble");
// 或者同时设置绑定的作用域
py::bind_vector<std::vector<int>>(m, "VectorInt", py::module_local(false));

Bonfire

Cannot move my eyes away from the bonfire

It was his will

It is a dance of life and death

It would be ashes of bodies

It is not an illusion as long as it burns

But I had stared at it for too long

So I lit my hand

And walked into the dark

As my lights going dimmer and dimmer

How I wish it is a cool autumn

I would just fall down on the ground

And flare up the whole forest

TF graph_matcher用法及源码探究

graph_matcher是tensorflow.contrib中量化模块（quantize）的一个子模块，用于在计算图中描述和匹配特定的模式。配合对模式匹配后的处理，可以在python层面实现计算图的pass。

用法: 以Conv + BN融合为例

Conv+BN是CNN网络中常见的组合。如果我们观察两者的计算公式，可以发现两者可以融合为一个算子达到运行加速的效果（实质上，Conv计算上等价于MatMul，所以融合也适用于FC+BN等组合）。

融合原理

首先，分别观察Conv和BN的计算公式：

​ Conv的计算公式：$z = w * x + b$

​ BN的计算公式：$y = \frac{(z - \mu_B) * \gamma}{ \sigma_B} + \beta$

实际上，我们可以通过更新Conv的weight和bias直接在Conv中完成Conv + BN所需完成的计算。略去推导，直接给出新的weight和bias的计算公式如下：

​ $w^{\prime} = \frac{w * \gamma}{\sigma_B}$, $b^{\prime} = \frac{(b - \mu_B)\gamma}{\sigma_B} + \beta$

代入新的$w^{‘}$和$b^{‘}$，容易验证新的Conv计算等价于Conv + BN:

graph_matcher实现

tensorflow.contrib.quantize.python中包含了Conv + BN融合的实现，完整的代码较长，我们重点关注其中对Conv+BN模式描述的部分。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61

def _FindFusedBatchNorms(graph):
    """Finds all ops and tensors related to found FusedBatchNorms.
    	Args:
    		graph: Graph to inspect.
    	Returns:
    		_FusedBatchNormMatches.
    """
    input_pattern = graph_matcher.OpTypePattern('*')
    weight_pattern = graph_matcher.OpTypePattern('*')
    gamma_pattern = graph_matcher.OpTypePattern('*')
    beta_pattern = graph_matcher.OpTypePattern('*')
    mean_pattern = graph_matcher.OpTypePattern('*')
    variance_pattern = graph_matcher.OpTypePattern('*')
    
    moving_average_pattern = graph_matcher.OpTypePattern('*')
    bn_decay_pattern = graph_matcher.OpTypePattern('*')
    layer_pattern = graph_matcher.OpTypePattern('Conv2D|DepthwiseConv2dNative|MatMul', inputs=[input_pattern, weight_pattern])
    ...
    layer_output_pattern = graph_matcher.OneofPattern([layer_pattern_with_identity, layer_pattern, batch_to_space_pattern])
    ...
    bn_matcher = graph_matcher.GraphMatcher(graph_matcher.OneofPattern([matmul_bn_output_reshape_pattern, batch_norm_pattern]))
    ...
    def _GetLayerMatch(match_result):
        """Populates a layer match object containing ops/tensors for folding BNs.
        Args:
        	match_result: Matched result from graph matcher
        Returns:
        	layer_op: Matching conv/fc op prior to batch norm
        	BatchNormMatch: _BatchNormMatch containing all required batch norm parameters.
        """
        ...
    layer_matches = []
    matched_layer_set = set()
    for match_result in bn_identity_matcher.match_graph(graph):
        layer_op, layer_match = _GetLayerMatch(match_result)
        if layer_op is not None:
            if layer_op not in matched_layer_set:
                matched_layer_set.add(layer_op)
                layer_matches.append(layer_match)
    ...
    return layer_matches
def _FoldFusedBatchNorms(graph, is_training, freeze_batch_norm_delay):
    """Finds fused batch norm layers and folds them into preceding layers.
    Folding only affects the following layers: Conv2D, fully connected, depthwise convolution.
    Args:
    	graph: Graph to walk and modify.
    	is_training: Bool, true if training.
    	freeze_batch_norm_delay: How many steps to wait before freezing moving mean
      and variance and using them for batch normalization.
  Raises:
  	ValueError: When batch norm folding fails.
  	"""
    for match in _FindFusedBatchNorms(graph):
        scope, sep, _ = match.layer_op.name.rpartition('/')
        with graph.as_default(), graph.name_scope(scope + sep):
            with graph.name_scope(scope + sep + 'BatchNorm_Fold' + sep):
                # new weights = old weights * gamma / sqrt(variance + epsilon)
                # new biases = -mean * gamma / sqrt(variance + epsilon) + beta
                multiplier_tensor = match.gamma_tensor * math_ops.rsqrt(match.variance_tensor + match.bn_op.get_attr('epsilon'))
                bias_tensor = math_ops.subtract(match.beta_tensor, match.mean_tensor * multiplier_tensor, name='bias')
    ...

流程

使用GraphMatcher进行pattern描述、匹配、替换主要分为以下几步：

1. 使用OpTypePattern及OneofPattern（语法糖）自底向上构建目标pattern
2. 用目标pattern构造GraphMatcher
3. 向GraphMatcher的match_graph方法传入要匹配的图，获得match_result
4. 从match_result中取出需要复用的pattern中的节点，构造新的节点替换pattern

graph_matcher实现

graph_matcher的实现主要包括三个部分：Pattern, GraphMatcher, MatchResult。

Pattern实现

Pattern作为一个抽象类，要求子类必须实现match方法。match方法接收两个参数：op；tensor。

Pattern类有两个子类：OpTypePattern类可以限定节点的类型、输入，可以描述一个类型树；Oneof

Pattern作为语法糖用于描述one-of关系，也就是匹配输入多个子模式之一即可。

NOTE: 当前实现了的模式中，match方法中的tensor只是占位用，没有实际使用到。

OpTypePattern

构造函数（def __init__(self, op_type, name=None, inputs=None, ordered_inputs=True)）通过限定节点的类型、输入来描述一个类型树；对应的，match中也会递归地对输入节点调用match函数。NOTE: OpType的匹配是使用字符串来完成的。

MatchResult

保存match的结果，可以从Pattern实例映射到对应的匹配到的op和tensor。

GraphMatcher

GraphMatcher中会保存一个pattern，提供方法来检验输入的op或者graph是否和pattern匹配，主要方法有：

1. match_op
2. match_ops
3. match_graph

拓展思考

TensorFlow中能够轻松的在python中操作图主要得益于图数据结构对python的暴露。当前MindSpore要在python中支持图pass(图中模式的匹配和替换)，可以对比两种思路：

1. python向C++注册pass，python中对模式和要替换的目标进行描述，C++中运行pass
   1. 优点：可以复用部分优化器部分的代码；执行效率较高
   2. 缺点：python pass中的pattern、target与C++通信较复杂
2. C++向python暴露图接口，直接在python中完成改图
   1. 优点：对图修改的逻辑全部包含在python中
   2. 缺点：效率较低，但此类任务通常较低频，性能要求不高

综合考虑，C++新增向python暴露图接口，直接在python中完成改图较合理。

SNN学习笔记5：LCA

稀疏编码

实验表明人脑对于外界刺激采取一种稀疏的内在表示，例如自然图像只需要用稀疏词典中一个很小的子集及合适的对应系数来进行稀疏近似（sparse approximation）。

稀疏近似

稀疏近似的数学表述如下：

​ 给定一个N维刺激$s \in \mathbb{R}^N$, 找到一个基于由M个向量$\{\phi_m\}$组成的词典$D$的表示。当词典$D$是overcomplete时（i.e M > N），我们可以有无穷多种方式来选取词典中向量对应的稀疏$\{a_m\}$ 来表示$s$: $s = \sum_{m = 1}^Ma_m\phi_m$。

在最优稀疏近似中，我们希望尽可能少的使用D中的向量，也就是系数不为0的向量尽可能的少：

上式中$||a||_0$表示$\ell^0$ norm, 也就是 $a = [a_1, a_2, ..., a_M]$中非零元素的个数。需要注意的是，这个组合优化问题是NP-hard的。

Basis Pursuit目标

对$\ell_0$的优化是NP-hard的，BP目标函数尝试将优化目标改为最小化系数向量的$\ell_1$ norm：

BP目标函数在信号$s$相对稀疏时也可以得到最优稀疏近似。

BPDN：重建误差

实际操作中，由于$s$中存在噪音，我们不应该要求完美重建。BPDN（Basis Pursuit De-Noising）目标函数在BP的基础上引入了MSE重建误差来平衡正则项与重建精度：

公式中的$\lambda$正是用来权衡重建误差与正则项的。

MP算法

在信号处理社区，常用MP(Matching Pursuit)算法来求解BPDN。MPs算法本质上是一种贪心算法，流程如下：

1. 将残差初始化为$s$: $r_0 = s$
2. 在第k次迭代，通过$\theta_k = argmax_m|\langle r_{k-1}, \phi_m\rangle|$找到词典M中的索引$\theta_k$
3. 更新残差：$r_k = r_{k - 1} - \phi_{\theta_k}d_k$

K次迭代后得到一个$s$的稀疏近似：$\hat{s} = \sum_{k = 1}^K \phi_{\theta_k}d_k$。

LCA

LCA(Locally Competitive Algorithm)是一种稀疏编码算法，相比MP算法，不仅考虑到了选取最稀疏表示的目标，也考虑了选取最能表征信号特性的向量的目标。同时，LCA对随时间变化的信号的处理进行了优化，LCA不用每一步都从头进行稀疏近似，而是基于上一步的表征向量进行更新。

框架

LCA中，词典中的每个向量$\phi_m \in D$都被关联到一个神经元。神经元中维护自己的膜电位$u_m(t)$，神经元的输入电流为输入与神经元的感受野的匹配度：$b_m(t) = \langle\phi_m, s(t)\rangle$。当神经元m的膜电位超过阈值$\lambda$时，输出一个激活信号$a_m = T_{\lambda}(u_m)$并向周边神经元n发射抑制信号$a_m G_{m, n}$其中$G_{m, n} = \langle \phi_m, \phi_n \rangle$。

NOTE: 从抑制信号的公式中可以看出，一个神经元的激活越强，对周边神经元的抑制越强；一个神经元与周围的神经元越相似，对周围神经元的抑制越强。这种机制会导致匹配度最高的神经元得到最大的电流输入，然后抑制周边神经元得到输入及进行反向抑制，以此达到获取稀疏表示的效果（WTA：winner takes all）。

上面的膜电位变化机制可以用下面的常微分方程来描述：

Demo实现

github上的一个参考实现：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56

import nengo
import numpy as np
import matplotlib.pyplot as plt

def LCA(d, n_neurons, dt):
    k = 1.
    beta = 1.
    tau_model = 0.1
    tau_actual = 0.1

    I = np.eye(d)
    inhibit = 1 - I
    B = 1. / tau_model
    A = (-k * I - beta * inhibit) / tau_model

    with nengo.Network(label="LCA") as net:
        net.input = nengo.Node(size_in=d)
        # array of ensembles: d ensembles, each with n_neurons neurons
        x = nengo.networks.EnsembleArray(
            n_neurons, d,
            eval_points=nengo.dists.Uniform(0., 1.),
            intercepts=nengo.dists.Uniform(0., 1.),
            encoders=nengo.dists.Choice([[1.]]),
            label="state")
        # transform: linear transformation mapping the pre output to the post input
        # synapse: synapse model for filtering
        nengo.Connection(x.output, x.input, transform=tau_actual * A + I, synapse=tau_actual)
        nengo.Connection(net.input, x.input, transform=tau_actual*B, synapse=tau_actual)
        net.output = x.output
    return net

def main():
    dt = 0.001
    with nengo.Network(seed=42) as model:
        # winner takes all
        wta = LCA(3, 200, dt)
        stimulus = nengo.Node([0.8, 0.7, 0.6])
        nengo.Connection(stimulus, wta.input, synapse=True)

        p_stimulus = nengo.Probe(stimulus, synapse=None)
        p_output = nengo.Probe(wta.output, synapse=0.01)
    with nengo.Simulator(model, dt=dt) as sim:
        sim.run(1.)

    fig = plt.figure()
    plt.plot(sim.trange(), sim.data[p_output])
    plt.title("(a)LCA")
    plt.xlabel("Time [s]")
    plt.ylabel("Decoded output")
    plt.locator_params(axis='y', nbins=5)
    plt.tight_layout()
    plt.show()
        

if __name__ == "__main__":
    main()

效果如下：

可以看到3个ensemble中只有得到最强输入的一个保留了下来。

SNN学习笔记4： Nengo

Nengo 是一个用于神经建模框架，其扩展NengoDL支持混用包含了生物细节的神经模型和现在流行的深度学习框架（例如：TensorFlow）。

安装

nengo安装命令：pip install nengo nengo-gui

测试是否安装成功可以尝试运行nengo-gui界面：$: nengo

nengo-dl安装命令：

1. 安装依赖的tensorflow：conda install tensorflow
2. 安装nengo-dl: pip install nengo-dl

使用

nengo有两种使用模式：GUI和Python解释器。Python解释器模式下，nengo的表现就是一个普通的Python库，因此下面仅介绍GUI模式的使用方式及限制。

GUI模式

直接在命令行运行$:nengo即可在网页中运行图形界面。此时左侧会展示当前图的结构（方形表示Nodes，圆形表示Ensembles, 圆角矩形表示Networks），右侧展示对应的代码。

可以点击左上角的文件夹图标可以运行很多内建的例子，例如/built-in examples/tutorial/15-lorenz.py可以运行一个洛伦兹吸引子的例子，效果很酷。

需要注意的是，GUI模式下代码有如下限制：

1. 顶层网络必须叫model
2. 不能构建Simulator对象
3. 不能使用Matplotlib绘图

核心概念

上图Nengo Core主要包含五个核心Nengo对象和一个基于Numpy的模拟器。五个对象如下：

1. nengo.Network： 一个网络可以包含ensembles、nodes、connections和其它网络
2. nengo.Ensemble：一组神经元，用于表征一个向量
   1. nengo.ensemble.Neurons: 用于连接ensemble中特定神经元的接口
3. Node：用于提供输入以及处理输出
4. Connection：连接两个对象 NOTE: 和TensorFlow等不同，连接作为一个独立的对象，方便进行独立的设置
   1. nengo.connection.LearningRule：为连接制定学习规则
5. Probe：用于将对象的数据在模拟器运行时取出

Nengo-DL Demo: Relu VS Spiking neurons

下面的例子会展示Relu作为激活神经元和脉冲神经元的差异。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46

import matplotlib.pyplot as plt
import nengo
from nengo.utils.matplotlib import rasterplot
import numpy as np
import nengo_dl

with nengo.Network() as net:
    # 输入节点，周期为1s的正弦波
    a = nengo.Node(lambda t: np.sin(2 * np.pi * t))

    # rate神经元，功能和Relu一致
    b_rate = nengo.Ensemble(10, 1, neuron_type=nengo.RectifiedLinear(), seed=2)
    nengo.Connection(a, b_rate)

    # spiking神经元
    b_spike = nengo.Ensemble(10, 1, neuron_type=nengo.SpikingRectifiedLinear(), seed=2)
    nengo.Connection(a, b_spike)

    # 模拟时取出输入输出数据
    p_a = nengo.Probe(a)
    p_rate = nengo.Probe(b_rate.neurons)
    p_spike = nengo.Probe(b_spike.neurons)

with nengo_dl.Simulator(net) as sim:
    # 运行模拟1S，上述Probe数据会被保存在sim.data字典中
    sim.run_steps(1000)

plt.figure()
plt.plot(sim.trange(), sim.data[p_a])
plt.xlabel("time")
plt.ylabel("input value")
plt.title("a")

plt.figure()
plt.plot(sim.trange(), sim.data[p_rate])
plt.xlabel("time")
plt.ylabel("firing rate")
plt.title("b_rate")

plt.figure()
# 时间栅格图，Spiking神经元的的脉冲事件出现与否以栅格形式展示
rasterplot(sim.trange(), sim.data[p_spike])
plt.xlabel("time")
plt.ylabel("neuron")
plt.title("b_spike")
plt.show()

上面代码运行的结果如下：

可以看到每个神经元的初始连接权重及bias不同，因此对输入信号的相应略有不同；spiking neurons会在电压超过0时产生脉冲发射事件，注意图二和图三中的颜色对应，我们还可以观察到电压值越高，对应的脉冲发射频率越高。

SNN学习笔记3：关联学习

关联学习（Associative Learning）认为想法（神经元激活模式）和想法、想法和经历（外界刺激）会互相关联、互相强化。

关联在学习和认知中有着核心的地位，记忆实际上就是一种关联，很多认知功能的本质就是一层或多层的关联。

赫布理论

突触（synaptic）用于连接多个神经元，其连接强度具有可以调整，这种可以调整的特性叫做突触可塑性（synaptic plasticity）。突触可塑性是记忆和学习的基础，赫布理论认为突触前的神经元对突触后神经元的反复刺激可以增加突触的传递效能，也就是强化了这两个神经元之间的连接。

赫布理论强调细胞A的激活导致了细胞B的激活以及两者之间连接的强化，这种先后/因果关系的理论也叫STDP(spike-timing-dependent plasticity)。

STDP

STDP基于一个神经元的输入脉冲及输出脉冲间的相对时间来调整神经元间连接的强度。如果一个神经元的输入倾向于刚好出现在输出脉冲前，那么此神经元与输入神经元之间的连接倾向加强；如果一个神经元的输入倾向于刚好出现在输出脉冲之后，那么此神经元和输入神经元之间的连接倾向减弱。

通过上面的调整机制，可能是当前神经元激活的原因的输入会被强化；不是当前神经元激活原因的输入会被弱化。这个调整的过程最后会收敛到只有一部分连接保留下来，另一部分的连接强度降低到0，以此达到一种稀疏连接的效果。

SNN学习笔记2：神经编码

大脑中的感知细胞在受到光、声音等外界刺激时，其动作电位的激活序列会呈现出一定的时序模式。神经编码认为大脑中的感知、认知等信息由神经元的激活表征，这些激活的各方面特征（时序 、强度等）不仅能够编码数字信号也能编码模拟信号。

一些概念

ISI

ISI(interspike intervals)，激活间间隙，表示两次激活中间间隔的时间的长度。虽然每次激活的持续时间、幅度和形状可能都有差异，但通常被处理为出现/不出现的点事件（point events）。

神经编码

神经编码（Neural Encoding）将外界输入刺激映射到神经元的反应，主要关注的是理解神经元如何对刺激作出反应，并建模来尝试预测神经元对其它刺激的反应。

神经解码

神经解码（Neural Decoding）关注的是编码的反向映射，也就是通过观察神经元的活动推导出对应的外界刺激。

一些神经编码的理论

一系列的神经脉冲中包含了丰富的信息，不同编码理论侧重有所不同，在信息精度/浓度上也各有取舍。对不同功能的神经元，适用的编码也会不同。例如，对于控制肌肉收缩的运动细胞，基本只关心脉冲的发射频率（firing rate）。而对于处理复杂认知任务（例如，视觉、听觉）的神经细胞，每个脉冲出现的精确时间都包含了信息。

频率编码（Rate coding）

频率编码模型将外界刺激的强度编码为脉冲的发射频率（firing rate），也就是外界的刺激越强，对应神经元的脉冲发射频率越高（通常是非线性变化）。频率编码假设外界刺激的绝大多数信息都包含在神经元的发射频率中，这是一种早期的编码方案，实验表明脉冲准确的时间中也包含了大量信息。

NOTE：目前发射频率没有一个公认的定义，常见的定义有：1）随时间平均 2）多次实验平均。

频率编码：脉冲计数码率

脉冲计数码率（spike-count rate, a.k.a temporal average）由对一次实验中的脉冲数目进行计数，然后除以实验持续的时间得到。显然，这种平均只适合恒定或者变化较慢的外界刺激，对于变化较快的外界刺激意义不大。

频率编码：时间相关发射率

时间相关发射率（time-dependent firing rate）定义为$\frac{C_{t}^{t + \Delta t}spikes}{\Delta t}$,其中$C_t^{t + \Delta t}spikes$表示$t$到$\Delta t$之间的脉冲计数。和脉冲计数码率不同，时间相关发射率不仅可以处理常量刺激，也可以处理时间相关的刺激。

NOTE: 时间相关发射率依赖于有多个独立神经元，每个神经元接收同一种刺激的假设。

时间编码（Temporal coding）

与频率编码不同，在时间编码中，脉冲出现的准确时间或者发射频率的波动被认为是携带信息的。

时间编码：二元编码

用二元符号来标记单位时间内是否有脉冲，1表示有脉冲，0表示没有。利用二元编码，我们可以区分频率编码中无法区分的序列，例如：0001110011和1110001100两者虽然频率一致，但是脉冲的时序显然不一样。

时间编码：ISI

ISI利用激活间的区间长度来编码激活序列。

时间编码：稀疏编码

对神经元的每次强激活单独编码，通常用线性生成模型来描述（Linear Generative Model）。

SNN学习笔记1：LIF

什么是SNN？

SNN（Spiking neural network，脉冲神经网络）号称是第三代神经网络, 与当前流行的神经网络的主要区别是将神经脉冲传播的动态过程纳入学习和推理中。SNN中的神经元不会在每次传播中都激活，而是只有当神经元的膜电位超过阈值时才激活，在激活时发射电脉冲（spikes，神经科学中常称作动作电位），这些电脉冲通过轴突传递给其它的神经元。

因此， 神经元的膜电位的变化的描述变得十分重要。下面介绍几种描述膜电位变化的模型。

Integrate-and-fire

Integrate-and-fire是最早用来描述膜电位变化的模型（1907年！），神经元膜电位的变化由下面的公式表示:

$C_m$ 表示神经元的电容, 我们容易看出上式只是电容公式Q=CV两边同时对时间求导。

此模型有如下性质：当有电流输入时，膜电位将会升高。

不应期

通过增加不应期（refractory period）$t_{ref}$可以使Integrate-and-fire模型更加精准。所谓不应期，就是在此期间神经元无法激活，这一现象在实际的神经元中已经被观察到，微观上可以由钠钾离子通道的状态来解释。

不应期限制了神经元激活的频率，激活频率与不应期$t_{ref}$的关系如下：

缺点

一个显著的缺点是，上述模型没有实现时间相关记忆（time-dependent memory）,即如果这个模型收到一个远超阈值的信号时，会永远将此信号记在自己的膜电位中。

LIF

LIF(Leaky integrate-and-fire)通过增加“漏电（Leaky）”项解决了Integrate-and-fire模型中缺少时间相关记忆的问题。模型的公式如下：

“漏电”项$\frac{V_m(t)}{R_m}$中的$R_m$表示膜电阻。

公式蕴含了下面这些有趣的性质：

1. 激活神经元的输入电流必须超过$I_{th} = \frac{V_{th}}{R_m}$, 否则漏电项会导致膜电位泄露
2. 大电流输入时，模型趋近于Integrate-and-fire+不应期

LIF: Nengo实现

Nengo是一个SNN相关的库，利用Nengo可以很容易实现上述LIF模型。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41

import matplotlib.pyplot as plt
import numpy as np
import nengo
from nengo.utils.matplotlib import rasterplot

from nengo.dists import Uniform

model = nengo.Network(label="A Single Neuron")
with model:
    neuron = nengo.Ensemble(1, dimensions=1, intercepts=Uniform(-.5, -.5), max_rates=Uniform(100, 100), encoders=[[1]])
    # input node
    cos = nengo.Node(lambda t: np.cos(8*t))
    nengo.Connection(cos, neuron)
    cos_probe = nengo.Probe(cos)
    # 神经元原始脉冲输出
    spikes = nengo.Probe(neuron.neurons)
    # 细胞体电压
    voltage = nengo.Probe(neuron.neurons, "voltage")
    # spikes filtered by a 10ms post-synaptic filter
    filtered = nengo.Probe(neuron, synapse=0.01)

with nengo.Simulator(model) as sim:
    sim.run(1)
    # plot the decoded output of the ensemble
    plt.figure()
    plt.plot(sim.trange(), sim.data[filtered])
    plt.plot(sim.trange(), sim.data[cos_probe])
    plt.xlim(0, 1)

    # plot the spiking output of the ensemble
    plt.figure(figsize=(10, 8))
    plt.subplot(221)
    rasterplot(sim.trange(), sim.data[spikes])
    plt.ylabel("Neuron")
    plt.xlim(0, 1)

    # plot the soma voltages of the neurons
    plt.subplot(222)
    plt.plot(sim.trange(), sim.data[voltage][:,0], 'r')
    plt.xlim(0, 1)
    plt.show()

右值引用和emplace

这篇文章初衷是好奇push_back和emplace_back的区别，了解之后发现绕不开右值引用，在此一并记录一下。

右值引用

右值引用（rvalue reference）是C++11中引入的新特性，显然是与C++11之前普通左值引用相对的一个概念。下面的右值引用的介绍很多参考自这篇文章。

左值与右值

一个粗略的定义如下：

左值是一个可以出现在赋值符号（=）左边或者右边的表达式, 可以理解为对一块内存的引用；

右值是一个只能出现在赋值符号右边的表达式， 注意右值不是对内存的引用，因此不能进行取地址操作。

如：

1
2
3
4
5
6
7
8
9
10
11

int a = 42;
int b = 43;
// a, b 均为左值
a = b;
b = a;
// a + b 为右值
int c = a  + b;
a + b = 42; // error!
// 不能对右值取地址
a  = foo(); // ok, foo() is a rvalue
int* p = &foo(); // invalid, 不能对右值取地址

为什么要右值引用？

move语义

假设类X中包括一个指向资源的指针m_pResource，我们想实现一个接收临时对象作为参数的拷贝赋值操作符，其实现可能如下：

1
2
3
4
5

X& X::operator=(const X &rhs) {
  //1. 拷贝rhs.m_pResource
  //2. 析构rhs.m_pResource指向的资源
  //3. 将拷贝的资源赋给self.m_pResource
}

不难看出，上面对资源m_pResource的拷贝和析构十分低效，我们可以直接和临时实例交换指针（此所谓move语义）。另外，对于非临时对象的拷贝，我们可能不想析构其资源。因此我们需要一个类型来标识这样的临时对象，对它进行单独的处理，例如：

1
2
3

X& X::operator=( <Desired type>rhs) {
  //1. swap this->m_pResource and rhs.m_pResource
}

其实右值引用就是我们想要的类型，上面代码中的<Desired type>我们可以用X&&替代，表示是X的右值引用。

emplace_back or push_back?

一个小实验

猜猜看下面的代码会进行几次复制操作？

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

#include <iostream>
#include <vector>

struct Point {
  float x, y;
  Point(float x, float y):x(x), y(y){}

  Point(const Point& point):x(point.x), y(point.y) {
    std::cout << "copying!" << std::endl;
  }
};

int main() {
  std::vector<Point> points;
  points.push_back(Point(1, 2));
  points.push_back(Point(3, 4));
  points.push_back(Point(5, 6));
  return 0;
}

答案是六次，其中三次是将临时Point对象拷贝至容器，有一次是容器容量从1到2的拷贝，有两次是容器容量从2到4的拷贝。

消除扩容拷贝

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

#include <iostream>
#include <vector>

struct Point {
  float x, y;
  Point(float x, float y):x(x), y(y){}

  Point(const Point& point):x(point.x), y(point.y) {
    std::cout << "copying!" << std::endl;
  }
};

int main() {
  std::vector<Point> points;
  // 预分配内存 NOTE: 和std::vector<Point> points(3)的区别，reserve只分配内存不调用构造函数
  points.reserve(3);
  points.push_back(Point(1, 2));
  points.push_back(Point(3, 4));
  points.push_back(Point(5, 6));
  return 0;
}

现在只用进行三次插入容器时的拷贝

move语义消除插入容器时的拷贝

直接将参数forward给容器，直接在容器中构造。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

#include <iostream>
#include <vector>

struct Point {
  float x, y;
  Point(float x, float y):x(x), y(y){}

  Point(const Point& point):x(point.x), y(point.y) {
    std::cout << "copying!" << std::endl;
  }
};

int main() {
  std::vector<Point> points;
  points.emplace_back(1, 2);
  points.emplace_back(3, 4);
  points.emplace_back(5, 6);
  return 0;
}

现在拷贝的次数为0次！

接口

首先看看C++11中两者的接口：

1
2
3
4
5

template< class... Args >
void emplace_back( Args&&... args );

void std::vector<T, Allocator>::push_back( const T& value );
void std::vector<T, Allocator>::push_back( T&& value );

注意到接收右值的接口的差别，push_back只能接收Vector中存储类型的右值作为参数，而emplace可以接收变长模板作为参数，尝试为变长模板找到最合适的构造函数直接在容器中构建。

NOTE：实验中的例子，emplace_back的变长参数也可以接收Point对象，此时会先调用Point的默认构造函数然后调用复制构造函数。对应代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

#include <iostream>
#include <vector>

struct Point {
  float x, y;
  Point(float x, float y):x(x), y(y){}

  Point(const Point& point):x(point.x), y(point.y) {
    std::cout << "copying!" << std::endl;
  }
};

int main() {
  std::vector<Point> points;
  points.reserve(3);
  points.emplace_back(Point(1, 2));
  points.emplace_back(Point(3, 4));
  points.emplace_back(Point(5, 6));
  return 0;
}

Pybind11 安装及简单使用

之前项目里用到了pybind11，效果强大&很好用。但是由于那个项目整个工程直接包含了pybind11的头文件和构建脚本，因此无需自己动手折腾pybind11的环境和构建。最近自己想用pybind11做些POC的小实验，在此记录一些搭建环境的过程。

安装

1. 依赖安装：sudo apt-get install python-dev cmake
2. pybind11 python: conda install pybind11
3. pybind11 安装：
   1. 下载repo：git clone https://github.com/pybind/pybind11.git
   2. cd pybind11
   3. mkdir build
   4. cd build
   5. cmake ..
   6. make check -j8

使用

尝试编译官网的玩具样例，代码如下（保存到文件toy.cc）:

1
2
3
4
5
6
7
8
9
10

#include <pybind11/pybind11.h>

int add(int i, int j) {
  return i + j;
}

PYBIND11_MODULE(example, m) {
  m.doc() = "pybind11 example plugin";
  m.def("add", &add, "A function which adds two numbers");
}

编译命令：

1

c++ -O3 -Wall -shared -std=c++11 -fPIC `python3 -m pybind11 --includes` toy.cc -o shared`python3-config --extension-suffix`

NOTE: 上述命令初看很唬人，我们尝试运行一下python3 -m pybind11 --includes和python3-config --extension-suffix, 得到的结果如下：-I/home/xxx/anaconda3/envs/mindspore/include/python3.7m -I/home/xxx/anaconda3/envs/mindspore/include, .cpython-37m-x86_64-linux-gnu.so。