import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize

class SVM:
    """SVMのスクラッチ実装（二次計画法ベース）"""
    def __init__(self, C=1.0, kernel='linear', gamma=1.0):
        self.C = C
        self.kernel_type = kernel
        self.gamma = gamma
        self.alpha = None
        self.b = None
        self.support_vectors = None
        self.support_vector_labels = None
        self.support_vector_alphas = None

    def _kernel(self, x1, x2):
        """カーネル関数"""
        if self.kernel_type == 'linear':
            return np.dot(x1, x2)
        elif self.kernel_type == 'rbf':
            return np.exp(-self.gamma * np.linalg.norm(x1 - x2) ** 2)
        elif self.kernel_type == 'polynomial':
            return (1 + np.dot(x1, x2)) ** 3

    def _kernel_matrix(self, X):
        """グラム行列の計算"""
        N = X.shape[0]
        K = np.zeros((N, N))
        for i in range(N):
            for j in range(N):
                K[i, j] = self._kernel(X[i], X[j])
        return K

    def fit(self, X, y):
        """学習"""
        N = X.shape[0]
        K = self._kernel_matrix(X)

        # 双対問題を最小化問題に変換: min -Σα_i + (1/2)ΣΣα_i α_j y_i y_j K_ij
        # scipy.optimize.minimize を使用
        def objective(alpha):
            return 0.5 * np.sum((alpha * y)[:, None] * (alpha * y)[None, :] * K) - np.sum(alpha)

        def gradient(alpha):
            return np.sum((alpha * y)[None, :] * K * y[:, None], axis=1) - 1

        # 制約: Σα_i y_i = 0
        constraints = {'type': 'eq', 'fun': lambda a: np.dot(a, y)}
        # 境界: 0 <= α_i <= C
        bounds = [(0, self.C) for _ in range(N)]

        # 最適化
        alpha0 = np.zeros(N)
        result = minimize(objective, alpha0, jac=gradient,
                         bounds=bounds, constraints=constraints,
                         method='SLSQP', options={'maxiter': 1000, 'ftol': 1e-10})
        self.alpha = result.x

        # サポートベクターの抽出（α > 閾値）
        sv_threshold = 1e-5
        sv_mask = self.alpha > sv_threshold
        self.support_vectors = X[sv_mask]
        self.support_vector_labels = y[sv_mask]
        self.support_vector_alphas = self.alpha[sv_mask]

        # バイアス項の計算（freeサポートベクターを使用）
        free_sv_mask = (self.alpha > sv_threshold) & (self.alpha < self.C - sv_threshold)
        if np.any(free_sv_mask):
            free_indices = np.where(free_sv_mask)[0]
            b_values = []
            for idx in free_indices:
                b_val = y[idx]
                for j in range(len(self.support_vectors)):
                    b_val -= self.support_vector_alphas[j] * self.support_vector_labels[j] * \
                             self._kernel(self.support_vectors[j], X[idx])
                b_values.append(b_val)
            self.b = np.mean(b_values)
        else:
            # 全てのサポートベクターで平均
            sv_indices = np.where(sv_mask)[0]
            b_values = []
            for idx in sv_indices:
                b_val = y[idx]
                for j in range(len(self.support_vectors)):
                    b_val -= self.support_vector_alphas[j] * self.support_vector_labels[j] * \
                             self._kernel(self.support_vectors[j], X[idx])
                b_values.append(b_val)
            self.b = np.mean(b_values)

    def decision_function(self, X):
        """決定関数 f(x) = Σα_i y_i k(x_i, x) + b"""
        N = X.shape[0]
        f = np.zeros(N)
        for i in range(N):
            for j in range(len(self.support_vectors)):
                f[i] += self.support_vector_alphas[j] * self.support_vector_labels[j] * \
                         self._kernel(self.support_vectors[j], X[i])
            f[i] += self.b
        return f

    def predict(self, X):
        """予測"""
        return np.sign(self.decision_function(X))

