Exercise 3 — MLP

"""MLP escrito do zero em NumPy, usado nos Exercises 2, 3 e 4.

O mesmo código de forward, backward e update serve para classificação binária
(saída sigmoide + entropia cruzada binária) e multiclasse (saída softmax +
entropia cruzada categórica). Só mudam os argumentos: tamanhos das camadas,
tipo de saída, taxa de aprendizado e gerador aleatório.

Notação da aula: em cada camada, u = W a + b, com W no formato
(neurônios que recebem, neurônios que enviam).
"""
import numpy as np


class MLP:
    def __init__(self, layer_sizes, output="sigmoid", lr=0.1, rng=None):
        """layer_sizes: lista com o tamanho de cada camada, ex. [2, 8, 1].
        output: "sigmoid" (binário, 1 neurônio de saída) ou "softmax" (multiclasse).
        lr: taxa de aprendizado do gradiente descendente.
        rng: np.random.Generator usado para inicializar os pesos.
        """
        self.layer_sizes = layer_sizes
        self.output = output
        self.lr = lr
        # Pesos com desvio 1/sqrt(n_entradas) (inicialização de Xavier, boa para tanh); vieses em zero
        self.W = [rng.normal(0.0, np.sqrt(1.0 / n_in), size=(n_out, n_in))
                  for n_in, n_out in zip(layer_sizes[:-1], layer_sizes[1:])]
        self.b = [np.zeros(n_out) for n_out in layer_sizes[1:]]

    # ------------------------------------------------------------ ativações
    @staticmethod
    def tanh(u):
        return np.tanh(u)

    @staticmethod
    def sigmoid(u):
        return 1.0 / (1.0 + np.exp(-np.clip(u, -500, 500)))

    @staticmethod
    def softmax(u):
        e = np.exp(u - u.max(axis=1, keepdims=True))   # subtrai o máximo para não estourar
        return e / e.sum(axis=1, keepdims=True)

    # ------------------------------------------------------------ forward
    def forward(self, X):
        """X: (N, n_entradas). Guarda a saída de cada camada para usar no backward."""
        self.activations = [X]
        for i, (W, b) in enumerate(zip(self.W, self.b)):
            u = self.activations[-1] @ W.T + b          # pré-ativação da camada
            if i < len(self.W) - 1:
                a = self.tanh(u)                        # camadas ocultas: tanh
            elif self.output == "sigmoid":
                a = self.sigmoid(u)
            else:
                a = self.softmax(u)
            self.activations.append(a)
        return self.activations[-1]

    # ------------------------------------------------------------ perda
    def loss(self, Y_hat, Y):
        """Entropia cruzada média. Y tem o mesmo formato de Y_hat (one-hot no caso multiclasse)."""
        eps = 1e-12
        if self.output == "sigmoid":
            return float(-np.mean(Y * np.log(Y_hat + eps) + (1 - Y) * np.log(1 - Y_hat + eps)))
        return float(-np.mean(np.sum(Y * np.log(Y_hat + eps), axis=1)))

    # ------------------------------------------------------------ backward
    def backward(self, Y):
        """Calcula os gradientes de todos os pesos e vieses (regra da cadeia, sem autodiff).

        Tanto sigmoide + BCE quanto softmax + entropia cruzada categórica dão
        dL/du = (Y_hat - Y) / N na saída, então o mesmo código serve aos dois casos.
        """
        N = Y.shape[0]
        delta = (self.activations[-1] - Y) / N          # dL/du da camada de saída
        self.grad_W = [None] * len(self.W)
        self.grad_b = [None] * len(self.b)
        for i in reversed(range(len(self.W))):
            self.grad_W[i] = delta.T @ self.activations[i]   # dL/dW = delta^T a_anterior
            self.grad_b[i] = delta.sum(axis=0)               # dL/db = soma dos deltas
            if i > 0:
                # volta uma camada: dL/da = delta W, e multiplica pela derivada da tanh (1 - a^2)
                delta = (delta @ self.W[i]) * (1 - self.activations[i] ** 2)

    # ------------------------------------------------------------ update
    def update(self):
        """Gradiente descendente: p <- p - lr * dL/dp."""
        for i in range(len(self.W)):
            self.W[i] -= self.lr * self.grad_W[i]
            self.b[i] -= self.lr * self.grad_b[i]

    # ------------------------------------------------------------ treino e previsão
    def fit(self, X, Y, epochs):
        """Treino em lote completo. Devolve a perda de treino de cada época."""
        history = []
        for _ in range(epochs):
            Y_hat = self.forward(X)
            history.append(self.loss(Y_hat, Y))
            self.backward(Y)
            self.update()
        return np.array(history)

    def predict(self, X):
        Y_hat = self.forward(X)
        if self.output == "sigmoid":
            return (Y_hat[:, 0] >= 0.5).astype(int)
        return Y_hat.argmax(axis=1)

    def n_params(self):
        return int(sum(W.size + b.size for W, b in zip(self.W, self.b)))

