文献链接：https://doi.org/10.1038/s42256-021-00430-y
源码文件下载链接：https://github.com/ykubo82/bioCHL

基本设置

导入的工具包

import numpy as np
from keras.datasets import mnist

import torch

import os
from sklearn.preprocessing import OneHotEncoder
from sklearn import utils  
from warnings import filterwarnings
filterwarnings('ignore')
import time

1. ‘os’:这是Python标准库中的一个模块，提供了与操作系统交互的功能。通过os模块，您可以执行各种文件和目录操作，例如创建、删除、移动、重命名文件，获取当前工作目录，执行系统命令等。它是Python中常用的一个模块，尤其在处理文件和目录时非常有用。
2. ‘sklearn.preprocessing.OneHotEncoder’:这是scikit-learn（简称sklearn）库中的一个类，用于执行独热编码（One-Hot Encoding）操作。独热编码是将分类变量转换为二进制向量的过程，使得每个类别都在向量中表示为一个单独的二进制特征。它通常用于处理机器学习算法中需要处理分类数据的情况，例如逻辑回归、支持向量机等。
3. ‘sklearn.utils’:这是scikit-learn库中的一个模块，提供了一些实用工具函数。其中的函数可以用于数据集的操作、交叉验证、类别平衡处理等。它帮助简化了在机器学习任务中的一些常见操作。
4. ‘warnings.filterwarnings’:这是Python标准库中warnings模块的filterwarnings函数。通过调用filterwarnings函数并传递’ignore’参数，我们可以在运行代码时忽略警告信息。警告通常用于提醒开发者潜在的问题或不推荐的用法，但有时可能会干扰代码的输出。在特定情况下，可以使用**filterwarnings(‘ignore’)**来暂时忽略这些警告。
5. ‘time’:这是Python标准库中的一个模块，提供与时间相关的功能。通过time模块，您可以测量程序的执行时间、进行时间延迟等。它在需要对代码执行时间进行评估或进行时间相关操作时非常有用。

指定GPU及文件保存地址

os.environ["CUDA_VISIBLE_DEVICES"] = '0'
torch.set_default_tensor_type(torch.cuda.FloatTensor)

## just in case, please create the directory before runnning
directory = 'results'

# list for AdaGrad
dy_squared  = []
dy_squared.append(None)
dy_squared.append(None)

torch.set_default_tensor_type(torch.cuda.FloatTensor)：

• 这行代码将PyTorch的默认张量类型设置为torch.cuda.FloatTensor，即GPU上的浮点张量类型。
• 设置默认张量类型为GPU上的浮点张量类型可以在不指定张量类型的情况下，确保大多数张量在GPU上进行计算，从而充分利用GPU的并行计算能力，加速运算。

directory = ‘results’：

• 这行代码创建一个名为"results"的目录，并将其赋值给变量directory。
• 如果代码中要使用名为"results"的目录，需要在使用之前确保该目录已经被创建。否则，在尝试写入或读取该目录时，会出现目录不存在的错误

dy_squared = []：

• 这行代码创建了一个空列表dy_squared。
• 列表是Python中的一种数据结构，用于存储一系列元素。在这里，dy_squared列表被用于存储AdaGrad算法的相关数据。

dy_squared.append(None)和dy_squared.append(None)：

• 这两行代码将None添加到dy_squared列表中两次。
• 根据代码的注释“list for AdaGrad”，这个列表是用来存储AdaGrad算法的数据的。dy_squared列表中的None元素可能是用作占位符，以便稍后将其他数据添加到列表中。

**分析：**这段代码的作用是设置PyTorch在GPU上运行的环境，创建一个名为"results"的目录，并初始化一个用于AdaGrad算法的列表。在后续的代码中，dy_squared列表可能会被用于存储AdaGrad算法的中间结果或相关信息。

初始化神经网络权重与偏置

def initialize_weights(node_sizes):
  w, b = [],[]
  for i in range(len(node_sizes)-1):
    w.append(torch.tensor(np.random.rand(node_sizes[i],node_sizes[i+1]) * np.sqrt(6. /(node_sizes[i]*node_sizes[i+1]))).float().cuda())
    b.append(torch.tensor(np.random.normal(size=node_sizes[i+1]) - 0.5).float().cuda())
  return w, b

主要步骤：

1. 创建空的权重列表w和偏置列表b。
2. 使用for循环遍历node_sizes列表（除了最后一个元素）的索引。
3. 在循环中，对于每个索引i，执行以下操作：a. 生成一个随机数矩阵作为权重，其形状为**(node_sizes[i], node_sizes[i+1])。使用np.random.rand生成0到1之间的随机数，乘以一个标量np.sqrt(6. /(node_sizes[i]*node_sizes[i+1]))进行缩放。然后将其转换为torch.tensor类型，并使用.float()将其转换为浮点类型。最后，通过调用.cuda()将其移动到GPU上（如果可用）。b. 生成一个正态分布的随机数作为偏置，其形状为node_sizes[i+1]。使用np.random.normal生成正态分布随机数，减去0.5以进行平移。然后将其转换为torch.tensor类型，并使用.float()将其转换为浮点类型。最后，通过调用.cuda()将其移动到GPU上（如果可用）。c. 将权重和偏置分别添加到权重列表w和偏置列表b**中。
4. 循环结束后，返回权重列表w和偏置列表b作为函数的输出。

数据集预处理

def preprocess_data():
  train_size = 60000
  test_size  = 10000  
  (x_train, y_train), (x_test, y_test) = mnist.load_data()
  
  x_train = np.reshape(x_train, (train_size, -1)) / 255.0
  x_test  = np.reshape(x_test,  (test_size, -1)) / 255.0

  # one hot encoding
  enc = OneHotEncoder()
  train_y = enc.fit_transform(y_train[:, np.newaxis]).toarray()
  test_y = enc.fit_transform(y_test[:, np.newaxis]).toarray()

  return torch.from_numpy(x_train), torch.from_numpy(x_test), torch.from_numpy(train_y), torch.from_numpy(test_y)

