"""
相干伊辛机（CIM）模拟器
基于测量反馈和量子噪声的离散自旋演化算法。
适用于求解伊辛模型基态（最小能量）问题。
"""
import numpy as np
import matplotlib.pyplot as plt
from typing import Optional, Tuple, List, Callable
class CoherentIsingMachine:
"""
相干伊辛机模拟器
参数：
n_spins: 自旋数量
J: 耦合矩阵 (n_spins x n_spins)，对称且对角元为 0
h: 外磁场 (n_spins,)
noise_amplitude: 初始噪声幅度（模拟量子涨落）
annealing_steps: 退火步数
"""
def __init__(self, n_spins: int, J: np.ndarray, h: Optional[np.ndarray] = None,
noise_amplitude: float = 1.0, annealing_steps: int = 1000):
self.n_spins = n_spins
self.J = J
self.h = h if h is not None else np.zeros(n_spins)
self.noise_amplitude = noise_amplitude
self.annealing_steps = annealing_steps
self.spins = np.random.choice([-1, 1], size=n_spins)
self.energy_history = []
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
def energy(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 is None:
spins = self.spins
return -0.5 * np.sum(self.J * np.outer(spins, spins)) - np.sum(self.h * spins)
def update_step(self, local_fields: np.ndarray, noise: float) -> np.ndarray:
"""
更新规则：
计算每个自旋的有效场，加上噪声，然后通过符号函数决定新自旋。
这里模拟 CIM 的测量反馈过程：局部场加上高斯噪声（模拟量子涨落），
再通过符号函数决定新自旋。实际运行时会遇到零值情况，需特殊处理。
"""
noisy_field = local_fields + noise * np.random.randn(self.n_spins)
new_spins = np.sign(noisy_field)
new_spins[new_spins == 0] = 1
return new_spins
def anneal(self, verbose: bool = True):
"""执行模拟退火过程（对应 CIM 中的噪声逐渐降低）"""
for step in range(self.annealing_steps):
noise = self.noise_amplitude * (1.0 - step / self.annealing_steps)
local_fields = self._local_fields(self.spins)
self.spins = self.update_step(local_fields, noise)
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}")
def get_best_state(self) -> Tuple[np.ndarray, float]:
"""返回最低能量状态及其能量"""
idx = np.argmin(self.energy_history)
return self.spins.copy(), self.energy()
if __name__ == "__main__":
print("=== Random Ising Model ===")
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}")
print("\n=== Max-Cut Problem ===")
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
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})")
