LeaveClusterOut¶

import colormaps
import numpy as np
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from sklearn.cluster import HDBSCAN

from geovalidate.cv import LeaveClusterOut

DATA_PATH = '/Users/ljwolf/Dropbox/work/projects/spquadtree/code/kc_house_data.csv'

df = pd.read_csv(DATA_PATH)
gdf = gpd.GeoDataFrame(
    df,
    geometry=gpd.points_from_xy(df['long'], df['lat']),
    crs='EPSG:4326',
)
print(str(len(gdf)) + ' house sales')

# haversine expects (lat, lon) in radians
coords_rad = np.radians(np.column_stack([gdf.geometry.y, gdf.geometry.x]))

best_mcs = 100
for mcs in [100, 150, 200, 250, 300, 400, 500]:
    labels_ = HDBSCAN(min_cluster_size=mcs, metric='haversine').fit(coords_rad).labels_
    n = (np.unique(labels_) >= 0).sum()
    n_noise = (labels_ == -1).sum()
    print('  min_cluster_size=' + str(mcs) + ' -> ' + str(n) + ' clusters, '
          + str(n_noise) + ' noise (' + str(round(100 * n_noise / len(labels_), 1)) + '%)')
    if n >= 8:
        best_mcs = mcs

print('\nUsing min_cluster_size=' + str(best_mcs))
hdb = HDBSCAN(min_cluster_size=best_mcs, metric='haversine').fit(coords_rad)

# Pass the pre-fitted HDBSCAN; LeaveClusterOut detects labels_ and skips refitting
lco = LeaveClusterOut(hdb, noise='train_only')
splits_to = list(lco.split(gdf))

labels = lco.labels_
n_noise = (labels == -1).sum()
print('Clusters (folds) : ' + str(lco.n_clusters_))
print('Noise points     : ' + str(n_noise) + ' (' + str(round(100 * n_noise / len(gdf), 1)) + '%)')

fig, ax = plt.subplots(figsize=(7, 8))
noise_mask = labels == -1

ax.scatter(
    gdf.geometry.x[noise_mask], gdf.geometry.y[noise_mask],
    color='lightgrey', s=0.5, zorder=1,
)
cluster_labels = labels[~noise_mask]
sc = ax.scatter(
    gdf.geometry.x[~noise_mask], gdf.geometry.y[~noise_mask],
    c=cluster_labels, cmap=colormaps.safe[:lco.n_clusters_], s=1.5, zorder=2,
    vmin=0, vmax=lco.n_clusters_ - 1,
)
handles = [
    mpatches.Patch(color=sc.cmap(sc.norm(i)), label='cluster ' + str(i))
    for i in range(lco.n_clusters_)
]
handles.append(mpatches.Patch(color='lightgrey', label='noise (' + str(n_noise) + ')'))
ax.legend(handles=handles, loc='lower left', fontsize=8)
ax.set_title(
    'HDBSCAN -- ' + str(lco.n_clusters_) + ' clusters, '
    + str(n_noise) + ' noise points (min_cluster_size=' + str(best_mcs) + ')'
)
ax.set_xlabel('longitude')
ax.set_ylabel('latitude')
plt.tight_layout()
plt.show()

lco_drop    = LeaveClusterOut(lco.clusterer_, noise='drop')
lco_nearest = LeaveClusterOut(lco.clusterer_, noise='nearest')

splits_drop    = list(lco_drop.split(gdf))
splits_nearest = list(lco_nearest.split(gdf))

stats = []
for i, ((tr_to, te), (tr_dr, _), (tr_nr, _)) in enumerate(
    zip(splits_to, splits_drop, splits_nearest)
):
    stats.append({
        'fold': i,
        'test size': len(te),
        'train (train_only)': len(tr_to),
        'train (drop)': len(tr_dr),
        'train (nearest)': len(tr_nr),
        'noise excluded (drop)': len(tr_to) - len(tr_dr),
    })
pd.DataFrame(stats)

fold = 0
tr_to, te     = splits_to[fold]
tr_dr, _      = splits_drop[fold]
tr_nr, te_nr  = splits_nearest[fold]

noise_in_train   = np.setdiff1d(tr_to, tr_dr)          # noise visible in train_only
noise_in_test_nr = np.intersect1d(te_nr, np.where(lco.labels_ == -1)[0])  # noise -> test in nearest

fig, axes = plt.subplots(1, 3, figsize=(18, 6), sharex=True, sharey=True)

configs = [
    (tr_to, te,    "noise='train_only'"),
    (tr_dr, te,    "noise='drop'"),
    (tr_nr, te_nr, "noise='nearest'"),
]

for ax, (train, test, title) in zip(axes, configs):
    ax.scatter(
        gdf.geometry.x.iloc[train], gdf.geometry.y.iloc[train],
        color='steelblue', s=0.4, alpha=0.4,
    )
    ax.scatter(
        gdf.geometry.x.iloc[test], gdf.geometry.y.iloc[test],
        color='firebrick', s=3, zorder=3,
    )
    ax.set_title(title, fontsize=12)
    ax.set_axis_off()

# highlight noise points on the train_only panel
axes[0].scatter(
    gdf.geometry.x.iloc[noise_in_train], gdf.geometry.y.iloc[noise_in_train],
    color='orange', s=1.5, zorder=2, alpha=0.8,
)

# highlight reassigned noise in the nearest panel (split between test and train)
noise_in_train_nr = np.intersect1d(tr_nr, np.where(lco.labels_ == -1)[0])
axes[2].scatter(
    gdf.geometry.x.iloc[noise_in_test_nr], gdf.geometry.y.iloc[noise_in_test_nr],
    color='darkorange', s=4, zorder=4, alpha=0.9,
)
axes[2].scatter(
    gdf.geometry.x.iloc[noise_in_train_nr], gdf.geometry.y.iloc[noise_in_train_nr],
    color='orange', s=1.5, zorder=2, alpha=0.8,
)

train_patch      = mpatches.Patch(color='steelblue',  label='train (cluster pts)')
test_patch       = mpatches.Patch(color='firebrick',  label='test cluster')
noise_train_patch = mpatches.Patch(color='orange',    label='noise -> train')
noise_test_patch  = mpatches.Patch(color='darkorange', label='noise -> test (nearest)')

axes[0].legend(handles=[train_patch, test_patch, noise_train_patch], loc='lower left', fontsize=8)
axes[1].legend(handles=[train_patch, test_patch], loc='lower left', fontsize=8)
axes[2].legend(handles=[train_patch, test_patch, noise_train_patch, noise_test_patch],
               loc='lower left', fontsize=8)

fig.suptitle('LeaveClusterOut fold ' + str(fold) + ': noise handling comparison', fontsize=13)
plt.tight_layout()
plt.show()