model_input.py

# IMPORT THE MAIN PACKAGES NEEDED FOR THE PROGRAM TO RUN
# Standard packages
import numpy as np
from scipy.interpolate import interp1d
import pandas as pd
from typing import Callable, Dict, List, Any, Tuple, Union
import logging
# Custom packages
from general_utils.save_df import*
from general_utils.read_file import filter_and_convert_to_time
'''
Author: Davide Carecci
Date: January 2025
Institution: Politecnico di Milano
Module for preparing model inputs for simulation.
Contains functions to convert input compositions (from agri-AcoDM to AM2HN-like format), evaluate flow rates, and prepare model inputs for simulation.
'''
def convert_input_composition(input_name: str, data: np.ndarray, inorganicC: float = 0, 
                              output_path = os.getcwd(), save_to_csv: bool = False) -> Tuple[pd.DataFrame, Dict[str, List[Callable[[float], float]]]]:
    '''
    Convert input composition data from agri-AcoDM format (CSV or CombiTimeTable) to AM2HN-like format.
    Following the nomenclature in the manual of BIOGoAlS.Twin and the author PhD thesis, the output are the $\boldsymbol{\theta_{in,i}}$ vectors,
    for the i-th co-feedstock in the input diet. 
    In practice, the outputs are interpolating functions for each component of the influent composition (X_h, S_1, S_2, Z, C, N),
    because it is considered that such compositions can vary over time.
    Parameters:
        input_name (str): Name of the input feedstock (e.g., 'maize', 'cowslurry', 'tomatosauce').
        data (np.ndarray): Input composition data read from CSV or CombiTimeTable. The first column must be time in days.
        inorganicC (float, optional): Inorganic dissolved carbon content. Defaults to 0.
        output_path (str, optional): Path to save output CSV file. Defaults to current working directory.
        save_to_csv (bool, optional): Whether to save the converted composition to a CSV file. Defaults to False.
    Returns:
        df (pd.DataFrame): DataFrame containing the converted composition data (Timestamp and composition components as columns).
        dict: Dictionary containing interp1d functions for each component of the influent composition.

    Note: data are the same of CombiTimeTables created from CSV files (see the BIOGoAlS.agri-AcoDM and BIOGoAlS.Twin manuals for details). 
        Must have the following columns in order:
        Time [days], TS [gTS/kgFM], VS [%TS], pH [-], TAC [mgCaCO3/L], TAN [mgN/L], ortoP [mgP/L], Ca [mol/L], Mg [mol/L],
        Acetate [gCOD/L], Propionate [gCOD/L], Butirate [gCOD/L], Valerate [gCOD/L], Sugars [%TS], Aminoacids [%TS], Long-Ch.Fatty Acids [%TS], th_BMP [NmL_ch4/g_VS], 
        Crude Protein [%TS], Lipid [%TS], Cellulose [%TS], Hemicell. [%TS], Lignin [%TS], Starch [%TS], BD_pr [%], BD_li [%], BD_cell [%], BD_hemicell [%].
    Note: 'inorganicC' must be set manually from function call coherent with those computed from the Modelica simulation of each feedstock's conversion model,
            where a change balance given the pH data is enforced to derive the dissolved inorganic carbon content.
    '''
    # Convert composition from agri-AcoDM to AM2HN-like format
    comp = [data[:,1]/100*(data[:,17]*data[:,23]/100 + data[:,18]*data[:,24]/100 + data[:,19]*data[:,25]/100 + data[:,20]*data[:,26]/100 + data[:,22]),
            data[:,1]/100*(1.07*data[:,13] + 1.53*data[:,14] + 2.878*data[:,15]),
            data[:,9]/64*1000+data[:,10]/112*1000+data[:,11]/160*1000+data[:,12]/208*1000,
            data[:,4]/100*2, # Note: no TAN/14 here! # NO Z OF MAIZE IS SEEN BY ADM1 BECCOUSE LIQUID=FALSE, BUT COMBI != 0
            inorganicC*np.ones(len(data[:,0])),
            data[:,5]/14
            ]
    # Create DataFrame and save to CSV if needed
    df = pd.DataFrame(np.column_stack((data[:,0], np.matrix(comp).T)), columns = ['Time','Xh', 'S1', 'S2', 'Z', 'C', 'N'])
    if save_to_csv:
        df.to_csv(os.path.join(output_path, input_name+'_read_composition.csv'), index=False)

    # Create interp1d functions for each component
    Xh_in = interp1d(data[:,0], comp[0], kind='previous', fill_value=(comp[0][0], comp[0][-1]), bounds_error=False)
    S1_in = interp1d(data[:,0], comp[1], kind='previous', fill_value=(comp[1][0], comp[1][-1]), bounds_error=False)
    S2_in = interp1d(data[:,0], comp[2], kind='previous', fill_value=(comp[2][0], comp[2][-1]), bounds_error=False)
    Z_in = interp1d(data[:,0], comp[3], kind='previous', fill_value=(comp[3][0], comp[3][-1]), bounds_error=False)
    C_in = interp1d(data[:,0], comp[4], kind='previous', fill_value=(comp[4][0], comp[4][-1]), bounds_error=False)
    N_in = interp1d(data[:,0], comp[5], kind='previous', fill_value=(comp[5][0], comp[5][-1]), bounds_error=False)

    return df, dict({input_name: [Xh_in, S1_in, S2_in, Z_in, C_in, N_in]})

