Lesson 0

The first model

The original notebook with its real outputs: CV 0.758, leaderboard 0.70.

Left: what & whyRight: the code and its real output

Lesson 0: the first model (0.70 on Kaggle)

This is the original notebook, exactly as it ran, with its real outputs. It's a strong pipeline: about 350 features, out-of-fold target encodings, two LightGBM models, a blend, and macro-F1 class weights.

Its own CV said 0.758; Kaggle said 0.70. Read it as a good model with one hidden flaw: the validation (step 1).

Step 1: load the data and make folds

StratifiedKFold makes 5 folds with the same tier mix, but rows are assigned at random. One operator's listings land in both training and validation folds.

⚠️ This is the flaw. Kaggle's test set contains operators that never appear in train, so random folds measure an easier task than the real one.

Step 2: look at the target

The tiers are roughly balanced, with tier 2 the largest. Balanced classes make macro F1 a fair metric: each tier counts equally.

train["price_tier"].value_counts().sort_index().plot.bar(title="price_tier distribution");

Step 3: basic features

yes/no flags → 1/0; dates → number of days
operator_ref, zone, borough → frequency encoding (how common the value is)
facility codes → one-hot columns (about 150 of them)
coordinates → rotations, distance to metro/zone centres, k-means map clusters
remaining text columns → pandas category (LightGBM handles them natively)

All of this is computed on train+test together, which is fine because no labels are used.

def basic_features(train, test):
    raw = pd.concat([train.drop(columns="price_tier"), test], ignore_index=True)
    df = raw.copy()

    for c in ["operator_certified", "operator_has_avatar", "operator_id_verified", "bookable"]:
        df[c] = df[c].map({"t": 1, "f": 0})
    for c in ["snapshot_start_date", "snapshot_end_date", "first_review_date", "last_review_date"]:
        df[c] = (pd.to_datetime(df[c]) - pd.Timestamp("2000-01-01")).dt.days
    for c in ["operator_ref", "zone", "borough"]:
        df[c] = raw[c].map(raw[c].value_counts())
    df["n_facilities"] = df["facilities"].str.count("FCL-")

    codes = raw["facilities"].fillna("[]").map(json.loads).explode().dropna()
    fac = pd.crosstab(codes.index, codes).reindex(raw.index, fill_value=0).clip(upper=1)
    fac.columns = ["f_" + c[4:] for c in fac.columns]

    lat, lon = raw["geo_lat"], raw["geo_lon"]
    geo = pd.DataFrame({"lat_p_lon": lat + lon, "lat_m_lon": lat - lon})
    for a in (15, 30, 60):
        r = np.deg2rad(a)
        geo[f"rot{a}_x"] = lat * np.cos(r) - lon * np.sin(r)
        geo[f"rot{a}_y"] = lat * np.sin(r) + lon * np.cos(r)
    for k in ("metro", "zone"):
        center = raw.groupby(k)[["geo_lat", "geo_lon"]].transform("median")
        geo[f"dist_{k}_center"] = np.hypot(lat - center["geo_lat"], lon - center["geo_lon"])
    km = MiniBatchKMeans(100, random_state=0, n_init=3, batch_size=4096).fit(raw[["geo_lat", "geo_lon"]])
    geo["geo_cluster"] = km.labels_
    geo["geo_cluster_dist"] = km.transform(raw[["geo_lat", "geo_lon"]]).min(1)

    df = df.drop(columns=["facilities", "space_id"])
    for c in df.select_dtypes(exclude="number").columns:
        df[c] = df[c].astype("category")
    X = pd.concat([df, fac, geo], axis=1)
    return X.iloc[: len(train)].reset_index(drop=True), X.iloc[len(train):].reset_index(drop=True)


X_train, X_test = basic_features(train, test)
X_train.shape

(279746, 215)

Step 4: out-of-fold target features

For groups such as zone, zone + seats, operator, operator + zone: what's the tier mix of that group in the other folds? Plus residuals vs. an (access type, seats) baseline, and tiers of the nearest map neighbours.

Out-of-fold is the right idea: no row sees its own label. ⚠️ But groups keyed on operator_ref (op, op_zone, op_sc, op_seats, dup, dup2) mostly identify the operator. With random folds they look brilliant. On the real test set the operator is new, or is a different operator with the same name, so they mislead.

