app.py

# --- INSTALACIÓN DE DEPENDENCIAS ADICIONALES ---
import os
import sys
import subprocess
os.system("pip install gradio==5.38.1")
import os
import io
import tempfile
import traceback
import zipfile
from typing import List, Tuple, Dict, Any, Optional, Union
from abc import ABC, abstractmethod
from unittest.mock import MagicMock
from dataclasses import dataclass
from enum import Enum
import json

from PIL import Image
import gradio as gr
import plotly.graph_objects as go
from plotly.subplots import make_subplots
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.integrate import odeint
from scipy.optimize import curve_fit, differential_evolution
from sklearn.metrics import mean_squared_error, r2_score
from docx import Document
from docx.shared import Inches
from fpdf import FPDF
from fpdf.enums import XPos, YPos

# --- SISTEMA DE INTERNACIONALIZACIÓN ---
class Language(Enum):
    ES = "Español"
    EN = "English"
    PT = "Português"
    FR = "Français"
    DE = "Deutsch"
    ZH = "中文"
    JA = "日本語"

TRANSLATIONS = {
    Language.ES: {
        "title": "🔬 Analizador de Cinéticas de Bioprocesos",
        "subtitle": "Análisis avanzado de modelos matemáticos biotecnológicos",
        "upload": "Sube tu archivo Excel (.xlsx)",
        "select_models": "Modelos a Probar",
        "analyze": "Analizar y Graficar",
        "results": "Resultados",
        "download": "Descargar",
        "biomass": "Biomasa",
        "substrate": "Sustrato",
        "product": "Producto",
        "time": "Tiempo",
        "parameters": "Parámetros",
        "model_comparison": "Comparación de Modelos",
        "dark_mode": "Modo Oscuro",
        "light_mode": "Modo Claro",
        "language": "Idioma",
        "theory": "Teoría y Modelos",
    },
    Language.EN: {
        "title": "🔬 Bioprocess Kinetics Analyzer",
        "subtitle": "Advanced analysis of biotechnological mathematical models",
        "upload": "Upload your Excel file (.xlsx)",
        "select_models": "Models to Test",
        "analyze": "Analyze and Plot",
        "results": "Results",
        "download": "Download",
        "biomass": "Biomass",
        "substrate": "Substrate",
        "product": "Product",
        "time": "Time",
        "parameters": "Parameters",
        "model_comparison": "Model Comparison",
        "dark_mode": "Dark Mode",
        "light_mode": "Light Mode",
        "language": "Language",
        "theory": "Theory and Models",
    },
}

# --- CONSTANTES MEJORADAS ---
C_TIME = 'tiempo'
C_BIOMASS = 'biomass'
C_SUBSTRATE = 'substrate'
C_PRODUCT = 'product'
COMPONENTS = [C_BIOMASS, C_SUBSTRATE, C_PRODUCT]

# --- SISTEMA DE TEMAS ---
THEMES = {
    "light": gr.themes.Soft(
        primary_hue="blue",
        secondary_hue="sky",
        neutral_hue="gray",
        font=[gr.themes.GoogleFont("Inter"), "ui-sans-serif", "sans-serif"]
    ),
    "dark": gr.themes.Base(
        primary_hue="blue",
        secondary_hue="cyan",
        neutral_hue="slate",
        font=[gr.themes.GoogleFont("Inter"), "ui-sans-serif", "sans-serif"]
    ).set(
        body_background_fill="*neutral_950",
        body_background_fill_dark="*neutral_950",
        button_primary_background_fill="*primary_600",
        button_primary_background_fill_hover="*primary_700",
    )
}

# --- MODELOS CINÉTICOS COMPLETOS ---

class KineticModel(ABC):
    def __init__(self, name: str, display_name: str, param_names: List[str],
                 description: str = "", equation: str = "", reference: str = ""):
        self.name = name
        self.display_name = display_name
        self.param_names = param_names
        self.num_params = len(param_names)
        self.description = description
        self.equation = equation
        self.reference = reference

    @abstractmethod
    def model_function(self, t: np.ndarray, *params: float) -> np.ndarray:
        pass

    def diff_function(self, X: float, t: float, params: List[float]) -> float:
        return 0.0

    @abstractmethod
    def get_initial_params(self, time: np.ndarray, biomass: np.ndarray) -> List[float]:
        pass

    @abstractmethod
    def get_param_bounds(self, time: np.ndarray, biomass: np.ndarray) -> Tuple[List[float], List[float]]:
        pass

# Modelo Logístico
class LogisticModel(KineticModel):
    def __init__(self):
        super().__init__(
            "logistic",
            "Logístico",
            ["X0", "Xm", "μm"],
            "Modelo de crecimiento logístico clásico para poblaciones limitadas",
            r"X(t) = \frac{X_0 X_m e^{\mu_m t}}{X_m - X_0 + X_0 e^{\mu_m t}}",
            "Verhulst (1838)"
        )

    def model_function(self, t: np.ndarray, *params: float) -> np.ndarray:
        X0, Xm, um = params
        if Xm <= 0 or X0 <= 0 or Xm < X0:
            return np.full_like(t, np.nan)
        exp_arg = np.clip(um * t, -700, 700)
        term_exp = np.exp(exp_arg)
        denominator = Xm - X0 + X0 * term_exp
        denominator = np.where(denominator == 0, 1e-9, denominator)
        return (X0 * term_exp * Xm) / denominator

    def diff_function(self, X: float, t: float, params: List[float]) -> float:
        _, Xm, um = params
        return um * X * (1 - X / Xm) if Xm > 0 else 0.0

    def get_initial_params(self, time: np.ndarray, biomass: np.ndarray) -> List[float]:
        return [
            biomass[0] if len(biomass) > 0 and biomass[0] > 1e-6 else 1e-3,
            max(biomass) if len(biomass) > 0 else 1.0,
            0.1
        ]

    def get_param_bounds(self, time: np.ndarray, biomass: np.ndarray) -> Tuple[List[float], List[float]]:
        initial_biomass = biomass[0] if len(biomass) > 0 else 1e-9
        max_biomass = max(biomass) if len(biomass) > 0 else 1.0
        return ([1e-9, initial_biomass, 1e-9], [max_biomass * 1.2, max_biomass * 5, np.inf])

