""" 相干伊辛机（CIM）模拟器 基于测量反馈和量子噪声的离散自旋演化算法。 适用于求解伊辛模型基态（最小能量）问题。 """import numpy as np import matplotlib.pyplot as plt from typing import Optional, Tuple, List, Callable classCoherentIsingMachine:""" 相干伊辛机模拟器 参数： n_spins: 自旋数量 J: 耦合矩阵 (n_spins x n_spins)，对称且对角元为0 h: 外磁场 (n_spins,) noise_amplitude: 初始噪声幅度（模拟量子涨落） annealing_steps: 退火步数 update_rule: 更新规则函数，默认使用带有噪声的阈值更新 """def__init__( self, n_spins:int, J: np.ndarray, h: Optional[np.ndarray]=None, noise_amplitude:float=1.0, annealing_steps:int=1000, update_rule: Optional[Callable]=None,): self.n_spins = n_spins self.J = J self.h = h if h isnotNoneelse np.zeros(n_spins) self.noise_amplitude = noise_amplitude self.annealing_steps = annealing_steps self.update_rule = update_rule if update_rule else self._default_update # 随机初始化自旋（+1/-1） self.spins = np.random.choice([-1,1], size=n_spins) self.energy_history =[]def_default_update(self, spins: np.ndarray, local_fields: np.ndarray, noise:float)-> np.ndarray:""" 默认更新规则： 计算每个自旋的有效场，加上噪声，然后通过符号函数决定新自旋。 新自旋 = sign(local_field + noise) """ noisy_field = local_fields + noise * np.random.randn(self.n_spins) new_spins = np.sign(noisy_field)# 避免零值 new_spins[new_spins ==0]=1return new_spins def_local_fields(self, spins: np.ndarray)-> np.ndarray:"""计算每个自旋感受到的局部场：h_i_eff = sum_j J_ij s_j + h_i"""return self.J @ spins + self.h defenergy(self, spins: Optional[np.ndarray]=None)->float:"""计算伊辛能量：E = -0.5 * sum_ij J_ij s_i s_j - sum_i h_i s_i"""if spins isNone: spins = self.spins return-0.5* np.sum(self.J * np.outer(spins, spins))- np.sum(self.h * spins)defanneal(self, verbose:bool=True):""" 执行模拟退火过程（对应CIM中的噪声逐渐降低） """for step inrange(self.annealing_steps):# 当前噪声幅度，线性退火 noise = self.noise_amplitude *(1.0- step / self.annealing_steps)# 计算局部场 local_fields = self._local_fields(self.spins)# 更新自旋（并行） new_spins = self.update_rule(self.spins, local_fields, noise)# 接受新状态（确定性更新，因为目标是下降） self.spins = new_spins # 记录能量 e = self.energy() self.energy_history.append(e)if verbose and(step +1)%(self.annealing_steps //10)==0:print(f"Step {step+1}/{self.annealing_steps}, Energy: {e:.4f}")defget_best_state(self)-> Tuple[np.ndarray,float]:"""返回最低能量状态及其能量""" idx = np.argmin(self.energy_history)# 注意：我们记录的是每一步的能量，但状态只保留了最后一步# 为了得到真正的最优状态，应该在每次更新后保存状态# 这里简化处理：重新模拟一次并跟踪最佳状态（略）# 演示中直接返回最终状态（模拟退火通常最终状态接近最优）return self.spins.copy(), self.energy()defexample_random_ising():"""示例1：随机伊辛模型""" n =50 np.random.seed(42) J = np.random.randn(n, n) J =(J + J.T)/2# 对称化 np.fill_diagonal(J,0)# 无自相互作用 h = np.random.randn(n) cim = CoherentIsingMachine(n_spins=n, J=J, h=h, noise_amplitude=2.0, annealing_steps=500) cim.anneal(verbose=True) plt.figure(figsize=(10,4)) plt.plot(cim.energy_history) plt.xlabel("Iteration") plt.ylabel("Energy") plt.title("Energy evolution for random Ising model") plt.grid(True) plt.show()print(f"Final energy: {cim.energy():.4f}")defexample_max_cut():""" 示例2：最大割问题（Max-Cut） 将无向图的最大割问题映射为伊辛模型： 割的权重 = 0.5 * sum_{i<j} w_ij (1 - s_i s_j) 能量最小化等价于最大化割权重。 因此，设置 J_ij = -w_ij, h=0。 """# 构造一个简单的图（如三角形加一个额外节点）# 图权重矩阵（对称） W = np.array([[0,1,2,0],[1,0,1,3],[2,1,0,2],[0,3,2,0]]) n = W.shape[0] J =-W # 伊辛耦合为负权重# 确保对角为0 np.fill_diagonal(J,0) cim = CoherentIsingMachine(n_spins=n, J=J, h=None, noise_amplitude=1.5, annealing_steps=300) cim.anneal(verbose=True)# 输出结果 partition = cim.spins cut_weight =0.5* np.sum(W *(1- np.outer(partition, partition)))print(f"Final partition: {partition}")print(f"Cut weight: {cut_weight:.2f} (Max possible: {np.sum(W)/2:.2f})") plt.figure(figsize=(10,4)) plt.plot(cim.energy_history) plt.xlabel("Iteration") plt.ylabel("Energy") plt.title("Energy evolution for Max-Cut problem") plt.grid(True) plt.show()if __name__ =="__main__":print("=== Random Ising Model ===") example_random_ising()print("\n=== Max-Cut Problem ===") example_max_cut()
