1# Author: Simon-Pierre Boucher — simon-pierre.boucher@uqo.ca2"""Column contract of the analytic layer (``data/processed/analytic.parquet`` and ``data/synthetic/analytic.parquet``).34The analytic table is the *only* input of the feature, comparable, agent and evaluation layers; it is identical5for real and synthetic data. ``ANALYTIC_COLUMNS`` lists the 24 columns produced by the synthetic generator6(``data/synthetic.py``) followed by 9 derived columns added here:78* ``region_group`` — Montréal CMA / Québec CMA / Outaouais / other (deterministic rule on region + distance);9* ``x_m``, ``y_m`` — EPSG:32198 (NAD83 / Quebec Lambert) metres, used for all distance computations;10* ``grid1k``, ``grid5k`` — square cell identifiers ``"<floor(x/size)>_<floor(y/size)>"``;11* ``prior_sale_price``, ``prior_sale_date`` — the most recent *strictly earlier* sale of the same ``parcel_id``12 inside the table (as-of information, null when none);13* ``price_tercile`` — 1/2/3 within the full analytic sample, by ``property_class`` (used only for stratification);14* ``log_price`` — natural log of the price (the prediction target).1516``validate_analytic`` raises on contract violations; ``coerce_analytic`` casts/reorders a frame into the contract.17"""18from __future__ import annotations1920import math21from collections import OrderedDict22from typing import Final2324import numpy as np25import polars as pl2627# ----------------------------------------------------------------------------------------------------------------28# Contract29# ----------------------------------------------------------------------------------------------------------------3031SYNTHETIC_COLUMNS: Final[OrderedDict[str, pl.DataType]] = OrderedDict(32 [33 ("tx_id", pl.String),34 ("parcel_id", pl.String),35 ("date", pl.Date),36 ("price", pl.Float64),37 ("lat", pl.Float64),38 ("lng", pl.Float64),39 ("region_name", pl.String),40 ("mun_code", pl.String),41 ("mun_name", pl.String),42 ("borough_name", pl.String),43 ("urban_class", pl.String),44 ("property_class", pl.String),45 ("subtype", pl.String),46 ("floor_area_m2", pl.Float64),47 ("land_area_m2", pl.Float64),48 ("year_built", pl.Int64),49 ("storeys", pl.Float64),50 ("n_units", pl.Int64),51 ("roll_land_value", pl.Float64),52 ("roll_building_value", pl.Float64),53 ("roll_total_value", pl.Float64),54 ("roll_reference_date", pl.String),55 ("new_construction", pl.Boolean),56 ("exclude_reason", pl.String),57 ]58)5960DERIVED_COLUMNS: Final[OrderedDict[str, pl.DataType]] = OrderedDict(61 [62 ("region_group", pl.String),63 ("x_m", pl.Float64),64 ("y_m", pl.Float64),65 ("grid1k", pl.String),66 ("grid5k", pl.String),67 ("prior_sale_price", pl.Float64),68 ("prior_sale_date", pl.Date),69 ("price_tercile", pl.Int8),70 ("log_price", pl.Float64),71 ]72)7374ANALYTIC_COLUMNS: Final[OrderedDict[str, pl.DataType]] = OrderedDict(list(SYNTHETIC_COLUMNS.items()) + list(DERIVED_COLUMNS.items()))7576#: columns that must never be null in an analytic table77REQUIRED_NON_NULL: Final[tuple[str, ...]] = (78 "tx_id", "parcel_id", "date", "price", "lat", "lng", "region_name", "mun_code", "mun_name", "urban_class",79 "property_class", "region_group", "x_m", "y_m", "grid1k", "grid5k", "price_tercile", "log_price",80)8182PROPERTY_CLASSES: Final[tuple[str, ...]] = ("house", "condo", "plex")83SUBTYPES: Final[tuple[str, ...]] = ("detached", "semi", "row", "integrated", "unknown")84URBAN_CLASSES: Final[tuple[str, ...]] = ("metro", "suburban", "rural")85REGION_GROUPS: Final[tuple[str, ...]] = ("Montréal CMA", "Québec CMA", "Outaouais", "other")8687#: centres used by the region-group rule (lat, lng) and radius in km88MONTREAL_CENTRE: Final[tuple[float, float]] = (45.5017, -73.5673)89QUEBEC_CENTRE: Final[tuple[float, float]] = (46.8139, -71.2080)90MONTREAL_RADIUS_KM: Final[float] = 45.091QUEBEC_RADIUS_KM: Final[float] = 35.092MONTREAL_REGIONS: Final[frozenset[str]] = frozenset({"Montréal", "Laval", "Montérégie", "Laurentides", "Lanaudière"})9394EARTH_RADIUS_KM: Final[float] = 6371.0088959697# ----------------------------------------------------------------------------------------------------------------98# Geometry helpers (shared by features)99# ----------------------------------------------------------------------------------------------------------------100101def haversine_km(lat1: np.ndarray | float, lng1: np.ndarray | float, lat2: np.ndarray | float, lng2: np.ndarray | float) -> np.ndarray | float:102 """Great-circle distance in km (vectorised)."""103 p1, p2 = np.radians(lat1), np.radians(lat2)104 dphi = p2 - p1105 dlmb = np.radians(lng2) - np.radians(lng1)106 a = np.sin(dphi / 2) ** 2 + np.cos(p1) * np.cos(p2) * np.sin(dlmb / 2) ** 2107 return 2 * EARTH_RADIUS_KM * np.arcsin(np.sqrt(a))108109110def haversine_expr(lat_col: str, lng_col: str, lat0: float, lng0: float) -> pl.Expr:111 """Polars expression: haversine distance (km) from (lat_col, lng_col) to a fixed point."""112 p1 = pl.col(lat_col).radians()113 p2 = math.radians(lat0)114 dphi = p2 - p1115 dlmb = math.radians(lng0) - pl.col(lng_col).radians()116 a = (dphi / 2).sin() ** 2 + p1.cos() * math.cos(p2) * (dlmb / 2).sin() ** 2117 return 2 * EARTH_RADIUS_KM * a.sqrt().arcsin()118119120def project_xy(lat: np.ndarray, lng: np.ndarray) -> tuple[np.ndarray, np.ndarray]:121 """WGS84 → EPSG:32198 (NAD83 / Quebec Lambert) metres. Uses pyproj when available; otherwise an122 equirectangular approximation centred on the projection's origin (adequate for distances < 50 km)."""123 lat = np.asarray(lat, dtype=float)124 lng = np.asarray(lng, dtype=float)125 try:126 from pyproj import Transformer127128 tr = Transformer.from_crs("EPSG:4326", "EPSG:32198", always_xy=True)129 x, y = tr.transform(lng, lat)130 return np.asarray(x, dtype=float), np.asarray(y, dtype=float)131 except Exception: # pragma: no cover - pyproj is installed in the reference environment132 lat0, lng0 = 44.0, -68.5 # EPSG:32198 latitude/longitude of origin133 x = np.radians(lng - lng0) * EARTH_RADIUS_KM * 1000 * np.cos(np.radians(lat))134 y = np.radians(lat - lat0) * EARTH_RADIUS_KM * 1000135 return x, y136137138def grid_id_expr(x_col: str, y_col: str, size_m: int) -> pl.Expr:139 """Cell id ``"<floor(x/size)>_<floor(y/size)>"`` — same convention as the WP 22 grids."""140 gx = (pl.col(x_col) / size_m).floor().cast(pl.Int64).cast(pl.String)141 gy = (pl.col(y_col) / size_m).floor().cast(pl.Int64).cast(pl.String)142 return gx + "_" + gy143144145def region_group_expr() -> pl.Expr:146 """Deterministic rule mapping (region_name, lat, lng) to the four analysis groups."""147 d_mtl = haversine_expr("lat", "lng", *MONTREAL_CENTRE)148 d_qc = haversine_expr("lat", "lng", *QUEBEC_CENTRE)149 return (150 pl.when(pl.col("region_name").is_in(list(MONTREAL_REGIONS)) & (d_mtl <= MONTREAL_RADIUS_KM)).then(pl.lit("Montréal CMA"))151 .when((pl.col("region_name") == "Capitale-Nationale") & (d_qc <= QUEBEC_RADIUS_KM)).then(pl.lit("Québec CMA"))152 .when(pl.col("region_name") == "Outaouais").then(pl.lit("Outaouais"))153 .otherwise(pl.lit("other"))154 )155156157# ----------------------------------------------------------------------------------------------------------------158# Derived columns159# ----------------------------------------------------------------------------------------------------------------160161def add_prior_sale(df: pl.DataFrame) -> pl.DataFrame:162 """Most recent sale of the same ``parcel_id`` with ``date < date`` (strictly earlier; same-day transfers are163 not visible). Ties on the earlier date are broken by the highest price (a deterministic, conservative choice)."""164 hist = (165 df.select("parcel_id", pl.col("date").alias("_pdate"), pl.col("price").alias("_pprice"))166 .sort(["parcel_id", "_pdate", "_pprice"])167 )168 left = df.with_row_index("_row").sort(["parcel_id", "date"])169 out = left.join_asof(hist, left_on="date", right_on="_pdate", by="parcel_id", strategy="backward", allow_exact_matches=False, check_sortedness=False)170 out = out.sort("_row").drop("_row")171 return out.rename({"_pprice": "prior_sale_price", "_pdate": "prior_sale_date"})172173174def add_price_tercile(df: pl.DataFrame) -> pl.DataFrame:175 """Tercile (1 = lowest) of ``price`` within ``property_class`` over the whole table (rank-based, ties by tx_id)."""176 rank = pl.struct(["price", "tx_id"]).rank(method="ordinal").over("property_class")177 n = pl.len().over("property_class")178 terc = ((rank - 1) * 3 // n + 1).clip(1, 3).cast(pl.Int8)179 return df.with_columns(terc.alias("price_tercile"))180181182def add_derived_columns(df: pl.DataFrame, recompute_xy: bool = False) -> pl.DataFrame:183 """Add/refresh the 9 derived columns. ``x_m``/``y_m`` are kept when present (WP 22 computed them with the same184 projection) unless ``recompute_xy``."""185 if recompute_xy or "x_m" not in df.columns or "y_m" not in df.columns:186 x, y = project_xy(df["lat"].to_numpy(), df["lng"].to_numpy())187 df = df.with_columns(pl.Series("x_m", x, dtype=pl.Float64), pl.Series("y_m", y, dtype=pl.Float64))188 df = df.with_columns(189 region_group_expr().alias("region_group"),190 grid_id_expr("x_m", "y_m", 1000).alias("grid1k"),191 grid_id_expr("x_m", "y_m", 5000).alias("grid5k"),192 pl.col("price").log().alias("log_price"),193 )194 df = df.drop([c for c in ("prior_sale_price", "prior_sale_date") if c in df.columns])195 df = add_prior_sale(df)196 df = add_price_tercile(df)197 return df198199200# ----------------------------------------------------------------------------------------------------------------201# Coercion and validation202# ----------------------------------------------------------------------------------------------------------------203204def coerce_analytic(df: pl.DataFrame) -> pl.DataFrame:205 """Cast every contract column to its dtype, add missing optional columns as null, drop extras, reorder."""206 exprs: list[pl.Expr] = []207 for name, dtype in ANALYTIC_COLUMNS.items():208 if name in df.columns:209 col = pl.col(name)210 if dtype == pl.Date and df.schema[name] == pl.String:211 col = col.str.to_date("%Y-%m-%d")212 col = col.cast(dtype, strict=True)213 if dtype == pl.Float64:214 col = col.fill_nan(None)215 exprs.append(col.alias(name))216 else:217 exprs.append(pl.lit(None, dtype=dtype).alias(name))218 out = df.select(exprs)219 return out.sort(["date", "tx_id"])220221222def validate_analytic(df: pl.DataFrame, *, allow_extra: bool = False) -> None:223 """Raise ``ValueError`` if ``df`` violates the analytic contract (columns, dtypes, nulls, domains, keys)."""224 problems: list[str] = []225 missing = [c for c in ANALYTIC_COLUMNS if c not in df.columns]226 if missing:227 problems.append(f"missing columns: {missing}")228 extra = [c for c in df.columns if c not in ANALYTIC_COLUMNS]229 if extra and not allow_extra:230 problems.append(f"unexpected columns: {extra}")231 for name, dtype in ANALYTIC_COLUMNS.items():232 if name in df.columns and df.schema[name] != dtype:233 problems.append(f"{name}: dtype {df.schema[name]} != {dtype}")234 if problems:235 raise ValueError("analytic contract violated: " + "; ".join(problems))236 for name in REQUIRED_NON_NULL:237 n = df[name].null_count()238 if n:239 problems.append(f"{name}: {n} nulls")240 if df["tx_id"].n_unique() != df.height:241 problems.append("tx_id not unique")242 if df.height:243 bad_cls = set(df["property_class"].unique().to_list()) - set(PROPERTY_CLASSES)244 if bad_cls:245 problems.append(f"property_class outside {PROPERTY_CLASSES}: {bad_cls}")246 bad_sub = set(df["subtype"].drop_nulls().unique().to_list()) - set(SUBTYPES)247 if bad_sub:248 problems.append(f"subtype outside {SUBTYPES}: {bad_sub}")249 bad_urb = set(df["urban_class"].unique().to_list()) - set(URBAN_CLASSES)250 if bad_urb:251 problems.append(f"urban_class outside {URBAN_CLASSES}: {bad_urb}")252 bad_rg = set(df["region_group"].unique().to_list()) - set(REGION_GROUPS)253 if bad_rg:254 problems.append(f"region_group outside {REGION_GROUPS}: {bad_rg}")255 if (df["price"] <= 0).any():256 problems.append("non-positive prices")257 if ((df["log_price"] - df["price"].log()).abs() > 1e-9).any():258 problems.append("log_price != log(price)")259 if df["price_tercile"].min() < 1 or df["price_tercile"].max() > 3:260 problems.append("price_tercile outside 1..3")261 ps = df.filter(pl.col("prior_sale_date").is_not_null())262 if ps.height and (ps["prior_sale_date"] >= ps["date"]).any():263 problems.append("prior_sale_date not strictly earlier than date")264 if (df["lat"].abs() > 90).any() or (df["lng"].abs() > 180).any():265 problems.append("coordinates out of range")266 if problems:267 raise ValueError("analytic contract violated: " + "; ".join(problems))268