文章参考：

《神经计算建模实战——基于brainpy》 吴思

【brainpy学习笔记】基础知识2(动力学模型的编程基础)-CSDN博客

Brainpy手册

前一篇文章介绍了JIT编译环境基本要素：数据操作和控制流，结合面向对象的思想，可以进行Brainpy的编程。

如何运行动力学模型，怎么搭建一个框架，并逐渐填充各个模块的内容。

了解微分方程及其数值和数值积分器的使用方式，如何定义一个动力学模型，使用BrainPy提供的各类运行器进行模拟训练和分析。

积分器：

介绍积分器的使用方法。求解微分方程是神经动力学模拟的核心。支持常微分方程（ODE）随机微分方程（SDE），分数阶微分方程（FDE）和延迟微分方程（DDE）。

定义ODE函数和对应的 数值计分方法。

定义ODE函数

（1）定义ODE函数

ODE描述系统状态随时间的演变，例如经典的单变量方程了解！让我们从常微分方程（ODE）开始。使用BrainPy实现ODE，首先需要了解如何用数值方法求解ODE以及BrainPy的基本用法。

在BrainPy中，常见的步骤如下：

def diff(x,y,t,p1,p2):

  dx = f1(x,t,y,p1)

  dy = g1(x,t,y,p2)

  return dx,dy

可以将参数分为3个部分：t表示当前时刻，在时刻t前传递的x和y表示动态变量，在时刻t后传递的p1和p2表示系统需要的参数（不随时间变化）。

在函数主体中，f1和g1根据用户需求定制，dx和dy按照与函数参数中对应变量的顺序返回。用户实际定义微分方程时，需要注意动态变量必须写在参数t之前。静态变量写在参数之后。

数值积分方法

（2）数值积分方法。将积分器Brainpy.odeint作为一个Python装饰器放在微分函数diff上，就可以完成对该函数的数值积分方法的定义。

python@bp.odeint

def diff(x,y,t,p1,p2):

  dx = f1(x,t,y,p1)

  dy = g1(x,t,y,p2)

  return dx,dy

包装后的diff的数据类型是ODEIntegartor的实例。brainpy.odeint的参数如下

method:字符串类型，用于指定对ODE函数数值积分的方法,如 ‘rk4’、‘adams’ 等,默认为 欧拉函数。

dt：浮点数类型，用于设置数值积分步长，默认为0.1.

如果扩展的话有很多函数，但是装饰器直接封装和设定了这些函数值。

• ode_func: 表示微分方程组的函数,需要用户自定义。
• y0: 微分方程组的初始条件。
• t: 待求解的时间点。
• args: 传递给 ode_func 的其他参数,可选。
• method: 数值积分的方法,如 ‘rk4’、‘adams’ 等,默认为 ‘rk4’。
• dt: 数值积分的步长,默认会自动选择。
• events: 检测微分方程解中事件的函数,可选。
• dense_output: 是否返回连续的解,默认为 False。
• atol、rtol: 绝对和相对误差tolerance,用于控制数值积分精度。

@bp.odeint(method='euler',dt=0.01)

def diff(x,y,t,p1,p2):

  dx = f1(x,t,y,p1)

  dy = g1(x,t,y,p2)

  return dx,dy