# 線形分離可能なデータの生成
np.random.seed(42)
N = 40
X_pos = np.random.randn(N // 2, 2) + np.array([2, 2])
X_neg = np.random.randn(N // 2, 2) + np.array([-2, -2])
X = np.vstack([X_pos, X_neg])
y = np.array([1] * (N // 2) + [-1] * (N // 2), dtype=float)

# ハードマージンSVM（C を大きくしてハードマージンを近似）
svm_hard = SVM(C=1000.0, kernel='linear')
svm_hard.fit(X, y)

# 決定境界の可視化
fig, ax = plt.subplots(figsize=(8, 8))

# 決定境界のメッシュ
x_min, x_max = X[:, 0].min() - 1, X[:, 0].max() + 1
y_min, y_max = X[:, 1].min() - 1, X[:, 1].max() + 1
xx, yy = np.meshgrid(np.linspace(x_min, x_max, 200),
                      np.linspace(y_min, y_max, 200))
X_grid = np.c_[xx.ravel(), yy.ravel()]
Z = svm_hard.decision_function(X_grid).reshape(xx.shape)

# 決定境界とマージン
ax.contour(xx, yy, Z, levels=[-1, 0, 1], colors=['blue', 'black', 'red'],
           linestyles=['--', '-', '--'], linewidths=[1.5, 2, 1.5])
ax.contourf(xx, yy, Z, levels=[-100, 0, 100], alpha=0.1, colors=['blue', 'red'])

# データ点
ax.scatter(X[y == 1, 0], X[y == 1, 1], c='red', edgecolors='k', s=50, label='$y=+1$', zorder=5)
ax.scatter(X[y == -1, 0], X[y == -1, 1], c='blue', edgecolors='k', s=50, label='$y=-1$', zorder=5)

# サポートベクターの強調
ax.scatter(svm_hard.support_vectors[:, 0], svm_hard.support_vectors[:, 1],
           s=200, facecolors='none', edgecolors='green', linewidths=2, label='Support vectors', zorder=6)

ax.set_xlabel('$x_1$', fontsize=12)
ax.set_ylabel('$x_2$', fontsize=12)
ax.set_title(f'Hard Margin SVM (# support vectors = {len(svm_hard.support_vectors)})', fontsize=14)
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('svm_hard_margin.png', dpi=150, bbox_inches='tight')
plt.show()

# ノイズを含むデータの生成
np.random.seed(42)
N = 100
X_pos = np.random.randn(N // 2, 2) + np.array([1, 1])
X_neg = np.random.randn(N // 2, 2) + np.array([-1, -1])
X = np.vstack([X_pos, X_neg])
y = np.array([1] * (N // 2) + [-1] * (N // 2), dtype=float)

fig, axes = plt.subplots(1, 3, figsize=(18, 5))
C_values = [0.1, 1.0, 100.0]

for ax, C in zip(axes, C_values):
    svm = SVM(C=C, kernel='linear')
    svm.fit(X, y)

    # 決定境界
    x_min, x_max = X[:, 0].min() - 1.5, X[:, 0].max() + 1.5
    y_min, y_max = X[:, 1].min() - 1.5, X[:, 1].max() + 1.5
    xx, yy = np.meshgrid(np.linspace(x_min, x_max, 150),
                          np.linspace(y_min, y_max, 150))
    X_grid = np.c_[xx.ravel(), yy.ravel()]
    Z = svm.decision_function(X_grid).reshape(xx.shape)

    ax.contour(xx, yy, Z, levels=[-1, 0, 1], colors=['blue', 'black', 'red'],
               linestyles=['--', '-', '--'], linewidths=[1, 2, 1])
    ax.contourf(xx, yy, Z, levels=[-100, 0, 100], alpha=0.1, colors=['blue', 'red'])

    ax.scatter(X[y == 1, 0], X[y == 1, 1], c='red', edgecolors='k', s=30, zorder=5)
    ax.scatter(X[y == -1, 0], X[y == -1, 1], c='blue', edgecolors='k', s=30, zorder=5)
    ax.scatter(svm.support_vectors[:, 0], svm.support_vectors[:, 1],
               s=120, facecolors='none', edgecolors='green', linewidths=2, zorder=6)

    ax.set_title(f'C = {C} (# SV = {len(svm.support_vectors)})', fontsize=13)
    ax.set_xlim(x_min, x_max)
    ax.set_ylim(y_min, y_max)
    ax.grid(True, alpha=0.3)

plt.suptitle('Soft Margin SVM: Effect of C Parameter', fontsize=15, y=1.02)
plt.tight_layout()
plt.savefig('svm_soft_margin_C.png', dpi=150, bbox_inches='tight')
plt.show()

# 非線形データの生成（月型データ）
from sklearn.datasets import make_moons

np.random.seed(42)
X, y_01 = make_moons(n_samples=100, noise=0.15, random_state=42)
y = 2 * y_01 - 1  # {0,1} -> {-1,+1}

fig, axes = plt.subplots(1, 3, figsize=(18, 5))

kernel_configs = [
    ('linear', 1.0, 'Linear Kernel'),
    ('rbf', 0.5, 'RBF Kernel ($\\gamma=0.5$)'),
    ('rbf', 5.0, 'RBF Kernel ($\\gamma=5.0$)'),
]

for ax, (kernel, gamma, title) in zip(axes, kernel_configs):
    svm = SVM(C=10.0, kernel=kernel, gamma=gamma)
    svm.fit(X, y)

    # 決定境界
    x_min, x_max = X[:, 0].min() - 0.5, X[:, 0].max() + 0.5
    y_min, y_max = X[:, 1].min() - 0.5, X[:, 1].max() + 0.5
    xx, yy = np.meshgrid(np.linspace(x_min, x_max, 150),
                          np.linspace(y_min, y_max, 150))
    X_grid = np.c_[xx.ravel(), yy.ravel()]
    Z = svm.decision_function(X_grid).reshape(xx.shape)

    ax.contourf(xx, yy, Z, levels=[-100, 0, 100], alpha=0.2, colors=['blue', 'red'])
    ax.contour(xx, yy, Z, levels=[0], colors='black', linewidths=2)

    ax.scatter(X[y == 1, 0], X[y == 1, 1], c='red', edgecolors='k', s=30, zorder=5)
    ax.scatter(X[y == -1, 0], X[y == -1, 1], c='blue', edgecolors='k', s=30, zorder=5)
    ax.scatter(svm.support_vectors[:, 0], svm.support_vectors[:, 1],
               s=120, facecolors='none', edgecolors='green', linewidths=2, zorder=6)

    ax.set_title(f'{title}\n(# SV = {len(svm.support_vectors)})', fontsize=12)
    ax.grid(True, alpha=0.3)

plt.suptitle('Kernel SVM: Nonlinear Classification', fontsize=15, y=1.05)
plt.tight_layout()
plt.savefig('svm_kernel.png', dpi=150, bbox_inches='tight')
plt.show()

# ヒンジ損失と他の損失関数の比較
fig, ax = plt.subplots(figsize=(10, 6))
yf = np.linspace(-2, 3, 500)

# ヒンジ損失
hinge = np.maximum(0, 1 - yf)
# 指数損失（AdaBoost）
exp_loss = np.exp(-yf)
# 対数損失（ロジスティック回帰）
log_loss = np.log(1 + np.exp(-yf)) / np.log(2)
# 0-1損失
zero_one = (yf < 0).astype(float)

ax.plot(yf, hinge, 'b-', linewidth=2.5, label='Hinge loss: $\\max(0, 1-yf)$')
ax.plot(yf, exp_loss, 'r-', linewidth=2, label='Exponential loss: $e^{-yf}$')
ax.plot(yf, log_loss, 'g-', linewidth=2, label='Log loss: $\\log_2(1+e^{-yf})$')
ax.plot(yf, zero_one, 'k--', linewidth=2, label='0-1 loss')

ax.set_xlabel('$yf(x)$ (margin)', fontsize=13)
ax.set_ylabel('Loss', fontsize=13)
ax.set_title('Comparison of Loss Functions for Classification', fontsize=14)
ax.legend(fontsize=11, loc='upper right')
ax.set_xlim(-2, 3)
ax.set_ylim(-0.1, 4)
ax.grid(True, alpha=0.3)
ax.axvline(x=0, color='gray', linestyle=':', alpha=0.5)
ax.axvline(x=1, color='gray', linestyle=':', alpha=0.5)

plt.tight_layout()
plt.savefig('loss_functions_comparison.png', dpi=150, bbox_inches='tight')
plt.show()