shuffle函数

def shuffle_data(X, y):
  return utils.shuffle(X, y, random_state=1234)

神经元动态计算模块（动力学计算，返回每层神经元的状态）

## calculating dynamics
## return states of neurons at each layer
def calculate_dynamics(input, time, delay, dt, batch_size, node_sizes, gamma, w, b, target=None, binarize=False):
    # initialization of activations
    activations     = [torch.zeros((batch_size,size)) for size in node_sizes]
    activations_new = [torch.zeros((batch_size,size)) for size in node_sizes]
     
    # clamped input
    activations[0] =  input
    length = len(node_sizes)

    store_all_activations     = [torch.zeros((size, batch_size, time)) for size in node_sizes]

    # similations for the free or clamped phase start
    for t in range(time):
      length = len(node_sizes)
      if (target is not None) and (t >= delay):
        activations[-1] = target
        length -= 1
      for j in range(1,length):
        if (target is None and j == length -1) or (t < delay and target is not None and j == length -1)  : # length + 1 does not exist
            activations_new[j] =  activations[j] + dt * (- activations[j] + sigmoid(torch.mm(activations[j-1].float().cuda(), w[j-1])  + b[j-1]))
        else: 
           activations_new[j] =  activations[j] + dt * (- activations[j] + sigmoid(torch.mm(activations[j-1].float().cuda(), w[j-1])  + gamma*torch.mm(activations[j+1].float().cuda(), torch.transpose(w[j], 0, 1))  + b[j-1]))

      # t -> t +1
      for k in range(1, length):
        activations[k] = activations_new[k]
        store_all_activations[k][:,:, t] = torch.transpose(activations[k], 0, 1)

    return activations, store_all_activations

1. 输入参数：

• input：输入数据的张量，形状为**(batch_size, input_size)，其中batch_size**表示批量大小，input_size表示输入层的节点数量。
• time：模拟的时间步数。
• delay：延迟的时间步数。
• dt：时间步长。
• batch_size：批量大小。
• node_sizes：包含每层节点数量的整数列表。
• gamma：用于计算下一层到当前层的反馈权重的因子。
• w：权重列表，包含每层之间的连接权重。
• b：偏置列表，包含每层的偏置项。
• target：目标输出的张量，形状为**(batch_size, output_size)，默认为None**，表示没有目标输出。
• binarize：一个布尔值，指示是否对神经元的输出进行二值化，默认为False。

2. 函数的返回值：

两个列表：

• activations：包含每层神经元状态的列表。每个元素是一个形状为**(batch_size, size)**的张量，表示对应层的神经元状态。
• store_all_activations：包含每层神经元状态的张量的列表。每个元素是一个形状为**(size, batch_size, time)**的张量，表示对应层神经元在每个时间步的状态。

3. 神经元状态列表（重点）

• activations是一个包含每层神经元状态张量的列表，初始状态为全零。
• activations_new是一个与activations具有相同结构的列表，用于存储更新后的神经元状态。

4. 创建用于存储每层神经元状态的张量列表

5. 随时间步的神经动力学循环（难点）

 # similations for the free or clamped phase start
    for t in range(time):
      length = len(node_sizes)
      if (target is not None) and (t >= delay):
        activations[-1] = target  #将目标输出值赋值给最后一层神经元的激活状态activations[-1]
        length -= 1
      for j in range(1,length): #除了输入层和输出层的每一层神经元
        if (target is None and j == length -1) or (t < delay and target is not None and j == length -1)  : # length + 1 does not exist
            activations_new[j] =  activations[j] + dt * (- activations[j] + sigmoid(torch.mm(activations[j-1].float().cuda(), w[j-1])  + b[j-1]))
        else: 
           activations_new[j] =  activations[j] + dt * (- activations[j] + sigmoid(torch.mm(activations[j-1].float().cuda(), w[j-1])  + gamma*torch.mm(activations[j+1].float().cuda(), torch.transpose(w[j], 0, 1))  + b[j-1]))

• for t in range(time):：通过循环来模拟神经网络的动力学，执行 time 个时间步。
• length = len(node_sizes)：获取神经网络的层数，即神经网络中的层的数量。
• if (target is not None) and (t >= delay):：这是一个条件判断语句，检查是否存在目标输出 target 并且当前时间步 t 大于或等于 delay。
  • 如果满足条件，表示模拟处于夹持阶段。将目标输出 target 赋值给最后一层神经元的激活状态 activations[-1]。这样，模拟时最后一层的神经元状态将被强制设为目标输出值。
  • 同时将length减去1，因为在夹持阶段不需要更新最后一层的状态。
• for j in range(1, length):：通过循环遍历除了输入层和输出层的每一层神经元，即遍历隐藏层。
  • if (target is None and j == length - 1) or (t < delay and target is not None and j == length - 1):：这是一个条件判断语句，检查当前层是否为输出层。
    • 如果满足条件，表示当前时间步处于自由阶段或夹持阶段的输出层（length + 1不存在）。在这种情况下，使用更新规则更新当前层的状态：
      • 计算当前层和前一层的加权输入，使用torch.mm函数执行矩阵乘法，然后加上偏置项 b[j-1]。
      • 通过 sigmoid 函数作为激活函数，计算更新后的状态。

神经元更新函数：

• 如果不满足上述条件，表示当前层为隐藏层。在这种情况下，使用更新规则更新当前层的状态：
  • 计算当前层和前一层的加权输入，使用 torch.mm 函数执行矩阵乘法，然后加上偏置项 b[j-1]。
  • 同时，还考虑下一层到当前层的反馈（反向权重）：gamma * torch.mm(activations[j+1].float().cuda(), torch.transpose(w[j], 0, 1))。
  • 最后，通过 sigmoid 函数作为激活函数，计算更新后的状态。