# Modelo Gompertz
class GompertzModel(KineticModel):
    def __init__(self):
        super().__init__(
            "gompertz",
            "Gompertz",
            ["Xm", "μm", "λ"],
            "Modelo de crecimiento asimétrico con fase lag",
            r"X(t) = X_m \exp\left(-\exp\left(\frac{\mu_m e}{X_m}(\lambda-t)+1\right)\right)",
            "Gompertz (1825)"
        )

    def model_function(self, t: np.ndarray, *params: float) -> np.ndarray:
        Xm, um, lag = params
        if Xm <= 0 or um <= 0:
            return np.full_like(t, np.nan)
        exp_term = (um * np.e / Xm) * (lag - t) + 1
        exp_term_clipped = np.clip(exp_term, -700, 700)
        return Xm * np.exp(-np.exp(exp_term_clipped))

    def diff_function(self, X: float, t: float, params: List[float]) -> float:
        Xm, um, lag = params
        k_val = um * np.e / Xm
        u_val = k_val * (lag - t) + 1
        u_val_clipped = np.clip(u_val, -np.inf, 700)
        return X * k_val * np.exp(u_val_clipped) if Xm > 0 and X > 0 else 0.0

    def get_initial_params(self, time: np.ndarray, biomass: np.ndarray) -> List[float]:
        return [
            max(biomass) if len(biomass) > 0 else 1.0,
            0.1,
            time[np.argmax(np.gradient(biomass))] if len(biomass) > 1 else 0
        ]

    def get_param_bounds(self, time: np.ndarray, biomass: np.ndarray) -> Tuple[List[float], List[float]]:
        initial_biomass = min(biomass) if len(biomass) > 0 else 1e-9
        max_biomass = max(biomass) if len(biomass) > 0 else 1.0
        return ([max(1e-9, initial_biomass), 1e-9, 0], [max_biomass * 5, np.inf, max(time) if len(time) > 0 else 1])

# Modelo Moser
class MoserModel(KineticModel):
    def __init__(self):
        super().__init__(
            "moser",
            "Moser",
            ["Xm", "μm", "Ks"],
            "Modelo exponencial simple de Moser",
            r"X(t) = X_m (1 - e^{-\mu_m (t - K_s)})",
            "Moser (1958)"
        )

    def model_function(self, t: np.ndarray, *params: float) -> np.ndarray:
        Xm, um, Ks = params
        return Xm * (1 - np.exp(-um * (t - Ks))) if Xm > 0 and um > 0 else np.full_like(t, np.nan)

    def diff_function(self, X: float, t: float, params: List[float]) -> float:
        Xm, um, _ = params
        return um * (Xm - X) if Xm > 0 else 0.0

    def get_initial_params(self, time: np.ndarray, biomass: np.ndarray) -> List[float]:
        return [max(biomass) if len(biomass) > 0 else 1.0, 0.1, 0]

    def get_param_bounds(self, time: np.ndarray, biomass: np.ndarray) -> Tuple[List[float], List[float]]:
        initial_biomass = min(biomass) if len(biomass) > 0 else 1e-9
        max_biomass = max(biomass) if len(biomass) > 0 else 1.0
        return ([max(1e-9, initial_biomass), 1e-9, -np.inf], [max_biomass * 5, np.inf, np.inf])

# Modelo Baranyi
class BaranyiModel(KineticModel):
    def __init__(self):
        super().__init__(
            "baranyi",
            "Baranyi",
            ["X0", "Xm", "μm", "λ"],
            "Modelo de Baranyi con fase lag explícita",
            r"X(t) = X_m / [1 + ((X_m/X_0) - 1) \exp(-\mu_m A(t))]",
            "Baranyi & Roberts (1994)"
        )

    def model_function(self, t: np.ndarray, *params: float) -> np.ndarray:
        X0, Xm, um, lag = params
        if X0 <= 0 or Xm <= X0 or um <= 0 or lag < 0:
            return np.full_like(t, np.nan)
        A_t = t + (1 / um) * np.log(np.exp(-um * t) + np.exp(-um * lag) - np.exp(-um * (t + lag)))
        exp_um_At = np.exp(np.clip(um * A_t, -700, 700))
        numerator = Xm
        denominator = 1 + ((Xm / X0) - 1) * (1 / exp_um_At)
        return numerator / np.where(denominator == 0, 1e-9, denominator)

    def get_initial_params(self, time: np.ndarray, biomass: np.ndarray) -> List[float]:
        return [
            biomass[0] if len(biomass) > 0 and biomass[0] > 1e-6 else 1e-3,
            max(biomass) if len(biomass) > 0 else 1.0,
            0.1,
            time[np.argmax(np.gradient(biomass))] if len(biomass) > 1 else 0.0
        ]

    def get_param_bounds(self, time: np.ndarray, biomass: np.ndarray) -> Tuple[List[float], List[float]]:
        initial_biomass = biomass[0] if len(biomass) > 0 else 1e-9
        max_biomass = max(biomass) if len(biomass) > 0 else 1.0
        return ([1e-9, max(1e-9, initial_biomass), 1e-9, 0], [max_biomass * 1.2, max_biomass * 10, np.inf, max(time) if len(time) > 0 else 1])

