#!/usr/bin/env python3
"""Chapitre 41 — Plus haute et meilleure utilisation, terrains et développement : calculs de l'exemple intégral.

Sujet : le terrain de 620 m² du cas 1 (cottage de 1998, 165 m², secteur Hull de Gatineau ; DONNÉES SYNTHÉTIQUES,
data/synthetic/cas1_gatineau_*). Le zonage retenu pour l'exemple est FICTIF (zone « H-2f » : 1 à 3 logements).

Quatre usages sont confrontés :
  S0 — maison actuelle conservée (valeur dans l'usage actuel : conclusion du ch. 20, 600 000 $) ;
  U  — démolition et maison neuve de 165 m² (prix de sortie : modèle hédonique du cas 1, âge 0) ;
  D  — démolition et duplex (2 × 120 m², vendus en copropriété divise) ;
  T  — démolition et triplex (3 × 100 m²).
Pour chacun : quatre tests de la PHMU, méthode résiduelle de développement avec calendrier (financement intérimaire
explicite via traite.finance.interets_intermediaires, actualisation du résidu, contribution pour fins de parcs),
sensibilité, puis option réelle de redéveloppement (arbre binomial à deux périodes) et seuil de Capozza–Helsley.

Réutilise : traite.cas1 (hedonique, predire), traite.finance (interets_intermediaires, marge_sur_prix, marge_sur_cout),
ch21_cout_neuf (devis du cas 1 : D + I), ch08_zonage (densite_admissible, contribution_parcs, prime_option,
valeur_residuelle_redeveloppement pour la version statique).

Usage : python3 analysis/python/cases/ch41_phmu.py
"""
from __future__ import annotations
import sys, pathlib, math
import numpy as np
import pandas as pd

ROOT = pathlib.Path(__file__).resolve().parents[3]
sys.path.insert(0, str(ROOT / "analysis/python")); sys.path.insert(0, str(ROOT / "analysis/python/cases"))
from traite import cas1, finance as fin  # noqa: E402
import ch21_cout_neuf as c21              # noqa: E402
import ch08_zonage as c08                 # noqa: E402

# ----------------------------------------------------------------------------- constantes partagées (provenance déclarée)
V_ACT = 600_000.0        # ch. 20 : valeur la plus probable conclue pour la maison du cas 1 (fourchette 575 000 – 640 000 $)
V_L_MARCHE = 277_000.0   # ch. 24/25 : valeur du terrain du sujet retenue par comparaison/extraction (fourchette 240 000 – 285 000 $)
IPCB_MAISON_T3_2026 = 116.5    # Statistique Canada 18-10-0289-01, Montréal, maison individuelle, T3 2026 (SRC-E-017)
IPCB_APPART_T3_2026 = 118.5    # idem, immeuble d'appartements de moins de cinq étages
TAUX_PRET = 0.0619       # taux hypothécaire affiché 5 ans, prêt conventionnel, 28 octobre 2026 (Banque du Canada, série V80691335 ; ch. 23)
PART_PRET = 0.75         # prêt de construction : 75 % des coûts (hypothèse, comme au ch. 23)
EFFICACITE = 0.85        # superficie habitable / aire brute de plancher (hypothèse, comme au ch. 23)

# zone FICTIVE « H-2f » : hypothèses de zonage déclarées (grille calquée sur la colonne H2 de Val-Fictif, ch. 8, plafond ramené à 3 logements)
GRILLE = {"s_min_lot": 450.0, "s_min_par_logement": 150.0, "logements_max": 3, "cos_max": 0.60, "hauteur_max_etages": 3,
          "stationnement_par_logement": 1, "marge_avant": 6.0, "marge_lat": 2.0, "marge_arr": 9.0}

# hypothèses communes du promoteur (déclarées dans le texte ; sensibilité à la fin)
HYP = {"demolition": 30_000.0,        # démolition, désamiantage sommaire, préparation du site (hypothèse, idem ch. 8)
       "frais_vente": 0.04,           # courtage, notaire, mise en marché : 4 % du produit de la vente (hypothèse)
       "marge": 0.15, "marge_base": "prix",   # marge du promoteur-marchand : 15 % du produit de la vente (hypothèse, idem ch. 8)
       "taux_parcs": 0.10,            # contribution pour fins de parcs exigée au redéveloppement : 10 % de la valeur du site (LAU 117.4, plafond)
       "taux_pret": TAUX_PRET, "part_pret": PART_PRET, "r_portage": TAUX_PRET}

