#!/usr/bin/env python3
"""Chapitre 21 — coût neuf (remplacement) du sujet du cas 1 (Gatineau, DONNÉES SYNTHÉTIQUES) et indexation IPCB (données réelles).

Imprime tous les chiffres cités dans le chapitre 21 :
  * devis par composantes (data/synthetic/cas1_gatineau_couts.csv, 38 postes) regroupé par corps de métier ;
  * coûts directs, indirects, incitatif de l'entrepreneur (traite.cout.cout_neuf), coût neuf total, $/m² et $/pi² ;
  * comparaison unitaire avec la fourchette Altus 2026 (155-215 $/pi², Montréal, unifamiliale) ;
  * indexation par ratio d'indices IPCB Montréal (Statistique Canada 18-10-0289-01, données réelles dans data/raw) ;
  * ajustement de taille (élasticité) ;
  * analyse de sensibilité (taux d'indirects, incitatif, indice, élasticité).

Usage (racine du dépôt) : PYTHONPATH=analysis/python python3 analysis/python/cases/ch21_cout_neuf.py
Toutes les fonctions financières proviennent de analysis/python/traite/cout.py (aucune réécriture)."""
from __future__ import annotations
import sys, json, csv, pathlib
import numpy as np, pandas as pd
sys.path.insert(0, str(pathlib.Path(__file__).resolve().parents[1]))
from traite import cout as ct

ROOT = pathlib.Path(__file__).resolve().parents[3]
D = ROOT / "data/synthetic"
RAW = ROOT / "data/raw/statcan_18100289/18100289.csv"

PI2_PAR_M2 = 10.7639
TAUX_INCITATIF = 0.10          # hypothèse pédagogique (voir texte, §21.6) ; assiette : directs + indirects
ALTUS_MIN, ALTUS_MAX = 155.0, 215.0   # $/pi², Altus Canadian Cost Guide 2026, Montréal, unifamiliale (tel que relayé)

# Regroupement des 38 postes du devis en sept corps de métier (convention du chapitre)
CORPS = {
    "Excavation et remblai": "Fondations et site", "Fondations béton (semelles, murs)": "Fondations et site",
    "Dalle de sous-sol": "Fondations et site", "Drain français et imperméabilisation": "Fondations et site",
    "Raccordements aqueduc/égout": "Fondations et site",
    "Charpente murs et planchers": "Structure et charpente", "Fermes de toit": "Structure et charpente",
    "Escaliers": "Structure et charpente",
    "Couverture bardeaux d'asphalte": "Enveloppe", "Revêtement extérieur (brique 40 % / fibrociment 60 %)": "Enveloppe",
    "Isolation et pare-vapeur": "Enveloppe", "Fenêtres (PVC double vitrage)": "Enveloppe", "Portes extérieures": "Enveloppe",
    "Porte de garage": "Enveloppe", "Gouttières et solins": "Enveloppe",
    "Gypse et tirage de joints": "Finition intérieure", "Peinture intérieure": "Finition intérieure",
    "Planchers bois franc": "Finition intérieure", "Planchers céramique": "Finition intérieure",
    "Planchers flottants (sous-sol)": "Finition intérieure", "Portes intérieures et moulures": "Finition intérieure",
    "Armoires de cuisine et comptoirs": "Finition intérieure", "Vanités et accessoires (2 sdb)": "Finition intérieure",
    "Plomberie brute et finition": "Mécanique et électricité", "Chauffe-eau": "Mécanique et électricité",
    "Électricité (entrée 200 A, filage, appareils)": "Mécanique et électricité", "Ventilation (VRC) et conduits": "Mécanique et électricité",
    "Thermopompe centrale et plinthes": "Mécanique et électricité", "Éclairage et finition électrique": "Mécanique et électricité",
    "Balcon et perron": "Aménagements extérieurs", "Entrée asphaltée": "Aménagements extérieurs",
    "Aménagement paysager de base": "Aménagements extérieurs", "Nettoyage et conteneurs": "Aménagements extérieurs",
    "Permis et études (arpentage, sol)": "Coûts indirects", "Supervision de chantier (salaires)": "Coûts indirects",
    "Assurances chantier et cautionnement": "Coûts indirects", "Plans et honoraires professionnels": "Coûts indirects",
    "Financement intérimaire": "Coûts indirects",
}
ORDRE_CORPS = ["Fondations et site", "Structure et charpente", "Enveloppe", "Finition intérieure",
               "Mécanique et électricité", "Aménagements extérieurs", "Coûts indirects"]


