"""Estudio de reposición de Arume (noviembre de 2026): simulación reproducible.

Compara tres formas de reponer el stock de muchos puntos de venta desde un almacén:
  1. Mínimo fijo: si el stock baja de un mínimo, se repone hasta un máximo (dos mínimos).
  2. Cobertura fija: cada semana se repone hasta N semanas de la venta media de las últimas 4.
  3. Previsión: cada semana se repone hasta la previsión del plazo más un colchón calculado
     con el error real de esa previsión (y la temporada del año anterior).
Con demanda sintética de tres tipos (estable, estacional e intermitente), venta perdida
cuando no hay stock, revisión semanal y un plazo de una semana que a veces son dos.

Para cada política se recorre una escala de «generosidad» (k) y se mide, en el segundo año,
la tasa de servicio (unidades servidas / demanda) y el stock medio en semanas de venta.
Así se compara con justicia: cuánto stock necesita cada política para servir lo mismo.

Uso: python sim.py <carpeta_salida>   (semilla fija: siempre da los mismos números)
"""
import json, sys, math
import numpy as np

OUT = sys.argv[1]
SEM = 20261102
rng = np.random.default_rng(SEM)
N = 400          # series producto-punto por tipo de demanda
T = 104          # semanas: la 1.ª mitad es histórico, se evalúa la 2.ª
P_PLAZO2 = 0.2   # probabilidad de que el envío tarde dos semanas en vez de una


def demanda(tipo):
    t = np.arange(T)
    if tipo == 'estable':
        lam = rng.uniform(8, 60, size=(N, 1)) * np.ones((1, T))
        return rng.poisson(lam)
    if tipo == 'estacional':
        base = rng.uniform(6, 40, size=(N, 1))
        pico = rng.uniform(22, 32, size=(N, 1))  # semana del pico (verano)
        amp = rng.uniform(0.4, 0.8, size=(N, 1))
        lam = base * (1 + amp * np.cos(2 * np.pi * (t[None, :] - pico) / 52))
        # sobredispersión: binomial negativa con la misma media
        r = 6.0
        return rng.negative_binomial(r, r / (r + np.maximum(lam, 0.05)))
    if tipo == 'intermitente':
        p = rng.uniform(0.12, 0.4, size=(N, 1))
        hay = rng.random((N, T)) < p
        tam = 1 + rng.poisson(rng.uniform(1, 4, size=(N, 1)) * np.ones((1, T)))
        return hay * tam
    raise ValueError(tipo)


def simula(D, politica, k, plazos, H):
    """Devuelve (servicio, stock medio en semanas, envíos por serie y año)."""
    n = D.shape[0]
    stock = np.zeros(n)
    llegadas = np.zeros((n, T + 8))
    ventas = np.zeros(D.shape)
    serv = dem = st_acum = env = 0.0
    prev_antes = {}
    errores = np.full((n, T), np.nan)
    filas = np.arange(n)
    # El primer año es histórico conocido: las ventas fueron la demanda. Se arranca con dos
    # semanas de stock y se simula el segundo año con cada política.
    ventas[:, :52] = D[:, :52]
    stock = 2 * D[:, :52].mean(axis=1)
    if politica == 'prevision':
        for w in range(2, 52):
            if w - H - 7 >= 0:
                prev_h = ventas[:, w - H - 7:w - H + 1].mean(axis=1) * H
                errores[:, w] = ventas[:, w - H + 1:w + 1].sum(axis=1) - prev_h
    for w in range(52, T):
        stock += llegadas[:, w]
        vend = np.minimum(stock, D[:, w])
        ventas[:, w] = vend
        stock -= vend
        if w >= 52:
            serv += vend.sum(); dem += D[:, w].sum(); st_acum += stock.sum()
        pos = stock + llegadas[:, w + 1:].sum(axis=1)
        if politica == 'minimo':
            media = ventas[:, max(0, w - 51):w + 1].mean(axis=1)
            minimo = np.ceil(k * H * media)
            objetivo = np.where(pos < minimo, 2 * minimo, pos)
        elif politica == 'cobertura':
            objetivo = np.ceil(k * H * ventas[:, max(0, w - 3):w + 1].mean(axis=1))
        else:
            nivel = ventas[:, max(0, w - 7):w + 1].mean(axis=1)
            if w >= 60:
                ref = ventas[:, w - 51:w - 51 + H].mean(axis=1)      # el año pasado, las próximas semanas
                ult = ventas[:, w - 59:w - 51].mean(axis=1)          # el año pasado, estas 8 semanas
                factor = np.clip(np.where(ult > 0.5, ref / np.maximum(ult, 0.5), 1.0), 0.4, 2.5)
            else:
                factor = np.ones(n)
            prev = nivel * factor * H
            if w >= 2:
                errores[:, w] = ventas[:, w - H + 1:w + 1].sum(axis=1) - prev_antes[w - H + 1] if w - H + 1 >= 52 else errores[:, w]
            prev_antes[w] = prev
            ventana = errores[:, max(0, w - 25):w + 1]
            cuenta = np.sum(~np.isnan(ventana), axis=1)
            sigma = np.where(cuenta >= 6, np.nanstd(np.where(cuenta[:, None] > 0, ventana, 0), axis=1), np.sqrt(np.maximum(prev, 1)))
            objetivo = np.ceil(prev + k * sigma)
        # nadie deja a cero algo que vende: como mínimo, una unidad en cada sitio
        vende = ventas[:, max(0, w - 51):w + 1].sum(axis=1) > 0
        objetivo = np.where(vende, np.maximum(objetivo, 1), objetivo)
        pedido = np.maximum(0, objetivo - pos)
        hay = pedido > 0
        llegadas[filas[hay], w + plazos[hay, w]] += pedido[hay]
        if w >= 52:
            env += hay.sum()
    venta_media = dem / (n * (T - 52))
    stock_semanas = (st_acum / (n * (T - 52))) / max(venta_media, 1e-9)
    return serv / dem, stock_semanas, env / n