TE_GROUPS = {
    "op": ["operator_ref"], "zone": ["zone"], "boro": ["borough"], "sc": ["space_class"],
    "metro_sc": ["metro", "space_class"], "zone_acc": ["zone", "access_type"],
    "zone_seats": ["zone", "max_seats"], "op_zone": ["operator_ref", "zone"],
    "zone_sc": ["zone", "space_class"], "op_sc": ["operator_ref", "space_class"],
    "op_seats": ["operator_ref", "max_seats"], "geo": ["geo_lat", "geo_lon"],
    "dup": ["operator_ref", "zone", "space_class", "max_seats"],
    "dup2": ["operator_ref", "space_class", "max_seats", "access_type"],
    "metro_acc_seats": ["metro", "access_type", "max_seats"],
}
RES_GROUPS = ["op", "zone", "geo", "op_zone", "zone_sc", "sc", "op_sc"]
GEO_KS = (1, 5, 20, 50)

both = pd.concat([train, test], ignore_index=True)
group_ids = {name: both[cols].astype(str).groupby(cols, sort=False, dropna=False).ngroup().to_numpy()
             for name, cols in TE_GROUPS.items()}
base_ids = both[["access_type", "max_seats"]].astype(str).groupby(["access_type", "max_seats"], dropna=False).ngroup().to_numpy()
coords = both[["geo_lat", "geo_lon"]].to_numpy()
prior = np.bincount(y, minlength=K) / len(y)
tiers = np.arange(1, K + 1)


def target_encoding(fit, apply):
    cols = {}
    for name, g in group_ids.items():
        n = g.max() + 1
        counts = np.zeros((n, K))
        np.add.at(counts, (g[fit], y[fit]), 1)
        total = counts.sum(1)
        p = ((counts + 20 * prior) / (total[:, None] + 20))[g[apply]]
        for c in range(K):
            cols[f"te_{name}_p{c + 1}"] = p[:, c]
        cols[f"te_{name}_mean"] = p @ tiers
        cols[f"te_{name}_n"] = np.log1p(total[g[apply]])
    return cols


def residual_encoding(fit, apply):
    yt = y[fit] + 1.0
    nb = base_ids.max() + 1
    base_mean = (np.bincount(base_ids[fit], yt, nb) + 20 * yt.mean()) / (np.bincount(base_ids[fit], minlength=nb) + 20)
    resid = yt - base_mean[base_ids[fit]]
    cols = {}
    for name in RES_GROUPS:
        g = group_ids[name]
        n = g.max() + 1
        cnt = np.bincount(g[fit], minlength=n)[g[apply]] + 10
        mean = np.bincount(g[fit], resid, n)[g[apply]] / cnt
        cols[f"res_{name}_mean"] = mean
        cols[f"res_{name}_rms"] = np.sqrt(np.bincount(g[fit], resid ** 2, n)[g[apply]] / cnt)
        cols[f"res_{name}_pred"] = base_mean[base_ids[apply]] + mean
    return cols


def geo_neighbours(fit, apply):
    points, inv = np.unique(coords[fit], axis=0, return_inverse=True)
    counts = np.zeros((len(points), K))
    np.add.at(counts, (inv.ravel(), y[fit]), 1)
    dist, idx = NearestNeighbors(n_neighbors=max(GEO_KS)).fit(points).kneighbors(coords[apply])
    cols = {}
    for k in GEO_KS:
        c = counts[idx[:, :k]].sum(1)
        p = (c + 5 * prior) / (c.sum(1)[:, None] + 5)
        for j in range(K):
            cols[f"geo{k}_p{j + 1}"] = p[:, j]
        cols[f"geo{k}_mean"] = p @ tiers
        cols[f"geo{k}_dist"] = dist[:, k - 1]
    return cols


def target_features(fit, apply):
    return pd.DataFrame({**target_encoding(fit, apply), **residual_encoding(fit, apply), **geo_neighbours(fit, apply)})


T_train = pd.concat([target_features(fit, val).set_index(val) for fit, val in folds]).sort_index()
T_test = target_features(np.arange(len(train)), np.arange(len(train), len(both)))

X_train = pd.concat([X_train, T_train], axis=1)
X_test = pd.concat([X_test, T_test], axis=1)
X_train.shape

(279746, 350)

Step 5: two LightGBM models

multiclass: 4 probabilities directly.
ordinal: 3 yes/no models, is tier > 1, > 2, > 3?, turned into 4 probabilities. This uses the order of the tiers.

Early stopping picks the number of trees on each validation fold. The fold scores (~0.755) are high because of the random folds.