"""Exercise 2 — classificação binária que uma reta não resolve.

Gera figures/fig1.png, fig2.png, fig3.png e imprime os números do relatório.
Roda de qualquer diretório: python ex2_binary.py
"""
from pathlib import Path

import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from sklearn.datasets import make_blobs
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression

from mlp import MLP

rng = np.random.default_rng(42)
FIG = Path(__file__).resolve().parent.parent / "figures"
FIG.mkdir(exist_ok=True)
COLORS = {0: "tab:blue", 1: "tab:orange"}

# ---------------------------------------------------------------- A — dados
# 500 pontos no centro (0, 0) e 250 em cada um dos centros (3, 3) e (-3, -3)
X, cluster = make_blobs(n_samples=[500, 250, 250], centers=[(0, 0), (3, 3), (-3, -3)],
                        cluster_std=1.2, random_state=42)
y = (cluster > 0).astype(int)                 # classe 1 = os dois clusters de fora
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42, stratify=y)
print(f"A) treino {X_train.shape}, teste {X_test.shape}, classe 1 no teste: {y_test.mean():.2f}")

fig, ax = plt.subplots(figsize=(6, 5.5))
for c in (0, 1):
    ax.scatter(*X[y == c].T, s=10, alpha=0.6, color=COLORS[c], label=f"classe {c}")
ax.set_xlabel("$x_1$"); ax.set_ylabel("$x_2$")
ax.set_title("Figure 1 — Exercise 2: os 1000 pontos, coloridos por classe")
ax.legend(); ax.set_aspect("equal")
fig.savefig(FIG / "fig1.png", bbox_inches="tight"); plt.close(fig)

# ---------------------------------------------------------------- B — baseline linear
logreg = LogisticRegression().fit(X_train, y_train)
base_acc = logreg.score(X_test, y_test)
print(f"B) regressão logística: acurácia no teste = {base_acc:.4f}")
print(f"   coeficientes = {logreg.coef_[0].round(4)}, intercepto = {logreg.intercept_[0]:.4f}")
print(f"   fração do teste prevista como classe 1 = {logreg.predict(X_test).mean():.3f}")

# ---------------------------------------------------------------- C — MLP
SIZES, LR, EPOCHS = [2, 8, 1], 0.5, 2000
model = MLP(SIZES, output="sigmoid", lr=LR, rng=rng)
history = model.fit(X_train, y_train[:, None], EPOCHS)
acc_train = (model.predict(X_train) == y_train).mean()
acc_test = (model.predict(X_test) == y_test).mean()
print(f"C) arquitetura {SIZES} (tanh na oculta, sigmoide na saída), lr = {LR}, {EPOCHS} épocas, "
      f"{model.n_params()} parâmetros")
print(f"   perda de treino: inicial {history[0]:.4f}, final {history[-1]:.4f}")
print(f"   acurácia: treino {acc_train:.4f}, teste {acc_test:.4f}")

fig, ax = plt.subplots(figsize=(7, 4))
ax.plot(np.arange(1, EPOCHS + 1), history, label="perda de treino (BCE)")
ax.set_xlabel("época"); ax.set_ylabel("entropia cruzada binária")
ax.set_title("Figure 2 — Exercise 2: perda de treino por época")
ax.legend()
fig.savefig(FIG / "fig2.png", bbox_inches="tight"); plt.close(fig)