$$\{\begin{array}{l}\tau \dot{w}=v+a-bw\\ \\ \dot{v}=v-\frac{v^3}{3}-w+I_{\mathrm{ext}}\end{array}$$

$v$,$w$为变量，其余为参数。

下面这个是把参数赋值的公式。

$$\begin{array}{rl}& {\mathbf{v}}^{\prime }=\mathbf{v}-{\mathbf{v}}^{\wedge }3/3-\mathbf{w}+1\\ & {\mathbf{w}}^{\prime }=0.08(\mathbf{v}+0.7-0.8\mathbf{w})\end{array}$$

import brainpy as bp

def FitzHugh_Nagumo(state, t, I_ext, a=0.7, b=0.8, tau=0.08):
    """
    FitzHugh-Nagumo neuron model.

    Parameters:
    state (list): the state variables [v, w]
    t (float): the time point
    I_ext (float): the external current input
    a (float): constant parameter
    b (float): constant parameter 
    tau (float): time constant

    Returns:
    list: the derivatives of the state variables [dv/dt, dw/dt]
    """
    v, w = state
    dv = v - v**3/3 - w + I_ext
    dw = (v + a - b*w) / tau
    return dv, dw

@bp.odeint(method='euler',dt=0.01)

#V,w是系统的状态变量。t是时间变量，Iext是外部输入，a,b,tau是系统参数

def integral(V,w,t,Iext,a,b,tau):

  dv = v - v**3/3 - w + Iext #定义状态变量v的导数

  dw = (v + a - b*w) / tau  #状态变量w的导数

  return dv, dw #返回导数

如果直接执行的话，和t无关。执行结果。

原因的说明：传入的t没有被使用。只是初始状态的计算。

def integral(V, w, t, Iext, a, b, tau):
    dv = v - v**3/3 - w + Iext 
    dw = (v + a - b*w) / tau
    return dv, dw

1. 它接收当前状态 (V, w) 作为输入
2. 计算在这个状态下的导数 (dv, dw)
3. 返回这些导数值

所以无论你传入什么 t 值，结果都只取决于：

• 状态变量 V 和 w 的当前值
• 参数 Iext, a, b, tau 的值

真实的方程
dV/dt = V - V³/3 - w + Iext
dw/dt = (V + a - b*w) / tau

告诉我们在某个特定状态点 (V,w) 下系统将如何变化，但要获得随时间演化的完整解，我们需要：

1. 从初始点开始
2. 使用数值积分方法（如龙格库塔法）
3. 用小的时间步长逐步积分
4. 记录每个时间点的状态

这就是为什么我们需要使用 IntegratorRunner 或其他积分器来获得完整的时间演化解。

• 它包装了积分过程
• 自动处理时间步进
• 存储中间结果
• 提供了方便的监视和绘图功能

在一个动力可能有多个变量随时间动态变化。有时这些变量是互相联系的，更新一个变量时需要输入其他变量。为了达到更高的积分精度。

可以使用brainpy.JointEq.

a,b=0.02,0.20

dV = lambda V,t,w,Iext:0.04*V*V+5*V+140-w+Iext #四个参数，一个输出的函数

dw = lambda w,t,V:a*(b*V-w)  #三个参数，

joint_eq = bp.JointEq([dV,dw])

integra12 = bp.odeint(joint_eq,method ='rk2')

使用积分求解器来计算出一个微分方程的数值解

a=0.7; b=0.8; tau=12.5; Iext=1.
runner = bp.integrators.IntegratorRunner(
    integral,         #包装好的方程组
    monitors=['V'], #监视记录方程组中的动态变量
    inits=dict(V=0.,w=0.), #以字典的方式初始化动态变量
    args=dict(a=a,b=b,tau=tau,Iext=Iext), #以字典的方式传递所需要的参数
    dt=0.01 #步长为0.01
)
使用积分运行器
runner.run(100.)
 
plt.plot(runner.mon.ts , runner.mon.V) #作图
plt.show()

除了调用，此处给出ai求解积分的代码实现。

import numpy as np

class SimpleIntegrator:
    def __init__(self, func, dt=0.01):
        """
        初始化积分器
        func: 微分方程函数,返回导数
        dt: 时间步长
        """
        self.func = func
        self.dt = dt
        
    def integrate(self, t_span, initial_conditions, args=()):
        """
        执行数值积分
        t_span: (start_time, end_time)
        initial_conditions: 初始状态值的列表或数组
        args: 传递给微分方程的额外参数
        """
        t_start, t_end = t_span
        n_steps = int((t_end - t_start) / self.dt)
        
        # 创建时间数组
        t = np.linspace(t_start, t_end, n_steps)
        
        # 初始化结果数组
        n_vars = len(initial_conditions)
        results = np.zeros((n_steps, n_vars))
        results[0] = initial_conditions
        
        # 使用四阶Runge-Kutta方法进行积分
        for i in range(1, n_steps):
            y = results[i-1]
            
            # RK4步骤
            k1 = np.array(self.func(*y, t[i-1], *args))
            k2 = np.array(self.func(*(y + k1 * self.dt/2), t[i-1] + self.dt/2, *args))
            k3 = np.array(self.func(*(y + k2 * self.dt/2), t[i-1] + self.dt/2, *args))
            k4 = np.array(self.func(*(y + k3 * self.dt), t[i], *args))
            
            # 更新状态
            results[i] = y + (k1 + 2*k2 + 2*k3 + k4) * self.dt/6
            
        return t, results

# 使用示例
def use_simple_integrator(V0, w0, t_span, Iext, a, b, tau):
    integrator = SimpleIntegrator(integral, dt=0.01)
    t, results = integrator.integrate(
        t_span=(0, t_span),
        initial_conditions=[V0, w0],
        args=(Iext, a, b, tau)
    )
    return t, results
# 使用简化版积分器
t, results = use_simple_integrator(
    V0=-1.0, 
    w0=0.0, 
    t_span=100, 
    Iext=1., 
    a=0.7, 
    b=0.8, 
    tau=12.5
)

# 绘图
plt.plot(t, results[:, 0])  # 画V的变化
plt.plot(t, results[:, 1]) 
plt.show()

同样有运行的结果。

在文档中搜索Numerical integration。

我们可以发现，关于数值积分求解的都在brainpy.integrators的库中。

toolbox处有相关的对应的微分方程的计算和求解

让我们看一个简单的例子，比如一个指数增长模型 $\frac{dy}{dt}=ay$，其中 (a) 是常数。你可以用如下代码实现它：

这些文档给出了计算式，还有简单的运行实例。修改方程或参数，尝试其他ODE模型。

搜索中，还有对应的求解算法的使用过程。

更新函数和动力系统计算介绍

所有这些支持都基于一个共同的概念： Dynamical System通过brainpy.DynamicalSystem进行脑的动力学仿真和方程的求解

更新函数部分的章节来自于动力学系统介绍：

该章节主要分为三个部分

什么是brainpy.DynamicalSystem？

如何定义brainpy.DynamicalSystem？(更新函数是在这个部分的。)

如何运行brainpy.DynamicalSystem？

本章组织，把课本内容和文档进行交叉对照

上一节，通过brainpy.odeint等积分器定义了模型的连续积分过程，并使用IntegratorRunner获得了一段时间内的数值积分结果。然而，一个动力学模型往往不仅包含连续的积分部分，还包含不连续的更新步骤。因此BrianPy提供了一个通用的DynamicalSystem类，以定义各种动力学模型。

什么是brainpy.DynamicalSystem？

所有用于大脑模拟和大脑启发计算的模型都是动力学系统（DynamicalSystem)