神经元更新函数：

• 更新当前层的状态 activations_new[j]。
• 在循环结束后，返回更新后的神经元状态列表 activations 和存储所有时间步的神经元状态列表 store_all_activations。

6. 时间步循进

# t -> t +1
      for k in range(1, length):
        activations[k] = activations_new[k]
        store_all_activations[k][:,:, t] = torch.transpose(activations[k], 0, 1)

用于在时间步 t 到 t + 1 时更新神经网络每个层的激活状态，并将当前时间步的激活状态记录在 store_all_activations 列表中的相应张量中。

7. 为什么前面两种判断条件下一个需要对w转置一个不需要？

对于隐藏层（非输出层）：

• 输入数据 activations[j-1] 的形状是 (batch_size, size_prev_layer)，其中 size_prev_layer 表示前一层的神经元数量。
• 权重矩阵 w[j-1] 的形状是 (size_prev_layer, size)，其中 size 表示当前层的神经元数量。
• 在矩阵乘法 torch.mm(activations[j-1].float().cuda(), w[j-1]) 中，前一层的输出作为矩阵乘法的左操作数，而权重矩阵作为右操作数。这样的矩阵乘法是合法的，因为左操作数的列数等于右操作数的行数，即 (batch_size, size_prev_layer) x (size_prev_layer, size)，结果得到形状 (batch_size, size)。
• 在这种情况下，不需要对权重矩阵进行转置。

对于输出层：

• 输入数据 activations[j-1] 的形状是 (batch_size, size_prev_layer)，其中 size_prev_layer 表示前一层的神经元数量。
• 权重矩阵 w[j-1] 的形状是 (size, size_prev_layer)，这里与隐藏层不同，size 表示当前层的神经元数量。
• 在矩阵乘法 torch.mm(activations[j-1].float().cuda(), w[j-1]) 中，前一层的输出作为矩阵乘法的左操作数，而权重矩阵作为右操作数。这样的矩阵乘法要求左操作数的列数等于右操作数的行数，即 (batch_size, size_prev_layer) x (size, size_prev_layer)，结果得到形状 (batch_size, size)。
• 在这种情况下，由于权重矩阵 w[j-1] 的形状与隐藏层不同，需要将其转置为 (size_prev_layer, size)，使得矩阵乘法的维度匹配。

8. dt的作用：

dt 表示时间步长，用于调整在每个时间步内更新神经元状态的幅度。乘以 dt 的作用是控制动力学模拟的步幅大小，从而影响神经元状态在每个时间步中的变化。
具体来说，神经网络的动力学模拟涉及到在每个时间步更新神经元的状态，以模拟网络在不同时间点的响应行为。在更新神经元状态时，使用了类似欧拉法（Euler’s method）的数值积分方法。这种方法通过对微分方程进行离散化，将微分方程转换为差分方程，然后通过在每个时间步上进行差分运算来更新状态。
在这里，乘以 dt 的作用是将微分方程转换为差分方程时引入的步长因子。通过调整 dt 的值，可以控制每个时间步中神经元状态的变化速度。较小的 dt 值会使得状态更新更加细致和平滑，因为每个时间步的更新幅度较小。相反，较大的 dt 值会导致状态更新较快，因为每个时间步的更新幅度较大。
总的来说，dt 的取值会影响神经网络动力学模拟的精度和速度。较小的 dt 值可能会更准确地模拟网络的动态行为，但需要更多的计算时间。较大的 dt 值可以加快计算速度，但可能会牺牲一些模拟的准确性。因此，在使用这段代码时，根据模拟的需要和计算资源的限制，可以适当调整 dt 的值。

突触权重（与偏置）更新模块（AdaGrad with Contrastive Hebbian Learning (AdaCHL)算法）

