import sys, csv, math, collections
import numpy as np
import pandas as pd
import pycountry
import geonamescache
from scipy.spatial import cKDTree
from scipy.optimize import linprog
from scipy.sparse import lil_matrix
from global_land_mask import globe
D = sys.argv[1] if len(sys.argv) > 1 else "."
SOUPLE = len(sys.argv) > 2 and sys.argv[2] == "souple"
ECRIRE = len(sys.argv) <= 2 or SOUPLE
NOM = "plan_mondial_artesia_5_v1_recalcul" + ("_souple" if SOUPLE else "")
MAX_KM, MAX_MER, DIST_MAX = 5000.0, 1500.0, 400.0
SOLEIL = [(1800, 871.3), (1900, 917.4), (2000, 959.8), (2100, 998.5), (2200, 1033.7), (2300, 1065.7), (2400, 1094.8), (2500, 1121.5), (2600, 1145.6), (2700, 1167.1)]
MIXTE = [(1800, 940.1), (2100, 1002.3), (2400, 1054.8), (2700, 1096.3)]
VENT = [(5.5, 662.8), (5.75, 715.1), (6.0, 766.0), (6.25, 815.2), (6.5, 862.7), (6.75, 908.1), (7.0, 951.0), (7.5, 1030.3), (8.0, 1099.1), (8.5, 1155.5), (9.0, 1201.0)]
CALAGE_V14 = 885.0 / 899.4
COUT_SOLEIL, COUT_MIXTE, COUT_CRETE, EAU_APPOINT = 465.0, 461.0, 364.0, 15.0
PACHA, GOLD, FRAC_GOLD, PART_V21, PRUDENCE_VENT = (155.6, 31.1), (56.0, 6.4), 0.47, 2.0 / 3.0, 0.9
FIXE = {"Titan Claudia": (1500, 43.0, 49.0), "Pacha existantes": (2700, 420.0, 84.0), "Gold existantes": (1126, 63.0, 7.2), "Sky existantes": (356, 0.0, 45.1), "O Sky sur les Sky": (583, 56.0, 29.3)}
HYBRIDE = 810.0
TOTAL_PLAN = 17360.0
FLOTTANT = (21.6, 85.2, 1100.0, 73.1)
LF = 0.5
GW = 1.0 / (8.76 * LF)
AC, POSTE, POSTE_EX, ARRIVEE, CONV, TERRE, MER, DETOUR, SEUIL_DOM, RAYON = 0.60, 50.0, 20.0, 100.0, 200.0, 0.30, 1.60, 1.15, 300.0, 50.0
DEM = np.array([0.85] * 6 + [1.0] * 4 + ([1.5] * 6 if SOUPLE else [1.0] * 6) + [1.0] + [1.05 if SOUPLE else 1.25] * 6 + [0.95])
OMEGA = float(DEM[[h for h in range(24) if h >= 17 or h < 5]].sum() / DEM.sum())
REND, HYDRO_DEP, EXT_MAX = 0.78, 0.3, 2.0
def hyd(k):
    return 19.1 * k ** 0.8 + 13.3 * k ** 0.75 + 67.8 * k ** 0.9
EXT_MEUR_GWH = [(hyd(3 * m) - hyd(m)) / (2 * m * 0.588) for m in (1, 2, 4, 6)]
NUIT = {"soleil": 0.18, "poste": 0.18, "mixte": 0.2 + 0.8 * 0.57 * 0.5, "crete": 0.19 + 0.81 * 0.5, "v21": 0.18, "v21v": 0.19 + 0.81 * 0.5, "flottant": 0.1, "gold4": 0.5}
STOCKE = {"soleil": 0.18, "poste": 0.18, "mixte": 0.2, "crete": 0.19}
MAXR = {"soleil": (1.0, 0.955, 0.918, 0.897), "mixte": (1.0, 0.957, 0.920, 0.900), "poste": (1.0, 0.952, 0.910, 0.889), "crete": (1.0, 0.968, 0.940, 0.929)}
TAILLES = ((6, 3), (4, 2), (2, 1))
PRIX_TYPE = {"soleil": "v14_soleil", "poste": "v14_poste", "mixte": "v14_mixte", "crete": "v14_crete", "v21": "v21_vent", "v21v": "v21_vent", "flottant": "hybrides_iles", "gold4": "hydraulique", "v15": "v15"}
def interp_bas(t, h):
    if h < t[0][0]:
        return t[0][1] - (t[1][1] - t[0][1]) * (t[0][0] - h) / (t[1][0] - t[0][0])
    return interp(t, h)