# Fronteira de decisão: avalia a rede numa grade e pinta as duas regiões
x1, x2 = np.meshgrid(np.linspace(X[:, 0].min() - 1, X[:, 0].max() + 1, 300),
                     np.linspace(X[:, 1].min() - 1, X[:, 1].max() + 1, 300))
grid_pred = model.predict(np.c_[x1.ravel(), x2.ravel()]).reshape(x1.shape)
fig, ax = plt.subplots(figsize=(6, 5.5))
ax.contourf(x1, x2, grid_pred, levels=[-0.5, 0.5, 1.5], colors=[COLORS[0], COLORS[1]], alpha=0.2)
for c in (0, 1):
    ax.scatter(*X_test[y_test == c].T, s=18, color=COLORS[c], edgecolors="k", linewidths=0.3,
               label=f"classe {c} (teste)")
ax.set_xlabel("$x_1$"); ax.set_ylabel("$x_2$")
ax.set_title(f"Figure 3 — Exercise 2: fronteira do MLP {SIZES} (teste: {acc_test:.1%})")
ax.legend(); ax.set_aspect("equal")
fig.savefig(FIG / "fig3.png", bbox_inches="tight"); plt.close(fig)

# ---------------------------------------------------------------- D — forma da região da classe 1
# Onde a região da classe 1 fica: proporção da grade em cada lado da diagonal x1 + x2 = 0
region1 = grid_pred == 1
s = x1 + x2
print(f"D) região da classe 1: {region1.mean():.1%} da grade; "
      f"acima de x1+x2=0: {region1[s > 0].mean():.1%}, abaixo: {region1[s < 0].mean():.1%}")

"""Dados dos Exercises 3 e 4 (os dois usam exatamente o mesmo conjunto e a mesma divisão)."""
import numpy as np
from sklearn.datasets import make_classification
from sklearn.model_selection import train_test_split


def load_multiclass():
    """1500 amostras, 4 features informativas, 3 classes, 2 clusters por classe; divisão 80/20 estratificada."""
    X, y = make_classification(n_samples=1500, n_features=4, n_informative=4, n_redundant=0,
                               n_repeated=0, n_classes=3, n_clusters_per_class=2,
                               class_sep=1.0, random_state=42)
    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42, stratify=y)
    return X_train, X_test, y_train, y_test


def one_hot(y, n_classes=3):
    return np.eye(n_classes)[y]

"""Exercise 3 — multiclasse com a mesma rede do Exercise 2.

Gera figures/fig4.png e fig5.png e imprime os números do relatório.
Roda de qualquer diretório: python ex3_multiclass.py
"""
from pathlib import Path

import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import confusion_matrix

from mlp import MLP                      # mesma classe do Exercise 2, sem nenhuma alteração
from multiclass_data import load_multiclass, one_hot

rng = np.random.default_rng(42)
FIG = Path(__file__).resolve().parent.parent / "figures"
FIG.mkdir(exist_ok=True)

# ---------------------------------------------------------------- A — dados
X_train, X_test, y_train, y_test = load_multiclass()
print(f"A) treino {X_train.shape}, teste {X_test.shape}, amostras por classe no teste: {np.bincount(y_test)}")

# ---------------------------------------------------------------- C — treino (a rede é a do Exercise 2)
SIZES, LR, EPOCHS = [4, 16, 3], 0.5, 2000
model = MLP(SIZES, output="softmax", lr=LR, rng=rng)
history = model.fit(X_train, one_hot(y_train), EPOCHS)
pred = model.predict(X_test)
acc_train = (model.predict(X_train) == y_train).mean()
acc_test = (pred == y_test).mean()
print(f"C) arquitetura {SIZES} (tanh na oculta, softmax na saída), lr = {LR}, {EPOCHS} épocas, "
      f"{model.n_params()} parâmetros")
print(f"   perda de treino: inicial {history[0]:.4f}, final {history[-1]:.4f}")
print(f"   acurácia: treino {acc_train:.4f}, teste {acc_test:.4f}")