# Modelo Monod
class MonodModel(KineticModel):
    def __init__(self):
        super().__init__(
            "monod",
            "Monod",
            ["μmax", "Ks", "Y", "m"],
            "Modelo de Monod con mantenimiento celular",
            r"\mu = \frac{\mu_{max} \cdot S}{K_s + S} - m",
            "Monod (1949)"
        )

    def model_function(self, t: np.ndarray, *params: float) -> np.ndarray:
        return np.full_like(t, np.nan)

    def diff_function(self, X: float, t: float, params: List[float]) -> float:
        μmax, Ks, Y, m = params
        S = 10.0
        μ = (μmax * S / (Ks + S)) - m
        return μ * X

    def get_initial_params(self, time: np.ndarray, biomass: np.ndarray) -> List[float]:
        return [0.5, 0.1, 0.5, 0.01]

    def get_param_bounds(self, time: np.ndarray, biomass: np.ndarray) -> Tuple[List[float], List[float]]:
        return ([0.01, 0.001, 0.1, 0.0], [2.0, 5.0, 1.0, 0.1])

# Modelo Contois
class ContoisModel(KineticModel):
    def __init__(self):
        super().__init__(
            "contois",
            "Contois",
            ["μmax", "Ksx", "Y", "m"],
            "Modelo de Contois para alta densidad celular",
            r"\mu = \frac{\mu_{max} \cdot S}{K_{sx} \cdot X + S} - m",
            "Contois (1959)"
        )

    def model_function(self, t: np.ndarray, *params: float) -> np.ndarray:
        return np.full_like(t, np.nan)

    def diff_function(self, X: float, t: float, params: List[float]) -> float:
        μmax, Ksx, Y, m = params
        S = 10.0
        μ = (μmax * S / (Ksx * X + S)) - m
        return μ * X

    def get_initial_params(self, time: np.ndarray, biomass: np.ndarray) -> List[float]:
        return [0.5, 0.5, 0.5, 0.01]

    def get_param_bounds(self, time: np.ndarray, biomass: np.ndarray) -> Tuple[List[float], List[float]]:
        return ([0.01, 0.01, 0.1, 0.0], [2.0, 10.0, 1.0, 0.1])

# Modelo Andrews
class AndrewsModel(KineticModel):
    def __init__(self):
        super().__init__(
            "andrews",
            "Andrews (Haldane)",
            ["μmax", "Ks", "Ki", "Y", "m"],
            "Modelo de inhibición por sustrato",
            r"\mu = \frac{\mu_{max} \cdot S}{K_s + S + \frac{S^2}{K_i}} - m",
            "Andrews (1968)"
        )

    def model_function(self, t: np.ndarray, *params: float) -> np.ndarray:
        return np.full_like(t, np.nan)

    def diff_function(self, X: float, t: float, params: List[float]) -> float:
        μmax, Ks, Ki, Y, m = params
        S = 10.0
        μ = (μmax * S / (Ks + S + S**2/Ki)) - m
        return μ * X

    def get_initial_params(self, time: np.ndarray, biomass: np.ndarray) -> List[float]:
        return [0.5, 0.1, 50.0, 0.5, 0.01]

    def get_param_bounds(self, time: np.ndarray, biomass: np.ndarray) -> Tuple[List[float], List[float]]:
        return ([0.01, 0.001, 1.0, 0.1, 0.0], [2.0, 5.0, 200.0, 1.0, 0.1])

# Modelo Tessier
class TessierModel(KineticModel):
    def __init__(self):
        super().__init__(
            "tessier",
            "Tessier",
            ["μmax", "Ks", "X0"],
            "Modelo exponencial de Tessier",
            r"\mu = \mu_{max} \cdot (1 - e^{-S/K_s})",
            "Tessier (1942)"
        )

    def model_function(self, t: np.ndarray, *params: float) -> np.ndarray:
        μmax, Ks, X0 = params
        return X0 * np.exp(μmax * t * 0.5)

    def diff_function(self, X: float, t: float, params: List[float]) -> float:
        μmax, Ks, X0 = params
        return μmax * X * 0.5

    def get_initial_params(self, time: np.ndarray, biomass: np.ndarray) -> List[float]:
        return [0.5, 1.0, biomass[0] if len(biomass) > 0 else 0.1]

    def get_param_bounds(self, time: np.ndarray, biomass: np.ndarray) -> Tuple[List[float], List[float]]:
        return ([0.01, 0.1, 1e-9], [2.0, 10.0, 1.0])