DynamicalSystem 是 BrainPyOject 的子类。 因此，它支持使用前面教程中所述的面向对象转换。

动态系统（DynamicalSystem）定义了模型在单个时间步的更新规则。

1.对于有状态的模型，DynamicalSystem 定义了从$t$到$t+dt$，$S(t+dt)=F(S(t),x,t,dt)$

• S是状态
• x是输入
• t是时间
• dt是时间步

大脑模拟中广泛使用的递归神经网络（如 GRU、LSTM）、神经元模型（如 HH、LIF）或突触模型就是这种情况。

然而，对于深度学习中的模型，如卷积和全连接线性层，DynamicalSystem 定义了输入到输出映射$y=F(x,t)$。

此处文档和书有一定的差异了。

所有的DynamicalSystem类,都需要定义一个更新函数update() .该函数规定了模型的状态从当前时刻t更新到下一个时刻t+dt的过程。适用于多种场景。

如何定义brainpy.DynamicalSystem？

首先，所有 DynamicalSystem 都应实现 .update() 函数，该函数接收两个参数：

class YourModel(bp.DynamicalSystem):
  def update(self, s, x):
    pass

update接受两类参数：公共参数和私有参数。

公共参数指一个网络中所有节点或所有该类实例都可以获取的参数。

私有参数指某个节点或DynamicalSystem实例独有的参数。

• s (or named as others): A dict, to indicate shared arguments across all nodes/layers in the network, like

  • the current time t, or当前时刻t
  • the current running index i, or当前迭代次数i
  • the current time step dt, or单步积分时长
  • the current phase of training or testing fit=True/False.当前是训练还是测试阶段

• x (or named as others): The individual input for this node/layer.（私有参数）该节点/层的单独输入

我们之所以称 s 为共享参数，是因为它们对所有节点/层来说都是相同和共享的。 相反，不同的节点/层有不同的输入 x。

update函数可以不接受私有参数，但公共参数是每个DynamicalSystem类都需要有的。公共参数往往以字典的方式存储。

  @bp.odeint(method='euler',dt=0.01)

#V,w是系统的状态变量。t是时间变量，Iext是外部输入，a,b,tau是系统参数

