CSP——脉冲神经网络

题目背景

在本题中，你需要实现一个 SNN（spiking neural network，脉冲神经网络）的模拟器。一个 SNN 由以下几部分组成：

1. 神经元：按照一定的公式更新内部状态，接受脉冲并可以发放脉冲
2. 脉冲源：在特定的时间发放脉冲
3. 突触：连接神经元-神经元或者脉冲源-神经元，负责传递脉冲

题目描述

神经元会按照一定的规则更新自己的内部状态。本题中，我们对时间进行离散化处理，即设置一个时间间隔 Δt，仅考虑时间间隔整数倍的时刻 t=kΔt(k∈Z+)，按照下面的公式，从 k−1 时刻的取值计算 k 时刻的变量的取值：

vk=vk−1+Δt(0.04vk−12+5vk−1+140−uk−1)+Ik

uk=uk−1+Δta(bvk−1−uk−1)

其中 v 和 u 是神经元内部的变量，会随着时间而变化，a 和 b 是常量，不会随着时间变化；其中 Ik 表示该神经元在 k 时刻接受到的所有脉冲输入的强度之和，如果没有接受到脉冲，那么 Ik=0。当进行上面的计算后，如果满足 vk≥30，神经元会发放一个脉冲，脉冲经过突触传播到其他神经元；同时，vk 设为 c 并且 uk 设为 uk+d，其中 c 和 d 也是常量。图 1 展示了一个神经元 v 变量随时间变化的曲线。

图1: 神经元 v 变量随时间变化的曲线

突触表示的是神经元-神经元、脉冲源-神经元的连接关系，包含一个入结点和一个出结点（可能出现自环和重边）。当突触的入结点（神经元或者脉冲源）在 k 时刻发放一个脉冲，那么在传播延迟 D(D>0) 个时刻以后，也就是在 k+D 时刻突触的出结点（神经元）会接受到一个强度为 w 的脉冲。

脉冲源在每个时刻以一定的概率发放一个脉冲，为了模拟这个过程，每个脉冲源有一个参数 0<r≤32,767，并统一采用以下的伪随机函数：

C++ 版本：

static unsigned long next = 1;
 
/* RAND_MAX assumed to be 32767 */
int myrand(void) {
    next = next * 1103515245 + 12345;
    return((unsigned)(next/65536) % 32768);
}

Python 版本：

