import numpy as np import pandas as pd import statsmodels.formula.api as smf import statsmodels.api as sm import plotly.graph_objects as go from plotly.subplots import make_subplots from scipy.optimize import minimize import plotly.express as px from scipy.stats import t, f import gradio as gr import io import os from zipfile import ZipFile class RSM_BoxBehnken: def __init__(self, data, x1_name, x2_name, x3_name, y_name, x1_levels, x2_levels, x3_levels): # ... (El código de la clase RSM_BoxBehnken se mantiene igual, solo se modifican las funciones que generan dataframes o strings) self.data = data.copy() self.model = None self.model_simplified = None self.optimized_results = None self.optimal_levels = None self.x1_name = x1_name self.x2_name = x2_name self.x3_name = x3_name self.y_name = y_name # Niveles originales de las variables self.x1_levels = x1_levels self.x2_levels = x2_levels self.x3_levels = x3_levels def get_levels(self, variable_name): if variable_name == self.x1_name: return self.x1_levels elif variable_name == self.x2_name: return self.x2_levels elif variable_name == self.x3_name: return self.x3_levels else: raise ValueError(f"Variable desconocida: {variable_name}") def fit_model(self): formula = f'{self.y_name} ~ {self.x1_name} + {self.x2_name} + {self.x3_name} + ' \ f'I({self.x1_name}**2) + I({self.x2_name}**2) + I({self.x3_name}**2) + ' \ f'{self.x1_name}:{self.x2_name} + {self.x1_name}:{self.x3_name} + {self.x2_name}:{self.x3_name}' self.model = smf.ols(formula, data=self.data).fit() print("Modelo Completo:") print(self.model.summary()) return self.model, self.pareto_chart(self.model, "Pareto - Modelo Completo") def fit_simplified_model(self): formula = f'{self.y_name} ~ {self.x1_name} + {self.x2_name} + ' \ f'I({self.x1_name}**2) + I({self.x2_name}**2) + I({self.x3_name}**2)' self.model_simplified = smf.ols(formula, data=self.data).fit() print("\nModelo Simplificado:") print(self.model_simplified.summary()) return self.model_simplified, self.pareto_chart(self.model_simplified, "Pareto - Modelo Simplificado") def optimize(self, method='Nelder-Mead'): if self.model_simplified is None: print("Error: Ajusta el modelo simplificado primero.") return def objective_function(x): return -self.model_simplified.predict(pd.DataFrame({self.x1_name: [x[0]], self.x2_name: [x[1]], self.x3_name: [x[2]]})) bounds = [(-1, 1), (-1, 1), (-1, 1)] x0 = [0, 0, 0] self.optimized_results = minimize(objective_function, x0, method=method, bounds=bounds) self.optimal_levels = self.optimized_results.x optimal_levels_natural = [ round(self.coded_to_natural(self.optimal_levels[0], self.x1_name), 3), round(self.coded_to_natural(self.optimal_levels[1], self.x2_name), 3), round(self.coded_to_natural(self.optimal_levels[2], self.x3_name), 3) ] optimization_table = pd.DataFrame({ 'Variable': [self.x1_name, self.x2_name, self.x3_name], 'Nivel Óptimo (Natural)': optimal_levels_natural, 'Nivel Óptimo (Codificado)': [round(x, 3) for x in self.optimal_levels] }) return optimization_table def plot_rsm_individual(self, fixed_variable, fixed_level): if self.model_simplified is None: print("Error: Ajusta el modelo simplificado primero.") return None varying_variables = [var for var in [self.x1_name, self.x2_name, self.x3_name] if var != fixed_variable] x_natural_levels = self.get_levels(varying_variables[0]) y_natural_levels = self.get_levels(varying_variables[1]) x_range_natural = np.linspace(x_natural_levels[0], x_natural_levels[-1], 100) y_range_natural = np.linspace(y_natural_levels[0], y_natural_levels[-1], 100) x_grid_natural, y_grid_natural = np.meshgrid(x_range_natural, y_range_natural) x_grid_coded = self.natural_to_coded(x_grid_natural, varying_variables[0]) y_grid_coded = self.natural_to_coded(y_grid_natural, varying_variables[1]) prediction_data = pd.DataFrame({ varying_variables[0]: x_grid_coded.flatten(), varying_variables[1]: y_grid_coded.flatten(), }) prediction_data[fixed_variable] = self.natural_to_coded(fixed_level, fixed_variable) z_pred = self.model_simplified.predict(prediction_data).values.reshape(x_grid_coded.shape) varying_variables = [var for var in [self.x1_name, self.x2_name, self.x3_name] if var != fixed_variable] fixed_level_coded = self.natural_to_coded(fixed_level, fixed_variable) subset_data = self.data[np.isclose(self.data[fixed_variable], fixed_level_coded)] valid_levels = [-1, 0, 1] experiments_data = subset_data[ subset_data[varying_variables[0]].isin(valid_levels) & subset_data[varying_variables[1]].isin(valid_levels) ] experiments_x_natural = experiments_data[varying_variables[0]].apply(lambda x: self.coded_to_natural(x, varying_variables[0])) experiments_y_natural = experiments_data[varying_variables[1]].apply(lambda x: self.coded_to_natural(x, varying_variables[1])) fig = go.Figure(data=[go.Surface(z=z_pred, x=x_grid_natural, y=y_grid_natural, colorscale='Viridis', opacity=0.7, showscale=True)]) for i in range(x_grid_natural.shape[0]): fig.add_trace(go.Scatter3d( x=x_grid_natural[i, :], y=y_grid_natural[i, :], z=z_pred[i, :], mode='lines', line=dict(color='gray', width=2), showlegend=False, hoverinfo='skip' )) for j in range(x_grid_natural.shape[1]): fig.add_trace(go.Scatter3d( x=x_grid_natural[:, j], y=y_grid_natural[:, j], z=z_pred[:, j], mode='lines', line=dict(color='gray', width=2), showlegend=False, hoverinfo='skip' )) colors = ['red', 'blue', 'green', 'purple', 'orange', 'yellow', 'cyan', 'magenta'] point_labels = [] for i, row in experiments_data.iterrows(): point_labels.append(f"{row[self.y_name]:.2f}") fig.add_trace(go.Scatter3d( x=experiments_x_natural, y=experiments_y_natural, z=experiments_data[self.y_name], mode='markers+text', marker=dict(size=4, color=colors[:len(experiments_x_natural)]), text=point_labels, textposition='top center', name='Experimentos' )) fig.update_layout( scene=dict( xaxis_title=varying_variables[0] + " (g/L)", yaxis_title=varying_variables[1] + " (g/L)", zaxis_title=self.y_name, ), title=f"{self.y_name} vs {varying_variables[0]} y {varying_variables[1]}
{fixed_variable} fijo en {fixed_level:.2f} (g/L) (Modelo Simplificado)", height=800, width=1000, showlegend=True ) return fig def generate_all_plots(self): if self.model_simplified is None: print("Error: Ajusta el modelo simplificado primero.") return levels_to_plot_natural = { self.x1_name: self.x1_levels, self.x2_name: self.x2_levels, self.x3_name: self.x3_levels } figs = [] for fixed_variable in [self.x1_name, self.x2_name, self.x3_name]: for level in levels_to_plot_natural[fixed_variable]: fig = self.plot_rsm_individual(fixed_variable, level) if fig is not None: figs.append(fig) return figs def coded_to_natural(self, coded_value, variable_name): levels = self.get_levels(variable_name) return levels[0] + (coded_value + 1) * (levels[-1] - levels[0]) / 2 def natural_to_coded(self, natural_value, variable_name): levels = self.get_levels(variable_name) return -1 + 2 * (natural_value - levels[0]) / (levels[-1] - levels[0]) def pareto_chart(self, model, title): tvalues = model.tvalues[1:] abs_tvalues = np.abs(tvalues) sorted_idx = np.argsort(abs_tvalues)[::-1] sorted_tvalues = abs_tvalues[sorted_idx] sorted_names = tvalues.index[sorted_idx] alpha = 0.05 dof = model.df_resid t_critical = t.ppf(1 - alpha / 2, dof) fig = px.bar( x=sorted_tvalues, y=sorted_names, orientation='h', labels={'x': 'Efecto Estandarizado', 'y': 'Término'}, title=title ) fig.update_yaxes(autorange="reversed") fig.add_vline(x=t_critical, line_dash="dot", annotation_text=f"t crítico = {t_critical:.2f}", annotation_position="bottom right") return fig def get_simplified_equation(self): if self.model_simplified is None: print("Error: Ajusta el modelo simplificado primero.") return None coefficients = self.model_simplified.params equation = f"{self.y_name} = {coefficients['Intercept']:.3f}" for term, coef in coefficients.items(): if term != 'Intercept': if term == f'{self.x1_name}': equation += f" + {coef:.3f}*{self.x1_name}" elif term == f'{self.x2_name}': equation += f" + {coef:.3f}*{self.x2_name}" elif term == f'{self.x3_name}': equation += f" + {coef:.3f}*{self.x3_name}" elif term == f'I({self.x1_name} ** 2)': equation += f" + {coef:.3f}*{self.x1_name}^2" elif term == f'I({self.x2_name} ** 2)': equation += f" + {coef:.3f}*{self.x2_name}^2" elif term == f'I({self.x3_name} ** 2)': equation += f" + {coef:.3f}*{self.x3_name}^2" return equation def generate_prediction_table(self): if self.model_simplified is None: print("Error: Ajusta el modelo simplificado primero.") return None self.data['Predicho'] = self.model_simplified.predict(self.data) self.data['Residual'] = self.data[self.y_name] - self.data['Predicho'] prediction_table = self.data[[self.y_name, 'Predicho', 'Residual']].copy() prediction_table[self.y_name] = prediction_table[self.y_name].round(3) prediction_table['Predicho'] = prediction_table['Predicho'].round(3) prediction_table['Residual'] = prediction_table['Residual'].round(3) return prediction_table def calculate_contribution_percentage(self): if self.model_simplified is None: print("Error: Ajusta el modelo simplificado primero.") return None anova_table = sm.stats.anova_lm(self.model_simplified, typ=2) ss_total = anova_table['sum_sq'].sum() contribution_table = pd.DataFrame({ 'Factor': [], 'Suma de Cuadrados': [], '% Contribución': [] }) for index, row in anova_table.iterrows(): if index != 'Residual': factor_name = index if factor_name == f'I({self.x1_name} ** 2)': factor_name = f'{self.x1_name}^2' elif factor_name == f'I({self.x2_name} ** 2)': factor_name = f'{self.x2_name}^2' elif factor_name == f'I({self.x3_name} ** 2)': factor_name = f'{self.x3_name}^2' ss_factor = row['sum_sq'] contribution_percentage = (ss_factor / ss_total) * 100 contribution_table = pd.concat([contribution_table, pd.DataFrame({ 'Factor': [factor_name], 'Suma de Cuadrados': [round(ss_factor, 3)], '% Contribución': [round(contribution_percentage, 3)] })], ignore_index=True) return contribution_table def calculate_detailed_anova(self): if self.model_simplified is None: print("Error: Ajusta el modelo simplificado primero.") return None formula_reduced = f'{self.y_name} ~ {self.x1_name} + {self.x2_name} + {self.x3_name} + ' \ f'I({self.x1_name}**2) + I({self.x2_name}**2) + I({self.x3_name}**2)' model_reduced = smf.ols(formula_reduced, data=self.data).fit() anova_reduced = sm.stats.anova_lm(model_reduced, typ=2) ss_total = np.sum((self.data[self.y_name] - self.data[self.y_name].mean())**2) df_total = len(self.data) - 1 ss_regression = anova_reduced['sum_sq'][:-1].sum() df_regression = len(anova_reduced) - 1 ss_residual = self.model_simplified.ssr df_residual = self.model_simplified.df_resid replicas = self.data[self.data.duplicated(subset=[self.x1_name, self.x2_name, self.x3_name], keep=False)] ss_pure_error = replicas.groupby([self.x1_name, self.x2_name, self.x3_name])[self.y_name].var().sum() df_pure_error = len(replicas) - len(replicas.groupby([self.x1_name, self.x2_name, self.x3_name])) ss_lack_of_fit = ss_residual - ss_pure_error df_lack_of_fit = df_residual - df_pure_error ms_regression = ss_regression / df_regression ms_residual = ss_residual / df_residual ms_lack_of_fit = ss_lack_of_fit / df_lack_of_fit ms_pure_error = ss_pure_error / df_pure_error f_lack_of_fit = ms_lack_of_fit / ms_pure_error p_lack_of_fit = 1 - f.cdf(f_lack_of_fit, df_lack_of_fit, df_pure_error) detailed_anova_table = pd.DataFrame({ 'Fuente de Variación': ['Regresión', 'Residual', 'Falta de Ajuste', 'Error Puro', 'Total'], 'Suma de Cuadrados': [round(ss_regression, 3), round(ss_residual, 3), round(ss_lack_of_fit, 3), round(ss_pure_error, 3), round(ss_total, 3)], 'Grados de Libertad': [df_regression, df_residual, df_lack_of_fit, df_pure_error, df_total], 'Cuadrado Medio': [round(ms_regression, 3), round(ms_residual, 3), round(ms_lack_of_fit, 3), round(ms_pure_error, 3), np.nan], 'F': [np.nan, np.nan, round(f_lack_of_fit, 3), np.nan, np.nan], 'Valor p': [np.nan, np.nan, round(p_lack_of_fit, 3), np.nan, np.nan] }) ss_curvature = anova_reduced['sum_sq'][f'I({self.x1_name} ** 2)'] + anova_reduced['sum_sq'][f'I({self.x2_name} ** 2)'] + anova_reduced['sum_sq'][f'I({self.x3_name} ** 2)'] df_curvature = 3 detailed_anova_table.loc[len(detailed_anova_table)] = ['Curvatura', round(ss_curvature, 3), df_curvature, round(ss_curvature / df_curvature, 3), np.nan, np.nan] detailed_anova_table = detailed_anova_table.reindex([0, 5, 1, 2, 3, 4]) detailed_anova_table = detailed_anova_table.reset_index(drop=True) return detailed_anova_table # --- Funciones para la interfaz de Gradio --- def load_data(x1_name, x2_name, x3_name, y_name, x1_levels_str, x2_levels_str, x3_levels_str, data_str): try: x1_levels = [float(x.strip()) for x in x1_levels_str.split(',')] x2_levels = [float(x.strip()) for x in x2_levels_str.split(',')] x3_levels = [float(x.strip()) for x in x3_levels_str.split(',')] data_list = [row.split(',') for row in data_str.strip().split('\n')] column_names = ['Exp.', x1_name, x2_name, x3_name, y_name] data = pd.DataFrame(data_list, columns=column_names) data = data.apply(pd.to_numeric, errors='coerce') if not all(col in data.columns for col in column_names): raise ValueError("El formato de los datos no es correcto.") global rsm rsm = RSM_BoxBehnken(data, x1_name, x2_name, x3_name, y_name, x1_levels, x2_levels, x3_levels) return data, x1_name, x2_name, x3_name, y_name, x1_levels, x2_levels, x3_levels, gr.update(visible=True) except Exception as e: return None, "", "", "", "", [], [], [], gr.update(visible=False), f"Error: {e}" def fit_and_optimize_model(): if 'rsm' not in globals(): return None, None, None, None, None, None, "Error: Carga los datos primero." model_completo, pareto_completo = rsm.fit_model() model_simplificado, pareto_simplificado = rsm.fit_simplified_model() optimization_table = rsm.optimize() equation = rsm.get_simplified_equation() prediction_table = rsm.generate_prediction_table() contribution_table = rsm.calculate_contribution_percentage() anova_table = rsm.calculate_detailed_anova() equation_formatted = equation.replace(" + ", "
+ ").replace(" ** ", "^").replace("*", " × ") equation_formatted = f"### Ecuación del Modelo Simplificado:
{equation_formatted}" return model_completo.summary().as_html(), pareto_completo, model_simplificado.summary().as_html(), pareto_simplificado, equation_formatted, optimization_table, prediction_table, contribution_table, anova_table current_plot_index = 0 plot_images = [] def generate_rsm_plot(fixed_variable, fixed_level, request: gr.Request): global current_plot_index, plot_images if 'rsm' not in globals(): return None, "Error: Carga los datos primero.", None, None if not plot_images: plot_images = rsm.generate_all_plots() if not plot_images: return None, "Error: No se pudieron generar los gráficos.", None, None current_plot_index = (current_plot_index) % len(plot_images) fig = plot_images[current_plot_index] img_bytes = fig.to_image(format="png") # Crear un archivo temporal para guardar la imagen temp_file = os.path.join(request.kwargs['temp_dir'], f"plot_{current_plot_index}.png") with open(temp_file, "wb") as f: f.write(img_bytes) return fig, "", temp_file, gr.update(visible=True) def download_excel(): if 'rsm' not in globals(): return None, "Error: Carga los datos primero." output = io.BytesIO() with pd.ExcelWriter(output, engine='xlsxwriter') as writer: rsm.data.to_excel(writer, sheet_name='Datos', index=False) rsm.generate_prediction_table().to_excel(writer, sheet_name='Predicciones', index=False) rsm.optimize().to_excel(writer, sheet_name='Optimizacion', index=False) rsm.calculate_contribution_percentage().to_excel(writer, sheet_name='Contribucion', index=False) rsm.calculate_detailed_anova().to_excel(writer, sheet_name='ANOVA', index=False) output.seek(0) return gr.File(value=output, visible=True, filename="resultados_rsm.xlsx") def download_images(request: gr.Request): global plot_images if 'rsm' not in globals(): return None, "Error: Carga los datos primero." if not plot_images: return None, "Error: No se han generado gráficos." zip_filename = "graficos_rsm.zip" zip_path = os.path.join(request.kwargs['temp_dir'], zip_filename) with ZipFile(zip_path, 'w') as zipf: for i, fig in enumerate(plot_images): img_bytes = fig.to_image(format="png") img_path = os.path.join(request.kwargs['temp_dir'], f"plot_{i}.png") with open(img_path, "wb") as f: f.write(img_bytes) zipf.write(img_path, f"plot_{i}.png") return gr.File(value=zip_path, visible=True, filename=zip_filename) def next_plot(): global current_plot_index current_plot_index += 1 return current_plot_index def prev_plot(): global current_plot_index current_plot_index -= 1 return current_plot_index # --- Crear la interfaz de Gradio --- with gr.Blocks() as demo: gr.Markdown("# Optimización de la producción de AIA usando RSM Box-Behnken") with gr.Row(): with gr.Column(): gr.Markdown("## Configuración del Diseño") x1_name_input = gr.Textbox(label="Nombre de la Variable X1 (ej. Glucosa)", value="Glucosa") x2_name_input = gr.Textbox(label="Nombre de la Variable X2 (ej. Extracto de Levadura)", value="Extracto_de_Levadura") x3_name_input = gr.Textbox(label="Nombre de la Variable X3 (ej. Triptófano)", value="Triptofano") y_name_input = gr.Textbox(label="Nombre de la Variable Dependiente (ej. AIA (ppm))", value="AIA_ppm") x1_levels_input = gr.Textbox(label="Niveles de X1 (separados por comas)", value="1, 3.5, 5.5") x2_levels_input = gr.Textbox(label="Niveles de X2 (separados por comas)", value="0.03, 0.2, 0.3") x3_levels_input = gr.Textbox(label="Niveles de X3 (separados por comas)", value="0.4, 0.65, 0.9") data_input = gr.Textbox(label="Datos del Experimento (formato CSV)", lines=5, value="""1,-1,-1,0,166.594 2,1,-1,0,177.557 3,-1,1,0,127.261 4,1,1,0,147.573 5,-1,0,-1,188.883 6,1,0,-1,224.527 7,-1,0,1,190.238 8,1,0,1,226.483 9,0,-1,-1,195.550 10,0,1,-1,149.493 11,0,-1,1,187.683 12,0,1,1,148.621 13,0,0,0,278.951 14,0,0,0,297.238 15,0,0,0,280.896""") load_button = gr.Button("Cargar Datos") with gr.Column(): gr.Markdown("## Datos Cargados") data_output = gr.Dataframe(label="Tabla de Datos") # Hacer que la sección de análisis sea visible solo después de cargar los datos with gr.Row(visible=False) as analysis_row: with gr.Column(): fit_button = gr.Button("Ajustar Modelo y Optimizar") download_excel_button = gr.Button("Descargar Tablas en Excel") gr.Markdown("**Modelo Completo**") model_completo_output = gr.HTML() pareto_completo_output = gr.Plot() gr.Markdown("**Modelo Simplificado**") model_simplificado_output = gr.HTML() pareto_simplificado_output = gr.Plot() equation_output = gr.HTML() optimization_table_output = gr.Dataframe(label="Tabla de Optimización", headers=["Variable", "Nivel Óptimo (Natural)", "Nivel Óptimo (Codificado)"]) prediction_table_output = gr.Dataframe(label="Tabla de Predicciones") contribution_table_output = gr.Dataframe(label="Tabla de % de Contribución") anova_table_output = gr.Dataframe(label="Tabla ANOVA Detallada") with gr.Column(): gr.Markdown("## Generar Gráficos de Superficie de Respuesta") fixed_variable_input = gr.Dropdown(label="Variable Fija", choices=["Glucosa", "Extracto_de_Levadura", "Triptofano"], value="Glucosa") fixed_level_input = gr.Slider(label="Nivel de Variable Fija", minimum=0, maximum=1, step=0.01, value=0.5) with gr.Row(): plot_button = gr.Button("Generar Gráfico") download_images_button = gr.Button("Descargar Gráficos en ZIP") prev_plot_button = gr.Button("<") next_plot_button = gr.Button(">") rsm_plot_output = gr.Plot() download_plot_button = gr.Button("Descargar Gráfico Actual") plot_image_output = gr.File(label="Gráfico Actual", visible=False) load_button.click( load_data, inputs=[x1_name_input, x2_name_input, x3_name_input, y_name_input, x1_levels_input, x2_levels_input, x3_levels_input, data_input], outputs=[data_output, x1_name_input, x2_name_input, x3_name_input, y_name_input, x1_levels_input, x2_levels_input, x3_levels_input, analysis_row] ) fit_button.click(fit_and_optimize_model, outputs=[model_completo_output, pareto_completo_output, model_simplificado_output, pareto_simplificado_output, equation_output, optimization_table_output, prediction_table_output, contribution_table_output, anova_table_output]) plot_button.click(generate_rsm_plot, inputs=[fixed_variable_input, fixed_level_input], outputs=[rsm_plot_output, equation_output, plot_image_output, download_plot_button]) download_excel_button.click(download_excel, outputs=download_excel_button, api_name="download_excel") download_images_button.click(download_images, outputs=download_images_button, api_name="download_images") download_plot_button.click(lambda x: x, inputs=[plot_image_output], outputs=[plot_image_output], api_name="download_plot") prev_plot_button.click(prev_plot, outputs=prev_plot_button).then(generate_rsm_plot, inputs=[fixed_variable_input, fixed_level_input], outputs=[rsm_plot_output, equation_output, plot_image_output, download_plot_button]) next_plot_button.click(next_plot, outputs=next_plot_button).then(generate_rsm_plot, inputs=[fixed_variable_input, fixed_level_input], outputs=[rsm_plot_output, equation_output, plot_image_output, download_plot_button]) # Ejemplo de uso gr.Markdown("## Ejemplo de uso") gr.Markdown("1. Introduce los nombres de las variables y sus niveles en las cajas de texto correspondientes.") gr.Markdown("2. Copia y pega los datos del experimento en la caja de texto 'Datos del Experimento'.") gr.Markdown("3. Haz clic en 'Cargar Datos' para cargar los datos en la tabla.") gr.Markdown("4. Haz clic en 'Ajustar Modelo y Optimizar' para ajustar el modelo y encontrar los niveles óptimos de los factores.") gr.Markdown("5. Selecciona una variable fija y su nivel en los controles deslizantes.") gr.Markdown("6. Haz clic en 'Generar Gráfico' para generar un gráfico de superficie de respuesta.") gr.Markdown("7. Usa '<' y '>' para navegar entre los gráficos generados.") gr.Markdown("8. Haz clic en 'Descargar Tablas en Excel' para obtener un archivo Excel con todas las tablas generadas.") gr.Markdown("9. Haz clic en 'Descargar Gráfico Actual' para descargar la imagen del gráfico actual en formato PNG.") gr.Markdown("10. Haz clic en 'Descargar Gráficos en ZIP' para descargar todas las imágenes de los gráficos en un archivo ZIP.") demo.launch()