def charger() -> tuple[pd.DataFrame, dict]:
    devis = pd.read_csv(D / "cas1_gatineau_couts.csv")
    devis["corps"] = devis["poste"].map(CORPS)
    assert devis["corps"].notna().all(), "poste sans corps de métier"
    sujet = json.loads((D / "cas1_gatineau_sujet.json").read_text())
    return devis, sujet


def regrouper(devis: pd.DataFrame) -> pd.DataFrame:
    g = devis.groupby("corps", sort=False)["montant"].agg(["sum", "count"]).reindex(ORDRE_CORPS)
    g.columns = ["montant", "postes"]
    g["part_direct_indirect"] = g["montant"] / g["montant"].sum()
    return g


def cout_neuf_sujet(devis: pd.DataFrame, taux_incitatif: float = TAUX_INCITATIF, base: str = "direct+indirect") -> dict:
    """Coût neuf du sujet : les indirects viennent du devis (5 postes), l'incitatif est appliqué par traite.cout.cout_neuf."""
    direct = float(devis.loc[~devis["indirect"], "montant"].sum())
    indirects = float(devis.loc[devis["indirect"], "montant"].sum())
    taux_ind = indirects / direct
    cn = ct.cout_neuf(direct, taux_ind, taux_incitatif, base)
    cn["taux_indirects"] = taux_ind
    cn["taux_incitatif"] = taux_incitatif
    return cn


def superficies(sujet: dict, devis: pd.DataFrame) -> dict:
    """Superficie habitable (JSON du sujet) et aire brute des étages (hypothèse : empreinte = aire de la dalle de sous-sol × 2 étages)."""
    hab = float(sujet["superficie_habitable_m2"])
    empreinte = float(devis.loc[devis["poste"] == "Dalle de sous-sol", "quantite"].iloc[0])
    brute = 2.0 * empreinte     # cottage à deux étages hors sol (type « à étages »)
    return {"habitable_m2": hab, "habitable_pi2": hab * PI2_PAR_M2, "empreinte_m2": empreinte,
            "brute_etages_m2": brute, "brute_etages_pi2": brute * PI2_PAR_M2}


# ----------------------------------------------------------------------------- IPCB (données réelles)

def ipcb_montreal(types: tuple[str, ...] = ("Bâtiments résidentiels [621]", "Maison individuelle [621112]",
                                           "Immeuble à hauteur restreinte (moins de cinq étages) [62111]",
                                           "Tour d'habitations (cinq étages ou plus) [621111]",
                                           "Bâtiments non résidentiels [622]")) -> pd.DataFrame:
    """Séries trimestrielles IPCB (2023 = 100), Montréal, agrégat des divisions, lues dans le CSV brut de Statistique Canada."""
    rows = []
    with open(RAW, encoding="utf-8-sig") as f:
        for r in csv.DictReader(f, delimiter=";"):
            if r["GÉO"] == "Montréal, Québec" and r["Division"] == "Agrégat des divisions" and r["Type d'immeuble"] in types and r["VALEUR"]:
                rows.append({"periode": r["PÉRIODE DE RÉFÉRENCE"], "type": r["Type d'immeuble"], "indice": float(r["VALEUR"]), "vecteur": r["VECTEUR"]})
    df = pd.DataFrame(rows)
    df["date"] = pd.PeriodIndex(df["periode"].str.replace("-01", "Q1").str.replace("-04", "Q2").str.replace("-07", "Q3").str.replace("-10", "Q4"), freq="Q").to_timestamp()
    return df.sort_values(["type", "date"]).reset_index(drop=True)


def indice(df: pd.DataFrame, type_: str, periode: str) -> float:
    v = df[(df["type"] == type_) & (df["periode"] == periode)]["indice"]
    if v.empty:
        raise KeyError(f"{type_} {periode} absent")
    return float(v.iloc[0])


# ----------------------------------------------------------------------------- sensibilité

def grille_sensibilite(direct: float, taux_ind: np.ndarray, taux_inc: np.ndarray) -> np.ndarray:
    """Coût neuf (assiette directs + indirects) pour chaque couple (taux d'indirects, taux d'incitatif)."""
    return np.array([[ct.cout_neuf(direct, ti, tc)["total"] for tc in taux_inc] for ti in taux_ind])