# Modelo Richards
class RichardsModel(KineticModel):
    def __init__(self):
        super().__init__(
            "richards",
            "Richards",
            ["A", "μm", "λ", "ν", "X0"],
            "Modelo generalizado de Richards",
            r"X(t) = A \cdot [1 + \nu \cdot e^{-\mu_m(t-\lambda)}]^{-1/\nu}",
            "Richards (1959)"
        )

    def model_function(self, t: np.ndarray, *params: float) -> np.ndarray:
        A, μm, λ, ν, X0 = params
        if A <= 0 or μm <= 0 or ν <= 0:
            return np.full_like(t, np.nan)
        exp_term = np.exp(-μm * (t - λ))
        return A * (1 + ν * exp_term) ** (-1/ν)

    def get_initial_params(self, time: np.ndarray, biomass: np.ndarray) -> List[float]:
        return [
            max(biomass) if len(biomass) > 0 else 1.0,
            0.5,
            time[len(time)//4] if len(time) > 0 else 1.0,
            1.0,
            biomass[0] if len(biomass) > 0 else 0.1
        ]

    def get_param_bounds(self, time: np.ndarray, biomass: np.ndarray) -> Tuple[List[float], List[float]]:
        max_biomass = max(biomass) if len(biomass) > 0 else 10.0
        max_time = max(time) if len(time) > 0 else 100.0
        return (
            [0.1, 0.01, 0.0, 0.1, 1e-9],
            [max_biomass * 2, 5.0, max_time, 10.0, max_biomass]
        )

# Modelo Stannard
class StannardModel(KineticModel):
    def __init__(self):
        super().__init__(
            "stannard",
            "Stannard",
            ["Xm", "μm", "λ", "α"],
            "Modelo de Stannard modificado",
            r"X(t) = X_m \cdot [1 - e^{-\mu_m(t-\lambda)^\alpha}]",
            "Stannard et al. (1985)"
        )

    def model_function(self, t: np.ndarray, *params: float) -> np.ndarray:
        Xm, μm, λ, α = params
        if Xm <= 0 or μm <= 0 or α <= 0:
            return np.full_like(t, np.nan)
        t_shifted = np.maximum(t - λ, 0)
        return Xm * (1 - np.exp(-μm * t_shifted ** α))

    def get_initial_params(self, time: np.ndarray, biomass: np.ndarray) -> List[float]:
        return [
            max(biomass) if len(biomass) > 0 else 1.0,
            0.5,
            0.0,
            1.0
        ]

    def get_param_bounds(self, time: np.ndarray, biomass: np.ndarray) -> Tuple[List[float], List[float]]:
        max_biomass = max(biomass) if len(biomass) > 0 else 10.0
        max_time = max(time) if len(time) > 0 else 100.0
        return ([0.1, 0.01, -max_time/10, 0.1], [max_biomass * 2, 5.0, max_time/2, 3.0])

# Modelo Huang
class HuangModel(KineticModel):
    def __init__(self):
        super().__init__(
            "huang",
            "Huang",
            ["Xm", "μm", "λ", "n", "m"],
            "Modelo de Huang para fase lag variable",
            r"X(t) = X_m \cdot \frac{1}{1 + e^{-\mu_m(t-\lambda-m/n)}}",
            "Huang (2008)"
        )

    def model_function(self, t: np.ndarray, *params: float) -> np.ndarray:
        Xm, μm, λ, n, m = params
        if Xm <= 0 or μm <= 0 or n <= 0:
            return np.full_like(t, np.nan)
        return Xm / (1 + np.exp(-μm * (t - λ - m/n)))

    def get_initial_params(self, time: np.ndarray, biomass: np.ndarray) -> List[float]:
        return [
            max(biomass) if len(biomass) > 0 else 1.0,
            0.5,
            time[len(time)//4] if len(time) > 0 else 1.0,
            1.0,
            0.5
        ]

    def get_param_bounds(self, time: np.ndarray, biomass: np.ndarray) -> Tuple[List[float], List[float]]:
        max_biomass = max(biomass) if len(biomass) > 0 else 10.0
        max_time = max(time) if len(time) > 0 else 100.0
        return (
            [0.1, 0.01, 0.0, 0.1, 0.0],
            [max_biomass * 2, 5.0, max_time/2, 10.0, 5.0]
        )

# --- REGISTRO ACTUALIZADO DE MODELOS ---
AVAILABLE_MODELS: Dict[str, KineticModel] = {
    model.name: model for model in [
        LogisticModel(),
        GompertzModel(),
        MoserModel(),
        BaranyiModel(),
        MonodModel(),
        ContoisModel(),
        AndrewsModel(),
        TessierModel(),
        RichardsModel(),
        StannardModel(),
        HuangModel()
    ]
}

# --- CLASE MEJORADA DE AJUSTE ---
class BioprocessFitter:
    def __init__(self, kinetic_model: KineticModel, maxfev: int = 50000,
                 use_differential_evolution: bool = False):
        self.model = kinetic_model
        self.maxfev = maxfev
        self.use_differential_evolution = use_differential_evolution
        self.params: Dict[str, Dict[str, float]] = {c: {} for c in COMPONENTS}
        self.r2: Dict[str, float] = {}
        self.rmse: Dict[str, float] = {}
        self.mae: Dict[str, float] = {}
        self.aic: Dict[str, float] = {}
        self.bic: Dict[str, float] = {}
        self.data_time: Optional[np.ndarray] = None
        self.data_means: Dict[str, Optional[np.ndarray]] = {c: None for c in COMPONENTS}
        self.data_stds: Dict[str, Optional[np.ndarray]] = {c: None for c in COMPONENTS}

    def _get_biomass_at_t(self, t: np.ndarray, p: List[float]) -> np.ndarray:
        return self.model.model_function(t, *p)

    def _get_initial_biomass(self, p: List[float]) -> float:
        if not p: return 0.0
        if any(k in self.model.param_names for k in ["Xo", "X0"]):
            try:
                idx = self.model.param_names.index("Xo") if "Xo" in self.model.param_names else self.model.param_names.index("X0")
                return p[idx]
            except (ValueError, IndexError): pass
        return float(self.model.model_function(np.array([0]), *p)[0])

    def _calc_integral(self, t: np.ndarray, p: List[float]) -> Tuple[np.ndarray, np.ndarray]:
        X_t = self._get_biomass_at_t(t, p)
        if np.any(np.isnan(X_t)): return np.full_like(t, np.nan), np.full_like(t, np.nan)
        integral_X = np.zeros_like(X_t)
        if len(t) > 1:
            dt = np.diff(t, prepend=t[0] - (t[1] - t[0] if len(t) > 1 else 1))
            integral_X = np.cumsum(X_t * dt)
        return integral_X, X_t

    def substrate(self, t: np.ndarray, so: float, p_c: float, q: float, bio_p: List[float]) -> np.ndarray:
        integral, X_t = self._calc_integral(t, bio_p)
        X0 = self._get_initial_biomass(bio_p)
        return so - p_c * (X_t - X0) - q * integral

    def product(self, t: np.ndarray, po: float, alpha: float, beta: float, bio_p: List[float]) -> np.ndarray:
        integral, X_t = self._calc_integral(t, bio_p)
        X0 = self._get_initial_biomass(bio_p)
        return po + alpha * (X_t - X0) + beta * integral

    def process_data_from_df(self, df: pd.DataFrame) -> None:
        try:
            time_col = [c for c in df.columns if c[1].strip().lower() == C_TIME][0]
            self.data_time = df[time_col].dropna().to_numpy()
            min_len = len(self.data_time)

            def extract(name: str) -> Tuple[np.ndarray, np.ndarray]:
                cols = [c for c in df.columns if c[1].strip().lower() == name.lower()]
                if not cols: return np.array([]), np.array([])
                reps = [df[c].dropna().values[:min_len] for c in cols]
                reps = [r for r in reps if len(r) == min_len]
                if not reps: return np.array([]), np.array([])
                arr = np.array(reps)
                mean = np.mean(arr, axis=0)
                std = np.std(arr, axis=0, ddof=1) if arr.shape[0] > 1 else np.zeros_like(mean)
                return mean, std

            self.data_means[C_BIOMASS], self.data_stds[C_BIOMASS] = extract('Biomasa')
            self.data_means[C_SUBSTRATE], self.data_stds[C_SUBSTRATE] = extract('Sustrato')
            self.data_means[C_PRODUCT], self.data_stds[C_PRODUCT] = extract('Producto')
        except (IndexError, KeyError) as e:
            raise ValueError(f"Estructura de DataFrame inválida. Error: {e}")

    def _calculate_metrics(self, y_true: np.ndarray, y_pred: np.ndarray,
                          n_params: int) -> Dict[str, float]:
        n = len(y_true)
        residuals = y_true - y_pred
        ss_res = np.sum(residuals**2)
        ss_tot = np.sum((y_true - np.mean(y_true))**2)
        r2 = 1 - (ss_res / ss_tot) if ss_tot > 0 else 0
        rmse = np.sqrt(ss_res / n)
        mae = np.mean(np.abs(residuals))
        if n > n_params + 1:
            aic = n * np.log(ss_res/n) + 2 * n_params
            bic = n * np.log(ss_res/n) + n_params * np.log(n)
        else:
            aic = bic = np.inf
        return {'r2': r2, 'rmse': rmse, 'mae': mae, 'aic': aic, 'bic': bic}

    def _fit_component_de(self, func, t, data, bounds, *args):
        def objective(params):
            try:
                pred = func(t, *params, *args)
                if np.any(np.isnan(pred)): return 1e10
                return np.sum((data - pred)**2)
            except:
                return 1e10
        result = differential_evolution(objective, bounds=list(zip(*bounds)), maxiter=1000, seed=42)
        if result.success:
            popt = result.x
            pred = func(t, *popt, *args)
            metrics = self._calculate_metrics(data, pred, len(popt))
            return list(popt), metrics
        return None, {'r2': np.nan, 'rmse': np.nan, 'mae': np.nan, 'aic': np.nan, 'bic': np.nan}

    def _fit_component(self, func, t, data, p0, bounds, sigma=None, *args):
        try:
            if self.use_differential_evolution:
                return self._fit_component_de(func, t, data, bounds, *args)
            if sigma is not None:
                sigma = np.where(sigma == 0, 1e-9, sigma)
            popt, _ = curve_fit(func, t, data, p0, bounds=bounds, maxfev=self.maxfev, ftol=1e-9, xtol=1e-9, sigma=sigma, absolute_sigma=bool(sigma is not None))
            pred = func(t, *popt, *args)
            if np.any(np.isnan(pred)):
                return None, {'r2': np.nan, 'rmse': np.nan, 'mae': np.nan, 'aic': np.nan, 'bic': np.nan}
            metrics = self._calculate_metrics(data, pred, len(popt))
            return list(popt), metrics
        except (RuntimeError, ValueError):
            return None, {'r2': np.nan, 'rmse': np.nan, 'mae': np.nan, 'aic': np.nan, 'bic': np.nan}

    def fit_all_models(self) -> None:
        t, bio_m, bio_s = self.data_time, self.data_means[C_BIOMASS], self.data_stds[C_BIOMASS]
        if t is None or bio_m is None or len(bio_m) == 0: return
        popt_bio = self._fit_biomass_model(t, bio_m, bio_s)
        if popt_bio:
            bio_p = list(self.params[C_BIOMASS].values())
            if self.data_means[C_SUBSTRATE] is not None and len(self.data_means[C_SUBSTRATE]) > 0:
                self._fit_substrate_model(t, self.data_means[C_SUBSTRATE], self.data_stds[C_SUBSTRATE], bio_p)
            if self.data_means[C_PRODUCT] is not None and len(self.data_means[C_PRODUCT]) > 0:
                self._fit_product_model(t, self.data_means[C_PRODUCT], self.data_stds[C_PRODUCT], bio_p)

    def _fit_biomass_model(self, t, data, std):
        p0, bounds = self.model.get_initial_params(t, data), self.model.get_param_bounds(t, data)
        popt, metrics = self._fit_component(self.model.model_function, t, data, p0, bounds, std)
        if popt:
            self.params[C_BIOMASS] = dict(zip(self.model.param_names, popt))
            self.r2[C_BIOMASS], self.rmse[C_BIOMASS], self.mae[C_BIOMASS], self.aic[C_BIOMASS], self.bic[C_BIOMASS] = metrics['r2'], metrics['rmse'], metrics['mae'], metrics['aic'], metrics['bic']
        return popt

    def _fit_substrate_model(self, t, data, std, bio_p):
        p0, b = [data[0], 0.1, 0.01], ([0, -np.inf, -np.inf], [np.inf, np.inf, np.inf])
        popt, metrics = self._fit_component(lambda t, so, p, q: self.substrate(t, so, p, q, bio_p), t, data, p0, b, std)
        if popt:
            self.params[C_SUBSTRATE] = {'So': popt[0], 'p': popt[1], 'q': popt[2]}
            self.r2[C_SUBSTRATE], self.rmse[C_SUBSTRATE], self.mae[C_SUBSTRATE], self.aic[C_SUBSTRATE], self.bic[C_SUBSTRATE] = metrics['r2'], metrics['rmse'], metrics['mae'], metrics['aic'], metrics['bic']

    def _fit_product_model(self, t, data, std, bio_p):
        p0, b = [data[0] if len(data)>0 else 0, 0.1, 0.01], ([0, -np.inf, -np.inf], [np.inf, np.inf, np.inf])
        popt, metrics = self._fit_component(lambda t, po, a, b: self.product(t, po, a, b, bio_p), t, data, p0, b, std)
        if popt:
            self.params[C_PRODUCT] = {'Po': popt[0], 'alpha': popt[1], 'beta': popt[2]}
            self.r2[C_PRODUCT], self.rmse[C_PRODUCT], self.mae[C_PRODUCT], self.aic[C_PRODUCT], self.bic[C_PRODUCT] = metrics['r2'], metrics['rmse'], metrics['mae'], metrics['aic'], metrics['bic']

    def system_ode(self, y, t, bio_p, sub_p, prod_p):
        X, _, _ = y
        dXdt = self.model.diff_function(X, t, bio_p)
        return [dXdt, -sub_p.get('p',0)*dXdt - sub_p.get('q',0)*X, prod_p.get('alpha',0)*dXdt + prod_p.get('beta',0)*X]

    def solve_odes(self, t_fine):
        p = self.params
        bio_d, sub_d, prod_d = p[C_BIOMASS], p[C_SUBSTRATE], p[C_PRODUCT]
        if not bio_d: return None, None, None
        try:
            bio_p = list(bio_d.values())
            y0 = [self._get_initial_biomass(bio_p), sub_d.get('So',0), prod_d.get('Po',0)]
            sol = odeint(self.system_ode, y0, t_fine, args=(bio_p, sub_d, prod_d))
            return sol[:, 0], sol[:, 1], sol[:, 2]
        except:
            return None, None, None

    def _generate_fine_time_grid(self, t_exp):
        return np.linspace(min(t_exp), max(t_exp), 500) if t_exp is not None and len(t_exp) > 1 else np.array([])

    def get_model_curves_for_plot(self, t_fine, use_diff):
        if use_diff and self.model.diff_function(1, 1, [1]*self.model.num_params) != 0:
            return self.solve_odes(t_fine)
        X, S, P = None, None, None
        if self.params[C_BIOMASS]:
            bio_p = list(self.params[C_BIOMASS].values())
            X = self.model.model_function(t_fine, *bio_p)
            if self.params[C_SUBSTRATE]:
                S = self.substrate(t_fine, *list(self.params[C_SUBSTRATE].values()), bio_p)
            if self.params[C_PRODUCT]:
                P = self.product(t_fine, *list(self.params[C_PRODUCT].values()), bio_p)
        return X, S, P

# --- FUNCIONES AUXILIARES ---
def format_number(value: Any, decimals: int) -> str:
    if not isinstance(value, (int, float, np.number)) or pd.isna(value):
        return "" if pd.isna(value) else str(value)
    decimals = int(decimals)
    if decimals == 0:
        if 0 < abs(value) < 1:
            return f"{value:.2e}"
        else:
            return str(int(round(value, 0)))
    return str(round(value, decimals))

# --- FUNCIONES DE PLOTEO MEJORADAS CON PLOTLY ---
def create_interactive_plot(plot_config: Dict, models_results: List[Dict],
                          selected_component: str = "all") -> go.Figure:
    time_exp = plot_config['time_exp']
    time_fine = np.linspace(min(time_exp), max(time_exp), 500)
    if selected_component == "all":
        fig = make_subplots(rows=3, cols=1, subplot_titles=('Biomasa', 'Sustrato', 'Producto'), vertical_spacing=0.08, shared_xaxes=True)
        components_to_plot, rows = [C_BIOMASS, C_SUBSTRATE, C_PRODUCT], [1, 2, 3]
    else:
        fig, components_to_plot, rows = go.Figure(), [selected_component], [None]
    colors = ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728', '#9467bd', '#8c564b', '#e377c2', '#7f7f7f', '#bcbd22', '#17becf']
    for comp, row in zip(components_to_plot, rows):
        data_exp, data_std = plot_config.get(f'{comp}_exp'), plot_config.get(f'{comp}_std')
        if data_exp is not None:
            error_y = dict(type='data', array=data_std, visible=True) if data_std is not None and np.any(data_std > 0) else None
            trace = go.Scatter(x=time_exp, y=data_exp, mode='markers', name=f'{comp.capitalize()} (Experimental)', marker=dict(size=10, symbol='circle'), error_y=error_y, legendgroup=comp, showlegend=True)
            if selected_component == "all": fig.add_trace(trace, row=row, col=1)
            else: fig.add_trace(trace)
    for i, res in enumerate(models_results):
        color, model_name = colors[i % len(colors)], AVAILABLE_MODELS[res["name"]].display_name
        for comp, row, key in zip(components_to_plot, rows, ['X', 'S', 'P']):
            if res.get(key) is not None:
                trace = go.Scatter(x=time_fine, y=res[key], mode='lines', name=f'{model_name} - {comp.capitalize()}', line=dict(color=color, width=2), legendgroup=f'{res["name"]}_{comp}', showlegend=True)
                if selected_component == "all": fig.add_trace(trace, row=row, col=1)
                else: fig.add_trace(trace)
    theme, template = plot_config.get('theme', 'light'), "plotly_white" if plot_config.get('theme', 'light') == 'light' else "plotly_dark"
    fig.update_layout(title=f"Análisis de Cinéticas: {plot_config.get('exp_name', '')}", template=template, hovermode='x unified', legend=dict(orientation="v", yanchor="middle", y=0.5, xanchor="left", x=1.02), margin=dict(l=80, r=250, t=100, b=80))
    if selected_component == "all":
        fig.update_xaxes(title_text="Tiempo", row=3, col=1)
        fig.update_yaxes(title_text="Biomasa (g/L)", row=1, col=1)
        fig.update_yaxes(title_text="Sustrato (g/L)", row=2, col=1)
        fig.update_yaxes(title_text="Producto (g/L)", row=3, col=1)
    else:
        fig.update_xaxes(title_text="Tiempo")
        labels = {C_BIOMASS: "Biomasa (g/L)", C_SUBSTRATE: "Sustrato (g/L)", C_PRODUCT: "Producto (g/L)"}
        fig.update_yaxes(title_text=labels.get(selected_component, "Valor"))
    return fig

# --- FUNCIÓN PRINCIPAL DE ANÁLISIS ---
def run_analysis(file, model_names, component, use_de, maxfev, exp_names, theme='light'):
    if not file: return None, pd.DataFrame(), "Error: Sube un archivo Excel."
    if not model_names: return None, pd.DataFrame(), "Error: Selecciona un modelo."
    try:
        xls = pd.ExcelFile(file.name)
    except Exception as e:
        return None, pd.DataFrame(), f"Error al leer archivo: {e}"
    results_data, msgs, models_results = [], [], []
    exp_list = [n.strip() for n in exp_names.split('\n') if n.strip()] if exp_names else []
    for i, sheet in enumerate(xls.sheet_names):
        exp_name = exp_list[i] if i < len(exp_list) else f"Hoja '{sheet}'"
        try:
            df = pd.read_excel(xls, sheet_name=sheet, header=[0,1])
            reader = BioprocessFitter(list(AVAILABLE_MODELS.values())[0])
            reader.process_data_from_df(df)
            if reader.data_time is None:
                msgs.append(f"WARN: Sin datos de tiempo en '{sheet}'.")
                continue
            plot_config = {'exp_name': exp_name, 'time_exp': reader.data_time, 'theme': theme}
            for c in COMPONENTS:
                plot_config[f'{c}_exp'], plot_config[f'{c}_std'] = reader.data_means[c], reader.data_stds[c]
            t_fine = reader._generate_fine_time_grid(reader.data_time)
            for m_name in model_names:
                if m_name not in AVAILABLE_MODELS:
                    msgs.append(f"WARN: Modelo '{m_name}' no disponible.")
                    continue
                fitter = BioprocessFitter(AVAILABLE_MODELS[m_name], maxfev=int(maxfev), use_differential_evolution=use_de)
                fitter.data_time, fitter.data_means, fitter.data_stds = reader.data_time, reader.data_means, reader.data_stds
                fitter.fit_all_models()
                row = {'Experimento': exp_name, 'Modelo': fitter.model.display_name}
                for c in COMPONENTS:
                    if fitter.params[c]:
                        row.update({f'{c.capitalize()}_{k}': v for k, v in fitter.params[c].items()})
                    row[f'R2_{c.capitalize()}'], row[f'RMSE_{c.capitalize()}'], row[f'MAE_{c.capitalize()}'], row[f'AIC_{c.capitalize()}'], row[f'BIC_{c.capitalize()}'] = fitter.r2.get(c), fitter.rmse.get(c), fitter.mae.get(c), fitter.aic.get(c), fitter.bic.get(c)
                results_data.append(row)
                X, S, P = fitter.get_model_curves_for_plot(t_fine, False)
                models_results.append({'name': m_name, 'X': X, 'S': S, 'P': P, 'params': fitter.params, 'r2': fitter.r2, 'rmse': fitter.rmse})
        except Exception as e:
            msgs.append(f"ERROR en '{sheet}': {e}")
            traceback.print_exc()
    msg = "Análisis completado." + ("\n" + "\n".join(msgs) if msgs else "")
    df_res = pd.DataFrame(results_data).dropna(axis=1, how='all')
    fig = None
    if models_results and reader.data_time is not None:
        fig = create_interactive_plot(plot_config, models_results, component)
    return fig, df_res, msg

# --- INTERFAZ GRADIO MEJORADA ---
def create_gradio_interface() -> gr.Blocks:
    def change_language(lang_key: str) -> Dict:
        lang = Language[lang_key]
        trans = TRANSLATIONS.get(lang, TRANSLATIONS[Language.ES])
        return trans["title"], trans["subtitle"]

    MODEL_CHOICES = [(model.display_name, model.name) for model in AVAILABLE_MODELS.values()]
    DEFAULT_MODELS = [m.name for m in list(AVAILABLE_MODELS.values())[:4]]

    with gr.Blocks(theme=THEMES["light"], css="""
        .gradio-container {font-family: 'Inter', sans-serif;}
        .theory-box {background-color: #f0f9ff; padding: 20px; border-radius: 10px; margin: 10px 0;}
        .dark .theory-box {background-color: #1e293b;}
        .model-card {border: 1px solid #e5e7eb; padding: 15px; border-radius: 8px; margin: 10px 0;}
        .dark .model-card {border-color: #374151;}
    """) as demo:
        current_theme = gr.State("light")
        current_language = gr.State("ES")
        with gr.Row():
            with gr.Column(scale=8):
                title_text = gr.Markdown("# 🔬 Analizador de Cinéticas de Bioprocesos")
                subtitle_text = gr.Markdown("Análisis avanzado de modelos matemáticos biotecnológicos")
            with gr.Column(scale=2):
                with gr.Row():
                    theme_toggle = gr.Checkbox(label="🌙 Modo Oscuro", value=False)
                    language_select = gr.Dropdown(choices=[(lang.value, lang.name) for lang in Language], value="ES", label="🌐 Idioma")
        with gr.Tabs() as tabs:
            with gr.TabItem("📚 Teoría y Modelos"):
                gr.Markdown("## Introducción a los Modelos Cinéticos\nLos modelos cinéticos en biotecnología describen el comportamiento dinámico de los microorganismos.")
                for model_name, model in AVAILABLE_MODELS.items():
                    with gr.Accordion(f"📊 {model.display_name}", open=False):
                        with gr.Row():
                            with gr.Column(scale=3):
                                gr.Markdown(f"**Descripción**: {model.description}\n\n**Ecuación**: ${model.equation}$\n\n**Parámetros**: {', '.join(model.param_names)}\n\n**Referencia**: {model.reference}")
                            with gr.Column(scale=1):
                                gr.Markdown(f"**Características**:\n- Parámetros: {model.num_params}\n- Complejidad: {'⭐' * min(model.num_params, 5)}")
            with gr.TabItem("🔬 Análisis"):
                with gr.Row():
                    with gr.Column(scale=1):
                        file_input = gr.File(label="📁 Sube tu archivo Excel (.xlsx)", file_types=['.xlsx'])
                        exp_names_input = gr.Textbox(label="🏷️ Nombres de Experimentos", placeholder="Experimento 1\nExperimento 2\n...", lines=3)
                        model_selection_input = gr.CheckboxGroup(choices=MODEL_CHOICES, label="📊 Modelos a Probar", value=DEFAULT_MODELS)
                        with gr.Accordion("⚙️ Opciones Avanzadas", open=False):
                            use_de_input = gr.Checkbox(label="Usar Evolución Diferencial", value=False, info="Optimización global más robusta pero más lenta")
                            maxfev_input = gr.Number(label="Iteraciones máximas", value=50000)
                    with gr.Column(scale=2):
                        component_selector = gr.Dropdown(choices=[("Todos los componentes", "all"), ("Solo Biomasa", C_BIOMASS), ("Solo Sustrato", C_SUBSTRATE), ("Solo Producto", C_PRODUCT)], value="all", label="📈 Componente a visualizar")
                        plot_output = gr.Plot(label="Visualización Interactiva")
                analyze_button = gr.Button("🚀 Analizar y Graficar", variant="primary")
            with gr.TabItem("📊 Resultados"):
                status_output = gr.Textbox(label="Estado del Análisis", interactive=False)
                results_table = gr.DataFrame(label="Tabla de Resultados", wrap=True)
                with gr.Row():
                    download_excel = gr.Button("📥 Descargar Excel")
                    download_json = gr.Button("📥 Descargar JSON")
                download_file = gr.File(label="Archivo descargado")
        def run_analysis_wrapper(file, models, component, use_de, maxfev, exp_names, theme):
            try:
                return run_analysis(file, models, component, use_de, maxfev, exp_names, 'dark' if theme else 'light')
            except Exception as e:
                print(f"--- ERROR EN ANÁLISIS ---\n{traceback.format_exc()}")
                return None, pd.DataFrame(), f"Error: {str(e)}"
        analyze_button.click(fn=run_analysis_wrapper, inputs=[file_input, model_selection_input, component_selector, use_de_input, maxfev_input, exp_names_input, theme_toggle], outputs=[plot_output, results_table, status_output])
        language_select.change(fn=change_language, inputs=[language_select], outputs=[title_text, subtitle_text])
        def apply_theme(is_dark):
            return gr.Info("Tema cambiado. Los gráficos nuevos usarán el tema seleccionado.")
        theme_toggle.change(fn=apply_theme, inputs=[theme_toggle], outputs=[])
        def download_results_excel(df):
            if df is None or df.empty:
                gr.Warning("No hay datos para descargar")
                return None
            with tempfile.NamedTemporaryFile(delete=False, suffix=".xlsx") as tmp:
                df.to_excel(tmp.name, index=False)
                return tmp.name
        def download_results_json(df):
            if df is None or df.empty:
                gr.Warning("No hay datos para descargar")
                return None
            with tempfile.NamedTemporaryFile(delete=False, suffix=".json") as tmp:
                df.to_json(tmp.name, orient='records', indent=2)
                return tmp.name
        download_excel.click(fn=download_results_excel, inputs=[results_table], outputs=[download_file])
        download_json.click(fn=download_results_json, inputs=[results_table], outputs=[download_file])
    return demo

# --- PUNTO DE ENTRADA ---
if __name__ == '__main__':
    gradio_app = create_gradio_interface()
    gradio_app.launch()