np.prod()用法解析

frame = np.frombuffer( binary_data[offset:offset + image_size_bytes], dtype=np.uint16 ).reshape((height, width))

import pcl import numpy as np import open3d as o3d def ransac_line_fitting(input_path, output_path): # 读取PCD文件 cloud = pcl.load(input_path) # 转换点云为numpy数组便于处理 points = cloud.to_array() # RANSAC直线拟合 model = pcl.SampleConsensusModelLine(cloud) ransac = pcl.RandomSampleConsensus(model) ransac.set_DistanceThreshold(0.01) ransac.computeModel() coeffs = ransac.getModelCoefficients() # 格式: [x0,y0,z0,dx,dy,dz] # 提取直线参数 line_point = np.array(coeffs[:3]) line_dir = np.array(coeffs[3:6]) line_dir /= np.linalg.norm(line_dir) # 归一化方向向量 # 确保方向向量X分量为正 if line_dir[0] < 0: line_dir *= -1 # 找到原始点云中X最小值点 min_x_idx = np.argmin(points[:, 0]) min_x_point = points[min_x_idx] # 将X最小值点投影到直线上 vec = min_x_point - line_point t = np.dot(vec, line_dir) # 投影参数计算~ start_point = line_point + t * line_dir # 计算总长度和步长 max_x_idx = np.argmax(points[:, 0]) max_x_point = points[max_x_idx] total_length = np.linalg.norm(max_x_point - min_x_point) step = total_length / (len(points) - 1) # 生成等间距点云 new_points = [] for i in range(len(points)): new_p = start_point + i * step * line_dir new_points.append(new_p) # 创建并保存新点云 new_cloud = pcl.PointCloud() new_cloud.from_array(np.array(new_points, dtype=np.float32)) pcl.save(new_cloud, output_path) if __name__ == "__main__": ransac_line_fitting("input.pcd", "output.pcd") 请把其中的pcl用open3d替代