def integral(V,w,t,Iext,a,b,tau):

  dv = v - v**3/3 - w + Iext #定义状态变量v的导数

  dw = (v + a - b*w) / tau  #状态变量w的导数

  return dv, dw #返回导数

课本内容

class FitzHughNagumo(bp.DynamicalSystem):
    def __init__(self,size,a=0.7,b=0.8,tau=12.5):
        super(FitzHughNagumo,self).__init__()
        #参数
        self.a=a
        self.b=b
        self.tau=tau
        #变量
        self.V=bm.Variable(bm.ones(size))
        self.W=bm.Variable(bm.zeros(size))
        self.input = bm.Variable(bm.zeros(size))
        #积分函数
        self.integral = bp.odeint(bp.JointEq(self.dV,self.dw))
    def dV(self,V,t,w,I):#V的微分方程
        return V - V * V * V/3 -w + I #I是外部输入
    def dw(self,w,t,V): #w的微分方程
        return (V+a -b*w)/tau
    def update(self,tdi):#更新函数
        #integral函数在这一步完成了从t到t+dt的积分。
		self.V.value,self.w.value=self.integral(self.V,self.w,tdi.t,self.input,tdi.dt)
        #积分完成后input被置零
        self.input[:]=0.
        

runner = bp.integrators.IntegratorRunner(
    integral,         #包装好的方程组
    monitors=['V'], #监视记录方程组中的动态变量
    inits=dict(V=0.,w=0.), #以字典的方式初始化动态变量
    args=dict(a=a,b=b,tau=tau,Iext=Iext), #以字典的方式传递所需要的参数
    dt=0.01 #步长为0.01
)

另一个例子：

实例： 用于大脑模拟的 LIF 神经元模型

在此，我们使用 “Leaky Integrate-and-Fire”（LIF）模型来说明动态系统的第一个约束条件。

LIF 模型是在大脑模拟中首次提出的神经元动力学建模方法。 其方程为
$$\begin{gathered} {\tau }_m\frac{dV}{dt}=-(V(t)-V_{rest})+I(t) \\ \mathrm{if}V(t)>V_{th},V(t)=V_{rest} \end{gathered}$$
${\tau }_m$：膜时间常数,决定了膜电位对输入信号的响应速度。它表示膜电位从某个值衰减到剩下63.2%所需的时间。(确定的参数)

$\frac{dV}{dt}$: 膜电位随时间的变化率。这个是方程左边的导数项。

$V(t)-V_{rest}$: 膜电位相对于静息电位的偏移量。

$I(t)$: 注入到神经元的电流输入。

如果膜电位 $V(t)$ 大于阈值 $V_{th}$,则膜电位重置为静息电位 $V_{rest}$。这模拟了动作电位的产生和重置。

综合起来,这个方程描述了膜电位如何随时间变化:

1. 膜电位会向静息电位 $V_{rest}$ 衰减,衰减速度由 ${\tau }_m$ 决定。

2. 注入的电流输入 $I(t)$ 会使膜电位上升。

3. 一旦膜电位超过阈值 $V_{th}$,就会重置为静息电位 $V_{rest}$,模拟动作电位的产生。

class LIF_for_BrainSimulation(bp.DynamicalSystem):
  def __init__(self, size, V_rest=0., V_th=1., tau=5., mode=None):
    super().__init__(mode=mode)

    # this model only supports non-batching mode
    # 此模型仅支持非批处理模式，数据会被实时地、逐个地处理。
    bp.check.is_subclass(self.mode, bm.NonBatchingMode)

    # parameters
    #固定参数
    self.size = size
    self.V_rest = V_rest
    self.V_th = V_th
    self.tau = tau

    # variables
    # 变量
    self.V = bm.Variable(bm.ones(size) * V_rest)
    self.spike = bm.Variable(bm.zeros(size, dtype=bool))

    # integrate differential equation with exponential euler method
    #用指数欧拉法整合微分方程
    self.integral = bp.odeint(f=lambda V, t, I: (-V + V_rest + I)/tau, method='exp_auto')

  def update(self, s, x):
    # define how the model states update
    # according to the external input
    # s是一个包含当前时间t和时间步长dt的字典
    # x 是给定的外部输入电流
    t, dt = s.get('t'), s.get('dt')
    # 变量的计算
    V = self.integral(self.V, t, x, dt=dt)
    #是否激活
    spike = V >= self.V_th
    #判断语句，根据这个，返回静息态，还是现有的电压
    self.V.value = bm.where(spike, self.V_rest, V)
    self.spike.value = spike
    return spike

