Optimización de hiperparámetros con Hyperopt

pip install torch torchvision hyperopt matplotlib numpy

import torch
import torch.nn as nn
import torch.optim as optim
from torch.utils.data import DataLoader, random_split
from torchvision import datasets, transforms
import numpy as np
import matplotlib.pyplot as plt
import time
import warnings
warnings.filterwarnings('ignore')

from hyperopt import fmin, tpe, hp, rand, Trials, STATUS_OK, space_eval

SEED = 42
torch.manual_seed(SEED)
np.random.seed(SEED)

device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
print(f"Device: {device}")

transform = transforms.Compose([
    transforms.ToTensor(),
    transforms.Normalize((0.2860,), (0.3530,)),  # media y std de FashionMNIST
])

full_train = datasets.FashionMNIST('./data', train=True,  download=True, transform=transform)
test_set   = datasets.FashionMNIST('./data', train=False, download=True, transform=transform)

# Split train/val (54K / 6K)
train_set, val_set = random_split(
    full_train, [54000, 6000],
    generator=torch.Generator().manual_seed(SEED)
)

print(f"Train: {len(train_set)} | Val: {len(val_set)} | Test: {len(test_set)}")
print(f"Clases: {full_train.classes}")
print(f"Input shape: {full_train[0][0].shape}")  # → [1, 28, 28]

class ConfigurableCNN(nn.Module):
    """
    CNN con arquitectura configurable para HPO.
    
    Parámetros configurables:
      n_conv_blocks: número de bloques conv (2-4)
      n_filters:     filtros base por bloque (16, 32, 64)
      activation:    función de activación ('relu', 'leaky_relu', 'gelu')
      fc_dim:        neuronas en la capa FC oculta (64, 128, 256, 512)
      dropout:       dropout rate en el clasificador (0.0 - 0.5)
      use_batchnorm: usar BatchNorm (True/False)
    """
    def __init__(self, n_conv_blocks=3, n_filters=32, activation='relu',
                 fc_dim=128, dropout=0.3, use_batchnorm=True):
        super().__init__()
        
        act_fn = {
            'relu': nn.ReLU, 'leaky_relu': nn.LeakyReLU, 'gelu': nn.GELU
        }[activation]
        
        # ── Bloques convolucionales ──
        layers = []
        in_ch = 1  # FashionMNIST: 1 canal
        for i in range(n_conv_blocks):
            out_ch = n_filters * (2 ** i)  # Doblar filtros cada bloque
            layers.append(nn.Conv2d(in_ch, out_ch, 3, padding=1))
            if use_batchnorm:
                layers.append(nn.BatchNorm2d(out_ch))
            layers.append(act_fn())
            layers.append(nn.Conv2d(out_ch, out_ch, 3, padding=1))
            if use_batchnorm:
                layers.append(nn.BatchNorm2d(out_ch))
            layers.append(act_fn())
            layers.append(nn.MaxPool2d(2))
            in_ch = out_ch
        
        self.features = nn.Sequential(*layers)
        
        # Calcular tamaño del feature map tras las convoluciones
        # 28 → 14 → 7 → 3 → 1 (con cada MaxPool2d(2))
        feat_size = 28 // (2 ** n_conv_blocks)
        flat_dim = in_ch * feat_size * feat_size
        
        # ── Clasificador MLP ──
        self.classifier = nn.Sequential(
            nn.Flatten(),
            nn.Linear(flat_dim, fc_dim),
            nn.ReLU(),
            nn.Dropout(dropout),
            nn.Linear(fc_dim, 10),
        )
        
        self._init_weights()
    
    def _init_weights(self):
        for m in self.modules():
            if isinstance(m, nn.Conv2d):
                nn.init.kaiming_normal_(m.weight, nonlinearity='relu')
            elif isinstance(m, nn.Linear):
                nn.init.kaiming_normal_(m.weight, nonlinearity='relu')
                nn.init.zeros_(m.bias)
    
    def forward(self, x):
        return self.classifier(self.features(x))