# variante « construction sur commande » (acheteur-occupant qui fait bâtir) : la convention du ch. 21 — incitatif de 10 % des coûts de
# construction (directs, indirects, démolition, financement), pas de frais de vente (rien n'est vendu), pas de contribution pour fins de parcs
HYP_OCCUPANT = {**HYP, "frais_vente": 0.0, "marge": 0.10, "marge_base": "cout", "taux_parcs": 0.0}

# programmes de construction (hypothèses) ; facteur_prix = décote déclarée d'une unité superposée/contiguë vendue en copropriété divise
# par rapport à la maison isolée neuve de même superficie que prédit le modèle hédonique (1,00 = maison isolée)
PROGRAMMES = {
    "S0": {"nom": "Maison actuelle conservée", "n": 1, "s_unite": 165.0, "garage": 1, "sous_sol": 1, "facteur_prix": 1.0, "duree_ans": 0.0, "neuf": False},
    "U": {"nom": "Maison neuve de 165 m²", "n": 1, "s_unite": 165.0, "garage": 1, "sous_sol": 1, "facteur_prix": 1.00, "duree_ans": 1.00, "neuf": True},
    "D": {"nom": "Duplex, 2 × 120 m²", "n": 2, "s_unite": 120.0, "garage": 0, "sous_sol": 1, "facteur_prix": 0.90, "duree_ans": 1.25, "neuf": True},
    "T": {"nom": "Triplex, 3 × 100 m²", "n": 3, "s_unite": 100.0, "garage": 0, "sous_sol": 1, "facteur_prix": 0.85, "duree_ans": 1.50, "neuf": True},
    "Q": {"nom": "Quadruplex, 4 × 85 m²", "n": 4, "s_unite": 85.0, "garage": 0, "sous_sol": 1, "facteur_prix": 0.80, "duree_ans": 1.75, "neuf": True},
}


# ----------------------------------------------------------------------------- intrants : coûts (ch. 21) et prix de sortie (cas 1)

def couts_unitaires() -> dict:
    """Coût direct + indirect par m² habitable du devis du cas 1 (ch. 21), et transposition aux plex par le ratio IPCB (ch. 23)."""
    devis, sujet = c21.charger(); cn = c21.cout_neuf_sujet(devis)
    hab = float(sujet["superficie_habitable_m2"])
    di_m2 = (cn["direct"] + cn["indirects"]) / hab
    ratio = IPCB_APPART_T3_2026 / IPCB_MAISON_T3_2026
    return {"direct": cn["direct"], "indirects": cn["indirects"], "incitatif": cn["incitatif"], "cout_neuf": cn["total"],
            "di_m2_maison": di_m2, "ratio_ipcb_plex": ratio, "di_m2_plex": di_m2 * ratio, "cn_m2_maison": cn["total"] / hab, "sujet": sujet}


def modele_cas1() -> dict:
    d = cas1.charger(); m = cas1.hedonique(d["ventes"]); return {"modele": m, "sujet": d["sujet"]}


def prix_unite(m: dict, sujet: dict, s_unite: float, n: int, garage: int, sous_sol: int, neuf: bool = True) -> float:
    """Prix prédit par le modèle hédonique du cas 1 pour une habitation isolée neuve de s_unite m² sur 620/n m² de terrain,
    à Hull, à la date d'évaluation (correction de Duan incluse). La décote d'unité superposée est appliquée ensuite (facteur_prix)."""
    x = dict(sujet); x.update({"superficie_habitable_m2": s_unite, "superficie_terrain_m2": sujet["superficie_terrain_m2"] / n,
                               "garage_places": garage, "sous_sol_fini": sous_sol})
    if neuf:
        x["age"] = 0; x["annee_construction"] = 2026
    return float(cas1.predire(m, cas1.sujet_df(x))[0])


def produit_vente(m: dict, sujet: dict, prog: dict) -> dict:
    p_isolee = prix_unite(m, sujet, prog["s_unite"], prog["n"], prog["garage"], prog["sous_sol"], prog["neuf"])
    p_unite = p_isolee * prog["facteur_prix"]
    return {"prix_isolee_equiv": p_isolee, "prix_unite": p_unite, "gdv": p_unite * prog["n"], "s_total": prog["s_unite"] * prog["n"]}


# ----------------------------------------------------------------------------- les quatre tests

