# Configuración general y librerías

import warnings
warnings.filterwarnings('ignore')

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

from sklearn.datasets import fetch_covtype
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.cluster import KMeans, AgglomerativeClustering, DBSCAN
from sklearn.mixture import GaussianMixture
from sklearn.metrics import (
    silhouette_score,
    calinski_harabasz_score,
    davies_bouldin_score,
    adjusted_rand_score,
    normalized_mutual_info_score,
    v_measure_score
)

sns.set_theme(style='whitegrid', context='notebook')
plt.rcParams['figure.figsize'] = (9, 5)

SEED = 42
np.random.seed(SEED)

# Cargamos CoverType
cov = fetch_covtype(as_frame=True)
df = cov.frame.copy()

# Convertimos etiquetas a [0..6] para comodidad en visualización/métricas externas
df['Cover_Type'] = df['Cover_Type'] - 1

print('Shape total:', df.shape)
print('Número de clases reales:', df['Cover_Type'].nunique())
df.head()

# Revisión estructural y nulos
print(df.info())
print('\nTotal de nulos en dataset:', int(df.isnull().sum().sum()))

# Distribución de clases reales
class_counts = df['Cover_Type'].value_counts().sort_index()
class_props = (class_counts / class_counts.sum()).round(4)

display(pd.DataFrame({'count': class_counts, 'proportion': class_props}))

plt.figure(figsize=(8, 4))
sns.barplot(x=class_counts.index, y=class_counts.values, palette='viridis')
plt.title('Distribución de clases reales (solo referencia)')
plt.xlabel('Clase real')
plt.ylabel('Número de muestras')
plt.tight_layout()
plt.show()

# EDA visual sobre una muestra para no saturar gráficos
eda_sample = df.sample(6000, random_state=SEED)

num_cols = [
    'Elevation', 'Aspect', 'Slope',
    'Horizontal_Distance_To_Hydrology',
    'Vertical_Distance_To_Hydrology',
    'Horizontal_Distance_To_Roadways',
    'Hillshade_9am', 'Hillshade_Noon', 'Hillshade_3pm',
    'Horizontal_Distance_To_Fire_Points'
]

_ = eda_sample[num_cols + ['Cover_Type']].hist(bins=30, figsize=(14, 10), edgecolor='black')
plt.suptitle('Distribuciones (muestra EDA)', y=1.02)
plt.tight_layout()
plt.show()

# Muestreo estratificado para trabajo práctico
N_SAMPLE = 30_000

sample_df = (
    df.groupby('Cover_Type', group_keys=False)
      .apply(lambda x: x.sample(max(1, int(np.floor(len(x) / len(df) * N_SAMPLE))), random_state=SEED))
)

# Ajuste fino para tener tamaño exacto
if len(sample_df) > N_SAMPLE:
    sample_df = sample_df.sample(N_SAMPLE, random_state=SEED)
elif len(sample_df) < N_SAMPLE:
    faltan = N_SAMPLE - len(sample_df)
    extra = df.drop(sample_df.index).sample(faltan, random_state=SEED)
    sample_df = pd.concat([sample_df, extra], axis=0)

sample_df = sample_df.sample(frac=1, random_state=SEED).reset_index(drop=True)

X = sample_df.drop(columns=['Cover_Type']).values
y_true = sample_df['Cover_Type'].values.astype(int)

print('Shape muestra de trabajo:', X.shape)
print('Clases reales en muestra:', np.unique(y_true))

# Estandarización (clave para distancias y DBSCAN)
scaler = StandardScaler()
X_sc = scaler.fit_transform(X)

# PCA para visualización 2D
pca2 = PCA(n_components=2, random_state=SEED)
X_pca2 = pca2.fit_transform(X_sc)

# PCA para reducir ruido/dimensión en algunos experimentos densidad
pca20 = PCA(n_components=20, random_state=SEED)
X_pca20 = pca20.fit_transform(X_sc)

print('Varianza explicada PCA 2D:', pca2.explained_variance_ratio_.sum().round(4))
print('Varianza explicada PCA 20D:', pca20.explained_variance_ratio_.sum().round(4))

# Visualización base en 2D por clase real (solo referencia)
plt.figure(figsize=(8, 6))
scatter = plt.scatter(X_pca2[:, 0], X_pca2[:, 1], c=y_true, s=8, alpha=0.45, cmap='tab10')
plt.title('PCA 2D coloreado por clase real (diagnóstico)')
plt.xlabel('PC1')
plt.ylabel('PC2')
plt.colorbar(scatter, label='Clase real')
plt.tight_layout()
plt.show()

