# Ising Model（伊辛模型）介紹

伊辛模型（Ising Model）是統計物理學中的經典模型，用於描述磁性物質中的自旋系統。由物理學家 **Wilhelm Lenz** 提出，其學生 **Ernst Ising** 在 1924 年進一步發展。該模型在研究相變、臨界現象以及分布式系統中具有廣泛應用。

---

## 1. 模型基礎

### 1.1 自旋系統
- **系統組成**：由一個 $d$ 維格子上的 **自旋變量** 組成。
- **自旋取值**：
- 每個格點 $i$ 的自旋 $s_i$ 取值為 $+1$（向上自旋）或 $-1$（向下自旋）。
- 自旋表示微觀粒子的磁矩方向。

---

### 1.2 哈密頓量（Hamiltonian）
- **總能量公式**：
$$
H = -J \sum_{\langle i, j \rangle} s_i s_j - h \sum_i s_i
$$
其中：
- $J$：**相互作用強度**，描述鄰居自旋之間的耦合程度。
- 當 $J > 0$：鄰居自旋傾向一致（鐵磁性）。
- 當 $J < 0$：鄰居自旋傾向相反（反鐵磁性）。
- $\langle i, j \rangle$：表示格子中所有相鄰自旋對。
- $h$：**外部磁場強度**。
- $s_i$：格點 $i$ 的自旋狀態（取值 $+1$ 或 $-1$）。

---

### 1.3 相互作用機制
- **局部交互**：
- 每個格點的自旋狀態受鄰居自旋和外部磁場的共同影響。
- **能量最小化**：
- 系統傾向於選擇自旋分佈，使總能量 $H$ 最小。

---
# 基於伊辛模型的自旋格子強化學習

## 程式碼

import torch
import numpy as np
import matplotlib.pyplot as plt
from collections import deque
from random import shuffle

device = torch.device("cuda" if torch.cuda.is_available() else "cpu")

# Softmax 策略
def softmax_policy(qvals, temp=0.9):
    soft = torch.exp(qvals / temp) / torch.sum(torch.exp(qvals / temp))
    action = torch.multinomial(soft, 1)
    return action

# 獲取格點在二維格子中的座標
def get_coords(grid, j):
    x = int(np.floor(j / grid.shape[0]))
    y = int(j - x * grid.shape[0])
    return x, y

# 計算二維格子的獎勵
def get_reward_2d(action, action_mean):
    r = (action * (action_mean - action / 2)).sum() / action.sum()
    return torch.tanh(5 * r)

# 獲取自旋格點的子狀態
def get_substate(b):
    s = torch.zeros(2).to(device)
    if b > 0:
        s[1] = 1
    else:
        s[0] = 1
    return s

# 計算格點的平均行動
def mean_action(grid, j):
    x, y = get_coords(grid, j)
    action_mean = torch.zeros(2).to(device)
    for i in [-1, 0, 1]:
        for k in [-1, 0, 1]:
            if i == k == 0:
                continue
            x_, y_ = x + i, y + k
            x_ = x_ if x_ >= 0 else grid.shape[0] - 1
            y_ = y_ if y_ >= 0 else grid.shape[1] - 1
            x_ = x_ if x_ < grid.shape[0] else 0
            y_ = y_ if y_ < grid.shape[1] else 0
            cur_n = grid[x_, y_]
            s = get_substate(cur_n)
            action_mean += s
    action_mean /= action_mean.sum()
    return action_mean

# 初始化自旋格子
def init_grid(size=(10, 10)):
    grid = torch.randint(0, 2, size).to(device)
    return grid

# 生成參數的函數
def gen_params(N, size):
    ret = []
    for i in range(N):
        vec = torch.randn(size).to(device) / 10.
        vec.requires_grad = True
        ret.append(vec)
    return ret

# 假設的 q 函數，用來估算 Q 值
def qfunc(state, params, layers=[(2, 10), (10, 2)], afn=torch.tanh):
    l1n = layers[0]
    l1s = np.prod(l1n)
    theta_1 = params[0:l1s].reshape(l1n)
    l2n = layers[1]
    l2s = np.prod(l2n)
    theta_2 = params[l1s:l2s + l1s].reshape(l2n)
    bias = torch.ones((1, theta_1.shape[1])).to(device)
    l1 = state @ theta_1 + bias
    l1 = torch.nn.functional.elu(l1)
    l2 = afn(l1 @ theta_2)
    return l2.flatten()