def test_legal(prog: dict, grille: dict = GRILLE, s_terrain: float = 620.0) -> dict:
    g = c08.densite_admissible(s_terrain, grille["s_min_par_logement"], grille["logements_max"], grille["cos_max"], grille["s_min_lot"])
    brut = prog["n"] * prog["s_unite"] / EFFICACITE
    ok_n = prog["n"] <= g["logements_admissibles"]; ok_cos = brut <= g["aire_plancher_max_m2"]
    return {"logements_admissibles": g["logements_admissibles"], "aire_max_m2": g["aire_plancher_max_m2"], "aire_brute_m2": brut,
            "ok_logements": ok_n, "ok_cos": ok_cos, "permis": bool(g["lot_conforme"] and ok_n and ok_cos)}


def test_physique(prog: dict, grille: dict = GRILLE, facade_m: float = 20.0, profondeur_m: float = 31.0) -> dict:
    """Enveloppe constructible après marges (façade et profondeur du lot : hypothèses, le dossier du cas ne les précise pas) ;
    stationnement : 1 case par logement de 2,6 m × 5,5 m en cour avant/latérale (hypothèse)."""
    largeur = facade_m - 2 * grille["marge_lat"]; prof = profondeur_m - grille["marge_avant"] - grille["marge_arr"]
    empreinte_max = largeur * prof
    brut = prog["n"] * prog["s_unite"] / EFFICACITE
    etages = math.ceil(brut / empreinte_max) if empreinte_max > 0 else 99
    ok_etages = etages <= grille["hauteur_max_etages"]
    cases = prog["n"] * grille["stationnement_par_logement"]; ok_stat = cases * 2.6 <= largeur + 2 * grille["marge_lat"]  # cases côte à côte en façade
    return {"empreinte_max_m2": empreinte_max, "etages_requis": etages, "ok_etages": ok_etages, "cases": cases, "ok_stationnement": ok_stat,
            "possible": bool(ok_etages and ok_stat)}


# ----------------------------------------------------------------------------- méthode résiduelle de développement (dynamique)

def residuel_developpement(gdv: float, s_total: float, di_m2: float, duree_ans: float, hyp: dict = HYP, demolir: bool = True) -> dict:
    """Valeur résiduelle du terrain avec calendrier :
        V_L = [GDV(1 − c_v) − C − D − F − m·GDV] / [(1 + τ_parcs)(1 + r)^T]
    où C = coûts directs + indirects (di_m2 × s_total), D = démolition, F = intérêts intérimaires du prêt de construction
    (tirage uniforme, mise de fonds d'abord), m = marge sur le prix, τ_parcs = contribution pour fins de parcs (assise sur
    la valeur du site), r = taux de portage du terrain, T = durée totale (chantier + vente). r = 0 redonne la version
    statique du ch. 8 (contribution incluse) ; r = 0 et τ = 0 redonnent le résidu « de manuel »."""
    C = di_m2 * s_total; D = hyp["demolition"] if demolir else 0.0
    couts_fin = C + D
    pret = hyp["part_pret"] * couts_fin
    F = fin.interets_intermediaires(pret, hyp["taux_pret"], duree_ans, mode="continu", mise_de_fonds=couts_fin - pret, cout_total=couts_fin)["interets"] if duree_ans > 0 else 0.0
    frais = hyp["frais_vente"] * gdv
    marge = hyp["marge"] * gdv if hyp.get("marge_base", "prix") == "prix" else hyp["marge"] * (C + D + F)
    reste = gdv - frais - C - D - F - marge
    facteur_temps = (1.0 + hyp["r_portage"]) ** duree_ans
    v_l = reste / ((1.0 + hyp["taux_parcs"]) * facteur_temps)
    return {"gdv": gdv, "frais_vente": frais, "couts_construction": C, "demolition": D, "pret": pret, "financement": F, "marge": marge,
            "reste": reste, "facteur_temps": facteur_temps, "valeur_residuelle": v_l, "contribution_parcs": v_l * hyp["taux_parcs"],
            "reste_statique": reste, "v_l_statique": reste / (1.0 + hyp["taux_parcs"]), "cout_total_hors_terrain": C + D + F + frais + marge,
            "marge_sur_cout": marge / (C + D + F + frais + v_l) if (C + D + F + frais + v_l) > 0 else float("nan")}


def seuil_gdv(res: dict, v_cible: float, hyp: dict = HYP) -> float:
    """Produit de vente qui rend la valeur résiduelle égale à v_cible (résolution de la forme linéaire en GDV)."""
    a, b = _a_b(res, hyp)
    return (v_cible + b) / a