# Quick test con valores por defecto
model = ConfigurableCNN().to(device)
n_params = sum(p.numel() for p in model.parameters())
print(f"ConfigurableCNN (default): {n_params:,} params")

# Verificar forward pass
x = torch.randn(4, 1, 28, 28).to(device)
out = model(x)
print(f"Output shape: {out.shape}")  # → [4, 10]

def quick_train(model, train_loader, val_loader, lr, epochs=5):
    """Entrena rápido y devuelve val accuracy."""
    optimizer = optim.Adam(model.parameters(), lr=lr)
    criterion = nn.CrossEntropyLoss()
    
    model.train()
    for epoch in range(epochs):
        for images, targets in train_loader:
            images, targets = images.to(device), targets.to(device)
            optimizer.zero_grad()
            loss = criterion(model(images), targets)
            loss.backward()
            optimizer.step()
    
    # Evaluar
    model.eval()
    correct, total = 0, 0
    with torch.no_grad():
        for images, targets in val_loader:
            images, targets = images.to(device), targets.to(device)
            correct += (model(images).argmax(1) == targets).sum().item()
            total += targets.size(0)
    return correct / total

# DataLoaders con batch_size fijo para el estudio
study_train = DataLoader(train_set, batch_size=128, shuffle=True, num_workers=2)
study_val   = DataLoader(val_set, batch_size=256, shuffle=False, num_workers=2)

# ── Test 1: Sensibilidad al Learning Rate ──────────────────
print("=== Sensibilidad al Learning Rate ===")
for lr in [1e-1, 1e-2, 1e-3, 3e-4, 1e-4]:
    model = ConfigurableCNN(n_conv_blocks=2, n_filters=32).to(device)
    acc = quick_train(model, study_train, study_val, lr=lr, epochs=5)
    print(f"  lr={lr:.0e}  →  val_acc={acc:.2%}")

# ── Test 2: Sensibilidad a la arquitectura ─────────────────
print("\n=== Sensibilidad a la Arquitectura ===")
for blocks, filters in [(2, 16), (2, 32), (3, 32), (3, 64)]:
    model = ConfigurableCNN(n_conv_blocks=blocks, n_filters=filters).to(device)
    n = sum(p.numel() for p in model.parameters())
    acc = quick_train(model, study_train, study_val, lr=1e-3, epochs=5)
    print(f"  blocks={blocks}, filters={filters} ({n:>7,} params)  →  val_acc={acc:.2%}")

# ── Test 3: Sensibilidad al Dropout ────────────────────────
print("\n=== Sensibilidad al Dropout ===")
for drop in [0.0, 0.1, 0.3, 0.5]:
    model = ConfigurableCNN(n_conv_blocks=2, n_filters=32, dropout=drop).to(device)
    acc = quick_train(model, study_train, study_val, lr=1e-3, epochs=5)
    print(f"  dropout={drop}  →  val_acc={acc:.2%}")

def build_model(params):
    """Construye un ConfigurableCNN a partir de un diccionario de HPs."""
    model = ConfigurableCNN(
        n_conv_blocks=params['n_conv_blocks'],
        n_filters=params['n_filters'],
        activation=params.get('activation', 'relu'),
        fc_dim=params['fc_dim'],
        dropout=params['dropout'],
        use_batchnorm=True,  # Fijo: siempre BatchNorm
    )
    return model.to(device)

# Test
test_params = {
    'n_conv_blocks': 3, 'n_filters': 32, 'activation': 'relu',
    'fc_dim': 128, 'dropout': 0.3,
}
m = build_model(test_params)
print(f"Modelo con {sum(p.numel() for p in m.parameters()):,} params")

# hp.choice(label, options)
#   → Elige un elemento de una lista finita. Para categorías y discretos.
#   Ejemplo: hp.choice('optimizer', ['adam', 'sgd'])

# hp.uniform(label, low, high)
#   → Valor continuo uniforme en [low, high]. Para rangos lineales.
#   Ejemplo: hp.uniform('dropout', 0.0, 0.5)