ESCALAS = {
    'minimo': [0.4, 0.55, 0.7, 0.85, 1.0, 1.2, 1.45, 1.75, 2.1, 2.5, 3.0, 3.6, 4.3, 5.2],
    'cobertura': [0.8, 1.0, 1.2, 1.4, 1.6, 1.8, 2.1, 2.5, 3.0, 3.6, 4.3, 5.2, 6.2],
    'prevision': [0.0, 0.5, 1.0, 1.3, 1.65, 2.0, 2.4, 2.9, 3.5, 4.2, 5.0, 6.0],
}
NOMBRES = {'minimo': 'Mínimo fijo', 'cobertura': 'Cobertura fija', 'prevision': 'Previsión'}
# Plazo de entrega: corto (el almacén repone a la tienda en una semana, a veces dos) y largo
# (compra a fábrica o importación: cuatro semanas, a veces tres o cinco).
PLAZOS = {
    'corto': {'H': 2, 'sorteo': lambda shape: (rng.random(shape) < P_PLAZO2).astype(int) + 1},
    'largo': {'H': 5, 'sorteo': lambda shape: rng.choice([3, 4, 5], size=shape, p=[0.25, 0.5, 0.25])},
}

res = {'semilla': SEM, 'series_por_tipo': N, 'semanas': T, 'escenarios': {}}
DEMANDAS = {tipo: demanda(tipo) for tipo in ['estable', 'estacional', 'intermitente']}
for (plazo, cfg), tipo in [((p, c), t) for p, c in PLAZOS.items() for t in DEMANDAS]:
    D = DEMANDAS[tipo]
    H = cfg['H']
    plazos = cfg['sorteo'](D.shape)  # el mismo sorteo para todas las políticas
    curvas = {}
    for pol, ks in ESCALAS.items():
        pts = []
        for k in ks:
            s, st, e = simula(D, pol, k, plazos, H)
            pts.append({'k': k, 'servicio': round(s, 4), 'stock_semanas': round(st, 3), 'envios_anio': round(e, 1)})
        curvas[pol] = pts
        print(plazo, tipo, pol, [(p['servicio'], p['stock_semanas']) for p in pts], flush=True)

    def stock_para(pol, objetivo):
        pts = sorted(curvas[pol], key=lambda p: p['servicio'])
        for a, b in zip(pts, pts[1:]):
            if a['servicio'] <= objetivo <= b['servicio'] and b['servicio'] > a['servicio']:
                f = (objetivo - a['servicio']) / (b['servicio'] - a['servicio'])
                return round(a['stock_semanas'] + f * (b['stock_semanas'] - a['stock_semanas']), 2)
        return None

    def servicio_con(pol, stock):
        pts = sorted(curvas[pol], key=lambda p: p['stock_semanas'])
        for a, b in zip(pts, pts[1:]):
            if a['stock_semanas'] <= stock <= b['stock_semanas'] and b['stock_semanas'] > a['stock_semanas']:
                f = (stock - a['stock_semanas']) / (b['stock_semanas'] - a['stock_semanas'])
                return round(a['servicio'] + f * (b['servicio'] - a['servicio']), 4)
        return None

    resumen = {}
    for obj in (0.95, 0.98):
        resumen[f'stock_para_{int(obj*100)}'] = {pol: stock_para(pol, obj) for pol in ESCALAS}
    for st in (1.5, 2.0):
        resumen[f'servicio_con_{st}_semanas'] = {pol: servicio_con(pol, st) for pol in ESCALAS}
    res['escenarios'][f'{tipo}|{plazo}'] = {'demanda': tipo, 'plazo': plazo, 'curvas': curvas, 'resumen': resumen, 'venta_media_semanal': round(float(D[:, 52:].mean()), 2)}
    print(plazo, tipo, json.dumps(resumen), flush=True)

json.dump(res, open(f'{OUT}/estudio-reposicion-2026.json', 'w'), ensure_ascii=False, indent=1)
with open(f'{OUT}/estudio-reposicion-2026.csv', 'w') as f:
    f.write('tipo_demanda,plazo,politica,k,tasa_servicio,stock_medio_semanas,envios_por_serie_y_anio\n')
    for clave, d in res['escenarios'].items():
        for pol, pts in d['curvas'].items():
            for p in pts:
                f.write(f"{d['demanda']},{d['plazo']},{NOMBRES[pol]},{p['k']},{p['servicio']},{p['stock_semanas']},{p['envios_anio']}\n")
print('ok')