# 主程式開始
size = (10, 10)
J = np.prod(size)  # 計算格子的總數量
hid_layer = 10
layers = [(2, hid_layer), (hid_layer, 2)]
params = gen_params(1, 2 * hid_layer + hid_layer * 2)  # 生成隨機參數
grid = init_grid(size=size)
grid_ = grid.clone()
grid__ = grid.clone()

# 顯示初始狀態並保存為圖片
plt.imshow(grid.cpu())
plt.title("Initial Spin Lattice")
plt.colorbar()
plt.savefig("initial_grid.png")
plt.show()
print("Initial grid sum:", grid.sum())

# 開始訓練過程
epochs = 75
lr = 0.0001
num_iter = 3  # A
losses = [[] for i in range(J)]  # B
replay_size = 50  # C
replay = deque(maxlen=replay_size)  # D
batch_size = 10  # E
gamma = 0.9  # F

for i in range(epochs): 
    act_means = torch.zeros((J, 2)).to(device)  # G
    q_next = torch.zeros(J).to(device)  # H
    for m in range(num_iter):  # I
        for j in range(J):  # J
            action_mean = mean_action(grid_, j).detach()
            act_means[j] = action_mean.clone()
            qvals = qfunc(action_mean.detach(), params[0], layers=layers)
            action = softmax_policy(qvals.detach(), temp=0.5)
            grid__[get_coords(grid_, j)] = action
            q_next[j] = torch.max(qvals).detach()
        grid_.data = grid__.data
    grid.data = grid_.data

    # 計算當前動作和獎勵
    actions = torch.stack([get_substate(a.item()) for a in grid.flatten()])
    rewards = torch.stack([get_reward_2d(actions[j], act_means[j]) for j in range(J)])
    exp = (actions, rewards, act_means, q_next)  # K
    replay.append(exp)
    shuffle(replay)

    # 使用批量回放進行更新
    if len(replay) > batch_size:  # L
        ids = np.random.randint(low=0, high=len(replay), size=batch_size)  # M
        exps = [replay[idx] for idx in ids]
        for j in range(J):
            jacts = torch.stack([ex[0][j] for ex in exps]).detach()
            jrewards = torch.stack([ex[1][j] for ex in exps]).detach()
            jmeans = torch.stack([ex[2][j] for ex in exps]).detach()
            vs = torch.stack([ex[3][j] for ex in exps]).detach()
            qvals = torch.stack([qfunc(jmeans[h].detach(), params[0], layers=layers) for h in range(batch_size)])
            target = qvals.clone().detach()
            target[:, torch.argmax(jacts, dim=1)] = jrewards + gamma * vs
            loss = torch.sum(torch.pow(qvals - target.detach(), 2))
            losses[j].append(loss.item())
            loss.backward()
            with torch.no_grad():
                params[0] = params[0] - lr * params[0].grad
            params[0].requires_grad = True

# 可視化損失值和最終格子狀態
fig, ax = plt.subplots(2, 1)
fig.set_size_inches(10, 10)

# 損失值的平均變化
mean_losses = np.mean([np.array(l) for l in losses if len(l) > 0], axis=0)
ax[0].plot(mean_losses)
ax[0].set_title("Average Loss over Epochs")

# 顯示最終自旋格子狀態
ax[1].imshow(grid.cpu())
ax[1].set_title("Final Spin Lattice")
plt.savefig('training_results.png')  # 保存圖片
plt.show()

# Softmax 策略

def softmax_policy(qvals, temp=0.9):
    soft = torch.exp(qvals / temp) / torch.sum(torch.exp(qvals / temp))
    action = torch.multinomial(soft, 1)
    return action

# Softmax 函數公式

$$
P(a_i) = \frac{\exp(Q(a_i) / \text{temp})}{\sum_j \exp(Q(a_j) / \text{temp})}
$$

- **$Q(a_i)$**：行為 $a_i$ 的 Q 值。
- **$\text{temp}$**：溫度參數，用於調節 Softmax 函數對 Q 值的響應程度。

---

## 溫度參數對指數函數的影響

### 1. 高溫（大 $\text{temp}$）

當 $\text{temp} \to \infty$ 時：

$$
\frac{Q(a_i)}{\text{temp}} \to 0
$$