PARAMS = dict(learning_rate=0.05, num_leaves=127, min_child_samples=20, feature_fraction=0.8,
              bagging_fraction=0.8, bagging_freq=1, verbose=-1, seed=SEED)


def fit_predict(params, X_fit, y_fit, X_val, y_val, X_test):
    model = lgb.train(params, lgb.Dataset(X_fit, y_fit), 5000,
                      valid_sets=[lgb.Dataset(X_val, y_val)],
                      callbacks=[lgb.early_stopping(100, verbose=False)])
    return model.predict(X_val, num_iteration=model.best_iteration), model.predict(X_test, num_iteration=model.best_iteration)


def cumulative_to_classes(p_above):
    p_above = np.minimum.accumulate(p_above, axis=1)
    edges = np.column_stack([np.ones(len(p_above)), p_above, np.zeros(len(p_above))])
    return np.clip(edges[:, :-1] - edges[:, 1:], 1e-6, 1)


oof = {m: np.zeros((len(train), K)) for m in ["multiclass", "ordinal"]}
pred = {m: np.zeros((len(test), K)) for m in ["multiclass", "ordinal"]}

for i, (fit, val) in enumerate(folds):
    Xf, Xv = X_train.iloc[fit], X_train.iloc[val]

    p_val, p_test = fit_predict({**PARAMS, "objective": "multiclass", "num_class": K},
                                Xf, y[fit], Xv, y[val], X_test)
    oof["multiclass"][val], pred["multiclass"] = p_val, pred["multiclass"] + p_test / FOLDS

    above_val, above_test = [], []
    for k in range(K - 1):
        v, t = fit_predict({**PARAMS, "objective": "binary"}, Xf, y[fit] > k, Xv, y[val] > k, X_test)
        above_val.append(v)
        above_test.append(t)
    oof["ordinal"][val] = cumulative_to_classes(np.column_stack(above_val))
    pred["ordinal"] += cumulative_to_classes(np.column_stack(above_test)) / FOLDS

    print(f"fold {i}", {m: round(f1_score(y[val], oof[m][val].argmax(1), average="macro"), 4) for m in oof})

{m: round(f1_score(y, p.argmax(1), average="macro"), 4) for m, p in oof.items()}

fold 0 {'multiclass': 0.7523, 'ordinal': 0.7557}

fold 1 {'multiclass': 0.753, 'ordinal': 0.7551}

fold 2 {'multiclass': 0.7515, 'ordinal': 0.7517}

fold 3 {'multiclass': 0.7509, 'ordinal': 0.7541}

fold 4 {'multiclass': 0.7522, 'ordinal': 0.757}

{'multiclass': 0.752, 'ordinal': 0.7547}

Step 6: blend and class weights

Average the two models, then multiply each tier's probability by a tuned weight to maximise macro F1. CV: 0.7577.

In the final notebook these weights are tuned with cross-fitting (tuned on 4 folds, scored on the 5th), so the reported score isn't flattered by the tuning.

def macro_f1(p, w):
    return f1_score(y, (p * w).argmax(1), average="macro")


def tune_class_weights(p, rounds=3):
    w = np.ones(K)
    for _ in range(rounds):
        for k in range(1, K):
            grid = w[k] * np.exp(np.linspace(-0.3, 0.3, 31))
            scores = [macro_f1(p, np.r_[w[:k], g, w[k + 1:]]) for g in grid]
            w[k] = grid[int(np.argmax(scores))]
    return w


blend_oof = (oof["multiclass"] + oof["ordinal"]) / 2
blend_test = (pred["multiclass"] + pred["ordinal"]) / 2
weights = tune_class_weights(blend_oof)

print("blend CV macro F1:", round(macro_f1(blend_oof, np.ones(K)), 4))
print("with class weights:", round(macro_f1(blend_oof, weights), 4), weights.round(3))

blend CV macro F1: 0.7559
with class weights: 0.7577 [1.    0.98  1.15  0.852]

Step 7: submission

Folds' test predictions are averaged, weighted and turned into tiers. Leaderboard: 0.70.

→ Next: the story of how we found the problem, and Lesson 1, the final notebook.

submission = pd.DataFrame({"ID": test["space_id"], "price_tier": (blend_test * weights).argmax(1) + 1})
submission.to_csv("submission.csv", index=False)
submission["price_tier"].value_counts().sort_index()

price_tier
1    25359
2    37553
3    27169
4    24457
Name: count, dtype: int64