import numpy as np import matplotlib.pyplot as plt import torch from torch import nn, optim import torch.nn.functional as F from scipy.stats import entropy from sklearn.neighbors import KernelDensity # 设置字体 plt.rcParams["font.family"] = "Times New Roman" plt.rcParams["axes.labelsize"] = 12 plt.rcParams["axes.titlesize"] = 14 # 辅助函数：将整数转换为one-hot编码向量 def one_hot(a, M): onehot = np.zeros(M) onehot[a] = 1 return onehot # 模值归一化函数 def magnitude_normalization(x): """ 模值归一化：将星座点缩放到单位圆内 确保所有点的模值不超过1 """ # 计算每个点的模值 (sqrt(real^2 + imag^2)) magnitudes = torch.sqrt(torch.sum(x**2, dim=1)) # 找到最大模值 max_magnitude = torch.max(magnitudes) # 避免除以零 if max_magnitude < 1e-6: return x # 归一化所有点 return x / max_magnitude.unsqueeze(1) def awgn(x, sigma2): noise_t = np.sqrt(sigma2) * torch.randn(x.shape) return torch.add(x, noise_t) # 计算互信息 def calculate_mutual_information(x, y, n_bins=50): """ 计算两个连续变量x和y之间的互信息 使用直方图方法估计联合分布和边际分布 """ # 将PyTorch张量转换为NumPy数组 if isinstance(x, torch.Tensor): x = x.detach().numpy() if isinstance(y, torch.Tensor): y = y.detach().numpy() # 确保输入是二维数组 if x.ndim == 1: x = x.reshape(-1, 1) if y.ndim == 1: y = y.reshape(-1, 1) # 合并数据 joint_data = np.hstack((x, y)) # 计算联合直方图 hist_result = np.histogramdd( joint_data, bins=n_bins, density=True ) # 解析返回值 if len(hist_result) == 2: # 2D情况: (H, edges) joint_hist = hist_result[0] edges = hist_result[1] x_edges = edges[0] y_edges = edges[1] elif len(hist_result) == 3: # 3D情况: (H, x_edges, y_edges) joint_hist, x_edges, y_edges = hist_result else: raise ValueError("Unsupported histogramdd return format") # 确保直方图是二维的 if joint_hist.ndim != 2: joint_hist = joint_hist.reshape(n_bins, n_bins) # 计算边际直方图 - 确保返回一维数组 x_hist, _ = np.histogram(x.flatten(), bins=x_edges, density=True) y_hist, _ = np.histogram(y.flatten(), bins=y_edges, density=True) # 确保边际直方图是一维的 x_hist = x_hist.flatten() y_hist = y_hist.flatten() # 计算互信息 mi = 0.0 # 计算箱宽（假设等间距） dx = np.diff(x_edges)[0] if len(x_edges) > 1 else 1.0 dy = np.diff(y_edges)[0] if len(y_edges) > 1 else 1.0 # 遍历所有箱子 for i in range(joint_hist.shape[0]): for j in range(joint_hist.shape[1]): # 获取标量值 p_xy_val = joint_hist[i, j] * dx * dy # 确保索引在范围内 if i < len(x_hist) and j < len(y_hist): p_x_val = x_hist[i] * dx p_y_val = y_hist[j] * dy else: # 如果索引超出范围，跳过 continue # 避免除以零和log(0) - 现在所有值都是标量 if p_xy_val > 1e-10 and p_x_val > 1e-10 and p_y_val > 1e-10: mi += p_xy_val * np.log2(p_xy_val / (p_x_val * p_y_val)) return mi # 生成标准32QAM星座图（均匀分布） def generate_uniform_32qam(): # 创建6x6网格 levels = np.array([-1.5, -0.9, -0.3, 0.3, 0.9, 1.5]) x, y = np.meshgrid(levels, levels) # 展平网格 points = np.vstack((x.flatten(), y.flatten())).T # 移除四个角点（32QAM标准） corners = [ [0, 0], # (-1.5, -1.5) [0, 5], # (-1.5, 1.5) [5, 0], # (1.5, -1.5) [5, 5] # (1.5, 1.5) ] # 移除角点 mask = np.ones(len(points), dtype=bool) for corner in corners: idx = corner[0] * 6 + corner[1] mask[idx] = False uniform_32qam = points[mask] # 模值归一化 magnitudes = np.sqrt(np.sum(uniform_32qam**2, axis=1)) max_magnitude = np.max(magnitudes) uniform_32qam = uniform_32qam / max_magnitude return uniform_32qam # 计算平均功率 def calculate_average_power(constellation): """ 计算星座图的平均功率 平均功率 = (1/N) * Σ |x_i|^2 """ magnitudes_sq = np.sum(constellation**2, axis=1) return np.mean(magnitudes_sq) class Mapper(nn.Module): """ 映射器网络（带隐藏层） """ def __init__(self): super().__init__() self.lin1 = nn.Linear(32, 64) self.lin2 = nn.Linear(64, 64) self.lin3 = nn.Linear(64, 2) def forward(self, y): y = F.relu(self.lin1(y)) y = F.relu(self.lin2(y)) return self.lin3(y) class Demapper(nn.Module): """ 解映射器网络（带隐藏层） """ def __init__(self): super().__init__() self.lin1 = nn.Linear(2, 64) self.lin2 = nn.Linear(64, 64) self.lin3 = nn.Linear(64, 32) def forward(self, y): y = F.relu(self.lin1(y)) y = F.relu(self.lin2(y)) return self.lin3(y) # 生成训练数据 M = 32 n = 10000 a = np.random.choice(range(M), size=n) onehot = np.array([one_hot(a[i], M) for i in range(n)]) # 转换为PyTorch张量 onehot_t = torch.tensor(onehot).float() a_t = torch.tensor(a).long() # 直接转换为长整型 # 整个星座图的one-hot编码（用于计算最小距离） onehot_plot = np.array([one_hot(i, M) for i in range(M)]) onehot_plot_t = torch.tensor(onehot_plot).float() # 损失函数 loss_fn = nn.CrossEntropyLoss() # 生成均匀星座图并计算其平均功率 uniform_32qam = generate_uniform_32qam() P_uniform = calculate_average_power(uniform_32qam) # 均匀星座图的平均功率 print(f"Uniform constellation average power: {P_uniform:.4f}") SNRdBs = np.array([5, 10, 15]) SNRs = 10 ** (SNRdBs / 10) # 初始化网络 mapper = Mapper() demap = Demapper() # 存储互信息结果 mi_results = [] # 训练循环（针对不同SNR） for (k, snr) in enumerate(SNRs): # 使用均匀星座图的平均功率计算噪声方差 sigma2 = P_uniform / snr # 噪声方差 print(f'\n---SNR = {SNRdBs[k]} dB---') print(f'Using uniform constellation average power: {P_uniform:.4f}') print(f'Noise variance: {sigma2:.6f}') # 优化器 optimizer = optim.Adam(list(mapper.parameters()) + list(demap.parameters()), lr=0.005) # 训练循环 for j in range(10001): optimizer.zero_grad() # 1. 生成训练批次的星座点 xhat = mapper(onehot_t) # 2. 生成整个星座图的星座点（用于计算最小距离） x_plot = mapper(onehot_plot_t) # 3. 模值归一化（使用整个星座图的最大模值） magnitudes = torch.sqrt(torch.sum(x_plot**2, dim=1)) max_mag = torch.max(magnitudes) if max_mag > 0: xhat_norm = xhat / max_mag x_plot_norm = x_plot / max_mag else: xhat_norm = xhat x_plot_norm = x_plot # 4. 通过AWGN信道 yhat = awgn(xhat_norm, sigma2) # 5. 解映射 l = demap(yhat) ce_loss = loss_fn(l, a_t) # 6. 安全计算最小距离（避免就地操作） # 计算所有点对之间的欧氏距离 dist_matrix = torch.cdist(x_plot_norm, x_plot_norm, p=2) # 创建对角线掩码（避免就地操作） eye_mask = torch.eye(M, dtype=torch.bool, device=dist_matrix.device) # 将对角线位置设为无穷大（非就地操作） dist_matrix_masked = dist_matrix.masked_fill(eye_mask, float('inf')) # 找到最小距离 min_dist = torch.min(dist_matrix_masked) # 7. 计算正则化项（惩罚过小的最小距离） d_min_threshold = 0.3 # 最小距离阈值（模值归一化后需要更大阈值） reg_loss = 10.0 * F.relu(d_min_threshold - min_dist) # 使用ReLU避免条件判断 total_loss = ce_loss + reg_loss # 反向传播 total_loss.backward() optimizer.step() # 每500次迭代打印损失 if j % 500 == 0: ce_bits = ce_loss.item() / np.log(2) reg_val = reg_loss.item() # 计算当前星座点平均模值 avg_magnitude = torch.mean(magnitudes).item() # 计算当前星座图平均功率 with torch.no_grad(): avg_power = torch.mean(torch.sum(x_plot_norm**2, dim=1)).item() print(f'Epoch {j}: CE Loss = {ce_bits:.4f} bits, Reg Loss = {reg_val:.4f}, ' f'Min Dist = {min_dist.item():.4f}, Avg Mag = {avg_magnitude:.4f}, ' f'Avg Power = {avg_power:.4f}') # 提前终止条件 if total_loss < 1e-3 and min_dist > d_min_threshold: print(f"Early stopping at epoch {j}") break # 可视化学习到的星座图 with torch.no_grad(): # 生成星座点 learned_x = mapper(onehot_plot_t) # 模值归一化 magnitudes = torch.sqrt(torch.sum(learned_x**2, dim=1)) max_mag = torch.max(magnitudes) if max_mag > 0: learned_x = learned_x / max_mag learned_x = learned_x.numpy() # 计算模值 learned_magnitudes = np.sqrt(learned_x[:, 0]**2 + learned_x[:, 1]**2) # 计算平均功率 learned_avg_power = np.mean(learned_magnitudes**2) plt.figure(figsize=(10, 8)) # 绘制星座图 plt.scatter(learned_x[:, 0], learned_x[:, 1], s=30, c=learned_magnitudes, cmap='viridis') # 绘制单位圆 theta = np.linspace(0, 2*np.pi, 100) plt.plot(np.cos(theta), np.sin(theta), 'r--', alpha=0.5) # 标记星座点索引 for i, (x, y) in enumerate(learned_x): plt.text(x, y, str(i), fontsize=8, ha='center', va='center') plt.title(f'32QAM, SNR = {SNRdBs[k]} dB') plt.xlabel('In-phase Component') plt.ylabel('Quadrature Component') plt.grid(True, linestyle='--', alpha=0.6) plt.axis('equal') plt.xlim(-1.2, 1.2) plt.ylim(-1.2, 1.2) # 添加颜色条显示模值 cbar = plt.colorbar() cbar.set_label('Magnitude') plt.tight_layout() # plt.savefig(f'optimized_constellation_snr_{SNRdBs[k]}dB.png', dpi=300) plt.show() # 打印模值统计信息 print(f"Magnitude statistics for SNR={SNRdBs[k]}dB:") print(f" Max magnitude: {np.max(learned_magnitudes):.4f}") print(f" Min magnitude: {np.min(learned_magnitudes):.4f}") print(f" Avg magnitude: {np.mean(learned_magnitudes):.4f}") print(f" Std magnitude: {np.std(learned_magnitudes):.4f}") print(f" Average power: {learned_avg_power:.4f} (Uniform: {P_uniform:.4f})") # ================= 互信息计算 ================= # 生成测试数据（10000个样本） test_a = np.random.choice(range(M), size=10000) test_onehot = np.array([one_hot(test_a[i], M) for i in range(10000)]) test_onehot_t = torch.tensor(test_onehot).float() with torch.no_grad(): # 优化后星座图的发送信号 test_xhat = mapper(test_onehot_t) # 模值归一化 test_x_plot = mapper(onehot_plot_t) test_magnitudes = torch.sqrt(torch.sum(test_x_plot**2, dim=1)) test_max_mag = torch.max(test_magnitudes) if test_max_mag > 0: test_xhat_norm = test_xhat / test_max_mag else: test_xhat_norm = test_xhat # 通过AWGN信道 test_yhat = awgn(test_xhat_norm, sigma2) # 计算优化后星座图的互信息 mi_optimized = calculate_mutual_information( test_xhat_norm, test_yhat, n_bins=50 ) # 计算标准均匀分布星座图的互信息 # 为测试数据选择星座点 uniform_x = uniform_32qam[test_a] # 添加AWGN噪声 - 使用相同的噪声方差 uniform_y = uniform_x + np.sqrt(sigma2) * np.random.randn(*uniform_x.shape) # 计算均匀分布星座图的互信息 mi_uniform = calculate_mutual_information( uniform_x, uniform_y, n_bins=50 ) # 存储结果 mi_results.append({ 'SNR_dB': SNRdBs[k], 'mi_uniform': mi_uniform, 'mi_optimized': mi_optimized, 'optimized_power': learned_avg_power, 'uniform_power': P_uniform }) # 打印互信息结果 print(f"\nMutual Information (SNR={SNRdBs[k]}dB):") print(f" Uniform Constellation: {mi_uniform:.4f} bits (Avg Power: {P_uniform:.4f})") print(f" Optimized Constellation: {mi_optimized:.4f} bits (Avg Power: {learned_avg_power:.4f})") print(f" Improvement: {(mi_optimized - mi_uniform):.4f} bits ({100*(mi_optimized - mi_uniform)/mi_uniform:.2f}%)") # 绘制互信息随SNR变化的曲线 plt.figure(figsize=(10, 6)) snrs = [res['SNR_dB'] for res in mi_results] mi_uni = [res['mi_uniform'] for res in mi_results] mi_opt = [res['mi_optimized'] for res in mi_results] plt.plot(snrs, mi_uni, 'bo-', label=f'Uniform Constellation (Avg Power: {P_uniform:.4f})') plt.plot(snrs, mi_opt, 'ro-', label='Optimized Constellation') plt.xlabel('SNR (dB)') plt.ylabel('Mutual Information (bits)') plt.title('Mutual Information Comparison (Using Uniform Constellation Power)') plt.grid(True, linestyle='--', alpha=0.6) plt.legend() plt.tight_layout() # plt.savefig('mutual_information_comparison_uniform_power.png', dpi=300) plt.show() # 绘制功率比较 plt.figure(figsize=(10, 6)) opt_power = [res['optimized_power'] for res in mi_results] uni_power = [res['uniform_power'] for res in mi_results] plt.bar([s - 0.2 for s in snrs], uni_power, width=0.4, label='Uniform Constellation', alpha=0.7) plt.bar([s + 0.2 for s in snrs], opt_power, width=0.4, label='Optimized Constellation', alpha=0.7) plt.xlabel('SNR (dB)') plt.ylabel('Average Power') plt.title('Average Power Comparison') plt.grid(True, axis='y', linestyle='--', alpha=0.6) plt.legend() plt.tight_layout() # plt.savefig('power_comparison.png', dpi=300) plt.show() # 保存互信息结果 # import pandas as pd # mi_df = pd.DataFrame(mi_results) # mi_df.to_csv('mutual_information_results_uniform_power.csv', index=False) # print("\nMutual information results saved to 'mutual_information_results_uniform_power.csv'") 代码运行之后显示： --------------------------------------------------------------------------- ValueError Traceback (most recent call last) ~\AppData\Local\Temp\ipykernel_9024\1337079151.py in <module> 381 test_xhat_norm, 382 test_yhat, --> 383 n_bins=50 384 ) 385 ~\AppData\Local\Temp\ipykernel_9024\1337079151.py in calculate_mutual_information(x, y, n_bins) 84 # 确保直方图是二维的 85 if joint_hist.ndim != 2: ---> 86 joint_hist = joint_hist.reshape(n_bins, n_bins) 87 88 # 计算边际直方图 - 确保返回一维数组 ValueError: cannot reshape array of size 6250000 into shape (50,50)