def _a_b(res: dict, hyp: dict) -> tuple[float, float]:
    """Forme linéaire du résidu : V_L = a·GDV − b. La marge entre dans a si elle est assise sur le prix, dans b si elle est assise sur le coût."""
    den = (1.0 + hyp["taux_parcs"]) * res["facteur_temps"]
    fixes = res["couts_construction"] + res["demolition"] + res["financement"]
    if hyp.get("marge_base", "prix") == "prix":
        return (1.0 - hyp["frais_vente"] - hyp["marge"]) / den, fixes / den
    return (1.0 - hyp["frais_vente"]) / den, fixes * (1.0 + hyp["marge"]) / den


# ----------------------------------------------------------------------------- option réelle : arbre binomial

def parametres_binomiaux(sigma: float, dt: float, r: float, y: float) -> dict:
    """u = e^{σ√Δt}, d = 1/u ; probabilité risque-neutre d'un actif à rendement de détention y :
    q = (F − d)/(u − d), F = [(1 + r)/(1 + y)]^Δt (facteur de croissance à terme)."""
    u = math.exp(sigma * math.sqrt(dt)); d = 1.0 / u
    F = ((1.0 + r) / (1.0 + y)) ** dt
    q = (F - d) / (u - d)
    return {"u": u, "d": d, "F": F, "q": q, "disc": (1.0 + r) ** (-dt)}


def option_redeveloppement(gdv0: float, a: float, b: float, sigma: float, dt: float, r: float, y: float, n: int = 2) -> dict:
    """Option américaine de redéveloppement sur n périodes de Δt années. À chaque nœud, le propriétaire peut vendre à un
    promoteur (ou se faire promoteur) pour V_L(GDV) = a·GDV − b' ; la prime par rapport à la valeur en usage V_act est
    max(a·GDV − b, 0) avec b = b' + V_act. Rend la prime (américaine et européenne), la règle VAN immédiate et la politique
    d'exercice par nœud."""
    p = parametres_binomiaux(sigma, dt, r, y); u, d, q, disc = p["u"], p["d"], p["q"], p["disc"]
    payoff = lambda g: max(a * g - b, 0.0)
    noeuds = {}
    V = [payoff(gdv0 * u ** j * d ** (n - j)) for j in range(n + 1)]
    for j in range(n + 1):
        noeuds[(n, j)] = {"gdv": gdv0 * u ** j * d ** (n - j), "exercice": V[j], "continuation": None, "valeur": V[j]}
    V_eur = list(V)
    for t in range(n - 1, -1, -1):
        Vn, Vn_eur = [], []
        for j in range(t + 1):
            g = gdv0 * u ** j * d ** (t - j)
            cont = disc * (q * V[j + 1] + (1 - q) * V[j]); ex = payoff(g)
            Vn.append(max(ex, cont)); Vn_eur.append(disc * (q * V_eur[j + 1] + (1 - q) * V_eur[j]))
            noeuds[(t, j)] = {"gdv": g, "exercice": ex, "continuation": cont, "valeur": max(ex, cont), "exercer": ex > cont and ex > 0}
        V, V_eur = Vn, Vn_eur
    return {**p, "prime_americaine": V[0], "prime_europeenne": V_eur[0], "van_immediate": payoff(gdv0), "noeuds": noeuds,
            "gdv_seuil_van": b / a, "a": a, "b": b, "n": n, "dt": dt, "sigma": sigma, "r": r, "y": y}


def coefficients_option(res: dict, v_act: float, hyp: dict = HYP) -> tuple[float, float]:
    """V_L(GDV) = a·GDV − b' (forme linéaire du résidu) ; prime d'exercice = a·GDV − (b' + V_act)."""
    a, b_prime = _a_b(res, hyp)
    return a, b_prime + v_act


def frontiere_exercice_t1(a: float, b: float, sigma: float, dt: float, r: float, y: float) -> float:
    """Plus petit GDV au nœud t = 1 (arbre à deux périodes) où l'exercice immédiat vaut au moins la continuation (bissection)."""
    p = parametres_binomiaux(sigma, dt, r, y); u, d, q, disc = p["u"], p["d"], p["q"], p["disc"]
    f = lambda g: (a * g - b) - disc * (q * max(a * g * u - b, 0.0) + (1 - q) * max(a * g * d - b, 0.0))
    lo, hi = b / a, 50 * b / a
    if f(hi) < 0:
        return float("inf")
    for _ in range(200):
        mid = 0.5 * (lo + hi)
        if f(mid) >= 0:
            hi = mid
        else:
            lo = mid
    return hi