指數函數趨近於 1：

$$
\exp\left(\frac{Q(a_i)}{\text{temp}}\right) \approx 1
$$

因此所有 Q 值的影響力變得相似，機率分布接近均勻分布：

$$
P(a_i) \approx \frac{1}{\text{總行為數}}
$$

- **原因**：當 $\text{temp}$ 很大時，Q 值被過度平滑化，無法有效區分高 Q 值和低 Q 值行為，最終所有行為的機率幾乎相等。

---

### 2. 低溫（小 $\text{temp}$）

當 $\text{temp} \to 0^+$ 時：

$$
\frac{Q(a_i)}{\text{temp}} \to \infty \quad (\text{若 } Q(a_i) \text{ 最大})
$$

此時，Q 值較大的行為對指數函數的影響變得極其顯著，而其他較小 Q 值的行為影響趨近於零。因此，Softmax 分布會強烈偏向 Q 值最大的行為。

- **數學解釋**：

若 $Q(a_k) > Q(a_j)$ 且 $\text{temp} \to 0^+$，則：

$$
\frac{\exp(Q(a_k) / \text{temp})}{\exp(Q(a_j) / \text{temp})} \to \infty
$$

導致機率分布中：

$$
P(a_k) \to 1 \quad \text{且} \quad P(a_j) \to 0 \quad (\text{若 } j \neq k)
$$

# Softmax 計算範例：高溫與低溫

## 1. 高溫：$$\text{temp} = 2$$

### 指數值計算
$$
\exp\left(\frac{Q}{\text{temp}}\right) = \exp\left(\frac{[2, 1, 0]}{2}\right) = [e^1, e^{0.5}, e^0] \approx [2.718, 1.649, 1.0]
$$

### 分母計算（正規化因子）
$$
Z = \sum \exp\left(\frac{Q}{\text{temp}}\right) = 2.718 + 1.649 + 1.0 = 5.367
$$

### Softmax 分布
$$
\text{soft} = \left[\frac{2.718}{5.367}, \frac{1.649}{5.367}, \frac{1.0}{5.367}\right] \approx [0.506, 0.307, 0.186]
$$

### 解釋
- 在高溫下，所有行為的機率相對均勻。
- 行為 $0$ 的主導性被削弱，低 Q 值行為（行為 $2$）仍有一定機率被選中。

---

## 2. 低溫：$$\text{temp} = 0.1$$

### 指數值計算
$$
\exp\left(\frac{Q}{\text{temp}}\right) = \exp\left(\frac{[2, 1, 0]}{0.1}\right) = [e^{20}, e^{10}, e^0]
$$

由於指數函數的快速增長特性：
$$
\exp(20) \gg \exp(10) \gg \exp(0)
$$

### 分母計算（正規化因子）
$$
Z \approx \exp(20) + \exp(10) + 1
$$

### Softmax 分布
$$
\frac{\exp(20)}{Z} \approx 1.0, \quad \frac{\exp(10)}{Z} \approx 0, \quad \frac{1}{Z} \approx 0
$$

### 解釋
- 在低溫下，Softmax 分布強烈偏向於 Q 值最大的行為。
- 行為 $0$ 的機率接近 $1.0$，而其他行為的機率接近 $0$。
# `get_reward_2d` 函數

此函數用於計算基於行動和行動均值的獎勵，並應用雙曲正切函數（`tanh`）對結果進行縮放。

## 函數實現

def get_reward_2d(action, action_mean):
    r = (action * (action_mean - action / 2)).sum() / action.sum()
    return torch.tanh(5 * r)

# 格點行為獎勵計算範例

## 1. 假設數據

### 1.1 格點的行為分布
格點 $j$ 的行為概率（action）：
$$
\text{action} = [0.7, 0.3]
$$
- $0.7$：格點 $j$ 以 70% 的概率選擇行為 $0$。
- $0.3$：格點 $j$ 以 30% 的概率選擇行為 $1$。

---

### 1.2 鄰居的平均行為偏好
鄰居行為的平均分布（action_mean）：
$$
\text{action_mean} = [0.6, 0.4]
$$
- $0.6$：鄰居偏好行為 $0$。
- $0.4$：鄰居偏好行為 $1$。

---

## 2. 運算過程

