# Author: Simon-Pierre Boucher — simon-pierre.boucher@uqo.ca """Column contract of the analytic layer (``data/processed/analytic.parquet`` and ``data/synthetic/analytic.parquet``). The analytic table is the *only* input of the feature, comparable, agent and evaluation layers; it is identical for real and synthetic data. ``ANALYTIC_COLUMNS`` lists the 24 columns produced by the synthetic generator (``data/synthetic.py``) followed by 9 derived columns added here: * ``region_group`` — Montréal CMA / Québec CMA / Outaouais / other (deterministic rule on region + distance); * ``x_m``, ``y_m`` — EPSG:32198 (NAD83 / Quebec Lambert) metres, used for all distance computations; * ``grid1k``, ``grid5k`` — square cell identifiers ``"_"``; * ``prior_sale_price``, ``prior_sale_date`` — the most recent *strictly earlier* sale of the same ``parcel_id`` inside the table (as-of information, null when none); * ``price_tercile`` — 1/2/3 within the full analytic sample, by ``property_class`` (used only for stratification); * ``log_price`` — natural log of the price (the prediction target). ``validate_analytic`` raises on contract violations; ``coerce_analytic`` casts/reorders a frame into the contract. """ from __future__ import annotations import math from collections import OrderedDict from typing import Final import numpy as np import polars as pl # ---------------------------------------------------------------------------------------------------------------- # Contract # ---------------------------------------------------------------------------------------------------------------- SYNTHETIC_COLUMNS: Final[OrderedDict[str, pl.DataType]] = OrderedDict( [ ("tx_id", pl.String), ("parcel_id", pl.String), ("date", pl.Date), ("price", pl.Float64), ("lat", pl.Float64), ("lng", pl.Float64), ("region_name", pl.String), ("mun_code", pl.String), ("mun_name", pl.String), ("borough_name", pl.String), ("urban_class", pl.String), ("property_class", pl.String), ("subtype", pl.String), ("floor_area_m2", pl.Float64), ("land_area_m2", pl.Float64), ("year_built", pl.Int64), ("storeys", pl.Float64), ("n_units", pl.Int64), ("roll_land_value", pl.Float64), ("roll_building_value", pl.Float64), ("roll_total_value", pl.Float64), ("roll_reference_date", pl.String), ("new_construction", pl.Boolean), ("exclude_reason", pl.String), ] ) DERIVED_COLUMNS: Final[OrderedDict[str, pl.DataType]] = OrderedDict( [ ("region_group", pl.String), ("x_m", pl.Float64), ("y_m", pl.Float64), ("grid1k", pl.String), ("grid5k", pl.String), ("prior_sale_price", pl.Float64), ("prior_sale_date", pl.Date), ("price_tercile", pl.Int8), ("log_price", pl.Float64), ] ) ANALYTIC_COLUMNS: Final[OrderedDict[str, pl.DataType]] = OrderedDict(list(SYNTHETIC_COLUMNS.items()) + list(DERIVED_COLUMNS.items())) #: columns that must never be null in an analytic table REQUIRED_NON_NULL: Final[tuple[str, ...]] = ( "tx_id", "parcel_id", "date", "price", "lat", "lng", "region_name", "mun_code", "mun_name", "urban_class", "property_class", "region_group", "x_m", "y_m", "grid1k", "grid5k", "price_tercile", "log_price", ) PROPERTY_CLASSES: Final[tuple[str, ...]] = ("house", "condo", "plex") SUBTYPES: Final[tuple[str, ...]] = ("detached", "semi", "row", "integrated", "unknown") URBAN_CLASSES: Final[tuple[str, ...]] = ("metro", "suburban", "rural") REGION_GROUPS: Final[tuple[str, ...]] = ("Montréal CMA", "Québec CMA", "Outaouais", "other") #: centres used by the region-group rule (lat, lng) and radius in km MONTREAL_CENTRE: Final[tuple[float, float]] = (45.5017, -73.5673) QUEBEC_CENTRE: Final[tuple[float, float]] = (46.8139, -71.2080) MONTREAL_RADIUS_KM: Final[float] = 45.0 QUEBEC_RADIUS_KM: Final[float] = 35.0 MONTREAL_REGIONS: Final[frozenset[str]] = frozenset({"Montréal", "Laval", "Montérégie", "Laurentides", "Lanaudière"}) EARTH_RADIUS_KM: Final[float] = 6371.0088 # ---------------------------------------------------------------------------------------------------------------- # Geometry helpers (shared by features) # ---------------------------------------------------------------------------------------------------------------- def haversine_km(lat1: np.ndarray | float, lng1: np.ndarray | float, lat2: np.ndarray | float, lng2: np.ndarray | float) -> np.ndarray | float: """Great-circle distance in km (vectorised).""" p1, p2 = np.radians(lat1), np.radians(lat2) dphi = p2 - p1 dlmb = np.radians(lng2) - np.radians(lng1) a = np.sin(dphi / 2) ** 2 + np.cos(p1) * np.cos(p2) * np.sin(dlmb / 2) ** 2 return 2 * EARTH_RADIUS_KM * np.arcsin(np.sqrt(a)) def haversine_expr(lat_col: str, lng_col: str, lat0: float, lng0: float) -> pl.Expr: """Polars expression: haversine distance (km) from (lat_col, lng_col) to a fixed point.""" p1 = pl.col(lat_col).radians() p2 = math.radians(lat0) dphi = p2 - p1 dlmb = math.radians(lng0) - pl.col(lng_col).radians() a = (dphi / 2).sin() ** 2 + p1.cos() * math.cos(p2) * (dlmb / 2).sin() ** 2 return 2 * EARTH_RADIUS_KM * a.sqrt().arcsin() def project_xy(lat: np.ndarray, lng: np.ndarray) -> tuple[np.ndarray, np.ndarray]: """WGS84 → EPSG:32198 (NAD83 / Quebec Lambert) metres. Uses pyproj when available; otherwise an equirectangular approximation centred on the projection's origin (adequate for distances < 50 km).""" lat = np.asarray(lat, dtype=float) lng = np.asarray(lng, dtype=float) try: from pyproj import Transformer tr = Transformer.from_crs("EPSG:4326", "EPSG:32198", always_xy=True) x, y = tr.transform(lng, lat) return np.asarray(x, dtype=float), np.asarray(y, dtype=float) except Exception: # pragma: no cover - pyproj is installed in the reference environment lat0, lng0 = 44.0, -68.5 # EPSG:32198 latitude/longitude of origin x = np.radians(lng - lng0) * EARTH_RADIUS_KM * 1000 * np.cos(np.radians(lat)) y = np.radians(lat - lat0) * EARTH_RADIUS_KM * 1000 return x, y def grid_id_expr(x_col: str, y_col: str, size_m: int) -> pl.Expr: """Cell id ``"_"`` — same convention as the WP 22 grids.""" gx = (pl.col(x_col) / size_m).floor().cast(pl.Int64).cast(pl.String) gy = (pl.col(y_col) / size_m).floor().cast(pl.Int64).cast(pl.String) return gx + "_" + gy def region_group_expr() -> pl.Expr: """Deterministic rule mapping (region_name, lat, lng) to the four analysis groups.""" d_mtl = haversine_expr("lat", "lng", *MONTREAL_CENTRE) d_qc = haversine_expr("lat", "lng", *QUEBEC_CENTRE) return ( pl.when(pl.col("region_name").is_in(list(MONTREAL_REGIONS)) & (d_mtl <= MONTREAL_RADIUS_KM)).then(pl.lit("Montréal CMA")) .when((pl.col("region_name") == "Capitale-Nationale") & (d_qc <= QUEBEC_RADIUS_KM)).then(pl.lit("Québec CMA")) .when(pl.col("region_name") == "Outaouais").then(pl.lit("Outaouais")) .otherwise(pl.lit("other")) ) # ---------------------------------------------------------------------------------------------------------------- # Derived columns # ---------------------------------------------------------------------------------------------------------------- def add_prior_sale(df: pl.DataFrame) -> pl.DataFrame: """Most recent sale of the same ``parcel_id`` with ``date < date`` (strictly earlier; same-day transfers are not visible). Ties on the earlier date are broken by the highest price (a deterministic, conservative choice).""" hist = ( df.select("parcel_id", pl.col("date").alias("_pdate"), pl.col("price").alias("_pprice")) .sort(["parcel_id", "_pdate", "_pprice"]) ) left = df.with_row_index("_row").sort(["parcel_id", "date"]) out = left.join_asof(hist, left_on="date", right_on="_pdate", by="parcel_id", strategy="backward", allow_exact_matches=False, check_sortedness=False) out = out.sort("_row").drop("_row") return out.rename({"_pprice": "prior_sale_price", "_pdate": "prior_sale_date"}) def add_price_tercile(df: pl.DataFrame) -> pl.DataFrame: """Tercile (1 = lowest) of ``price`` within ``property_class`` over the whole table (rank-based, ties by tx_id).""" rank = pl.struct(["price", "tx_id"]).rank(method="ordinal").over("property_class") n = pl.len().over("property_class") terc = ((rank - 1) * 3 // n + 1).clip(1, 3).cast(pl.Int8) return df.with_columns(terc.alias("price_tercile")) def add_derived_columns(df: pl.DataFrame, recompute_xy: bool = False) -> pl.DataFrame: """Add/refresh the 9 derived columns. ``x_m``/``y_m`` are kept when present (WP 22 computed them with the same projection) unless ``recompute_xy``.""" if recompute_xy or "x_m" not in df.columns or "y_m" not in df.columns: x, y = project_xy(df["lat"].to_numpy(), df["lng"].to_numpy()) df = df.with_columns(pl.Series("x_m", x, dtype=pl.Float64), pl.Series("y_m", y, dtype=pl.Float64)) df = df.with_columns( region_group_expr().alias("region_group"), grid_id_expr("x_m", "y_m", 1000).alias("grid1k"), grid_id_expr("x_m", "y_m", 5000).alias("grid5k"), pl.col("price").log().alias("log_price"), ) df = df.drop([c for c in ("prior_sale_price", "prior_sale_date") if c in df.columns]) df = add_prior_sale(df) df = add_price_tercile(df) return df # ---------------------------------------------------------------------------------------------------------------- # Coercion and validation # ---------------------------------------------------------------------------------------------------------------- def coerce_analytic(df: pl.DataFrame) -> pl.DataFrame: """Cast every contract column to its dtype, add missing optional columns as null, drop extras, reorder.""" exprs: list[pl.Expr] = [] for name, dtype in ANALYTIC_COLUMNS.items(): if name in df.columns: col = pl.col(name) if dtype == pl.Date and df.schema[name] == pl.String: col = col.str.to_date("%Y-%m-%d") col = col.cast(dtype, strict=True) if dtype == pl.Float64: col = col.fill_nan(None) exprs.append(col.alias(name)) else: exprs.append(pl.lit(None, dtype=dtype).alias(name)) out = df.select(exprs) return out.sort(["date", "tx_id"]) def validate_analytic(df: pl.DataFrame, *, allow_extra: bool = False) -> None: """Raise ``ValueError`` if ``df`` violates the analytic contract (columns, dtypes, nulls, domains, keys).""" problems: list[str] = [] missing = [c for c in ANALYTIC_COLUMNS if c not in df.columns] if missing: problems.append(f"missing columns: {missing}") extra = [c for c in df.columns if c not in ANALYTIC_COLUMNS] if extra and not allow_extra: problems.append(f"unexpected columns: {extra}") for name, dtype in ANALYTIC_COLUMNS.items(): if name in df.columns and df.schema[name] != dtype: problems.append(f"{name}: dtype {df.schema[name]} != {dtype}") if problems: raise ValueError("analytic contract violated: " + "; ".join(problems)) for name in REQUIRED_NON_NULL: n = df[name].null_count() if n: problems.append(f"{name}: {n} nulls") if df["tx_id"].n_unique() != df.height: problems.append("tx_id not unique") if df.height: bad_cls = set(df["property_class"].unique().to_list()) - set(PROPERTY_CLASSES) if bad_cls: problems.append(f"property_class outside {PROPERTY_CLASSES}: {bad_cls}") bad_sub = set(df["subtype"].drop_nulls().unique().to_list()) - set(SUBTYPES) if bad_sub: problems.append(f"subtype outside {SUBTYPES}: {bad_sub}") bad_urb = set(df["urban_class"].unique().to_list()) - set(URBAN_CLASSES) if bad_urb: problems.append(f"urban_class outside {URBAN_CLASSES}: {bad_urb}") bad_rg = set(df["region_group"].unique().to_list()) - set(REGION_GROUPS) if bad_rg: problems.append(f"region_group outside {REGION_GROUPS}: {bad_rg}") if (df["price"] <= 0).any(): problems.append("non-positive prices") if ((df["log_price"] - df["price"].log()).abs() > 1e-9).any(): problems.append("log_price != log(price)") if df["price_tercile"].min() < 1 or df["price_tercile"].max() > 3: problems.append("price_tercile outside 1..3") ps = df.filter(pl.col("prior_sale_date").is_not_null()) if ps.height and (ps["prior_sale_date"] >= ps["date"]).any(): problems.append("prior_sale_date not strictly earlier than date") if (df["lat"].abs() > 90).any() or (df["lng"].abs() > 180).any(): problems.append("coordinates out of range") if problems: raise ValueError("analytic contract violated: " + "; ".join(problems))