fig, ax = plt.subplots(figsize=(7, 4))
ax.plot(np.arange(1, EPOCHS + 1), history, label="perda de treino (entropia cruzada)")
ax.set_xlabel("época"); ax.set_ylabel("entropia cruzada categórica")
ax.set_title(f"Figure 4 — Exercise 3: perda de treino por época, MLP {SIZES}")
ax.legend()
fig.savefig(FIG / "fig4.png", bbox_inches="tight"); plt.close(fig)

cm = confusion_matrix(y_test, pred)
print("   matriz de confusão (linha = real, coluna = prevista):")
print(cm)
fig, ax = plt.subplots(figsize=(5, 4.5))
im = ax.imshow(cm, cmap="Blues")
for i in range(3):
    for j in range(3):
        ax.text(j, i, cm[i, j], ha="center", va="center",
                color="white" if cm[i, j] > cm.max() / 2 else "black")
ax.set_xticks(range(3)); ax.set_yticks(range(3))
ax.set_xticklabels([f"classe {k}" for k in range(3)])
ax.set_yticklabels([f"classe {k}" for k in range(3)])
ax.set_xlabel("classe prevista"); ax.set_ylabel("classe real")
ax.set_title(f"Figure 5 — Exercise 3: matriz de confusão\n(teste, acurácia {acc_test:.1%})")
fig.colorbar(im, ax=ax, label="nº de amostras")
fig.savefig(FIG / "fig5.png", bbox_inches="tight"); plt.close(fig)

# ---------------------------------------------------------------- D — par mais confundido
pairs = {(i, j): cm[i, j] + cm[j, i] for i in range(3) for j in range(i + 1, 3)}
worst = max(pairs, key=pairs.get)
print(f"D) erros por par de classes (nos dois sentidos): {pairs}; mais confundido: {worst}")


# Estrutura dos clusters: com shuffle=False o make_classification devolve os mesmos pontos e rótulos
# (só a ordem das linhas e das colunas muda, o que não altera distâncias), agrupados em blocos de 250
# por cluster (6 clusters = 2 por classe). Isso permite medir onde cada cluster está.
from sklearn.datasets import make_classification
X_all, y_all = make_classification(n_samples=1500, n_features=4, n_informative=4, n_redundant=0,
                                   n_repeated=0, n_classes=3, n_clusters_per_class=2,
                                   class_sep=1.0, random_state=42, shuffle=False)
cluster_id = np.repeat(np.arange(6), 250)
centers = np.vstack([X_all[cluster_id == k].mean(0) for k in range(6)])
cluster_class = [int(np.bincount(y_all[cluster_id == k]).argmax()) for k in range(6)]
print("   centros dos 6 clusters (classe: centro):")
for k in range(6):
    print(f"   cluster {k} (classe {cluster_class[k]}): {centers[k].round(2)}")
print("   pares de clusters de classes diferentes a menos de 2.5 de distância (vizinhos no hipercubo):")
for i in range(6):
    for j in range(i + 1, 6):
        d = np.linalg.norm(centers[i] - centers[j])
        if cluster_class[i] != cluster_class[j] and d < 2.5:
            print(f"   cluster {i} (classe {cluster_class[i]}) - cluster {j} (classe {cluster_class[j]}): {d:.2f}")
for (i, j) in pairs:
    d = min(np.linalg.norm(centers[a] - centers[b]) for a in range(6) for b in range(6)
            if {cluster_class[a], cluster_class[b]} == {i, j})
    print(f"   classes {i} e {j}: clusters mais próximos a {d:.2f}")
print(f"   rótulos trocados pelo flip_y padrão (1%): {(y_all != np.array(cluster_class)[cluster_id]).sum()} de 1500")
spread = [np.linalg.norm(X_all[cluster_id == k] - centers[k], axis=1).mean() for k in range(6)]
print(f"   distância média de um ponto ao centro do próprio cluster: {np.round(spread, 2).tolist()}")
mid_14 = (centers[1] + centers[4]) / 2
print(f"   cluster 2 (classe 2) até o ponto médio dos clusters 1 e 4 (classe 1): {np.linalg.norm(centers[2] - mid_14):.2f}; "
      f"metade da distância entre 1 e 4: {np.linalg.norm(centers[1] - centers[4]) / 2:.2f}")

