# -*- coding: utf-8 -*-
"""Trayectoria optima (X, Y) de la celda por correlacion G_POA ~ Isc.

Metodo (ver conversacion 2026-10-18):
  1. Base de datos: todos los logs de logs_nuevos/ con dia (G > 20 W/m2),
     ventana robusta P0.5-P99.5 para descartar outliers de la trayectoria.
  2. Normalizacion solar: rho = Isc / Isc_max(dia). Validacion del proxy:
     Pearson r y Spearman rho de Isc vs G por dia (r >= 0.7 => G es proxy).
  3. Cielo claro: indice de claridad CI = integral(G) / integral(G_clear),
     donde G_clear es la envolvente superior (max por bin de 10 min de hora
     solar) entre todos los dias. Definen la trayectoria los dias con
     CI >= 0.85*CI_max (top 3 por CI); con envolvente patchwork de pocos
     dias el CI integral no alcanza 0.8 absoluto.
  4. Tabla operativa T: h -> (X*, Y*). Por bin de 15 min de hora solar,
     mediana de X e Y de los dias claros; relleno de huecos + suavizado.
  5. Ruta espacial: Dijkstra sobre la rejilla (costo = -score, 8 vecinos)
     desde el primer punto de la tabla hasta el ultimo.
  6. Backtest por dia (incluye nublados): E_real, E_open (tabla en bucle
     abierto, muestreo bilineal de la rejilla del dia) y E_hc (tabla +
     hill-climbing de 1 paso = max del bloque 3x3). Para los dias claros la
     tabla se construye SIN el dia evaluado (leave-one-out); los nublados
     se evaluan contra la tabla global (referencia). Metrica: dE %.
  7. Salida: figura de 3 paneles (mapa + rutas, barras dE %, curva de la
     tabla) y CSV: tabla_trayectoria.csv (10 min) y resumen_backtest.csv.

Mismo mapeo y calibracion que heatmap_extremos_p80.py (adc_manager.cpp:
C=canal1 Isc, P=canal2 piranometro on-board, V=canal3 Voc; ADS122C04
24-bit, Vref 2.048 V, canal P en complemento a 2).

Uso:  py -3.12 trayectoria_optima.py
"""

import glob
import os
import re
from datetime import datetime

import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.colors import Normalize
from matplotlib.patches import FancyArrowPatch

ROOT = os.path.dirname(os.path.abspath(__file__))
LOGS_DIR = os.path.join(ROOT, "logs_nuevos")
OUT_DIR = os.path.join(ROOT, "figures", "trayectoria_optima")

ADC_MAX = 8388607.0
GAIN_SHUNT = 20.0
R_SHUNT = 0.075
DEG_PER_USTEP = 80.0 / 53000.0

DAY_G_MIN = 20.0
PAD_PCT = 0.5
T0, T1 = 7.5, 17.5          # ventana de hora solar de la tabla
BIN_TAB = 0.25              # bin de la tabla (h)
TAB_STEP = 0.1              # paso del CSV (h)
NX, NY = 61, 46             # rejilla comun

RE = re.compile(
    r"LOG:\[(?P<t>[^\]]+)\]\s*"
    r"X:(?P<x>-?\d+),Y:(?P<y>-?\d+),C:(?P<c>-?\d+),V:(?P<v>-?\d+),"
    r"P:(?P<p>-?\d+),S:(?P<s>-?\d+),er:(?P<er>-?\d+)"
)


def trapz(y, x):
    f = getattr(np, "trapezoid", None)
    if f is None:
        f = np.trapz
    return f(y, x)


def sign24(raw):
    raw &= 0xFFFFFF
    return raw - 0x1000000 if raw & 0x800000 else raw


def isc_ma(c_raw):
    vsh_mV = c_raw / ADC_MAX * 2048.0 / GAIN_SHUNT
    return (vsh_mV / 1000.0) / R_SHUNT * 1000.0


def robust_limits(arr, lo=PAD_PCT, hi=100.0 - PAD_PCT):
    a = np.asarray(arr, dtype=float)
    return float(np.percentile(a, lo)), float(np.percentile(a, hi))