def interp(t, h):
    h = min(max(h, t[0][0]), t[-1][0])
    for (a, x), (b, y) in zip(t, t[1:]):
        if a <= h <= b:
            return x + (y - x) * (h - a) / (b - a)
    return t[-1][1]
def iso3(c):
    if c == "XK":
        return "OWID_KOS"
    if len(c) == 3:
        return c
    p = pycountry.countries.get(alpha_2=c)
    return p.alpha_3 if p else None
def xyz(lat, lon):
    la, lo = np.radians(lat), np.radians(lon)
    return np.c_[np.cos(la) * np.cos(lo), np.cos(la) * np.sin(lo), np.sin(la)]
def km(a, b, c, d):
    a, b, c, d = map(math.radians, (a, b, c, d))
    return 6371.0 * 2 * math.asin(math.sqrt(min(1.0, math.sin((c - a) / 2) ** 2 + math.cos(a) * math.cos(c) * math.sin((d - b) / 2) ** 2)))
def route(a, b, c, d):
    L = km(a, b, c, d)
    n = max(int(L / 25), 2)
    p = np.linspace(0, 1, n)
    v1, v2 = xyz(np.array([a]), np.array([b]))[0], xyz(np.array([c]), np.array([d]))[0]
    om = math.acos(max(-1.0, min(1.0, float(v1 @ v2))))
    if om < 1e-9:
        return 0.0, 0.0
    q = (np.sin((1 - p) * om)[:, None] * v1 + np.sin(p * om)[:, None] * v2) / math.sin(om)
    la, lo = np.degrees(np.arcsin(np.clip(q[:, 2], -1, 1))), np.degrees(np.arctan2(q[:, 1], q[:, 0]))
    mer = L * (1 - globe.is_land(np.clip(la, -89.99, 89.99), lo).mean())
    return (L - mer) * DETOUR, mer
def centre(lats, lons, w):
    v = (xyz(np.asarray(lats, float), np.asarray(lons, float)) * np.asarray(w, float)[:, None]).sum(0)
    return math.degrees(math.atan2(v[2], math.hypot(v[0], v[1]))), math.degrees(math.atan2(v[1], v[0]))
def groupes(n):
    out, r = [], int(math.floor(n + 1e-9))
    for t, k in TAILLES:
        while r >= t:
            out.append((t, k))
            r -= t
    return out, n - sum(t for t, k in out)
pays = pd.read_csv(f"{D}/electricite_pays_prix_part_fossile.csv").set_index("iso")
F = pays.fossile_twh
besoin = F * TOTAL_PLAN / F.sum()
nonloc_twh = sum(v[1] for v in FIXE.values())
nonloc = nonloc_twh * F / F.sum()
PC = pd.read_csv(f"{D}/prix_centrales_pays.csv", keep_default_na=False).set_index("iso")
PX = pd.read_csv(f"{D}/prix_lignes_electriques_pays.csv", keep_default_na=False).set_index("iso")
FA = PX.facteur.astype(float).to_dict()
FD = (PX.continu_terre_meur_km_gw.astype(float) / TERRE).to_dict()
RY = PX.rayon_km.astype(float).to_dict()
def fp(k, t):
    return float(PC.loc[k, PRIX_TYPE[t]]) if k in PC.index else 1.0
def fh(k):
    return float(PC.loc[k, "hydraulique"]) if k in PC.index else 1.0