logreg = LogisticRegression(max_iter=1000).fit(X_train, y_train)
print(f"D) regressão logística (baseline linear): acurácia no teste = {logreg.score(X_test, y_test):.4f}")

"""Exercise 4 — uma segunda camada oculta ajuda?

Mesmos dados, mesma divisão, mesma taxa de aprendizado e mesmo número de épocas do
Exercise 3. A única mudança é a profundidade. Gera figures/fig6.png e imprime a tabela.
Roda de qualquer diretório: python ex4_depth.py
"""
from pathlib import Path

import numpy as np
import pandas as pd
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

from mlp import MLP
from multiclass_data import load_multiclass, one_hot

FIG = Path(__file__).resolve().parent.parent / "figures"
FIG.mkdir(exist_ok=True)

X_train, X_test, y_train, y_test = load_multiclass()
Y_train = one_hot(y_train)

LR, EPOCHS = 0.5, 2000                       # mesmo orçamento do Exercise 3
ARCHS = {"Exercise 3 (1 camada oculta)": [4, 16, 3],
         "Exercise 4 (2 camadas ocultas)": [4, 16, 16, 3]}
SEEDS = [42, 43, 44]                         # seed 42 é a execução principal; 43 e 44 entram na média


def run(sizes, seed):
    """Treina uma rede com pesos iniciais vindos de default_rng(seed)."""
    rng = np.random.default_rng(seed)
    model = MLP(sizes, output="softmax", lr=LR, rng=rng)
    history = model.fit(X_train, Y_train, EPOCHS)
    return model, history, (model.predict(X_test) == y_test).mean(), (model.predict(X_train) == y_train).mean()


rows, curves = [], {}
for name, sizes in ARCHS.items():
    accs, losses, accs_train = [], [], []
    for seed in SEEDS:
        model, history, acc, acc_train = run(sizes, seed)
        accs.append(acc); losses.append(history[-1]); accs_train.append(acc_train)
        if seed == 42:
            curves[name] = history
            print(f"A) {name}: {sizes}, seed 42 -> perda final {history[-1]:.4f}, acurácia no teste {acc:.4f}")
    print(f"   {name}: acurácias por seed {np.round(accs, 4).tolist()}, perdas finais {np.round(losses, 4).tolist()}")
    rows.append({"modelo": name, "camadas ocultas": len(sizes) - 2, "parâmetros": model.n_params(),
                 "perda final de treino (seed 42)": losses[0], "acurácia no teste (seed 42)": accs[0],
                 "perda final de treino (média 3 seeds)": np.mean(losses),
                 "acurácia no treino (média 3 seeds)": np.mean(accs_train),
                 "acurácia no teste (média 3 seeds)": np.mean(accs),
                 "acurácia no teste (desvio 3 seeds)": np.std(accs)})

table = pd.DataFrame(rows).set_index("modelo")
print("B) tabela de comparação:")
print(table.round(4).to_string())
diff = rows[1]["acurácia no teste (média 3 seeds)"] - rows[0]["acurácia no teste (média 3 seeds)"]
print(f"B) diferença de acurácia (Ex. 4 - Ex. 3), média de 3 seeds: {diff * 100:+.2f} p.p.")

# Estabilidade do treino (seed 42): em quantas épocas a perda subiu em vez de cair
for name, history in curves.items():
    rises = np.diff(history)
    msg = f"C) {name}: a perda subiu em {(rises > 0).sum()} de {EPOCHS} épocas"
    if (rises > 0).any():
        k = rises.argmax()                      # maior subida: da época k+1 para a época k+2
        msg += f"; maior subida na época {k + 2}: {history[k]:.3f} -> pico {history[k:k + 5].max():.3f}"
    print(msg)

fig, ax = plt.subplots(figsize=(8, 4.5))
for name, history in curves.items():
    ax.plot(np.arange(1, EPOCHS + 1), history, label=f"{name}: {ARCHS[name]}")
ax.set_xlabel("época"); ax.set_ylabel("entropia cruzada categórica")
ax.set_title("Figure 6 — Exercises 3 e 4: perda de treino por época (seed 42)")
ax.legend()
fig.savefig(FIG / "fig6.png", bbox_inches="tight"); plt.close(fig)