# ----------------------------------------------------------------------------- Capozza–Helsley (version simplifiée)

def seuil_capozza_helsley(R_A: float, C: float, r: float, g: float, sigma: float) -> dict:
    """Rente de réserve R* déclenchant la conversion lorsque la rente du nouvel usage suit dR = g dt + σ dz :
        R* = R_A + rC − g/r + 1/α,  α = [−g + sqrt(g² + 2 r σ²)] / σ²   (σ → 0 : R* → R_A + rC).
    Prime d'irréversibilité = 1/α − g/r (nulle sans incertitude)."""
    if sigma <= 0:
        alpha = r / g if g > 0 else float("inf"); prime = 0.0
    else:
        alpha = 2.0 * r / (g + math.sqrt(g * g + 2.0 * r * sigma * sigma))   # forme rationalisée de [−g + √(g² + 2rσ²)]/σ², stable quand σ → 0
        prime = 1.0 / alpha - g / r
    R_cert = R_A + r * C
    return {"alpha": alpha, "R_certitude": R_cert, "prime_irreversibilite": prime, "R_seuil": R_cert + prime}


# ----------------------------------------------------------------------------- exemple intégral

OPTION = {"sigma": 0.10, "dt": 5.0, "r": 0.02, "y": 0.045, "n": 2}   # volatilité réelle annuelle du prix des unités, pas de 5 ans, taux réel, rendement locatif net


def prix_logement_plex_existant() -> dict:
    """Prix par logement des triplex existants du cas 2 (Québec, Limoilou ; traite.cas2) — contrepoint « vente en bloc »."""
    from traite import cas2
    d = cas2.charger(); c = cas2.comparables_plex(d["ventes"], d["sujet"]); ind = cas2.indicateurs_marche(c)
    return {"prix_par_logement_triplex": ind["prix_par_logement_triplex"], "mrb_triplex": ind["mrb_median_triplex"], "n": ind["n"]}


