從 Mendeley Food FTIR spectra 開始學 2T2D-COS

資料來源是 Data for: Spectra Data Classification with Kernel Extreme Learning。它不是專門為 2T2D-COS 建立的資料集，但因為每個樣本都有完整 FTIR spectrum 與類別標籤，非常適合用來教「從一維光譜生成 2D correlation image」。

以下圖形直接由 Mendeley coffee FTIR CSV 重新計算產生，不是示意圖。可作為參考答案，回到 Python / Orange 重算即可確認是否得到相同趨勢。

下面程式下載 Mendeley zip、讀取 coffee FTIR CSV、整理成 X: samples × wavenumbers，做 SNV 前處理，再用與 SpectraView 一致的 Noda 公式算真 2T2D 與廣義 2D-COS（不是會讓異步退化為 0 的減平均版）。

import io, zipfile, requests
import numpy as np
import pandas as pd
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.model_selection import StratifiedKFold, cross_val_score
from sklearn.svm import SVC

url = "https://data.mendeley.com/public-api/zip/frrv2yd9rg/download/1"
r = requests.get(url, headers={"User-Agent": "Mozilla/5.0"})
r.raise_for_status()
z = zipfile.ZipFile(io.BytesIO(r.content))

coffee_file = [f for f in z.namelist() if "instant_coffee" in f and f.endswith(".csv")][0]
raw = pd.read_csv(z.open(coffee_file), header=None)

sample_id = raw.iloc[0, 1:].astype(str).to_numpy()
group_code = raw.iloc[1, 1:].astype(str).to_numpy()
y = raw.iloc[2, 1:].astype(str).to_numpy()      # Arabica / Robusta
wavenumber = pd.to_numeric(raw.iloc[3:, 0]).to_numpy()
X = raw.iloc[3:, 1:].T.astype(float).to_numpy() # samples × wavenumbers

print(X.shape, y[:5], wavenumber.min(), wavenumber.max())

Xz = StandardScaler().fit_transform(X)
pca = PCA(n_components=2).fit_transform(Xz)

clf = SVC(kernel="linear")
cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=7)
scores = cross_val_score(clf, Xz, y, cv=cv)
print("CV accuracy:", scores.mean(), scores.std())

# Noda 2018 two-trace 2D correlation.
# 不要對兩條光譜減平均：只有兩條時減平均會得到 ya-yb 與 yb-ya，
# 使異步圖恆為 0、同步圖塌成 1D 差異光譜的外積（rank-1）。
# 另外避免變數名用 async（Python 3.7+ 是保留字）。
def two_trace_2dcos(ya, yb):
    sync = 0.5 * (np.outer(ya, ya) + np.outer(yb, yb))
    asyn = 0.5 * (np.outer(ya, yb) - np.outer(yb, ya))
    return sync, asyn

# SNV 前處理，再算兩類平均的真 2T2D
Xs = (X - X.mean(1, keepdims=True)) / X.std(1, keepdims=True)
arab, robu = Xs[y == "Arabica"].mean(0), Xs[y == "Robusta"].mean(0)
sync, asyn = two_trace_2dcos(arab, robu)
ratio = np.abs(asyn).max() / np.abs(sync).max()
print(sync.shape, "max|async|/max|sync| =", round(ratio, 3))  # ~0.15 非退化

# 廣義 2D-COS：給一條「隨擾動變化」的光譜序列 M (m_spectra × n_wavenumbers)。
def noda_matrix(m):
    N = np.zeros((m, m))
    for j in range(m):
        for k in range(m):
            if j != k:
                N[j, k] = 1.0 / (np.pi * (k - j))
    return N

def generalized_2dcos(M, reference="mean"):
    D = M - M.mean(axis=0) if reference == "mean" else M
    m = D.shape[0]
    sync = D.T @ D / (m - 1)
    asyn = D.T @ (noda_matrix(m) @ D) / (m - 1)
    return sync, asyn

# 把不同烘焙時間 / 摻配比例的光譜疊成 M，再 sync, asyn = generalized_2dcos(snv(M))。
# Arabica/Robusta 分類沒有這種擾動軸，因此只能做兩跡 2T2D；
# 要展現 band-to-band covariation 與先後順序，請改用真正的擾動序列。

除了 Python 版 2T2D-COS，Orange Data Mining 的 Spectroscopy add-on 也有 2D Correlation Plot widget。建議先用 Orange 做互動式探索，看懂「一維光譜如何變成二維相關圖」；再用 Python 重現流程、輸出圖檔與建立機器學習特徵。