def load_log(path):
    ts, xs, ys, cs, ps = [], [], [], [], []
    with open(path, "r", encoding="utf-8", errors="replace") as f:
        for line in f:
            m = RE.search(line)
            if m:
                ts.append(m.group("t"))
                xs.append(int(m.group("x")))
                ys.append(int(m.group("y")))
                cs.append(int(m.group("c")))
                ps.append(sign24(int(m.group("p"))))
    return (np.array(ts), np.array(xs, dtype=float), np.array(ys, dtype=float),
            np.array(cs, dtype=float), np.array(ps, dtype=float))


def parse_t(t_strs):
    out = np.empty(len(t_strs), dtype=float)
    for i, s in enumerate(t_strs):
        dt = datetime.strptime(s, "%Y-%m-%d %H:%M:%S")
        out[i] = dt.hour + dt.minute / 60.0 + dt.second / 3600.0
    return out


def bin_stat(x, y, vals, xbins, ybins):
    xi = np.clip(np.searchsorted(xbins, x, side="right") - 1, 0, len(xbins) - 2)
    yi = np.clip(np.searchsorted(ybins, y, side="right") - 1, 0, len(ybins) - 2)
    shape = (len(ybins) - 1, len(xbins) - 1)
    flat_sum = np.zeros(shape[0] * shape[1])
    flat_cnt = np.zeros(shape[0] * shape[1])
    for a, b, v in zip(yi, xi, vals):
        k = a * shape[1] + b
        flat_sum[k] += v
        flat_cnt[k] += 1.0
    with np.errstate(invalid="ignore"):
        flat = flat_sum / np.where(flat_cnt == 0, np.nan, flat_cnt)
    return flat.reshape(shape)


def bilinear(grid, x, y, xbins, ybins):
    """Muestreo bilineal de la rejilla en posiciones (x, y) arbitraras."""
    nx, ny = grid.shape[1], grid.shape[0]
    fx = (np.asarray(x, dtype=float) - xbins[0]) / (xbins[-1] - xbins[0]) * nx
    fy = (np.asarray(y, dtype=float) - ybins[0]) / (ybins[-1] - ybins[0]) * ny
    fx = np.clip(fx, 0, nx - 1)
    fy = np.clip(fy, 0, ny - 1)
    j0 = np.clip(fx.astype(int), 0, nx - 2)
    i0 = np.clip(fy.astype(int), 0, ny - 2)
    wx = fx - j0
    wy = fy - i0
    v = (grid[i0, j0] * (1 - wx) * (1 - wy) + grid[i0, j0 + 1] * wx * (1 - wy)
         + grid[i0 + 1, j0] * (1 - wx) * wy + grid[i0 + 1, j0 + 1] * wx * wy)
    return v


def nan_max_pool3(grid):
    """Maximo del bloque 3x3 (nan-aware): aproxima un hill-climbing de 1 paso."""
    import warnings
    p = np.pad(grid, 1, mode="edge")
    n, m = grid.shape
    with np.errstate(invalid="ignore"), warnings.catch_warnings():
        warnings.simplefilter("ignore", RuntimeWarning)
        out = np.nanmax(
            np.stack([p[di:di + n, dj:dj + m] for di in range(3) for dj in range(3)]),
            axis=0)
    return out