# hp.loguniform(label, low, high)
#   → exp(uniform(low, high)). Para valores que varían en órdenes de magnitud.
#   Ejemplo: hp.loguniform('lr', np.log(1e-4), np.log(1e-1))
#   ⚠️ low y high son el LOG del rango deseado

# hp.quniform(label, low, high, q)
#   → round(uniform(low, high) / q) * q. Para discretos con paso.
#   Ejemplo: hp.quniform('neurons', 32, 512, 16)

# hp.normal(label, mu, sigma)
#   → Distribución normal. Para centrar la búsqueda alrededor de un valor conocido.

# hp.lognormal(label, mu, sigma)
#   → exp(normal(mu, sigma)). Para LR cuando tienes un prior fuerte.

search_space = {
    # ── Optimización ──
    'lr': hp.loguniform('lr', np.log(1e-4), np.log(1e-1)),
    'weight_decay': hp.loguniform('weight_decay', np.log(1e-6), np.log(1e-2)),
    'optimizer': hp.choice('optimizer', ['adam', 'sgd']),
    'batch_size': hp.choice('batch_size', [32, 64, 128, 256]),

    # ── Arquitectura ──
    'n_conv_blocks': hp.choice('n_conv_blocks', [2, 3]),
    'n_filters': hp.choice('n_filters', [16, 32, 64]),
    'fc_dim': hp.choice('fc_dim', [64, 128, 256, 512]),
    'activation': hp.choice('activation', ['relu', 'leaky_relu', 'gelu']),

    # ── Regularización ──
    'dropout': hp.uniform('dropout', 0.0, 0.5),
}

# Verificar el espacio
from hyperopt.pyll.stochastic import sample
example = sample(search_space)
print("Ejemplo de punto muestreado:")
for k, v in sorted(example.items()):
    print(f"  {k:>16s}: {v}")

NUM_EPOCHS_HPO = 10  # Épocas por trial (suficiente para discriminar)

def objective(params):
    """
    Función objetivo para Hyperopt.
    Entrena un modelo con los HPs dados y devuelve -val_accuracy.
    (Hyperopt MINIMIZA, así que negamos la accuracy.)
    """
    try:
        # 1. Construir modelo
        model = build_model(params)

        # 2. DataLoaders con batch_size del HP
        bs = params['batch_size']
        tr_loader = DataLoader(train_set, batch_size=bs, shuffle=True,
                               num_workers=2, pin_memory=True)
        vl_loader = DataLoader(val_set, batch_size=256, shuffle=False,
                               num_workers=2, pin_memory=True)

        # 3. Optimizer
        if params['optimizer'] == 'adam':
            opt = optim.Adam(model.parameters(),
                             lr=params['lr'],
                             weight_decay=params['weight_decay'])
        else:  # sgd
            opt = optim.SGD(model.parameters(),
                            lr=params['lr'],
                            momentum=0.9,
                            weight_decay=params['weight_decay'])

        criterion = nn.CrossEntropyLoss()

        # 4. Entrenar
        model.train()
        for epoch in range(NUM_EPOCHS_HPO):
            for images, targets in tr_loader:
                images, targets = images.to(device), targets.to(device)
                opt.zero_grad()
                loss = criterion(model(images), targets)
                loss.backward()
                nn.utils.clip_grad_norm_(model.parameters(), max_norm=1.0)
                opt.step()

        # 5. Evaluar en validación
        model.eval()
        correct, total = 0, 0
        with torch.no_grad():
            for images, targets in vl_loader:
                images, targets = images.to(device), targets.to(device)
                preds = model(images).argmax(dim=1)
                correct += (preds == targets).sum().item()
                total += targets.size(0)

        val_acc = correct / total

        # 6. Retornar resultado
        return {
            'loss': -val_acc,       # Hyperopt minimiza → negamos la accuracy
            'status': STATUS_OK,
            'val_acc': val_acc,     # Metadata extra
            'n_params': sum(p.numel() for p in model.parameters()),
        }

    except Exception as e:
        # Si algo falla (OOM, dimensiones incompatibles...), retornar loss mala
        print(f"  ❌ Trial fallido: {e}")
        return {'loss': 0.0, 'status': STATUS_OK, 'val_acc': 0.0}