1. AdaGrad(Adaptive Gradient Algorithm, 自适应梯度算法）

• AdaGrad是一种自适应学习率算法，旨在根据每个参数在训练过程中的历史梯度来自动调整学习率。
• 对于每个参数，AdaGrad维护一个梯度累积的平方和，用于动态调整学习率。梯度累积的平方和会随着训练的进行而增大，导致学习率逐渐减小。这样，较少更新的参数会得到更大的学习率，而较频繁更新的参数会得到较小的学习率，从而在训练初期更大幅度地更新参数，在训练后期则更小幅度地更新参数，从而更好地适应训练数据。

2. CHL（Contrastive Hebbian Learning, 对比Hebbian学习）

• 对比性海比学习是一种基于突触权重的学习规则，旨在模拟生物神经元之间的突触学习过程。
• 在对比性海比学习中**，突触权重根据突触前后神经元的激活状态差异进行调整**。如果两个神经元同时激活，说明它们之间可能存在相关性，此时增强对应的突触权重；反之，如果两个神经元的激活状态差异较大，则减弱对应的突触权重。

AdaGrad with Contrastive Hebbian Learning (AdaCHL) 将以上两种方法结合起来，实现对神经网络权重的优化。在AdaCHL算法中，对于每层的权重更新，首先计算梯度变化 dy，然后使用AdaGrad的思想调整学习率，并结合对比性海比学习的概念，根据当前层与前一层之间的激活状态差异对权重进行更新。

**Pytorch代码实现：

**

def update_weights(w, b, learning_rate, gamma, free_act, clamped_act, length=2, batch_size=32):
    for i in range(1, length):       
      
      # AdaptiveGrad with Contrastive Hebbian Learning (AdaCHL)
      # Reffered to AdaGrad
      #首先计算梯度变化dy
      dy =  (torch.mm(torch.transpose(clamped_act[i-1].float().cuda(), 0, 1), clamped_act[i].float().cuda())  - torch.mm(torch.transpose(clamped_act[i-1].float().cuda(), 0, 1), free_act[i].float().cuda()))/float(batch_size)
      global dy_squared  #声明dy_squared是全局变量
      if dy_squared[i-1] is None:
        dy_squared[i-1] = dy * dy
      else:
        dy_squared[i-1] += dy * dy  #历史梯度的平方和
    
      dy_update  = dy/(torch.sqrt(dy_squared[i-1]) + 1e-7)  #自适应调整学习率因子
      
      w[i-1] += learning_rate[i-1]*(dy_update) #更新权重
      b[i-1] += learning_rate[i-1]*((clamped_act[i].float().cuda() -  free_act[i].float().cuda())[0])/float(batch_size) 
      #[0] 表示取张量的第一个元素，这里之所以用 [0] 是因为偏置项 b[i-1] 的形状通常是 (size,)，即一维的向量。
    return w, b

1. 计算梯度变化 dy：通过夹持阶段的激活状态 clamped_act 与自由阶段的激活状态 free_act 之间的差异得到 dy。
2. AdaGrad调整学习率：将 dy 除以历史梯度的平方和来适应性地调整学习率，使得对梯度较大的参数使用较小的学习率。

（自适应调整学习率缩放比率）

3. 结合对比Hebbian学习：根据 dy 调整后的学习率，对当前层与前一层之间的连接权重 w 进行更新。

论文给的原方程：

4. 更新偏置项：同时根据梯度变化 dy 和学习率对当前层的偏置项 b 进行更新。

5. global dy_squared 的作用是声明 dy_squared 是一个全局变量，从而使得在函数内部的代码中可以访问并修改这个全局变量的值。

当在函数内部使用变量时，默认情况下，Python会将其视为局部变量，即使全局变量具有相同的名称。如果你想在函数内部修改全局变量的值，必须使用 global 关键字来明确告诉Python该变量是全局变量。
在这段代码中，dy_squared 是一个列表，用于存储每一层的 dy * dy 的累积和。在函数 update_weights 中，首先通过 global dy_squared 声明了 dy_squared 是一个全局变量，然后在后续的代码中可以对 dy_squared 进行读取和修改操作。
具体来说，dy_squared 列表的每个元素 dy_squared[i-1] 存储了第 i 层的 dy * dy 的累积和。在每次计算 dy 的更新时，会根据是否为第一次计算来决定是初始化 dy_squared[i-1] 还是累加到原有值上。
全局变量 dy_squared 的使用是为了在多次调用 update_weights 函数时保持 dy_squared 的状态，确保能够正确地累积每层的 dy * dy 的值，并在每次更新权重时使用该累积值来计算学习率调整，这是 AdaGrad 算法的核心机制。

神经元未来状态预测函数（最小二乘线性回归模型）

def predict_dynamics(store_free_all_activations, prediction_inp_size, node_size, update_data_idx, train_ls_idx, length=5):
  #length = np.shape(store_free_all_activations)[0]
  pred_all_activations = []
  
  # for each layer
  for i in range(1,length):
    layer_l     = store_free_all_activations[i]
    node_size_l = np.shape(layer_l)[0] 
    activaitons = []
    
    # for each neuron
    for j in range(node_size_l):
      ## training data for the prediction
      one_neuron_train_data     =  layer_l[j, train_ls_idx, :prediction_inp_size].cpu().numpy()
      
      ## testing data (prediction data)
      one_neuron_test_data      = layer_l[j, update_data_idx, :prediction_inp_size].cpu().numpy()
  
      shape_train                      = np.shape(one_neuron_train_data)
      shape_test                       = np.shape(one_neuron_test_data)
      
      ## adding offset for trainig and prediction data
      one_neuron_input_offset_train    = np.ones((shape_train[0], shape_train[1]+1))
      one_neuron_input_offset_test     = np.ones((shape_test[0], shape_test[1]+1))     
      one_neuron_input_offset_train[:, :-1] = one_neuron_train_data
      one_neuron_input_offset_test[:, :-1]  = one_neuron_test_data
      
      ## targets for traininig 
      one_neuron_train_target               = layer_l[j, train_ls_idx, -1].cpu().numpy()
      
      ## training for linear regresssion
      pred_activation                       = np.linalg.lstsq(one_neuron_input_offset_train, one_neuron_train_target, rcond=None)[0]
      
      # prediction
      #pred_negative_acts              = one_neuron_input_offset_test @ pred_activation # or
      pred_negative_acts              = np.dot(one_neuron_input_offset_test,  pred_activation)
      
      # if values are negative, they will be 0
      pred_negative_acts              = np.clip(pred_negative_acts, a_min=0, a_max=None)    

      activaitons.append(torch.from_numpy(pred_negative_acts))
    pred_all_activations.append(torch.transpose(torch.stack(activaitons), 0,1))
  
  return pred_all_activations

函数输入参数

• store_free_all_activations：自由阶段所有神经元的激活状态存储的列表。
• prediction_inp_size：预测输入的大小（数据点的数量）。
• node_size：每层的神经元数量。
• update_data_idx：用于预测的测试数据的索引。
• train_ls_idx：用于训练的数据的索引。
• length：神经网络的层数，默认为5层。

备注：

layer_l	->	[node_size_l, num_samples, time_steps]
#node_size_l	:当前层的神经元数量
#num_samples	:样本数(批量大小）
#time_steps		:时间步数

执行步骤

1. 初始化 pred_all_activations 列表，用于存储每一层的预测激活状态。
2. 对于每一层（从第1层到第length-1层）：
   • 获取该层的激活状态 layer_l，并获取该层的神经元数量 node_size_l。
   • 对于每个神经元（从第0个到第node_size_l-1个）：

这句代码的作用是从神经网络某一层 layer_l 的第 j 个神经元中提取用于训练的数据，包含 train_ls_idx 所指示的样本，并取这些样本前 prediction_inp_size 个时间步的激活状态，并将其转换为NumPy数组，以备后续的线性回归训练。

  - 获取训练数据和测试数据，用于线性回归预测。
  - 对训练数据和测试数据进行偏移处理（添加偏置项）。（下文解释，难点）
  - 获取训练目标（最后一个时间步的激活状态）。
  - 使用线性回归（最小二乘法）对训练数据和训练目标进行拟合，得到预测用的线性回归系数 **pred_activation**。
  - 对测试数据进行预测，得到未来时间步的预测激活状态 **pred_negative_acts**。
  - 将预测激活状态中小于0的值截断为0，确保预测的状态为非负值。
  - 将每个神经元的预测结果添加到 **activaitons** 列表中。

3. 将每一层的所有神经元的预测激活状态（即 activaitons 列表）通过 torch.stack 和 torch.transpose 进行整理，并添加到 pred_all_activations 列表中。
4. 返回 pred_all_activations，即每一层未来时间步的预测激活状态。
   

训练与预测神经元传入数据区分开来：

引入偏移量（小难点）

one_neuron_input_offset_train = np.ones((shape_train[0], shape_train[1]+1))：创建一个由全为1的元素组成的数组，其形状为 (shape_train[0], shape_train[1]+1)。这里 shape_train[0] 表示训练样本的数量，shape_train[1] 表示训练样本每个时间步的维度（特征数），shape_train[1]+1 是在每个样本的最后一列添加一个偏置项的位置。
one_neuron_input_offset_test = np.ones((shape_test[0], shape_test[1]+1))：类似地，创建一个由全为1的元素组成的数组，形状为 (shape_test[0], shape_test[1]+1)。这里 shape_test[0] 表示测试（预测）样本的数量，shape_test[1] 表示测试（预测）样本每个时间步的维度（特征数），shape_test[1]+1 是在每个样本的最后一列添加一个偏置项的位置。
one_neuron_input_offset_train[:, :-1] = one_neuron_train_data：将 one_neuron_train_data 的数据（除了最后一列）复制到 one_neuron_input_offset_train 的每一行中，这样，one_neuron_input_offset_train 中的前 shape_train[1] 列将和 one_neuron_train_data 中的数据相同。
one_neuron_input_offset_test[:, :-1] = one_neuron_test_data：类似地，将 one_neuron_test_data 的数据（除了最后一列）复制到 one_neuron_input_offset_test 的每一行中，这样，one_neuron_input_offset_test 中的前 shape_test[1] 列将和 one_neuron_test_data 中的数据相同。
作用：通过以上步骤，one_neuron_input_offset_train 和 one_neuron_input_offset_test 现在包含了相应的训练和测试（预测）数据，并在每个样本的最后一列添加了全为1的偏置项，以备后续进行线性回归（最小二乘法）计算。这样一来，在线性回归过程中，就可以同时学习合适的斜率（系数）和截距（偏置），以更好地拟合数据。

获取用于训练的目标值（label）

• layer_l[j]: 选择神经网络中的第 i 层的第 j 个神经元的所有时间步的激活状态。
• layer_l[j, train_ls_idx]: 从上一步得到的结果中选择用于训练的样本，其中 train_ls_idx 表示用于训练的数据的索引，即选择的样本的索引。
• layer_l[j, train_ls_idx, -1]: 从上一步得到的结果中选择每个训练样本的最后一个时间步的激活状态，这就是用于训练的目标值（标签）。
  

线性回归模型训练（线性回归模型参数）

1. 线性回归模型：

2. 对应代码：

np.linalg.lstsq()函数：
np.linalg.lstsq 来计算线性回归模型的参数。该函数接收两个参数，一个是输入数据 one_neuron_input_offset_train，一个是目标值 one_neuron_train_target。rcond 参数用于控制奇异矩阵的容差，默认值为 None，表示使用默认容差。
[0]: 上一步得到的结果是一个元组，其中包含了多个返回值。而我们对线性回归的结果只关心第一个返回值，即线性回归的参数。通过 [0] 取得这个参数值。（原因不明）

3. rcond参数介绍（了解）：

在线性代数中，rcond（相对条件数）是用来评估矩阵的条件数的一个指标。条件数是用来衡量矩阵的稳定性和数值精度的一个度量。具体来说，条件数越大，矩阵越接近奇异（singular），也就是说，矩阵的行列式接近于零，它在数值计算中容易出现问题，例如出现数值不稳定或舍入误差等。

线性回归运算

负值置零（np.clip()函数)