def dijkstra(grid, i0, j0, i1, j1):
    """Ruta de maximo score (costo = -score) con 8 vecinos en la rejilla."""
    import heapq
    n, m = grid.shape
    INF = np.inf
    best = np.full((n, m), INF)
    prev = np.full((n, m), -1, dtype=int)
    best[i0, j0] = 0.0
    pq = [(0.0, i0, j0)]
    nbrs = [(-1, -1, 1.0), (-1, 0, 1.414), (-1, 1, 1.414), (0, -1, 1.0),
            (0, 1, 1.0), (1, -1, 1.414), (1, 0, 1.0), (1, 1, 1.414)]
    while pq:
        d, i, j = heapq.heappop(pq)
        if d > best[i, j] + 1e-12:
            continue
        if i == i1 and j == j1:
            break
        for di, dj, w in nbrs:
            ni, nj = i + di, j + dj
            if ni < 0 or ni >= n or nj < 0 or nj >= m:
                continue
            if not np.isfinite(grid[ni, nj]):
                continue
            nd = d + w * (0.0 if grid[ni, nj] >= 1 else 1.0 - grid[ni, nj])
            if nd < best[ni, nj] - 1e-12:
                best[ni, nj] = nd
                prev[ni, nj] = i * m + j
                heapq.heappush(pq, (nd, ni, nj))
    if not np.isfinite(best[i1, j1]):
        return None
    path = []
    cur = i1 * m + j1
    while cur != -1:
        i, j = divmod(cur, m)
        path.append((i, j))
        cur = prev[i, j]
    path.reverse()
    return path


def build_table(days_subset, tb):
    """Tabla T: centro de cada bin de tb -> (X*, Y*). Mediana por bin +
    relleno de huecos + suavizado. Devuelve None si hay menos de 4 bins
    con datos."""
    tab_x = np.full(len(tb) - 1, np.nan)
    tab_y = np.full(len(tb) - 1, np.nan)
    for k in range(len(tb) - 1):
        pts_x, pts_y = [], []
        for d in days_subset:
            m = (d["t_solar"] >= tb[k]) & (d["t_solar"] < tb[k + 1])
            pts_x.extend(d["x"][m].tolist())
            pts_y.extend(d["y"][m].tolist())
        if len(pts_x) >= 3:
            tab_x[k] = np.median(pts_x)
            tab_y[k] = np.median(pts_y)
    tf = (tb[:-1] + tb[1:]) / 2
    fin = np.isfinite(tab_x)
    if fin.sum() < 4:
        return None
    tab_x = np.interp(tf, tf[fin], tab_x[fin])
    tab_y = np.interp(tf, tf[fin], tab_y[fin])
    # suavizado con nucleo [1,2,3,2,1]/9
    ker = np.array([1, 2, 3, 2, 1], dtype=float) / 9
    tab_x = np.convolve(tab_x, ker, mode="same")
    tab_y = np.convolve(tab_y, ker, mode="same")
    return tf, tab_x, tab_y


def cell_index(x, y, xbins, ybins):
    nx, ny = len(xbins) - 1, len(ybins) - 1
    j = int(np.clip((x - xbins[0]) / (xbins[-1] - xbins[0]) * nx, 0, nx - 1))
    i = int(np.clip((y - ybins[0]) / (ybins[-1] - ybins[0]) * ny, 0, ny - 1))
    return i, j


