library(shiny)

# ============================================================
#  INTEGRACIÓN NUMÉRICA
# ============================================================
integrar <- function(f, a, b, n = 5000) {
  h  <- (b - a) / n
  xs <- a + (0:n) * h
  ys <- f(xs)
  ys[c(1, length(ys))] <- ys[c(1, length(ys))] / 2
  sum(ys) * h
}

# ============================================================
#  MODELOS  (l0 = 100 fijo en todos)
# ============================================================
L0 <- 100

S_moivre  <- function(x, omega) 1 - x / omega
mu_moivre <- function(x, omega) 1 / (omega - x)

S_exp  <- function(x, mu) exp(-mu * x)
mu_exp <- function(x, mu) mu

S_cuad  <- function(x, omega) 1 - (x / omega)^2
mu_cuad <- function(x, omega) (2 * x) / (omega^2 - x^2)

construir_modelo <- function(tipo, omega, mu) {
  switch(tipo,
    base = list(
      nombre    = "Base: l(x) = 100 - x",
      S         = function(x) S_moivre(x, 100),
      mu        = function(x) mu_moivre(x, 100),
      omega     = 100,
      paso1_tex = "$$ S(x) = 1 - \\frac{x}{100}, \\qquad \\mu(x) = \\frac{1}{100 - x} $$"
    ),
    moivre = list(
      nombre    = paste0("Moivre: 100·(1 − x/", omega, ")"),
      S         = function(x) S_moivre(x, omega),
      mu        = function(x) mu_moivre(x, omega),
      omega     = omega,
      paso1_tex = paste0("$$ S(x) = 1 - \\frac{x}{", omega,
                         "}, \\qquad \\mu(x) = \\frac{1}{", omega, " - x} $$")
    ),
    exp = list(
      nombre    = paste0("Exponencial: 100·e^(−", sprintf("%.4f", mu), "x)"),
      S         = function(x) S_exp(x, mu),
      mu        = function(x) mu_exp(x, mu),
      omega     = Inf,
      paso1_tex = paste0("$$ S(x) = e^{-", sprintf("%.4f", mu),
                         "x}, \\qquad \\mu(x) = ", sprintf("%.4f", mu), " $$")
    ),
    cuad = list(
      nombre    = paste0("Cuadrático: 100·(1 − x²/", omega, "²)"),
      S         = function(x) S_cuad(x, omega),
      mu        = function(x) mu_cuad(x, omega),
      omega     = omega,
      paso1_tex = paste0("$$ S(x) = 1 - \\frac{x^2}{", omega,
                         "^2}, \\qquad \\mu(x) = \\frac{2x}{", omega,
                         "^2 - x^2} $$")
    )
  )
}