# Quick test con los HPs por defecto
result = objective(test_params | {'lr': 1e-3, 'weight_decay': 1e-4,
                                   'optimizer': 'adam', 'batch_size': 128})
print(f"Test objective: val_acc={result['val_acc']:.2%}, loss={result['loss']:.4f}")

N_RANDOM_TRIALS = 30  # Número de configuraciones a probar

# Trials almacena todo el historial de la búsqueda
random_trials = Trials()

print(f"🎲 Iniciando Random Search ({N_RANDOM_TRIALS} trials)...")
print("=" * 60)

best_random = fmin(
    fn=objective,
    space=search_space,
    algo=rand.suggest,       # ← Random search
    max_evals=N_RANDOM_TRIALS,
    trials=random_trials,
    rstate=np.random.default_rng(SEED),
    verbose=False,
)

# Convertir índices a valores reales
best_random_params = space_eval(search_space, best_random)

print("\n" + "=" * 60)
print("🏆 Mejor configuración (Random Search):")
for k, v in sorted(best_random_params.items()):
    print(f"  {k:>16s}: {v}")

best_random_acc = -random_trials.best_trial['result']['loss']
print(f"\n  → Val accuracy: {best_random_acc:.2%}")

# Extraer accuracies de todos los trials
random_accs = [-t['result']['loss'] for t in random_trials.trials]

# Estadísticas
print(f"Trials completados: {len(random_accs)}")
print(f"Mejor accuracy:     {max(random_accs):.2%}")
print(f"Peor accuracy:      {min(random_accs):.2%}")
print(f"Media:              {np.mean(random_accs):.2%}")
print(f"Std:                {np.std(random_accs):.2%}")

# Distribución de accuracies
fig, axes = plt.subplots(1, 2, figsize=(12, 4))

# Histograma
axes[0].hist(random_accs, bins=15, color='#6c5ce7', alpha=0.7, edgecolor='white')
axes[0].axvline(max(random_accs), color='#00b894', linestyle='--', label=f'Best: {max(random_accs):.2%}')
axes[0].set_xlabel('Val Accuracy'); axes[0].set_ylabel('Nº Trials')
axes[0].set_title('Distribución de Accuracies (Random Search)')
axes[0].legend()

# Accuracy vs trial number (¿mejora con más trials?)
axes[1].scatter(range(1, len(random_accs)+1), random_accs,
                c='#6c5ce7', alpha=0.6, s=30)
best_so_far = [max(random_accs[:i+1]) for i in range(len(random_accs))]
axes[1].plot(range(1, len(random_accs)+1), best_so_far,
             color='#00b894', linewidth=2, label='Best so far')
axes[1].set_xlabel('Trial'); axes[1].set_ylabel('Val Accuracy')
axes[1].set_title('Convergencia (Random Search)')
axes[1].legend()

plt.tight_layout()
plt.savefig('random_search_results.png', dpi=150)
plt.show()

# Extraer HPs de cada trial
lrs = [t['misc']['vals']['lr'][0] for t in random_trials.trials]

fig, axes = plt.subplots(1, 3, figsize=(14, 4))

# LR vs Accuracy (el más revelador)
axes[0].scatter(lrs, random_accs, c=random_accs, cmap='viridis', s=40, alpha=0.8)
axes[0].set_xscale('log')
axes[0].set_xlabel('Learning Rate (log)'); axes[0].set_ylabel('Val Accuracy')
axes[0].set_title('LR vs Accuracy')

# Dropout vs Accuracy
dropouts = [t['misc']['vals']['dropout'][0] for t in random_trials.trials]
axes[1].scatter(dropouts, random_accs, c=random_accs, cmap='viridis', s=40, alpha=0.8)
axes[1].set_xlabel('Dropout'); axes[1].set_ylabel('Val Accuracy')
axes[1].set_title('Dropout vs Accuracy')