### 2.1 相容性項計算
點乘的第一部分計算當前行為與鄰居行為的匹配度：
$$
\text{compatibility} = \text{action} \cdot \text{action_mean}
$$
對應計算：
$$
\text{compatibility} = (0.7 \times 0.6) + (0.3 \times 0.4) = 0.42 + 0.12 = 0.54
$$
- **解釋**：
- $0.42$：當前行為 $0$ 與鄰居偏好 $0$ 的匹配程度。
- $0.12$：當前行為 $1$ 與鄰居偏好 $1$ 的匹配程度。
- **總和** $0.54$：行為的總相容性，越高表示與鄰居越一致。

---

### 2.2 自抑制項計算
點乘的第二部分計算行為的極端性懲罰：
$$
\text{self-suppression} = \text{action} \cdot \frac{\text{action}}{2}
$$
對應計算：
$$
\text{self-suppression} = (0.7 \times 0.35) + (0.3 \times 0.15) = 0.245 + 0.045 = 0.29
$$
- **解釋**：
- $0.245$：行為 $0$ 的極端性懲罰。
- $0.045$：行為 $1$ 的極端性懲罰。
- **總和** $0.29$：自抑制懲罰，越高表示行為分布越極端。

---

### 2.3 最終點乘結果
將兩部分結合：
$$
r = \text{compatibility} - \text{self-suppression}
$$
計算：
$$
r = 0.54 - 0.29 = 0.25
$$

---

### 2.4 正規化與非線性壓縮
最終獎勵值進行非線性處理：
$$
\text{reward} = \tanh(5 \cdot r)
$$
對應計算：
$$
\text{reward} = \tanh(5 \cdot 0.25) = \tanh(1.25) \approx 0.848
$$

---

## 3. 結果解釋

- **相容性高**：當前行為 $\text{action} = [0.7, 0.3]$ 與鄰居行為 $\text{action_mean} = [0.6, 0.4]$ 有較高的匹配度（$0.54$）。
- **行為懲罰**：$0.7$ 和 $0.3$ 的分布較平衡，因此自抑制項 $0.29$ 不算太高。
- **最終獎勵**：經過壓縮，獎勵值為 $\approx 0.848$，說明該行為與鄰居偏好相對一致，獲得了較高的激勵。

---

## 4. 不同情況比較

### 4.1 完全一致的行為
若當前行為與鄰居偏好完全一致：
$$
\text{action} = [0.6, 0.4], \quad \text{action_mean} = [0.6, 0.4]
$$
計算：
$$
\text{compatibility} = (0.6 \times 0.6) + (0.4 \times 0.4) = 0.36 + 0.16 = 0.52
$$
$$
\text{self-suppression} = (0.6 \times 0.3) + (0.4 \times 0.2) = 0.18 + 0.08 = 0.26
$$
$$
r = 0.52 - 0.26 = 0.26
$$
$$
\text{reward} = \tanh(5 \cdot 0.26) \approx 0.856
$$
- **結果**：行為完全一致時，獎勵值達到更高，為 $\approx 0.856$。

---

### 4.2 完全不一致的行為
若當前行為與鄰居偏好完全不一致：
$$
\text{action} = [0.4, 0.6], \quad \text{action_mean} = [0.6, 0.4]
$$
計算：
$$
\text{compatibility} = (0.4 \times 0.6) + (0.6 \times 0.4) = 0.24 + 0.24 = 0.48
$$
$$
\text{self-suppression} = (0.4 \times 0.2) + (0.6 \times 0.3) = 0.08 + 0.18 = 0.26
$$
$$
r = 0.48 - 0.26 = 0.22
$$
$$
\text{reward} = \tanh(5 \cdot 0.22) \approx 0.804
$$
- **結果**：行為不一致時，獎勵值降低，但仍有一定的正值。

---

### 4.3 極端行為
若行為極端化，例如完全偏向一個行為：
$$
\text{action} = [1.0, 0.0], \quad \text{action_mean} = [0.6, 0.4]
$$
計算：
$$
\text{compatibility} = (1.0 \times 0.6) + (0.0 \times 0.4) = 0.6 + 0 = 0.6
$$
$$
\text{self-suppression} = (1.0 \times 0.5) + (0.0 \times 0) = 0.5 + 0 = 0.5
$$
$$
r = 0.6 - 0.5 = 0.1
$$
$$
\text{reward} = \tanh(5 \cdot 0.1) \approx 0.462
$$
- **結果**：行為極端化時，自抑制項懲罰變高，獎勵值顯著降低。