def safe_internal_metrics(X_eval, labels):
    """Métricas internas robustas a casos degenerados (1 cluster o todo ruido)."""
    unique = np.unique(labels)
    # Ignoramos ruido -1 de DBSCAN para métricas internas cuando exista
    mask = labels != -1

    if len(unique) <= 1 or mask.sum() < 5 or len(np.unique(labels[mask])) <= 1:
        return {'Silhouette': np.nan, 'CalinskiHarabasz': np.nan, 'DaviesBouldin': np.nan}

    X_m = X_eval[mask]
    y_m = labels[mask]

    return {
        'Silhouette': silhouette_score(X_m, y_m),
        'CalinskiHarabasz': calinski_harabasz_score(X_m, y_m),
        'DaviesBouldin': davies_bouldin_score(X_m, y_m)
    }


def external_metrics(y_true, labels):
    """Métricas externas usando clases reales solo para diagnóstico."""
    return {
        'ARI': adjusted_rand_score(y_true, labels),
        'NMI': normalized_mutual_info_score(y_true, labels),
        'VMeasure': v_measure_score(y_true, labels)
    }


def summarize_model(name, X_eval, y_true, labels):
    row = {'Modelo': name, 'N_clusters': len(np.unique(labels[labels != -1]))}
    row.update(safe_internal_metrics(X_eval, labels))
    row.update(external_metrics(y_true, labels))
    row['Noise_ratio'] = float(np.mean(labels == -1))
    return row


def plot_cluster_map(X2d, labels, title):
    plt.figure(figsize=(8, 6))
    plt.scatter(X2d[:, 0], X2d[:, 1], c=labels, s=8, alpha=0.5, cmap='tab20')
    plt.title(title)
    plt.xlabel('PC1')
    plt.ylabel('PC2')
    plt.tight_layout()
    plt.show()

k_range = range(2, 13)
km_inertia, km_silhouette = [], []

for k in k_range:
    km_tmp = KMeans(n_clusters=k, n_init=20, random_state=SEED)
    labels_tmp = km_tmp.fit_predict(X_sc)

    km_inertia.append(km_tmp.inertia_)
    km_silhouette.append(silhouette_score(X_sc, labels_tmp))

fig, axes = plt.subplots(1, 2, figsize=(14, 4.5))
axes[0].plot(list(k_range), km_inertia, marker='o')
axes[0].set_title('K-Means: método del codo (inercia)')
axes[0].set_xlabel('K')
axes[0].set_ylabel('Inercia')

axes[1].plot(list(k_range), km_silhouette, marker='o', color='darkgreen')
axes[1].set_title('K-Means: silhouette vs K')
axes[1].set_xlabel('K')
axes[1].set_ylabel('Silhouette')

plt.tight_layout()
plt.show()

k_best = int(list(k_range)[np.argmax(km_silhouette)])
print('K sugerido por silhouette:', k_best)

kmeans = KMeans(n_clusters=k_best, n_init=30, random_state=SEED)
kmeans_labels = kmeans.fit_predict(X_sc)

plot_cluster_map(X_pca2, kmeans_labels, f'K-Means (K={k_best}) en PCA2D')

kmeans_row = summarize_model(f'KMeans (K={k_best})', X_sc, y_true, kmeans_labels)
kmeans_row

# Submuestra para búsqueda de hiperparámetros jerárquicos
idx_sub = np.random.RandomState(SEED).choice(len(X_sc), size=6000, replace=False)
X_h_sub = X_sc[idx_sub]
y_h_sub = y_true[idx_sub]

linkages = ['ward', 'complete', 'average']
k_range_h = range(2, 13)
rows_h = []

for link in linkages:
    for k in k_range_h:
        # Ward requiere distancia euclídea y funciona bien en datos estandarizados
        model = AgglomerativeClustering(n_clusters=k, linkage=link)
        labels_h = model.fit_predict(X_h_sub)

        row = {'linkage': link, 'k': k}
        row.update(safe_internal_metrics(X_h_sub, labels_h))
        row.update(external_metrics(y_h_sub, labels_h))
        rows_h.append(row)

ag_search = pd.DataFrame(rows_h)
ag_search.head()