# Weight decay vs Accuracy
wds = [t['misc']['vals']['weight_decay'][0] for t in random_trials.trials]
axes[2].scatter(wds, random_accs, c=random_accs, cmap='viridis', s=40, alpha=0.8)
axes[2].set_xscale('log')
axes[2].set_xlabel('Weight Decay (log)'); axes[2].set_ylabel('Val Accuracy')
axes[2].set_title('Weight Decay vs Accuracy')

plt.tight_layout()
plt.savefig('hp_analysis_random.png', dpi=150)
plt.show()

N_TPE_TRIALS = 50  # Presupuesto total (incluye warm-start)

# Crear nuevo objeto Trials para TPE
tpe_trials = Trials()

print(f"🧠 Iniciando TPE ({N_TPE_TRIALS} trials)...")
print("=" * 60)

best_tpe = fmin(
    fn=objective,
    space=search_space,
    algo=tpe.suggest,          # ← TPE bayesiano
    max_evals=N_TPE_TRIALS,
    trials=tpe_trials,
    rstate=np.random.default_rng(SEED),
    verbose=False,
)

best_tpe_params = space_eval(search_space, best_tpe)

print("\n" + "=" * 60)
print("🏆 Mejor configuración (TPE):")
for k, v in sorted(best_tpe_params.items()):
    print(f"  {k:>16s}: {v}")

best_tpe_acc = -tpe_trials.best_trial['result']['loss']
print(f"\n  → Val accuracy: {best_tpe_acc:.2%}")

random_accs = [-t['result']['loss'] for t in random_trials.trials]
tpe_accs    = [-t['result']['loss'] for t in tpe_trials.trials]

# Best-so-far curves
random_best = [max(random_accs[:i+1]) for i in range(len(random_accs))]
tpe_best    = [max(tpe_accs[:i+1]) for i in range(len(tpe_accs))]

fig, axes = plt.subplots(1, 2, figsize=(13, 5))

# ── Convergencia ──
axes[0].plot(range(1, len(random_best)+1), random_best,
             color='#fdcb6e', linewidth=2, label=f'Random ({max(random_accs):.2%})')
axes[0].plot(range(1, len(tpe_best)+1), tpe_best,
             color='#6c5ce7', linewidth=2, label=f'TPE ({max(tpe_accs):.2%})')
axes[0].set_xlabel('Trial')
axes[0].set_ylabel('Best Val Accuracy')
axes[0].set_title('Convergencia: Random Search vs TPE')
axes[0].legend()
axes[0].grid(alpha=0.2)

# ── Distribuciones ──
axes[1].hist(random_accs, bins=12, color='#fdcb6e', alpha=0.5, label='Random', edgecolor='white')
axes[1].hist(tpe_accs, bins=12, color='#6c5ce7', alpha=0.5, label='TPE', edgecolor='white')
axes[1].set_xlabel('Val Accuracy')
axes[1].set_ylabel('Nº Trials')
axes[1].set_title('Distribución de accuracies')
axes[1].legend()

plt.tight_layout()
plt.savefig('random_vs_tpe.png', dpi=150)
plt.show()

# ── Resumen ──
print("\n📊 Resumen comparativo:")
print(f"{'':>20s} {'Random':>10s}  {'TPE':>10s}")
print(f"{'Trials':>20s} {len(random_accs):>10d}  {len(tpe_accs):>10d}")
print(f"{'Mejor accuracy':>20s} {max(random_accs):>10.2%}  {max(tpe_accs):>10.2%}")
print(f"{'Media accuracy':>20s} {np.mean(random_accs):>10.2%}  {np.mean(tpe_accs):>10.2%}")
print(f"{'Std':>20s} {np.std(random_accs):>10.2%}  {np.std(tpe_accs):>10.2%}")

# Reentrenar el mejor modelo con más épocas
print("🔄 Reentrenando mejor configuración con 30 épocas...")
best_model = build_model(best_tpe_params)

opt = optim.Adam(best_model.parameters(),
                 lr=best_tpe_params['lr'],
                 weight_decay=best_tpe_params['weight_decay'])