计算模式

其次，明确考虑动态系统支持哪种计算模式。

大脑仿真通常在没有批处理维度的情况下构建模型（我们称之为非批处理模式，如上述 LIF 模型所示），而脑启发计算则通过批处理数据（批处理模式或训练模式）来训练模型。

因此，要编写一个适用于国外大脑仿真和脑启发计算应用的模型，你需要考虑你的模型支持哪种模式，是其中之一，还是两者都支持。

实例： 用于大脑模拟和大脑启发计算的 LIF 神经元模型

在考虑计算模式时，我们可以为大脑模拟和大脑启发计算编程一个通用的 LIF 模型。

为了克服大脑模拟 LIF 模型中尖峰的非差分特性，

spike = V >= self.V_th

脑启发计算中使用的 LIF 模型使用代梯度函数计算尖峰状态。 通常，我们用一个平滑函数代替尖峰的后向梯度，如
$$g^{\prime }(x)=\frac{1}{(\alpha \ast |x|+1.{)}^2}$$

脉冲神经网络阶跃函数(激活，或者不激活的0，1阶跃函数)，存在不可微的问题，所以，指定一些替代函数，可以解决该问题。

def inv_square_grad(x, alpha=2.):
    return 1. / (1. + alpha * x * x)

class LIF(bp.DynamicalSystem):
  def __init__(self, size, f_surrogate=None, V_rest=0., V_th=1., tau=5.,mode=None):
    super().__init__(mode=mode)
    #设置有非批次模式，批次模式，训练模式
    bp.check.is_subclass(self.mode, [bm.NonBatchingMode, bm.BatchingMode, bm.TrainingMode])

    # Parameters
    self.size = size
    self.num = bp.tools.size2num(size)
    self.V_rest = V_rest
    self.V_th = V_th
    self.tau = tau#tau 是膜电位的时间常数，表示膜电位衰减到初始值的 63.2% 所需的时间。
    # 一个替代函数，用于替代非微分函数Heaviside。
    # 替代梯度函数，用于解决脉冲神经网络中阶跃函数不可微的问题。
    if f_surrogate is None:
      f_surrogate = bm.surrogate.inv_square_grad
    self.f_surrogate = f_surrogate
	# 添加 varshape 属性
    c
    # integrate differential equation with exponential euler method
    # 表示膜电位衰减到初始值的 63.2% (1-1/e) 所需的时间
    # f是返回的微分
    self.integral = bp.odeint(f=lambda V, t, I: (-V + V_rest + I)/tau, method='exp_auto')

    # Initialize a Variable:
    # - if non-batching mode, batch axis of V is None
    # - if batching mode,     batch axis of V is 0
    # 按照不同的mode，来初始化V
    self.V = bp.init.variable_(bm.zeros, self.size, self.mode)
    self.V[:] = self.V_rest
    self.spike = bp.init.variable_(bm.zeros, self.size, self.mode)
  #重置神经元的状态，根据批次大小重置膜电位和动作电位
  def reset_state(self, batch_size=None):
    self.V.value = bp.init.variable_(bm.ones, self.size, batch_size) * self.V_rest
    self.spike.value = bp.init.variable_(bm.zeros, self.size, batch_size)
  #
  def update(self, s, x):
    t, dt = s.get('t'), s.get('dt', bm.dt)
    V = self.integral(self.V, t, x, dt=dt)
    # replace non-differential heaviside function
    # with a surrogate gradient function
    spike = self.f_surrogate(V - self.V_th)
    # reset membrane potential
    # 重置V和激活的变量
    self.V.value = (1. - spike) * V + spike * self.V_rest
    self.spike.value = spike
    return spike

tau是膜时间常数，描述了神经元膜电位变化的时间尺度。

模型构成

我们上面定义的 LIF 模型可以递归地组成大脑仿真和大脑启发计算中的网络。

下面的代码片段利用 LIF 模型构建了 E/I 平衡网络 EINet，这是大脑仿真中的一个经典网络模型。

运行下面的代码，需要在LIF中添加下面这个属性，以满足突触函数设置的条件。

 # 添加 varshape 属性      
 self.varshape = self.size

EINet（DynamicalSystem）