def tornade(direct: float, taux_ind0: float, taux_inc0: float, ratio_idx: tuple[float, float], elast: tuple[float, float],
            d_ind: float = 0.05, d_inc: float = 0.05) -> list[dict]:
    """Variation du coût neuf quand un seul paramètre bouge (bas / haut), les autres restant à leur valeur centrale."""
    base = ct.cout_neuf(direct, taux_ind0, taux_inc0)["total"]
    out = [
        {"param": "Taux d'indirects", "bas": ct.cout_neuf(direct, taux_ind0 - d_ind, taux_inc0)["total"], "haut": ct.cout_neuf(direct, taux_ind0 + d_ind, taux_inc0)["total"],
         "etiquette": f"{(taux_ind0-d_ind)*100:.1f} % → {(taux_ind0+d_ind)*100:.1f} %"},
        {"param": "Taux d'incitatif", "bas": ct.cout_neuf(direct, taux_ind0, taux_inc0 - d_inc)["total"], "haut": ct.cout_neuf(direct, taux_ind0, taux_inc0 + d_inc)["total"],
         "etiquette": f"{(taux_inc0-d_inc)*100:.0f} % → {(taux_inc0+d_inc)*100:.0f} %"},
        {"param": "Indice de coûts (date)", "bas": ct.indexer_cout(base, 1.0, ratio_idx[0]), "haut": ct.indexer_cout(base, 1.0, ratio_idx[1]),
         "etiquette": f"×{ratio_idx[0]:.3f} → ×{ratio_idx[1]:.3f}"},
        {"param": "Élasticité de taille", "bas": base * ct.cout_par_superficie_ajuste(1.0, 140.0, 165.0, elast[0]) / ct.cout_par_superficie_ajuste(1.0, 140.0, 165.0, -0.10),
         "haut": base * ct.cout_par_superficie_ajuste(1.0, 140.0, 165.0, elast[1]) / ct.cout_par_superficie_ajuste(1.0, 140.0, 165.0, -0.10),
         "etiquette": f"{elast[0]:.2f} → {elast[1]:.2f}"},
    ]
    for o in out:
        o["base"] = base
    return out