criterion = nn.CrossEntropyLoss()
bs = best_tpe_params['batch_size']
final_train = DataLoader(train_set, batch_size=bs, shuffle=True, num_workers=2)
test_loader = DataLoader(test_set, batch_size=256, shuffle=False, num_workers=2)

best_model.train()
for epoch in range(30):
    for images, targets in final_train:
        images, targets = images.to(device), targets.to(device)
        opt.zero_grad()
        loss = criterion(best_model(images), targets)
        loss.backward()
        nn.utils.clip_grad_norm_(best_model.parameters(), max_norm=1.0)
        opt.step()

# Evaluar en TEST (la primera y única vez)
best_model.eval()
correct, total = 0, 0
with torch.no_grad():
    for images, targets in test_loader:
        images, targets = images.to(device), targets.to(device)
        correct += (best_model(images).argmax(1) == targets).sum().item()
        total += targets.size(0)

test_acc = correct / total
print(f"\n📊 Test accuracy (mejor modelo, 30 épocas): {test_acc:.2%}")
print(f"   (Val accuracy durante HPO fue: {best_tpe_acc:.2%})")

import pandas as pd

# Construir tabla de todos los trials TPE
rows = []
for t in tpe_trials.trials:
    vals = t['misc']['vals']
    row = {
        'val_acc': -t['result']['loss'],
        'lr': vals['lr'][0],
        'weight_decay': vals['weight_decay'][0],
        'dropout': vals['dropout'][0],
        'optimizer': ['adam', 'sgd'][vals['optimizer'][0]],
        'batch_size': [32, 64, 128, 256][vals['batch_size'][0]],
        'n_conv_blocks': [2, 3][vals['n_conv_blocks'][0]],
        'n_filters': [16, 32, 64][vals['n_filters'][0]],
        'fc_dim': [64, 128, 256, 512][vals['fc_dim'][0]],
        'activation': ['relu', 'leaky_relu', 'gelu'][vals['activation'][0]],
    }
    rows.append(row)

df = pd.DataFrame(rows)

# Correlación con val_acc (para HPs numéricos)
numeric_cols = ['lr', 'weight_decay', 'dropout', 'n_conv_blocks', 'n_filters', 'fc_dim']
correlations = df[numeric_cols + ['val_acc']].corr()['val_acc'].drop('val_acc').abs()
correlations = correlations.sort_values(ascending=False)

print("📊 Correlación absoluta con val_acc:")
for hp_name, corr in correlations.items():
    bar = '█' * int(corr * 40)
    print(f"  {hp_name:>16s}: {corr:.3f} {bar}")

# Top 5 mejores trials
print("\n🏆 Top 5 configuraciones:")
top5 = df.nlargest(5, 'val_acc')
print(top5[['val_acc', 'lr', 'n_filters', 'n_conv_blocks', 'fc_dim',
            'dropout', 'optimizer']].to_string(index=False))

fig, ax = plt.subplots(figsize=(8, 5))

scatter = ax.scatter(
    df['lr'], df['val_acc'],
    c=range(len(df)), cmap='plasma',
    s=60, alpha=0.7, edgecolors='white', linewidth=0.5,
)
ax.set_xscale('log')
ax.set_xlabel('Learning Rate (log scale)', fontsize=12)
ax.set_ylabel('Val Accuracy', fontsize=12)
ax.set_title('TPE Exploration: LR vs Accuracy (color = trial order)', fontsize=13)
cbar = plt.colorbar(scatter, ax=ax, label='Trial #')

# Marcar el mejor
best_idx = df['val_acc'].idxmax()
ax.scatter(df.loc[best_idx, 'lr'], df.loc[best_idx, 'val_acc'],
           s=200, facecolors='none', edgecolors='#00b894', linewidth=3,
           label=f'Best: {df.loc[best_idx, "val_acc"]:.2%}')
ax.legend(fontsize=10)
plt.tight_layout()
plt.savefig('tpe_exploration.png', dpi=150)
plt.show()