class EINet(bp.DynamicalSystem):
  def __init__(self, num_exc, num_inh):
    super().__init__()
    self.E = LIF(num_exc, V_rest=-55, V_th=-50., tau=20.)
    self.I = LIF(num_inh, V_rest=-55, V_th=-50., tau=20.)
    self.E2E = bp.synapses.Exponential(self.E, self.E, bp.conn.FixedProb(0.02),
                                       g_max=1.62, tau=5., output=None)
    self.E2I = bp.synapses.Exponential(self.E, self.I, bp.conn.FixedProb(0.02),
                                       g_max=1.62, tau=5., output=None)
    self.I2E = bp.synapses.Exponential(self.I, self.E, bp.conn.FixedProb(0.02),
                                       g_max=-9.0, tau=10., output=None)
    self.I2I = bp.synapses.Exponential(self.I, self.I, bp.conn.FixedProb(0.02),
                                       g_max=-9.0, tau=10., output=None)

  def update(self, s, x):
    # x is the background input
    e2e = self.E2E(s)
    e2i = self.E2I(s)
    i2e = self.I2E(s)
    i2i = self.I2I(s)
    self.E(s, e2e + i2e + x)
    self.I(s, e2i + i2i + x)

with bm.environment(mode=bm.nonbatching_mode):
  net1 = EINet(3200, 800)

此外，我们的 LIF 模型还可用于脑启发计算场景。 下面的 AINet 使用 LIF 模型构建了一个人工智能训练模型。

# This network can be used in AI applications

class AINet(bp.DynamicalSystem):
  def __init__(self, sizes):
    super().__init__()
    self.neu1 = LIF(sizes[0])
    self.syn1 = bp.layers.Dense(sizes[0], sizes[1])
    self.neu2 = LIF(sizes[1])
    self.syn2 = bp.layers.Dense(sizes[1], sizes[2])
    self.neu3 = LIF(sizes[2])

  def update(self, s, x):
    x = self.neu1(s, x)
    x = self.syn1(s, x)
    x = self.neu2(s, x)
    x = self.syn2(s, x)
    x = self.neu3(s, x)
    return x

with bm.environment(mode=bm.training_mode):
  net2 = AINet([100, 50, 10])

如何运行brainpy.DynamicalSystem？

在构建上述的底层神经元，神经元网络，之后，该如何运行上面的结构呢

和第一部分积分一样，如何运行求解。

如上所述，DynamicalSystem 只定义了单个时间步的更新规则，因此要长期运行一个 DynamicalSystem 实例，我们需要一个 for 循环机制。

brainpy.math.for_loop

for_loop 是一个结构性控制流应用程序接口，用于运行一个对输入进行循环的函数。

此外，该应用程序接口还能及时将循环过程编译为机器代码。 假设我们有 200 个时间步，步长为 0.1，我们可以用以下方法运行模型：

with bm.environment(dt=0.1):
  # construct a set of shared argument with the given time steps
  shared = bm.shared_args_over_time(num_step=200)
  # construct the inputs with shape of (time, batch, feature)
  currents = bm.random.rand(200, 10, 100)

  # run the model
  net2.reset_state(batch_size=10)
  out = bm.for_loop(net2, (shared, currents))

out.shape

brainpy.DSRunner

在 BrainPy 中运行模型的另一种方法是使用结构运行对象 DSRunner 和 DSTrainer。 它们为监控动态系统中的变量提供了更灵活的方式。 详细内容请参阅 DSRunner 教程。

with bm.environment(dt=0.1):
  runner = bp.DSRunner(net1, monitors={'E.spike': net1.E.spike, 'I.sike': net1.I.spike})
  runner.run(inputs=bm.ones(1000) * 20.)

bp.visualize.raster_plot(runner.mon['ts'], runner.mon['E.spike'])

突触计算（课本主线）

在实际网络模型中，突触往往占据大部分内存和计算时间。因此，研究如何高效地存储和计算突触是加快网络运行的必要步骤。

根据突触的特点来提高计算效率。

突触连接

稀疏连接下。

突触计算通常涉及突触前神经元(“pre”)、突触后神经元(“post”)和突触(“syn”)三种数据之间的转换。

一个突触对应一个从突触前神经元到突触后神经元的唯一映射。

如何存储这样的突触连接方式呢？brainpy.connect.Connector提供了一种统一的格式来生成任意的突触连接存储格式。先实例化一个Connector,再调用Connector.require(),以生成所需要的连接结构。

例如一个概率连接的Connector类定义如下：

