{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6f7421c0",
   "metadata": {},
   "source": [
    "# Desarrollo de producto alimenticio — línea de una startup de snacks saludables (Monte Carlo de una cola M/M/1)\n",
    "\n",
    "**Objetivo.** Simular por eventos discretos una cola M/M/1, generando los números pseudoaleatorios con un generador congruencial lineal (LCG) propio, estimar sus métricas de desempeño por Monte Carlo (muchas réplicas) y compararlas contra la solución analítica.\n",
    "\n",
    "Curso: Modelado de Sistemas bajo Incertidumbre · Universidad de los Andes · 2026-20\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e7ab4530",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "894bd64c",
   "metadata": {},
   "source": [
    "## Datos del problema\n",
    "\n",
    "Una startup de snacks saludables tiene **una sola** línea de horneado y empaque. Los pedidos llegan Poisson con $\\lambda = 6$ pedidos/hora; el tiempo de horneado + empaque es exponencial con media 8 minutos ($\\mu = 60/8 = 7.5$ pedidos/hora)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "01f4c481",
   "metadata": {},
   "source": [
    "## Generador de números pseudoaleatorios (LCG) y transformada inversa exponencial\n",
    "\n",
    "\n",
    "\n",
    "Parámetros $a=1664525$, $c=1013904223$, $m=2^{32}$, semilla $X_0=2026$. El mismo estado del LCG se usa para **todos** los números que hacen falta en la simulación completa (llegadas y servicios, de todas las réplicas), en una sola secuencia continua."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df0a7efc",
   "metadata": {},
   "outputs": [],
   "source": [
    "A_LCG, C_LCG, M_LCG = 1664525, 1013904223, 2**32\n",
    "SEMILLA_LCG = 2026\n",
    "\n",
    "def siguiente_uniforme(estado):\n",
    "    nuevo_estado = (A_LCG * estado + C_LCG) % M_LCG\n",
    "    u = nuevo_estado / M_LCG\n",
    "    return u, nuevo_estado\n",
    "\n",
    "def siguiente_exponencial(estado, tasa):\n",
    "    # TODO 1: u, \n",
    "    # TODO 2: x = \n",
    "    # TODO 3: retorna (x, estado)\n",
    "    raise NotImplementedError(\"Implementa siguiente_exponencial\")\n",
    "\n",
    "# --- prueba rápida (no la modifiques) ---\n",
    "x1, e1 = siguiente_exponencial(SEMILLA_LCG, tasa=6.0)\n",
    "x2, e2 = siguiente_exponencial(e1, tasa=6.0)\n",
    "assert x1 > 0 and x2 > 0\n",
    "assert e1 != SEMILLA_LCG and e2 != e1\n",
    "print(\"OK: siguiente_exponencial pasa la prueba básica. x1 =\", round(x1, 4), \" x2 =\", round(x2, 4))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a53732f",
   "metadata": {},
   "source": [
    "## Simular un día\n",
    "\n",
    "Implementa `simular_dia(lam, mu, horizonte, estado)` siguiendo el pseudocódigo: cada vez que necesites un número, pide uno nuevo con `siguiente_exponencial`, pasando y actualizando `estado`. Debe retornar `(llegadas, inicios, salidas, estado)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36dce7e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "def simular_dia(lam, mu, horizonte, estado):\n",
    "    \"\"\"Simula un día de cola M/M/1 usando el LCG. Retorna (llegadas, inicios, salidas, estado).\"\"\"\n",
    "    # TODO \n",
    "    # TODO \n",
    "    #         \n",
    "    # TODO 4:\n",
    "    # TODO 6: \n",
    "\n",
    "    # TODO 7: retorna (llegadas, inicios, salidas, estado)\n",
    "    raise NotImplementedError(\"Implementa simular_dia\")\n",
    "\n",
    "# --- prueba rápida (no la modifiques) ---\n",
    "llegadas_t, inicios_t, salidas_t, estado_t = simular_dia(lam=6.0, mu=7.5, horizonte=10.0, estado=SEMILLA_LCG)\n",
    "assert len(llegadas_t) == len(inicios_t) == len(salidas_t)\n",
    "assert np.all(inicios_t >= llegadas_t)      # nadie empieza antes de llegar\n",
    "assert np.all(salidas_t >= inicios_t)       # nadie sale antes de empezar\n",
    "assert estado_t != SEMILLA_LCG              # el estado del LCG avanzó\n",
    "print(f\"OK: se simularon {len(llegadas_t)} pedidos en un día de prueba.\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc7f5690",
   "metadata": {},
   "source": [
    "## Monte Carlo: repetir el día M veces\n",
    "\n",
    "Corre `simular_dia` $M=500$ veces (jornadas de 10 horas), **encadenando el estado del LCG de un día al siguiente** (no lo reinicies). De cada día, calcula el $W$ promedio y el $W_q$ promedio (en minutos)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c2e0dbe1",
   "metadata": {},
   "source": [
    "**Fórmulas:**\n",
    "\n",
    "$$\\bar W_k=\\frac{1}{n_k}\\sum_{i=1}^{n_k}(\\text{salida}_i-\\text{llegada}_i) \\qquad \\bar{Wq}_k=\\frac{1}{n_k}\\sum_{i=1}^{n_k}(\\text{inicio}_i-\\text{llegada}_i)$$\n",
    "\n",
    "$$\\hat W=\\frac{1}{M}\\sum_{k=1}^{M}\\bar W_k \\qquad SE_{\\hat W}=\\frac{s_W}{\\sqrt M} \\qquad IC_{95\\%}(\\hat W)=\\hat W\\pm1.96\\, SE_{\\hat W}$$\n",
    "\n",
    "(análogo para $\\hat{Wq}$)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aa456e3f",
   "metadata": {},
   "outputs": [],
   "source": [
    "M = 500\n",
    "horizonte = 10.0  # horas\n",
    "\n",
    "W_dias = np.empty(M)\n",
    "Wq_dias = np.empty(M)\n",
    "\n",
    "estado = SEMILLA_LCG\n",
    "# TODO: \n",
    "\n",
    "\n",
    "# TODO: W_hat\n",
    "# TODO: SE_W \n",
    "# TODO: SE_Wq \n",
    "\n",
    "# TODO: imprime W_hat y Wq_hat con su IC 95% (± 1.96*SE), y compáralos con los valores\n",
    "#       analíticos (W=40.00 min, Wq=32.00 min)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2ba3279a",
   "metadata": {},
   "source": [
    "## ¿Por qué no coinciden exactamente?\n",
    "\n",
    "Vas a notar que $W$ y $W_q$ simulados quedan **por debajo** de los valores analíticos. La fórmula de M/M/1 asume un sistema en operación continua (estado estacionario); tu simulación, en cambio, **reinicia vacía cada día** — los primeros pedidos de cada jornada casi nunca esperan. Para comprobarlo, corremos una sola simulación **larga y continua** (sin reiniciar, siguiendo con el mismo `estado` del LCG), descartando el arranque, y comparamos de nuevo."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6bc3ebc5",
   "metadata": {},
   "outputs": [],
   "source": [
    "llegadas_l, inicios_l, salidas_l, estado = simular_dia(lam=6.0, mu=7.5, horizonte=4000.0, estado=estado)\n",
    "corte = int(0.1 * len(llegadas_l))  # descarta el primer 10% como calentamiento\n",
    "\n",
    "W_largo = ((salidas_l - llegadas_l) * 60)[corte:]\n",
    "Wq_largo = ((inicios_l - llegadas_l) * 60)[corte:]\n",
    "\n",
    "print(f\"W promedio (corrida larga, sin reinicios):  {W_largo.mean():.2f} min   analítico: 40.00 min\")\n",
    "print(f\"Wq promedio (corrida larga, sin reinicios):  {Wq_largo.mean():.2f} min   analítico: 32.00 min\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "89f5cd98",
   "metadata": {},
   "source": [
    "## Conclusión\n",
    "\n",
    "_¿A qué se debe la diferencia entre la simulación de jornadas de 10 horas y la fórmula analítica? ¿Cuál de las dos crees que describe mejor un negocio real que cierra cada noche?_"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}