def main():
    os.makedirs(OUT_DIR, exist_ok=True)
    paths = sorted(glob.glob(os.path.join(LOGS_DIR, "log_*.txt")),
                   key=lambda p: re.search(r"(\d{4}-\d{2}-\d{2})", p).group(1))
    print("Logs disponibles:")
    days = []
    for p in paths:
        ts, xs, ys, cs, ps = load_log(p)
        g = np.maximum(0.0, ps) / ADC_MAX * 2048.0 / 2.0
        isc = isc_ma(cs)
        day = g > DAY_G_MIN
        if day.sum() < 50:
            print("   %-28s sin dia suficiente, se omite" % os.path.basename(p))
            continue
        x0, x1 = robust_limits(xs[day])
        y0, y1 = robust_limits(ys[day])
        ins = day & (xs >= x0) & (xs <= x1) & (ys >= y0) & (ys <= y1)
        if ins.sum() < 50:
            print("   %-28s se omite" % os.path.basename(p))
            continue
        t = parse_t(ts[ins])
        # mediodia solar: maximo de G con media movil de 1 h
        k = max(3, int(np.round(1.0 / max(1e-6, np.median(np.diff(t))))))
        kwin = 2 * k + 1
        kern = np.ones(kwin) / kwin
        gg_pad = np.pad(g[ins], k, mode="edge")
        gg_s = np.convolve(gg_pad, kern, mode="valid")[:len(g[ins])]
        t_noon = t[int(np.argmax(gg_s))]
        t_solar = t - t_noon + 12.0
        days.append(dict(
            name=os.path.splitext(os.path.basename(p))[0],
            x=xs[ins], y=ys[ins], c=isc[ins], g=g[ins],
            t=t, t_solar=t_solar,
            x0=x0, x1=x1, y0=y0, y1=y1,
        ))
        print("   %-28s n_dia=%d" % (os.path.basename(p), int(ins.sum())))
    if len(days) < 2:
        print("No hay suficientes dias de datos.")
        return

    # rejilla comun sobre la union de las ventanas de los dias
    xlo, xhi = min(d["x0"] for d in days), max(d["x1"] for d in days)
    ylo, yhi = min(d["y0"] for d in days), max(d["y1"] for d in days)
    xbins = np.linspace(xlo, xhi, NX)
    ybins = np.linspace(ylo, yhi, NY)

    # --- paso 1: metricas por dia y rejillas -------------------------------
    bins_clear = np.arange(7.0, 18.0 + 1e-9, 1.0 / 6.0)  # bins de 10 min
    clear_vals = np.full(len(bins_clear) - 1, np.nan)
    for d in days:
        d["rho"] = d["c"] / d["c"].max()
        d["grid_rho"] = bin_stat(d["x"], d["y"], d["rho"], xbins, ybins)
        d["grid_g"] = bin_stat(d["x"], d["y"], d["g"], xbins, ybins)
        sel = (d["t_solar"] >= 7.0) & (d["t_solar"] <= 18.0)
        for k in range(len(bins_clear) - 1):
            m = sel & (d["t_solar"] >= bins_clear[k]) & (d["t_solar"] < bins_clear[k + 1])
            if m.any():
                clear_vals[k] = np.nanmax([clear_vals[k], d["g"][m].max()])
    # interpolacion de huecos de la envolvente
    fin = np.isfinite(clear_vals)
    clear_vals = np.interp(np.arange(len(clear_vals)), np.where(fin)[0], clear_vals[fin])
    tc = (bins_clear[:-1] + bins_clear[1:]) / 2
    for d in days:
        g_clear_d = np.interp(d["t_solar"], tc, clear_vals)
        # CI por integral: dias con cielo claro capturan mas energia total
        with np.errstate(invalid="ignore"):
            d["CI"] = float(trapz(d["g"], d["t_solar"]) /
                             max(1e-9, trapz(g_clear_d, d["t_solar"])))
        if d["c"].std() > 0 and d["g"].std() > 0:
            d["r"] = float(np.corrcoef(d["c"], d["g"])[0, 1])
        else:
            d["r"] = np.nan
        cr = np.argsort(np.argsort(d["c"]))
        gr = np.argsort(np.argsort(d["g"]))
        d["rho_s"] = float(np.corrcoef(cr, gr)[0, 1])
    ci_max = max(d["CI"] for d in days)
    # umbral relativo: la envolvente es un patchwork con pocos dias, el CI
    # integral nunca llega a 0.8 absoluto; el dia mas claro queda ~0.7-0.9
    clear = sorted([d for d in days if d["CI"] >= 0.85 * ci_max],
                   key=lambda d: d["CI"], reverse=True)[:3]
    print("\nDias que definen la trayectoria (CI >= 0.85*CI_max = %.2f): %s"
          % (0.85 * ci_max,
             ", ".join("%s(%.2f)" % (d["name"], d["CI"]) for d in clear)))
    if not clear:
        print("Ningun dia claro; se usa el de mayor CI.")
        clear = [max(days, key=lambda d: d["CI"])]

    # --- paso 2: tabla T: h -> (X*, Y*) ------------------------------------
    tb = np.arange(T0, T1 + 1e-9, BIN_TAB)
    res = build_table(clear, tb)
    if res is None:
        print("Tabla sin datos suficientes.")
        return
    tab_t, tab_x, tab_y = res

    # tabla a paso de 10 min para el CSV
    csv_t = np.arange(tab_t[0], tab_t[-1] + 1e-9, TAB_STEP)
    csv_x = np.interp(csv_t, tab_t, tab_x)
    csv_y = np.interp(csv_t, tab_t, tab_y)
    csv_path = os.path.join(OUT_DIR, "tabla_trayectoria.csv")
    with open(csv_path, "w", encoding="utf-8") as f:
        f.write("hora_solar,X_ref_upaso,Y_ref_upaso,X_ref_deg,Y_ref_deg\n")
        for a, b, c in zip(csv_t, csv_x, csv_y):
            f.write("%.1f,%.0f,%.0f,%.2f,%.2f\n"
                    % (a, b, c, b * DEG_PER_USTEP, c * DEG_PER_USTEP))

    # --- paso 3: rejilla pool y Dijkstra ------------------------------------
    import warnings
    with warnings.catch_warnings():
        warnings.simplefilter("ignore", RuntimeWarning)
        pool_rho = np.nanmean(np.stack([d["grid_rho"] for d in clear]), axis=0)
        pool_g = np.nanmean(np.stack([d["grid_g"] for d in clear]), axis=0)
        pool_rho_n = pool_rho / np.nanmax(pool_rho)
    i0, j0 = cell_index(tab_x[0], tab_y[0], xbins, ybins)
    i1, j1 = cell_index(tab_x[-1], tab_y[-1], xbins, ybins)
    path = dijkstra(pool_rho_n, i0, j0, i1, j1)
    if path:
        pcx = (xbins[np.array([j for i, j in path])] +
               xbins[np.array([j for i, j in path]) + 1]) / 2
        pcy = (ybins[np.array([i for i, j in path])] +
               ybins[np.array([i for i, j in path]) + 1]) / 2
    else:
        pcx, pcy = np.array([]), np.array([])

    # --- paso 4: backtest por dia -------------------------------------------
    rows = []
    for d in days:
        is_clear = any(x is d for x in clear)
        g_hc = nan_max_pool3(d["grid_rho"])
        m = (d["t_solar"] >= tab_t[0]) & (d["t_solar"] <= tab_t[-1])
        tt = d["t_solar"][m]
        e_real = float(trapz(d["rho"][m], tt)) if len(tt) > 1 else 0.0
        if is_clear and len(clear) >= 2:
            # leave-one-out: tabla construida sin el dia evaluado
            r_loo = build_table([x for x in clear if x is not d], tb)
            if r_loo is None:
                r_loo = (tab_t, tab_x, tab_y)
                tipo = "self-fit"
            else:
                tipo = "LOO"
        elif is_clear:
            r_loo = (tab_t, tab_x, tab_y)
            tipo = "self-fit"
        else:
            r_loo = (tab_t, tab_x, tab_y)
            tipo = "referencia"
        u_t, u_x, u_y = r_loo
        lo = max(tt.min(), u_t[0]) if len(tt) > 0 else None
        hi = min(tt.max(), u_t[-1]) if len(tt) > 0 else None
        if len(tt) < 2 or lo is None or lo >= hi:
            e_open = 0.0
            e_hc = 0.0
        else:
            tx_t = u_t[(u_t >= lo) & (u_t <= hi)]
            tx_x = np.interp(tx_t, u_t, u_x)
            tx_y = np.interp(tx_t, u_t, u_y)
            g_num = np.nan_to_num(d["grid_rho"], nan=0.0)
            g_hc_num = np.nan_to_num(g_hc, nan=0.0)
            if len(tx_t) < 2:
                e_open = 0.0
                e_hc = 0.0
            else:
                e_open = float(trapz(bilinear(g_num, tx_x, tx_y, xbins, ybins), tx_t))
                e_hc = float(trapz(bilinear(g_hc_num, tx_x, tx_y, xbins, ybins), tx_t))
        de_open = (e_open - e_real) / e_real * 100.0 if e_real > 0 else np.nan
        de_hc = (e_hc - e_real) / e_real * 100.0 if e_real > 0 else np.nan
        rows.append((d["name"], d["CI"], d["r"], d["rho_s"], e_real,
                     de_open, de_hc, "claro" if is_clear else "", tipo))
        print("%-22s CI=%.2f  r=%.3f  rho_s=%.3f  dE_open=%+.1f%%  dE_hc=%+.1f%%  [%s] %s"
              % (d["name"].replace("log_", ""), d["CI"], d["r"], d["rho_s"],
                 de_open, de_hc, tipo, "<-- claro" if is_clear else ""))
    res_path = os.path.join(OUT_DIR, "resumen_backtest.csv")
    with open(res_path, "w", encoding="utf-8") as f:
        f.write("log,CI,r_isc_g,spearman,E_real,dE_open_pct,dE_hc_pct,rol,tipo\n")
        for r in rows:
            f.write("%s,%.3f,%.3f,%.3f,%.4f,%.1f,%.1f,%s,%s\n"
                    % (r[0], r[1], r[2], r[3], r[4], r[5], r[6], r[7], r[8]))

    # --- paso 5: figura de 3 paneles -----------------------------------------
    fig = plt.figure(figsize=(18.5, 11.0))
    gs = fig.add_gridspec(2, 2, width_ratios=[1.7, 1.0],
                          left=0.05, right=0.985, top=0.845, bottom=0.09,
                          wspace=0.22, hspace=0.30)
    axA = fig.add_subplot(gs[:, 0])
    axB = fig.add_subplot(gs[0, 1])
    axC = fig.add_subplot(gs[1, 1])

    im = axA.imshow(
        np.ma.masked_invalid(pool_rho),
        extent=[xbins[0], xbins[-1], ybins[0], ybins[-1]],
        origin="lower", aspect="auto", cmap="inferno",
        norm=Normalize(vmin=0, vmax=np.nanmax(pool_rho)),
        interpolation="nearest")
    cb = fig.colorbar(im, ax=axA, pad=0.02)
    cb.set_label("Isc media por celda, normalizada por dia (pool de dias claros)",
                 fontsize=9)
    cb.ax.tick_params(labelsize=8)
    cont_g = axA.contour((xbins[:-1] + xbins[1:]) / 2, (ybins[:-1] + ybins[1:]) / 2,
                         pool_g, levels=5, colors="white", linewidths=0.7, alpha=0.6)
    axA.clabel(cont_g, inline=True, fontsize=6.5, fmt="%.0f")

    cmap = plt.get_cmap("tab10")
    for k, d in enumerate(days):
        axA.plot(d["x"], d["y"], ".", ms=1.4, color=cmap(k % 10), alpha=0.22,
                 zorder=2, label=d["name"].replace("log_", "10-"))

    axA.plot(tab_x, tab_y, "-", color="white", lw=2.6, zorder=5,
             label="Tabla operativa T: h -> (X*, Y*)")
    step = max(1, len(tab_x) // 12)
    for k in range(0, len(tab_x) - 1, step):
        axA.add_patch(FancyArrowPatch((tab_x[k], tab_y[k]), (tab_x[k + 1], tab_y[k + 1]),
                                      arrowstyle="-|>", mutation_scale=14,
                                      color="white", lw=1.2, alpha=0.9, zorder=5))
    for h in range(8, 18):
        if h < tab_t[0] or h > tab_t[-1]:
            continue
        k = int(np.argmin(np.abs(tab_t - h)))
        axA.plot(tab_x[k], tab_y[k], "o", ms=7, mfc="gold", mec="black",
                 mew=0.8, zorder=7)
        axA.annotate(str(h), (tab_x[k], tab_y[k]), textcoords="offset points",
                     xytext=(5, 5), fontsize=7.5, color="gold", zorder=7)
    if len(pcx):
        axA.plot(pcx, pcy, "--", color="yellow", lw=1.8, zorder=4,
                 label="Ruta de maximo score (Dijkstra)")
    axA.set_title("A) Trayectoria optima sobre el mapa de Isc (pool de dias claros)",
                  fontsize=11)
    axA.set_xlabel("X (micro-pasos)")
    axA.set_ylabel("Y (micro-pasos)")
    axA.legend(loc="upper center", bbox_to_anchor=(0.5, -0.12), ncol=2,
               fontsize=7.5, framealpha=0.95, borderpad=0.4, handletextpad=0.5,
               labelspacing=0.4, columnspacing=1.0)
    tx = axA.secondary_xaxis("top", functions=(lambda s: s * DEG_PER_USTEP,
                                               lambda d: d / DEG_PER_USTEP))
    tx.set_xlabel("X [deg]", fontsize=8, labelpad=2)
    tx.tick_params(labelsize=7)
    ty = axA.secondary_yaxis("right", functions=(lambda s: s * DEG_PER_USTEP,
                                                 lambda d: d / DEG_PER_USTEP))
    ty.set_ylabel("Y [deg]", fontsize=8, labelpad=2, color="dimgray")
    ty.tick_params(labelsize=7, colors="dimgray")

    names = [r[0].replace("log_", "") for r in rows]
    xx = np.arange(len(names))
    w = 0.38
    axB.bar(xx - w / 2, [r[5] for r in rows], w, color="steelblue",
            label="tabla, bucle abierto")
    axB.bar(xx + w / 2, [r[6] for r in rows], w, color="darkorange",
            label="tabla + hill-climbing 1 paso")
    axB.axhline(0, color="black", lw=0.9)
    for k, r in enumerate(rows):
        axB.annotate("CI=%.2f r=%.2f" % (r[1], r[2]),
                     (k, max(2.0, max(r[5], r[6])) + 2.0),
                     ha="center", fontsize=7)
    axB.set_xticks(xx)
    axB.set_xticklabels(names, rotation=30, ha="right", fontsize=8)
    axB.set_ylabel("dE vs registrado [%]")
    axB.set_title("B) Backtest por dia (positivo = la tabla captura mas)\n"
                  "dias claros evaluados con leave-one-out",
                  fontsize=11)
    axB.grid(axis="y", ls="--", alpha=0.35)
    axB.legend(fontsize=8, loc="upper left")

    axC.plot(csv_t, csv_x, "-", color="steelblue", lw=2.0, label="X*(h) [upaso]")
    axC.set_xlabel("hora solar")
    axC.set_ylabel("X* [upaso]", color="steelblue")
    axC.tick_params(axis="y", labelcolor="steelblue")
    axC2 = axC.twinx()
    axC2.plot(csv_t, csv_y, "--", color="darkorange", lw=2.0, label="Y*(h) [upaso]")
    axC2.set_ylabel("Y* [upaso]", color="darkorange")
    axC2.tick_params(axis="y", labelcolor="darkorange")
    vxp = float(np.abs(np.gradient(csv_x)).max()) * 60.0
    vyp = float(np.abs(np.gradient(csv_y)).max()) * 60.0
    axC.set_title("C) Tabla operativa (para firmware, 10 min)\n"
                  "velocidad maxima: X %.0f  Y %.0f upaso/min" % (vxp, vyp),
                  fontsize=11)
    axC.grid(True, ls="--", alpha=0.35)
    h1, l1 = axC.get_legend_handles_labels()
    h2, l2 = axC2.get_legend_handles_labels()
    axC.legend(h1 + h2, l1 + l2, fontsize=8, loc="upper left")

    fig.suptitle(
        "Trayectoria optima de la celda por correlacion G_POA ~ Isc  -  "
        "%d logs, %d dias claros (CI >= 0.85*CI_max)\n"
        "logica hibrida por hora del dia: setpoint de la tabla T + hill-climbing "
        "sobre Isc solo con cielo claro; sin busqueda local con nubes"
        % (len(days), len(clear)),
        fontsize=13, y=0.99)
    out = os.path.join(OUT_DIR, "trayectoria_optima_resumen.png")
    fig.savefig(out, dpi=150, bbox_inches="tight")
    plt.close(fig)
    print("\nTabla    -> %s" % os.path.basename(csv_path))
    print("Figura   -> %s" % os.path.basename(out))
    print("Backtest -> %s" % os.path.basename(res_path))


if __name__ == "__main__":
    main()