def main() -> dict:
    devis, sujet = charger()
    g = regrouper(devis)
    cn = cout_neuf_sujet(devis)
    sup = superficies(sujet, devis)
    f = lambda x: f"{x:,.2f}".replace(",", " ")
    print("== DONNÉES SYNTHÉTIQUES — sujet cas 1 :", sujet["type"], sujet["annee_construction"], sujet["superficie_habitable_m2"], "m² hab.,", sujet["secteur"])
    print(f"== Devis : {len(devis)} postes ; {int((~devis['indirect']).sum())} directs, {int(devis['indirect'].sum())} indirects")
    print("== Regroupement par corps de métier ($) :")
    for k, r in g.iterrows():
        print(f"   {k:28s} {int(r['postes']):2d} postes  {f(r['montant']):>14s}  {r['part_direct_indirect']*100:5.1f} %")
    print(f"== Coûts directs D        = {f(cn['direct'])} $")
    print(f"== Coûts indirects I      = {f(cn['indirects'])} $  (I/D = {cn['taux_indirects']*100:.2f} % ; I/(D+I) = {cn['indirects']/(cn['direct']+cn['indirects'])*100:.2f} %)")
    print(f"== Sous-total D + I       = {f(cn['direct']+cn['indirects'])} $")
    print(f"== Incitatif {TAUX_INCITATIF*100:.0f} % (D+I)   = {f(cn['incitatif'])} $")
    print(f"== COÛT NEUF C_N          = {f(cn['total'])} $  (arrondi {round(cn['total']/1000)*1000:,.0f} $)".replace(",", " "))
    print(f"   part directs {cn['direct']/cn['total']*100:.1f} %, indirects {cn['indirects']/cn['total']*100:.1f} %, incitatif {cn['incitatif']/cn['total']*100:.1f} %")
    alt = cout_neuf_sujet(devis, TAUX_INCITATIF, "direct")
    print(f"== Variante assiette 'direct' : incitatif {f(alt['incitatif'])} $, C_N {f(alt['total'])} $ (écart {f(cn['total']-alt['total'])} $)")
    print(f"== Superficies : habitable {sup['habitable_m2']:.0f} m² = {sup['habitable_pi2']:.0f} pi² ; empreinte {sup['empreinte_m2']:.0f} m² ; aire brute des étages {sup['brute_etages_m2']:.0f} m² = {sup['brute_etages_pi2']:.0f} pi²")
    um_hab = cn["total"] / sup["habitable_m2"]; um_brut = cn["total"] / sup["brute_etages_m2"]
    print(f"== C_N / m² habitable = {um_hab:,.0f} $/m² = {um_hab/PI2_PAR_M2:,.1f} $/pi²".replace(",", " "))
    print(f"== C_N / m² brut      = {um_brut:,.0f} $/m² = {um_brut/PI2_PAR_M2:,.1f} $/pi²".replace(",", " "))
    di = cn["direct"] + cn["indirects"]
    print(f"== (D+I) / m² habitable = {di/sup['habitable_m2']:,.0f} $/m² = {di/sup['habitable_m2']/PI2_PAR_M2:,.1f} $/pi² ; (D+I)/m² brut = {di/sup['brute_etages_m2']:,.0f} $/m² = {di/sup['brute_etages_m2']/PI2_PAR_M2:,.1f} $/pi²".replace(",", " "))
    # comparaison unitaire Altus
    alt_mid = 0.5 * (ALTUS_MIN + ALTUS_MAX)
    print(f"== Comparaison unitaire Altus 2026 ({ALTUS_MIN:.0f}-{ALTUS_MAX:.0f} $/pi², milieu {alt_mid:.0f}) :")
    for nom, s_pi2 in [("aire brute des étages", sup["brute_etages_pi2"]), ("superficie habitable", sup["habitable_pi2"])]:
        print(f"   sur {nom} ({s_pi2:.0f} pi²) : {ALTUS_MIN*s_pi2:,.0f} – {ALTUS_MAX*s_pi2:,.0f} $, milieu {alt_mid*s_pi2:,.0f} $".replace(",", " "))
    print(f"   milieu Altus en $/m² : {alt_mid*PI2_PAR_M2:,.0f} $/m²".replace(",", " "))
    # indexation IPCB
    ipcb = ipcb_montreal()
    res, mi = "Bâtiments résidentiels [621]", "Maison individuelle [621112]"
    print("== IPCB Montréal (2023 = 100), données réelles Statistique Canada 18-10-0289-01 :")
    for t in ipcb["type"].unique():
        s = ipcb[ipcb["type"] == t]
        print(f"   {t:60s} {len(s)} trimestres, {s['periode'].iloc[0]} → {s['periode'].iloc[-1]} ; premier {s['indice'].iloc[0]}, dernier {s['indice'].iloc[-1]} ({s['vecteur'].iloc[0]})")
    pts = {p: indice(ipcb, mi, p) for p in ["2017-01", "2019-07", "2021-07", "2023-07", "2025-07", "2026-04", "2026-07"]}
    print("   Maison individuelle :", pts)
    ptr = {p: indice(ipcb, res, p) for p in ["2017-01", "2019-07", "2021-07", "2023-07", "2025-07", "2026-07"]}
    print("   Bâtiments résidentiels :", ptr)
    var_an = indice(ipcb, res, "2026-07") / indice(ipcb, res, "2025-07") - 1
    var_an_mi = indice(ipcb, mi, "2026-07") / indice(ipcb, mi, "2025-07") - 1
    print(f"   Variation annuelle T3 2025 → T3 2026 : résidentiel {var_an*100:.2f} %, maison individuelle {var_an_mi*100:.2f} %")
    r_2017 = indice(ipcb, mi, "2026-07") / indice(ipcb, mi, "2017-01")
    r_2019 = indice(ipcb, mi, "2026-07") / indice(ipcb, mi, "2019-07")
    r_2021 = indice(ipcb, mi, "2026-07") / indice(ipcb, mi, "2021-07")
    print(f"   Ratios maison individuelle : T3 2026 / T1 2017 = {r_2017:.4f} ; / T3 2019 = {r_2019:.4f} ; / T3 2021 = {r_2021:.4f}")
    # taux annuel composé 2017T1→2026T3 (9,5 ans)
    n_ans = (2026 + 6 / 12) - (2017 + 0 / 12)
    print(f"   Croissance annuelle composée T1 2017 → T3 2026 ({n_ans:.1f} ans) : {(r_2017**(1/n_ans)-1)*100:.2f} %/an")
    # exemple d'indexation : coût du sujet exprimé en dollars du T3 2019 et du T3 2021 (rétropolation)
    c_2019 = ct.indexer_cout(cn["total"], indice(ipcb, mi, "2026-07"), indice(ipcb, mi, "2019-07"))
    c_2021 = ct.indexer_cout(cn["total"], indice(ipcb, mi, "2026-07"), indice(ipcb, mi, "2021-07"))
    print(f"   C_N du sujet en dollars du T3 2019 : {c_2019:,.0f} $ ; du T3 2021 : {c_2021:,.0f} $".replace(",", " "))
    # exemple direct : un devis de 2019T3 de 290 000 $ porté au T3 2026
    c0 = 290_000.0
    print(f"   Exemple indexation : {c0:,.0f} $ (T3 2019) × {indice(ipcb, mi, '2026-07')}/{indice(ipcb, mi, '2019-07')} = {ct.indexer_cout(c0, indice(ipcb, mi, '2019-07'), indice(ipcb, mi, '2026-07')):,.0f} $".replace(",", " "))
    # erreur de série : non résidentiel et tours comme mauvais substituts
    for t in ["Bâtiments non résidentiels [622]", "Tour d'habitations (cinq étages ou plus) [621111]", "Immeuble à hauteur restreinte (moins de cinq étages) [62111]"]:
        r = indice(ipcb, t, "2026-07") / indice(ipcb, t, "2019-07")
        print(f"   Ratio T3 2026 / T3 2019 avec {t} : {r:.4f} (écart vs maison individuelle {(r/r_2019-1)*100:+.1f} %)")
    # ajustement de taille
    print("== Ajustement de taille (élasticité) :")
    for e in [-0.05, -0.10, -0.15, -0.20]:
        u = ct.cout_par_superficie_ajuste(um_hab, 165.0, 140.0, e)   # coût unitaire d'une maison de 140 m² déduit du sujet
        u2 = ct.cout_par_superficie_ajuste(um_hab, 165.0, 200.0, e)
        print(f"   e = {e:+.2f} : 140 m² → {u:,.0f} $/m² ({u*140:,.0f} $) ; 200 m² → {u2:,.0f} $/m² ({u2*200:,.0f} $)".replace(",", " "))
    # exemple du texte : coût unitaire de référence 2 300 $/m² observé sur 200 m², appliqué à 165 m²
    u_ref = 2300.0
    for e in [0.0, -0.10]:
        u = ct.cout_par_superficie_ajuste(u_ref, 200.0, 165.0, e)
        print(f"   Référence {u_ref:.0f} $/m² (200 m²) → 165 m² avec e = {e:+.2f} : {u:,.1f} $/m², coût {u*165:,.0f} $".replace(",", " "))
    # sensibilité
    print("== Sensibilité (assiette directs + indirects) :")
    ti = np.array([0.12, 0.15, 0.18, 0.2095, 0.24, 0.27]); tc = np.array([0.05, 0.08, 0.10, 0.12, 0.15])
    G = grille_sensibilite(cn["direct"], ti, tc)
    print("   indirects \\ incitatif", [f"{x*100:.0f} %" for x in tc])
    for i, x in enumerate(ti):
        print(f"   {x*100:5.1f} %  ", [f"{v:,.0f}".replace(",", " ") for v in G[i]])
    print(f"   Étendue de la grille : {G.min():,.0f} – {G.max():,.0f} $ ({(G.max()/G.min()-1)*100:.1f} %)".replace(",", " "))
    torn = tornade(cn["direct"], cn["taux_indirects"], TAUX_INCITATIF, (1 / (1 + var_an_mi), 1 + var_an_mi), (-0.05, -0.20))
    for o in torn:
        print(f"   {o['param']:24s} {o['etiquette']:22s} bas {o['bas']:,.0f}  haut {o['haut']:,.0f}  (base {o['base']:,.0f}) ; étendue {(o['haut']-o['bas'])/o['base']*100:+.1f} %".replace(",", " "))
    return {"devis": devis, "sujet": sujet, "groupes": g, "cout_neuf": cn, "superficies": sup, "ipcb": ipcb, "grille": (ti, tc, G), "tornade": torn}


if __name__ == "__main__":
    main()