• pred_negative_acts: 这是一个NumPy数组，其中包含了线性回归模型对测试数据的预测值。
• np.clip(pred_negative_acts, a_min=0, a_max=None): 这是使用 NumPy 提供的 np.clip() 函数。该函数用于对数组中的元素进行限制（截断）操作，将超出指定范围的值进行裁剪。
• a_min=0: 这是 np.clip() 函数的参数，表示截断的下限。所有小于等于 a_min 的元素将被截断为 a_min。
• a_max=None: 这是 np.clip() 函数的参数，表示截断的上限。所有大于等于 a_max 的元素将被截断为 a_max。由于 a_max 设置为 None，意味着没有上限，即不对大于 a_min 的值进行截断。

Usage: 通过这样的操作，将线性回归模型预测得到的负值都截断为零，确保了预测结果只包含非负值。这在某些应用场景中是很有用的，例如对于预测非负的实数值或概率值时。(神经元刺激状态最小为0）

保存预测神经元状态

1. activaitons.append(torch.from_numpy(pred_negative_acts))：在每个循环迭代中，将预测得到的神经元激活状态 pred_negative_acts 转换为 PyTorch 的张量（tensor）对象，并将其添加到名为 activaitons 的列表中。
2. pred_all_activations.append(torch.transpose(torch.stack(activaitons), 0, 1))：在每个循环迭代中，将 activaitons 列表中的元素堆叠（stack）起来，得到一个新的张量。然后，通过 torch.transpose() 函数进行张量的转置，将不同神经元的预测结果排列在张量的第一维上，不同时间步的预测结果排列在张量的第二维上。最后，将这个转置后的张量添加到名为 pred_all_activations 的列表中。

准确率计算模块

def check_accuracy(data_x, data_y, data_size, batch_size, node_sizes, free_time, clamped_time, delay, dt, gamma, w, b, n_activations):
  accs = []
  test_size = 10000 
  index = int(test_size/float(batch_size))
  for i in range(index): 
    x = torch.reshape(data_x[i*batch_size:(i+1)*batch_size], (batch_size, node_sizes[0]))
    y = torch.reshape(data_y[i*batch_size:(i+1)*batch_size], (batch_size, node_sizes[-1]))
      
    free_act, _          = calculate_dynamics(x, free_time, delay, dt, batch_size, node_sizes, gamma, w, b)


    acc =  torch.argmax(free_act[-1].float().cuda(), dim=1) == torch.argmax(y.float().cuda(), dim=1)  
    accs.append(acc)
  return torch.mean(torch.stack(accs).float().cuda())

输入参数

1. check_accuracy: 这是函数的名称，用于检查神经网络在给定数据上的准确性。
2. data_x: 输入数据（特征），是一个PyTorch张量，包含用于测试的输入样本。
3. data_y: 输出数据（标签），是一个PyTorch张量，包含与输入数据对应的真实标签。
4. data_size: 数据集大小，即总样本数量。
5. batch_size: 批处理大小，用于指定每个批次中的样本数量。
6. node_sizes: 包含神经网络每一层神经元数量的列表。
7. free_time: 自由相的时间步数，用于指定自由相的时间长度。
8. clamped_time: 固定相的时间步数，用于指定固定相的时间长度。
9. delay: 固定相后到自由相开始之间的时间延迟。
10. dt: 时间步长。
11. gamma: 反馈强度参数。
12. w: 包含每层神经元权重的列表。
13. b: 包含每层神经元偏置的列表。
14. n_activations: 记录每层神经元状态的时间步数。
    

输入数据与标签数据提取

返回神经元状态张量列表（见上文）

计算准确率

1. torch.argmax(free_act[-1].float().cuda(), dim=1): 这部分代码首先将 free_act 中的最后一个时间步的激活状态提取出来，并执行以下操作：
   • free_act[-1]: 提取 free_act 最后一个时间步的激活状态。在这里，free_act 是一个包含每个时间步激活状态的列表，[-1] 表示获取列表的最后一个元素，即最后一个时间步的激活状态。
   • .float(): 将激活状态转换为浮点数类型，以便进行后续的计算。
   • .cuda(): 将激活状态移动到GPU上，如果GPU可用的话。
2. torch.argmax(y.float().cuda(), dim=1): 这部分代码用于计算真实标签数据 y 在每个样本中的最大值索引，并执行以下操作：
   • y.float(): 将真实标签数据转换为浮点数类型，以便进行后续的计算。
   • .cuda(): 将真实标签数据移动到GPU上，如果GPU可用的话。
3. ==: 这是一个逻辑运算符，用于比较两个张量的元素是否相等。它会对两个张量进行逐元素比较，并返回一个布尔类型的张量，其中元素为True表示对应位置的元素相等，元素为False表示对应位置的元素不相等。
4. acc = torch.argmax(free_act[-1].float().cuda(), dim=1) == torch.argmax(y.float().cuda(), dim=1): 这行代码将上述两个逻辑张量相等性比较的结果赋值给变量 acc。acc 是一个布尔类型的张量，它记录了在给定输入数据 x 下神经网络的自由相激活状态与真实标签数据 y 的对应位置是否相等，即记录了预测结果是否正确。
   

模型训练模块（重难点）

源码

def train_model(epoch, w, b, learning_rate, gamma, batch_size, minibatch_size, free_time, clamped_time, delay, dt, node_sizes, n_activations, prediction_inp_size, pred=False):
  train_accs   = []
  test_accs    = []
  train_x, test_x, train_y, test_y = preprocess_data()
  train_size = np.shape(train_x)[0]
  test_size  = np.shape(test_x)[0]
  epoch_train_size = int(train_size/batch_size)
  
  ## check accuracies before training   
  ## check training accuracy (mean)
  train_acc  = check_accuracy(train_x, train_y, train_size, batch_size, node_sizes, free_time, clamped_time, delay, dt, gamma, w, b, n_activations)      

  ## check testing accuracy (mean)
  test_acc   = check_accuracy(test_x, test_y, test_size, batch_size, node_sizes, free_time, clamped_time, delay, dt, gamma, w, b, n_activations)

  print('epoch:' + str(0))
  print('accuracy for training: ' + str(train_acc.cpu().numpy()))
  print('accuracy for testing: '  + str(test_acc.cpu().numpy()))
      
  f = None
  if os.path.isfile(directory + '/log.txt'):
    f = open(directory + '/log.txt', 'a')
  else:
    os.mkdir(directory)    
    f = open(directory + '/log.txt', 'w')

  np.save(directory + '/w_epoch_' + str(0) + '.npy', w)    
  np.save(directory + '/b_epoch_' + str(0) + '.npy', b)        
  np.save(directory + '/dy_squared_epoch_' + str(0) + '.npy', dy_squared)
        
  f.write("Epoch: " + str(0) + '\n')
  f.write("accuracy for training: " + str(train_acc.cpu().numpy()) + '\n')
  f.write("accuracy for testing: " + str(test_acc.cpu().numpy()) + '\n')
  f.close()  
  for i in range(epoch):

    start = time.time()
    
    train_x, train_y = shuffle_data(train_x, train_y)

    for j in range(epoch_train_size):
      one_x = torch.reshape(train_x[j*batch_size:(j+1)*batch_size], (batch_size, node_sizes[0]))
      one_y = torch.reshape(train_y[j*batch_size:(j+1)*batch_size], (batch_size, node_sizes[-1]))
      
      ## the free phase
      free_act,    store_free_all_activations    = calculate_dynamics(one_x, free_time, delay, dt, batch_size, node_sizes, gamma, w, b)
      
      ## randomly picked up data indices for the prediction and clamped phase data to update the weights
      update_data_idx                            = np.random.choice(batch_size, size=minibatch_size, replace=False)
      
      ## these indices are for training LS model to predict the activations  
      train_ls_idx                               = [k for k in range(batch_size) if k not in update_data_idx]

      if pred:
          ## predict the dynamics for both hidden and output
          free_pred_acts = predict_dynamics(store_free_all_activations, prediction_inp_size, node_sizes, update_data_idx, train_ls_idx, length=len(node_sizes))
          
          ## store predicted dynamics into the free phase activations
          input_act = free_act[0]
          del free_act
          
          free_act = []
          free_act.append(input_act[update_data_idx, :])  
          free_act.append(free_pred_acts[0])
          free_act.append(free_pred_acts[1])

      ## the clamped phase 
      clamped_act, store_clamped_all_activations = calculate_dynamics(one_x[update_data_idx,:], clamped_time, delay, dt, minibatch_size, node_sizes, gamma, w, b, target=one_y[update_data_idx,:])
            
      ## update the weights
      w, b = update_weights(w, b, learning_rate, gamma, free_act, clamped_act, length=len(node_sizes), batch_size=minibatch_size)

    ## after every xxx epoch, check the accuracies for training and testing  
    if i % 1 == 0:
      ## check training accuracy (mean)
      train_acc  = check_accuracy(train_x, train_y, train_size, batch_size, node_sizes, free_time, clamped_time, delay, dt, gamma, w, b, n_activations)      

      ## check testing accuracy (mean)
      test_acc   = check_accuracy(test_x, test_y, test_size, batch_size, node_sizes, free_time, clamped_time, delay, dt, gamma, w, b, n_activations)
    
      print('epoch:' + str(i+1))
      print('accuracy for training: ' + str(train_acc.cpu().numpy()))
      print('accuracy for testing: '  + str(test_acc.cpu().numpy()))

      np.save(directory + '/w_epoch_' + str(i+1) + '.npy', w)    
      np.save(directory + '/b_epoch_' + str(i+1) + '.npy', b)         
      
      f = None
      if os.path.isfile(directory + '/log.txt'):
        f = open(directory + '/log.txt', 'a')
      else:
        f = open(directory + '/log.txt', 'w')
        
      f.write("Epoch: " + str(i+1) + '\n')
      f.write("accuracy for training: " + str(train_acc.cpu().numpy()) + '\n')
      f.write("accuracy for testing: " + str(test_acc.cpu().numpy()) + '\n')
      f.close()        

      train_accs.append(train_acc)
      test_accs.append(test_acc)
      end = time.time()
      print(end - start)

  return train_accs, test_accs

数据集预处理

train_x, test_x, train_y, test_y = preprocess_data()  #数据预处理

结果记录文档保存操作

f = None
  if os.path.isfile(directory + '/log.txt'):
    f = open(directory + '/log.txt', 'a')
  else:
    os.mkdir(directory)    
    f = open(directory + '/log.txt', 'w')

  np.save(directory + '/w_epoch_' + str(0) + '.npy', w)    
  np.save(directory + '/b_epoch_' + str(0) + '.npy', b)        
  np.save(directory + '/dy_squared_epoch_' + str(0) + '.npy', dy_squared)
        
  f.write("Epoch: " + str(0) + '\n')
  f.write("accuracy for training: " + str(train_acc.cpu().numpy()) + '\n')
  f.write("accuracy for testing: " + str(test_acc.cpu().numpy()) + '\n')
  f.close()  

在指定文件夹内新建log.txt文本文件，若文件已存在，则设定为追加模式（a)，否则设定为写入模式（w)赋值给f。