H = pd.read_csv(f"{D}/owid-energy-data.csv", usecols=["iso_code", "year", "hydro_electricity"])
H = H[H.iso_code.notna()].dropna(subset=["hydro_electricity"]).sort_values("year").groupby("iso_code").tail(1).set_index("iso_code").hydro_electricity
reb = {(round(float(r["lat"]), 3), round(float(r["lon"]), 3)): float(r["rebord_km_6ms"]) for r in csv.DictReader(open(f"{D}/plateaux_rebord_vent.csv", encoding="utf-8"))}
def pacha(p):
    return (PACHA[0] + FRAC_GOLD * GOLD[0], PACHA[1] + FRAC_GOLD * GOLD[1]) if p >= 1500 else (0.0, 0.0)
C = []
def ajoute(**d):
    d.setdefault("gp", 0.0)
    d.setdefault("kp", 0.0)
    C.append(d)
for r in csv.DictReader(open(f"{D}/plateaux_ceinture_soleil.csv", encoding="utf-8")):
    if r["classe"] == "C" or r["occupation"] == "occupe" or int(r["dist_grande_ville_km"]) > DIST_MAX:
        continue
    n = int(r["sites"]) if r["occupation"] == "libre" else int(r["sites"]) // 2
    if n == 0:
        continue
    h, p = int(r["heures_tracker"]), int(r["pluie_mm"])
    kmr = reb.get((round(float(r["lat"]), 3), round(float(r["lon"]), 3)), 0.0)
    nmix = min(n, int(kmr // 12) if p >= 200 else 0)
    gp, kp = pacha(p)
    gs, ks = interp(SOLEIL, h), COUT_SOLEIL + (EAU_APPOINT if p < 75 else 0.0)
    gm = interp(MIXTE, h)
    if nmix and COUT_MIXTE / gm < ks / gs:
        ajoute(type="mixte", n=nmix, g=gm, gp=gp, k=COUT_MIXTE, kp=kp, pacha=p >= 1500, lat=float(r["lat"]), lon=float(r["lon"]), pays=r["pays"], ville=r["ville"], dv=float(r["dist_grande_ville_km"]))
        n -= nmix
    if n:
        ajoute(type="soleil", n=n, g=gs, gp=gp, k=ks, kp=kp, pacha=p >= 1500, lat=float(r["lat"]), lon=float(r["lon"]), pays=r["pays"], ville=r["ville"], dv=float(r["dist_grande_ville_km"]))
for r in csv.DictReader(open(f"{D}/cretes_vent_v14.csv", encoding="utf-8")):
    if int(r["empl_AB"]) == 0 or r["occupation"] == "occupe" or int(r["dist_grande_ville_km"]) > DIST_MAX or r["vent_150m_moyen"] == "":
        continue
    n = int(r["empl_AB"]) if r["occupation"] == "libre" else int(r["empl_AB"]) // 2
    if n == 0:
        continue
    v, p = float(r["vent_150m_moyen"]) * PRUDENCE_VENT, int(r["pluie_mm"])
    gp, kp = pacha(p)
    ajoute(type="crete", n=n, g=interp(VENT, v) * CALAGE_V14, gp=gp, k=COUT_CRETE, kp=kp, pacha=p >= 1500, lat=float(r["lat"]), lon=float(r["lon"]), pays=r["pays"], ville=r["ville"], dv=float(r["dist_grande_ville_km"]))
for r in csv.DictReader(open(f"{D}/v21_retenues_1mw.csv", encoding="utf-8")):
    if r["type"] != "V21 Soleil" or float(r["gwh_v21"]) <= 0:
        continue
    ajoute(type="v21", n=PART_V21, g=float(r["gwh_v21"]), k=float(r["meur_v21"]), pacha=False, lat=float(r["lat"]), lon=float(r["lon"]), pays=r["pays"], ville=r["nom"], dv=0.0)
for r in csv.DictReader(open(f"{D}/vent_barrages_temperes_v21.csv", encoding="utf-8")):
    if float(r["gwh"]) <= 0:
        continue
    ajoute(type="v21v", n=PART_V21, g=float(r["gwh"]), k=float(r["meur"]), pacha=False, lat=float(r["lat"]), lon=float(r["lon"]), pays=r["pays"], ville=r["nom"], dv=0.0)
for r in csv.DictReader(open(f"{D}/postes_fossiles_relief_v14.csv", encoding="utf-8")):
    n = int(r["v14_possibles"])
    if n == 0 or r["heures_tracker"] == "":
        continue
    gp, kp = pacha(int(r["pluie_mm"]))
    ajoute(type="poste", n=n, g=interp_bas(SOLEIL, float(r["heures_tracker"])), gp=gp, k=COUT_SOLEIL + 10.0, kp=kp, pacha=int(r["pluie_mm"]) >= 1500, lat=float(r["latitude"]), lon=float(r["longitude"]), pays=r["country"], ville=r["name"], dv=float(r["dist_relief_km"]))
for r in csv.DictReader(open(f"{D}/retenues_soleil_flottant.csv", encoding="utf-8")):
    if r["pvout_kwh_kwp"] == "" or float(r["poste_mw"]) < 20:
        continue
    n = max(1, min(5, int(float(r["poste_mw"]) // 200)))
    ajoute(type="flottant", n=n, g=FLOTTANT[0] + FLOTTANT[1] * float(r["pvout_kwh_kwp"]) / FLOTTANT[2], k=FLOTTANT[3], pacha=False, lat=float(r["lat"]), lon=float(r["lon"]), pays=r["pays"], ville=r["nom"], dv=5.0)
c = pd.DataFrame(C)
c["iso"] = c.pays.map(iso3)
c = c[c.iso.isin(pays.index)].reset_index(drop=True)
c["gtot"] = c.g + c.gp
c["twh"] = c.gtot / 1000.0
c["gw"] = c.twh * GW
c["fprix"] = [fp(k, t) for k, t in zip(c.iso, c.type)]
c["kloc"] = c.k * c.fprix + c.kp * [fh(k) for k in c.iso]
c["nuit"] = [(NUIT[t] * g + 0.5 * gp) / (g + gp) for t, g, gp in zip(c.type, c.g, c.gp)]
c["stocke"] = [STOCKE.get(t, 0.0) * g / (g + gp) for t, g, gp in zip(c.type, c.g, c.gp)]
g = pd.read_csv(f"{D}/global_power_plant_database.csv", low_memory=False)
g["fossile"] = g.primary_fuel.isin({"Coal", "Gas", "Oil", "Petcoke"})
g = g[(g.fossile & (g.capacity_mw >= 50)) | (~g.fossile & (g.capacity_mw >= 100))].reset_index(drop=True)
libre = np.where(g.fossile, 0.95 * g.capacity_mw / 1000.0, 0.0)
arbre = cKDTree(xyz(g.latitude.values, g.longitude.values))
P = xyz(c.lat.values, c.lon.values)
vois = arbre.query_ball_point(P, 2 * math.sin(RAYON / 6371.0 / 2))
GX = xyz(g.latitude.values, g.longitude.values)
sorte, dpo = ["ville"] * len(c), [0.0] * len(c)
for i in np.argsort((c.kloc / c.gtot).values):
    t = c.type[i]
    if t in ("v21", "v21v"):
        sorte[i] = "barrage"
        continue
    if t == "poste":
        sorte[i] = "poste_v14"
        continue
    if t == "flottant":
        sorte[i] = "barrage_flottant"
        continue
    vs = vois[i]
    if not vs:
        continue
    ds = 6371.0 * 2 * np.arcsin(np.clip(np.linalg.norm(GX[vs] - P[i], axis=1) / 2, 0, 1))
    besoin_gw = c.gw[i] * c.n[i]
    pris = None
    for k in np.argsort(ds):
        j = vs[k]
        if g.fossile[j] and libre[j] >= besoin_gw:
            libre[j] -= besoin_gw
            pris = ("fossile", ds[k])
            break
    if pris is None:
        pris = ("autre", float(ds.min()))
    sorte[i], dpo[i] = pris[0], pris[1] * DETOUR
c["sorte"], c["dpo"] = sorte, dpo
villes = collections.defaultdict(list)
for v in geonamescache.GeonamesCache(min_city_population=15000).get_cities().values():
    k = iso3(v["countrycode"])
    if k:
        villes[k].append((v["latitude"], v["longitude"], v["population"]))
gros = {k: [(a, b) for a, b, p in v if p >= 1000000] or [(a, b) for a, b, p in sorted(v, key=lambda t: -t[2])[:1]] for k, v in villes.items()}
charge = {k: centre([a for a, b, p in v], [b for a, b, p in v], [p for a, b, p in v]) for k, v in villes.items()}
c["dgros"] = [min((km(r.lat, r.lon, a, b) for a, b in gros.get(r.iso, [])), default=0.0) for r in c.itertuples()]
def raccord(r):
    if r.sorte == "barrage":
        return r.gw * ARRIVEE
    if r.type == "poste":
        return r.gw * (r.dv * DETOUR * AC + POSTE_EX)
    if r.type == "flottant":
        return r.gw * (r.dv * AC + POSTE_EX + ARRIVEE)
    if r.sorte == "fossile":
        return r.gw * (r.dpo * AC + POSTE_EX)
    d = r.dpo if r.sorte == "autre" else r.dv
    base = r.gw * (d * AC + (POSTE_EX if r.sorte == "autre" else POSTE)) if d <= SEUIL_DOM else r.gw * (CONV + d * DETOUR * TERRE + POSTE)
    dom = r.gw * (CONV + r.dgros * DETOUR * TERRE) if r.dgros > SEUIL_DOM else 0.0
    return base + dom + r.gw * ARRIVEE
c["racc"] = c.apply(raccord, axis=1) * c.iso.map(lambda k: FA.get(k, 1.0))
V = []
for r in csv.DictReader(open(f"{D}/remparts_v15.csv", encoding="utf-8")):
    k = iso3(r["pays"])
    if k not in pays.index:
        continue
    o = float(r["gwh_pointe"]) / 1000.0
    V.append(dict(iso=k, lat=float(r["lat"]), lon=float(r["lon"]), ville=r["lieu"], rempart=int(r["rempart"]), out=o, draw=o / REND - 0.172 * o, k=float(r["meur"]) * fp(k, "v15"), racc=0.3 * o / 0.986 * ARRIVEE * FA.get(k, 1.0), stock=float(r["stockage_mwh"])))
v = pd.DataFrame(V)
for r in csv.DictReader(open(f"{D}/gold4_sources_artesia_5.csv", encoding="utf-8")):
    if float(r["gwh"]) < 5:
        continue
    k = iso3(r["pays"])
    if k in pays.index:
        c.loc[len(c)] = dict(type="gold4", n=1.0, g=float(r["gwh"]), gp=0.0, k=float(r["meur"]), kp=0.0, pacha=False, lat=0.0, lon=0.0, pays=r["pays"], ville=r["source"], dv=0.0, iso=k, gtot=float(r["gwh"]), twh=float(r["gwh"]) / 1000, gw=float(r["gwh"]) / 1000 * GW, fprix=fh(k), kloc=float(r["meur"]) * fh(k), nuit=0.5, stocke=0.0, sorte="source", dpo=0.0, dgros=0.0, racc=0.0)
T = []
for i, r in c.iterrows():
    if r.type in MAXR and r.n >= 2:
        gr, reste = groupes(r.n)
        par = collections.Counter(gr)
        for (t, k), m in par.items():
            T.append((i, t * m, MAXR[r.type][k], k))
        if reste > 1e-9:
            T.append((i, reste, 1.0, 0))
    else:
        T.append((i, r.n, 1.0, 0))
tr = pd.DataFrame(T, columns=["site", "n", "red", "classe"])
tr["cout"] = c.kloc.values[tr.site] * tr.red + c.racc.values[tr.site]
prod = {k: centre(t.lat, t.lon, t.gtot * t.n) for k, t in c[c.type != "gold4"].groupby("iso")}
isos = list(pays.index)
ix = {k: i for i, k in enumerate(isos)}
arcs = []
for a in isos:
    if a not in prod:
        continue
    for b in isos:
        if b == a or besoin[b] < 0.05 or b not in charge:
            continue
        t, m = route(prod[a][0], prod[a][1], charge[b][0], charge[b][1])
        if t + m <= MAX_KM and m <= MAX_MER:
            fa, fb = FA.get(a, 1.0), FA.get(b, 1.0)
            ta, tb = min(t / 2, RY.get(a, 300)), min(t / 2, RY.get(b, 300))
            terre = ta * FD.get(a, fa) + tb * FD.get(b, fb) + (t - ta - tb) * (fa + fb) / 2
            arcs.append((a, b, t, m, GW * (CONV * (0.6 + 0.2 * (fa + fb)) + terre * TERRE + m * 1.61 * (0.6 + 0.2 * (fa + fb))), 0.015 + 0.03 * (t + m) / 1000.0))
nT, nA, nV, nH = len(tr), len(arcs), len(v), len(isos)
o_A, o_V, o_Y, o_E = nT, nT + nA, nT + nA + nV, nT + nA + nV + nH
o_Hd = o_E + 4 * nH
o_Md, o_Mn = o_Hd + nH, o_Hd + 2 * nH
N = o_Hd + 3 * nH
cout = np.zeros(N)
cout[:nT] = tr.cout.values
cout[o_A:o_V] = [a[4] for a in arcs]
cout[o_V:o_Y] = v.k.values + v.racc.values
cout[o_Y:o_E] = [HYBRIDE * fp(k, "flottant") + GW * (ARRIVEE + 30 * AC + POSTE) * FA.get(k, 1.0) for k in isos]
cout[o_E:o_Hd] = [EXT_MEUR_GWH[j] * 1000 / 365 * fh(k) for j in range(4) for k in isos]
cout[o_Md:N] = 5000.0
A = lil_matrix((6 * nH, N))
bub = np.zeros(6 * nH)
for q, (i, cl) in enumerate(zip(tr.site.values, tr.classe.values)):
    k = ix[c.iso[i]]
    tw = c.twh[i]
    A[k, q] = -tw * (1 - c.nuit[i])
    A[nH + k, q] = -tw * c.nuit[i]
    A[(2 + cl) * nH + k, q] = -EXT_MAX * tw * c.stocke[i]
for q, a in enumerate(arcs):
    for row, f in ((ix[a[0]], 0.5), (nH + ix[a[0]], 0.5)):
        A[row, o_A + q] += f
    for row, f in ((ix[a[1]], 0.5), (nH + ix[a[1]], 0.5)):
        A[row, o_A + q] -= f * (1 - a[5])
for q, r in enumerate(v.itertuples()):
    k = ix[r.iso]
    A[k, o_V + q] = r.draw
    A[nH + k, o_V + q] = -r.out
for h in range(nH):
    A[h, o_Y + h] = -0.9
    A[nH + h, o_Y + h] = -0.1
    for j in range(4):
        A[h, o_E + j * nH + h] = 1.0 / REND
        A[nH + h, o_E + j * nH + h] = -1.0
        A[(2 + j) * nH + h, o_E + j * nH + h] = 1.0
    A[h, o_Hd + h] = 1.0
    A[nH + h, o_Hd + h] = -1.0
    A[h, o_Md + h] = -1.0
    A[nH + h, o_Mn + h] = -1.0
for k in isos:
    h = ix[k]
    bub[h] = -((1 - OMEGA) * besoin[k] - 0.5 * nonloc[k])
    bub[nH + h] = -(OMEGA * besoin[k] - 0.5 * nonloc[k])
bornes = [(0, n) for n in tr.n] + [(0, None)] * nA + [(0, 1)] * nV + [(0, None)] * nH + [(0, None)] * (4 * nH) + [(0, HYDRO_DEP * float(H.get(k, 0.0))) for k in isos] + [(0, None)] * (2 * nH)
res = linprog(cout, A_ub=A.tocsr(), b_ub=bub, bounds=bornes, method="highs")
x = res.x
tr["x"] = x[:nT]
f = x[o_A:o_V]
v["x"] = x[o_V:o_Y]
hy, hd, md, mn = x[o_Y:o_E], x[o_Hd:o_Md], x[o_Md:o_Mn], x[o_Mn:N]
ext4 = x[o_E:o_Hd].reshape(4, nH)
ext = ext4.sum(0)
c["pris"] = tr.groupby("site").x.sum().reindex(c.index).fillna(0.0)
def cout_site(r):
    if r.type in MAXR and r.pris >= 2 - 1e-6:
        gr, reste = groupes(r.pris)
        return sum(t * MAXR[r.type][k] for t, k in gr) * r.kloc + reste * r.kloc, len(gr), sum(t for t, k in gr)
    return r.pris * r.kloc, 0, 0
u = c[c.pris > 1e-6].copy()
cs = u.apply(cout_site, axis=1, result_type="expand")
u["k_tot"], u["n_max"], u["v14_en_max"] = cs[0], cs[1], cs[2]
mx = collections.Counter()
for r in u.itertuples():
    if r.type in MAXR and r.pris >= 2 - 1e-6:
        for t, k in groupes(r.pris)[0]:
            mx[(r.type, k)] += 1
md_cent = u.k_tot.sum() / 1000
md_cent_fr = (u.k * u.pris + u.kp * u.pris).sum() / 1000
md_racc = (u.racc * u.pris).sum() / 1000
md_exp = sum(a[4] * q for a, q in zip(arcs, f)) / 1000
vu = v[v.x > 1e-6]
md_v15 = ((vu.k + vu.racc) * vu.x).sum() / 1000
md_hyb = sum(hy[ix[k]] * HYBRIDE * fp(k, "flottant") for k in isos) / 1000
md_hyb_res = sum(hy[ix[k]] * GW * (ARRIVEE + 30 * AC + POSTE) * FA.get(k, 1.0) for k in isos) / 1000
md_ext = float((x[o_E:o_Hd] * cout[o_E:o_Hd]).sum()) / 1000
n_pacha = (u.pris * u.pacha).sum()
n_gold = FRAC_GOLD * n_pacha
n_sky = round(n_gold * 225.0 / 1126)
md_sky = n_sky * 87.2 / 1000
md_fixe = sum(q[2] for q in FIXE.values())
md_res_fixe = nonloc_twh * GW * (ARRIVEE + 30 * AC + POSTE) / 1000
total = md_cent + md_racc + md_exp + md_v15 + md_hyb + md_hyb_res + md_ext + md_sky + md_fixe + md_res_fixe
twh_sites = (u.twh * u.pris).sum()
nb = u.groupby("type").pris.sum().round().to_dict()
tw = (u.twh * u.pris).groupby(u.type).sum().round().to_dict()
stock_v15 = (vu.stock * vu.x).sum() / 1000
nuit_v14 = (u.twh * u.pris * u.stocke).sum() * 1000 / 365
nuit_v15 = (vu.out * vu.x).sum() * 1000 / 365
nuit_ext = ext.sum() * 1000 / 365
nuit_hyd = hd.sum() * 1000 / 365
print(f"statut {res.status} omega {OMEGA:.3f} manque jour {md.sum():.0f} nuit {mn.sum():.0f} TWh", {isos[i]: round(md[i] + mn[i]) for i in np.argsort(-(md + mn))[:6] if md[i] + mn[i] > 1})
print(f"centrales {nb} TWh {tw} ; sites {twh_sites:.0f} TWh ; hybrides {hy.sum():.0f} TWh ; export {f.sum():.0f} TWh")
print(f"V15 {vu.x.sum():.0f} ({stock_v15:.0f} GWh de stockage) ; extension bassins V14 {nuit_ext:.0f} GWh/j ; barrages deplaces {nuit_hyd:.0f} GWh/j ; stockage V14 {nuit_v14:.0f} GWh/j ; V15 {nuit_v15:.0f} GWh/j")
print(f"MAX {dict(mx)}")
print(f"Md€ centrales {md_cent:.0f} (au prix francais sans MAX {md_cent_fr:.0f}) ; raccordement {md_racc:.0f} ; export {md_exp:.0f} ; V15 {md_v15:.0f} ; hybrides {md_hyb:.0f} + reseau {md_hyb_res:.0f} ; bassins agrandis {md_ext:.0f} ; sky {md_sky:.0f} ; fixes {md_fixe:.0f} + reseau {md_res_fixe:.0f} ; TOTAL {total:.0f}")
R = dict(total=total, cent=md_cent, cent_prix_fr_sans_max=md_cent_fr, racc=md_racc, exp=md_exp, v15=md_v15, hyb=md_hyb, hyb_res=md_hyb_res, ext=md_ext, sky=md_sky, fixe=md_fixe, res_fixe=md_res_fixe, twh_exp=f.sum(), pertes=sum(a[5] * q for a, q in zip(arcs, f)), hyb_twh=hy.sum(), n_v15=vu.x.sum(), stock_v15_gwh=stock_v15, nuit_v14=nuit_v14, nuit_v15=nuit_v15, nuit_ext=nuit_ext, nuit_hyd=nuit_hyd, nuit_ext_seules=ext4[0].sum() * 1000 / 365, nuit_ext_max=ext4[1].sum() * 1000 / 365, nuit_ext_max2=ext4[2].sum() * 1000 / 365, nuit_ext_max3=ext4[3].sum() * 1000 / 365, manque=md.sum() + mn.sum(), n_pacha=n_pacha, n_gold=n_gold, n_sky=n_sky, omega=OMEGA)
for (t, k), m in mx.items():
    R[f"max{k}_{t}"] = m
if ECRIRE:
    pd.Series(R).round(2).to_csv(f"{D}/{NOM}_totaux.csv", header=["valeur"], lineterminator="\n")
if ECRIRE and not SOUPLE:
    o = u[["type", "lat", "lon", "pays", "ville", "pris", "twh", "k", "kloc", "k_tot", "racc", "sorte", "n_max", "v14_en_max"]].copy()
    o["twh_pris"] = o.twh * o.pris
    o.round(3).to_csv(f"{D}/plan_mondial_artesia_5_v1_recalcul_sites.csv", index=False, lineterminator="\n")
    vu[["iso", "ville", "lat", "lon", "rempart", "stock", "out", "k", "x"]].round(3).to_csv(f"{D}/plan_mondial_artesia_5_v1_recalcul_v15.csv", index=False, lineterminator="\n")
    out = collections.defaultdict(float)
    inn = collections.defaultdict(float)
    with open(f"{D}/plan_mondial_artesia_5_v1_recalcul_lignes.csv", "w", newline="", encoding="utf-8") as fh_:
        w = csv.writer(fh_, lineterminator="\n")
        w.writerow(["de", "vers", "twh_par_an", "gw", "km_terre", "km_mer", "md_euros", "lat_de", "lon_de", "lat_vers", "lon_vers"])
        for a, q in sorted(zip(arcs, f), key=lambda t: -t[1]):
            if q > 0.05:
                out[a[0]] += q
                inn[a[1]] += q * (1 - a[5])
                w.writerow([pays.pays[a[0]], pays.pays[a[1]], round(q, 1), round(q * GW, 1), round(a[2]), round(a[3]), round(a[4] * q / 1000, 1), round(prod[a[0]][0], 2), round(prod[a[0]][1], 2), round(charge[a[1]][0], 2), round(charge[a[1]][1], 2)])
    pp = (u.twh * u.pris).groupby(u.iso).sum()
    nn = u.pris.groupby(u.iso).sum()
    nv = vu.x.groupby(vu.iso).sum()
    with open(f"{D}/plan_mondial_artesia_5_v1_recalcul_pays.csv", "w", newline="", encoding="utf-8") as fh_:
        w = csv.writer(fh_, lineterminator="\n")
        w.writerow(["iso", "pays", "fossile_twh", "besoin_plan_twh", "centrales_du_plan", "production_sites_twh", "centrales_fixes_twh", "panneaux_hybrides_twh", "exporte_twh", "importe_twh", "v15", "bassins_agrandis_gwh_jour", "barrages_deplaces_gwh_jour", "manque_twh"])
        for k in pays.sort_values("fossile_twh", ascending=False).index:
            if besoin[k] < 0.05 and pp.get(k, 0) < 0.05:
                continue
            h = ix[k]
            w.writerow([k, pays.pays[k], round(F[k], 1), round(besoin[k], 1), round(nn.get(k, 0)), round(pp.get(k, 0), 1), round(nonloc[k], 1), round(hy[h], 1), round(out[k], 1), round(inn[k], 1), round(nv.get(k, 0), 1), round(ext[h] * 1000 / 365, 1), round(hd[h] * 1000 / 365, 1), round(md[h] + mn[h], 1)])