def u(t: np.ndarray, items_dict: list, mult_factors: list) -> list:
    """
    Generalized function to evaluate flow rates for a dynamic list of items.
    Args:
        t (float): The time at which to evaluate the functions.
        items_dict (dict): A dictionary where keys are item names and values are lists of interp1d functions.
        mult_factors (list): List of multiplicative factors for each item.
    Returns:
        list: Dilution rates computed for each item.
    """
    return [mult_factor * item_functions[0](t) for mult_factor,item_functions in zip(mult_factors,items_dict.values())] 

def param_u(t: np.ndarray, items_dict: list) -> list:
    """
    Generalized function to compute a matrix of parameters for a dynamic list of items.
    Args:
        t (float): The time at which to evaluate the functions.
        items_dict (dict): A dictionary where keys are item names and values are lists of interp1d functions.
    Returns:
        list: A matrix (list of lists) where each row corresponds to a parameter across items.
    """
    # Number of parameters is assumed to be the length of each value in the dictionary
    num_params = len(next(iter(items_dict.values())))  # Get the length of the first list of functions
    return [
        [items_dict[item][param_index](t) for item in items_dict]  # Evaluate each parameter for all items
        for param_index in range(num_params)
    ]

def prepare_model_inputs(modelname: str, maize: pd.DataFrame, cowslurry: pd.DataFrame, tomatosauce: pd.DataFrame, 
                         d_flow: pd.DataFrame, u_flow: pd.DataFrame, 
                         x: pd.DataFrame, 
                         integrator_parameters: Dict, mult_factors: List = [1,1,1], t_span_end: float = None,
                         save_to_csv: bool = False, output_path: str = os.getcwd(), start=None, end=None, fill_option: List=["zeros", "zeros", "zeros"], 
                         sic: List[float] = [0,193,0], log: bool = False) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
    '''
    Function to prepare model inputs before simulation. 
    This function is very case-specific and bug-prone if inputs are not carefully prepared before its call. 
    However, it wraps all the necessary steps to prepare the inputs (i.e. parameters, initial state, influent characterization and flow rates) before simulation.
    Nominally, should be used on the outputs of the 'read_all' function in 'model_read_file.py'.
    From the dataframes and dictionaries read from files, it returns all the necessary inputs in numpy arrays for simulation.
    It filter and convert inputs from agri-AcoDM to AM2HN-like format. 
    It defines initial state, u_values, param_u_values, time_points, timestamp_points, delta_t, for an integration coherent with the 'integrator_parameters'.
    
    Parameters:
        modelname (str): Name of the model (for file saving files).
        maize, cowslurry, tomatosauce, d_flow, u_flow, x (pd.DataFrame): DataFrames containing input data read from files.
        integrator_parameters (Dict): Dictionary containing integrator parameters (keys: 'start_timestamp', 'end_timestamp', 'mesh_finesse').
        mult_factors (List, optional): Multiplicative factors for flow rates. Defaults to [1,1,1].
        t_span_end (float, optional): End time for simulation span. Defaults to None. To be used only if I don't want to simulate until the end of the input data.
        save_to_csv (bool, optional): Whether to save prepared inputs to CSV files. Defaults to False.
        output_path (str, optional): Path to save output files. Defaults to current working directory.
        start (pd.Timestamp, optional): Start timestamp for filtering input data. Defaults to None.
        end (pd.Timestamp, optional): End timestamp for filtering input data. Defaults to None.
        fill_option (List[str], optional): Options for filling missing data in filtering. Defaults to ["zeros", "zeros", "zeros"]. 
            Check 'filter_and_convert_to_time' function for details.
        sic (List[float], optional): Inorganic dissolved carbon content for maize, cowslurry, and tomatosauce. Defaults to [0,193,0].
        log (bool, optional): Whether to log information during preparation. Defaults to False.
    Returns:
        x0 (np.ndarray): Initial state values.
        timestamp_points_discrete (np.ndarray): Timestamps for discrete time points.
        time_points_discrete (np.ndarray): Discrete time points in days.
        delta_step_discrete (np.ndarray): Time step differences between discrete time points in days.
        u_values_discrete (np.ndarray): Discrete flow rate values.
        param_u_values_discrete (np.ndarray): Discrete parameter matrix values.
    
    Note: if not directly available as data, 'sic' values must be coherent with those computed from the Modelica simulation of each feedstock's conversion model,
        where a change balance given the pH data is enforced to derive the dissolved inorganic carbon content.
        (see BIOGoAlS.agri-AcoDM and BIOGoAlS.Twin manuals for details).
    Note: see also 'prepare_model_inputs' method in the 'model_input.py' file.
    Note: 'mesh_finesse' in 'integrator_parameters' is used to define a mesh over which the inputs are evaluated. It is expressed in number of points per hour.
        It is strongly recommended to have values inside the input flow rate and feedstock's composition file that does not have any row that has a Timestamp 
        not multiple of (3600/mesh_finesse) seconds from the start timestamp, otherwise unexpected behaviors may occur.
    '''
    logging.info('## Preparing inputs... ##')
    start = start if start is not None else integrator_parameters['start_timestamp']
    end = end if end is not None else integrator_parameters['end_timestamp']
    
    # Cut input data (and convert the Timestamp column to time)
    d_flow = filter_and_convert_to_time(file_path='/d_flow', df = d_flow,
                                    start_timestamp=start, 
                                    end_timestamp=end,fill_option=fill_option[0], log = log)
    u_flow = filter_and_convert_to_time(file_path='/u_flow', df = u_flow,
                                    start_timestamp=start, 
                                    end_timestamp=end,fill_option=fill_option[0], log = log)
    maize = filter_and_convert_to_time(file_path='/maize', df = maize,
                                    start_timestamp=start, 
                                    end_timestamp=end,fill_option=fill_option[1], log = log)
    cowslurry = filter_and_convert_to_time(file_path='/cowslurry', df = cowslurry,
                                    start_timestamp=start, 
                                    end_timestamp=end,fill_option=fill_option[1], log = log)
    tomatosauce = filter_and_convert_to_time(file_path='/tomatosauce', df = tomatosauce,
                                    start_timestamp=start, 
                                    end_timestamp=end,fill_option=fill_option[1], log = log)
    
    # Note: important if values are less than 1!!
    d_flow[1][:,0] = np.round(d_flow[1][:,0])
    u_flow[1][:,0] = np.round(u_flow[1][:,0])

    # Cut state data and extract initial states values
    x = filter_and_convert_to_time(file_path='/x', df=x,
                                    start_timestamp=start, 
                                    end_timestamp=end,fill_option=fill_option[2])
    x0 = x[0][x[0]['Timestamp']==0.0].values[0][1:10]
    logging.info(f'Initial condition is {dict(zip(x[0].keys()[1:], x0))}')

    # Convert input compositions from agri-AcoDM to AM2HN-like format
    maize_comp = convert_input_composition('maize', maize[1], sic[0], output_path=output_path, save_to_csv=save_to_csv)
    slurry_comp = convert_input_composition('cowslurry', cowslurry[1], sic[1], output_path=output_path, save_to_csv=save_to_csv)
    tomato_comp = convert_input_composition('tomatosauce', tomatosauce[1], sic[2], output_path=output_path, save_to_csv=save_to_csv)
    comps = dict()
    comps.update(maize_comp[1])
    comps.update(slurry_comp[1])
    comps.update(tomato_comp[1])

    # Interpolate the input flow rate signal (gFM/day)
    maize_flow = interp1d(d_flow[1][:,0], d_flow[1][:,1], kind='previous', fill_value=(d_flow[1][0,0], d_flow[1][-1,1]), bounds_error=False)
    slurry_flow = interp1d(d_flow[1][:,0], d_flow[1][:,2], kind='previous', fill_value=(d_flow[1][0,0], d_flow[1][-1,2]), bounds_error=False)
    tomato_flow = interp1d(u_flow[1][:,0], u_flow[1][:,1], kind='previous', fill_value=(u_flow[1][0,0], u_flow[1][-1,1]), bounds_error=False)
    flows = dict()
    flows.update({'maize':[maize_flow]})
    flows.update({'cowslurry': [slurry_flow]})
    flows.update({'tomatosauce': [tomato_flow]})

    # Set simulation (discrete integration points) coherent with integrator parameters
    if t_span_end == None:
        t_span = (0, (end-start).total_seconds()/86400)                                                                                     # Note: time span of the simulation
    else: 
        t_span = (0, t_span_end)
    if log:
        logging.info(f'The timespan of the simulation is {t_span}, but the flow inputs exist for {(max(d_flow[1][0,0]/86400, u_flow[1][0,0]/86400), min(d_flow[1][-1,0]/86400, u_flow[1][-1,0]/86400))} days')
    mesh = np.arange(0,np.round(t_span[1]*86400,0)+3600/integrator_parameters['mesh_finesse'],3600/integrator_parameters['mesh_finesse'])   # Note: it does not take the initial point...always 3600. 
    time_points_discrete = np.union1d(d_flow[1][:,0], mesh)                                                                                 # Note: union between a 'mesh_finesse' mesh and maize/slurry event points
    time_points_discrete = np.union1d(time_points_discrete, u_flow[1][:,0])                                                                 # Note: union between a 'mesh_finesse' mesh and tomatosauce event points
    time_points_discrete = time_points_discrete[time_points_discrete <= np.round(t_span[1]*86400,0)]                                        # Note: up to end of simulation
    delta_step_discrete = np.diff(time_points_discrete)                                                                                     # Note: difference between subsequent points in the 'mesh_finesse' mesh and maize/slurry event points
    # Verify that delta_step_discrete has only equidistant values
    if log:
        logging.info(f'Indices with different delta step {np.where(delta_step_discrete != delta_step_discrete[0])[0]}')
    
    # Compute the input values over the discrete time points
    u_test_discrete = u(time_points_discrete, flows, mult_factors)                                                                          # Note: it still includes TAN. To exclude TAN take only u(time_points_discrete,V)[:-1]
    param_u_test_discrete = param_u(time_points_discrete, comps)                                                                            # Note: to be transposed
    u_values_discrete = np.vstack(u_test_discrete).T[:-1,:]                                                                                 # Note 16.01.2025: don't take anymore the last value!!
    param_u_test_discrete = np.array(param_u_test_discrete).T[:-1,:]                                                                        # Note 16.01.2025: don't take anymore the last value!!
    param_u_values_discrete = param_u_test_discrete
    timestamp_points_discrete = (start + pd.to_timedelta(time_points_discrete-time_points_discrete[0], unit='S')).round('S')
    # Move from seconds to days
    time_points_discrete = time_points_discrete/86400
    delta_step_discrete = delta_step_discrete/86400

    # Save the input values to CSV if needed
    if save_to_csv:
        logging.info(f'Saving to csv u and param_u files...')
        df = pd.DataFrame(np.column_stack((time_points_discrete[:-1], u_values_discrete)),columns=['Time','maize','cowslurry','tomatosauce'])
        df.insert(0, 'Timestamp', timestamp_points_discrete[:-1])
        df.to_csv(os.path.join(output_path, f'{modelname}_u_values.csv'), index=False)
        # save_df_with_check(df, 'u_test_discrete.csv')
        df = pd.DataFrame(np.column_stack((time_points_discrete[:-1], param_u_values_discrete[:,0,:])),columns=['Time','Xh', 'S1', 'S2', 'Z', 'C', 'N'])
        df.insert(0, 'Timestamp', timestamp_points_discrete[:-1])
        df.to_csv(os.path.join(output_path,f'{modelname}_param_u_values_maize.csv'), index=False)
        # save_df_with_check(df, 'param_u_test_discrete_maize.csv')
        df = pd.DataFrame(np.column_stack((time_points_discrete[:-1], param_u_values_discrete[:,1,:])),columns=['Time','Xh', 'S1', 'S2', 'Z', 'C', 'N'])
        df.insert(0, 'Timestamp', timestamp_points_discrete[:-1])
        df.to_csv(os.path.join(output_path,f'{modelname}_param_u_values_cowslurry.csv'), index=False)
        # save_df_with_check(df, 'param_u_test_discrete_cowslurry.csv')
        df = pd.DataFrame(np.column_stack((time_points_discrete[:-1], param_u_values_discrete[:,2,:])),columns=['Time','Xh', 'S1', 'S2', 'Z', 'C', 'N'])
        df.insert(0, 'Timestamp', timestamp_points_discrete[:-1])
        df.to_csv(os.path.join(output_path,f'{modelname}_param_u_values_tomatosauce.csv'), index=False)

    return x0, timestamp_points_discrete, time_points_discrete, delta_step_discrete, u_values_discrete, param_u_values_discrete