# Curvas de silhouette por linkage (análogo de curva de aprendizaje para elección de k)
plt.figure(figsize=(9, 5))
for link in linkages:
    sub = ag_search[ag_search['linkage'] == link].sort_values('k')
    plt.plot(sub['k'], sub['Silhouette'], marker='o', label=link)

plt.title('Agglomerative: silhouette vs número de clusters')
plt.xlabel('n_clusters')
plt.ylabel('Silhouette')
plt.legend()
plt.tight_layout()
plt.show()

best_ag = ag_search.sort_values('Silhouette', ascending=False).iloc[0]
best_link, best_k = best_ag['linkage'], int(best_ag['k'])
print('Mejor configuración jerárquica (submuestra):', best_link, '| k =', best_k)

# Ajuste final en todo el conjunto de trabajo
agg = AgglomerativeClustering(n_clusters=best_k, linkage=best_link)
agg_labels = agg.fit_predict(X_sc)

plot_cluster_map(X_pca2, agg_labels, f'Agglomerative ({best_link}, k={best_k}) en PCA2D')

agg_row = summarize_model(f'Agglomerative ({best_link}, k={best_k})', X_sc, y_true, agg_labels)
agg_row

eps_values = np.linspace(1.2, 3.2, 9)
min_samples = 12

db_rows = []
for eps in eps_values:
    db = DBSCAN(eps=float(eps), min_samples=min_samples, n_jobs=-1)
    labels_db = db.fit_predict(X_pca20)

    n_clusters = len(np.unique(labels_db[labels_db != -1]))
    noise_ratio = np.mean(labels_db == -1)
    internal = safe_internal_metrics(X_pca20, labels_db)

    db_rows.append({
        'eps': float(eps),
        'n_clusters': n_clusters,
        'noise_ratio': noise_ratio,
        'Silhouette': internal['Silhouette']
    })

db_scan = pd.DataFrame(db_rows)
db_scan

fig, axes = plt.subplots(1, 2, figsize=(14, 4.5))

axes[0].plot(db_scan['eps'], db_scan['n_clusters'], marker='o')
axes[0].set_title('DBSCAN: clusters detectados vs eps')
axes[0].set_xlabel('eps')
axes[0].set_ylabel('Nº clusters (sin ruido)')

axes[1].plot(db_scan['eps'], db_scan['noise_ratio'], marker='o', color='firebrick')
axes[1].set_title('DBSCAN: ratio de ruido vs eps')
axes[1].set_xlabel('eps')
axes[1].set_ylabel('Proporción de ruido')

plt.tight_layout()
plt.show()

# Elegimos eps con mejor silhouette válida
valid = db_scan.dropna(subset=['Silhouette'])
if len(valid) == 0:
    eps_best = 2.0
else:
    eps_best = float(valid.sort_values('Silhouette', ascending=False).iloc[0]['eps'])

print('eps elegido para DBSCAN:', eps_best)

dbscan = DBSCAN(eps=eps_best, min_samples=min_samples, n_jobs=-1)
db_labels = dbscan.fit_predict(X_pca20)

plot_cluster_map(X_pca2, db_labels, f'DBSCAN (eps={eps_best:.2f}, min_samples={min_samples}) en PCA2D')

db_row = summarize_model(f'DBSCAN (eps={eps_best:.2f})', X_pca20, y_true, db_labels)
db_row

comp_range = range(2, 13)
bic_vals, aic_vals = [], []

for c in comp_range:
    gmm_tmp = GaussianMixture(
        n_components=c,
        covariance_type='full',
        random_state=SEED,
        n_init=3,
        max_iter=200
    )
    gmm_tmp.fit(X_pca20)
    bic_vals.append(gmm_tmp.bic(X_pca20))
    aic_vals.append(gmm_tmp.aic(X_pca20))

fig, axes = plt.subplots(1, 2, figsize=(14, 4.5))
axes[0].plot(list(comp_range), bic_vals, marker='o')
axes[0].set_title('GMM: BIC vs nº componentes')
axes[0].set_xlabel('n_components')
axes[0].set_ylabel('BIC (menor es mejor)')

axes[1].plot(list(comp_range), aic_vals, marker='o', color='purple')
axes[1].set_title('GMM: AIC vs nº componentes')
axes[1].set_xlabel('n_components')
axes[1].set_ylabel('AIC (menor es mejor)')

plt.tight_layout()
plt.show()