conn = bp.connect.FixedProb(prob=0.3, include_self=False ,seed=134)
#prob:连接概率 #include_self:是否有自环 seed:随机数种子
print(conn)
>>>FixedProb(prob=0.3, pre_ratio=1.0, include_self=False, allow_multi_conn=False, seed=134)

在具有4个突触前神经元和四个突触后神经元的情况下，我们使用该连接器自动生成连接矩阵（conn_mat）、突触前连接的神经元索引(pre_ids)、突触后连接的神经元索引(post_ids)、突触前到突触后神经元连接的索引(pre2post)等

conn(pre_size=4,post_size=4)
#FixedProb(prob=0.3, pre_ratio=1.0, include_self=False, allow_multi_conn=False, seed=134)
res = conn.require('conn_mat','pre_ids','post_ids','pre2post','pre2syn') #事先请求所有需要用到的数据结构
res
(Array([[False, False, False,  True],
        [False, False, False, False],
        [ True, False, False, False],
        [ True, False, False, False]], dtype=bool),
 Array([0, 2, 3], dtype=int32),
 Array([3, 0, 0], dtype=int32),
 (Array([3, 0, 0], dtype=int32), Array([0, 1, 1, 2, 3], dtype=int32)))

用矩阵存储连接信息在稀疏连接的情况下显得非常冗余，因为在16个元素中，只有一个是True

。对于稀疏连接，连接矩阵显然不是最优突触连接存储方式。为了更好地降低内存占用率，我们可以用pre_ids和post_ids来存储索引。存储空间占用少，索引具有相同的元素数。

每对pre_ids[i]和post_ids[i]相连形成突触。具体来说，上述例子生成的：

 Array([0, 2, 3], dtype=int32),
 Array([3, 0, 0], dtype=int32),

pre_ids和post_ids是最常见的稀疏连接结构，但在事件驱动下其计算效率不是最高的。

pre2post是一种高效的稀疏连接结构。可以表示一个突触前神经元与哪些突触后神经元相连。

pre2post包含两个数组，第一个数组为：突触后神经元的序号；第二个数组为：与第i号（索引项）突触前神经元相连的那些突触后神经元在第一个数组中的开始和结束位置。所以，第二个数组会多一个。

对于上面的例子，进行了（0，3），（2，0）（3，0）这三对连接。

[3, 0, 0]突触后神经元的排序。因为没有2这个后神经元，因此，不在此列。

[0, 1, 1, 2, 3]

从第一个前神经元开始。取两项（0，1）只有3，即0和3，为一个组合（0，3）。

|3| 0| 0|

0 1 2 3

对于第二个 （1, 1）。无元素.即前神经元1无后突触。

依次类推。

pre2Post包含的数组有：突触后神经元数组（Post-Synaptic Neuron）和突触前神经元指针的指针数组（Pre-Index Vector）.

即图中的下面两行。

前两行是index的第一种保存方式。可以看到通过数组索引，我们把数组的大小进行了压缩。指示前神经元分布到多少。

基于pre2Post，事件驱动更新会非常高效，因为只需要沿着突触前神经元迭代，一旦遇到突触前神经元产生的脉冲发放，就可以取出与该突触前神经元相连的突触后神经元。，并对齐进行更新保证事件驱动更新的特性。

给一个参考文章作者做的图

除了pre2post,BrainPy还提供了一种适用于不同场景的数据结构:pre2syn.

定义与pre2post类似，也是由两个数组构成的元组，只是将突触后神经元数组换成了突触数组。

res[4]

(Array([0, 1, 2, 3], dtype=int32), Array([0, 1, 2, 3, 4], dtype=int32))

突触前神经元和突触连接的索引信息。还可以生成post2pre，post2syn,syn2post等数据结构。

post2pre: 一个 1D 数组,其长度等于突触后神经元的总数。post2pre[j]表示第 j 个突触后神经元是由哪个突触前神经元连接而来。

连接模式是在require()函数中创建的，而不是在初始化connector是创建的。对于概率恒为p的连接模式，当我们初始化Connector会保存构建突触连接的所有参数信息，不会直接构建对应的连接模式，防止空间损耗。

require（）函数只创建用户需要用到的数据，能够大大减小内存占用，这意味着用户调用一次require（），根据概率p生成随机的连接模式，并返回数据。由于随机性，如果需要使用什么，要直接申请所有的数据，从而保证相同的连接模式。

突触计算