shuffle输入数据

训练过程

  for i in range(epoch):

    start = time.time()  #记录本次epoch的训练时间
    
    train_x, train_y = shuffle_data(train_x, train_y)

    for j in range(epoch_train_size):
      one_x = torch.reshape(train_x[j*batch_size:(j+1)*batch_size], (batch_size, node_sizes[0]))  #输入端size
      one_y = torch.reshape(train_y[j*batch_size:(j+1)*batch_size], (batch_size, node_sizes[-1])) #输出端size
      
      ## the free phase
      free_act,    store_free_all_activations    = calculate_dynamics(one_x, free_time, delay, dt, batch_size, node_sizes, gamma, w, b)
      
      ## randomly picked up data indices for the prediction and clamped phase data to update the weights
      update_data_idx                            = np.random.choice(batch_size, size=minibatch_size, replace=False)
      
      ## these indices are for training LS model to predict the activations  
      train_ls_idx                               = [k for k in range(batch_size) if k not in update_data_idx]

返回当前自由相神经元状态

## the free phase  (返回当前自由相神经元状态）
      free_act,    store_free_all_activations    = calculate_dynamics(one_x, free_time, delay, dt, batch_size, node_sizes, gamma, w, b)

随机抽取子样本用于预测未来神经元状态（重点）