def exemple_complet() -> dict:
    cu = couts_unitaires(); mc = modele_cas1(); m, sujet = mc["modele"], mc["sujet"]
    rows, res, tests = [], {}, {}
    for k, prog in PROGRAMMES.items():
        pv = produit_vente(m, sujet, prog) if prog["neuf"] else {"prix_isolee_equiv": float(cas1.predire(m, cas1.sujet_df(sujet))[0]), "prix_unite": V_ACT, "gdv": V_ACT, "s_total": 165.0}
        tl = test_legal(prog); tp = test_physique(prog)
        di = cu["di_m2_maison"] if prog["n"] == 1 else cu["di_m2_plex"]
        if prog["neuf"]:
            r_prom = residuel_developpement(pv["gdv"], pv["s_total"], di, prog["duree_ans"], HYP)
            r_occ = residuel_developpement(pv["gdv"], pv["s_total"], di, prog["duree_ans"], HYP_OCCUPANT)
        else:
            r_prom = r_occ = None
        res[k] = {"programme": prog, "vente": pv, "legal": tl, "physique": tp, "promoteur": r_prom, "occupant": r_occ, "di_m2": di}
        rows.append({"usage": k, "nom": prog["nom"], "n": prog["n"], "s_total": pv["s_total"], "prix_unite": pv["prix_unite"], "gdv": pv["gdv"],
                     "permis": tl["permis"], "possible": tp["possible"],
                     "v_l_promoteur": r_prom["valeur_residuelle"] if r_prom else float("nan"),
                     "v_l_occupant": r_occ["valeur_residuelle"] if r_occ else float("nan")})
    tableau = pd.DataFrame(rows)
    # contrepoint : triplex vendu en bloc à un investisseur au prix par logement des triplex existants du cas 2
    bloc = prix_logement_plex_existant(); uT = res["T"]
    gdv_bloc = bloc["prix_par_logement_triplex"] * 3
    res["T_bloc"] = {"gdv": gdv_bloc, "vente": {"gdv": gdv_bloc}, "promoteur": residuel_developpement(gdv_bloc, uT["vente"]["s_total"], uT["di_m2"], uT["programme"]["duree_ans"], HYP), **bloc}
    # PHMU du terrain vacant : usages permis et possibles, résidu positif ; maximalement productif = plus grande valeur résiduelle
    # (marge du promoteur pour les plex ; marge du constructeur sur commande pour la maison, usage dont l'acheteur type est un occupant)
    cand = tableau[(tableau.usage != "S0") & tableau.permis & tableau.possible].copy()
    cand["v_l_retenue"] = np.where(cand.n == 1, cand.v_l_occupant, cand.v_l_promoteur)
    cand["faisable"] = cand["v_l_retenue"] > 0
    cand_f = cand[cand.faisable]
    phmu_vacant = cand_f.loc[cand_f.v_l_retenue.idxmax(), "usage"]
    phmu_vacant_promoteur = cand_f.loc[cand_f.v_l_promoteur.idxmax(), "usage"]
    # PHMU du bien amélioré : max(V_act, V_L − D) — la démolition est déjà incluse dans V_L
    meilleure_redev = cand["v_l_promoteur"].max()
    primes = {u: c08.prime_option(V_ACT, v) for u, v in zip(cand.usage, cand.v_l_promoteur)}
    phmu_ameliore = "S0" if V_ACT >= meilleure_redev else cand.loc[cand.v_l_promoteur.idxmax(), "usage"]
    contribution_batiment = V_ACT - V_L_MARCHE
    # seuils de bascule (triplex, convention promoteur)
    rT = res["T"]["promoteur"]; seuil_T = seuil_gdv(rT, V_ACT); hausse_T = seuil_T / rT["gdv"] - 1
    seuil_T_marche = seuil_gdv(rT, V_L_MARCHE); hausse_T_marche = seuil_T_marche / rT["gdv"] - 1
    # option réelle (triplex) : arbre binomial à deux périodes de 5 ans
    a, b = coefficients_option(rT, V_ACT); O = OPTION
    opt = option_redeveloppement(rT["gdv"], a, b, O["sigma"], O["dt"], O["r"], O["y"], O["n"])
    front = frontiere_exercice_t1(a, b, O["sigma"], O["dt"], O["r"], O["y"])
    # sensibilité de l'option : σ et V_act
    sens_opt = {}
    for v_act in (600_000.0, 450_000.0, 300_000.0):
        a_, b_ = coefficients_option(rT, v_act)
        sens_opt[v_act] = {s: option_redeveloppement(rT["gdv"], a_, b_, s, O["dt"], O["r"], O["y"], O["n"])["prime_americaine"] for s in (0.05, 0.10, 0.15, 0.20, 0.25)}
    # Capozza–Helsley : rente implicite de la maison R_A = y·V_act ; rente du triplex R = y·GDV ; coût de conversion C = coûts + démolition + financement
    ch = {"R_A": O["y"] * V_ACT, "R_T": O["y"] * rT["gdv"], "C": rT["couts_construction"] + rT["demolition"] + rT["financement"]}
    ch["g"] = 0.01 * ch["R_T"]; ch["sigma"] = 0.10 * ch["R_T"]
    ch["resultat"] = seuil_capozza_helsley(ch["R_A"], ch["C"], O["y"], ch["g"], ch["sigma"])
    ch["hausse_requise_certitude"] = ch["resultat"]["R_certitude"] / ch["R_T"] - 1; ch["hausse_requise_incertitude"] = ch["resultat"]["R_seuil"] / ch["R_T"] - 1
    # sensibilité de la valeur résiduelle du triplex (convention promoteur)
    sens = sensibilite_triplex(res["T"], cu)
    return {"couts": cu, "modele": m, "sujet": sujet, "usages": res, "tableau": tableau, "candidats": cand, "phmu_vacant": phmu_vacant,
            "phmu_vacant_promoteur": phmu_vacant_promoteur, "phmu_ameliore": phmu_ameliore, "primes": primes, "contribution_batiment": contribution_batiment,
            "seuil_gdv_T": seuil_T, "hausse_requise_T": hausse_T, "seuil_gdv_T_marche": seuil_T_marche, "hausse_requise_T_marche": hausse_T_marche,
            "option": opt, "a": a, "b": b, "frontiere_t1": front, "sens_option": sens_opt, "capozza": ch, "sensibilite": sens, "v_act": V_ACT, "v_l_marche": V_L_MARCHE}