基于突触连接，可以高效地进行突触计算。假设已知所有突触前神经原地脉冲事件events,突触连接的权重weight(标量，即所有突触连接的权重相同)，为了获得这些脉冲事件作用在突触后神经元上的电流，可以是共下面的计算方法。

没有电子板的，不想手打了，直接粘贴。

根据存储结构的不同，选择对应的突触计算函数。（事件，存储，权重）。

也不是很懂，往后过。

运行器

上述组件构建，脑动力学模型。这些模型的模拟，训练，分析需要依赖Brainpy提供的运行器。BrainPy提供的

brainpy.DSRunner可以帮助用户运行和模拟动力学模型。

brainpy.train.DSTrainer可以训练动力学模型

brainpy.analysis.DSAnalyzer可以进行动力学分析。

Brainpy运行器的基本声明方式：

runner = SomeRunner(
    target,
    inputs,
    monitors,
    dyn_vars,
    jit
)

在定义好运行器以后，用户可以调用.run()方法来运行，它的参数有：

前面也有EI模型，但是基类不太一样

class EINet(bp.dyn.Network):
    def __init__(self, scale=1.0, method='exp_auto'):
        super(EINet, self).__init__()
 
        num_exc = int(3200 * scale)
        num_inh = int(800 * scale)
 
        pars = dict(V_rest=-60., V_th=-50., V_reset=-60., tau=20., tau_ref=5.)
        self.E = bp.neurons.LIF(num_exc, **pars, method=method)
        self.I = bp.neurons.LIF(num_inh, **pars, method=method)
        self.E.V[:] = bm.random.randn(num_exc) * 2 - 55.
        self.I.V[:] = bm.random.randn(num_inh) * 2 - 55.
 
        prob = 0.1
        we = 0.6 / scale / (prob / 0.02)**2
        wi = 6.7 / scale / (prob / 0.02)**2
        self.E2E = bp.dyn.ExpCOBA(self.E, self.E, bp.conn.FixedProb(prob),
                                  E=0., g_max=we, tau=5., method=method)
        self.E2I = bp.dyn.ExpCOBA(self.E, self.I, bp.conn.FixedProb(prob),
                                  E=0., g_max=we, tau=5., method=method)
        self.I2E = bp.dyn.ExpCOBA(self.I, self.E, bp.conn.FixedProb(prob),
                                  E=-80., g_max=wi, tau=10., method=method)
        self.I2I = bp.dyn.ExpCOBA(self.I, self.I, bp.conn.FixedProb(prob),
                                  E=-80., g_max=wi, tau=10., method=method)
 
# 实例化该网络模型
net = EINet(scale=1., method='exp_auto')
# 初始化DSRunner
runner = bp.dyn.DSRunner(
    net, # 目标模型
    monitors={'E.spike': net.E.spike}, # 监视动态变量spike
    inputs=[(net.E.input, 20.), (net.I.input, 20.)] # 两群神经元的输入大小为恒定值20.
)
# 运行1000ms
runner.run(1000.)
# 可视化结果
bp.visualize.raster_plot(runner.mon.ts, runner.mon['E.spike'], show=True)

训练的调用方法和Pytorch和sklearn类似。

trainer = bp.DSTrainer(model)
trainer.fit([X,Y])
trainer.predict(X)

ob(prob),
E=0., g_max=we, tau=5., method=method)
self.I2E = bp.dyn.ExpCOBA(self.I, self.E, bp.conn.FixedProb(prob),
E=-80., g_max=wi, tau=10., method=method)
self.I2I = bp.dyn.ExpCOBA(self.I, self.I, bp.conn.FixedProb(prob),
E=-80., g_max=wi, tau=10., method=method)

实例化该网络模型

net = EINet(scale=1., method=‘exp_auto’)

初始化DSRunner

runner = bp.dyn.DSRunner(
net, # 目标模型
monitors={‘E.spike’: net.E.spike}, # 监视动态变量spike
inputs=[(net.E.input, 20.), (net.I.input, 20.)] # 两群神经元的输入大小为恒定值20.
)

运行1000ms

runner.run(1000.)

可视化结果

bp.visualize.raster_plot(runner.mon.ts, runner.mon[‘E.spike’], show=True)

[外链图片转存中...(img-M31nulIp-1730988902104)]

训练的调用方法和Pytorch和sklearn类似。

```python
trainer = bp.DSTrainer(model)
trainer.fit([X,Y])
trainer.predict(X)