"""
相干伊辛机（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: 退火步数
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 is not None else 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
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] = 1
return 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
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 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)
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}")
def get_best_state(self) -> Tuple[np.ndarray, float]:
"""返回最低能量状态及其能量"""
idx = np.argmin(self.energy_history)
return self.spins.copy(), self.energy()
def example_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}")
def example_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
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()
