🖋️ PyTorch 的自动微分机制在进行PINNs代码实战前，有必要对 PyTorch 的自动微分机制，也就是autograd 模块的实现进行复习。它使得神经网络中的梯度计算变得非常简单和高效。以下是复习自动微分机制的基本概念和步骤：张量及梯度追踪 (requires_grad=True)在 PyTorch 中，张量（Tensor）是主要的数据结构，类似于多维数组。自动微分机制会记录对张量的操作以构建计算图，但前提是张量需要启用梯度追踪。这可以通过设置 requires_grad=True 来完成：x = torch.tensor(2.0, requires_grad=True) t = torch.tensor(3.0, requires_grad=True)这段代码创建了两个张量 x 和 t，并开启了梯度追踪。这样在计算过程中，PyTorch 会自动记录对 x 和 t 的操作，以便在之后计算梯度。定义函数并构建计算图PyTorch 会基于操作来构建计算图。每个操作都创建一个节点，并将输出作为另一个节点的输入。这样，计算图会记录下从输入到输出的所有依赖关系：f = x ** 2 + 3 * x * t + t ** 2这定义了一个函数：$$ f(x, t) = x^2 + 3xt + t^2 $$由于 x 和 t 都开启了梯度追踪，f 也会有梯度属性，并且 PyTorch 会自动构建计算图。计算梯度为了得到 f 对 x 和 t 的偏导数，我们使用 torch.autograd.grad 函数。它的常用参数包括：outputs: 需要求导的输出张量，这里是 f。inputs: 需要求导的变量，这里是 [x, t]。grad_outputs: 当 outputs 是标量时通常设为 1，表示链式法则中的外部导数。create_graph: 是否保留计算图。通常在二阶或更高阶导数计算时需要设为 True。Tips / 提示grad_outputs 参数在 PyTorch 中是用于指定目标张量 outputs 的“外部导数”（或梯度）：在微积分中，链式法则用于计算复合函数的导数。假设我们有两个变量 $y$ 和 $z$ ，且它们之间存在如下依赖关系：$$ z = f(y), \quad y = g(x) $$那么，如果我们要计算$\frac{\partial z}{\partial x}$，可以使用链式法则得到：$$ \frac{\partial z}{\partial x} = \frac{\partial z}{\partial y} \cdot \frac{\partial y}{\partial x} $$在 PyTorch 的 autograd 中，如果我们想要计算 z 对 x 的梯度，grad_outputs 就提供了“外部导数” $\frac{\partial z}{\partial y}$，从而帮助 autograd 计算梯度。这里设置设置grad_outputs=torch.ones_like(f)其实就是另“外部导数” $\frac{\partial z}{\partial y}=1$.gradients = torch.autograd.grad( outputs=f, inputs=[x, t], grad_outputs=torch.ones_like(f), create_graph=True )这行代码会返回 f 对 x 和 t 的偏导数，并将结果存储在 gradients 中。获取偏导数结果print("df/dx:", gradients[0]) # df/dx: tensor(13., grad_fn=<AddBackward0>) print("df/dt:", gradients[1]) # df/dt: tensor(12., grad_fn=<AddBackward0>)这将显示 f 对 x 和 t 的偏导数。根据上面 $f(x, t)$ 的定义，当 $x=2, t=3$ 时，可得：$$ \frac{\partial f}{\partial x} = 2x+3t=13, \frac{\partial f}{\partial t} = 3x+2t=12 $$矩阵、向量求导在 PyTorch 中，对矩阵或张量求导的过程可以理解为对每个元素进行逐元素（element-wise）求导。假设我们有两个张量 X 和 u，其中 u 是 X 的函数。X 的形状是 (2, 2)，即：# 假设 x 和 t 是这些值，直接设置 requires_grad=True x = torch.tensor([2.0, 3.0], requires_grad=True) t = torch.tensor([0.0, 1.0], requires_grad=True) from torch.autograd import Variable # 创建叶子节点副本 X = Variable(torch.stack([x, t], dim=1), requires_grad=True)u 的形状与 X 相同，因为 u 是对 X 的逐元素操作：u = X ** 2 + 2 * X对于矩阵 u 中的每个元素 u[i, j]，它是 X[i, j] 的函数。因此，计算 u 对 X 的梯度相当于逐元素对 X 中的每个元素计算偏导数。这会生成一个与 X 相同形状的矩阵，每个元素对应于 u 的某个元素对 X 相应元素的偏导数。torch.autograd.grad 可以对 u 中的每个元素求 X 中每个元素的偏导。由于 u 和 X 是逐元素关联的，torch.autograd.grad 会计算每个 u[i, j] 对应的 X[i, j] 的梯度。在代码中，grad_outputs=torch.ones_like(u) 用来设置“外部导数”，这相当于在链式法则中使用单位向量来表示标量梯度。具体来说，对于 u[i, j] 对 X[i, j] 的偏导数，我们设定 grad_outputs[i, j] = 1，即每个元素都作为独立的标量。du_dX = torch.autograd.grad( inputs=X, # 计算 X 的梯度 outputs=u, # u 是基于 X 的函数 grad_outputs=torch.ones_like(u), # 设置外部导数为 1 retain_graph=True, create_graph=True, allow_unused=True )[0]du_dX[:, 0] 和 du_dX[:, 1] 分别给出了 u 对 X 的两个列的偏导：du_dX[:, 0] 表示 u 对 x 的偏导数。du_dX[:, 1] 表示 u 对 t 的偏导数。print("du/dx:", du_dX[:, 0]) # du/dx: tensor([6., 8.], grad_fn=<SelectBackward0>) print("du/dt:", du_dX[:, 1]) # du/dt: tensor([2., 4.], grad_fn=<SelectBackward0>)🔛 求解过程IntroBurgers方程是一类非线性偏微分方程，可以用于模拟流体动力学、非线性声学和交通流等。在论文中，PINNs使用自动微分计算网络输出对输入的导数，并将其替代到Burgers方程中形成损失函数。通过最小化损失函数，使得神经网络的输出逼近真实解。RAISSI M, PERDIKARIS P, KARNIADAKIS G 2017. Physics Informed Deep Learning (Part I): Data-driven Solutions of Nonlinear Partial Differential Equations.Burgers方程的具体形式为：$$ u_t + u u_x - \left(\frac{0.01}{\pi}\right) u_{xx} = 0, x \in [-1, 1], \, t \in [0, 1] $$生成空间-时间网格h, k = .1, .1 x = torch.arange(-1, 1 + h, h) # len(x) = 21, 空间维度 t = torch.arange(0, 1 + k, k) # len(t) = 11, 时间维度x 和 t 分别定义了空间和时间维度的点，通过步长 h 和 k 划分区间。这里 x 的范围为 [-1, 1]，而 t 的范围为 [0, 1]。X = torch.stack(torch.meshgrid(x, t, indexing='ij')).reshape(2, -1).T使用 torch.meshgrid 生成了一个二维的网格，它包含了空间和时间上所有组合的点。然后用 torch.stack 将 x 和 t 的网格组合成 [len(x) * len(t), 2] 的张量 X。该张量的每一行表示一个 (x, t) 组合，这样的结构便于后续计算。定义边界条件边界条件和初始条件的定义是PINNs的关键步骤。对于Burgers方程，训练数据包括空间上的边界条件和时间上的初始条件：bc1 = torch.stack(torch.meshgrid(x[0], t, indexing='ij')).reshape(2, -1).T # 下边界点 bc2 = torch.stack(torch.meshgrid(x[-1], t, indexing='ij')).reshape(2, -1).T # 上边界点 ic = torch.stack(torch.meshgrid(x, t[0], indexing='ij')).reshape(2, -1).T # 起始时刻点bc1 和 bc2 分别代表空间维度上边界的 (x, t) 点，其中 bc1 对应 x=-1 处的下边界，bc2 对应 x=1 处的上边界。ic 表示初始条件的采样点，即 t=0 时刻整个空间维度的所有点 (x, t=0)。合并训练数据将 bc1、bc2 和 ic 合并成 X_train，它包含了所有边界和初始条件的点，用于训练PINN模型使其满足这些物理边界和初始条件。X_train = torch.cat([bc1, bc2, ic]) # X_train定义目标值（真实解）y_bc1 = torch.zeros(len(bc1)) y_bc2 = torch.zeros(len(bc2)) y_ic = -torch.sin(math.pi * ic[:, 0]) y_train = torch.cat([y_bc1, y_bc2, y_ic]).unsqueeze(1)y_bc1 和 y_bc2 分别是上下边界的目标值，这里均为0，即满足边界条件 $u(t, -1) = u(t, 1) = 0$。y_ic 是初始条件的目标值，由 $u(0, x) = -\sin(\pi x)$ 确定。将 y_bc1、y_bc2 和 y_ic 合并成 y_train，以与 X_train 中的点对应，最终构成神经网络的目标输出。数据迁移到GPUdevice = torch.device("cuda") if torch.cuda.is_available() else torch.device("cpu") X = X.to(device) X.requires_grad = True X_train = X_train.to(device) y_train = y_train.to(device)将数据迁移到GPU（如果可用），加速计算。设置 X.requires_grad = True 以启用自动微分，因为Burgers方程中涉及到对网络输出 $u(t, x)$ 关于时间和空间的偏导数。多层感知机首先建立一个深度多层感知机用于近似非线性函数$u(t,x)$：class MultiLayerPerceptron(nn.Module): """ Multi-layer Perceptron. """ def __init__( self, input_size, hidden_size, output_size, depth, act=torch.nn.Tanh, ): """ :param input_size: input size. :param hidden_size: hidden size. :param output_size: output size. :param depth: depth of the network. :param act: activation function. """ super(MultiLayerPerceptron, self).__init__() layers = [('Input', torch.nn.Linear(input_size, hidden_size)), ('Input_Activation', act())] for i in range(depth): layers.append( ('Hidden_%d' % i, torch.nn.Linear(hidden_size, hidden_size)) ) layers.append(('Activation_%d' % i, act())) layers.append(('Output', torch.nn.Linear(hidden_size, output_size))) self.layers = torch.nn.Sequential(OrderedDict(layers)) def forward(self, x): out = self.layers(x) return out该 MLP 的层数和激活函数类型都可以灵活调整，便于适应不同的应用需求。建立PINNs网络class PINNForBurgers: def __init__(self, device, X, X_train, y_train, model): self.model = model self.X = X self.X_train = X_train self.y_train = y_train self.criterion = torch.nn.MSELoss() self.iter = 1 self.optimizer = torch.optim.LBFGS( self.model.parameters(), lr=1.0, max_iter=50000, max_eval=50000, history_size=50, tolerance_grad=1e-7, tolerance_change=1.0 * np.finfo(float).eps, line_search_fn="strong_wolfe", # better numerical stability ) self.adam = torch.optim.Adam(self.model.parameters()) def loss_func(self): # this is more like a not so elegant hack to zero grad both optimizers self.adam.zero_grad() self.optimizer.zero_grad() y_pred = self.model(self.X_train) loss_data = self.criterion(y_pred, self.y_train) u = self.model(self.X) du_dX = torch.autograd.grad( inputs=self.X, # The shape of X: [len(x) * len(t), 2], X的第一列是空间坐标x, 第二列是是时间坐标t. outputs=u, grad_outputs=torch.ones_like(u), retain_graph=True, create_graph=True )[0] du_dt = du_dX[:, 1] du_dx = du_dX[:, 0] du_dxx = torch.autograd.grad( inputs=self.X, outputs=du_dX, grad_outputs=torch.ones_like(du_dX), retain_graph=True, create_graph=True )[0][:, 0] loss_pde = self.criterion(du_dt + u.squeeze() * du_dx, 0.01 / math.pi * du_dxx) loss = loss_pde + loss_data loss.backward() if self.iter % 100 == 0: print(self.iter, loss.item()) self.iter = self.iter + 1 return loss def train(self): self.model.train() for i in range(1000): self.adam.step(self.loss_func) self.optimizer.step(self.loss_func) def eval_(self): self.model.eval()这里面最主要的就是损失函数 (loss_func)的建立，这里采用的是KUMAR V, GLEYZER L, KAHANA A, et al. MyCrunchGPT: A chatGPT assisted framework for scientific machine learning[C]//,202310.48550/arXiv.2306.15551.提到的方法：也就是在损失函数中加入物理损失来计算求解。这里的损失函数定义为 loss = loss_pde + loss_data。其中，loss_data也就是已知边界点的计算误差：self.criterion(y_pred, self.y_train)；loss_pde则是指物理损失 self.criterion(du_dt + u.squeeze() * du_dx, 0.01 / math.pi * du_dxx)，也就是方程 $u_t + u u_x - \left(\frac{0.01}{\pi}\right) u_{xx} = 0$ 的误差。训练结果从最后的结果可以看出，我们根据已知边界点（x符号）的数据信息以及物理方程的定义拟合出了非线性函数 $u(t,x)$ 的近似表达式。提取不同时刻下对应的函数解：🪛 方程常数项的学习如上面所述，Burgers方程的具体形式为：$$ u_t + u u_x - \left(\frac{0.01}{\pi}\right) u_{xx} = 0, x \in [-1, 1], \, t \in [0, 1] $$但是很多工程问题我们是并不能确切地使用完整的物理方程来进行描述的。例如，如果我们不确定 $u u_x$ 与 $u_{xx}$ 前的系数，那么我们也可以借助于神经网络强大的学习能力来进行参数“拟合”。class PhysicsInformedNN: """ The physics-guided neural network. """ def __init__(self, X, u, model, device): # Data self.x, self.t = self.handle_m(X) self.u = torch.tensor(u).float().to(device) # PED params self.lambda_1, self.lambda_2 = [torch.nn.Parameter( torch.tensor( [v], requires_grad=True ).to(device) ) for v in [0.0, -6.0]] self.nn = model # 是 lambda_1 与 lambda_2 成为模型的一部分, 从而被优化器自动管理, 并在训练过程中被更新. for name, param in {'lambda_1': self.lambda_1, 'lambda_2': self.lambda_2}.items(): self.nn.register_parameter(name, param) # ...Warning / 注意在 Burgers 方程里，$u_{xx}$ 前面的常数项被称为扩散系数（diffusion coefficient），其控制了扩散或扩散速率的强弱，且永远是正值，这是因为如果扩散系数是负值，会导致“负扩散”现象，即不均匀性反而增加，这是不符合物理实际的。那我们应该如何限制自定义参数永远为正值呢？f = u_t + self.lambda_1 * u * u_x - torch.exp(self.lambda_2) * u_xx其实，只要使用torch.exp()来将其包裹起来，根据指数函数的特点，无论self.lambda_2在训练过程中如何变化，$u_xx$ 前的系数torch.exp(self.lambda_2)永远都为正值。具体的实现方式就是将 $u u_x$ 与 $u_{xx}$ 的系数设置为可学习参数，将其绑定注册到神经网络中，即可在神经网络的训练过程中进行参数的“拟合”。对比结果最后的结果对比证明，当原始数据不存在噪音的情况下，PINN基本可以“完美地”拟合 $u u_x$ 与 $u_{xx}$ 的常数项，其真实值分别为 1 和 $\frac{0.01}{\pi}$ ，拟合的结果为0.99952和0.0031969，当数据存在1%的噪音时，拟合精度显著下降。💡 一个小Trick在PINNs-Torch: Enhancing Speed and Usability of Physics-Informed Neural Networks with PyTorch提供的代码中，作者使用了个小Trick——采用两个优化器（LBFGS 和 Adam）对训练过程中的PINN进行优化。定义如下：self.optimizer = torch.optim.LBFGS( self.nn.parameters(), lr=1.0, max_iter=50000, max_eval=50000, history_size=50, tolerance_grad=1e-5, tolerance_change=1.0 * np.finfo(float).eps, line_search_fn="strong_wolfe" # can be "strong_wolfe" ) self.optimizer_Adam = torch.optim.Adam(self.nn.parameters())在训练函数中的具体调用方式为：def train(self, epoch_num: int): self.nn.train() for n_iter in range(epoch_num): # ... loss = ... # Backward and optimize self.optimizer_Adam.zero_grad() loss.backward() self.optimizer_Adam.step() # ... # Backward and optimize self.optimizer.step(self.loss_func)使用两个优化器（LBFGS 和 Adam）是设计PINN网络常见的优化策略，目的是结合两种优化算法的优点，达到快速收敛和高精度优化的效果。两种优化器的特点：Adam 优化器优点快速收敛：Adam 是一种自适应优化算法，能够根据梯度的历史记录动态调整学习率，在优化初期表现优异。鲁棒性：适合处理非凸问题和梯度震荡，尤其是参数较多或损失函数复杂时。全局探索：Adam 能够快速探索参数空间，找到一个接近最优解的区域。不足:在优化后期，Adam 的更新步长可能过大，难以细致地调整参数，导致优化精度不足。LBFGS 优化器优点:高精度：LBFGS 是一种拟牛顿方法，利用二阶信息（Hessian 矩阵的近似）进行优化，收敛速度快且优化结果精确。适合平滑问题：在损失函数较平滑时（例如 PINNs 中的物理残差），LBFGS 可以很好地找到局部最优解。不足:计算开销大：由于涉及近似二阶导数的计算，LBFGS 的开销较大。对初始点敏感：如果初始点较差，LBFGS 可能陷入次优解。使用两个优化器是为了取长补短，结合两者的优势。🧰 半监督神经网络物理信息神经网络本质上是一种半监督学习方法，因为它结合了监督学习和无监督学习的特性：监督学习部分：使用已知的边界条件和初始条件数据进行监督学习，self.x_u 和 self.t_u: 边界条件点的位置，self.u: 这些点上已知的精确解。无监督学习部分：PINNs 在定义域内随机采样伪点这些点没有已知的解，但它们需要满足偏微分方程的物理约束，计算Loss的过程中，通过让物理方程的误差趋近于0来进行梯度优化。例如上图中，我们可以利用边界点（已知对应 $u$ 值的点）来进行计算神经网络部分的Loss，然后在平面内随机采样一些点（不知道对应 $u$ 值，但这些点预测的 $u$ 要符合物理方程）作为物理信息部分的Loss。
文章头图来自：CAI S, MAO Z, WANG Z, et al. 2021. Physics-informed neural networks (PINNs) for fluid mechanics: a review. Acta Mechanica Sinica [J], 37: 1727-1738.🥕 基本信息1、什么是物理信息神经网络？Expand / 拓展物理信息神经网络 (Physics-Informed Neural Networks, PINNs) 是一种结合物理学知识和深度学习的强大工具，特别适用于求解偏微分方程（PDE）等复杂物理问题。PINNs 将神经网络的泛化能力与物理定律结合，通过编码物理方程中的约束条件来指导神经网络的学习。这种方法特别适用于那些难以通过传统数值方法求解的复杂问题，比如流体力学、热传导、结构力学等领域。简单来说，可以理解为下图：开题报告2、PINNs的优势（最明显的两点）：Expand / 拓展总而言之，言而总之，PINNs有以下优点：少量数据需求：由于物理方程本身已经提供了大量的约束信息，PINNs 通常需要的训练数据较少，在小样本情况下仍能表现出较好的效果。可解释性：因为 PINNs 是基于物理方程的，解的可解释性比纯数据驱动的神经网络模型更强。🤔 初步认识目前主要有两种方式将物理信息嵌入到神经网络中。一种是将偏微分方程（Partial Differential Equations，PDE）的先验知识作为损失函数隐式地嵌入到神经网络中，称为软约束方式。另一种是将PDE的先验知识显式地嵌入到神经网络架构中，称为硬约束方式。方式思路特点典型代表缺点软约束将神经网络视为方程解的逼近器不依赖于标签数据、要求 PDE 具体形式和方程系数是完全已知的、通常用来求解单个方程实例PINNs难以平衡不同的损失函数项，目前主要采用网格搜索方法来确定每个损失函数项对应的权重系数，但对于损失函数项较多的PDE 需要繁琐的调参工作，效率很低。硬约束将物理方程的信息显式地嵌入到神经网络架构中利用少量得标签数据即可训练好神经网络-需要针对特定问题进行神经网络整体架构的设计，移植性、通用性较差下面，以软约束方式为例，来对PINNs有个初步的认识。介绍下图来自文章：KUMAR V, GLEYZER L, KAHANA A, et al. MyCrunchGPT: A chatGPT assisted framework for scientific machine learning[C]//,202310.48550/arXiv.2306.15551.，在我学习PINNs网络的初期，这个图片让我比较清晰地理解了PINNs到底是怎么基于神经网络来进行实现的。这个图片由三部分组成：神经网络：$\mathbf{NN}\left( w, b \right) $；偏微分方程：$\mathbf{PDE}:\mathcal{L} \left( \boldsymbol{u}\left( x,\,\,t \right) , \theta \right) =g$；损失函数：$\mathbf{LOSS}_{MSE}={\color[RGB]{0, 128, 255} \mathrm{MSE}_{\left\{ \mathrm{u}, \mathrm{BC}, \mathrm{IC} \right\}}}+{\color[RGB]{0, 255, 64} \mathrm{MSE}_R}$。下面将分别介绍这三个组成部分。注意：对于下面的内容，您需要有深度学习方面的基础。神经网络图片的左上部分，你可以看到输入项 $ x $ 和 $ t $（代表空间和时间的坐标）被输入到神经网络 $\mathbf{NN}$ 中，经过多层神经元和激活函数的处理，网络输出的是 $ u_{\text{NN}} $，也就是神经网络对未知函数 $ u $ 的预测解。这个过程可以理解为神经网络尝试拟合未知物理量，比如某个区域的温度、压力、速度等随时间和空间变化的分布。偏微分方程图片右上部分展示了如何使用自动微分（Automatic Differentiation）来获取神经网络输出 $ u_{\text{NN}} $ 相对于输入 $ x $ 和 $ t $ 的偏导数（例如 $ \frac{\partial u}{\partial t} $ 和 $ \frac{\partial^2 u}{\partial x^2} $）。自动微分是一种通过神经网络计算导数的高效方法，可以直接在网络中对输出变量进行求导，而不需要通过数值差分的方式。有了这些导数之后，我们可以将它们代入偏微分方程（PDE）中，检查网络解 $ u_{\text{NN}} $ 是否符合物理定律。偏微分方程通常表示为 $\mathcal{L} \left( \boldsymbol{u}\left( x,\,\,t \right) , \theta \right) =g$，其中 $ \mathcal{L} $ 是作用在 $ u $ 上的物理算子，$ g $ 是已知的源项或外部驱动力。通过代入导数，我们计算出方程的残差 （Residual $ R $），即：$$ R = \mathcal{L}(\boldsymbol{u}, \theta) - g $$这个残差反映了网络输出解和实际物理方程解之间的差异。损失函数神经网络的损失不仅仅依赖于标准的监督学习误差 ${\color[RGB]{0, 128, 255} \mathrm{MSE}_{\left\{ \mathrm{u}, \mathrm{BC}, \mathrm{IC} \right\}}}$。通常情况下，监督学习的损失计算的是神经网络预测值与已知数据之间的均方误差（MSE）。在 PINNs 中，神经网络会根据初始条件（IC）和边界条件（BC）对 $ u_{\text{NN}} $ 进行监督训练。也就是说，网络需要保证在空间和时间边界以及初始时刻，输出解是准确的。除了标准的监督损失外，PINNs 还会将物理方程的残差 $ R $ ，也就是 ${\color[RGB]{0, 255, 64} \mathrm{MSE}_R}$ 纳入到损失函数中。通过最小化这个残差，网络被强制去学习一个解，不仅能拟合给定的数据，还能在整个定义域内尽量满足物理方程的约束。因此，PINNs 的总损失可以写作：$$ \mathbf{LOSS}_{MSE}={\color[RGB]{0, 128, 255} \mathrm{MSE}_{\left\{ \mathrm{u}, \mathrm{BC}, \mathrm{IC} \right\}}}+{\color[RGB]{0, 255, 64} \mathrm{MSE}_R} $$这里的 ${\color[RGB]{0, 255, 64} \mathrm{MSE}_R}$ 是物理方程残差的均方误差，用来确保神经网络的解满足物理定律。神经网络最终优化的目标不仅是数据上的误差，还包括物理方程残差。Warning / 注意虽然标签数据可以作为训练神经网络时的额外损失函数项，但是PINNs方法的训练过程可以完全不依赖于标签数据。因此，PINNS方法也可以视作是一种无监督学习方法。工作流程综上所述，PINNs的工作流程可以总结为：输入与神经网络预测：输入空间和时间的坐标 $ x $ 和 $ t $，神经网络会预测出对应时刻和空间位置的物理量 $ u_{\text{NN}} $。自动微分计算导数：通过自动微分，计算 $ u_{\text{NN}} $ 对时间和空间的导数。这些导数是用于后续代入物理方程的关键。计算残差：将计算出的导数代入偏微分方程，得到物理方程残差 $ R $，并将其作为物理损失的一部分。优化损失函数：总损失由两部分组成：监督学习的均方误差（用于拟合初始条件、边界条件等已知点的数据）和物理方程的残差损失。通过不断迭代优化，使得网络解既能拟合数据，又符合物理方程。🏗️ 网络架构设计注意点网络层：得益于物理信息的融合，简单的全连接神经网络PINNs即可表现出较高的精度，反而越复杂的网络结构有可能导致难以收敛；目前，如何确定全连接网络中最优的隐藏层数目以及每个隐藏层的神经元数目，尚无数学分析的方法，主要依赖手动调参。激活函数：由于 PINNs 需要计算高阶导数，因此选择高阶导数存在且连续的激活函数很重要。Tanh 激活函数具有良好的数学特性，比如它的高阶导数是连续的。因此，原始的 PINNs 方法以及大多数研究者倾向于选择 Tanh 作为激活函数，以确保在计算高阶导数时模型的稳定性和精度。除此之外，也有相关研究表明 Sine 激活函数与自适应激活函数也可以显著地加快训练速度并提高精度。损失函数自适应加权：GradNorm；联合似然公式；学习率退火算法。模型权重初始化：随机初始化；如果目标任务与源任务相似度较高，也可以采用迁移学习的思路——将源任务中训练得到的模型权重作为目标任务的初始权重；元学习方法——通过多次相关任务的学习来改进学习算法本身，从而加快处理新任务的速度。📝 基本概念偏微分方程偏微分方程（Partial Differential Equation, PDE）：描述未知函数及其偏导数之间关系的方程，其解析解难以显式地表达出来，通过数值求解方法（数值求解的核心是将连续的数学问题离散化，将一个复杂的方程分解成许多简单的、可以通过迭代和逼近解决的小问题。）得到PDE的近似解是求解方程的首选。有限差分方法有限差分方法（Finite Difference Method，FDM）是一种传统的数值求解方法。它的基本思想是通过泰勒展开将控制方程的导数替换为离散网格节点处函数值的差商，然后通过求解差分方程近似求解原微分方程。举个例子🌰，设函数 $ y(x) = x^2 $，我们想在 $ x = 1 $ 处求其导数，步长 $ h = 0.01 $。真实导数为：$$ \frac{dy}{dx} = 2x $$在 $ x = 1 $ 处，真实导数为 $ 2 $。使用前向差分法：$$ \frac{dy}{dx} \approx \frac{y(x+h) - y(x)}{h} = \frac{y(1+0.01) - y(1)}{0.01} = \frac{(1.01)^2 - 1^2}{0.01} = \frac{1.0201 - 1}{0.01} = 2.01 $$这个结果与真实导数 $ 2 $ 很接近。使用中心差分法：$$ \frac{dy}{dx} \approx \frac{y(x+h) - y(x-h)}{2h}= \frac{y(1+0.01) - y(1-0.01)}{2 \times 0.01} = \frac{(1.01)^2 - (0.99)^2}{0.02} = \frac{1.0201 - 0.9801}{0.02} = 2 $$这种方法通常更精确，因为它对 $x$ 前后点的变化都进行了平均，中心差分法在这个例子中与真实导数完全相符。Data-driven solution/discovery of PDE这两个概念分别对应于使用数据来解决偏微分方程（PDE）以及从数据中自动推导出描述物理系统的偏微分方程。简单来说可以这样理解：基于数据的偏微分方程求解（Data-driven solution of PDE）：利用数据训练出一个“函数逼近器”，可以输入变量值，得到已知偏微分方程的解。基于数据的偏微分方程发现（Data-driven discovery of PDE）：根据数据推导出这个“函数逼近器”本身，也就是找到能够描述数据规律的偏微分方程。参数化和非线性偏微分方程的通用形式可以表示为：$$ u_t + \mathcal{N}[u; \lambda] = 0 $$其中：$ u_t $ 表示未知函数 $ u $ 对时间 $ t $ 的偏导数，即随时间变化的速率。$ \mathcal{N}[u; \lambda] $ 是一个非线性算子，作用在未知函数 $ u $ 上，可能包括 $ u $ 关于空间或其他变量的偏导数。$ \lambda $ 表示一组参数，控制着方程的特性或行为。参数可以是系统中的物理常数（如扩散系数、波速）或者环境条件。这种形式的偏微分方程通常用于描述随时间演化的物理过程，比如热传导、波动或流体流动等。这里的非线性算子 $ \mathcal{N}[u; \lambda] $ 可以有不同的具体形式，取决于具体应用，例如可能包含一阶或二阶导数，甚至更复杂的非线性项。Burgers方程是一类著名的非线性偏微分方程，常用于描述流体力学、交通流、湍流等物理现象。Burgers方程在一维空间中的标准形式为：$$ \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} = \nu \frac{\partial^2 u}{\partial x^2} $$其中：$ u = u(x, t) $ 表示一个随空间 $ x $ 和时间 $ t $ 变化的未知函数。$ \nu $ 是粘性系数（或扩散系数），表示粘性或扩散效应的强度。$ u \frac{\partial u}{\partial x} $ 是一个非线性项，描述对流（convective）现象。$ \nu \frac{\partial^2 u}{\partial x^2} $ 是一个扩散项，描述粘性或扩散效应。Burgers方程可以看作参数化和非线性偏微分方程通用形式的一个特例：对应到Burgers方程中，我们可以将上面表达式中的算子 $ \mathcal{N}[u; \lambda] $ 具体化为：$$ \mathcal{N}[u; \lambda] = -u \frac{\partial u}{\partial x} + \nu \frac{\partial^2 u}{\partial x^2} $$其中参数 $ \lambda $ 就是粘性系数 $ \nu $。这样，Burgers方程的形式可以表示为：$$ u_t + \left(-u \frac{\partial u}{\partial x} + \nu \frac{\partial^2 u}{\partial x^2}\right) = 0 $$这表明，Burgers方程是一种特定形式的带参数的非线性偏微分方程，其中的非线性项 $ u \frac{\partial u}{\partial x} $ 和扩散项 $ \nu \frac{\partial^2 u}{\partial x^2} $ 可以分别描述不同的物理效应。这种形式使得Burgers方程既有非线性行为（如冲击波）又有扩散行为（如平滑效应），因此在物理和工程中有广泛应用。🏙️ PINNs在城市水系统中的应用不同知识整合的层级：当现有知识较为全面时，可以通过数据驱动模型学习已有的数学方程描述的规律，从而成为完整机械模型的高效替代。当现有知识不完整时，可以将这些知识部分嵌入机械模型，并与数据驱动模型结合，形成混合模型（hybrid model）。如果知识无法用数学方程完整描述，也可以用于指导数据驱动模型的结构设计，减少需要训练的参数数量，从而提高模型性能。PINN的基本架构核心原理：PINN将物理定律（如偏微分方程或常微分方程）直接嵌入到神经网络中。输入/输出：网络接收输入（如时间 $t$、空间位置 $x$）并生成输出（如速度 $u, v$）。物理约束：通过自动微分计算输出相对于输入的导数，并将其与物理定律（如连续性方程、动量方程）进行对比。损失函数：总损失由三部分组成：物理损失（Physical Loss）：约束模型输出符合物理方程。边界条件损失（BC Loss）：确保边界条件满足。初始条件损失（IC Loss）：约束初始条件。训练方法：通过反向传播优化损失函数，从而得到物理一致的网络参数。混合建模结构串行结构（Serial Structure）：机械模型用于处理已知部分。数据驱动模型用于补充未知部分（黑箱）。适用场景：系统的部分过程物理规律清晰。并行结构（Parallel Structure）：机械模型和数据驱动模型并行运行，数据驱动模型负责补偿机械模型与实际观测之间的偏差。适用场景：机械模型存在明显预测误差，但仍能提供部分有价值的参考。知识驱动的甲烷排放建模背景：甲烷是污水管道的重要温室气体，其排放受水面角度、管道坡度和流量等影响。建模思路：利用物理几何知识（如曼宁公式）引导模型结构。数据驱动部分负责拟合参数 $k_r, k_\theta, \alpha, \beta, \gamma$，优化甲烷排放预测。知识驱动的土壤温度预测背景：埋地污水管道的热传导影响土壤温度，涉及复杂的物理过程。建模思路：物理组件：采用一维热传导方程的解析解作为基础。数据驱动组件：利用与曲率相关的变量 $\frac{\delta}{\delta + R}$ 描述复杂影响。模型优势：减少了参数数量，提高了训练效率和模型解释性。Summary of the pros and cons of mechanistic models, data-driven models and knowledge informed data-driven models.Modelling approachesMechanistic ModelsData-Driven ModelsKnowledge-Informed Data-Driven modelsStrength1. Based on established physical, chemical, and biological principles;2. Widely accepted and applied in the water industry; 3. Excellent interpretability - well-suited for understanding system dynamics.1. Leverages abundant data from advanced sensor technologies;2. Can capture system complexities and reduce computational costs;3. Self-adaptive to changing environmental and operational conditions.1. Combines the strengths of mechanistic and data-driven models;2. Reduces data requirements by incorporating mechanistic knowledge;3. Enhances model interpretability and efficiency.Limitations1. Requires comprehensive mechanistic knowledge;2. May be overly complex and/or computationally demanding for real-time applications;3. Modelling practice can be highly complicated and cannot capture time-varying processes.1. Data-hungry, requiring comprehensive datasets for training;2. Limited interpretability;3. Difficulties in acquiring high-quality datasets1. Complexity in model integration;2. Requires expertise in both mechanistic and data-driven modelling;3. May still face challenges in real-time applications.🗂️ 参考资料视频Steve Brunton, Physics Informed Neural Networks (PINNs) Physics Informed Machine Learning(YouTube).文献PINNs开山之作：RAISSI M, PERDIKARIS P, KARNIADAKIS G E 2019. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics [J], 378: 686-707.查文舒, 李道伦, 沈路航, 等. 基于神经网络的偏微分方程求解方法研究综述[J]. 力学学报, 2022, 54(03): 543-556.李野, 陈松灿. 基于物理信息的神经网络：最新进展与展望[J]. 计算机科学, 2022, 49(04): 254-262.黄翔. 基于物理信息嵌入神经网络的偏微分方程求解方法研究[D]. 中国科学技术大学, 2023.DUAN H, LI J, YUAN Z 2024. Making waves: Knowledge and data fusion in urban water modelling. Water Research X [J], 24: 100234.ASHRAF I, STROTHERM J, HERMES L, et al. 2024. Physics-Informed Graph Neural Networks for Water Distribution Systems. Proceedings of the AAAI Conference on Artificial Intelligence [J], 38: 21905-21913.代码https://github.com/rezaakb/pinns-torchhttps://github.com/jayroxis/PINNs
❓ 需求分析首先展示下我的问题，以便分析本篇博客中的哪些步骤适用于你的需要。已知条件：有某地区管网的.inp文件，并且有该地区水源节点历史两个月的余氯时间序列数据。求解结果：得到该管网结构中所有节点在历史两个月中的余氯值，更准确的结果如下表展示所示：Time节点 01节点 02……节点 m时刻 012.564.28……2.56时刻 023.564.56……2.56…………………………时刻 n3.443.89……2.56🫡 基础知识什么是WNTR？简单来说：WNTR（Water Network Tool for Resilience）是由美国国家标准与技术研究院（NIST）开发的一套用于分析水分配系统的开源工具，用于模拟和评估给水管网系统的性能和恢复力。主要作用包括：模拟水分配系统、恢复力分析、情景分析，以及优化调度等。简单来说，WNTR就是基于Python的高级版EPANET。EPANET中的Patterns是干什么的？在EPANET中，"Patterns"（模式）是用来表示随时间变化的水需求、压力或其他参数的变化模式。模式是一个数组或一系列因子，这些因子用于调整基础需求或其他输入参数，以模拟系统在不同时间段的行为。应用场景：水需求模式：模式通常用于表示水需求的日变化。例如，居民区的用水需求在早晨和傍晚会有高峰，而在夜间会较低。通过设置不同时间段的需求因子，可以精确模拟真实的用水情况。泵操作模式：模式还可以用于控制泵的操作。例如，可以设置泵在特定时间段的运行速度或开启/关闭状态，以模拟不同的运行策略。压力调节模式：在某些情况下，模式也可以用于调节系统中的压力，模拟压力管理策略。举个例子：假设有一个居民区的用水模式，表示一天24小时内的需求变化：时间步长（小时）: 1 模式因子: [0.5, 0.6, 0.7, 0.8, 1.0, 1.2, 1.5, 1.8, 2.0, 1.8, 1.6, 1.4, 1.2, 1.0, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.3, 0.4, 0.5]在这个模式中，0点到1点的需求因子是0.5，表示这段时间的需求是基础需求的50%。早晨6点到9点的需求因子是1.2到2.0，表示用水需求较高，为基础需求的120%到200%。🕹️ 具体步骤简要介绍在开始细致讲解前，先对整体的操作步骤进行一个大致总结，以便各位观众按需研究：以文本形式打开.inp文件，将管网水源节点的余氯时间序列添加为一个Pattern；在EPANET中打开这个.inp文件，将该Pattern与水源节点相关联，并进行相关设置；使用WNTR完成水力模拟并保存所有节点的水质数据。管网数据文件格式的调整在.inp文件中，可以将管网数据以文本的形式进行编辑，非常方便。那么第一步需要进行的操作就是：如果你的管网文件是.net格式，需要将其调整为.inp格式。这里有个注意点，将.net文件另存为.inp文件不能直接使用 File → Save As... 命令来操作，否则生成的文件会显示打开错误：正确的操作步骤为：File → Export → Network，这样导出的.inp文件才是完整并且可直接进行文本编辑的。将时间序列设置为Pattern完成上面的.inp文件格式转换后，就可以直接使用文本编辑工具（Notepad、VS Code等）打开并进行编辑了。在文件中直接搜索[PATTERNS]，定位到Pattern的存储位置，会有如下的字样：[PATTERNS] ;ID Multipliers ; 1 0.781 0.752 0.738 0.706 0.704 0.73 然后即可在后面进行Pattern的添加，例如：[PATTERNS] ;ID Multipliers ; 1 0.781 0.752 0.738 0.706 0.704 0.73 ; demo 1 0.986378898 0.974703668 0.976000916 0.980541283 0.97340642上面添加的Pattern名为demo，可以使用EPANET打开该文件进行查看：这里有个小问题，如果我的水源节点历史时间序列长度为2个月（60天）、每5分钟一个数据，那么一共会有60*24*(60/5)=17280条数据，也就是说我设置的Pattern长度也必须为17280。同时，还需要进行对应的时间设置：另外还需要将水源节点的余氯值与这个Pattern相关联：关于节点的两个属性：Initial Quality代表了水力模型起始运行那一刻的水质；Source Quality代表了水力模拟过程中，水源水质的变化情况。为什么要在EPANET中进行这些设置呢？就是因为我是用WNTR倒腾了几周还是会有各种各样的问题，一直没有解决，然后才采用了这种方法。在EPANET中设置的属性是会保存在.inp文件中的，但水力模拟的结果或者过程信息却是不会保存在文件中的。在WNTR中进行水力模拟其实，按照上面的步骤进行操作后，我们已经可以在EPANET中进行水力模拟，得到想要的结果了。但是问题就在于：在 Report → Table 的导出选项中，每次仅可以选择一个节点或者管道的属性进行导出：这如果一个网络拓扑结果中有成千上万个节点，那我得导出到猴年马月呢~😟因此，在这里我们借助于WNTR进行水力模拟，其可以得到Numpy数据格式的文件，非常方便进行数据分析与处理。使用WNTR进行水力模拟就非常简单了，直接上代码：# 导入WNTR包 import wntr # 加载EPANET模型 inp_file = r"C:\Users\myxc\Desktop\Anytown-2.inp" wn = wntr.network.WaterNetworkModel(inp_file) # 运行水力模拟(使用EPANET内核) sim = wntr.sim.EpanetSimulator(wn) results = sim.run_sim() # 获取所有节点的余氯数据 cl_concentration = results.node['quality'] * 1000 cl_concentration请注意：在WNTR中设置Pattern或者水源节点的一些具体属性会非常麻烦，但是应该也是可以进行的。本文作者折腾了很长时间，还是不能完美实现自己的需求。因此，才采用了在EPANET中进行属性设置、在WNTR中进行水力模拟的方法。其中，cl_concentration返回结果便是管网中所有节点在任意时刻的余氯监测值（模拟）：但是目前还有个疑问我没有搞明白：可以看到我上面对返回的数据矩阵进行了 * 1000操作，WNTR返回的数据单位需要 * 1000才与EPANET的相等。😕暂时还没找到原因~
介绍图（Graph）是一种表示对象之间关系的抽象数据结构。图由节点（Vertex）和边（Edge）组成，节点表示对象，边表示对象之间的关系。图可以用于建模各种实际问题，如社交网络、交通网络、电力网络等。NetworkX是一个用Python编写的库，专门用于创建、操作和研究复杂网络的结构、动态和功能。它提供了简单易用的接口来处理图论和网络结构。NetworkX适用于处理大型网络结构，并提供了许多内置的图算法，如路径寻找、图的构建和修改、节点属性操作等。NetworkX官方文档（网站）：https://networkx.org/；使用pip安装：pip install networkx并回车。基本概念无向图（Undirected Graph）import networkx as nx # 创建一个无向图 G = nx.Graph() # 添加节点 G.add_node(1) G.add_nodes_from([2, 3]) # 添加边 G.add_edge(1, 2) G.add_edges_from([(2, 3), (1, 3)]) # 查看图的节点和边 print("图的节点: ", G.nodes(), "; 图的边: ", G.edges(), '.') # 可视化 nx.draw(G, node_size=500, with_labels=True)控制台输出结果 - 无向图有向图（Directed Graph）有向图的创建方式很简单，只需要把上面无向图的对象：# 创建一个无向图 G = nx.Graph()换成：# 创建一个有向图 G = nx.DiGraph()即可。控制台输出结果 - 有向图有权图（Directed Graph）创建有权图时需要添加权重信息，且可视化的代码略有不同：import networkx as nx import matplotlib.pyplot as plt # 创建一个带权重的无向图 G = nx.Graph() # 添加带权重的边 G.add_edges_from([ (2, 3, {'diameter': 1.0,'length': 12.0}), (1, 3, {'diameter': 2.0,'length': 10.0}), (2, 1, {'diameter': 3.0,'length': 8.0}) ]) # 提取权重信息 edge_weights = nx.get_edge_attributes(G, 'length') # 可视化图 pos = nx.spring_layout(G) # 选择布局算法，这里使用弹簧布局 nx.draw(G, pos, with_labels=True, font_weight='bold', node_size=700, node_color='skyblue', font_color='black') # 绘制权重标签 nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_weights) # 显示图形 plt.show() # 查看图的节点和边 print("图的节点: ", G.nodes(), "; 图的边: ", G.edges(data=True), '.')控制台输出结果 - 有权图邻接矩阵邻接矩阵（Adjacency Matrix）： 邻接矩阵是一个二维矩阵，其中的行和列分别对应图中的节点。矩阵的元素表示节点之间是否存在边。对于无向图，邻接矩阵是对称的；对于有向图，邻接矩阵不一定对称。# 获取邻接矩阵(默认是稀疏矩阵格式) adj_matrix = nx.adjacency_matrix(G) # 将稀疏矩阵转换为密集矩阵(如果需要) dense_adj_matrix = adj_matrix.todense()请注意，如果你的图是有向图，你可以使用 nx.adjacency_matrix(G, directed=True) 来获取有向图的邻接矩阵。如果你想要自定义矩阵的表示方式，你可以使用 toarray() 方法将稀疏矩阵转换为 NumPy 数组。关联矩阵在 networkx 库中，nx.incidence_matrix(G) 用于生成图 G 的关联矩阵（Incidence Matrix）。关联矩阵和邻接矩阵不同，它是一种表示图中节点与边之间关系的矩阵。对于一个包含 $ N $ 个节点和 $ M $ 条边的图，其关联矩阵是一个 $ N \times M $ 的矩阵，其中矩阵的行对应图中的节点，列对应图中的边。关联矩阵中的元素定义如下：如果边 $ e $ 与节点 $ v $ 相连，则矩阵中对应位置的值为 1（或者对于有向图，出边为 -1，入边为 1）；如果边 $ e $ 与节点 $ v $ 不相连，则对应位置的值为 0。举个例子：# 创建一个图 G = nx.Graph() G.add_edges_from([(0, 1), (0, 2), (1, 2)]) # 获取关联矩阵 inc_matrix = nx.incidence_matrix(G) # 绘制节点和边 pos = nx.spring_layout(G) # 生成节点的布局 nx.draw(G, pos, node_size=500, with_labels=True) # 为每条边创建一个标签（例如，使用边的索引作为标签） edge_labels = {edge: i for i, edge in enumerate(G.edges())} # 将边的标签添加到图中 nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_labels, font_size=20)获取关联矩阵：inc_matrix.toarray()度、平均度、度分布、度矩阵以及常见复杂网络常见度分布下面的演示均以：G.add_edges_from([ # Fig. 7 (0, 1), (0, 2), (0, 7), (1, 0), (1, 2), (1, 3), (1, 4), (1, 7), (2, 0), (2, 1), (2, 3), (3, 1), (3, 2), (3, 4), (3, 5), (3, 7), (4, 1), (4, 3), (4, 5), (4, 6), (4, 7), (5, 3), (5, 4), (5, 6), (6, 4), (6, 5), (7, 0), (7, 1), (7, 3), (7, 4) ])演示图结构为示例。度度（Degree）的定义：对于无向图 G，节点 i 的度 $ \text{degree}(i) $ 是与节点 i 相连的边的数量。对于有向图 G，节点 i 的入度 $ \text{in-degree}(i) $ 是指向节点 i 的边的数量，出度 $ \text{out-degree}(i) $ 是从节点 i 出发的边的数量。NetworkX求度：# 获取节点的度 degrees = dict(G.degree()) print("节点的度：", degrees)节点的度： {0: 3, 1: 5, 2: 3, 7: 4, 3: 5, 4: 5, 5: 3, 6: 2}平均度平均度（Average Degree）是图中所有节点的度的平均值。对于无向图 G，平均度 $ \langle k \rangle $ 可以通过所有节点的度之和除以节点数得到。对于有向图 G，同样可以计算平均入度和平均出度。# 计算平均度 average_degree = sum(dict(G.degree()).values()) / len(G) print("平均度：", average_degree)平均度： 3.75度分布度分布（Degree Distribution）描述了图中节点的度的分布情况，即度为 k 的节点有多少个。度分布是图结构的一个重要特征，它可以帮助我们了解网络中节点的连接模式。# 获取度分布 degree_distribution = nx.degree_histogram(G) # 绘制度分布直方图 plt.bar(range(len(degree_distribution)), degree_distribution, tick_label=range(len(degree_distribution))) plt.xlabel("Degree") plt.ylabel("Frequency") plt.title("Degree Distribution") plt.show()度矩阵度矩阵（Degree Matrix）是一个表示图中节点度信息的对角矩阵。对于一个无向图，度矩阵的定义如下：对于无向图 G，其度矩阵 $ D $ 是一个 $ n \times n $ 的矩阵，其中 $ n $ 是图中的节点数。$ D $ 的对角线元素 $ D_{ii} $ 表示节点 $ i $ 的度，即与节点 $ i $ 相连的边的数量。# 获取度矩阵 np.diag(list(dict(G.degree()).values()))控制台输出结果 - 度矩阵复杂网络常见度分布在复杂网络理论中，度分布描述了网络中各节点的连接数（即度）的分布情况。不同类型的网络具有不同的度分布特性。典型复杂网络的度分布有：泊松分布、幂律分布（Power-Law Distribution）、指数分布等。泊松分布大多数节点的度数集中在较小的数值范围内，这是典型的泊松分布特征，符合Erdős-Rényi随机图模型的特点。import networkx as nx import matplotlib.pyplot as plt # 参数设置 n = 500 # 节点数 p = 0.05 # 每对节点间存在边的概率 # 使用 Erdős-Rényi 模型创建随机图 G = nx.erdos_renyi_graph(n, p) # 准备绘制网络结构图和度分布统计图 fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 6)) # 在左边绘制网络图 nx.draw(G, ax=ax1, node_size=30, with_labels=False) ax1.set_title("Erdős-Rényi Network (n={}, p={})".format(n, p)) # 在右边绘制度分布统计图 degrees = [G.degree(n) for n in G.nodes()] ax2.hist(degrees, bins=range(min(degrees), max(degrees) + 1), density=True) ax2.set_title("Degree Distribution") ax2.set_xlabel("Degree") ax2.set_ylabel("Frequency") plt.show()幂律分布Barabási-Albert 模型生成的网络展现了显著的幂律分布特征，即大多数节点拥有较少的连接，而少数节点（枢纽节点）拥有非常多的连接。m = 5 n = 200 # 节点数 G = nx.barabasi_albert_graph(n, m)指数分布大多数节点具有较低的度数，而只有少数节点具有较高的度数。这反映了指数分布的特点，即概率随着度数的增加而迅速减少。由于指数分布的特性，大多数节点集中在较低的度数区间，而高度数的节点数量迅速减少。import numpy as np import networkx as nx import matplotlib.pyplot as plt # 参数设置 m = 5 n = 1000 # 节点数 # 由于NetworkX没有直接生成指数分布网络的函数，我们可以手动创建一个具有指数分布度的网络 # 创建一个空的图 G = nx.Graph() # 添加节点 G.add_nodes_from(range(n)) # 生成一个符合指数分布的度序列 degree_sequence = np.random.exponential(scale=m, size=n).astype(int) # 为了确保总度数为偶数（每条边两个端点），如果是奇数则调整 if sum(degree_sequence) % 2 != 0: degree_sequence[0] += 1 # 生成具有指定度序列的随机图（可能会有多重边和自环） G = nx.configuration_model(degree_sequence) # 转换为简单图（无多重边和自环） G = nx.Graph(G) G.remove_edges_from(nx.selfloop_edges(G)) # 准备绘制网络结构图和度分布统计图 fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 6)) # 在左边绘制网络图 nx.draw(G, ax=ax1, node_size=30, with_labels=False) ax1.set_title("Exponential Distribution Network (n={}, scale={})".format(n, m)) # 在右边绘制度分布统计图 degrees = [d for n, d in G.degree()] ax2.hist(degrees, bins=range(min(degrees), max(degrees) + 1), density=True) ax2.set_title("Degree Distribution") ax2.set_xlabel("Degree") ax2.set_ylabel("Frequency") plt.show()Quote / 参考幂律分布与指数分布看起来没有什么差别啊？幂律分布：长尾特性：幂律分布的一个显著特点是它的长尾特性。这意味着即使是非常大的值也有可能出现，尽管其频率较低。对数坐标下的直线：在对数-对数坐标图中，幂律分布呈现为一条直线。现实世界的例子：社交网络中的连接分布、互联网链接的分布、城市的人口分布等。指数分布：迅速衰减：指数分布的特点是随着值的增加，其概率迅速减少。这意味着大值出现的可能性比幂律分布要小得多。半对数坐标下的直线：在半对数坐标图中（即纵坐标为对数坐标），指数分布表现为一条直线。现实世界的例子：无线电活动的持续时间、顾客在商店的停留时间等。平均近邻度平均近邻度（Average Nearest Neighbor Degree, Average Neighbor Degree）是指网络中一个节点的邻居节点的平均度数。这个指标可以用来判断网络中的节点是否倾向于与相似或不同度数的其他节点相连。平均近邻度的计算方法是：对于网络中的每个节点，计算其所有邻居节点的度数。对这些度数求平均值。重复此过程，计算所有节点的平均近邻度。公式为：$$ k_{nn,i} = \frac{1}{k_i} \sum_{j \in \mathcal{N}(i)} k_j $$其中，$ \mathcal{N}(i) $ 是节点 $ i $ 的邻居集合，$ k_i $ 是节点 $ i $ 的度数，$ k_j $ 是 $ i $ 的邻居节点 $ j $ 的度数。相关性(匹配)系数相关性(匹配)系数（Assortativity Coefficient, Degree Correlation Coefficient）衡量的是网络中高度数的节点倾向于与其他高度数的节点相连（正相关），还是与低度数的节点相连（负相关）。相关性系数的计算方法是：计算网络中每对连接节点的度数乘积。对这些乘积求平均值。计算度数的平均值的平方和度数平方的平均值。将这些值代入皮尔逊相关系数的公式中。公式为：$$ r = \frac{M^{-1} \sum_{i,j} A_{ij} k_i k_j - [M^{-1} \sum_{i} k_i]^2}{M^{-1} \sum_{i} k_i^2 - [M^{-1} \sum_{i} k_i]^2} $$其中，$ A_{ij} $ 是邻接矩阵中的元素，$ k_i $ 和 $ k_j $ 是节点的度数，$ M $ 是边的总数。下面举个例子对比分析下平均近邻度与相关性(匹配)系数：创建两个网络结构：# 创建两个不同的网络结构 # 网络1：经典的Erdős-Rényi随机图 G1 = nx.erdos_renyi_graph(n=100, p=0.05) # 网络2：Barabási-Albert无标度网络 G2 = nx.barabasi_albert_graph(n=100, m=2)计算两个网络的相关性（匹配）系数：G1_r = nx.degree_pearson_correlation_coefficient(G1) G2_r = nx.degree_pearson_correlation_coefficient(G2)计算两个网络的平均近邻度：average_neighbor_degree_G1 = nx.average_neighbor_degree(G1) average_neighbor_degree_G2 = nx.average_neighbor_degree(G2) # 计算两个网络的平均近邻度的平均值 G1_and =sum(average_neighbor_degree_G1.values()) / len(average_neighbor_degree_G1) G2_and = sum(average_neighbor_degree_G2.values()) / len(average_neighbor_degree_G2)绘图可视化观察：# 在左边绘制网络图 nx.draw(G1, ax=ax1, node_size=50, with_labels=False) nx.draw(G2, ax=ax2, node_size=50, with_labels=False) ax1.set_title(f"Graph A, r={round(G1_r, 3)}, and={round(G1_and, 3)}") ax2.set_title(f"Graph B, r={round(G2_r, 3)}, and={round(G2_and, 3)}") plt.show()相关性（匹配）系数：两个网络都显示出负相关性，但Barabási-Albert无标度网络（Graph B）的负相关性更强。这表明在无标度网络中，网络中的枢纽（高度数的节点）更可能与度数较低的节点相连，这是典型的无标度网络特性。平均近邻度：Barabási-Albert无标度网络的平均近邻度高于Erdős-Rényi随机图。这反映了无标度网络中节点连接的非均匀性，其中一些节点（即枢纽）与大量其他节点相连，从而提高了平均近邻度。总的来说，这两个网络展现出不同的结构特征，无标度网络的节点连通性更加集中和不均匀，而随机图则表现出更均匀的连接模式。这些特性在网络的稳定性和动态行为上有重要影响，例如在信息传播或疾病传播的研究中。拉普拉斯矩阵拉普拉斯矩阵（Laplacian Matrix）有多种定义方式，其中最常见的计算方法是使用度矩阵减去邻接矩阵。假设无向图 G 有 n 个节点，其邻接矩阵为 A，度矩阵为 D。标准拉普拉斯矩阵 L 的计算如下：计算度矩阵 D 的对角线元素，即每个节点的度：$ D_{ii} = \sum_{j} A_{ij} $；计算拉普拉斯矩阵 L：$ L = D - A $。# 获取邻接矩阵和度矩阵 adj_matrix = nx.adjacency_matrix(G).toarray() degree_matrix = np.diag(list(dict(G.degree()).values())) # 计算标准拉普拉斯矩阵 laplacian_matrix = degree_matrix - adj_matrix其他几种邻接矩阵：nx.laplacian_matrix(G)：返回G的拉普拉斯矩阵；nx.directed_laplacian_matrix(G)：返回有向图的拉普拉斯矩阵；nx.normalized_laplacian_matrix(G)：返回G的标准化拉普拉斯矩阵。路径和距离在图论中，路径和距离是描述图中节点之间连接关系和位置关系的重要概念。路径（Path）：在图中，路径是指图中的一系列节点，其中任意相邻两个节点之间都有边相连。路径的长度是指路径上边的数量。如果路径中的所有节点都是不同的，则路径是简单路径。距离（Distance）：在图中，两个节点之间的距离是指连接这两个节点的最短路径的长度。如果两个节点之间没有路径相连，则它们之间的距离通常被定义为无穷大。获取所有节点间的最短路径和距离获取图中的所有最短路径和距离：# 获取所有节点对之间的最短路径和距离 all_shortest_paths = dict(nx.all_pairs_shortest_path(G)) all_shortest_distances = dict(nx.all_pairs_shortest_path_length(G)) # 打印最短路径和距离 all_info = [] for source in all_shortest_paths: for target in all_shortest_paths[source]: paths = all_shortest_paths[source][target] distance = all_shortest_distances[source][target] print(f"最短路径从节点 {source} 到节点 {target}: {paths}, 距离: {distance}")节点间最短路径和距离获取节点间最短路径和距离：i, j = 1, 3 print(nx.shortest_path(G, i, j), nx.shortest_path_length(G, i, j))获取节点间的所有路径# 寻找所有简单路径 paths = list(nx.all_simple_paths(G, 6, 0)) print("所有简单路径:") for path in paths: print(path) # 独立路径的数量 print("独立路径的数量:", len(paths))这个代码计算了从节点 6 到节点 0 的所有简单路径。请注意，这些路径可能在某些较复杂的图结构中共享边。如果你的网络足够大或复杂，计算所有独立路径的数量可能会变得非常困难，甚至不可行，因为可能的路径数量随着网络的大小和复杂性指数级增长。离心率离心率（Eccentricity）用于描述一个节点到图中所有其他节点的最短路径中最长的那一个。对于图中的节点 v，其离心率定义为从 v 到图中所有其他节点的最短路径的最大长度。公式表示为：$$ \text{Eccentricity}(v) = \max_{u \in G} d(v, u) $$其中，d(v, u) 是节点 v 和 u 之间的最短路径长度。如果不指定节点 v，nx.eccentricity(G) 返回一个字典，包含图 G 中每个节点的离心率。如果指定节点 v，nx.eccentricity(G, v) 返回节点 v 的离心率。离心率仅适用于连通图。在不连通的图中，某些节点间可能不存在路径，因此无法计算离心率。在有向图中，应考虑路径的方向性，这可能影响离心率的计算结果。网络直径网络直径是图论中的一个重要概念，它指的是网络中所有节点对的最短路径中最长的那一个。对于连通图：# 计算直径 diameter = nx.diameter(G) print("网络的直径是:", diameter)如果图不连通（即存在不可达的节点对），这个函数将会抛出错误。在处理不连通图时，你可能需要先找出图中的连通分量，然后分别计算每个连通分量的直径。G = nx.Graph() G.add_edges_from([ (0, 1), (0, 2), (0, 7), (1, 2), (1, 3), (1, 4), (1, 7), (2, 3), (3, 4), (3, 5), (3, 7), (4, 5), (4, 6), (4, 7), (5, 6), (8,9),(9,10),(10,8) ]) nx.draw(G, node_size=500, with_labels=True) # 找出所有连通分量 connected_components = nx.connected_components(G) # 遍历每个连通分量，并计算其直径 for i, component in enumerate(connected_components): subgraph = G.subgraph(component) diameter = nx.diameter(subgraph) print(f"连通分量 {i+1} 的直径是: {diameter}")离心率（Eccentricity）和图的直径（Diameter）是图论中两个密切相关的概念，它们都与图中节点间距离的概念有关。离心率（Eccentricity）：公式表示为：$ \text{eccentricity}(v) = \max_{u \in G} d(v, u) $，其中 $ d(v, u) $ 是节点 $ v $ 到节点 $ u $ 的最短路径长度。图的直径（Diameter）：$ \text{diameter}(G) = \max_{v \in G} \text{eccentricity}(v) $。图的直径实际上是图中所有节点离心率的最大值。换句话说，直径是图中具有最大离心率的节点的离心率。一个节点的离心率告诉我们从该节点到达图中任何其他节点所需的最长距离。而图的直径则告诉我们在整个图中，任意两个节点之间的最长最短路径长度。在一个连通图中，最远的节点对的最短路径长度等于该图的直径。因此，离心率和直径都是衡量图中节点间距离的重要指标，但离心率着重于单个节点与其他所有节点的最远距离，而直径关注的是整个图中任意两个节点间的最远距离。效率、平均局部效率、全局效率节点对之间的效率（Efficiency between Node Pairs）:定义：节点对之间的效率是衡量这两个节点之间路径效率的指标。它是节点间最短路径的倒数。公式：$E_{ij} = \frac{1}{d(i, j)}$，其中，$ E_{ij} $ 是节点 $ i $ 和节点 $ j $ 之间的效率，$ d(i, j) $ 是它们之间的最短路径长度。NetworkX代码：nx.efficiency(G, 3, 7)。平均局部效率（Average Local Efficiency）:定义：平均局部效率是网络中每个节点的局部效率的平均值。局部效率是一个节点的邻居间效率的平均值。公式：$E_{local}(i) = \frac{1}{N_i} \sum_{j, h \in G_i, j \neq h} \frac{1}{d(j, h)}$，其中，$ G_i $ 是节点 $ i $ 的邻居子图，$ N_i $ 是子图中节点的数量，$ d(j, h) $ 是子图中节点 $ j $ 和 $ h $ 之间的最短路径长度。NetworkX代码：nx.local_efficiency(G)。全局效率（Global Efficiency）:定义：全局效率是网络中所有节点对效率的平均值。它反映了网络整体的交换或通信效率。公式：$E_{global} = \frac{1}{N(N-1)} \sum_{i \neq j \in G} \frac{1}{d(i, j)}$，其中，$ N $ 是图中的节点总数。NetworkX代码：nx.global_efficiency(G)。集聚系数、平均集聚系数与全局集聚系数演示实例：G = nx.Graph() G.add_edges_from([(1, 2), (1, 3), (1, 4), (1, 5), (2, 3), (2, 4), (2, 5), (3, 4), (3, 5), (4, 5)]) nx.draw(G, node_size=500, with_labels=True)在图论中，集聚系数（Clustering Coefficient）是用于度量节点周围邻居之间连接紧密程度的指标。它可以帮助我们了解图中的局部连接性。有三种主要的集聚系数：节点的集聚系数、平均集聚系数和全局集聚系数。节点的集聚系数是一个节点邻居之间实际存在的边数与可能存在的最大边数之比。对于一个无向图中的节点 $i$，其集聚系数 $C_i$ 的计算方式如下：$$ C_i = \frac{2 \times \text{实际存在的边数}}{\text{邻居节点数} \times (\text{邻居节点数} - 1)} $$这个值在 [0, 1] 范围内，表示节点的邻居之间连接紧密的程度。# 计算节点的集聚系数 node_clustering = nx.clustering(G) print("节点的集聚系数:", node_clustering)节点的集聚系数: {1: 1.0, 2: 1.0, 3: 1.0, 4: 1.0, 5: 1.0}以节点2为例，节点2的邻居节点是1、3、4、5。节点1、3、4、5可能存在的最大边数组合为$C_4^2 = \frac{4 \times 3}{2 \times 1} = 6$.$$ C_n^k = \frac{n!}{k!(n-k)!} $$节点1、3、4、5实际存在的连接情况是：1: [1, 3], [1, 4], [1, 5] 3: [3, 4], [3, 5] 4: [4, 5]共有6种。因此，节点2的聚集系数为1.平均集聚系数是图中所有节点的集聚系数的平均值。对于一个无向图 $G$，其平均集聚系数 $C$ 的计算方式如下：$$ C = \frac{1}{n} \sum_{i=1}^{n} C_i $$其中，$n$ 是图中的节点数。# 计算平均集聚系数 average_clustering = nx.average_clustering(G) print("平均集聚系数:", average_clustering)全局集聚系数是图中所有闭合三角形的数量与所有可能的闭合三角形的数量的比值。闭合三角形是由三个节点和它们之间的边组成的子图。对于一个无向图 $G$，其全局集聚系数 $C_{\text{global}}$ 的计算方式如下：$$ C_{\text{global}} = \frac{\text{3×闭合三元组的数量}}{\text{连接三元组的数量}} $$连接三元组（Open Triplet）是图论中的一个概念，它指的是在图中任选三个节点，其中至少有两个节点是相互连接的。这种三元组被称为“开放的”，因为它们不一定构成一个闭合的三角形。连接三元组是用来计算图的集聚系数的一个重要概念。在具体的定义中，连接三元组通常包含以下两种情况：闭合三元组（Closed Triplet）：这是图中的三个节点，它们之间的每一对节点都相互连接。换句话说，这三个节点形成了一个闭合的三角形。非闭合三元组：这也是图中的三个节点，但它们之间不是每一对节点都相互连接。这意味着虽然其中两个节点之间有边相连，但至少有一对节点之间没有直接的连线，因此不形成闭合的三角形。在计算图的全局集聚系数时，会考虑图中所有可能的连接三元组。全局集聚系数是闭合三元组数量与连接三元组总数量的比例。这个比例说明了在所有可能形成三角形的节点组合中，有多少实际形成了闭合的三角形。例如，在社交网络分析中，闭合三元组可能表示一种更强的社会关系，因为如果A认识B，B认识C，且A也认识C，这可能意味着这三个人之间有更紧密的社交联系。而非闭合的三元组则可能表示潜在的、未完全形成的社交联系。平均集聚系数关注的是每个节点的局部连接性，而全局集聚系数关注的是整个图中的全局连接性。# 计算全局集聚系数 global_clustering = nx.transitivity(G) print("全局集聚系数:", global_clustering)全局集聚系数: 1.0上面的示例图是一个完全图（任意两个顶点都是相连的），在完全图中，每一组三个节点都会形成一个闭合三角形，所以闭合三元组的数量等于连接三元组的数量。因此，全局集聚系数是 1.0。实际存在的闭合三角形的数量（闭合三元组的数量）：from itertools import combinations # 获取图中所有可能的闭合三角形数量 all_triangles_count = 0 for node in G.nodes(): # 遍历所有节点 neighbors = set(G.neighbors(node)) # 获取该节点的所有邻居节点 for pair in combinations(neighbors, 2): # 获取所有邻居节点的两两组合 if G.has_edge(pair[0], pair[1]): # 检查这个两两组合之间是否有边(因为它俩都是node的邻居节点，因此它俩肯定与node相连)，如果它俩相连，那么就说明node与pair[0]、pair[1]构成了三角形. all_triangles_count += 1上面的代码等同于：# 获取图中实际存在的闭合三角形数量(也需要除以3) actual_triangles_count = sum(nx.triangles(G).values()) 需要注意的是：all_triangles_count需要除以3得到的才是闭合三元组的数量，因为闭合三元组中的三个顶点会在不同的组合中进行重复计算。连接三元组的数量：open_triplets_count = 0 for node in G.nodes(): neighbors_count = len(list(G.neighbors(node))) open_triplets_count += (neighbors_count * (neighbors_count - 1)) // 2 # open_triplets_count 现在就是图中连接三元组的总数量从 $n$ 个数中随机挑选出 2 个的组合数可以用组合数学中的公式来计算，这个公式是 $ C(n, 2) $ 或者 $ \binom{n}{2} $，可以表示为：$$ C(n, 2) = \frac{n!}{2!(n-2)!} = \frac{n \times (n-1)}{2} $$图的连通性连通性描述的是图中节点之间是否存在路径相连的性质。一个图是连通的，意味着从图中的任意一个节点到另一个节点都存在路径。图结构连通性的简单判断G = nx.Graph() G.add_nodes_from([i for i in range(1, 8)]) G.add_edges_from([(1, 2), (1, 3), (2, 3), (4, 7), (5, 6), (5, 7), (6, 7)]) nx.draw(G, node_size=500, with_labels=True) # 查看图的连通性 nx.is_connected(G)什么叫做图的连通性？在无向图中，如果对于每一对不同的顶点 $u$ 和 $v$，都存在至少一条由边连接的路径从 $u$ 到 $v$，则该图是连通的。请注意这个概念并不等同于完全图的概念，完全图的概念是每一对不同的顶点 $u$ 和 $v$都直接相连，而连通图的每一对不同的顶点 $u$ 都可以达到 $v$，两者可以不直接相连。Fiedler值与图连通性的关系通过创建两个图来展示不同的连通性：# 创建具有较强连通性的图（低Fiedler值） G_strong = nx.Graph() G_strong.add_edges_from([(1, 2), (2, 3), (3, 4), (4, 1), (1, 3), (2, 4)]) # 创建具有较弱连通性的图（高Fiedler值） G_weak = nx.Graph() G_weak.add_edges_from([(1, 2), (2, 3), (3, 4)]) # 绘制这两个图 plt.figure(figsize=(12, 6)) # 绘制具有较强连通性的图 plt.subplot(1, 2, 1) nx.draw(G_strong, with_labels=True, node_color='lightblue', node_size=700) plt.title("Strongly Connected Graph") # 绘制具有较弱连通性的图 plt.subplot(1, 2, 2) nx.draw(G_weak, with_labels=True, node_color='lightgreen', node_size=700) plt.title("Weakly Connected Graph") plt.show()可以看出，左边的图具有更高的连通性，其应该具有更高的高Fiedler值，表明要将图分割成孤立的子图，需要切断更多的边。这通常表示图在整体上更加紧密连接，没有明显的弱连接点。右边的图具有更低的连通性，表明图可以通过切断较少的边来分割成不同的部分。这通常发生在图中存在一个或多个"瓶颈"区域，这些区域的边相对较少，是连接大的图区域的桥梁。G_strongG_weak通过对两个图结构的拉普拉斯矩阵进行特征值分解发现，左边图结构的Fiedler值为4.0，而右边的为0.586左右。因此，说明左边图结构的连通性更强。Thinking/ 思考两点自己的思考：对于一个管网系统来说，Fiedler似乎可以反应该管网的冗余性，Fiedler值越大说明该管网中网状结构越多。如果对一个管网进行切分，是否可以设计一种算法，来使得切出来的这个分区的Fiedler值最小，这样就保证了“切最少的管道”就可以得到了这个分区。遍历所有连通分量请查看上面“网络直径”一节中，如何遍历图中所有连通分量并获取其网络直径。节点对之间的连通性在图论中，节点连通性是指在不破坏图的连通性的情况下，可以从图中移除的最小节点数，以便特定的两个节点不再连通。换句话说，它是指在保持图的连通性的前提下，需要移除的节点数，才能使得这两个节点不再通过任何路径相连。获取所有节点对之间的连通性：from networkx.algorithms import approximation as approx # 计算所有节点对之间的节点连通性：独立路径数量 print(approx.all_pairs_node_connectivity(G))获取单个节点对之间的连通性：approx.local_node_connectivity(G, 6, 0)上面的代码返回结果是2，就代表图中节点6与节点0之间的连通性是2。也就说明，如果要破坏节点2与节点6之间的连通性，最少要移除两条边才可以。全局连通性全局连通性是指整个图的最小节点连通性，也就是在保持图连通的情况下，需要从图中移除的最少节点数，使得图变得不连通。approx.node_connectivity(G)函数将返回图 $G$ 的节点连通性，即使图变得不连通所需移除的最少节点数。这个值为整个图的一个全局度量，而不是针对特定节点对的。一个特殊的图结构——柏拉图式的八面体图：G2 = nx.octahedral_graph() # 柏拉图式的八面体图 nx.draw(G2, node_size=500, with_labels=True)print(approx.all_pairs_node_connectivity(G2)) print(approx.node_connectivity(G2))这个图结构所有节点间的连通性都是4。因此，其全局连通性也是4。NetworkX中关于连通分量的一些其他操作获取连通组件（子图）的数量：print(nx.number_connected_components(G))对所有联通分量按照节点数的多少排序，获取最大连通子图并绘图：Gcc = sorted(nx.connected_components(G), key=len, reverse=True) # 获取最大连通子图 largest_cc = G.subgraph(Gcc[0]) nx.draw(largest_cc, node_size=500, with_labels=True)复制连通分量，并进行修改：# largest_cc.add_edge(3,6) # 连通子图不支持修改：NetworkXError: Frozen graph can't be modified # 可以复制，在复制的图上修改 LCC = largest_cc.copy() LCC.add_edge(3,6) nx.draw(LCC, node_size=500, with_labels=True)链接密度Link Density（链接密度）通常是用于网络分析中的一个概念，特别是在图论和社交网络分析中。它描述了网络中链接（或边）的数量与网络中可能存在的最大链接数量之间的比例。换句话说，它衡量了网络的饱和度或连接紧密程度。$$ \text{Link Density} = \frac{2E}{N(N-1)} $$其中：$ E $ 是图中实际存在的边的数量。$ N $ 是图中节点的数量。对于一个有向图，由于边的方向性，计算公式略有不同：$$ \text{Link Density} = \frac{E}{N(N-1)} $$在这个公式中：$ E $ 仍然代表边的数量。$ N $ 代表节点的数量。因子“2”在无向图中用于考虑每条边连接两个不同的节点，而在有向图中则不需要。解释：当Link Density接近1时，表示网络中的节点之间几乎全部相互连接，这是一个高度紧密的网络。当Link Density接近0时，表示网络中的节点之间很少连接，这是一个稀疏网络。举例：import networkx as nx import matplotlib.pyplot as plt # 创建一个完全连接的网络（Link Density = 1） G1 = nx.complete_graph(5) # 创建一个稀疏连接的网络（较小的Link Density） G2 = nx.Graph() G2.add_edges_from([(0, 1), (1, 2), (3, 4)]) # 只添加几条边 # 绘制这两个网络 plt.figure(figsize=(12, 6)) # 绘制完全连接的网络 plt.subplot(121) nx.draw(G1, with_labels=True, node_color='lightblue', node_size=800) plt.title("Link Density = 1") # 绘制稀疏连接的网络 plt.subplot(122) nx.draw(G2, with_labels=True, node_color='lightgreen', node_size=800) plt.title("Link Density < 1") plt.show()在NetworkX中获取图结构的连接密度：print(nx.density(G1), nx.density(G2))其实可以发现，对于一个图，如果Link Density=1，那么其就是一个完全图。网格系数网格系数（Meshedness Coefficient）是图论中用于特别衡量平面图的一个指标。它测量图中有界面的数量与具有相同顶点数量的其他平面图可能拥有的面的数量之间的比例。这个系数有助于理解图中的网格程度或互连程度。网格系数为0表示树状结构，没有循环，代表最小的连通性。系数为1代表最大平面图，表示高度互连的结构，具有最大数量的面。在平面图中计算网格系数涉及计算有界面的数量（f）并将其与具有相同顶点数的图的最大可能面数进行比较。网格系数（M）的计算公式为：$$ M = \frac{f - 1}{2v - 5} $$其中：$ f $ 是图中有界面的数量；$ v $ 是图中顶点的数量。Tips / 提示什么叫做有界面的数量？在图论和网络分析中，"有界面的数量"（或简称为"面"）是指平面图中被边界包围的区域的数量。编写一个计算网格系数的函数：# 在NetworkX中计算网格系数 def calculate_meshedness(G): """计算网格系数""" # 计算有界面的数量：边的数量 - 顶点的数量 + 1 f = G.number_of_edges() - G.number_of_nodes() + 1 # 计算网格系数 v = G.number_of_nodes() if v > 2: M = (f - 1) / (2 * v - 5) else: M = 0 # 对于少于3个顶点的图，网格系数无意义 return M代数连通性代数连通性（Algebraic Connectivity）是图的拉普拉斯矩阵的第二小的特征值。在 NetworkX 中，可以使用 nx.algebraic_connectivity(G) 函数来计算图 $G$ 的代数连通性。这个函数内部会完成上述的计算步骤，并返回代数连通性的值。代数连通性通常也被称为Fiedler值。有关其更详细的解释，请看下面“特征值分解的意义”章节。点中心性点中心性（degree centrality）是最基本的中心性度量之一，它反映了节点的连接数相对于网络中可能的连接数的比例。在无向图中，点中心性简单地等于某个节点的邻接节点数（或者说是边数）除以可能的最大邻接节点数（即网络中节点总数减去1）。度中心性对于网络中的节点 $ v $，其度中心性（Degree Centrality） $ C_D(v) $ 可以通过以下公式计算：$$ C_D(v) = \frac{\text{节点 } v \text{ 的度数}}{N - 1} $$其中，$ N $ 是网络中节点的总数。在 NetworkX 中，可以使用 degree_centrality(G) 函数来计算所有节点的点中心性。这个函数返回一个字典，其中键是节点，值是对应节点的点中心性。其他几种中心性特征向量中心性（Eigenvector Centrality）:定义：特征向量中心性考虑了一个节点的邻居的重要性。一个节点的特征向量中心性高，如果它连接到许多其他具有高中心性的节点。公式：$ C_E(v) = \frac{1}{\lambda} \sum_{t \in M(v)} C_E(t) $其中，$ C_E(v) $ 是节点 $ v $ 的特征向量中心性，$ M(v) $ 是与 $ v $ 直接相连的节点集合，$ \lambda $ 是特征值。接近度中心性（Closeness Centrality）:定义：接近度中心性基于节点到网络中所有其他节点的平均最短路径长度。距离越短，接近度中心性越高。公式：$ C_C(v) = \frac{N-1}{\sum_{t \neq v} d(v, t)} $，其中$ C_C(v) $ 是节点 $ v $ 的接近度中心性，$ d(v, t) $ 是节点 $ v $ 到 $ t $ 的最短路径长度。介数中心性（Betweenness Centrality）:定义：介数中心性是基于节点在网络中所有最短路径中出现的频率。一个节点的介数中心性高，如果更多的最短路径通过它。公式：$ C_B(v) = \sum_{s \neq v \neq t} \frac{\sigma_{st}(v)}{\sigma_{st}} $，其中，$ C_B(v) $ 是节点 $ v $ 的介数中心性，$ \sigma_{st} $ 是节点 $ s $ 到 $ t $ 的所有最短路径的数量，$ \sigma_{st}(v) $ 是经过节点 $ v $ 的那些路径的数量。PageRank中心性（PageRank Centrality）:定义：PageRank是衡量网页重要性的算法，但也可用于网络中。一个节点的PageRank高，如果有许多其他高PageRank的节点指向它。公式：$ PR(v) = \frac{1-d}{N} + d \sum_{t \in M(v)} \frac{PR(t)}{L(t)} $，其中，$ PR(v) $ 是节点 $ v $ 的PageRank值，$ M(v) $ 是指向 $ v $ 的节点集合，$ L(t) $ 是节点 $ t $ 的出度，$ d $ 是阻尼系数（通常设为0.85）。如何在NetworkX中获得：# 度中心性（Degree Centrality） nx.degree_centrality(G) # 特征向量中心性（Eigenvector Centrality） nx.eigenvector_centrality(G) # 接近度中心性（Closeness Centrality） nx.closeness_centrality(G) # 介数中心性（Betweenness Centrality） nx.betweenness_centrality(G) # PageRank中心性（PageRank Centrality） nx.pagerank(G)中心点优势Central-point dominance（中心点优势）是网络分析中的一个度量，用来量化网络中最中心节点的控制力相对于网络中其他节点的控制力。它是由Linton Freeman在1977年提出的。这个度量是基于中心性的概念，特别是点中心性（point centrality）。在某个网络中，中心点优势是由网络中点中心性最高的节点（即最中心的节点）的点中心性值，减去网络中所有节点的点中心性平均值，再除以最大可能的差值来计算的。计算步骤如下：计算网络中每个节点的点中心性（可以是接近中心性、介数中心性或特征向量中心性等）。找到具有最高点中心性的节点。计算该节点的点中心性与所有节点的点中心性平均值的差值。除以最大可能的差值，这通常是最中心节点的点中心性减去最小可能的点中心性（通常是1/n，其中n是节点总数）。公式表示为：$$ CPD = \frac{C_{max} - \overline{C}}{Max(C_{max} - C_i)} $$其中：$ CPD $ 是中心点优势。$ C_{max} $ 是最高点中心性。$ \overline{C} $ 是所有节点点中心性的平均值。$ Max(C_{max} - C_i) $ 是最高点中心性与任意节点中心性之差的最大可能值。在 NetworkX 中没有直接计算中心点优势的函数，但可以使用已有的中心性计算函数和一些额外的计算步骤来求解。中心点优势的具体含义包括：网络权力集中度：CPD较高意味着网络中有一个或几个节点拥有远高于其他节点的中心性，表明网络中的权力或影响力可能非常集中。脆弱性和稳健性：当一个网络的CPD很高时，网络可能对最中心节点的移除或故障非常敏感。这可能会影响网络的整体功能和稳定性。信息流动和效率：一个中心点优势很高的网络可能会有更有效的信息流动，因为中心节点可以作为信息流动的枢纽。然而，这也可能导致信息过度集中，使得网络对于中心节点的依赖度增加。Thinking/ 思考在给水管网的图结构中，中心点优势（Central Point Dominance, CPD）可以有以下几种现实意义：关键节点的识别：高中心点优势意味着给水管网中存在一个或几个关键节点，这些节点对整个网络的运作至关重要。例如，这可能是主要的水源地、水处理中心或主要的分配点。网络脆弱性分析：如果给水管网的中心点优势很高，这可能表明网络对单个节点或几个节点的依赖性很大。这样的网络可能对于这些关键节点的故障或损坏特别敏感，从而影响整个网络的供水能力。优化和升级规划：通过分析中心点优势，可以识别那些对整个给水管网至关重要的节点，从而有助于优先考虑这些节点的维护、升级或加固，以提高网络的整体效率和可靠性。应急响应和风险管理：在应对突发情况（如管道破裂、污染事件等）时，了解中心点优势可以帮助决策者优先处理那些对整个网络影响最大的节点，从而更有效地减轻整体影响。提高供水效率：通过优化中心节点的运作，可以提高整个给水系统的供水效率和效果，确保水资源的合理分配和使用。因此，在给水管网的管理和优化中，理解和分析中心点优势是非常重要的，它有助于指导网络设计、运维、危机管理和改进策略。实际数据分析表明，CPD与管网节点数呈现略微正相关关系，这表明节点数越多的管网中，网络对单个节点或几个节点的依赖性越大。这是为什么呢？不科学啊？不应该是越大的管网抵抗风险的能力越强吗？关键比例"Critical Fraction"（关键比例）是网络科学中的一个概念，通常用于描述网络在遭受随机故障或有意攻击时维持功能的能力。具体来说，它指的是导致网络失去大规模连通性或功能的最小比例的节点或边的移除。这个概念在理解和评估网络的鲁棒性和脆弱性方面非常重要。举个🌰，创建两个随机图结构：import networkx as nx # 创建两个网络 G_A = nx.erdos_renyi_graph(100, 0.05) # 网络A G_B = nx.watts_strogatz_graph(100, 4, 0.1) # 网络B # 准备绘制网络结构图和度分布统计图 fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 6)) # 在左边绘制网络图 nx.draw(G_A, ax=ax1, node_size=50, with_labels=False) nx.draw(G_B, ax=ax2, node_size=50, with_labels=False) ax1.set_title("Graph A") ax2.set_title("Graph B") plt.show()两个随机图结构定义一个计算关键比例的函数：import random def calculate_critical_fraction(G: nx.Graph, iter_num: int = 100) -> float: """ 图结构关键比例计算函数 :param G: 图结构对象, nx.Graph格式 :param iter_num: 随机循环数量, 通过多次重复移除节点的实验, 并计算每次实验的关键比例, 然后取这些实验的平均值. :return: 多次计算关键比例的平均值,范围[0, 1] """ ratio = [] for _ in range(iter_num): # 复制原始网络 G_temp = G.copy() # 计算原始网络中最大连通分量的节点数 original_size = max(len(c) for c in nx.connected_components(G_temp)) for i in range(len(G_temp)): # 随机选择并移除一个节点 node_to_remove = random.choice(list(G_temp.nodes())) G_temp.remove_node(node_to_remove) # 计算移除节点后网络中最大连通分量的节点数 largest_cc_size = max(len(c) for c in nx.connected_components(G_temp)) if largest_cc_size <= original_size / 2: # 如果移除节点后, 网络中最大连通分量的节点数小于原始网络中的1/2 ratio.append(i) break # 返回平均关键比例 return sum(ratio) / len(ratio) / len(G)对比两个图结构的关键比例：critical_fraction_A = calculate_critical_fraction(G_A.copy()) critical_fraction_B = calculate_critical_fraction(G_B.copy()) print("Critical Fraction of Network A:", critical_fraction_A) print("Critical Fraction of Network B:", critical_fraction_B)Critical Fraction of Network A: 0.45640000000000003Critical Fraction of Network B: 0.3701比较两个网络的平均关键比例，发现$CF_A>CF_B$，说明网络B在面临随机故障时更加脆弱。Thinking/ 思考感觉“关键比例”这个概念不适合用于对比两个给水管网在面临随机故障时的脆弱性。因为给水管网中节点的重要性并是仅仅是看哪个节点的度比较大（给水管网中节点的度都在1\~4之间）。关键比例仅仅对比两个管网在随机减少节点后，具有最大连通分量的节点数量，这对于给水管网来说并不合理。K核值节点的 k-核（k-core）是一个网络分析中的概念，用于识别网络中的一个最大子图，在这个子图中，每个节点至少与 k 个其他节点相连。计算一个网络中所有节点的 k-核值的过程涉及到逐步移除度数小于 k 的节点，直到所有剩余的节点度数都大于或等于 k。举个例子就明白了，创建一个图结构：# 创建图 G = nx.Graph() G.add_node('W') G.add_edges_from([('A', 'B'), ('A', 'C'), ('A', 'D'), ('B', 'C'), ('C', 'D'), ('E', 'F'), ('B', 'D'), ('B', 'G'), ('G', 'C')]) nx.draw(G, node_size=500, with_labels=True)分别获取并绘制上面图结构中K核值为1、2、3的节点：import matplotlib.pyplot as plt plt.rcParams.update({ 'font.size': 20, 'font.family': ['Times New Roman', 'SimSun'] }) # 准备绘制网络结构图和度分布统计图 fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(16, 6)) k_1 = nx.k_core(G, k=1) k_2 = nx.k_core(G, k=2) k_3 = nx.k_core(G, k=3) # 在左边绘制网络图 nx.draw(k_1, ax=ax1, node_size=500, with_labels=True, node_color='r') nx.draw(k_2, ax=ax2, node_size=500, with_labels=True, node_color='g') nx.draw(k_3, ax=ax3, node_size=500, with_labels=True, node_color='b') ax1.set_title('K=1') ax2.set_title('K=2') ax3.set_title('K=3') plt.show()这个时候就很容易明白K核值的定义了：nx.k_core(G, k=1)返回的是图结构中度≥1的节点；nx.k_core(G, k=2)返回的是图结构中度≥2的节点；nx.k_core(G, k=3)返回的是图结构中度≥3的节点。打印单个节点的K核值：for node, k_core in k_core_values.items(): print("节点", node, "的K核值是", k_core)我们会发现打印结果与上面的分析结论是一样的。度同配性如果网络中度数相似的节点倾向于相互连接，即高度数的节点倾向于与其他高度数的节点相连，低度数的节点倾向于与其他低度数的节点相连，这种现象称为度同配（Assortative Mixing）。度同配性（Degree Assortativity）度同配性系数是一个介于 -1 和 1 之间的数值，用来量化度同配性的程度。系数接近 1 表示强同配性，即网络中类似度数的节点倾向于相互连接。系数接近 -1 表示强异配性，即网络中不同度数的节点倾向于相互连接。在 networkx 中，可以使用 nx.degree_assortativity_coefficient(G) 来计算图 G 的度同配性系数。这个函数会考虑图中每个节点的度数，并计算整个网络的度同配性。度同配性系数的计算涉及到评估网络中节点的度数（即每个节点的连接数）之间的相关性。具体计算方法基于皮尔逊相关系数，但是针对的是节点的度数而不是节点的特征或标签。计算步骤：准备数据：对于每条边，考虑其两个端点的度数。例如，如果一条边连接度为2的节点和度为3的节点，则这条边对应于（2,3）和（3,2）两组度数对。计算平均度数和方差：计算网络中所有节点度数的平均值和方差。计算协方差：对于所有边的度数对，计算它们的协方差。协方差测量两个变量（在这里是度数）之间的总体偏差。应用皮尔逊相关系数公式：使用下面的公式计算度同配性系数：$r = \frac{\sum_{xy} xy(e_{xy} - a_x a_y)}{\sigma_a \sigma_b}$。其中，$ e_{xy} $ 是连接度数为 $ x $ 和 $ y $ 的节点的边的比例，$ a_x $ 是度数为 $ x $ 的节点的比例，$ \sigma_a $ 和 $ \sigma_b $ 是度分布的标准差。桥连接（这个🌉名词实在是不知道怎么翻译……）在图论中，桥是指图中的一条边，移除它会导致图的连通性分量增加。换句话说，桥是连接图中两个部分的唯一边，移除它会使这两部分变得不连通。举个例子：G_strong = nx.Graph() G_strong.add_edges_from([ (0, 1), (0, 2), (0, 7), (1, 2), (1, 3), (1, 7), (2, 3), (3, 5), (3, 7), (4, 5),(4, 6), (5, 6) ]) G_weak = G_strong.copy() G_weak.remove_edge(3, 5) # 绘制这两个图 plt.figure(figsize=(12, 6)) plt.subplot(1, 2, 1) nx.draw(G_strong, with_labels=True, node_color='lightblue', node_size=700) plt.title(f"Graph with a bridge{list(nx.bridges(G_strong))[0]}") plt.subplot(1, 2, 2) nx.draw(G_weak, with_labels=True, node_color='lightgreen', node_size=700) plt.title("Graph with no bridges") plt.show()上面的左图中，边 (3, 5)就是一个Bridge，因为当切断这条边时，该图就会变成两个连通分量。桥连接分析网络的结构稳定性和脆弱性时非常有用，特别是在识别网络中的关键连接或潜在的脆弱点方面。在NetworkX 库中：nx.has_bridges(G)：返回一个Bool值，返回 True，则表示图 G 中至少有一条边是桥；如果返回 False，则表示图中没有桥。nx.bridges(G)：返回一个可迭代对象，包含了网络中的所有桥。补充个知识点，还有个东西叫做“本地桥”，在NetworkX 库中使用nx.local_bridges(G)即可获取。还是举个例子来说明吧：import networkx as nx # 创建一个图 G = nx.Graph() G.add_edges_from([ ('A', 'B'), ('B', 'C'), ('C', 'D'), ('D', 'A'), ]) # 找出所有桥 print("图中的桥:") for bridge in nx.bridges(G): print(bridge) # 找出所有本地桥 print("\n图中的本地桥:") for local_bridge in nx.local_bridges(G, with_span=True): print(local_bridge) nx.draw(G, with_labels=True, node_color='lightblue', node_size=700)在上面的图结构中，没有桥，因为断开任何一条边也不能让图变成两个连通分量。所谓“本地桥”是指，桥两端的节点在桥断开后，没有共同的节点相连。将上面的图结构中，添加一个点E，使其与点A、点D相连，那么桥 (A, D ) 就不是本地桥了：G.add_edges_from([ ('E', 'A'), ('E', 'D') ]) # 找出所有桥 print("图中的桥:") for bridge in nx.bridges(G): print(bridge) # 找出所有本地桥 print("\n图中的本地桥:") for local_bridge in nx.local_bridges(G, with_span=True): print(local_bridge) nx.draw(G, with_labels=True, node_color='lightblue', node_size=700)图的特征值分解下面的讲解以Introduction to Graph Signal Processing中P19页的Fig. 10（图1，连通图）、P14页的Fig. 7（图2，非连通图）为例。G_strong = nx.Graph() G_strong.add_edges_from([ (0, 1), (0, 2), (0, 7), (1, 2), (1, 3), (1, 4), (1, 7), (2, 3), (3, 4), (3, 5), (3, 7), (4, 5), (4, 6), (4, 7), (5, 6) ]) G_weak = nx.Graph() G_weak.add_edges_from([ (0, 1), (0, 2), (1, 2), (3, 4), (3, 5), (3, 7), (4, 5), (4, 6), (4, 7), (5, 6) ]) # 绘制这两个图 plt.figure(figsize=(12, 6)) plt.subplot(1, 2, 1) nx.draw(G_strong, with_labels=True, node_color='lightblue', node_size=700) plt.title("Strongly Connected Graph") plt.subplot(1, 2, 2) nx.draw(G_weak, with_labels=True, node_color='lightgreen', node_size=700) plt.title("Weakly Connected Graph") plt.show()特征值分解的意义邻接矩阵更倾向于揭示图的整体连通性，而拉普拉斯矩阵的特征值和特征向量则更多地用于研究图的局部结构特征，如社区结构或者图的连通分量。为什么拉普拉斯矩阵的第一个特征值总是0？拉普拉斯矩阵的第一个特征值总是0，其对应的特征向量是一个所有元素都相同的向量。这代表图中所有节点的均匀分布。拉普拉斯矩阵 $ L $ 的第一个特征值总是 0 的原因与拉普拉斯矩阵的定义和图的性质有关。拉普拉斯矩阵 $ L $ 定义为度矩阵 $ D $ 减去邻接矩阵 $ A $，即 $ L = D - A $。这里是为什么其第一个特征值总是 0：和为零的行（列）：在拉普拉斯矩阵中，每一行的和（以及每一列的和，因为 $ L $ 是对称的）都是 0。这是因为每一行的非对角元素（即 $-A$ 的部分）与对角线上的元素（即 $ D $ 的部分，它是节点的度数）相加抵消。换句话说，对于每个节点 $ i $，它的度数 $ D_{ii} $ 等于与它相连的边的数量，而 $ -A $ 中的元素表示与节点 $ i $ 相连的邻居节点，因此当你在一行中把这些值加起来时，结果是 0。恒等特征向量：由于每行的和为零，这意味着全 1 向量（所有元素都是 1 的向量）是拉普拉斯矩阵的一个特征向量，因为 $ L \times \textbf{1} = \textbf{0} $（其中 $\textbf{1}$ 是全1向量，$\textbf{0}$ 是零向量）。这意味着乘以全 1 向量后，每行的元素相加都得到 0，这与零向量相匹配。特征值 0 的意义：特征值 0 对应的特征向量表示图中所有节点的均匀分布，它反映了图的连通性。对于连通图，特征值 0 是唯一的，它的代数重数（特征值的重复次数）等于图的连通分量的数量。因此，对于连通图，特征值 0 的代数重数是 1。如果图不是完全连通的，特征值 0 的代数重数将等于图的连通分量数量。简而言之，拉普拉斯矩阵的每一行和每一列的和为零这个属性保证了第一个特征值必定是0。这与图的基本性质密切相关，特别是与图的连通性有关。其他特征值的意义：第二个特征值（Fiedler值）和特征向量：第二个最小的特征值通常被称为Fiedler值，它在图的谱聚类和社区检测中非常重要。Fiedler值的大小可以表示图的连通性：Fiedler值越小，图的连通性越弱。对应的Fiedler向量可以用来识别图中的社区或集群。通过分析Fiedler向量的组件，可以将图划分为不同的部分，其中每个部分相对内部紧密连接，而与其他部分的连接较少。其他特征值和特征向量：更高的特征值和对应的特征向量可以揭示图的更复杂的结构特征，如多层次的社区结构。在某些情况下，这些特征向量也用于嵌入技术，将图的节点映射到低维空间，以便于可视化和进一步分析。特征值分解下面对上面的示例图进行拉普拉斯矩阵的特征值分解。拉普拉斯矩阵首先，获取图结构的拉普拉斯矩阵：# 获取拉普拉斯矩阵 L = nx.laplacian_matrix(G).toarray() L$$ \begin{bmatrix} 3 & -1 & -1 & -1 & 0 & 0 & 0 & 0 \\ -1 & 5 & -1 & -1 & -1 & -1 & 0 & 0 \\ -1 & -1 & 3 & 0 & -1 & 0 & 0 & 0 \\ -1 & -1 & 0 & 4 & -1 & -1 & 0 & 0 \\ 0 & -1 & -1 & -1 & 5 & -1 & -1 & 0 \\ 0 & -1 & 0 & -1 & -1 & 5 & -1 & -1 \\ 0 & 0 & 0 & 0 & -1 & -1 & 3 & -1 \\ 0 & 0 & 0 & 0 & 0 & -1 & -1 & 2 \\ \end{bmatrix} VS \begin{bmatrix} 2 & -1 & -1 & 0 & 0 & 0 & 0 & 0 \\ -1 & 2 & -1 & 0 & 0 & 0 & 0 & 0 \\ -1 & -1 & 2 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 3 & -1 & -1 & -1 & 0 \\ 0 & 0 & 0 & -1 & 4 & -1 & -1 & -1 \\ 0 & 0 & 0 & -1 & -1 & 3 & 0 & -1 \\ 0 & 0 & 0 & -1 & -1 & 0 & 2 & 0 \\ 0 & 0 & 0 & 0 & -1 & -1 & 0 & 2 \\ \end{bmatrix} $$需要注意的是，返回的拉普拉斯矩阵的行列顺序并不与图中的顺序相同，矩阵中的行列数据是按照节点的添加顺序来的。如何查看节点的顺序：list(G.nodes()) # [0, 1, 2, 7, 3, 4, 5, 6]对于图1来说，因为节点7添加的早，所以排在节点3之前。也就是说，拉普拉斯矩阵中第4行代表的是第7个元素的连接情况。特征值#求矩阵特征值以及特征向量 eig_value, eig_vector = np.linalg.eig(L) eig_value特征值及其可视化：$$ [4.44089210e-16, 1.13357211e+00, 3.05427452e+00, 3.31740297e+00, 4.00000000e+00, 5.67905993e+00, 6.47332881e+00, 6.34236165e+00] $$def autolabel(rects): for rect in rects: height = rect.get_height() height_pos = height + .1 if height >= -0.1 else height - .1 plt.text( rect.get_x() + rect.get_width() / 2, height_pos, '{:.5f}'.format(height), size=10, family="Times new roman", horizontalalignment='center', verticalalignment='center' ) fig = plt.bar([i for i in range(len(eig_value))], np.sort(eig_value)) autolabel(fig) plt.title('Eigenvalues') plt.show()图1图2图1有一个接近零的特征值，表明图是连通的。图2特征值有两个接近于零的值，这与图中的两个连通分量相对应。特征值为0的数量恰好等于图的连通分量的数量。图2的 Fiedler 值（最小非0值）为1.58，该 Fiedler 值对应的特征向量是：$$ \begin{bmatrix} 0. \\ 0. \\ 0. \\ -0.4472136 \\ -0.4472136 \\ -0.4472136 \\ -0.4472136 \\ -0.4472136 \\ \end{bmatrix} $$说明，该 Fiedler 值代表的是3、4、5、6、7节点对应的连通分量的连通性。总结：图1的连通性更强，因为其特征值中仅有一个为0；图2包含两个连通分量，因为其特征值中包含两个0。图2中3、4、5、6、7节点组成的连通分量的连通性要高于图1整体的连通性。因为图2中3、4、5、6、7节点组成的连通分量的 Fiedler 值为1.58，大于图1整体的连通分量1.13。Fiedler 值也被称为代数连通性，它反映了图的稀疏性或稠密性。一个较高的Fiedler值意味着图具有较强的连通性，而较低的Fiedler 值则意味着图可能容易被分割成独立的组件。Thinking/ 思考管网数据分析表明：随着管网中节点数的增多，Fiedler 值（代数连通性）越来越小，说明越大的管网环状比例越小。是不是可以这样理解为：大部分管网新建的时候都是城区环状管网，随着城市的发展，开始逐渐向周围以支状管网扩散，这就导致了大型管网的环状比例都较小。特征向量特征向量（图1）：$$ \begin{bmatrix} 0.35355339 & 0.40809973 & -0.38534069 & -0.32322928 & 0.61237244 & 0.10317586 & -0.24182446 & 0.10660991 \\ 0.35355339 & 0.21821011 & 0.07059821 & -0.06640201 & -0.20412415 & 0.58748394 & 0.64631494 & 0.11603433 \\ 0.35355339 & 0.36261755 & -0.3221172 & 0.67753284 & -0.20412415 & -0.23709192 & -0.00834996 & -0.28766179 \\ 0.35355339 & 0.18086105 & 0.27243317 & -0.5085369 & -0.20412415 & -0.62680632 & 0.20197088 & -0.18470141 \\ 0.35355339 & 0.05048967 & 0.33222523 & 0.17458035 & -0.20412415 & -0.05547633 & -0.37548832 & 0.7388255 \\ 0.35355339 & -0.15837433 & 0.24016423 & -0.13207484 & -0.20412415 & 0.41726191 & -0.52854257 & -0.52883224 \\ 0.35355339 & -0.40809973 & 0.38534069 & 0.32322928 & 0.61237244 & -0.10317586 & 0.24182446 & -0.10660991 \\ 0.35355339 & -0.65380405 & -0.59330365 & -0.14509945 & -0.20412415 & -0.08537128 & 0.06409502 & 0.14633561 \\ \end{bmatrix} $$将特征值从小到大排序，然后可视化特征向量：def plot_eigenvectors(value, vector): vector = vector[:, np.argsort(value)] value = np.sort(value) # plt.figure(dpi=100) fig, ax_arr = plt.subplots(len(value), 1, sharex='col', sharey='row', figsize=(5, 2 * len(value))) fig.subplots_adjust(hspace=0.05, wspace=0.03) for i in range(len(value)): ax_arr[i].stem([i for i in range(len(value))], vector[:, i]) ax_arr[i].set_ylabel(f'$u_{i}}}(n)$') plot_eigenvectors(eig_value, eig_vector)图1图2上面两个图是按照特征值排序后的特征向量，下面依次来解释特征向量有着什么样的含义：第一个特征向量的含义：对于任何图，拉普拉斯矩阵的第一个特征值总是0。对于零特征值，其对应的特征向量（通常是一个所有元素都相等的向量）实际上表示了图中所有节点的一种“平均”状态或者均匀状态。对于连通图（图1），这反映了图中所有节点在某种意义上是等价的，没有任何节点是孤立的。对于非连通图（图2），对于非连通图，拉普拉斯矩阵的零特征值的个数等于图的连通分量的数量，每个零特征值的特征向量反映了对应连通分量中的节点关系。例如图2的第一个特征向量有三个值不为0，那么就说明图2的第一个连通分量中有三个元素；第二个特征向量有四个值不为0，那么就说明第二个连通分量中有四个元素。总的来说，连通图的第一个特征向量不反映任何局部结构或节点间的差异，而是显示了图作为一个整体的特性，即全局特征。Fiedler向量（第二个特征向量）：在Fiedler向量中，正值和负值分别代表了图中的不同群体，这种划分往往揭示了图的潜在社区结构，其中社区内部的节点比跨社区的节点更紧密地连接。例如图1中，Fiedler向量的前5个值为负值，后3个为正值，这说明图1更容易分为两个簇：前5个节点（0、1、2、7、3）、后3个节点（4、5、6）。Thinking / 思考除了可以根据Fiedler向量对图中节点进行简单的正负划分外，还可以根据不同值的接近程度进行更为细致的划分。例如，假设一个Fiedler向量为[-100,-99,-1,2,3,6,7,8]：基本的二分法：首先，根据Fiedler向量的正负符号，可以将图分为两个主要群体：第一个群体：包含节点1、2、3（对应于Fiedler向量中的值为-100, -99, -1）；第二个群体：包含节点4、5、6、7、8（对应于Fiedler向量中的值为2, 3, 6, 7, 8）。更细致的分组：在每个主要群体内，可以根据Fiedler向量值的接近程度进行更细致的分组：在第一个主要群体中，节点1和2（Fiedler向量值为-100和-99）更接近彼此，而节点3（值为-1）相对较远；在第二个主要群体中，节点4和5（Fiedler向量值为2和3）更接近彼此，而节点6、7、8（值为6, 7, 8）更接近彼此。如果Fiedler向量的这种特性，在管网的图结构中，是否可以对管网中的节点进行分类，分为不同的簇，进而达到简化管网的目的？NetworkX中获取特征值（谱）获取邻接矩阵的特征值：print(np.real(nx.adjacency_spectrum(G)))计算拉普拉斯矩阵的特征值：print(np.real(nx.laplacian_spectrum(G)))NetworkX中常见图（网络）生成器完全图（全连接图）# 生成包含n个节点的完全图（全连接图） n = 10 G1 = nx.complete_graph(10) pos = nx.circular_layout(G1) nx.draw(G1, pos, node_size=500, node_color='red', with_labels=True)星型图# 生成包含n个节点的星型图 n = 10 G2 = nx.star_graph(n-1) pos = nx.spring_layout(G2) nx.draw(G2, pos, node_size=500, node_color='red', with_labels=True)ER随机图# 生成包含n个节点，连边概率为p的ER随机图 n, p = 20, 0.2 G3 = nx.erdos_renyi_graph(n, p) nx.draw(G3, node_size=500, node_color='red', with_labels=True)WS小世界网络# 生成包含n个节点，重连概率为p的WS小世界网络 n, k, p = 20, 4, 0.2 # 当p=0时，便退化成了k近邻规则网络 G4 = nx.watts_strogatz_graph(n, k, p) pos = nx.circular_layout(G4) nx.draw(G4, pos, node_size=500, node_color='red', with_labels=True)BA无标度网络# 生成包含n个节点，参数m=2的BA无标度网络 n, m = 50, 2 G5 = nx.barabasi_albert_graph(n, m) pos = nx.spring_layout(G5) nx.draw(G5, pos, node_size=500, node_color='red', with_labels=True)外部变量（数据）创建图结构dict对象将字典数据转化为网络：d = {0: {1: {"weight": 1}, 2: {"weight": 2), 1: {2: {"weight": 3.5}}} G = nx.Graph(d) # 或者 # G = nx.from_dict_of_dicts(d) edge_width = nx.get_edge_attributes(G, "weight") # print(width) nx.draw(G, node_size=500, node_color='red', with_labels=True, width=list(edge_width.values()))反过来把网络转换成字典数据：print(nx.to_dict_of_dicts(G))list对象转化为list对象：edgelist = nx.to_edgelist(G) print(edgelist)从list格式创建：G = nx.from_edgelist(edgelist)array对象转化为array对象（邻接矩阵）格式：nx.adjacency_matrix(G)从array对象（邻接矩阵）格式转化：nx.from_numpy_array(A)DataFrame对象df = pd.read_excel("edges.xlsx") G = nx.from_pandas_edgelist(df, 'source', 'target','weight', create_using = nx.Graph()) # 若为有向图，create_using = nx.DiGraph() # 若为无权无向网络则：G = nx.from_pandas_edgelist(df, 'source', 'target', create_using = nx.Graph()) edge_width = nx.get_edge_attributes(G, "weight") nx.draw(G, node_size=500, node_color='red', with_labels=True, width=list(edge_width.values()))NetworkX可视化常用属性值设置节点大小与节点颜色：nx.draw(G, node_size=100, node_color="red")更加复杂的格式：# 当设置的属性较多时，可以将其保存在字典中，以**不定长参数传入 options = { 'pos': nx.spring_layout(G), 'node_size': 300, 'node_color': "red", 'edge_color': "gray", 'width': 1.0, # 连边粗细 'with_labels': True, } nx.draw(G, **options)布局方式networkx 库提供了多种布局算法来可视化图，每种布局算法都有其独特的特性和适用场景。以下是一些常用的布局算法：Spring Layout（Fruchterman-Reingold Force-directed Algorithm）:实现了所谓的“弹簧”布局算法（也称为Fruchterman-Reingold算法），它是一种力导向图绘制算法。使用方法：nx.spring_layout(G)Circular Layout:将节点放置在一个圆上，通常用于突出循环结构。使用方法：nx.circular_layout(G)Random Layout:将节点随机放置在画布上，用于快速生成节点位置。使用方法：nx.random_layout(G)Spectral Layout:基于图的拉普拉斯特征向量（使用拉普拉斯矩阵的第二小和第三小的特征值对应的特征向量，这两个特征向量用于在二维空间中定义节点的位置），尝试将连接的节点放置得更近。使用方法：nx.spectral_layout(G)Shell Layout:将图中的节点排列在同心圆上，适用于层级结构。使用方法：nx.shell_layout(G)Kamada-Kawai Layout:基于路径长度的力导向算法，旨在使路径长度在图上的表现与其在几何空间中的长度相匹配。使用方法：nx.kamada_kawai_layout(G)Bipartite Layout:专门用于二分图的布局，将节点分为两组，每组形成一列。使用方法：nx.bipartite_layout(G, nodes)Planar Layout:如果图是平面图，这个算法会找到一个平面布局（即让所有的边都不相交）。使用方法：nx.planar_layout(G)Circular LayoutSpectral LayoutBipartite LayoutPlanar Layout点定位与边权重的设置创建一个有权无向图：G = nx.Graph() edge_list = [(0, 1, 2), (0, 2, 8), (0, 3, 1), (1, 2, 6), (1, 4, 1), (2, 3, 7), (2, 4, 5), (2, 5, 1), (2, 6, 2), (3, 6, 9), (4, 5, 3), (4, 7, 8), (5, 6, 4), (5, 7, 6), (6, 7, 3)] G.add_weighted_edges_from(edge_list) nx.draw(G, node_size=500, with_labels=True)自定义节点的位置并将边的权重作为图中边的宽度进行绘制：# 自定义各个节点的坐标 pos = {0: (-2, 0), 1: (-1, 1), 2: (-1, 0), 3: (-1, -1), 4: (1, 1), 5: (1, 0), 6: (1, -1), 7: (2, 0)} # 提取边的权重并作为标签 e_labels = nx.get_edge_attributes(G, 'weight') e_width = [G.get_edge_data(*e)['weight'] for e in G.edges()] options = { 'pos': pos, 'node_size': 500, 'node_color': "red", 'edge_color': "gray", 'width': e_width, 'with_labels': True, } nx.draw(G, **options) nx.draw_networkx_edge_labels(G, pos, edge_labels=e_labels)以节点度区分颜色G = nx.barabasi_albert_graph(20, 2) # 绘制网络图，按度值大小设定节点大小和颜色 # 设置节点大小与度成正比 nodesize = [G.degree(i) * 100 for i in G.nodes()] node_colors = [G.degree(i) for i in G.nodes()] options = { 'pos': nx.spring_layout(G), 'node_size': nodesize, 'node_color': node_colors, 'cmap': plt.cm.cool, # 设置节点colormap 'edge_color': "gray", 'with_labels': True, } nx.draw(G, **options) plt.show()更高阶的玩法更多NetworkX可视化案例见：https://networkx.org/documentation/stable/auto_examples/index.html；Gephi是网络图可视化领域中无敌的软件：https://gephi.org/。实战——FackBook社交网络分析说明这个案例来自NetworkX官方网站，通过对10个FackBook朋友圈社交网络的分析，来复习下图论中的一些基本概念以及相关NetworkX包的使用。请注意：每个节点代表一个FackBook匿名用户，每条边表示两个用户之间的朋友关系；FackBook朋友圈社交网络为无权无向图；节点0, 107, 348, 414, 686, 698, 1684, 1912, 3437, 3980为spotlight nodes，也就是上面提到的10个朋友圈的主人，这也就意味着他们是图结构中的核心节点。引包import pandas as pd import numpy as np import networkx as nx import matplotlib.pyplot as plt from random import randint读取文件并创建图facebook_df = pd.read_csv( "./data/facebook_combined.txt.gz", compression="gzip", sep=" ", names=["start_node", "end_node"], ) facebook_df读取到的DataFrame共有88234行，每一行代表一个边（即一个连接关系）：根据DataFrame创建一个图结构：G = nx.from_pandas_edgelist(facebook_df, "start_node", "end_node")接下来就是进行图的可视化。因为默认情况下我们并不清楚图的数据结构是什么样的，因此选择随机布局（nx.random_layout(G)）来绘制图形：fig, ax = plt.subplots(figsize=(15, 9)) ax.axis("off") plot_options = {"node_size": 10, "with_labels": False, "width": 0.15} nx.draw_networkx(G, pos=nx.random_layout(G), ax=ax, **plot_options)显然，从上面的可视化结果中我们看不出来任何信息。因为默认的随机布局看起来杂乱无章。nx.spring_layout(G)可以根据节点和边的关系对图结构进行布局，但是考虑到本图结构中节点数众多，可以通过设置其参数iterations来避免算法迭代次数过多：pos = nx.spring_layout(G, iterations=15, seed=666) fig, ax = plt.subplots(figsize=(20, 20)) ax.axis("off") plot_options = {"node_size": 2, "with_labels": False, "width": 0.05} nx.draw_networkx(G, pos=pos, ax=ax, **plot_options)经过重新布局后，图中的结构似乎有那么一丝“拨云见日”了。①：图结构有聚类的性质，不同的节点分为几个大类；②：少数节点拥有较高的度，符合社交网络的度分布情况。图结构的基本拓扑属性查看节点数与边数：G.number_of_nodes(), G.number_of_edges()平均度（即图中每个节点平均有几个邻居）：np.mean([d for _, d in G.degree()])或者：sum(dict(G.degree()).values()) / len(G)图结构中有很多属性与路径有关，例如网络直径和平均最短路径。计算这两个属性值都需要涉及到计算图中所有节点对之间的最短路径，这一个非常消耗计算资源的步骤。shortest_path_lengths = dict(nx.all_pairs_shortest_path_length(G))shortest_path_lengths返回的结果是一个dict-of-dict格式的对象：shortest_path_lengths[42][0], shortest_path_lengths[1][32] # 节点42到节点0的最短路径、节点1到节点32的最短路径有上面离心率和直径的基本概念可以得之：$ \text{diameter}(G) = \max_{v \in G} \text{eccentricity}(v) $。因此，可以通过：diameter = max(nx.eccentricity(G, sp=shortest_path_lengths).values()) diameter求图中所有节点的离心率最大值来“曲线”求图的直径，通过参数sp=shortest_path_lengths来使用计算好的最短路径来求，避免了计算资源的消耗。同理，对于求平均最短路径而言，我们可以使用*nx.average_shortest_path_length()命令直接求，但是这样会造成计算资源的浪费。因此，仍可以通过上面计算的shortest_path_lengths来“曲线”求解：# Compute the average shortest path length for each node average_path_lengths = [ np.mean(list(spl.values())) for spl in shortest_path_lengths.values() ] # The average over all nodes np.mean(average_path_lengths)求得平均最短路径为3.6，这说明在社交网络中一份“信息”的平均传递人数是3.6。😂这说明，当你把一个“秘密”在互联网上告诉了另一个人，另一个人是有很大概率告诉另外一个人的。这也说明了“三人成虎”是有科学依据的。😂查看最短路径的直方图分布情况：plt.rcParams.update({ 'font.size': 20, 'font.family': ['Times New Roman', 'SimSun'] }) # We know the maximum shortest path length (the diameter), so create an array # to store values from 0 up to (and including) diameter path_lengths = np.zeros(diameter + 1, dtype=int) # Extract the frequency of shortest path lengths between two nodes for pls in shortest_path_lengths.values(): pl, cnts = np.unique(list(pls.values()), return_counts=True) path_lengths[pl] += cnts # Express frequency distribution as a percentage (ignoring path lengths of 0) freq_percent = 100 * path_lengths[1:] / path_lengths[1:].sum() # Plot the frequency distribution (ignoring path lengths of 0) as a percentage fig, ax = plt.subplots(figsize=(15, 8)) ax.bar(np.arange(1, diameter + 1), height=freq_percent) ax.set_title("Distribution of shortest path length in G") ax.set_xlabel("Shortest Path Length") ax.set_ylabel("Frequency (%)")查看网络密度：nx.density(G)网络密度为0.0108，远远小于1，说明社交网络也是一个稀疏网络。也很容易理解，中国有14亿人，你认识的有几个呢？获取连通组件（子图）的数量：nx.number_connected_components(G)返回结果为1。说明，整个社交网络都是连通的。FaceBook社交网络中心性分析度中心性采用度中心性进行计算：degree_centrality = nx.centrality.degree_centrality(G) # save results in a variable to use again, nx.degree_centrality(G)也可以 (sorted(degree_centrality.items(), key=lambda item: item[1], reverse=True))[:8]返回的结果中有一些有意思的信息：节点107、1684、1912、3437以及节点0都是spotlight nodes，因此其具有更好的中心性也无可厚非。但节点2543、2347、1888也就有非常高的度中心性，说明这三个人可能非常喜欢吃瓜🍉🍉🍉。获取度排名前8的选手：(sorted(G.degree, key=lambda item: item[1], reverse=True))[:8]获取网络中心性分布直方图：plt.figure(figsize=(15, 8)) plt.hist(degree_centrality.values(), bins=25) plt.xticks(ticks=[0, 0.025, 0.05, 0.1, 0.15, 0.2]) # set the x axis ticks plt.title("Degree Centrality Histogram ") plt.xlabel("Degree Centrality") plt.ylabel("Counts")以节点度为节点的大小和颜色区分，可视化网络结构：import matplotlib.pyplot as plt node_size = [ v**2 * 10000 for v in degree_centrality.values() ] plt.figure(figsize=(20, 20)) nx.draw_networkx(G, pos=pos, node_size=node_size, with_labels=False, width=0.15,node_color = node_size, cmap= plt.cm.cool) plt.axis("off")介数中心性对于一个节点的介数中心性，简单来讲就是一个网络中所有任意两节点之间最短路径中有几条经过了该节点。反映到社交网络中，介数中心性反映了该节点影响他人能力的大小。例如：当一个人有很高的介数中心性，说明社交网络中的很多信息都需要经过她/他来传播。计算介数中心性排名前八的用户：betweenness_centrality = nx.centrality.betweenness_centrality( G ) # save results in a variable to use again (sorted(betweenness_centrality.items(), key=lambda item: item[1], reverse=True))[:8]在介数中心性排名前八的用户中，节点0, 107, 1684, 1912, 3437的度中心性也在前八之中，说明这些节点不仅仅是spotlight nodes，其对于社交网络中信息的传播也十分重要。而对于节点567, 1085，其度中心性远远落后于前八的值：degree_centrality[567], degree_centrality[1085] # (0.015601783060921248, 0.016344725111441305)但是两者却拥有了很高的介数中心性，说明两者虽然在影响力上并不大，但是在信息传播过程中却非常重要（有点像吃瓜群众）。接近度中心性接近度中心性在控制社交网络中虚假信息的传播方面非常有研究价值。如果一个节点（用户）有很高的接近度中心性，那么说明其发布的信息更容易在短时间内就迅速影响到一大批人。计算接近度中心性排名前八的用户：closeness_centrality = nx.centrality.closeness_centrality( G ) # save results in a variable to use again (sorted(closeness_centrality.items(), key=lambda item: item[1], reverse=True))[:8]聚类系数查看图平均聚类系数：nx.clustering(G).values()所有节点聚类系数的直方图分布：plt.figure(figsize=(15, 8)) plt.hist(nx.clustering(G).values(), bins=50) plt.title("Clustering Coefficient Histogram ") plt.xlabel("Clustering Coefficient") plt.ylabel("Counts")桥检测nx.has_bridges(G)结果返回True，说明图结构中存在“易破坏的点”。bridges = list(nx.bridges(G)) print(len(bridges)) for i in bridges: print(i)结果返回75，说明图结构中存在75个桥。之所以存在这么多的桥是因为该社交网络仅使用了10个spotlight nodes，对于每个spotlight nodes的朋友节点的朋友圈并没有考虑。查看本地桥数量：local_bridges = list(nx.local_bridges(G, with_span=False)) len(local_bridges)可视化桥与本地桥：plt.figure(figsize=(20, 20)) nx.draw_networkx(G, pos=pos, node_size=10, with_labels=False, width=0.15) nx.draw_networkx_edges( G, pos, edgelist=local_bridges, width=0.5, edge_color="lawngreen" ) # green color for local bridges nx.draw_networkx_edges( G, pos, edgelist=bridges, width=0.5, edge_color="r" ) # red color for bridges plt.axis("off")度同配性计算社交网络的度同配性：nx.degree_assortativity_coefficient(G)或者：nx.degree_pearson_correlation_coefficient( G ) # use the potentially faster scipy.stats.pearsonr function.显示FB社交网络的度同配性为0.064，比较接近于0。说明：用户的社交关系在度数上既不显示同配也不显示异配的特征（non-assortative）。网络社区分类networkx 库提供了多种社区检测算法，其中 nx.community.label_propagation_communities(G) 和 nx.community.asyn_fluidc(G, k, seed=0) 是两种不同的社区发现算法。它们在社区检测的方法和理论基础上有所不同。标签传播算法（Label Propagation）:函数：nx.community.label_propagation_communities(G)原理：此算法基于网络中节点的标签传播。最初，每个节点被赋予一个唯一的标签，然后在迭代过程中，每个节点会采用其邻居中最常见的标签。最终，具有相同标签的节点被分为同一个社区。特点：这种算法简单且运行速度快，但它可能不稳定，即不同的运行可能产生不同的结果。此算法适用于大型网络，但可能不适合发现重叠社区或具有层次结构的社区。异步流体社区检测（Asynchronous Fluid Communities）:函数：nx.community.asyn_fluidc(G, k, seed=0)原理：这是一种基于流体动力学的社区检测方法。算法试图通过模拟 k 个流体的扩散来确定社区。每个流体源（源自一个节点）扩散到网络，逐渐形成社区，直到达到稳定状态。特点：这种算法适用于中等大小的网络，并且需要预先指定社区数量（k）。它通常能够找到比较均衡的社区结构，但对于非常大的网络，算法的效率和效果可能受限。两者的不同之处：算法原理：标签传播算法依赖于标签的传播和统一，而异步流体社区检测基于流体动力学的模型。社区数量：标签传播算法不需要预先指定社区的数量，它会自然形成，而异步流体社区检测需要预先指定社区数量 k。适用性：标签传播算法适用于大型网络且运行快速，但可能不够稳定；异步流体社区检测适用于中等规模的网络，需要预设社区数量，可能对社区大小更加均衡。标签传播算法：colors = ["" for x in range(G.number_of_nodes())] # initialize colors list counter = 0 for com in nx.community.label_propagation_communities(G): color = "#%06X" % randint(0, 0xFFFFFF) # creates random RGB color counter += 1 for node in list( com ): # fill colors list with the particular color for the community nodes colors[node] = color plt.figure(figsize=(15, 9)) plt.axis("off") nx.draw_networkx( G, pos=pos, node_size=10, with_labels=False, width=0.15, node_color=colors ) 异步流体社区检测可视化：colors = ["" for x in range(G.number_of_nodes())] for com in nx.community.asyn_fluidc(G, 8, seed=0): color = "#%06X" % randint(0, 0xFFFFFF) # creates random RGB color for node in list(com): colors[node] = color plt.figure(figsize=(15, 9)) plt.axis("off") nx.draw_networkx( G, pos=pos, node_size=10, with_labels=False, width=0.15, node_color=colors )
定义什么是欧拉公式？$$ e^{i\pi}+1=0 $$这个公式将：$e$：自然对数的底；$i$：虚数的单位；$\pi$：圆周率结合到了一起，优美巧妙，因此也被称为“上帝公式”。Quote / 参考维基百科：欧拉公式。复数i与自然常数e在推导与证明欧拉公式之前，有必要对其中两个重要组成部分：复数——$i$；自然常数——$e$进行复习。复数i复数$i$的定义其实很简单，就是：$$ i^2=-1 $$最初提出复数的那一帮人，就是因为解高次方程的时候会遇到这样的情况：$$ x^2+1=0 $$虽然在我们对于自然数的理解下，这样求解出来的数没有任何实际意义，但是其在数学中却是真实存在的。数学家们之间的解方程大赛对于复数这个在自然界不存在，但是在数学中却是存在的东西，笛卡尔将其命名为Imaginary number——想象中的数字。意义所以，复数在数学中的意义到底应该是什么呢？仅仅就是为了解方程吗🤔？我们又是如何来表达它呢？接下来，让我们通过一个小例子来认识下它：从自然数1开始，将其与$i$相乘，得到了：$$ 1×i=i $$在此基础上，然后再乘一个$i$：$$ 1×i×i=-1 $$可以看到，自然数1再经过两次乘$i$后，完成了反向。如果再乘上两个：$$ 1×i×i×i×i=1 $$之后会发现，1竟然又回到了最初的起点。看到了这里，你的心里可能已经猜到了😆——这个过程，像极了旋转：所以，这个时候我们很容易会发现复数虚部所在的轴与实部所在的轴构成了一个平面。在这个平面中，一个复数的实部为x轴，虚部$i$为y轴，两者张成了一个二维复平面。对于一个复数$4+3i$，其在复平面上的表达类似于一个向量：运算加法运算在复平面中，复数的运算法则遵循矢量运算法则（平行四边形法则，实部与实部相加，虚部与虚部相加），例如：$$ (-1+2i)+(2-1i)=1+i $$乘法运算令：$$ y=\begin{cases} 1, n=0\\ \left( 1+i \right) ^n, n>0\\ \end{cases} $$计算得：$$ y_{n=0}=1 \\ y_{n=1}=1+i \\ y_{n=2}=\left( 1+i \right) ^2=1+2i+\left( i \right) ^2=2i \\ y_{n=3}=\left( 1+i \right) ^3=2i\times \left( 1+i \right) =-2+2i \\ y_{n=4}=\left( 1+i \right) ^4=2i\times 2i=-4 \\ ... $$依次绘制出$y_{n=0}、y_{n=1}、y_{n=2}、...$，会发现：$y_{m-1}$乘上$1+i$后变成了$y_{m}$。相比之下，$y_{m}$由$y_{m-1}$所代表的向量旋转45°，且模长变为原来得$\sqrt{2}$倍得来。其中45°是$1+i$与x轴的夹角，$\sqrt{2}$是$1+i$在复平面的模长。由此，我们可以进行一个推广：在复平面上，一个复数$b=m+ni$乘上任意一个复数$a$，其几何含义为：让这个复数$a$所代表的向量进行旋转并缩放，旋转的角度（复数$b$的辐角，即$m+ni$与x轴的夹角）为$arctan\frac{n}{m}$，缩放倍数（复数$b$的模长）为$\sqrt{m^2+n^2}$。一个有趣的现象，如果我们在复平面上绘制出$(1+i)^n$所代表的点：参考Python绘图代码def get_value(n): x = (1 + 1j) ** n return x.real, x.imag def plot_arrow(x, y, ax=plt, lines=True): if lines: ax.vlines(x, 0, y, 'r', '--') ax.hlines(y, 0, x, 'r', '--') ax.arrow(0, 0, x, y, color='cyan', head_width=0.2, fc='y', length_includes_head=True ) def plot_text(x, y, ax=plt): valign, halign = None, None if y > 0: text_y = f'{y}i' if x != 0: text_y = '+' + text_y y += .2 valign = 'bottom' elif y == 0: text_y = '' y += .2 valign = 'bottom' else: text_y = f'{y}i' y -= .2 valign = 'top' if x > 0: text_x = f'{x}' x += .2 halign = 'left' elif x == 0: text_x = '' x += .2 halign = 'left' else: text_x = f'{x}' x -= .2 halign = 'right' ax.text(x=x, y=y, #文本x、y轴坐标 s=f'${text_x}{text_y}$', #文本内容 fontdict=dict(fontsize=12), #字体属性字典 verticalalignment=valign, horizontalalignment=halign, bbox={ #添加文字背景色 'facecolor': 'white', #填充色 'edgecolor': 'gray', #外框色 'alpha': 0.3, #框透明度 'pad': .1, #本文与框周围距离 'boxstyle': 'round' } ) # 绘图全局设置 plt.rcParams.update({ # 'font.size': 15, 'font.family': ['Times New Roman', 'SimSun'] }) limit_coor_num = 200 plt.hlines(0, -limit_coor_num, limit_coor_num, '#ccc', '--') plt.vlines(0, -limit_coor_num, limit_coor_num, '#ccc', '--') plt.xlim([-limit_coor_num, limit_coor_num]) plt.ylim([-limit_coor_num, limit_coor_num]) plt.yticks( [i for i in range(-limit_coor_num, limit_coor_num + 1) if i % 40 == 0], [f'{i}i' for i in range(-limit_coor_num, limit_coor_num + 1) if i % 40 == 0] ) plt.xlabel('Real') plt.ylabel('Imaginary') plt.title('Complex Plane') # plt.grid() # point_list = [(1, 0), (0, 1), (-1, 0), (0, -1)] # point_list = [(-1, 2), (2, -1), (1, 1)] # for x, y in point_list: # plot_arrow(x, y) # plot_text(x, y) # for i in range(8): # x, y = get_value(i) # plot_arrow(x, y, lines=False) # plot_text(x, y) points_plt = [get_value(i) for i in np.linspace(0, 100, 1000)] x, y = zip(*points_plt) plt.scatter(x, y) # plot_arrow(0, -2, lines=False) # plot_text(0, -2) plt.show()Quote / 参考Bilibili @隐姓埋名的刘老：复数，一个充满故事又奇特的数。自然常数e在数学中$e$被称为自然常数（Natural Constant）。$$ e=2.71828... $$可是，自然数有这么多，为什么偏偏$e$会被称为自然常数？先来看一个简单的案例，在经济学中，有一个概念叫做复利运算。假设银行一年的年利率为100%，你在银行存了1块钱，那么一年之后你的存款应该会变为：$$ 1×(1+100\%)=2 $$但其实银行的利率往往还可以进行拆分，比如你把1块钱存在银行半年，这个时候利率就变成了50%，在半年后你拿到的钱就是：$$ 1×(1+50\%)=1.5 $$然后你可以将你的一块五再存银行半年，再经过半年后，你的钱就会变成：$$ 1×(1+50\%)×(1+50\%)=1.5×(1+50\%)=2.25 $$这个时候，你会惊奇地发现，一年分两次存钱比一次存钱要多拿到$0.25$！！！那是不是意味着我分的次数越多就拿到的越多呢？假设我们分为n期，那经过n期后我们的利率就会变为：$$ (1+\frac{1}{n})^n $$问题就是，随着n的增大，上面的利率真的会无限增大吗？在高等数学中，有这么一个极限定义：$$ \underset{n\rightarrow \infty}{\lim}\left( 1+\frac{1}{n} \right) ^n=e $$也就是说，随着$n$的无限增大，$(1+\frac{1}{n})^n$并不是也可以无穷增大的，而是趋近于一个值$2.718281828459045...$。而这个值，就是自然常数$e$的值。上面这个分为n期的假设，其实非常符合自然界中万事万物发展的规律，一切事物的变化都不是突然来的，都是经过n期的逐渐积累，然后引发了质变。这也就是其名字“自然常数”的来历，因此其定义像极了世间万物发展的模式。Quote / 参考Bilibili @隐姓埋名的刘老：自然常数e这个数，怎么就自然了？；Bilibili @隐姓埋名的刘老：用几何直觉理解欧拉公式！【中学生也能懂|manim】。欧拉公式思想的转变由上面对于自然常数$e$的讲解可知：$$ \underset{n\rightarrow \infty}{\lim}\left( 1+\frac{1}{n} \right) ^n=e $$$e$的定义并不是一个“实实在在”存在的数字，那么对于$e^{i\pi}$，其也并不意味着$i\pi$个$e$相乘！！！那么，$e^{i\pi}$到底代表什么呢？对于上面提到的等式：$$ \underset{n\rightarrow \infty}{\lim}\left( 1+\frac{1}{n} \right) ^n=e $$还可以进行进一步的推广：$$ \underset{n\rightarrow \infty}{\lim}\left( 1+\frac{x}{n} \right) ^n=e^x $$那么欧拉公式就可以等价变换为：$$ e^{i\pi} = \underset{n\rightarrow \infty}{\lim}\left( 1+\frac{i\pi}{n} \right) ^n $$上面这个等式，才是$e^{i\pi}$的真正含义！！！可视化求解定义求极限的函数：def get_value(n, p): x = np.power((1 + 3.14j / n), p) return x.real, x.imag绘图：plt.xlim([-3, 3]) plt.ylim([0, 4]) plt.yticks( [i for i in range(0, 5) if i % 2 == 0], [f'{i}i' for i in range(0, 5) if i % 2 == 0] ) for i in [1, 2, 5, 10, 20, 10000]: point_coor = [get_value(i, j) for j in range(1, i + 1)] x, y = zip(*point_coor) plt.scatter(x, y, s=(1 / i) * 100, marker='*') plt.plot(x, y, label=f'n={i}') # x、y轴等比例 ax = plt.gca() ax.set_aspect(1) plt.xlabel('Real') plt.ylabel('Imaginary') plt.title(r'$e^{i\pi} = \underset{n\rightarrow \infty}{\lim}\left( 1+\frac{i\pi}{n} \right) ^n=-1$') plt.legend(ncol=3, loc='lower center', bbox_to_anchor=(0.5, -0.4)) plt.show()由上面可视化的结果可以看出，对于：$$ e^{i\pi} = \underset{n\rightarrow \infty}{\lim}\left( 1+\frac{i\pi}{n} \right) ^n $$当$n\rightarrow \infty$时，$\left( 1+\frac{i\pi}{n} \right)^n$趋近于$-1$。那么对于更一般的情况：$$ e^{i\theta} = \underset{n\rightarrow \infty}{\lim}\left( 1+\frac{i\theta}{n} \right) ^n $$应该为多少呢？首先将求解极限的函数改写为：def get_value(n, p, theta): x = np.power((1 + theta * 1j / n), p) return x.real, x.imag以便用于支持任意角度的输入。接下来以$\frac{\pi}{2}、\frac{2\pi}{3}、\frac{5\pi}{4}、\frac{7\pi}{4}$为例，可视化其结果：fig, ax_arr = plt.subplots(2, 2) fig.subplots_adjust(hspace=0.5) theta_list = [1 / 2, 2 / 3, 5 / 4, 7 / 4] sub_title = [r'\frac{1}{2}', r'\frac{2}{3}', r'\frac{5}{4}', r'\frac{7}{4}'] for i in range(ax_arr.shape[0]): for j in range(ax_arr.shape[1]): cur_ax = ax_arr[i, j] for k in [1, 2, 5, 10, 20, 10000]: point_coor = [get_value(k, v, theta_list[i * 2 + j] * np.pi) for v in range(1, k + 1)] x, y = zip(*point_coor) cur_ax.scatter(x, y, s=(1 / k) * 100, marker='*') cur_ax.plot(x, y, label=f'n={i}') cur_ax.axvline(0, color='#ccc', linestyle='--') cur_ax.axhline(0, color='#ccc', linestyle='--') # x、y轴等比例 ax = plt.gca() ax.set_aspect(1) if j == 0: cur_ax.set_ylabel('Imaginary') if i == 1: cur_ax.set_xlabel('Real') cur_ax.set_title(f'$' + sub_title[i * 2 + j] + '\pi$') lines, labels = fig.axes[-1].get_legend_handles_labels() # 获取最后一个子图的图例 fig.legend(lines, labels, ncol=3, loc='upper center', bbox_to_anchor=(0.5, 0.03)) plt.tight_layout() plt.show()当$n\rightarrow \infty$时，$\left( 1+\frac{i\theta}{n} \right)^n$无限趋近于点$(\cos\theta, \sin\theta)$，也就是说：$$ e^{i\theta} = \underset{n\rightarrow \infty}{\lim}\left( 1+\frac{i\theta}{n} \right) ^n = \cos\theta+i\sin\theta $$这，便是宇宙间最优美的公式——欧拉公式。
Introduction / 引言大学期间用来打发无聊时间学的Python没想到竟然在写毕业论文的时候用处这么大，整个硕士论文所做研究，从前期的数据整理、数据分析，到最后的数据可视化我基本上都使用Python来完成，这篇博客就来分享下我毕业论文课题中所做数据分析相关的Python代码。本博文所有相关的代码都上传在GitHub仓库：Data-Analysis-for-Thesis，如果帮到你了，记得给我来个Star☺️，也可以顺便去参观下我的个人博客亚灿网志。Python环境配置环境安装首先是在Python官网下载你计算机对应的Python软件，然后安装。安装过程基本都是傻瓜式，不做过多叙述，一路回车即可。Python官网IDE选择然后就是IDE（Integrated Development Environment）的配置，最常见的肯定是VSCode，特点是足够轻量化且免费。但是我选择的是PyCharm，我选择它的理由是以前学Python的时候就用的它，因此快捷键什么的都比较熟悉。PyCharm正版是收费的，教育版可以申请一年的使用权限，破解版的自行百度即可。VSCodePyCharmWarning / 注意Pycharm也是可以申请教育版嘀，但是教育版并不支持编写Jupyter Notebook的.ipynb文件。Tips / 提示关于IDE的选择大家没必要纠结，新手不太了解想简单上手的话就选择VSCode；用过PyCharm并且熟练怎么破解安装的那就选择PyCharm。包的安装数据矩阵分析及处理：Pandas、Numpy、Math、Scipy；绘图可视化：Matplotlib、Seaborn；其他包：hues可以在控制台打印出彩色的提示信息，用法也比较简单，官方手册。hues的简单用法数据整理及预处理数据介绍论文所做内容简单来说就是对不同样本点的多个指标进行相关性分析，具体的数据格式就是：样本点变量1变量2变量nD-H1122424...L-N2324232...............基本常量的定义IMG_TYPE = '.svg' # 出图格式，我选择矢量图svg DATA_PATH, EXPORT_PATH, DPI = './data/', './export/', 300 # 数据存储路径、出图写入路径、出图DPI ROW_NUM, COL_NUM = 5, 3 # 对于多子图的图片，定义默认布图网格为5×3文件读取读取csv文件需要使用pandas的pd.read_csv()方法，具体的参数有：index_col：设置行索引为哪一列，可以使用序号或者列名称；sep：csv文件中的分隔符，默认常见的用法都可以自动识别，不需要设置；header：设置表头，参数为None就是没有表头，设置为n就是把第n行读取为表头；names：设置列名称，参数为list；usecols：仅读取文件内某几列。Quote / 参考具体用法可以参考李庆辉所著《深入浅出Pandas——利用Python进行数据处理与分析》3.2章 读取CSV（PDF P89）。数据表合并首先遇到的第一个需求就是，所有样本点的列变量存储在不同的数据表中，比如，样本点的指标分为上覆水的指标与沉积物的指标两部分，分别存储在两个或者多个数据表中，那么如何将两个或者多个数据表进行合并呢？all_df_0 = pd.concat([gene_abundance, sediment_env_df, eems_df, overlying_water_env_df, sediment_size_df], axis=1) # 将五个数据表按照行索引合并Quote / 参考具体用法可以参考李庆辉所著《深入浅出Pandas——利用Python进行数据处理与分析》7.2章 数据连接pd.concat（PDF P274）。根据行索引为每个样本点设置分类行索引也就是每个样本点的标记名，分别为“D-H1”、“L-N3”之类的，其中第一个字符的值为“D”、“L”、“W”，分别代表枯水期（Dry Season）、平水期（Level Season）、丰水期（Wet Season）；第二个字符为一个分隔符“-”；第三个字符的值为“N”、“P”、“H”，分别代表三条不同的河流，南淝河、派河、杭埠河；最后的数字代表是第几个样本点。例如：“D-H1”代表枯水期杭埠河第一个样本点、“L-N3”代表平水期南淝河第三个样本点。那么问题来了，我想要为合并后的数据表新增两列“River”、“Period”，分别来反应这个样本点的属性，应该如何实现呢？样本点变量1变量2变量nRiverPeriodD-H1122424...Hangbu RiverDry SeasonL-N2324232...Nanfei RiverLevel Season..................思路其实也很简单，就是使用apply函数分别对每一行（也就是每一个样本点）进行处理，获取该行的行索引，然后对行索引的字符进行判断即可：all_df['Period'] = all_df.apply(lambda x: 'Dry' if 'D' in x.name else ('Wet' if 'W' in x.name else 'Level'), axis=1) + ' Season' all_df['River'] = all_df.apply(lambda x: 'Nanfei' if 'N' in x.name else ('Pai' if 'P' in x.name else 'Hangbu'), axis=1) + ' River'Tips / 提示这里使用了Python列表推导式相关的知识，具体的讲解可以看之前的博文：「Python」列表推导式。在对每一行的样本点添加River、Period变量后，会有一个问题，River、Period两列的数据都是Object字符串类型。这种数据类型有两个问题：如果数据矩阵有几十万行，那么这两列会占用很大的内存空间；对数据进行绘图过程中，我想把River变量按照Nanfei River、Pai River、Hangbu River的顺序排列，那么就很麻烦。因为字符串变量默认是按照首字母的顺序来进行排序的，默认排序是Hangbu River、Nanfei River、Pai River。为了解决这两个问题，我们可以将这两列的数据由原来的object类型转换为Category类型，Category的好处就是，当数据量较大时，可以显著减小数据所占用的内存；第二还可以对数据类型进行排序。具体的处理方法：# 对季节、河流两列进行排序，首先定义category类型顺序 river_order = CategoricalDtype( # 河流的顺序定义为南淝河、派河、杭埠河 ['Nanfei River', 'Pai River', 'Hangbu River'], ordered=True ) period_order = CategoricalDtype( # 时期的顺序定义为枯水期、平水期、丰水期 ['Dry Season', 'Level Season', 'Wet Season'], ordered=True ) # 将两列Object类型数据转换为category类型并排序 all_df['River'] = all_df['River'].astype(river_order) all_df['Period'] = all_df['Period'].astype(period_order) all_df = all_df.sort_values(by=['Period', 'River', 'SortNum'])Warning / 注意在将一列数据转化为Category对象后，如果数据表中没有某个Category，但是绘图的时候还是会占用一个位置，下面举例说明：这个数据表中的Period列已经不包含Level Season的数据，但是使用.value_counts()方法统计的时候还是会对其进行统计。如果这个时候进行绘图，会出现一个问题：中间数据明明是空的，但是还是会占位。解决这个问题的方法是使用：nanfei_dry_and_wet_gene['Period'].cat.remove_unused_categories()删除Period列中没有用到的Category，这个时候再次进行绘图就可以正常显示了。异常值处理缺失值的填充Pandas中缺失值的填充所用方法时pd.fillna()，具体的参数可以填写：In [16]: pd.DataFrame.fillna Out[16]: <function pandas.core.frame.DataFrame.fillna(self, value: 'object | ArrayLike | None' = None, method: 'FillnaOptions | None' = None, axis: 'Axis | None' = None, inplace: 'bool' = False, limit=None, downcast=None) -> 'DataFrame | None'>value：直接将缺失值填充为字符串或者数字；method：填充方式，method='ffill' 向前填充，method='bfill'向后填充，也就是说用前面的值来填充NA或用后面的值来填充NA。另外，在使用读取pd.read_csv()读取csv文件的时候，也可以通过参数：na_values=Nonekeep_default_na=Truena_filter=True的设置来对NA值进行过滤或者识别。删除缺失值使用pd.DataFrame.dropna()方法完成缺失值的删除：In [17]: pd.DataFrame.dropna Out[17]: <function pandas.core.frame.DataFrame.dropna(self, axis: 'Axis' = 0, how: 'str' = 'any', thresh=None, subset: 'IndexLabel' = None, inplace: 'bool' = False)>通过参数how的属性值来设置：any：当每一行有一个缺失值时就删除这一行；all：当一行所有的数据都时缺失值时再删除这一行。重复值的删除使用pd.DataFrame.drop_duplicates()方法完成缺失值的删除：In [18]: pd.DataFrame.drop_duplicates Out[18]: <function pandas.core.frame.DataFrame.drop_duplicates(self, subset: 'Hashable | Sequence[Hashable] | None' = None, keep: "Literal['first'] | Literal['last'] | Literal[False]" = 'first', inplace: 'bool' = False, ignore_index: 'bool' = False) -> 'DataFrame | None'>通过参数keep的属性值来设置：first：所有重复行删除，保留第一行；last：所有重复行删除，保留最后一行。数据处理与可视化绘图前的小准备画图格式的定义如何在Matplotlib中显示中文：plt.rcParams['font.sans-serif'] = ['SimHei'] # 用于显示中文 plt.rcParams['axes.unicode_minus'] = False如何定义Seaborn绘图的默认字体属性：sns.set( style='darkgrid', # 绘图风格 font='Times New Roman', # 默认字体 font_scale=3 # 默认字体比例 )如何实现子图编号CHAR = [chr(i) for i in range(97, 123)] # 获取26个英文字母，用于给子图编号定义一个26个英文字母的list，循环绘制子图的时候直接调用即可。重复代码的打包每次进行数据分析我都会新建一个.ipynb文件，而数据分析前都需要经过数据表合并、数据清洗等工作，那么最好的方式其实是将数据分析前的准备工作进行一个打包，然后在.ipynb文件的第一行引入包即可。例如：我新建一个ResearchMain.py文件，然后将所有数据表合并、数据清洗的代码都放在这个文件里：# 引入数据分析常用的包 ... # 读取文件 ... # 合并文件 ... # 数据表清洗 ...然后在每次新建.ipynb文件进行数据分析时，我都会在第一行使用：from ResearchMain import *来引入所有ResearchMain.py文件中定义的变量与方法。QPCR标准曲线可视化Tips / 提示什么是标准曲线？简单来说，自变量x与因变量y之间存在某种线性关系——$y=ax+b$，那么我们可以通过多次改变自变量x的值，然后观察y的值并记录，得到几组对应的$x_1$、$x_2$、$x_3$、$x_4$、$x_5$、...与$y_1$、$y_2$、$y_3$、$y_4$、$y_5$、...，那么我们就可以通过得到的这几组数据来对自变量x与因变量y进行线性拟合，从而得到一个标准曲线$y=ax+b$，有了标准曲线之后，我们就可以直接输入任意的自变量x值，计算出因变量y的值。在Numpy中，拟合函数直接有现成的，可以直接调用：In [37]: x, y = [1, 2, 3, 4, 5], [2, 4, 6, 8, 11] # 需要进行拟合的自变量与因变量 In [38]: np.polyfit(x, y, 1) # 对自变量x与因变量y进行拟合，且拟合为一次函数 Out[38]: array([ 2.2, -0.4]) # 拟合结果为y=2.2x-0.4拟合完毕如何拼接拟合方程？In [25]: fitting_arr = np.polyfit(x, y, 1) In [26]: fitting_equ = np.poly1d(fitting_arr) # 获取拟合方程 In [27]: fitting_equ([6, 7, 8, 9, 10]) # 利用拟合方程计算任意自变量对应的因变量 Out[27]: array([12.8, 15. , 17.2, 19.4, 21.6]) In [28]: fitting_equ(11) Out[28]: 23.799999999999997 # 23.799999999999997=11*2.2-0.4获取拟合的$R^2$：In [43]: np.corrcoef(y, fitting_equ(x))[0, 1] ** 2 Out[43]: 0.9918032786885246$R^2$的计算原理其实就是把自变量带入拟合方程，将计算出来的因变量与原始的因变量进行比较，计算其相关性。搞懂上面的原理之后，进行标准曲线的可视化其实就很简单了。图中可以看出，还生成了一个拼接的一元一次方程，方程的拼接可以直接用我写好的函数，函数的具体用法以及讲解已经在注释里说的很清楚了：Tips / 提示函数的主要作用就是传入np.polyfit(X, Y, DEG)返回的list对象，比如返回的是[2, 3, -4]，那么就输出2x^2+3x-4，可以自动识别n元函数。def make_fit_equ_str(paras): """ 拼接拟合方程 :param paras: 传入参数为np.polyfit(X, Y, DEG)返回的对象 :return: 返回一个拼接拟合方程的字符串(LaTex格式) """ fit_equ_str = '' for i in range(len(paras)): # 遍历常数项的过程中，主要需要进行两步“组装”：①如果该常数项不是第一个常数项，且该常数项大于0，需要转化为字符串并在前面添加一个“+”； if paras[i]: # 如果常数项不为0(为0直接跳过这一项) cur_item = str(round(paras[i], 4)) # 每个常数项都保留两位小数，并转化为字符串 # 如果不是第一个常数项，且该常数项大于0，需要转化为字符串前面添加一个“+” if paras[i] > 0 and fit_equ_str != '': # 注意两种特殊情况：[0, 0, -5, 7]、[0, -5, 7]如何排除 cur_item = '+' + cur_item # ②为每一个常数项添加x的幂指数 if i == (len(paras) - 1) and cur_item != '0': # 最后一个项只有一个常数，不需要添加x了（如果为0就不用+了） fit_equ_str += cur_item else: # 前面的每一项都需要在最后加上一个x^幂数 idempotent_num = len(paras) - 1 - i # 对应的幂指数 if idempotent_num == 1: # 拟合为一次多项式 fit_equ_str += cur_item + 'x' else: # 拟合为二次及二次以上的多项式 fit_equ_str += cur_item + 'x^' + str(idempotent_num) #方程拼接 # fit_equ_str = '$' + fit_equ_str + '$' # 将多项式转化为LaTex公式 return fit_equ_strQPCR数据处理上一步的操作已经制作出来的功能基因的标准曲线，接下来就是把仪器测得的CT值（自变量）转化为丰度值（因变量），就是带进拟合的标准曲线计算下就OK了。具体绘图代码直接看GitHub代码即可，没有什么难度。三维荧光数据可视化读取数据表使用日立F-7000荧光光谱仪对沉积物中溶解性有机质（Dissolved Organic Matter, DOM）结构特征和组成成分进行表征。仪器得到的数据是.txt格式，且有用的数据表是从Data Points这一行后面开始的。得到的txt文件所以说，我们要先读取.txt文件，循环读取每一行，直到读取到Data Points这一行，说明已经到数据表了。def get_skip_rows(path): """ 读取txt文件，并在文件中查找含有'Data Points'的行，数据矩阵就在这一行的下面 :param path: 文件路径 :return: 数据矩阵开始的行号 """ f = open(path) for index, line in enumerate(f.readlines()): if 'Data Points' in line: return index + 1这个函数的作用就是输入.txt文件的路径，然后会返回需要的数据表是在第几行开始的。然后使用pd.read_table()方法读取.txt，并通过设置skiprows的值，来跳过前面无用的数据。例如利用get_skip_rows()函数获取到.txt文件中数据表从第156行开始：df = pd.read_table(search_info['Path'], skiprows=156, index_col=0)这样的话就可以完美跳过.txt文件前面无用的数据，直接读取所需的数据表。消除瑞利散射瑞利散射消除前后对比瑞利散射的消除其实很简单，观察数据表就可以看出来，瑞利散射其实就是不该出现在某个区间内数据峰，我们只需要慢慢根据$E_x$与$E_m$的设置范围来进行消除就行了：for i in range(modify_df.shape[0]): # 遍历ex for j in range(modify_df.shape[1]): # 遍历em ex = 200 + i * 5 # ex范围为200~450，间隔为5 em = 280 + j * 5 # em范围为280~550，间隔为5 if ex + 455 > em > ex + 145: modify_df.iloc[i, j] = 0 # 不能使用0, 因为0是有意义的数据光谱图的分区直接使用：cur_ax.axvline(250, color='white', linestyle='--', linewidth=5) # 垂直线 cur_ax.axhline(330, color='white', linestyle='--', linewidth=5) # 水平线就可以在图上画线了。需要讲解的其实就这几个部分，别的直接看代码就行了。粒径数据可视化沉积物粒径百分比分布使用Malvern Mastersizer 2000型激光粒度仪进行分析。得到数据后需要手动整理为.csv格式。数据格式首先来讲解下数据格式，每一列代表一个样本，每一行代表对应粒径所占百分比。例如图中红方框所示就是代表D-N4样本点对应粒径为0.955 μm颗粒占比为0.03%。清除空行上图中可以看出，数据有很多空行，那么首先第一步就是清除掉这些空行：df.dropna(how='all', inplace=True) # 删除缺失值（行全为空） df.dropna(how='all', axis=1, inplace=True) # 删除缺失值（列全为空）绘图还是直接看代码吧。常规指标数据可视化# 子图摆放设置 ROW_NUM = 2 COL_NUM = int(len(df.columns) / ROW_NUM) # 绘图设置 sns.set( style='darkgrid', font='Times New Roman', font_scale=2 ) # 布局设置 plt.figure(dpi=DPI) fig, ax_arr = plt.subplots(ROW_NUM, COL_NUM, sharex='col', figsize=(20, 10)) fig.subplots_adjust(hspace=0.2, wspace=0.2) # 循环绘图 for index, s in enumerate(df.iteritems()): cur_df = pd.concat([s[1], affix_info], axis=1) # 计算行列数 row_num = index // COL_NUM col_num = index - row_num * COL_NUM cur_ax = ax_arr[row_num][col_num] cur_plot = sns.boxplot( # 绘图 x="Period", # X轴 y=s[0], # Y轴 hue="River", # 颜色分类 data=cur_df, # 数据表 dodge=True, palette="Set1", # 配色 ax=cur_ax # 绘图坐标轴 ) # 子图图名、x轴、y轴、图例 cur_ax.tick_params(axis='x', rotation=90) if len(df.columns) - (index + 1) >= COL_NUM: cur_plot.set_xlabel(None) cur_plot.set_title(f'({CHAR[index]}) ${s[0][2:]}$ ({unit[s[0]]})',fontproperties="Times New Roman", fontsize=20) cur_ax.get_legend().remove() # 移除子图图例 cur_plot.set_ylabel(None) # 子图的图例相同，获取最后一个子图的图例 lines, labels = fig.axes[-1].get_legend_handles_labels() fig.legend(lines, labels, ncol=3, loc='lower center', bbox_to_anchor=(0.5, -0.16)) # 保存图片 # fig.savefig(EXPORT_PATH + '沉积物环境因子时空变化' + IMG_TYPE, dpi=DPI, bbox_inches='tight' )数据分析主成分分析主成分分析（Prin­ci­pal Com­po­nent Analy­sis, PCA）的具体逻辑及Python实现方法可以看：PCA主成分析原理、理解和代码实现。聚类分析K-means聚类与层次聚类相关代码直接查看GitHub项目Data-Analysis-for-Thesis。Mantel Test见博文：Mantel Test算法原理讲解、代码实现、绘图可视化。
? PREPOSITIONS - at on in/since for/during whileLinks / 链接Video course link: https://www.youtube.com/watch?v=D5ucI2Tw7H0at/on/inat: Specific.on: 1 day, 1 date.in: month, year, season, century.Example / 举例We could meet at 3:00.See you on Friday.I will be graduating from this university in May.since/forsince: + point(9:00, March, 1995)for: + period(5 minutes, 5 decades)Example / 举例She's been working here since 2018.I've lived here for two years.during/whileduring: + noun(this week, this meeting)while: + verb(watching TV, I was watching TV)Example / 举例I was busy during the day.The phone rang while we were sleeping.12 Confusing English VerbsLinks / 链接Video course link: https://www.youtube.com/watch?v=TIUuTNoLPWEbring/takeBring this over here and taking this over there.Example / 举例Could you please bring this gift to our teacher for me?please take this over there.hear/listenWhen we hear something, it just means some auditory information is coming through our ears. listen means we are paying attention to something. When you use "listen", you should add the "to" end of the "listen".Example / 举例Did you hear that sound?Please, could you be quiet? I's just trying to listen to the news.remember/remindTo remember just means don't forget. But when you remind, you tell someone else or you help someone else not to forget.Example / 举例Remember to drink eight glasses of water a dayPlease remind me to call John.borrow/lendLend or Borrow?Example / 举例I borrow some books from the library.The library lends me textbooks.coming/goingExample / 举例Come here! I want to show you something.When are you going there?spent/wastedWhen we use "spent", it just means the feeling is more neutral or positive. But "wasted" means some thing is useless, it's negative.Example / 举例I spent lots of money on my car, because I really love it.Don't waste your money on that movie. It's really bad!Talking about FAMILYLinks / 链接Video course link: https://www.youtube.com/watch?v=O3h7bUFrP0QThe favorite/favourite is the child in the family or the person in the family that the parents like the best.A daddy's girl is a girl or woman who has a very close relationship with their father.The black sheep of the family is the person in the family who has embarrassed the family, or who is a troublemaker in the family.Sam gets along with his sister. This means Sam and his sister have a good relationship.The breadwinner is The person in the family who earns the most money."Diabetes is a disease that runs in the family." This means Many people in the family have diabetes.A stay-at-home parent is someone who stays with the kids, takes care of the kids, and doesn't work full-time outside the house.A mama's boy is a boy or man who has an unhealthy relationship with his mom. Specifically, he is too dependent on his mother for everything. Some people use this as an insult.Take after means to be similar to an older member of your family in appearance or character.The baby of the family is the youngest in the family.
AdviceWhat three pieces of advice will you give your children?Answer / 回答Actually, I'm not finishing my studies at this point. So I never thought about this question. Okay, Let me think... Firstly, I think my children should know that being a good man is more important than studying. Secondly, everything is easy said than done. You must become a brave person and step out of your comfort zone valiantly. Last but not least, Happiness is more important than money. Spend more time on things you enjoy. That's all, Thank you.Whose advice do you follow more, your parents' or your friends' advice?Answer / 回答Okay. About this question, I think it depends on different situations. Parents undoubtedly have more experience than us. But in some situations, parents aren't very clear about contemporary environment. So I think my friends will be more useful in these situations.If you could give the leader of your country some advice, what would it be?Answer / 回答In my daily life, I don't care about politics. So I think I will choose to say some questions about my major. You know my major is water environment science. I main job is to research how to make wastewater cleaner and more efficient. I hope the leader in my country should pay more attention to the water environment. Everyone should drink clean water every day. Water safety and water quality are important to our bodies.If you could go back in time and give yourself some advice, what would you tell your younger self?Answer / 回答This is an interesting question. If I could go back to my junior high school. I want to tell my younger self two things. First, pay more attention to your studies. Second, get more experience from your parents.Who was the wisest person in your country's history? What kind of advice did they give?Answer / 回答In my opinion, the wisest person in China's history is the leaders of the Communist Party. During World War II, these leaders guided the Chinese to defeat the invasion of Japan. Today, China is the world's second largest economy and everyone can have a better life. We should all appreciate the elder generation.Where do you go to get good advice?Answer / 回答Today, the Internet is really convenient for everyone. You can use it to get information from anywhere, anytime. Personally, I like to get new information from Google.How often do you give your friends advice? What advice do you give them?Answer / 回答When my friends make big decisions, I think I will tell them some advice from my experience.Tell your partner a problem you have and let them give you advice.Answer / 回答If you have a partner and you should complete a project together with her/him. I think if you guys communicate with each other regularly is really beneficial for you in completing this project.DrugsWhat illegal drugs can you name?Answer / 回答Okay. In my country, my government doesn't allow anyone near illegal drugs. So I don't know too much about illegal drugs. The only chance for me to know it is on the TV news or in some crime films. Maybe many years ago, there were some karaoke would sell ecstasy to people. Today, however my country doesn't have this phenomenon.What are some drugs that used to be legal but aren't anymore?Answer / 回答I don't know too much about illegal drugs. But I know there are some food additives that have been used in the past but aren't anymore. Because today in my country, people can earn more money. They will pay more attention to their body and food health. Most urban people want to have a house in the countryside. So, they can plant some healthy vegetables.Why are some dangerous drugs legal (alcohol, nicotine, caffeine, etc.) but other safer drugs are illegal? Do you use any legal drugs?Answer / 回答Alcohol is allowed to drink in China, and some elders like to use alcohol as a coping mechanism against hard times. But most contemporary people don't like alcohol.Do you think your country has a drug problem?Answer / 回答Okay, China is one of the safest countries in the world. And I think my government is right. We should keep away from drugs. In China, there many polices to protect us from illegal drugs. We should appreciate them.What would happen if all illegal drugs were legalized?Answer / 回答It's a craze. I will never allow this phenomenon to happen in my country. If the government allows all illegal drugs to become legal, thousands of teenagers will be deceived by some bad person. It's really bad for the development of our country.Where do most illegal drugs come from?Answer / 回答I watched some films about illegal drugs. In these films, most police always say China's illegal drugs come from the Golden Triangle.Why do people start taking illegal drugs?Answer / 回答Okay. I think most people are just deceived by bad person. When we meet these bad people for the first time. The bad people always say "It makes you relax, it's nothing.". And they're giving it away for free for the first time.Is addiction to drugs a disease or a crime?Answer / 回答Both. In my country, taking illegal drugs are being committed. In addition, it can make you addictive, and it's very harmful to your body.Eating HabitsHello, my teacher. Long time no see. I had the flu last week, so I canceled online courses. I'm so sorry.Talk about what you have eaten today.Answer / 回答Okay. I always like to eat buns and eggs for breakfast. On the noon, I would like to eat noodles every day. While discussing this topic. I want to share an interesting thing about myself. My hometown is in northern China where wheat is a traditional crop. I should eat wheat almost every day. In this situation, if my mother cooked rice in one day, It's really novel to me. Prior to university, I really liked to eat rice. However, when I went to my university in southern China. In there, rice is a traditional crop. So you should eat rice almost every day. Since I should eat rice every day, I became exactly hate it.What is the unhealthiest food you can think of?Answer / 回答Personally, I think takeout is the unhealthiest food in the world. When I return home, my mother never allow me to eat takeout. She always said that sellers would use inferior products to cook. And it's really harmful to my health. But you know, when I am at university, I don't have the conditioning to cook. So sometimes, I will eat takeout. Because it's really convenient.What do people need to eat more of?Answer / 回答Actually, I think everyone should eat healthy vegetables and fruits. That's because vegetables and fruits can bring vitamins to our bodies. It's really good for us to grow up, especially children and teenagers.Do you usually eat healthy food?Answer / 回答I really love fruits, especially sweet fruits. I eat fruits almost everyday. Such as apples and bananas.What are some different types of diets you have heard of?Answer / 回答Okay. I never have diets. In my opinion, if you aren't very fat, diets aren't necessary for you. Health is important, but happiness from eating is also very important. However, if you want to have a beautiful body and go fitness, I think you should have strict diets.Do you think your parents ate healthier food at your age than you do? How about your grandparents?Answer / 回答While my parents were at my age, my country is very poor. Most of teenagers liking my parents couldn't have enough food to eat. And during this period, they didn't have takeout.Do you eat a lot of small portions or a few large portions of food each day?Answer / 回答While my parents were at my age, my country is very poor. Most of teenagers, like my parents, couldn't have enough food to eat. And during this period, they didn't have takeout.Do you think people from your country usually eat healthy food?Answer / 回答No. food additives are abused in China. I think my government should pay more attention to this phenomenon and solve this issue quickly. Because it's really bad for our bodies.What are some things about your eating habits you want to change?Answer / 回答okay. I really want to change my eating habits. Because at this moment, I am at university, I should live in my dormitory. I can't have kitchen set to cook by myself. So the first thing before I change my eating habits is that I should have my own kitchen set. Maybe when I become a PhD, I can have a separate room. So I can cook by myself.Do you think that organic food is much better than normal food or are they about the same?Answer / 回答Okay. I really want to change my eating habits. Because I am at university right now, I should be living in my dormitory. I can't have kitchen set to cook by myself. So the first thing I did before I change my eating habits is that I should have my own kitchen set. Maybe when I become a PhD, I can have a separate room. So I can cook by myself.Why is organic food so much more expensive?Answer / 回答About this question, the first reason I think is that if you want to make your vegetables healthier, you must use less fertilizers and pesticides(化肥和农药). You should pay more attention to its growth. Secondly, organic food is smaller than food grown with fertilizers and pesticides. So it can't bring more benefit to sellers.What is traditional food in your country like? Is it healthy?Answer / 回答In my opinion, I think anyone in this world has 100% healthy food. You know framers use fertilizers and pesticides in each of their fields. In addition to, the air and water almost totally is polluted by factories and cars. So, my opinion is that nothing is healthy in this world.Do you think that a person's eating habits affect how long they will live?Answer / 回答Exactly. If a person always eat unhealthy food, it's really bad for his body health. Maybe he will get some strange illness. If you only eat snacks and drinks, you will become very fat, and fat can bring many illnesses to you.Electric CarsHello, teacher. Long time no see. How about you recently? Two weeks ago, I had the flu, and it caused me to have a fever. I felt so terrible. So I canceled my online lessons. This week I went to Shanghai on a business trip with my professor. I returned my university yesterday.Would you ever drive an electric car?Answer / 回答Exactly. There tons of electric cars in China. Especially in big cities, more than half of the cars are electric. However in the countryside, this proportion is smaller. Because it's not convenient to charge in the countryside area.What are some of the benefits and disadvantages of electric cars?Answer / 回答Okay. In my opinion, electricity is an eco-friendly form of energy. In China, almost every family has a car. If all of us use gasoline to drive cars. It will release tons of pollutants to the air and sky. In addition, electric cars can accelerate faster than conventional cars. In China, our government encourages people to buy electric cars. If you choose to buy an electric car, you will pay lower taxes and the government will give you some subsidies. However, electric cars have a fatal shortcoming. The performance of electric cars depends on the battery life. You know batteries is not durable in the winter. In northern China, there has a long cold season, so these places have fewer electric cars.When do you think the first electric car was made?Answer / 回答I don't have too many knowledge about cars. Electric cars began to appear when I was at my university. In my mind, the first brand making electric cars is Tesla. Tesla is a famous company in the world. The CEO of Tesla is Musk. He totally is a genius.What technology needs to be improved to make electric cars more popular?Answer / 回答Okay. As we all know, electric cars are eco-friendly to the environment. However, it's not very convenient for people living in the countryside to use. Because the electricity that electric cars need is not common, you should use professional equipment to charge. If you live or work in big cities, that's OK. But if you live in the countryside, it's not very convenient to charge electric cars. In addition, the batteries of electric cars are not durable in the cold season. So if you live or work in a city which has a long winter, electric cars are not suitable for you.Are electric cars really environmentally friendly if the electricity to charge them comes from coal power plants?Answer / 回答This question is very interesting. I had saw some news about this question. Some people would say "the electricity to charge carts comes from coal power plants", and the coal power plants will release polluting gas to the air. However, what I want to say is electric cars are the future. In the future, gasoline is nonrenewable. But we can use the power of solar and wind to generate electricity.Do you think electric cars, hydrogen cars, or something totally different will be the future of automobiles?Answer / 回答Okay. Firstly, I want to give you my answer to this question. I think we certainly can do it. Let me tell you why. A few days ago, I watched a video about airplanes. It's crazy what we have in just a little over 100 years since the Wright bros made the first airplane in the world. I Can't even imagine what we will have in another 100 years. Compared to airplanes, cars are so easy to make. I think we can do many things that we can't imagine at this time in the future. However, the first thing is to stay away from wars. Wars can destroy any things.How long does it take to charge an electric car?Answer / 回答If you live or work in a big city, it's very convenient for you to use special equipments to charge your electric car. It only takes a few hours to fully charge.Which country do you think is trying the hardest to make electric cars popular?Answer / 回答About this question, the two difficulties about making electric cars popular are these: First, you must have the special equipments to charge; Second, your living area should have a long winter.Do you think electric cars are popular in your country? Why or why not?Answer / 回答Absolutely. In my country, my government attaches great importance to the development of electric cars. Firstly, if you buy an electric car in China, you can pay less in taxes. In addition, it's very difficult for you to buy a car number for a fuel car. Because the government wants to encourage people to buy electric cars and reduce air pollution. Last but not least, the main manufacturers of fuel cars is in Europe and America, it's very difficult for Chinese companies to make engines and transmissions. However, electric cars don't need these components. If every Chinese person buys electric cars, it will be very beneficial to the development of Chinese companies.Do you think electric planes will become popular soon?Answer / 回答I have full confidence in human creativity. As long as we have enough time, let technology and science develop. I think we can do this in the future.
? “THANKS”: 5 Common ExpressionsLinks / 链接Video course link: https://www.youtube.com/watch?v=Io6-K5RrjvUExam link: https://www.engvid.com/thanks-5-common-expressions/Thank you.Usage: We use "thank you" when someone helps us or does something for us.Example / 举例Thank you for your help.Thank you for the gift.Thank you for coming to my party.Tips / 提示When always use:You're welcome.It's no problem.to answer the "thank you".Thanks.Usage: Sometimes, we want to speak a little bit more informally, especially if we're talking with our friends or family. So, another way to say "thank you" is to say "thanks".Example / 举例Thanks.Thanks for your help.Warning / 注意You should never use it in writing business emails, because it's informal.Thanks anyway.Usage: we use "Thanks anyway" or "Thanks anyways" when we ask someone for help, but they can't give us the help that we need.Example / 举例When you are in a supermarket.You: Hi, I'd like to buy some apples. Where are the apples?Waiter: Sorry, we're out of apples. No apples today.You: Thanks anyway.You are talking with your friend Jill.You: Can you drive me to the airport?Jill: I can't. My car is at the mechanic.You: Thanks anyway.Tips / 提示What's the different between "Thanks anyway" or "Thanks anyways"?"Thanks anyway" is more formal. But "Thanks anyways" is more slang.Thanks though.Usage: we use it when someone want to help us or do something for us or offer something, but we don't want the help or the offer. This is a polite way to express I don't need your help.Example / 举例When you are in your office.Mickey: Can I buy you coffee?You: Thanks though.Your friend Barbie wants to invite you to the beach.Barbie: I'm going to the beach, want to come?You: I have my work. Thanks though for the invite.Your friend Tony wants to invite you to dinner.Tony: Can I give you a hand with dinner?You: It's already cooking. Thanks though.Thanks, but no thanks.Usage: It has multiple meanings.Impolite: strong rejection.When telling a story to show strong rejection.You're not interest.Example / 举例When someone want to play a trick on you.Jack: I cooked a dish with insects. want some?Thanks, but no thanks.When you want to complain that your teacher assigned too much homework during the holiday.My teacher wants us to study on the Spring Festival. Thanks, but no thanks.Your mother want to give you some salad.Mother: Yacan, do you want a little bit more salad?Yacan: Ugh. Thanks, but no thanks.? Top 10 Important Phrasal Verbs for Your Daily RoutineLinks / 链接Video course link: https://www.youtube.com/watch?v=z1KVcYIg0HUExam link: https://www.engvid.com/10-phrasal-verbs-daily-routine/sleep inUsage: When you stay in bed latter than usual.Example / 举例On the weekend, I sleep in until 11 am.get upUsage: When you start to leave your bed.Example / 举例I get up at 7 am to take my dog for a walk.head/go outUsage: When you start to leave your house.Example / 举例I head out at 8 am and go to my office.go inUsage: When you return your house.Example / 举例Last night, I went to a party and I got in at about midnight.heat upUsage: When you want to make food becoming warm or hot.Example / 举例At lunch, I take my lunch out of the refrigerator, I go to the microwave, and I heat up my lunch.work outUsage: It means going to exercise.Example / 举例I worked out almost every day when I was in my university.warm upUsage: It means that you do some light exercise to prepare for an activity.Example / 举例I warm up before working out.meet upUsage: It means that you meet your friends informally.Example / 举例I meet up with my friends at a cafe.tidy upUsage: It means organizing your things in the right place.Example / 举例Your room is so messy, and you must tidy it up before our friends come to visit.?️ 6 ways to pronounce ‘-OUGH’ words in EnglishLinks / 链接Video course link: https://www.youtube.com/watch?v=iwl8JnrYwB8Exam link: https://www.engvid.com/6-ways-to-pronounce-ough-words-in-english/Tips / 提示How to learn wards pronouncing ‘-OUGH’?Google pronunciation.Just learning key words.how to pronounce boughtawboughtUsage: "bought" is the past tense of the verb "buy". It speaks liking "bawt".Example / 举例I like to buy books. Yesterday, I bought a book about American history.taughtUsage: "taught" is the past tense of the verb "teach". It speaks liking "tawt".Example / 举例Yesterday, I taught some students English grammar.foughtUsage: "fought" is the past tense of the verb "fight". It speaks liking "fawt".Example / 举例Knights fought in armour.thoughtUsage: "thought" is the past tense of the verb "think". It speaks liking "thawt".Example / 举例Yesterday, I thought about some questions about English grammar.soughtUsage: "sought" is the past tense of the verb "seek". It speaks liking "sawt".Example / 举例Yesterday, I sought about some questions about English grammar.Tips / 提示This is a bit of formal word. We use it more less in our daily life.oughtUsage: "ought" means "should", but it's not a such common word.Example / 举例You ought to help your mother doing some housework.o/oethoughalthoughthorough(thur·ow)dough(dow)Tips / 提示thorough: When we doing some things thoroughly, it means we do them very carefully and diligently.dough: When we make the bread or cookies, we often talk about dough. It's the bread or cookies before we bake it.ufftoughenoughroughow/cowdroughtUsage: Some where has a long time no rain.Pronunciation: drowt.Example / 举例There is three month no raining. It's so drought.ploughUsage: It's framing word. It's a type of machine we use on farmland. It speaks liking "plaw".Example / 举例We always use ploughs to make our farmland becoming more loose.awfcoughUsage: A cough is a reflex action to clear your airways of mucus and irritants such as dust or smoke.Pronunciation: kaaf.Example / 举例Are you sure you don't want to take a chest X-ray?" she inquired, followed by a strategic cough.? Tricky English Question StructuresLinks / 链接Video course link: https://www.youtube.com/watch?v=hKloP6atKXsExam link: https://www.youtube.com/watch?v=hKloP6atKXsExample 01You like pizza, don't you?You don't like pizza, do you?Warning / 注意It's very important for you to understand what the meaning of the question is, not the structure of the question.You just answer the question "Do you like pizza", it's OK. Answer the meaning of the question, not the two verbs inside of it.Example 02Question: Would/Do you mind if I smoked?Your answer is "Yes"(you mind smoking) or "No"(go ahead, I don't care about smoking).Example 04Are you in? - means are you interested, do you want to come or do you want to join us?Example / 举例Someone suggests a few of us are going to the bar for some drinks. He/she will ask you a question "Are you in"?Are you into ...? - When you want to give more information, you can use it.Example / 举例Are you into coming to the bar with us?Are you down with ...? -Example / 举例Are you down with putting in $5 to buy lottery tickets?Are you cool with ...? - means are you agree with that ...?Example / 举例I'm going to take your homework and give it to the teacher, are you cool with that?? Avoid these mistakes talking about moneyLinks / 链接Video course link: https://www.youtube.com/watch?v=oMLn8KwyjjMExam link: https://www.engvid.com/learning-english-avoid-these-mistakes-talking-about-money/Money is uncountableWarning / 注意The first thing that we should know is that "money" is a uncountable noun, so we should never use an "s" on the word "money".In English, we can count dollars, cents, euros, bills. But "money" is a category, we can't count "money".Pronunciation: muh·neeSome usage about moneyWords[ ] Money;[x] Monies;[ ] Much money;[x] Many money;[ ] dollars;[ ] cents;[ ] Pesos;[ ] Euros;[ ] bills(paper money).Phrase[ ] A lot of money;[ ] earn a salary/money;[x] win a salary/money.Something about moneymortgage: mor·guhj(Don't pronounce the "t")bought: bawt(Don't pronounce the "gh")budget: buh·juht(Don't pronounce the "d")cheque/check: chekMaking a killing"Making a killing" means somebody earns a bunch of money.Currency"Currency" should only be used in some situations where we are talking about money from other countries.Spender/Saver"Spender" is a person who like to spend money freely."Saver" is the opposite of "Spender", and is one who likes to save their money.Cost an arm and a legWhen we say "cost an arm and a leg". It means something is really expensive.Warning / 注意You should never use "cost a leg and an arm".? 8 common ways to describe your coworkersLinks / 链接Video course link: https://www.youtube.com/watch?v=JncspyYVjlsExam link: https://www.engvid.com/learn-english-8-common-ways-to-describe-coworkers/A coworker or colleague means somebody you work with.Positive traitsGo-to(noun)"Go-to" means somebody who are most reliable or knowledgeable in your office.Example / 举例She is my go-to when I need help.Going places"Going places" means somebody who is probably to be very successful in the future.Example / 举例Alain is going places.Team player"team player" means somebody who cares about the other members of you team, and focus on helping everyone.Example / 举例Shelly is a team player. I really love doing project with her.Work wife/husband/spouse"Work wife/husband/spouse" means a coworker you are very close to.Example / 举例She is my work wife.Negative traitsNot pulling his/her weight"Not pulling his/her weight" means a coworker who isn't working hard. As result, other members of you office should work harder.Example / 举例Becky is not pulling her weight here. She's lazy.Slacker"Slacker" means somebody who is very lazy and doesn't work hard.Example / 举例Becky is a slacker. She's lazy.Two-faced"Two-faced" means somebody who appears one way when they talk to you, but then they do something differently when they talk to other people.Example / 举例Watch out! Peter can't be trusted. He is two-faced.Stress case"Stress case" means somebody who always panics or overreacts.Example / 举例Jenny is a stress case. She has called me 8 times since lunch.?️ Silent Letters in English: The 10 Most Common Errorssilent letterswordstranslationpronunciationpspsychologyn. 心理学；心理；心理特征；心理影响sai·kaa·luh·jee psychiatristn. 精神科医生；精神病学家sai·kai·uh·truhstknkneeln. 跪vt. 跪；跪着neel known. 知识；知晓，知情v. 知道，知悉，了解now knifen. 刀vt. 用刀伤害(或杀害)nife knockn. 敲击；撞击（声）v. 敲；击naakwrwritev. 写；编写；写信rite wristn. 手腕；腕关节rist warpv. 弯曲；(使)扭曲，变形/ˈwɔɹp/ wrestlen. 搏斗；摔跤；力战v. 摔跤；奋力对付；努力处理reh·slhhourn. 小时；约一小时的时间au·ur honestadj. 诚实的，老实的adv. 的确，确实，真的aa·nuhst honorn. 荣誉；荣幸v. 尊敬；尊重aa·nrdgjudgen. 法官；审判员v. 判断；断定；认为juhj pledgen. 誓言；保证；抵押v. 抵押；保证给予(或做)plej dodgen. 推脱的计策；逃避的诡计v. 闪避；躲开；避开；闪开daajtfastenv. 系牢；(使两部分)扎牢fa·sn listenn. 听vi. （注意地）听；听从li·sn balletn. 芭蕾；芭蕾舞团/ˈbæleɪ/ buffetn. 自助餐；(火车)饮食柜台vt. 对…打来打去；将…推来搡去adj. 快餐式的buh-FAYmbcombn. 梳子；篦子；压发梳v. 梳，梳理（头发）kowm climbn. 攀登；攀缘v. 攀登；(尤指吃力地向某处)爬klime lambn. 羔羊；小羊lam thumbn. 拇指thuhmlkwalkn. 步行；行走；散步v. 步行；走；行走waak talkn. 交谈；谈话；讨论v. 说话；讲话；谈话taak chalkn. 粉笔chaak folkn. 人们；各位；大伙儿；亲属；家属fowktlecastlen. 城堡；堡垒ka·sl whistlen. 哨子；口哨v. 吹口哨；打呼哨wi·slmncolumnn. 柱；(通常为)圆形石柱kaa·luhm autumnn. 秋，秋季aa·tm hymnn. 赞美诗；圣歌(尤指基督徒唱的颂扬上帝的歌)him solemnadj. 冷峻的；表情严肃的；庄严的saa·luhm
How to be diplomatic (soften your language).In some situations, you should be diplomatic so that you don't offend anyone.Soften the language.Using modal verbs.Using negative questions.1、This isn't good.I think this could be improved.Why don't we make some changes to this?Wouldn't this be better if ...?2、I don't like it.I'm not a big fan of ...I'm not too keen on ......isn't really to my taste....isn't my cup of tea....isn't my favourite.3、You're wrong.I'm not sure that's entirely accurate.4、I’m not satisfied with ...This isn't exactly what I had in my mind.I have some reservations.5、Would you like to ...I'm afraid I can't at the moment. Actually, I need to ...6、Can I speak to you?I'm a bit busy right now, could you give me a minute?That's not a bad idea, but .... How about ...?Do/Did you Know ...?What's the difference between "Did you know" and "Do you know"?Do you know means we don't know something(the answer) but we want to get some infoamation.Examples:Do you know where I can buy nice shoes?Do you know what time it is?Did you know means I know this and I want to tell you a fact. I want to see you know the fact.Examples:Did you know Mt.everst is the world's tallest mountain?Did you know the most common letter in English is "e"?Phrasal Verbs: Add “OUT” to change the meaning of these 8 verbsCourse Links: https://www.youtube.com/watch?v=Cr4RYby4L6QLooking out is used to mean be careful(Look out! The car is coming!). On the other hand, Look out means paying attention to do something(I will look out for you!).Staying out means to remain out for a long time(He stayed out all night).Pointing out means to draw your attention to something(My teacher always points out my mistakes carefully).Working out means exercising or solving questions.Wearing out means something can't be used anymore it's no longer usable. On the other hand, it means you're so tired(My T-shirt is wearing out. / Today I'm so tired, I am wearing out.).Ruling out means you eliminate the possibility(You can rule out possibilities or answers).Standing out means something or some people are noticeableFinding out means discovering or learning(I need to find out when the next bus leaves).Preposition tips秒懂英语里的常见介词～Psychology of investingrational, 理性的. emotional: 情绪化的. fear, adj.害怕的. greed, adj.贪婪的.dip one's toes into something: 正尝试做某事. I heard my favorite friend is dipping his toes into investing.I suggest .... as well: 我建议……. I suggest you can learn some new skills about your work as well.In my humble opinion, 在我看来. In my humble opinion, they're ultimately driven by emotions like fear and greed.a pretty significant edge, 一个很大的优势. Understanding this will give you a pretty significant edge.playing the long game, 做长期打算. Emotional states can make us forget we're playing the long game.Looks like I've got some homework to do, 看来我还有一些功课需要去做.Children with AutismWordsautism: n.自闭症coping mechanisms: 应对机制. A lot of people use alcohol as a coping mechanism against hard times.overload: n. 超载, 过重的负担. Nobody says that money is overload. / When I get into this meeting, people dump lots of information on me. I can't process it, so it's information overload for me.stimulation: n. 刺激; 鼓励; 鼓舞. stimming: n. 自我刺激行为.Phraseinformation gap: There is a huge information gap out there. And when we're dealing with little kids at school, sometimes an information gap can be really dangerous thing.Adjusting to a new environment: 适应新环境.somebody's first encounter with...: 第一次遇到…….That's great advice!