# models.py

import numpy as np
from scipy.optimize import curve_fit
from sympy import symbols, sympify, lambdify
import warnings

class BioprocessModel:
    def __init__(self):
        self.params = {}
        self.models = {}
        self.r2 = {}
        self.rmse = {}

    def set_model(self, model_type, equation_str, param_str):
        equation_str = equation_str.strip()
        if '=' in equation_str:
            equation_str = equation_str.split('=', 1)[1].strip()

        params = [param.strip() for param in param_str.split(',')]
        self.models[model_type] = {
            'equation_str': equation_str,
            'params': params
        }

        t = symbols('t')
        param_symbols = symbols(params)
        expr = sympify(equation_str)
        func = lambdify((t, *param_symbols), expr, 'numpy')
        self.models[model_type]['function'] = func

    def fit_model(self, model_type, time, data, bounds):
        func = self.models[model_type]['function']
        params = self.models[model_type]['params']

        p0 = np.ones(len(params))
        lower_bounds, upper_bounds = bounds

        lower_bounds = np.array(lower_bounds)
        upper_bounds = np.array(upper_bounds)

        if len(lower_bounds) != len(params):
            lower_bounds = np.full(len(params), -np.inf)
        if len(upper_bounds) != len(params):
            upper_bounds = np.full(len(params), np.inf)

        with warnings.catch_warnings():
            warnings.simplefilter("ignore")
            popt, _ = curve_fit(func, time, data, p0=p0, bounds=(lower_bounds, upper_bounds), maxfev=10000)

        self.params[model_type] = dict(zip(params, popt))
        y_pred = func(time, *popt)
        ss_res = np.sum((data - y_pred) ** 2)
        ss_tot = np.sum((data - np.mean(data)) ** 2)
        self.r2[model_type] = 1 - (ss_res / ss_tot)
        self.rmse[model_type] = np.sqrt(np.mean((data - y_pred) ** 2))
        return y_pred