next = 1
def myrand():
    global next
    next = (next * 1103515245 + 12345) % (2 ** 64)
    return (next // 65536) % 32768

Java 版本：

long next = 1;
int myrand() {
    next = next * 1103515245 + 12345;
    return (int)((Long.divideUnsigned(next, 65536)) % 32768);
}

在每个时间刻，按照编号顺序从小到大，每个脉冲源调用一次上述的伪随机函数，当 r>myrand() 时，在当前时间刻发放一次脉冲，并通过突触传播到神经元。

进行仿真的时候，已知 0 时间刻各个神经元的状态，从 1 时间刻开始按照上述规则进行计算，直到完成 T 时刻的计算，再输出 T 时刻神经元的 v 值和发放的脉冲次数分别的最小值和最大值。

规定输入数据中结点按如下方式顺序编号：[0,N−1] 为神经元的编号，[N,N+P−1] 为脉冲源的编号。

代码中请使用双精度浮点类型。

输入格式

从标准输入读入数据。

输入的第一行包括四个以空格分隔的正整数 N S P T，表示一共有 N 个神经元，S 个突触和 P 个脉冲源，输出时间刻 T 时神经元的 v 值。

输入的第二行是一个正实数 Δt，表示时间间隔。

输入接下来的若干行，每行有以空格分隔的一个正整数 RN 和六个实数 v u a b c d，按顺序每一行对应 RN 个具有相同初始状态和常量的神经元：其中 v u 表示神经元在时刻 0 时的变量取值；a b c d 为该神经元微分方程里的四个常量。保证所有的 RN 加起来等于 N。它们从前向后按编号顺序描述神经元，每行对应一段连续编号的神经元的信息。

输入接下来的 P 行，每行是一个正整数 r，按顺序每一行对应一个脉冲源的 r 参数。

输入接下来的 S 行，每行有以空格分隔的两个整数 s(0≤s<N+P)、t(0≤t<N) 、一个实数 w(w≥0) 和一个正整数 D，其中 s 和 t 分别是入结点和出结点的编号；w 和 D 分别表示脉冲强度和传播延迟。

输出格式

输出到标准输出。

输出共有两行，第一行由两个近似保留 3 位小数的实数组成，分别是所有神经元在时刻 T 时变量 v 的取值的最小值和最大值。第二行由两个整数组成，分别是所有神经元在整个模拟过程中发放脉冲次数的最小值和最大值。

只要按照题目要求正确实现就能通过，不会因为计算精度的问题而得到错误答案。

样例1输入

1 1 1 10
0.1
1 -70.0 -14.0 0.02 0.2 -65.0 2.0
30000
1 0 30.0 2

样例1输出

-35.608 -35.608
2 2

样例1解释

该样例有 1 个神经元、1 个突触和 1 个脉冲源，时间间隔 Δt=0.1。唯一的脉冲源通过脉冲强度为 30.0、传播延迟为 2 的突触传播到唯一的神经元。

该样例一共进行 10 个时间步的模拟，随机数生成器生成 10 次随机数如下：

16838
5758
10113
17515
31051
5627
23010
7419
16212
4086

因此唯一的脉冲源在时刻 1-4 和 6-10 发放脉冲。在时间刻从 1 到 10 时，唯一的神经元的 v 取值分别为：

-70.000
-70.000
-40.000
-8.200
-65.000
-35.404
-32.895
0.181
-65.000
-35.608

该神经元在时刻 5 和时刻 9 发放，最终得到的 v=−35.608 。

样例2输入

2 4 2 10
0.1
1 -70.0 -14.0 0.02 0.2 -65.0 2.0
1 -69.0 -13.0 0.04 0.1 -60.0 1.0
30000
20000
2 0 15.0 1
3 1 20.0 1
1 0 10.0 2
0 1 40.0 3

样例2输出

-60.000 -22.092
1 2

子任务

读完题目，理解题目的意思后就知道是道模拟题，按照题目的意思，一步步模拟实现就好了，当然要注意一下代码的时间和空间复杂度。先放上AC代码稳定军心。

AC代码：

#include<iostream>
#include<vector>
using namespace std;
static unsigned long Next = 1;//题目的next要改一下，不然会报错 
/* RAND_MAX assumed to be 32767 */
int myrand(void) 
{
    Next = Next * 1103515245 + 12345;
    return((unsigned)(Next/65536) % 32768);
}
struct edg//突触结构体 
{
	int end;
	double w;
	int D;
};
double v[1001];//神经元的状态参数，v,u,a,b,c,d 
double u[1001];
double a[1001], b[1001], c[1001], d[1001]; 
int r[1001]={0};//脉冲源的r参数 
vector<edg>e[2001];//突触，下标为起始发送源，存储的是该发送源能到达的目标神经元 
double ik[1001][1001]={0};//行为时间点，列为目标神经元 
int sum[1001]={0};//神经元的发送次数 
int main()
{
	int n=0, s=0, p=0, T=0;
	scanf("%d%d%d%d", &n, &s, &p, &T);
	double dt=0.0;
	scanf("%lf", &dt);
	int sum_n=0, rn=0, g=0;
	while(sum_n<n)//输入神经元 
	{
		scanf("%d%lf%lf%lf%lf%lf%lf", &rn, &v[g], &u[g], &a[g], &b[g], &c[g], &d[g]);
		g++;
		for(int i=1; i<rn; i++)
		{
			v[g]=v[g-1];
			u[g]=u[g-1];
			a[g]=a[g-1];
			b[g]=b[g-1];
			c[g]=c[g-1];
			d[g]=d[g-1];
			g++;
		}
		sum_n=sum_n+rn;
	}
	for(int i=0; i<p; i++)//输入脉冲源 
	{
		scanf("%d", &r[i]);
	}
	for(int i=0; i<s; i++)//输入突触 
	{
		edg p;
		int start=0;
		scanf("%d%d%lf%d", &start, &p.end, &p.w, &p.D);
		e[start].push_back(p);
	}
	//开始模拟 
	for(int t=1; t<=T; t++)
	{
		int temp_t=t%1000;//T最大为1000，后面的时刻可以把前面用过的t给覆盖掉 
		for(int i=0; i<n; i++)//神经元遍历 
		{
			double v_pre=v[i];
			v[i]=v[i]+dt*(0.04*v[i]*v[i]+5*v[i]+140-u[i])+ik[temp_t][i];
			ik[temp_t][i]=0;//用完之后清零，防止%1000的时候再次读到前面的脉冲 
			u[i]=u[i]+dt*a[i]*(b[i]*v_pre-u[i]);
			if(v[i]>=30)
			{
				sum[i]++;
				v[i]=c[i];
				u[i]=u[i]+d[i];
				for(int j=0; j<e[i].size(); j++)
				{
					ik[(t+e[i][j].D)%1000][e[i][j].end]+=e[i][j].w;
				}	
			}
		}
		for(int i=0; i<p; i++)//脉冲源遍历 
		{
			if(r[i]>myrand())
			{
				for(int j=0; j<e[i+n].size(); j++)
				{
					//脉冲源的编号是在神经元之后，即（i+n）
					ik[(t+e[i+n][j].D)%1000][e[i+n][j].end]+=e[i+n][j].w;
				}		
			}
		}
	}
	double v_min=v[0], v_max=v[0];
	int sum_min=sum[0], sum_max=sum[0];
	for(int i=1; i<n; i++)
	{
		if(v_min>v[i])v_min=v[i];
		if(v_max<v[i])v_max=v[i];
		if(sum_min>sum[i])sum_min=sum[i];
		if(sum_max<sum[i])sum_max=sum[i];
	}
	printf("%.3lf %.3lf\n%d %d", v_min, v_max, sum_min, sum_max);
	return 0;
}

下面就来讲解一下代码以及几点疑惑：

1、为什么v,u,a,b,c,d不直接弄个神经元结构体，而是分别用数组存储？在空间上，结构体存储和各个参数分开存储没有区别，但是在获取数据的方式上略微不同，将导致耗费的时间不同。例如，要获取神经元结构体中的参数v，计算机要先在神经元结构体数组中找到该神经元，然后再在该神经元结构体这块内存中再次寻找参数v（类比间接寻址），而用数组存储不同，计算机直接在数组v中顺序查找，直接获取对应参数v（类比直接寻址）。因为本题卡时长，所以每一步优化都很关键。

2、设立vector<edg> e[2001]，类似于邻接表，数组e的下标是发送脉冲的脉冲源或者神经元，每个数组元素都是一个vector，表示该脉冲源或者神经元将影响到哪些神经元以及脉冲w为多大。因为脉冲源和神经元都是是1000，所以数组e要设2001。

3、ik二维数组为什么是行为时间点，列为目标神经元，而不是行为目标神经元，列为时间点？在模拟的时候，最外围for循环是时间，然后依次遍历所有的神经元在该时刻所受的脉冲之和，即ik。我们知道二维数组的写入和读取都是一行一行的执行，先执行完一行再执行下一行，当ik中行为目标神经元，列为时间点，在每一时刻，每次访问神经元对应ik都要一行一行的执行，就会造成时间上的浪费。而行为时间点，列为目标神经元时，当访问该时刻时，一行都是接下来要访问的神经元，所以缩短了访问时间。

4、为什么时间要取模1000？因为直接数组e二维数组开10^5*10^3，

空间复杂度度：（4+8+4）*10^5*10^3=1600MB，超内存了，所以不能开10^5。取模1000，因为D最大是1000，当前面时刻的ik用完之后，可以清零直接覆盖掉，不会产生影响。

其它需要注意的点，注释上也有说明，如果还有不清楚的可以留言哦！