## randomly picked up data indices for the prediction and clamped phase data to update the weights
      update_data_idx                            = np.random.choice(batch_size, size=minibatch_size, replace=False)
      
      ## these indices are for training LS model to predict the activations  
      train_ls_idx                               = [k for k in range(batch_size) if k not in update_data_idx]

      if pred:
          ## predict the dynamics for both hidden and output
          free_pred_acts = predict_dynamics(store_free_all_activations, prediction_inp_size, node_sizes, update_data_idx, train_ls_idx, length=len(node_sizes))
          
          ## store predicted dynamics into the free phase activations
          input_act = free_act[0]
          del free_act

1. update_data_idx 是通过 np.random.choice 函数随机从 batch_size 个数据索引中选择 minibatch_size 个索引。这里 batch_size 表示一个批次中的样本数量，而 minibatch_size 表示每次选择的样本数量。这样做的目的是为了从当前批次中随机选择一部分样本，用于进行预测相（clamped phase）和权重更新。
2. train_ls_idx 则是在当前批次的所有数据索引中排除了 update_data_idx 中的索引。这样，train_ls_idx 就包含了那些没有被用于预测相的数据索引，用于训练线性回归模型来预测激活状态。
3. 当 pred 为 True 时，表示在预测相过程中，需要进行激活状态的预测。此时，调用 predict_dynamics 函数，传入相关参数 store_free_all_activations、prediction_inp_size、node_sizes、update_data_idx、train_ls_idx 以及 length，来预测隐藏层和输出层的激活状态。
4. 预测的激活状态 free_pred_acts 存储在变量中。
5. 最后，将 free_act 的第一个元素（即隐藏层的激活状态）存储在 input_act 中，并使用 del free_act 删除 free_act 的引用，释放内存。
   