import numpy as np import itertools def read_lammps_dump(file_path): """ 读取 LAMMPS dump 文件并返回原子坐标和盒子尺寸 (修正盒子尺寸处理，支持非零起始坐标) """ with open(file_path, 'r') as f: lines = f.readlines() box_bounds = np.zeros((3, 2)) for i, line in enumerate(lines): if 'ITEM: NUMBER OF ATOMS' in line: num_atoms = int(lines[i + 1].strip()) elif 'ITEM: BOX BOUNDS' in line: for dim in range(3): lo, hi = map(float, lines[i+1+dim].split()[:2]) box_bounds[dim] = [lo, hi] box_size = box_bounds[:, 1] - box_bounds[:, 0] volume = np.prod(box_size) elif 'ITEM: ATOMS' in line: start_index = i + 1 break # 读取原子坐标 (修正为读取任意列位置) coordinates = np.zeros((num_atoms, 3)) for i in range(num_atoms): line = lines[start_index + i] data = line.split() # 查找x,y,z的列位置（假设最后三列是坐标） coordinates[i] = np.array([float(data[-3]), float(data[-2]), float(data[-1])]) return coordinates, box_size, volume def compute_rdf(coordinates, box_size, volume, num_bins=100, r_max=10.0): """ 计算径向分布函数 (RDF) (关键修正：周期性边界条件 + 正确归一化) """ num_atoms = len(coordinates) dr = r_max / num_bins rdf_hist = np.zeros(num_bins) # 预计算所有原子坐标的周期性修正 coordinates = coordinates % box_size # 确保原子在盒子内 # 计算所有原子对之间的距离（带周期性边界条件） for i, j in itertools.combinations(range(num_atoms), 2): delta = coordinates[j] - coordinates[i] # 应用最小镜像约定 delta -= box_size * np.round(delta / box_size) dist = np.linalg.norm(delta) bin_idx = int(dist / dr) if bin_idx < num_bins: rdf_hist[bin_idx] += 1 # 修正：每个原子对只统计一次 # 计算归一化因子 r = np.linspace(dr/2, r_max - dr/2, num_bins) # bin中心点坐标 shell_volumes = 4/3 * np.pi * ((r + dr/2)**3 - (r - dr/2)**3) density = num_atoms / volume # 修正归一化公式 (考虑原子对总数) total_pairs = num_atoms * (num_atoms - 1) / 2 rdf = rdf_hist / (shell_volumes * density * total_pairs) return r, rdf # 使用示例 -------------------------------------------------- # 读取文件（注意现在返回box_size和volume） coordinates, box_size, volume = read_lammps_dump('Cu_npt.xyz') # 计算RDF（保持与LAMMPS相同的参数设置） r, g_r = compute_rdf(coordinates, box_size, volume, num_bins=100, r_max=10.0) # 保存结果（添加注释说明） with open('rdf_py_corrected.data', 'w') as f: f.write(f"# 使用参数：bins=100, r_max=10.0, 原子数={len(coordinates)}\n") f.write('# BinID Radius(A) g(r)\n') for i in range(len(r)): f.write(f'{i+1} {r[i]:.4f} {g_r[i]:.6f}\n') print("修正版RDF计算完成，结果已保存到 rdf_py_corrected.data")这个脚本rdf峰值不到一，明显不对，而且峰太多。请大刀阔斧地修改，使它的数据能与lammps计算结果对得上