---

- **高相容性**：行為與鄰居偏好一致時，獲得更高的獎勵。
- **平衡行為**：自抑制項防止極端行為，保證行為分布更平衡。
- **極端行為懲罰**：行為極端化會導致懲罰增加，獎勵值顯著下降。
- **獎勵設計特點**：
- 鼓勵行為與鄰居一致（高相容性）。
- 懲罰極端行為，促進穩定性與協同效應。

# `Get_substate`

def get_substate(b):
    s = torch.zeros(2).to(device)
    if b > 0:
        s[1] = 1
    else:
        s[0] = 1
    return s

### 自旋狀態轉換函數：`get_substate`

該函數用於將輸入值 $b$ 轉換為二維向量，用以表示二元自旋系統的狀態。

---

### **數學描述**

若 $b > 0$，則輸出：
$$
[0, 1]
$$
表示「向上自旋」。

若 $b \leq 0$，則輸出：
$$
[1, 0]
$$
表示「向下自旋」。

---

### **程式碼實現**

def get_substate(b):
    s = torch.zeros(2).to(device)
    if b > 0:
        s[1] = 1
    else:
        s[0] = 1
    return s

### 範例

# 假設 b 為正數
b = 1
state = get_substate(b)
print(state)  # 輸出: tensor([0., 1.], device='cuda:0')

# 假設 b 為負數
b = -1
state = get_substate(b)
print(state)  # 輸出: tensor([1., 0.], device='cuda:0')

# `Mean action`

def mean_action(grid, j):
    x, y = get_coords(grid, j)
    action_mean = torch.zeros(2).to(device)
    for i in [-1, 0, 1]:
        for k in [-1, 0, 1]:
            if i == k == 0:
                continue
            x_, y_ = x + i, y + k
            # 處理格點的環繞邊界
            x_ = x_ if x_ >= 0 else grid.shape[0] - 1
            y_ = y_ if y_ >= 0 else grid.shape[1] - 1
            x_ = x_ if x_ < grid.shape[0] else 0
            y_ = y_ if y_ < grid.shape[1] else 0
            cur_n = grid[x_, y_]
            s = get_substate(cur_n)
            action_mean += s
    action_mean /= action_mean.sum()
    return action_mean

# 範例：計算目標格點的平均行動

## 自旋格子設定
假設一個 $5 \times 5$ 的自旋格子，定義如下：

grid = torch.tensor([
    [0, 1, 0, 1, 0],
    [1, 0, 1, 0, 1],
    [0, 1, 0, 1, 0],
    [1, 0, 1, 0, 1],
    [0, 1, 0, 1, 0],
]).to(device)

## 目標格點平均行動計算範例

目標格點索引 $j=12$，對應的座標為 $(2, 2)$。

### 計算步驟

#### 1. 鄰居格點的自旋狀態
目標格點 $(2, 2)$ 的鄰居狀態為：

$$
\{0, 1, 0, 1, 1, 0, 1, 0\}
$$

---

#### 2. 轉換為向量形式
將鄰居的自旋狀態轉換為向量形式：

$$
\{
[1, 0], [0, 1], [1, 0], [0, 1], [0, 1], [1, 0], [0, 1], [1, 0]
\}
$$

---

#### 3. 累加結果
累加所有向量：

$$
\text{action_mean} = [5, 3]
$$

---

#### 4. 歸一化
將累加結果歸一化：

$$
\text{action_mean} = \left[\frac{5}{8}, \frac{3}{8}\right] = [0.625, 0.375]
$$

---

### 返回結果
目標格點的平均行動為：

$$
[0.625, 0.375]
$$

## 網路結構特點

### 1. 兩層全連接結構
- **第一層**：
- 輸入：$(2,)$
- 隱藏層輸出：$(10,)$
- **第二層**：
- 隱藏層輸入：$(10,)$
- 最終輸出：$(2,)$

---

### 2. 激活函數
- **ELU**：應用於第一層，增強非線性特性。
- **Tanh**：應用於第二層，進一步提升網路的非線性能力。

---

### 3. 參數共享
- **特點**：所有格點共享網路參數。
- **優勢**：大幅簡化模型設計，提升計算效率。
## 實驗結果