brainpy 动力学编程基础

文章参考：

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

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

Brainpy手册

文章目录

积分器：
- 定义ODE函数
- 数值积分方法
更新函数和动力系统计算介绍
- 什么是brainpy.DynamicalSystem？
- 如何定义brainpy.DynamicalSystem？
- - 课本内容
  - 实例：用于大脑模拟的 LIF 神经元模型
  - 计算模式
  - - 实例：用于大脑模拟和大脑启发计算的 LIF 神经元模型
  - 模型构成
  - EINet（DynamicalSystem）
- 如何运行brainpy.DynamicalSystem？
- - brainpy.math.for_loop
  - brainpy.DSRunner
突触计算（课本主线）
- 突触连接
- 突触计算
运行器
实例化该网络模型
初始化DSRunner
运行1000ms
可视化结果

前一篇文章介绍了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{cases}\tau\dot{w}=v+a-bw\\\\\dot{v}=v-\frac{v^{3}}{3}-w+I_{\mathrm{ext}}\end{cases}$

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

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

$\begin{aligned}&\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{aligned}$

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

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

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

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

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

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

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

这就是为什么我们需要使用 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 + d t$ ， $S (t + d t) = F (S (t), x, t, d t)$

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.

秒客网