# -*- coding: utf-8 -*- """04 — Moteur de rétro-évaluation. Pour chaque vente cible (2024 = validation, 2025-2026 = test), le prix est masqué ; seules les ventes conclues STRICTEMENT AVANT la date cible (et l'indice/les modèles connus à cette date) sont utilisées. On évalue d'un coup toutes les configurations (sélection × rayon × fenêtre × nombre × ajustement × pondération × réconciliation) sur exactement les mêmes cibles. Sorties : targets.parquet, configs_E1.csv/est_E1.npy (surface k,d,T), configs_E2.csv/est_E2.npy/sd_E2.npy (factoriel sélection × ajustement × pondération × réconciliation), n_avail.npy.""" import os for _v in ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS", "VECLIB_MAXIMUM_THREADS"): os.environ[_v] = "1" # un seul fil BLAS par processus (fork + BLAS multi-fils = interblocage sur macOS) import multiprocessing as mp, time from scipy.spatial import cKDTree from common import * df = load_parquet("sales.parquet") meta = load_json("index_meta.json"); regions = meta["regions"]; M_EVAL0 = meta["M_EVAL0"] rid_map = {r: i for i, r in enumerate(regions)} IDX = np.load(OUT / "index_asof.npy") ASOF = load_json("asof.json")["asof"] MODELS = {a: load_json(f"models_{a}.json") for a in ASOF} FEATS = MODELS[ASOF[0]]["feats"] GBM = {a: np.load(OUT / f"gbm_pred_{a}.npy") for a in ASOF} # ----------------------------------------------------------------- tableaux globaux (partagés par fork) N = len(df) DAY = df["day"].to_numpy(np.int32); MONTH = df["month"].to_numpy(np.int32) RID = df["region"].map(rid_map).to_numpy(np.int32) X = df["x_km"].to_numpy(np.float64); Y = df["y_km"].to_numpy(np.float64) LNP = df["ln_price"].to_numpy(np.float64); P = df["amount"].to_numpy(np.float64) FLOOR = df["floor_area"].to_numpy(np.float64); LAND = df["land_area"].to_numpy(np.float64) AGE = df["age"].to_numpy(np.float64); LIEN = df["lien"].to_numpy(np.int8); ETAGES = df["etages"].to_numpy(np.float64) GENRE = df["genre"].to_numpy(np.int8) PROV = pd.factorize(df["id_provinc"])[0].astype(np.int32) MUN = pd.factorize(df["code_mun"])[0].astype(np.int32) VINT = df["role_vintage"].to_numpy(np.int8); LNROLE = np.log(df["role_vente"].to_numpy(np.float64)) PART = df["part_terrain"].to_numpy(np.float64) NEUF = df["neuf"].to_numpy(np.float64) # matrice des caractéristiques hédoniques (même ordre que FEATS) XF = np.column_stack([np.log(FLOOR), np.log(LAND), AGE, AGE ** 2 / 100, (LIEN == 1), (LIEN >= 2), (GENRE == 1), (GENRE == 2), (GENRE == 3), ETAGES, NEUF]).astype(np.float64) BETA_G = {a: np.array([MODELS[a]["beta_global"][f] for f in FEATS]) for a in ASOF} BETA_R = {a: np.array([[MODELS[a]["beta_region"][r][f] for f in FEATS] for r in regions]) for a in ASOF} ASOF_ARR = np.array(ASOF) TREE = cKDTree(np.column_stack([X, Y])) D_MAX, T_MAX_DAYS = 25.0, 731 KS = [1, 3, 5, 7, 10, 15, 20, 30, 50] DS = [0.5, 1, 2, 5, 10, 25] TS = [3, 6, 12, 24] SELS = ["geo", "sim", "mix", "adj"] ADJS = ["raw", "none", "conv", "hedo", "hedo_reg", "local", "gbm", "role"] WRS = [("equal", "median"), ("equal", "mean"), ("dist", "median"), ("dist", "mean"), ("rec", "median"), ("rec", "mean"), ("sim", "median"), ("sim", "mean"), ("all", "median"), ("all", "mean")] ANCHORS = [(3, 2, 6), (5, 5, 12), (10, 10, 24)] E1 = [dict(sel=s, d=d, T=T, k=k, adj=a, w="equal", recon="median") for s in ("geo", "sim", "mix") for d in DS for T in TS for k in KS for a in ("none", "hedo")] E2 = [dict(sel=s, d=d, T=T, k=k, adj=a, w=w, recon=r) for (k, d, T) in ANCHORS for s in SELS for a in ADJS for (w, r) in WRS] LOCAL_COLS = [0, 1, 2, 4, 5, 9] # ln_floor, ln_land, age, lien1, lien2, etages LAMBDA_LOCAL = 50.0 def wmedian(v, w): o = np.argsort(v); v = v[o]; w = w[o] c = np.cumsum(w); return v[np.searchsorted(c, 0.5 * c[-1])] def recon(vals, w, how): ok = ~np.isnan(vals) if not ok.any(): return np.nan, np.nan v = vals[ok]; ww = w[ok] sd = v.std() if len(v) > 1 else np.nan if how == "median": return (np.median(v) if ww.min() == ww.max() else wmedian(v, ww)), sd return np.log(np.sum(ww * np.exp(v)) / ww.sum()), sd def evaluate(i): """Toutes les configurations pour la cible i. Renvoie (est_E1, est_E2, sd_E2, n_avail, info).""" period = int(ASOF_ARR[ASOF_ARR <= MONTH[i]].max()) # dernier trimestre de ré-estimation ≤ mois de la cible e1 = np.full(len(E1), np.nan, np.float32); e2 = np.full(len(E2), np.nan, np.float32) s2 = np.full(len(E2), np.nan, np.float32); nav = np.zeros(len(DS) * len(TS), np.int32) cand = np.asarray(TREE.query_ball_point((X[i], Y[i]), D_MAX), dtype=np.int64) keep = (DAY[cand] < DAY[i]) & (DAY[cand] >= DAY[i] - T_MAX_DAYS) & (PROV[cand] != PROV[i]) cand = cand[keep] info = dict(idx=int(i), asof=period, n_cand=int(len(cand))) if len(cand) == 0: return e1, e2, s2, nav, info dist = np.hypot(X[cand] - X[i], Y[cand] - Y[i]) mago = (DAY[i] - DAY[cand]) / 30.4375 m_eval = MONTH[i]; lvl = IDX[m_eval - M_EVAL0] tadj = np.where(MONTH[cand] == m_eval, 0.0, lvl[RID[cand], m_eval - 1] - lvl[RID[cand], MONTH[cand]]) tadj = np.nan_to_num(tadj, nan=0.0) lnT = LNP[cand] + tadj dX = XF[i] - XF[cand] # cible − comparable # --- ajustements bg = BETA_G[period]; br = BETA_R[period][RID[i]] adj = {} adj["raw"] = LNP[cand] adj["none"] = lnT Pt = np.exp(lnT) d_floor = (FLOOR[i] - FLOOR[cand]) * 0.5 * Pt / FLOOR[cand] d_land = np.clip((LAND[i] - LAND[cand]) * 0.5 * Pt * PART[cand] / LAND[cand], -0.15 * Pt, 0.15 * Pt) d_age = np.clip(-0.005 * Pt * (AGE[i] - AGE[cand]), -0.2 * Pt, 0.2 * Pt) d_lien = np.where((LIEN[i] == 0) & (LIEN[cand] > 0), 0.05 * Pt, np.where((LIEN[i] > 0) & (LIEN[cand] == 0), -0.05 * Pt, 0.0)) tot = np.maximum(d_floor + d_land + d_age + d_lien, -0.6 * Pt) adj["conv"] = np.log(Pt + tot) gross_conv = (np.abs(d_floor) + np.abs(d_land) + np.abs(d_age) + np.abs(d_lien)) / Pt adj["hedo"] = lnT + dX @ bg gross_hedo = np.abs(dX * bg).sum(axis=1) adj["hedo_reg"] = lnT + dX @ br loc = dist <= 10.0 if loc.sum() >= 40: Xl = XF[cand][loc][:, LOCAL_COLS]; yl = lnT[loc] b0 = br[LOCAL_COLS] Xc = Xl - Xl.mean(axis=0); yc = yl - yl.mean() A = Xc.T @ Xc + LAMBDA_LOCAL * np.eye(len(LOCAL_COLS)) bl = np.linalg.solve(A, Xc.T @ yc + LAMBDA_LOCAL * b0) adj["local"] = lnT + dX[:, LOCAL_COLS] @ bl else: adj["local"] = adj["hedo_reg"] f = GBM[period] adj["gbm"] = lnT + (f[i] - f[cand]) same = (MUN[cand] == MUN[i]) & (VINT[cand] == VINT[i]) adj["role"] = np.where(same, LNROLE[i] + (lnT - LNROLE[cand]), np.nan) # --- scores de sélection dissim = np.sqrt((dX[:, 0] / 0.25) ** 2 + (dX[:, 1] / 0.6) ** 2 + (dX[:, 2] / 15) ** 2 + (LIEN[cand] != LIEN[i]) * 1.0 + (ETAGES[cand] != ETAGES[i]) * 0.5) score = {"geo": dist, "sim": dissim + dist / 25.0, "mix": dissim + dist / 2.0, "adj": gross_hedo + dist / 25.0} # --- pondérations wts = {"equal": np.ones(len(cand)), "dist": 1.0 / (0.25 + dist), "rec": np.exp(-mago / 6.0), "sim": 1.0 / (0.1 + gross_hedo)} wts["all"] = wts["dist"] * wts["rec"] * wts["sim"] # --- masques (d, T) et classements top = {} for di, d in enumerate(DS): for ti, T in enumerate(TS): m = (dist <= d) & (mago <= T) nav[di * len(TS) + ti] = int(m.sum()) if not m.any(): continue ids = np.flatnonzero(m) for s in SELS: sc = score[s][ids] if len(ids) > 50: part = np.argpartition(sc, 50)[:50]; part = part[np.argsort(sc[part])] else: part = np.argsort(sc) top[(s, d, T)] = ids[part] # --- E1 for c, cfg in enumerate(E1): t = top.get((cfg["sel"], cfg["d"], cfg["T"])) if t is None: continue sub = t[:cfg["k"]] e1[c], _ = recon(adj[cfg["adj"]][sub], wts["equal"][sub], "median") # --- E2 for c, cfg in enumerate(E2): t = top.get((cfg["sel"], cfg["d"], cfg["T"])) if t is None: continue sub = t[:cfg["k"]] e2[c], s2[c] = recon(adj[cfg["adj"]][sub], wts[cfg["w"]][sub], cfg["recon"]) return e1, e2, s2, nav, info def worker(chunk): return [evaluate(i) for i in chunk] if __name__ == "__main__": targets = df.index[df["period"].isin(["val", "test"])].to_numpy() if os.environ.get("CV_SAMPLE"): rng = np.random.default_rng(1); targets = np.sort(rng.choice(targets, int(os.environ["CV_SAMPLE"]), replace=False)) SUF = "" if os.environ.get("CV_SHARD"): # répartition sur plusieurs nœuds : CV_SHARD=i/N → cibles i, i+N, i+2N, … i_sh, n_sh = map(int, os.environ["CV_SHARD"].split("/")); targets = targets[i_sh::n_sh]; SUF = f".shard{i_sh}of{n_sh}" log(f"{len(targets):,} cibles ; {len(E1)} configurations E1, {len(E2)} E2 ; {N:,} ventes dans le bassin") pd.DataFrame(E1).to_csv(OUT / "configs_E1.csv", index_label="cfg") pd.DataFrame(E2).to_csv(OUT / "configs_E2.csv", index_label="cfg") est1 = np.full((len(targets), len(E1)), np.nan, np.float32) est2 = np.full((len(targets), len(E2)), np.nan, np.float32) sd2 = np.full((len(targets), len(E2)), np.nan, np.float32) nav = np.zeros((len(targets), len(DS) * len(TS)), np.int32) infos = [] chunks = [targets[j:j + 400] for j in range(0, len(targets), 400)] t0 = time.time(); done = 0 with mp.get_context(os.environ.get("CV_MP", "fork")).Pool(int(os.environ.get("CV_PROCS", 20))) as pool: for res in pool.imap(worker, chunks): for (e1, e2, s2, nv, info) in res: est1[done] = e1; est2[done] = e2; sd2[done] = s2; nav[done] = nv; infos.append(info); done += 1 if done % 20000 < 400: log(f" {done:,}/{len(targets):,} cibles ({time.time() - t0:.0f} s)") np.save(OUT / f"est_E1{SUF}.npy", est1); np.save(OUT / f"est_E2{SUF}.npy", est2); np.save(OUT / f"sd_E2{SUF}.npy", sd2); np.save(OUT / f"n_avail{SUF}.npy", nav) T = df.loc[targets, ["idx", "id", "date", "period", "amount", "ln_price", "region", "region_nom", "mrc", "mrc_nom", "code_mun", "mun_nom", "unite_voisinage", "floor_area", "land_area", "age", "lien", "etages", "role_vente", "month", "lat", "lng"]].reset_index(drop=True) T["n_cand"] = [x["n_cand"] for x in infos]; T["asof"] = [x["asof"] for x in infos] T["n_2km_12m"] = nav[:, DS.index(2) * len(TS) + TS.index(12)] T["n_5km_12m"] = nav[:, DS.index(5) * len(TS) + TS.index(12)] save_parquet(T, f"targets{SUF}.parquet") pd.DataFrame({"DS": [DS], "TS": [TS], "KS": [KS]}).to_json(OUT / "grid_meta.json") log(f"rétro-évaluation terminée en {time.time() - t0:.0f} s")