def sensibilite_triplex(uT: dict, cu: dict) -> pd.DataFrame:
    """Variations une à une autour du cas de base (triplex, convention promoteur)."""
    prog = uT["programme"]; pv = uT["vente"]; base = uT["promoteur"]["valeur_residuelle"]
    def run(gdv=pv["gdv"], di=cu["di_m2_plex"], duree=prog["duree_ans"], **kw):
        h = {**HYP, **kw}; return residuel_developpement(gdv, pv["s_total"], di, duree, h)["valeur_residuelle"]
    lignes = [
        ("Prix des unités", "−15 %", run(gdv=pv["gdv"] * 0.85), "+15 %", run(gdv=pv["gdv"] * 1.15)),
        ("Coût de construction", "−10 %", run(di=cu["di_m2_plex"] * 0.90), "+10 %", run(di=cu["di_m2_plex"] * 1.10)),
        ("Marge du promoteur", "10 % du prix", run(marge=0.10), "20 % du prix", run(marge=0.20)),
        ("Durée totale", "12 mois", run(duree=1.0), "24 mois", run(duree=2.0)),
        ("Taux (prêt et portage)", "5,00 %", run(taux_pret=0.05, r_portage=0.05), "7,50 %", run(taux_pret=0.075, r_portage=0.075)),
        ("Contribution pour fins de parcs", "0 %", run(taux_parcs=0.0), "10 %", base),
        ("Frais de vente", "2 %", run(frais_vente=0.02), "6 %", run(frais_vente=0.06)),
    ]
    return pd.DataFrame(lignes, columns=["parametre", "bas_etiquette", "bas", "haut_etiquette", "haut"]).assign(base=base)


def courbe_densite(res: dict, cu: dict, facteurs_prix=(0.85, 1.0, 1.15)) -> pd.DataFrame:
    """Valeur résiduelle (convention promoteur) selon le nombre de logements, pour plusieurs niveaux de prix."""
    rows = []
    for k in ("U", "D", "T", "Q"):
        u = res["usages"][k]; prog = u["programme"]
        for f in facteurs_prix:
            r = residuel_developpement(u["vente"]["gdv"] * f, u["vente"]["s_total"], u["di_m2"], prog["duree_ans"], HYP)
            rows.append({"usage": k, "n": prog["n"], "facteur_prix": f, "v_l": r["valeur_residuelle"], "permis": u["legal"]["permis"]})
    return pd.DataFrame(rows)


def surface_prix_cout(uT: dict, cu: dict, fp=np.linspace(0.80, 1.40, 61), fc=np.linspace(0.80, 1.20, 41)) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    prog = uT["programme"]; pv = uT["vente"]
    Z = np.array([[residuel_developpement(pv["gdv"] * p, pv["s_total"], cu["di_m2_plex"] * c, prog["duree_ans"], HYP)["valeur_residuelle"] for p in fp] for c in fc])
    return fp, fc, Z


def _d(x): return f"{x:,.0f} $".replace(",", " ")