# ============================================================
#  UI
# ============================================================
ui <- fluidPage(
  tags$head(
    tags$link(rel = "stylesheet", type = "text/css", href = "estilo.css"),

    tags$script(HTML("
      window.MathJax = {
        tex: {
          inlineMath: [['$','$'], ['\\\\(','\\\\)']],
          displayMath: [['$$','$$'], ['\\\\[','\\\\]']]
        },
        svg: { fontCache: 'global' },
        startup: { typeset: false }
      };
    ")),
    tags$script(src = "https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-svg.js"),

    tags$script(HTML("
      (function() {
        function reTipografiar() {
          if (window.MathJax && MathJax.typesetPromise) {
            MathJax.typesetPromise().catch(function(){});
          }
        }
        function esperarMathJax() {
          if (window.MathJax && MathJax.startup && MathJax.startup.promise) {
            MathJax.startup.promise.then(function() {
              var obs = new MutationObserver(function(muts) {
                for (var i = 0; i < muts.length; i++) {
                  var t = muts[i].target;
                  if (t.closest && t.closest('mjx-container')) return;
                }
                clearTimeout(window.__mjxTimer);
                window.__mjxTimer = setTimeout(reTipografiar, 50);
              });
              obs.observe(document.body, { childList: true, subtree: true });
              reTipografiar();
            });
          } else {
            setTimeout(esperarMathJax, 100);
          }
        }
        esperarMathJax();
      })();
    "))
  ),

  div(class = "container",
    div(class = "header",
      h1(HTML("Esperanza de vida conjunta · \\(e_{xy}\\) y \\(e_{\\overline{xy}}\\)"),
         style = "font-size: 20px; margin: 8px 0; color: #1f6feb;"),
      h2("JuanMdeLejarza",
         style = "font-size: 10px; margin-top: 8px; margin-bottom: 8px; color: #1f6feb;"),
      actionButton("generar", HTML("➕ Nuevo par de edades"), class = "btn-generar")
    ),

    div(class = "enunciado",
      span(class = "contexto", "1. Elige el modelo de supervivencia  (l₀ = 100)"),
      fluidRow(
        column(6,
          selectInput("modelo", "Modelo:",
            choices = c("Base: l(x) = 100 - x"        = "base",
                        "Moivre: 100·(1 − x/ω)"       = "moivre",
                        "Exponencial: 100·e^(−μx)"    = "exp",
                        "Cuadrático: 100·(1 − x²/ω²)" = "cuad"),
            selected = "base")
        ),
        column(6,
          conditionalPanel("input.modelo == 'moivre' || input.modelo == 'cuad'",
            numericInput("omega", "ω (edad máxima):", value = 100, min = 2, step = 1)),
          conditionalPanel("input.modelo == 'exp'",
            numericInput("mu", "μ (tanto constante):", value = 0.02,
                         min = 0.0001, step = 0.001))
        )
      ),
      div(class = "nota-modelo",
        HTML(paste0(
          "<strong>Funciones por modelo:</strong><br>",
          "• <strong>Base / Moivre</strong>: ",
          "\\(S(x) = 1 - \\dfrac{x}{\\omega}\\), ",
          "\\(\\mu(x) = \\dfrac{1}{\\omega - x}\\)<br>",
          "• <strong>Exponencial</strong>: ",
          "\\(S(x) = e^{-\\mu x}\\), \\(\\mu(x) = \\mu\\)<br>",
          "• <strong>Cuadrático</strong>: ",
          "\\(S(x) = 1 - \\dfrac{x^2}{\\omega^2}\\), ",
          "\\(\\mu(x) = \\dfrac{2x}{\\omega^2 - x^2}\\)"
        ))
      )
    ),

    div(class = "enunciado",
      span(class = "contexto", "2. Edades de las dos personas"),
      fluidRow(
        column(6, numericInput("edad_x", "Edad x:", value = 40, min = 0, step = 1)),
        column(6, numericInput("edad_y", "Edad y:", value = 50, min = 0, step = 1))
      )
    ),

    div(class = "enunciado",
      span(class = "contexto", "3. ¿Qué esperanza quieres calcular?"),
      radioButtons("tipo", NULL,
        choices = c(
          "Hasta disolución (el primero que muere)  ·  \\(e_{xy}\\)" = "disolucion",
          "Hasta extinción (el último que muere)    ·  \\(e_{\\overline{xy}}\\)" = "extincion",
          "Ambas (con la identidad)" = "ambas"
        ),
        selected = "ambas", inline = FALSE)
    ),

    uiOutput("desarrollo"),

    div(class = "footer",
      HTML("Esperanza de vida conjunta · l₀ = 100"))
  )
)

# ============================================================
#  SERVER
# ============================================================
server <- function(input, output, session) {

  mod <- reactive({
    if (input$modelo == "exp") {
      validate(need(input$mu > 0, "μ debe ser positivo"))
      construir_modelo("exp", omega = 100, mu = input$mu)
    } else if (input$modelo == "moivre" || input$modelo == "cuad") {
      validate(need(input$omega > 1, "ω debe ser > 1"))
      construir_modelo(input$modelo, omega = input$omega, mu = 0.02)
    } else {
      construir_modelo("base", omega = 100, mu = 0.02)
    }
  })

  observeEvent(mod(), {
    m <- mod()
    if (!is.finite(m$omega)) return()
    if (!is.null(input$edad_x) && input$edad_x >= m$omega)
      updateNumericInput(session, "edad_x", value = max(0, floor(m$omega * 0.4)))
    if (!is.null(input$edad_y) && input$edad_y >= m$omega)
      updateNumericInput(session, "edad_y", value = max(0, floor(m$omega * 0.5)))
  }, ignoreInit = TRUE)

  observeEvent(input$generar, {
    m <- mod()
    if (is.finite(m$omega) && m$omega > 10) {
      max_edad <- max(1, floor(m$omega - 10))
      min_edad <- min(20, max_edad - 1)
      updateNumericInput(session, "edad_x", value = sample(min_edad:max_edad, 1))
      updateNumericInput(session, "edad_y", value = sample(min_edad:max_edad, 1))
    } else {
      updateNumericInput(session, "edad_x", value = sample(20:70, 1))
      updateNumericInput(session, "edad_y", value = sample(20:70, 1))
    }
  })

  output$desarrollo <- renderUI({
    x <- input$edad_x
    y <- input$edad_y
    m <- mod()
    req(x, y, x >= 0, y >= 0)
    if (is.finite(m$omega)) req(x < m$omega, y < m$omega)

    mu_val  <- if (!is.null(input$mu)) input$mu else 0.02

    S_fn  <- m$S
    mu_fn <- m$mu

    tmax <- if (is.finite(m$omega)) m$omega - max(x, y) else 80

    p_res  <- function(t, xx) S_fn(xx + t) / S_fn(xx)
    mu_res <- function(t, xx) mu_fn(xx + t)

    # ============================================================
    #  CÁLCULO DE ESPERANZAS
    # ============================================================
    if (input$modelo == "exp") {
      ex_ind  <- 1 / mu_val
      ey_ind  <- 1 / mu_val
      exy     <- 1 / (2 * mu_val)
    } else {
      ex_ind <- {
        tm <- if (is.finite(m$omega)) m$omega - x - 1e-9 else 80
        if (tm <= 0) NA_real_ else integrar(function(t) t * p_res(t, x) * mu_res(t, x), 0, tm, 5000)
      }
      ey_ind <- {
        tm <- if (is.finite(m$omega)) m$omega - y - 1e-9 else 80
        if (tm <= 0) NA_real_ else integrar(function(t) t * p_res(t, y) * mu_res(t, y), 0, tm, 5000)
      }
      exy <- {
        tm <- tmax - 1e-9
        if (tm <= 0) NA_real_ else {
          f <- function(t) {
            p_xy  <- p_res(t, x) * p_res(t, y)
            mu_xy <- mu_res(t, x) + mu_res(t, y)
            t * p_xy * mu_xy
          }
          integrar(f, 0, tm, 5000)
        }
      }
    }

    exy_ext <- ex_ind + ey_ind - exy

    px5  <- p_res(5, x);  py5  <- p_res(5, y)
    pxy5 <- px5 * py5
    mux5 <- mu_res(5, x); muy5 <- mu_res(5, y)

    # ============================================================
    #  PASOS 2, 3 y 4
    # ============================================================
    paso2_tex <- switch(input$modelo,
      base   = paste0(
        "$$ p_x(t) = \\frac{S(x+t)}{S(x)} = \\frac{1 - \\frac{x+t}{100}}{1 - \\frac{x}{100}} = 1 - \\frac{t}{", 100 - x, "} $$",
        "<br>$$ p_y(t) = 1 - \\frac{t}{", 100 - y, "} $$"),
      moivre = paste0(
        "$$ p_x(t) = 1 - \\frac{t}{", m$omega - x, "}, \\qquad p_y(t) = 1 - \\frac{t}{", m$omega - y, "} $$"),
      exp = paste0(
        "$$ p_x(t) = e^{-", sprintf("%.4f", mu_val), "t}, \\qquad p_y(t) = e^{-",
        sprintf("%.4f", mu_val), "t} $$"),
      cuad = paste0(
        "$$ p_x(t) = \\frac{1 - \\frac{(x+t)^2}{", m$omega, "^2}}{1 - \\frac{x^2}{", m$omega, "^2}}, \\qquad p_y(t) = \\frac{1 - \\frac{(y+t)^2}{", m$omega, "^2}}{1 - \\frac{y^2}{", m$omega, "^2}} $$")
    )

    paso3_tex <- switch(input$modelo,
      base   = paste0(
        "$$ \\mu_x(t) = \\frac{1}{", 100 - x, " - t}, \\qquad \\mu_y(t) = \\frac{1}{", 100 - y, " - t} $$"),
      moivre = paste0(
        "$$ \\mu_x(t) = \\frac{1}{", m$omega - x, " - t}, \\qquad \\mu_y(t) = \\frac{1}{", m$omega - y, " - t} $$"),
      exp = paste0(
        "$$ \\mu_x(t) = ", sprintf("%.4f", mu_val), ", \\qquad \\mu_y(t) = ", sprintf("%.4f", mu_val), " $$"),
      cuad = paste0(
        "$$ \\mu_x(t) = \\frac{2(x+t)}{", m$omega, "^2 - (x+t)^2}, \\qquad \\mu_y(t) = \\frac{2(y+t)}{", m$omega, "^2 - (y+t)^2} $$")
    )

    paso4_tex <- switch(input$modelo,
      base   = paste0(
        "$$ p_{xy}(t) = \\left(1 - \\frac{t}{", 100 - x, "}\\right)\\left(1 - \\frac{t}{", 100 - y, "}\\right) $$",
        "<br>$$ \\mu_{xy}(t) = \\frac{1}{", 100 - x, " - t} + \\frac{1}{", 100 - y, " - t} $$"),
      moivre = paste0(
        "$$ p_{xy}(t) = \\left(1 - \\frac{t}{", m$omega - x, "}\\right)\\left(1 - \\frac{t}{", m$omega - y, "}\\right) $$",
        "<br>$$ \\mu_{xy}(t) = \\frac{1}{", m$omega - x, " - t} + \\frac{1}{", m$omega - y, " - t} $$"),
      exp = paste0(
        "$$ p_{xy}(t) = e^{-", sprintf("%.4f", 2 * mu_val), "t} $$",
        "<br>$$ \\mu_{xy}(t) = ", sprintf("%.4f", 2 * mu_val), " $$"),
      cuad = paste0(
        "$$ p_{xy}(t) = p_x(t)\\,p_y(t), \\qquad \\mu_{xy}(t) = \\mu_x(t) + \\mu_y(t) $$")
    )

    # ============================================================
    #  PASO 5: ESPERANZAS INDIVIDUALES
    # ============================================================
    paso5_individuales <- switch(input$modelo,

      base = paste0(
        "La esperanza de vida individual de una persona de edad \\(z\\) es:",
        "$$ e_z = \\int_0^{100 - z} t \\cdot p_z(t)\\,\\mu_z(t)\\,dt $$",

        "<br><strong>(a) Para el modelo base</strong>, ",
        "\\(p_z(t) = 1 - \\dfrac{t}{100 - z}\\) y ",
        "\\(\\mu_z(t) = \\dfrac{1}{100 - z - t}\\). ",
        "Al multiplicarlas, el factor \\((100 - z - t)\\) se cancela:",
        "$$ p_z(t)\\,\\mu_z(t) = \\frac{100 - z - t}{100 - z}\\cdot\\frac{1}{100 - z - t} ",
        "= \\frac{1}{100 - z} $$",

        "<br><strong>(b) La integral es inmediata:</strong>",
        "$$ e_z = \\int_0^{100 - z} t \\cdot \\frac{1}{100 - z}\\,dt ",
        "= \\frac{1}{100 - z}\\cdot\\frac{t^2}{2}\\Bigg|_0^{100 - z} ",
        "= \\frac{100 - z}{2} $$",

        "<br><strong>(c) Aplicando a \\(x = ", x, "\\) e \\(y = ", y, "\\):</strong>",
        "$$ e_x = \\frac{100 - ", x, "}{2} = ", sprintf("%.6f", ex_ind), "\\;\\text{años}, \\qquad ",
        "e_y = \\frac{100 - ", y, "}{2} = ", sprintf("%.6f", ey_ind), "\\;\\text{años} $$"
      ),

      moivre = paste0(
        "La esperanza de vida individual de una persona de edad \\(z\\) es:",
        "$$ e_z = \\int_0^{\\omega - z} t \\cdot p_z(t)\\,\\mu_z(t)\\,dt $$",

        "<br><strong>(a) Para Moivre</strong> con \\(\\omega = ", m$omega, "\\), ",
        "\\(p_z(t) = 1 - \\dfrac{t}{\\omega - z}\\) y ",
        "\\(\\mu_z(t) = \\dfrac{1}{\\omega - z - t}\\). ",
        "Al multiplicarlas, el factor \\((\\omega - z - t)\\) se cancela:",
        "$$ p_z(t)\\,\\mu_z(t) = \\frac{\\omega - z - t}{\\omega - z}\\cdot\\frac{1}{\\omega - z - t} ",
        "= \\frac{1}{\\omega - z} $$",

        "<br><strong>(b) La integral es inmediata:</strong>",
        "$$ e_z = \\int_0^{\\omega - z} t \\cdot \\frac{1}{\\omega - z}\\,dt ",
        "= \\frac{1}{\\omega - z}\\cdot\\frac{t^2}{2}\\Bigg|_0^{\\omega - z} ",
        "= \\frac{\\omega - z}{2} $$",

        "<br><strong>(c) Aplicando a \\(x = ", x, "\\) e \\(y = ", y, "\\):</strong>",
        "$$ e_x = \\frac{", m$omega, " - ", x, "}{2} = ", sprintf("%.6f", ex_ind), "\\;\\text{años}, \\qquad ",
        "e_y = \\frac{", m$omega, " - ", y, "}{2} = ", sprintf("%.6f", ey_ind), "\\;\\text{años} $$"
      ),

      exp = {
        paste0(
          "La esperanza de vida individual de una persona de edad \\(z\\) es:",
          "$$ e_z = \\int_0^{\\infty} t \\cdot p_z(t)\\,\\mu_z(t)\\,dt $$",

          "<br><strong>(a) Para la exponencial</strong> con parámetro \\(\\mu = ",
          sprintf("%.4f", mu_val), "\\), ",
          "\\(p_z(t) = e^{-", sprintf("%.4f", mu_val), "t}\\) y ",
          "\\(\\mu_z(t) = ", sprintf("%.4f", mu_val), "\\). ",
          "Sustituyendo:",
          "$$ e_z = \\int_0^{\\infty} t \\cdot e^{-", sprintf("%.4f", mu_val), "t} \\cdot ",
          sprintf("%.4f", mu_val), "\\,dt ",
          "= ", sprintf("%.4f", mu_val), "\\int_0^{\\infty} t\\,e^{-",
          sprintf("%.4f", mu_val), "t}\\,dt $$",

          "<br><strong>(b) Integramos por partes</strong> con ",
          "\\(u = t\\), \\(dv = e^{-at}dt\\), \\(a = ", sprintf("%.4f", mu_val), "\\):<br>",
          "\\(du = dt\\), \\(v = -\\dfrac{1}{a}e^{-at}\\):",
          "$$ \\int_0^{\\infty} t\\,e^{-at}\\,dt ",
          "= \\underbrace{\\left[-\\frac{t}{a}e^{-at}\\right]_0^{\\infty}}_{=\\,0} ",
          "+ \\frac{1}{a}\\int_0^{\\infty} e^{-at}\\,dt $$",

          "<br><strong>(c) El corchete se anula</strong>: en \\(t\\to\\infty\\), ",
          "\\(t\\,e^{-at}\\to 0\\) (la exponencial decae más rápido que crece \\(t\\)); ",
          "en \\(t=0\\) vale \\(0\\). Queda:",
          "$$ \\int_0^{\\infty} t\\,e^{-at}\\,dt = \\frac{1}{a}\\int_0^{\\infty} e^{-at}\\,dt ",
          "= \\frac{1}{a}\\cdot\\frac{1}{a} = \\frac{1}{a^2} $$",

          "<br><strong>(d) Sustituimos</strong> \\(a = ", sprintf("%.4f", mu_val), "\\) ",
          "y multiplicamos por el factor \\(\\mu\\):",
          "$$ \\int_0^{\\infty} t\\,e^{-", sprintf("%.4f", mu_val), "t}\\,dt ",
          "= \\frac{1}{", sprintf("%.4f", mu_val), "^2} = ",
          sprintf("%.4f", 1/mu_val^2), " $$",
          "$$ e_z = ", sprintf("%.4f", mu_val), " \\cdot ",
          sprintf("%.4f", 1/mu_val^2), " = \\frac{1}{\\mu} ",
          "= ", sprintf("%.6f", 1/mu_val), "\\;\\text{años} $$",

          "<br><strong>(e) Para \\(x\\) e \\(y\\)</strong> el valor es el mismo ",
          "(falta de memoria de la exponencial):",
          "$$ e_x = e_y = ", sprintf("%.6f", ex_ind), "\\;\\text{años} $$"
        )
      },

      cuad = {
        w <- m$omega
        paste0(
          "La esperanza de vida individual de una persona de edad \\(z\\) es:",
          "$$ e_z = \\int_0^{\\omega - z} t \\cdot p_z(t)\\,\\mu_z(t)\\,dt $$",

          "<br><strong>(a) Para el modelo cuadrático</strong> con \\(\\omega = ", w, "\\), ",
          "\\(p_z(t) = \\dfrac{1 - \\frac{(z+t)^2}{", w, "^2}}{1 - \\frac{z^2}{", w, "^2}}\\) y ",
          "\\(\\mu_z(t) = \\dfrac{2(z+t)}{", w, "^2 - (z+t)^2}\\). ",
          "El factor \\(\\left(1 - \\frac{(z+t)^2}{", w, "^2}\\right)\\) del numerador de \\(p_z\\) ",
          "se cancela con el denominador de \\(\\mu_z\\):",
          "$$ p_z(t)\\,\\mu_z(t) = \\frac{2(z+t)}{", w, "^2 - z^2} $$",

          "<br><strong>(b) La integral</strong> es polinómica:",
          "$$ e_z = \\frac{2}{", w, "^2 - z^2}\\int_0^{", w, " - z} t\\,(z+t)\\,dt ",
          "= \\frac{2}{", w, "^2 - z^2}\\left(\\frac{z\\,(", w, " - z)^2}{2} + \\frac{(", w, " - z)^3}{3}\\right) $$",

          "<br><strong>(c) Aplicando a \\(x = ", x, "\\) e \\(y = ", y, "\\):</strong>",
          "$$ e_x = ", sprintf("%.6f", ex_ind), "\\;\\text{años}, \\qquad ",
          "e_y = ", sprintf("%.6f", ey_ind), "\\;\\text{años} $$"
        )
      }
    )

    # ============================================================
    #  PASO 6: DISOLUCIÓN
    # ============================================================
    paso6_disolucion <- switch(input$modelo,

      base = {
        dx <- 100 - x; dy <- 100 - y
        paste0(
          "$$ e_{", x, ",", y, "} = \\int_0^{", tmax, "} t \\cdot p_{xy}(t)\\,\\mu_{xy}(t)\\,dt $$",
          "<br><strong>(a) Sustituimos:</strong>",
          "$$ e_{", x, ",", y, "} = \\int_0^{", tmax, "} t \\cdot ",
          "\\left(1 - \\frac{t}{", dx, "}\\right)\\left(1 - \\frac{t}{", dy, "}\\right)",
          "\\left(\\frac{1}{", dx, " - t} + \\frac{1}{", dy, " - t}\\right)\\,dt $$",
          "<br><strong>(b) Simplificamos</strong> con \\(1 - \\frac{t}{", dx, "} = \\frac{", dx, " - t}{", dx, "}\\):",
          "$$ \\left(1 - \\frac{t}{", dx, "}\\right)\\frac{1}{", dx, " - t} = \\frac{1}{", dx, "}, \\qquad ",
          "\\left(1 - \\frac{t}{", dy, "}\\right)\\frac{1}{", dy, " - t} = \\frac{1}{", dy, "} $$",
          "<br>Integrando: \\(t \\cdot \\left(\\frac{1}{", dx, "} + \\frac{1}{", dy, "} - \\frac{2t}{", dx, "\\cdot", dy, "}\\right)\\)",
          "<br><strong>(c) Integramos:</strong>",
          "$$ e_{", x, ",", y, "} = \\left(\\frac{1}{", dx, "} + \\frac{1}{", dy, "}\\right)\\frac{t^2}{2} ",
          "- \\frac{2}{", dx, "\\cdot", dy, "}\\frac{t^3}{3} \\;\\Bigg|_0^{", tmax, "} $$",
          "<br><strong>(d) Evaluamos</strong> en \\(t = ", tmax, "\\):",
          "$$ e_{", x, ",", y, "} = \\left(\\frac{1}{", dx, "} + \\frac{1}{", dy, "}\\right)\\frac{", tmax, "^2}{2} ",
          "- \\frac{2}{", dx, "\\cdot", dy, "}\\frac{", tmax, "^3}{3} $$",
          "<br>$$ \\boxed{\\; e_{", x, ",", y, "} = ", sprintf("%.6f", exy), "\\;\\text{años} \\;} $$"
        )
      },

      moivre = {
        dx <- m$omega - x; dy <- m$omega - y
        paste0(
          "$$ e_{", x, ",", y, "} = \\int_0^{", tmax, "} t \\cdot p_{xy}(t)\\,\\mu_{xy}(t)\\,dt $$",
          "<br><strong>(a) Sustituimos:</strong>",
          "$$ e_{", x, ",", y, "} = \\int_0^{", tmax, "} t \\cdot ",
          "\\left(1 - \\frac{t}{", dx, "}\\right)\\left(1 - \\frac{t}{", dy, "}\\right)",
          "\\left(\\frac{1}{", dx, " - t} + \\frac{1}{", dy, " - t}\\right)\\,dt $$",
          "<br><strong>(b) Simplificamos:</strong>",
          "$$ \\left(1 - \\frac{t}{", dx, "}\\right)\\frac{1}{", dx, " - t} = \\frac{1}{", dx, "}, \\qquad ",
          "\\left(1 - \\frac{t}{", dy, "}\\right)\\frac{1}{", dy, " - t} = \\frac{1}{", dy, "} $$",
          "<br>Integrando: \\(t \\cdot \\left(\\frac{1}{", dx, "} + \\frac{1}{", dy, "} - \\frac{2t}{", dx, "\\cdot", dy, "}\\right)\\)",
          "<br><strong>(c) Integramos:</strong>",
          "$$ e_{", x, ",", y, "} = \\left(\\frac{1}{", dx, "} + \\frac{1}{", dy, "}\\right)\\frac{t^2}{2} ",
          "- \\frac{2}{", dx, "\\cdot", dy, "}\\frac{t^3}{3} \\;\\Bigg|_0^{", tmax, "} $$",
          "<br><strong>(d) Evaluamos</strong> en \\(t = ", tmax, "\\):",
          "$$ e_{", x, ",", y, "} = \\left(\\frac{1}{", dx, "} + \\frac{1}{", dy, "}\\right)\\frac{", tmax, "^2}{2} ",
          "- \\frac{2}{", dx, "\\cdot", dy, "}\\frac{", tmax, "^3}{3} $$",
          "<br>$$ \\boxed{\\; e_{", x, ",", y, "} = ", sprintf("%.6f", exy), "\\;\\text{años} \\;} $$"
        )
      },

      exp = {
        mu2 <- 2 * mu_val; a2 <- mu2^2; I <- 1/a2; prod <- mu2 * I
        paste0(
          "$$ e_{", x, ",", y, "} = \\int_0^{\\infty} t \\cdot p_{xy}(t)\\,\\mu_{xy}(t)\\,dt $$",
          "<br><strong>(a) Sustituimos</strong> \\(p_{xy} = e^{-", sprintf("%.4f", mu2), "t}\\), \\(\\mu_{xy} = ", sprintf("%.4f", mu2), "\\):",
          "$$ e_{", x, ",", y, "} = ", sprintf("%.4f", mu2), "\\int_0^{\\infty} t\\,e^{-", sprintf("%.4f", mu2), "t}\\,dt $$",
          "<br><strong>(b) Por partes</strong> con \\(u=t\\), \\(dv = e^{-at}dt\\), \\(a = ", sprintf("%.4f", mu2), "\\):",
          "$$ \\int_0^{\\infty} t\\,e^{-at}\\,dt = \\underbrace{\\left[-\\frac{t}{a}e^{-at}\\right]_0^{\\infty}}_{=\\,0} ",
          "+ \\frac{1}{a}\\int_0^{\\infty} e^{-at}\\,dt = \\frac{1}{a^2} $$",
          "<br><strong>(c) Sustituimos</strong> \\(a = ", sprintf("%.4f", mu2), "\\): \\(\\int_0^{\\infty} t\\,e^{-at}\\,dt = \\frac{1}{", sprintf("%.4f", a2), "} = ", sprintf("%.4f", I), "\\)",
          "<br><strong>(d) Multiplicamos:</strong> \\(e_{", x, ",", y, "} = ", sprintf("%.4f", mu2), "\\cdot", sprintf("%.4f", I), " = ", sprintf("%.4f", prod), "\\)",
          "<br>$$ \\boxed{\\; e_{", x, ",", y, "} = \\frac{1}{2\\mu} = ", sprintf("%.6f", exy), "\\;\\text{años} \\;} $$"
        )
      },

      cuad = {
        w <- m$omega
        paste0(
          "$$ e_{", x, ",", y, "} = \\int_0^{", tmax, "} t \\cdot p_{xy}(t)\\,\\mu_{xy}(t)\\,dt $$",
          "<br><strong>(a) Sustituimos:</strong>",
          "$$ p_x(t) = \\frac{1 - \\frac{(x+t)^2}{", w, "^2}}{1 - \\frac{x^2}{", w, "^2}}, \\quad \\mu_x(t) = \\frac{2(x+t)}{", w, "^2 - (x+t)^2} $$",
          "<br><strong>(b) Cancelamos</strong> \\(1 - \\frac{(x+t)^2}{", w, "^2}\\) con el denominador de \\(\\mu_x(t)\\):",
          "$$ p_x(t)\\,\\mu_x(t) = \\frac{2(x+t)}{", w, "^2 - x^2}, \\qquad p_y(t)\\,\\mu_y(t) = \\frac{2(y+t)}{", w, "^2 - y^2} $$",
          "<br><strong>(c) El integrando</strong> es un polinomio de grado 4 en \\(t\\).",
          "<br><strong>(d) Integramos</strong> numéricamente entre \\(0\\) y \\(", tmax, "\\):",
          "<br>$$ \\boxed{\\; e_{", x, ",", y, "} \\approx ", sprintf("%.6f", exy), "\\;\\text{años} \\;} $$"
        )
      }
    )

    # ============================================================
    #  PASO 7: EXTINCIÓN POR LA IDENTIDAD
    # ============================================================
    paso7_extincion <- paste0(
      "El grupo se extingue cuando muere el <strong>último</strong> de sus miembros: ",
      "\\(T_{\\overline{xy}} = \\max(T_x, T_y)\\).<br><br>",

      "<strong>Idea clave (unión e intersección de tiempos).</strong> ",
      "El tiempo total vivido por los dos es \\(T_x + T_y\\). ",
      "Ese total se reparte entre:<br>",
      "• el tiempo en que están <strong>los dos vivos</strong> (intersección) \\(= T_{xy} = \\min(T_x, T_y)\\),<br>",
      "• el tiempo en que está <strong>al menos uno vivo</strong> (unión) \\(= T_{\\overline{xy}} = \\max(T_x, T_y)\\).<br>",
      "Como el tiempo compartido se cuenta dos veces al sumar \\(T_x + T_y\\), hay que restarlo una vez:",

      "$$ T_{\\overline{xy}} = T_x + T_y - T_{xy} $$",

      "<br>Tomando esperanzas (linealidad):",

      "$$ e_{\\overline{xy}} = e_x + e_y - e_{xy} $$",

      "<br><strong>Aplicación numérica</strong> con los valores ya calculados:",

      "<br>\\(e_x = ", sprintf("%.6f", ex_ind), "\\), \\(",
      "e_y = ", sprintf("%.6f", ey_ind), "\\), \\(",
      "e_{xy} = ", sprintf("%.6f", exy), "\\)",

      "$$ e_{\\overline{", x, ",", y, "}} = ", sprintf("%.6f", ex_ind),
      " + ", sprintf("%.6f", ey_ind),
      " - ", sprintf("%.6f", exy),
      " = \\boxed{\\; ", sprintf("%.6f", exy_ext), "\\;\\text{años} \\;} $$",

      "<br><em>No hace falta una integral nueva: la extinción es la unión de las dos vidas, ",
      "y la unión es la suma menos la intersección.</em>"
    )

    # ============================================================
    #  TIPO DE CÁLCULO SELECCIONADO
    # ============================================================
    mostrar_disolucion <- input$tipo %in% c("disolucion", "ambas")
    mostrar_extincion  <- input$tipo %in% c("extincion", "ambas")
    mostrar_ambas      <- input$tipo == "ambas"

    # ============================================================
    #  UI FINAL
    # ============================================================
    tagList(
      div(class = "seccion fade-in",
        h3(HTML("Resolución paso a paso")),

        # PASO 1
        div(class = "paso",
          HTML("<strong>1. Funciones del modelo</strong>"),
          div(class = "teoria",
            HTML(paste0("Modelo: <strong>", m$nombre, "</strong>, con \\(l_0 = 100\\)."))),
          div(class = "deduccion",
            HTML(paste0("<br>", m$paso1_tex,
                        "<br>", if (is.finite(m$omega))
                          paste0("\\(\\omega =\\) ", m$omega, ".")
                        else
                          paste0("\\(\\mu =\\) ", sprintf("%.4f", mu_val), "."))))
        ),

        # PASO 2
        div(class = "paso",
          HTML("<strong>2. Supervivencias individuales residuales</strong>"),
          div(class = "teoria",
            HTML("\\(p_x(t) = \\dfrac{S(x+t)}{S(x)}\\), \\(p_y(t) = \\dfrac{S(y+t)}{S(y)}\\).")),
          div(class = "deduccion",
            HTML(paste0(
              paso2_tex,
              "<br><br>Ejemplo para \\(t=5\\): ",
              "\\(p_x(5) = ", sprintf("%.6f", px5), "\\), ",
              "\\(p_y(5) = ", sprintf("%.6f", py5), "\\)"
            )))
        ),

        # PASO 3
        div(class = "paso",
          HTML("<strong>3. Tanto instantáneo de mortalidad</strong>"),
          div(class = "teoria",
            HTML("Trasladado al tiempo \\(t\\):")),
          div(class = "deduccion",
            HTML(paste0(
              paso3_tex,
              "<br><br>Ejemplo para \\(t=5\\): ",
              "\\(\\mu_x(5) = ", sprintf("%.6f", mux5), "\\), ",
              "\\(\\mu_y(5) = ", sprintf("%.6f", muy5), "\\)"
            )))
        ),

        # PASO 4
        div(class = "paso",
          HTML("<strong>4. Supervivencia conjunta y tanto conjunto (disolución)</strong>"),
          div(class = "teoria",
            HTML("Para \\(T_{xy} = \\min(T_x, T_y)\\):")),
          div(class = "deduccion",
            HTML(paste0(
              paso4_tex,
              "<br><br>Ejemplo para \\(t=5\\): ",
              "\\(p_{xy}(5) = ", sprintf("%.6f", pxy5), "\\), ",
              "\\(\\mu_{xy}(5) = ", sprintf("%.6f", mux5 + muy5), "\\)",
              if (is.finite(m$omega))
                paste0("<br><br>Válido en \\(0 \\leq t \\leq ", tmax, "\\).")
              else "")))
        ),

        # PASO 5 - INDIVIDUALES
        div(class = "paso",
          HTML("<strong>5. Esperanzas de vida individuales · \\(e_x\\) y \\(e_y\\)</strong>"),
          div(class = "teoria",
            HTML("Antes de combinar, necesitamos cuánto vive por término medio cada persona por separado.")),
          div(class = "deduccion", HTML(paso5_individuales))
        ),

        # PASO 6 - DISOLUCIÓN
        if (mostrar_disolucion) {
          div(class = "paso",
            HTML("<strong>6. Esperanza de vida conjunta hasta disolución · \\(e_{xy}\\)</strong>"),
            div(class = "teoria",
              HTML(paste0(
                "Aplicamos \\(\\displaystyle e_{xy} = \\int_0^{t_{max}} t \\cdot p_{xy}(t)\\,\\mu_{xy}(t)\\,dt\\) con ",
                "\\(x = ", x, "\\), \\(y = ", y, "\\), \\(t_{max} = ", tmax, "\\):"
              ))),
            div(class = "deduccion", HTML(paso6_disolucion))
          )
        },

        # PASO 7 - EXTINCIÓN
        if (mostrar_extincion) {
          div(class = "paso",
            HTML("<strong>7. Esperanza de vida conjunta hasta extinción · \\(e_{\\overline{xy}}\\)</strong>"),
            div(class = "teoria",
              HTML("Ahora el grupo se extingue cuando muere el <strong>último</strong>. ",
                   "No integramos de nuevo: usamos la identidad con las esperanzas ya calculadas en los pasos 5 y 6.")),
            div(class = "deduccion", HTML(paso7_extincion))
          )
        },

        # RESULTADO
        div(class = "resultado-box",
          HTML(paste0(
            "<strong>Resultados:</strong><br>",
            if (mostrar_disolucion)
              paste0("$$ e_{", x, ",", y, "} = ", sprintf("%.6f", exy), "\\;\\text{años} \\quad \\text{(disolución)} $$"),
            if (mostrar_extincion)
              paste0("$$ e_{\\overline{", x, ",", y, "}} = ", sprintf("%.6f", exy_ext), "\\;\\text{años} \\quad \\text{(extinción)} $$"),
            if (mostrar_ambas) paste0(
              "<br><strong>Identidad:</strong> \\(e_{\\overline{xy}} = e_x + e_y - e_{xy}\\)",
              "<br>\\(e_x = ", sprintf("%.6f", ex_ind), "\\), ",
              "\\(e_y = ", sprintf("%.6f", ey_ind), "\\)",
              "<br>\\(e_x + e_y - e_{xy} = ", sprintf("%.6f", ex_ind + ey_ind - exy), "\\)",
              " = \\(e_{\\overline{xy}} = ", sprintf("%.6f", exy_ext), "\\)",
              "<br><em>La identidad se cumple por construcción: la extinción se ha calculado así.</em>"
            ),
            "<p>El grupo formado por una persona de <strong>", x,
            "</strong> años y otra de <strong>", y,
            "</strong> años:",
            if (mostrar_disolucion)
              paste0(" permanece unido una media de <strong>", sprintf("%.2f", exy), " años</strong> hasta la primera muerte"),
            if (mostrar_ambas) " ; ",
            if (mostrar_extincion)
              paste0(" se extingue (última muerte) tras una media de <strong>", sprintf("%.2f", exy_ext), " años</strong>"),
            ".</p>"
          ))
        )
      )
    )
  })
}

shinyApp(ui, server)