1# Author: Simon-Pierre Boucher — simon-pierre.boucher@uqo.ca2"""Feature matrix for the statistical / ML baselines (arms A–D).34Every feature is observable as of the valuation date (sale date − 1 day):5* subject attributes and the roll in force at the sale (public before the sale);6* calendar position (``month_index``) — the model learns the market level from *training* sales only;7* ``prior_sale_*`` (strictly earlier sale of the same parcel);8* static geo layers (``geo_features.parquet``);9* **as-of comparable summaries** computed on sales strictly before the sale date:10 - ``nb_*``: **all** same-class sales of other parcels within 3 km dated in [sale date − 365 d, sale date − 1]11 (exact; monthly KD-trees of candidate sales): count, median log price, median price/m2, median assessment12 ratio, mean distance;13 - ``mun_*``: municipality × class rolling medians over the 365 days strictly before the sale date14 (polars ``rolling`` with a left-closed window — same-day sales are excluded), exact and cheap.15The target is ``log_price``; ``price`` and ``date`` are carried for evaluation and never used as inputs.16"""17from __future__ import annotations1819import warnings20from typing import Any2122import numpy as np23import polars as pl24from scipy.spatial import cKDTree2526from ..data.analytic_schema import PROPERTY_CLASSES2728NB_RADIUS_M = 3000.029NB_WINDOW_DAYS = 36530MONTH0 = (2021, 1)3132ID_COLUMNS = ["tx_id", "split", "log_price", "price", "date", "parcel_id"]33CATEGORICAL = ["region_name", "region_group", "mun_code", "property_class", "subtype", "grid1k", "grid5k", "urban_class"]34GEO_NUMERIC = ["in_flood_zone", "transit_stops_within_1km", "da_population_density", "da_median_income", "da_pct_renters", "distance_to_cbd_km"]353637def _month_index(d: pl.Expr) -> pl.Expr:38 return (d.dt.year() - MONTH0[0]) * 12 + (d.dt.month() - MONTH0[1])394041def subject_features(analytic: pl.DataFrame) -> pl.DataFrame:42 """Row-wise features from the analytic table alone (no cross-row information)."""43 roll_ref = pl.col("roll_reference_date").str.to_date("%Y-%m-%d", strict=False)44 return analytic.with_columns(45 pl.col("floor_area_m2").log().alias("log_floor_area"),46 pl.col("land_area_m2").log().alias("log_land_area"),47 (pl.col("date").dt.year() - pl.col("year_built")).cast(pl.Float64).alias("age_years"),48 pl.col("storeys").alias("storeys_f"),49 pl.col("n_units").cast(pl.Float64).alias("n_units_f"),50 pl.when(pl.col("roll_total_value") > 0).then(pl.col("roll_total_value").log()).alias("log_roll_total"),51 pl.when(pl.col("roll_land_value") > 0).then(pl.col("roll_land_value").log()).alias("log_roll_land"),52 pl.when(pl.col("roll_building_value") > 0).then(pl.col("roll_building_value").log()).alias("log_roll_building"),53 pl.when(pl.col("roll_total_value") > 0).then(pl.col("roll_land_value") / pl.col("roll_total_value")).alias("roll_land_share"),54 pl.when((pl.col("roll_total_value") > 0) & (pl.col("floor_area_m2") > 0)).then((pl.col("roll_total_value") / pl.col("floor_area_m2")).log()).alias("log_roll_per_m2"),55 (((pl.col("date") - roll_ref).dt.total_days()) / 30.4375).alias("roll_age_months"),56 _month_index(pl.col("date")).cast(pl.Float64).alias("month_index"),57 pl.col("date").dt.month().cast(pl.Float64).alias("sale_month"),58 pl.col("prior_sale_price").log().alias("prior_sale_log_price"),59 (((pl.col("date") - pl.col("prior_sale_date")).dt.total_days()) / 30.4375).alias("months_since_prior_sale"),60 pl.col("new_construction").cast(pl.Float64).alias("new_construction_f"),61 )626364def municipality_rolling(analytic: pl.DataFrame) -> pl.DataFrame:65 """Municipality × class medians over the 365 days strictly before each sale date66 (tx_id, mun_n_12m, mun_med_log_price_12m, mun_med_ppm2_12m, mun_med_ratio_12m). Row-aligned ``rolling_*_by``67 expressions with a left-closed window: same-day and later sales are never in the window."""68 base = analytic.select(69 "tx_id", "mun_code", "property_class", "date",70 pl.col("log_price"),71 pl.when(pl.col("floor_area_m2") > 0).then(pl.col("price") / pl.col("floor_area_m2")).alias("ppm2"),72 pl.when(pl.col("roll_total_value") > 0).then(pl.col("price") / pl.col("roll_total_value")).alias("ratio"),73 ).sort(["mun_code", "property_class", "date", "tx_id"])74 grp = ["mun_code", "property_class"]75 win = dict(window_size=f"{NB_WINDOW_DAYS}d", closed="left", min_samples=1)76 out = base.with_columns(77 pl.col("log_price").is_not_null().cast(pl.Int32).rolling_sum_by("date", **win).over(grp).fill_null(0).alias("mun_n_12m"),78 pl.col("log_price").rolling_median_by("date", **win).over(grp).alias("mun_med_log_price_12m"),79 pl.col("ppm2").rolling_median_by("date", **win).over(grp).alias("mun_med_ppm2_12m"),80 pl.col("ratio").rolling_median_by("date", **win).over(grp).alias("mun_med_ratio_12m"),81 )82 return out.select("tx_id", "mun_n_12m", "mun_med_log_price_12m", "mun_med_ppm2_12m", "mun_med_ratio_12m")838485def _grouped_median(subj: np.ndarray, val: np.ndarray, n_subj: int) -> np.ndarray:86 """Median of ``val`` per subject id (NaN ignored), vectorised through sorting."""87 out = np.full(n_subj, np.nan)88 ok = ~np.isnan(val)89 subj, val = subj[ok], val[ok]90 if subj.size == 0:91 return out92 order = np.lexsort((val, subj))93 subj, val = subj[order], val[order]94 starts = np.flatnonzero(np.r_[True, subj[1:] != subj[:-1]])95 counts = np.diff(np.r_[starts, subj.size])96 lo = starts + (counts - 1) // 297 hi = starts + counts // 298 out[subj[starts]] = 0.5 * (val[lo] + val[hi])99 return out100101102def neighbour_summary(analytic: pl.DataFrame, radius_m: float = NB_RADIUS_M, window_days: int = NB_WINDOW_DAYS) -> pl.DataFrame:103 """Exact as-of neighbourhood summaries: for every sale, all same-class sales within ``radius_m`` whose date lies104 in [date − window_days, date − 1], other parcels only. Implementation: per class, one KD-tree per calendar105 month of *candidate* sales; subjects of month m query the trees of months m−12..m, then exact date masks and106 vectorised grouped medians. Because only earlier sales can enter a window, later sales never influence a107 feature (verified by ``leakage.check_feature_matrix_asof``)."""108 frames = []109 for cls in PROPERTY_CLASSES:110 sub = analytic.filter(pl.col("property_class") == cls).sort(["date", "tx_id"])111 n = sub.height112 if n == 0:113 continue114 x = sub["x_m"].to_numpy(); y = sub["y_m"].to_numpy()115 d = sub["date"].cast(pl.Int32).to_numpy().astype(np.int64)116 month = (sub["date"].dt.year() * 12 + sub["date"].dt.month() - 1).to_numpy().astype(np.int64)117 lp = sub["log_price"].to_numpy()118 fa = sub["floor_area_m2"].fill_null(float("nan")).to_numpy()119 price = sub["price"].to_numpy()120 ppm2 = np.where(fa > 0, price / np.where(fa > 0, fa, 1.0), np.nan)121 rt = sub["roll_total_value"].fill_null(float("nan")).to_numpy()122 ratio = np.where(rt > 0, price / np.where(rt > 0, rt, 1.0), np.nan)123 parcel = sub["parcel_id"].cast(pl.Categorical).to_physical().to_numpy().astype(np.int64)124 xy = np.column_stack([x, y])125 months = np.unique(month)126 trees: dict[int, tuple[cKDTree, np.ndarray]] = {}127 for m in months:128 rows = np.flatnonzero(month == m)129 trees[int(m)] = (cKDTree(xy[rows]), rows)130 n_12m = np.zeros(n, dtype=np.int32)131 med_lp = np.full(n, np.nan); med_ppm2 = np.full(n, np.nan); med_ratio = np.full(n, np.nan); mean_dist = np.full(n, np.nan)132 for m in months:133 subj_rows = np.flatnonzero(month == m)134 pairs_s: list[np.ndarray] = []135 pairs_c: list[np.ndarray] = []136 for mm in range(int(m) - 13, int(m) + 1):137 t = trees.get(mm)138 if t is None:139 continue140 tree, cand_rows = t141 lists = tree.query_ball_point(xy[subj_rows], r=radius_m, return_sorted=False)142 lens = np.fromiter((len(l) for l in lists), dtype=np.int64, count=len(lists))143 if lens.sum() == 0:144 continue145 cand = cand_rows[np.concatenate([np.asarray(l, dtype=np.int64) for l in lists if len(l)])]146 subj = np.repeat(subj_rows, lens)147 keep = (d[cand] < d[subj]) & (d[cand] >= d[subj] - window_days) & (parcel[cand] != parcel[subj])148 pairs_s.append(subj[keep]); pairs_c.append(cand[keep])149 if not pairs_s:150 continue151 S = np.concatenate(pairs_s); C = np.concatenate(pairs_c)152 cnt = np.bincount(S, minlength=n)153 n_12m += cnt.astype(np.int32)154 dist = np.hypot(x[C] - x[S], y[C] - y[S]) / 1000.0155 sum_dist = np.bincount(S, weights=dist, minlength=n)156 with np.errstate(all="ignore"):157 mean_dist = np.where(cnt > 0, sum_dist / np.maximum(cnt, 1), mean_dist)158 for arr, target in ((lp, med_lp), (ppm2, med_ppm2), (ratio, med_ratio)):159 med = _grouped_median(S, arr[C], n)160 target[subj_rows] = med[subj_rows]161 frames.append(pl.DataFrame({162 "tx_id": sub["tx_id"], "nb_n_12m": n_12m, "nb_med_log_price": med_lp, "nb_med_ppm2": med_ppm2,163 "nb_med_ratio": med_ratio, "nb_mean_dist_km": mean_dist,164 }))165 return pl.concat(frames)166167168def build_feature_matrix(analytic: pl.DataFrame, splits: pl.DataFrame, geo: pl.DataFrame | None = None, comps_summary: bool = True) -> pl.DataFrame:169 """Assemble the baseline design matrix (one row per sale in ``splits``)."""170 df = subject_features(analytic).join(splits.select("tx_id", "split"), on="tx_id", how="inner")171 numeric = [172 "log_floor_area", "log_land_area", "age_years", "storeys_f", "n_units_f", "log_roll_total", "log_roll_land", "log_roll_building",173 "roll_land_share", "log_roll_per_m2", "roll_age_months", "lat", "lng", "x_m", "y_m", "month_index", "sale_month",174 "prior_sale_log_price", "months_since_prior_sale", "new_construction_f",175 ]176 extra: list[str] = []177 if comps_summary:178 df = df.join(municipality_rolling(analytic), on="tx_id", how="left").join(neighbour_summary(analytic), on="tx_id", how="left")179 extra += ["mun_n_12m", "mun_med_log_price_12m", "mun_med_ppm2_12m", "mun_med_ratio_12m",180 "nb_n_12m", "nb_med_log_price", "nb_med_ppm2", "nb_med_ratio", "nb_mean_dist_km"]181 if geo is not None:182 g = geo.select(["tx_id", *[c for c in GEO_NUMERIC if c in geo.columns]])183 if "in_flood_zone" in g.columns:184 g = g.with_columns(pl.col("in_flood_zone").cast(pl.Float64))185 df = df.join(g, on="tx_id", how="left")186 extra += [c for c in GEO_NUMERIC if c in g.columns]187 cols = ID_COLUMNS + numeric + extra + CATEGORICAL188 return df.select(cols).with_columns([pl.col(c).cast(pl.Float64) for c in numeric + extra]).sort(["date", "tx_id"])189190191def feature_columns(matrix: pl.DataFrame) -> tuple[list[str], list[str]]:192 """(numeric feature names, categorical feature names) of a built matrix."""193 num = [c for c in matrix.columns if c not in ID_COLUMNS and c not in CATEGORICAL]194 cat = [c for c in CATEGORICAL if c in matrix.columns]195 return num, cat196