def main() -> dict:
    R = exemple_complet(); cu = R["couts"]
    print("=== Ch. 41 — PHMU du terrain du cas 1 (620 m², Hull ; DONNÉES SYNTHÉTIQUES ; zonage FICTIF H-2f) ===")
    print(f"Coûts ch. 21 : D = {_d(cu['direct'])} ; I = {_d(cu['indirects'])} ; D+I/m² hab. = {cu['di_m2_maison']:.2f} $/m² ; ratio IPCB plex {cu['ratio_ipcb_plex']:.4f} → {cu['di_m2_plex']:.2f} $/m²")
    print(f"V_act (ch. 20) = {_d(R['v_act'])} ; V_L marché (ch. 24/25) = {_d(R['v_l_marche'])} ; contribution du bâtiment = {_d(R['contribution_batiment'])}")
    print("\n--- Programmes, tests et produits de vente ---")
    for k, u in R["usages"].items():
        if "legal" not in u:
            continue
        tl, tp, pv = u["legal"], u["physique"], u["vente"]
        print(f"{k:>2} {u['programme']['nom']:<28} n={u['programme']['n']} S={pv['s_total']:.0f} m² brut={tl['aire_brute_m2']:.0f} (max {tl['aire_max_m2']:.0f}) "
              f"permis={tl['permis']} possible={tp['possible']} (étages {tp['etages_requis']}, cases {tp['cases']}) | prix isolée équiv. {_d(pv['prix_isolee_equiv'])} → unité {_d(pv['prix_unite'])} ; GDV {_d(pv['gdv'])}")
    print("\n--- Résidu de développement (convention promoteur : marge 15 % du prix, frais de vente 4 %, parcs 10 %, prêt 75 % à 6,19 %, portage 6,19 %) ---")
    for k in ("U", "D", "T", "Q"):
        r = R["usages"][k]["promoteur"]
        print(f"{k}: GDV {_d(r['gdv'])} − frais {_d(r['frais_vente'])} − C {_d(r['couts_construction'])} − démol. {_d(r['demolition'])} − F {_d(r['financement'])} − marge {_d(r['marge'])} = reste {_d(r['reste'])} ; "
              f"/(1,10 × {r['facteur_temps']:.4f}) → V_L = {_d(r['valeur_residuelle'])} (parcs {_d(r['contribution_parcs'])}) ; statique {_d(r['v_l_statique'])}")
    print("--- Convention constructeur sur commande (incitatif 10 % des coûts de construction, pas de frais de vente ni de parcs) ---")
    for k in ("U", "D", "T"):
        r = R["usages"][k]["occupant"]; print(f"{k}: V_L = {_d(r['valeur_residuelle'])} (marge {_d(r['marge'])}, reste {_d(r['reste'])}, F {_d(r['financement'])})")
    tb = R["usages"]["T_bloc"]; rb = tb["promoteur"]
    print(f"--- Triplex vendu en bloc (cas 2 : {_d(tb['prix_par_logement_triplex'])} par logement, n = {tb['n']} ventes) : GDV {_d(tb['gdv'])} → reste {_d(rb['reste'])} → V_L = {_d(rb['valeur_residuelle'])}")
    print(f"\nPHMU du terrain vacant : {R['phmu_vacant']} (par résidu promoteur seul : {R['phmu_vacant_promoteur']}) ; PHMU du bien amélioré : {R['phmu_ameliore']}")
    print("Primes d'option immédiates max(0, V_L − V_act) : " + ", ".join(f"{u} {_d(p)}" for u, p in R["primes"].items()))
    print(f"Seuil GDV triplex pour égaler V_act : {_d(R['seuil_gdv_T'])} (+{R['hausse_requise_T']*100:.1f} %) ; pour égaler V_L marché : {_d(R['seuil_gdv_T_marche'])} ({R['hausse_requise_T_marche']*100:+.1f} %)")
    o = R["option"]
    print(f"\n--- Option réelle (triplex) : σ = {o['sigma']:.2f}, Δt = {o['dt']:.0f} ans, r = {o['r']:.3f}, y = {o['y']:.3f} ; u = {o['u']:.4f}, d = {o['d']:.4f}, F = {o['F']:.4f}, q = {o['q']:.4f}, disc = {o['disc']:.4f}")
    print(f"a = {R['a']:.4f} ; b = {_d(R['b'])} ; GDV* (VAN nulle) = {_d(o['gdv_seuil_van'])} ; frontière d'exercice à t = 1 : {_d(R['frontiere_t1'])}")
    for (t, j), nd in sorted(o["noeuds"].items()):
        print(f"  nœud t={t} j={j}: GDV {_d(nd['gdv'])} ; exercice {_d(nd['exercice'])} ; continuation {('—' if nd['continuation'] is None else _d(nd['continuation']))} ; valeur {_d(nd['valeur'])}")
    print(f"Prime américaine {_d(o['prime_americaine'])} ({o['prime_americaine']/R['v_act']*100:.2f} % de V_act) ; européenne {_d(o['prime_europeenne'])} ; VAN immédiate {_d(o['van_immediate'])}")
    print("Sensibilité de la prime (lignes V_act, colonnes σ = 5/10/15/20/25 %) :")
    for v, d in R["sens_option"].items():
        print(f"  V_act {_d(v)} : " + "  ".join(_d(x) for x in d.values()))
    ch = R["capozza"]; cr = ch["resultat"]
    print(f"\n--- Capozza–Helsley : R_A {_d(ch['R_A'])}/an ; R_T {_d(ch['R_T'])}/an ; C {_d(ch['C'])} ; g {_d(ch['g'])}/an ; σ {_d(ch['sigma'])}/an ; r = y = {OPTION['y']}")
    print(f"α = {cr['alpha']:.6e} ; R* certitude {_d(cr['R_certitude'])} (+{ch['hausse_requise_certitude']*100:.1f} %) ; prime d'irréversibilité {_d(cr['prime_irreversibilite'])} ; R* {_d(cr['R_seuil'])} (+{ch['hausse_requise_incertitude']*100:.1f} %)")
    print("\n--- Sensibilité V_L triplex (promoteur) ---")
    print(R["sensibilite"].assign(bas=lambda d: d.bas.round(0), haut=lambda d: d.haut.round(0), base=lambda d: d.base.round(0)).to_string(index=False))
    return R


if __name__ == "__main__":
    main()