TeX 52.9%
Python 45.7%
Shell 1.4%
1# -*- coding: utf-8 -*-2"""04 — Moteur de rétro-évaluation.3Pour chaque vente cible (2024 = validation, 2025-2026 = test), le prix est masqué ; seules les ventes4conclues STRICTEMENT AVANT la date cible (et l'indice/les modèles connus à cette date) sont utilisées.5On évalue d'un coup toutes les configurations (sélection × rayon × fenêtre × nombre × ajustement ×6pondération × réconciliation) sur exactement les mêmes cibles.7Sorties : targets.parquet, configs_E1.csv/est_E1.npy (surface k,d,T), configs_E2.csv/est_E2.npy/sd_E2.npy8(factoriel sélection × ajustement × pondération × réconciliation), n_avail.npy."""9import os10for _v in ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS", "VECLIB_MAXIMUM_THREADS"):11 os.environ[_v] = "1" # un seul fil BLAS par processus (fork + BLAS multi-fils = interblocage sur macOS)12import multiprocessing as mp, time13from scipy.spatial import cKDTree14from common import *1516df = load_parquet("sales.parquet")17meta = load_json("index_meta.json"); regions = meta["regions"]; M_EVAL0 = meta["M_EVAL0"]18rid_map = {r: i for i, r in enumerate(regions)}19IDX = np.load(OUT / "index_asof.npy")20ASOF = load_json("asof.json")["asof"]21MODELS = {a: load_json(f"models_{a}.json") for a in ASOF}22FEATS = MODELS[ASOF[0]]["feats"]23GBM = {a: np.load(OUT / f"gbm_pred_{a}.npy") for a in ASOF}2425# ----------------------------------------------------------------- tableaux globaux (partagés par fork)26N = len(df)27DAY = df["day"].to_numpy(np.int32); MONTH = df["month"].to_numpy(np.int32)28RID = df["region"].map(rid_map).to_numpy(np.int32)29X = df["x_km"].to_numpy(np.float64); Y = df["y_km"].to_numpy(np.float64)30LNP = df["ln_price"].to_numpy(np.float64); P = df["amount"].to_numpy(np.float64)31FLOOR = df["floor_area"].to_numpy(np.float64); LAND = df["land_area"].to_numpy(np.float64)32AGE = df["age"].to_numpy(np.float64); LIEN = df["lien"].to_numpy(np.int8); ETAGES = df["etages"].to_numpy(np.float64)33GENRE = df["genre"].to_numpy(np.int8)34PROV = pd.factorize(df["id_provinc"])[0].astype(np.int32)35MUN = pd.factorize(df["code_mun"])[0].astype(np.int32)36VINT = df["role_vintage"].to_numpy(np.int8); LNROLE = np.log(df["role_vente"].to_numpy(np.float64))37PART = df["part_terrain"].to_numpy(np.float64)38NEUF = df["neuf"].to_numpy(np.float64)39# matrice des caractéristiques hédoniques (même ordre que FEATS)40XF = np.column_stack([np.log(FLOOR), np.log(LAND), AGE, AGE ** 2 / 100, (LIEN == 1), (LIEN >= 2),41 (GENRE == 1), (GENRE == 2), (GENRE == 3), ETAGES, NEUF]).astype(np.float64)42BETA_G = {a: np.array([MODELS[a]["beta_global"][f] for f in FEATS]) for a in ASOF}43BETA_R = {a: np.array([[MODELS[a]["beta_region"][r][f] for f in FEATS] for r in regions]) for a in ASOF}44ASOF_ARR = np.array(ASOF)45TREE = cKDTree(np.column_stack([X, Y]))4647D_MAX, T_MAX_DAYS = 25.0, 73148KS = [1, 3, 5, 7, 10, 15, 20, 30, 50]49DS = [0.5, 1, 2, 5, 10, 25]50TS = [3, 6, 12, 24]51SELS = ["geo", "sim", "mix", "adj"]52ADJS = ["raw", "none", "conv", "hedo", "hedo_reg", "local", "gbm", "role"]53WRS = [("equal", "median"), ("equal", "mean"), ("dist", "median"), ("dist", "mean"), ("rec", "median"), ("rec", "mean"),54 ("sim", "median"), ("sim", "mean"), ("all", "median"), ("all", "mean")]55ANCHORS = [(3, 2, 6), (5, 5, 12), (10, 10, 24)]56E1 = [dict(sel=s, d=d, T=T, k=k, adj=a, w="equal", recon="median")57 for s in ("geo", "sim", "mix") for d in DS for T in TS for k in KS for a in ("none", "hedo")]58E2 = [dict(sel=s, d=d, T=T, k=k, adj=a, w=w, recon=r)59 for (k, d, T) in ANCHORS for s in SELS for a in ADJS for (w, r) in WRS]60LOCAL_COLS = [0, 1, 2, 4, 5, 9] # ln_floor, ln_land, age, lien1, lien2, etages61LAMBDA_LOCAL = 50.06263def wmedian(v, w):64 o = np.argsort(v); v = v[o]; w = w[o]65 c = np.cumsum(w); return v[np.searchsorted(c, 0.5 * c[-1])]6667def recon(vals, w, how):68 ok = ~np.isnan(vals)69 if not ok.any(): return np.nan, np.nan70 v = vals[ok]; ww = w[ok]71 sd = v.std() if len(v) > 1 else np.nan72 if how == "median":73 return (np.median(v) if ww.min() == ww.max() else wmedian(v, ww)), sd74 return np.log(np.sum(ww * np.exp(v)) / ww.sum()), sd7576def evaluate(i):77 """Toutes les configurations pour la cible i. Renvoie (est_E1, est_E2, sd_E2, n_avail, info)."""78 period = int(ASOF_ARR[ASOF_ARR <= MONTH[i]].max()) # dernier trimestre de ré-estimation ≤ mois de la cible79 e1 = np.full(len(E1), np.nan, np.float32); e2 = np.full(len(E2), np.nan, np.float32)80 s2 = np.full(len(E2), np.nan, np.float32); nav = np.zeros(len(DS) * len(TS), np.int32)81 cand = np.asarray(TREE.query_ball_point((X[i], Y[i]), D_MAX), dtype=np.int64)82 keep = (DAY[cand] < DAY[i]) & (DAY[cand] >= DAY[i] - T_MAX_DAYS) & (PROV[cand] != PROV[i])83 cand = cand[keep]84 info = dict(idx=int(i), asof=period, n_cand=int(len(cand)))85 if len(cand) == 0: return e1, e2, s2, nav, info86 dist = np.hypot(X[cand] - X[i], Y[cand] - Y[i])87 mago = (DAY[i] - DAY[cand]) / 30.437588 m_eval = MONTH[i]; lvl = IDX[m_eval - M_EVAL0]89 tadj = np.where(MONTH[cand] == m_eval, 0.0, lvl[RID[cand], m_eval - 1] - lvl[RID[cand], MONTH[cand]])90 tadj = np.nan_to_num(tadj, nan=0.0)91 lnT = LNP[cand] + tadj92 dX = XF[i] - XF[cand] # cible − comparable93 # --- ajustements94 bg = BETA_G[period]; br = BETA_R[period][RID[i]]95 adj = {}96 adj["raw"] = LNP[cand]97 adj["none"] = lnT98 Pt = np.exp(lnT)99 d_floor = (FLOOR[i] - FLOOR[cand]) * 0.5 * Pt / FLOOR[cand]100 d_land = np.clip((LAND[i] - LAND[cand]) * 0.5 * Pt * PART[cand] / LAND[cand], -0.15 * Pt, 0.15 * Pt)101 d_age = np.clip(-0.005 * Pt * (AGE[i] - AGE[cand]), -0.2 * Pt, 0.2 * Pt)102 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))103 tot = np.maximum(d_floor + d_land + d_age + d_lien, -0.6 * Pt)104 adj["conv"] = np.log(Pt + tot)105 gross_conv = (np.abs(d_floor) + np.abs(d_land) + np.abs(d_age) + np.abs(d_lien)) / Pt106 adj["hedo"] = lnT + dX @ bg107 gross_hedo = np.abs(dX * bg).sum(axis=1)108 adj["hedo_reg"] = lnT + dX @ br109 loc = dist <= 10.0110 if loc.sum() >= 40:111 Xl = XF[cand][loc][:, LOCAL_COLS]; yl = lnT[loc]112 b0 = br[LOCAL_COLS]113 Xc = Xl - Xl.mean(axis=0); yc = yl - yl.mean()114 A = Xc.T @ Xc + LAMBDA_LOCAL * np.eye(len(LOCAL_COLS))115 bl = np.linalg.solve(A, Xc.T @ yc + LAMBDA_LOCAL * b0)116 adj["local"] = lnT + dX[:, LOCAL_COLS] @ bl117 else:118 adj["local"] = adj["hedo_reg"]119 f = GBM[period]120 adj["gbm"] = lnT + (f[i] - f[cand])121 same = (MUN[cand] == MUN[i]) & (VINT[cand] == VINT[i])122 adj["role"] = np.where(same, LNROLE[i] + (lnT - LNROLE[cand]), np.nan)123 # --- scores de sélection124 dissim = np.sqrt((dX[:, 0] / 0.25) ** 2 + (dX[:, 1] / 0.6) ** 2 + (dX[:, 2] / 15) ** 2125 + (LIEN[cand] != LIEN[i]) * 1.0 + (ETAGES[cand] != ETAGES[i]) * 0.5)126 score = {"geo": dist, "sim": dissim + dist / 25.0, "mix": dissim + dist / 2.0, "adj": gross_hedo + dist / 25.0}127 # --- pondérations128 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)}129 wts["all"] = wts["dist"] * wts["rec"] * wts["sim"]130 # --- masques (d, T) et classements131 top = {}132 for di, d in enumerate(DS):133 for ti, T in enumerate(TS):134 m = (dist <= d) & (mago <= T)135 nav[di * len(TS) + ti] = int(m.sum())136 if not m.any(): continue137 ids = np.flatnonzero(m)138 for s in SELS:139 sc = score[s][ids]140 if len(ids) > 50:141 part = np.argpartition(sc, 50)[:50]; part = part[np.argsort(sc[part])]142 else:143 part = np.argsort(sc)144 top[(s, d, T)] = ids[part]145 # --- E1146 for c, cfg in enumerate(E1):147 t = top.get((cfg["sel"], cfg["d"], cfg["T"]))148 if t is None: continue149 sub = t[:cfg["k"]]150 e1[c], _ = recon(adj[cfg["adj"]][sub], wts["equal"][sub], "median")151 # --- E2152 for c, cfg in enumerate(E2):153 t = top.get((cfg["sel"], cfg["d"], cfg["T"]))154 if t is None: continue155 sub = t[:cfg["k"]]156 e2[c], s2[c] = recon(adj[cfg["adj"]][sub], wts[cfg["w"]][sub], cfg["recon"])157 return e1, e2, s2, nav, info158159def worker(chunk):160 return [evaluate(i) for i in chunk]161162if __name__ == "__main__":163 targets = df.index[df["period"].isin(["val", "test"])].to_numpy()164 if os.environ.get("CV_SAMPLE"):165 rng = np.random.default_rng(1); targets = np.sort(rng.choice(targets, int(os.environ["CV_SAMPLE"]), replace=False))166 SUF = ""167 if os.environ.get("CV_SHARD"): # répartition sur plusieurs nœuds : CV_SHARD=i/N → cibles i, i+N, i+2N, …168 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}"169 log(f"{len(targets):,} cibles ; {len(E1)} configurations E1, {len(E2)} E2 ; {N:,} ventes dans le bassin")170 pd.DataFrame(E1).to_csv(OUT / "configs_E1.csv", index_label="cfg")171 pd.DataFrame(E2).to_csv(OUT / "configs_E2.csv", index_label="cfg")172 est1 = np.full((len(targets), len(E1)), np.nan, np.float32)173 est2 = np.full((len(targets), len(E2)), np.nan, np.float32)174 sd2 = np.full((len(targets), len(E2)), np.nan, np.float32)175 nav = np.zeros((len(targets), len(DS) * len(TS)), np.int32)176 infos = []177 chunks = [targets[j:j + 400] for j in range(0, len(targets), 400)]178 t0 = time.time(); done = 0179 with mp.get_context(os.environ.get("CV_MP", "fork")).Pool(int(os.environ.get("CV_PROCS", 20))) as pool:180 for res in pool.imap(worker, chunks):181 for (e1, e2, s2, nv, info) in res:182 est1[done] = e1; est2[done] = e2; sd2[done] = s2; nav[done] = nv; infos.append(info); done += 1183 if done % 20000 < 400:184 log(f" {done:,}/{len(targets):,} cibles ({time.time() - t0:.0f} s)")185 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)186 T = df.loc[targets, ["idx", "id", "date", "period", "amount", "ln_price", "region", "region_nom", "mrc", "mrc_nom", "code_mun",187 "mun_nom", "unite_voisinage", "floor_area", "land_area", "age", "lien", "etages", "role_vente", "month", "lat", "lng"]].reset_index(drop=True)188 T["n_cand"] = [x["n_cand"] for x in infos]; T["asof"] = [x["asof"] for x in infos]189 T["n_2km_12m"] = nav[:, DS.index(2) * len(TS) + TS.index(12)]190 T["n_5km_12m"] = nav[:, DS.index(5) * len(TS) + TS.index(12)]191 save_parquet(T, f"targets{SUF}.parquet")192 pd.DataFrame({"DS": [DS], "TS": [TS], "KS": [KS]}).to_json(OUT / "grid_meta.json")193 log(f"rétro-évaluation terminée en {time.time() - t0:.0f} s")194