gmm_best_c = int(list(comp_range)[np.argmin(bic_vals)])
print('n_components elegido por BIC:', gmm_best_c)

# Curva de convergencia aproximada de EM (lower bound) aumentando max_iter con warm_start
iters = list(range(5, 105, 5))

gmm_prog = GaussianMixture(
    n_components=gmm_best_c,
    covariance_type='full',
    random_state=SEED,
    n_init=1,
    warm_start=True,
    max_iter=1
)

lower_bounds = []
for it in iters:
    gmm_prog.set_params(max_iter=it)
    gmm_prog.fit(X_pca20)
    lower_bounds.append(gmm_prog.lower_bound_)

plt.figure(figsize=(8, 4.5))
plt.plot(iters, lower_bounds, marker='o')
plt.title('GMM: evolución del lower bound (EM)')
plt.xlabel('max_iter')
plt.ylabel('Lower bound (mayor es mejor)')
plt.tight_layout()
plt.show()

gmm = GaussianMixture(
    n_components=gmm_best_c,
    covariance_type='full',
    random_state=SEED,
    n_init=5,
    max_iter=200
)
gmm.fit(X_pca20)
gmm_labels = gmm.predict(X_pca20)

plot_cluster_map(X_pca2, gmm_labels, f'GMM (n_components={gmm_best_c}) en PCA2D')

gmm_row = summarize_model(f'GMM (C={gmm_best_c})', X_pca20, y_true, gmm_labels)
gmm_row

results = pd.DataFrame([kmeans_row, agg_row, db_row, gmm_row])

# Ordenamos por silhouette (si hay NaN, quedan al final)
results_sorted = results.sort_values('Silhouette', ascending=False, na_position='last').reset_index(drop=True)
results_sorted

fig, axes = plt.subplots(2, 2, figsize=(15, 10))

sns.barplot(data=results_sorted, x='Silhouette', y='Modelo', ax=axes[0, 0], palette='Blues')
axes[0, 0].set_title('Silhouette (mayor es mejor)')

sns.barplot(data=results_sorted, x='DaviesBouldin', y='Modelo', ax=axes[0, 1], palette='Reds_r')
axes[0, 1].set_title('Davies-Bouldin (menor es mejor)')

sns.barplot(data=results_sorted, x='ARI', y='Modelo', ax=axes[1, 0], palette='Greens')
axes[1, 0].set_title('ARI vs clases reales (diagnóstico)')

sns.barplot(data=results_sorted, x='Noise_ratio', y='Modelo', ax=axes[1, 1], palette='Purples')
axes[1, 1].set_title('Ratio de ruido detectado')

plt.tight_layout()
plt.show()

label_dict = {
    results_sorted.loc[0, 'Modelo']: None
}

# Mapeo explícito para evitar ambigüedad en el notebook
model_labels_map = {
    kmeans_row['Modelo']: kmeans_labels,
    agg_row['Modelo']: agg_labels,
    db_row['Modelo']: db_labels,
    gmm_row['Modelo']: gmm_labels
}

best_name = results_sorted.loc[0, 'Modelo']
best_labels = model_labels_map[best_name]

counts = pd.Series(best_labels).value_counts().sort_index()
print('Mejor modelo por silhouette:', best_name)
display(counts.to_frame('size'))

plt.figure(figsize=(8, 4))
sns.barplot(x=counts.index.astype(str), y=counts.values, palette='mako')
plt.title(f'Tamaño de clusters del mejor modelo: {best_name}')
plt.xlabel('ID cluster')
plt.ylabel('Número de muestras')
plt.tight_layout()
plt.show()

assert len(results_sorted) == 4, 'Deben compararse 4 modelos'
assert results_sorted['N_clusters'].notnull().all(), 'N_clusters no debe tener nulos'
assert (results_sorted['N_clusters'] >= 1).all(), 'Cada modelo debe tener al menos 1 cluster'
assert np.isfinite(results_sorted['ARI']).all(), 'ARI debe ser finito'
assert np.isfinite(results_sorted['NMI']).all(), 'NMI debe ser finito'
assert np.isfinite(results_sorted['VMeasure']).all(), 'VMeasure debe ser finito'

print('Silhouette máximo:', float(results_sorted['Silhouette'].dropna().max()))
print('Mejor modelo (silhouette):', results_sorted.iloc[0]['Modelo'])
print('✅ Sanity checks completados correctamente')