自由相神经元状态列表更新

clamped phase（重点）

## the clamped phase 
      clamped_act, store_clamped_all_activations = calculate_dynamics(one_x[update_data_idx,:], clamped_time, delay, dt, minibatch_size, node_sizes, gamma, w, b, target=one_y[update_data_idx,:])

calculate_dynamics()函数返回的2个对象：

传入的参数：

1. one_x[update_data_idx, :]: 这是一个张量，表示从 one_x 中选取在 update_data_idx 中指定的行的所有列作为输入。它是在固定相阶段用于计算神经元活动状态的输入数据。
2. clamped_time: 这是一个整数，表示固定相阶段的时间步数，即在固定相阶段网络将进行多少个时间步的计算。
3. delay: 这是一个浮点数，表示神经元之间的延迟（synaptic delay），即信号在神经元之间传播的时间。
4. dt: 这是一个浮点数，表示时间步长（time step），即每个时间步的持续时间。
5. minibatch_size: 这是一个整数，表示每个批次中的样本数量，用于在固定相阶段更新权重。
6. node_sizes: 这是一个列表，表示每层网络中的神经元数量。
7. gamma: 这是一个浮点数，表示神经元的学习速率（learning rate）。
8. w: 这是一个列表，表示网络的权重。
9. b: 这是一个列表，表示网络的偏置项。
10. target=one_y[update_data_idx, :]: 这是一个可选参数，表示在固定相阶段的目标输出。在一些情况下，可以将固定相阶段的输出与目标输出进行比较，从而更好地指导权重的更新。

函数 calculate_dynamics 将根据传入的参数计算在固定相阶段的神经元活动状态 clamped_act，并返回该状态以及其他可能有用的信息，如 store_clamped_all_activations。store_clamped_all_activations 是一个列表，用于存储在固定相阶段的每个时间步的神经元活动状态，以备后续使用。整个过程是为了进行权重更新，从而使网络逐渐逼近所期望的输出。

AdaCHL突触权重更新（重点）

## update the weights
      w, b = update_weights(w, b, learning_rate, gamma, free_act, clamped_act, length=len(node_sizes), batch_size=minibatch_size)

def update_weights(w, b, learning_rate, gamma, free_act, clamped_act, length=2, batch_size=32):
    for i in range(1, length):       
      
      # AdaptiveGrad with Contrastive Hebbian Learning (AdaCHL)
      # Reffered to AdaGrad
      #首先计算梯度变化dy
      dy =  (torch.mm(torch.transpose(clamped_act[i-1].float().cuda(), 0, 1), clamped_act[i].float().cuda())  - torch.mm(torch.transpose(clamped_act[i-1].float().cuda(), 0, 1), free_act[i].float().cuda()))/float(batch_size)
      global dy_squared  #声明dy_squared是全局变量
      if dy_squared[i-1] is None:
        dy_squared[i-1] = dy * dy
      else:
        dy_squared[i-1] += dy * dy  #历史梯度的平方和
    
      dy_update  = dy/(torch.sqrt(dy_squared[i-1]) + 1e-7)  #自适应调整学习率因子
      
      w[i-1] += learning_rate[i-1]*(dy_update) #更新权重
      b[i-1] += learning_rate[i-1]*((clamped_act[i].float().cuda() -  free_act[i].float().cuda())[0])/float(batch_size) 
      #[0] 表示取张量的第一个元素，这里之所以用 [0] 是因为偏置项 b[i-1] 的形状通常是 (size,)，即一维的向量。
    return w, b

保存ACC数据与训练时长

运行函数

源码

## running the code
## return training and testing accuraices
##
def run():

  ## setting up the hyper parameters 
  gamma          = 1.0                    # gamma for CHL 
  free_time      = 120                    # total simulation time for CHL
  clamped_time   = 120
  dt             = 0.1                    # time step for CHL
  epoch          = 3                      # training epoch
  learning_rate  = [0.03, 0.02]           # learning rate
  n_activations  = 100                    # how many examples should be used for activation maps
  pred           = True

  batch_size          = 500               # mini-batch size 
  minibatch_size      = 10                # update mini-batch size
  hidden_size         = [1000]            # hidden size for the network
  output_size         = 10                # output size for the network
  input_size          = 784               # input size for the network 28x28 = 784 mnist
  prediction_inp_size = 12                # how many inputs should be used for predictions of dynamics on the free phase
  delay               = 13                # time delay for the clamped phase 

  print('Delay: ' + str(delay))
  print('Pred inp size: ' + str(prediction_inp_size))  

  # all of the sizes for the network
  # at first, adding the first layer
  node_sizes = [input_size] 

  # adding the hidden layers 
  for i in hidden_size:
    node_sizes.append(i)

  # adding output layers
  node_sizes.append(output_size)

  # initilization of weights and biases for the model
  w, b = initialize_weights(node_sizes)

  return train_model(epoch, w, b, learning_rate, gamma,  batch_size, minibatch_size, free_time, clamped_time, delay, dt, node_sizes, n_activations, prediction_inp_size, pred=pred)

# run the code
train_accs, test_accs = run()