<jupyter_start><jupyter_text>Notebook C: Evaluation of ML Algorithms for Rate Prediction of Syngas Fermentation  Setup imports<jupyter_code>import pandas as pd
import numpy as np
import sklearn.preprocessing, sklearn.neural_network, sklearn.svm, sklearn.ensemble
import scipy.stats
import matplotlib.pyplot as plt<jupyter_output><empty_output><jupyter_text>Load data that was generated in notebook A<jupyter_code>rates_df = pd.read_csv(f'../data/rates_data.csv')
print(f'Shape of the rates data: {rates_df.shape[0]} rows by {rates_df.shape[1]} columns')<jupyter_output>Shape of the rates data: 836 rows by 18 columns<jupyter_text>Create train and test sets<jupyter_code>train_data = rates_df[rates_df.composition.isin([1,2,3,4,5,6,7])]
test_data = rates_df[rates_df.composition.isin([8,9,10])]
print(f'Shape of the training data: {train_data.shape[0]} rows by {train_data.shape[1]} columns')
print(f'Shape of the testing data: {test_data.shape[0]} rows by {test_data.shape[1]} columns')<jupyter_output>Shape of the training data: 657 rows by 18 columns
Shape of the testing data: 179 rows by 18 columns<jupyter_text>Define a function that generates the input and output arrays for scikit learn's API<jupyter_code>def get_X_y_arrays(data):
    
    # prevent set with copy error
    data_copy = data.copy()
    
    # ML input
    X = data_copy [[
        'biomass (g/L)', 'ethanol (mM)', 'acetate (mM)', 'butanol (mM)', 
         'butyrate (mM)', 'N2', 'CO', 'CO2', 'H2', 'flow rate (mL/min)'
    ]]
    
    # ML output
    y = data_copy [[
        'biomass rate', 'ethanol rate', 'acetate rate', 'butanol rate', 'butyrate rate'
    ]]
 
    return np.array(X), np.array(y)<jupyter_output><empty_output><jupyter_text>Split data into ML input and output<jupyter_code>X_train, y_train = get_X_y_arrays(train_data)
X_test, y_test = get_X_y_arrays(test_data)

print(f'Shape of the train X array: {X_train.shape[0]} rows by {X_train.shape[1]} columns')
print(f'Shape of the train y array: {y_train.shape[0]} rows by {y_train.shape[1]} columns')
print(f'Shape of the test X array: {X_test.shape[0]} rows by {X_test.shape[1]} columns')
print(f'Shape of the test y array: {y_test.shape[0]} rows by {y_test.shape[1]} columns')<jupyter_output>Shape of the train X array: 657 rows by 10 columns
Shape of the train y array: 657 rows by 5 columns
Shape of the test X array: 179 rows by 10 columns
Shape of the test y array: 179 rows by 5 columns<jupyter_text>Define a dictionary of models with optimized parameters from notebook B<jupyter_code>optimized_models = {
    "acetate": {
        "nn_fine": sklearn.neural_network.MLPRegressor(
            shuffle=True, 
            activation = 'tanh', 
            hidden_layer_sizes = [50, 40],
            max_iter = 5000
        ),
        "nn_coarse": sklearn.neural_network.MLPRegressor(
            shuffle=True, 
            activation = 'tanh', 
            hidden_layer_sizes = [40, 20],
            max_iter = 5000
        ),
        "svm_rbf": sklearn.svm.SVR(
            kernel = 'rbf', 
            C = 10000, 
            epsilon = 0.1, 
            gamma = 0.01
        ),
        'rf': sklearn.ensemble.RandomForestRegressor(
            max_depth = 32,
            # max_samples = 0.5,
            n_estimators = 130
        ),
        'en': sklearn.linear_model.ElasticNet(
            alpha =  0.1,
            l1_ratio = 0.4
        ),
        'lasso': sklearn.linear_model.Lasso(
            alpha = 0.1
        ),
        'knn': sklearn.neighbors.KNeighborsRegressor(
            algorithm = 'ball_tree',
            leaf_size = 5,
            n_neighbors = 4,
            weights = 'distance'
        ),
    },
    "biomass": {
        "nn_fine": sklearn.neural_network.MLPRegressor(
            shuffle=True, 
            activation = 'tanh', 
            hidden_layer_sizes = [100, 80, 60, 70],
            max_iter = 5000
        ),
        "nn_coarse": sklearn.neural_network.MLPRegressor(
            shuffle=True, 
            activation = 'tanh', 
            hidden_layer_sizes = [100, 100, 60, 80],
            max_iter = 5000
        ),
        "svm_rbf": sklearn.svm.SVR(
            kernel = 'rbf', 
            C = 10000,
            epsilon = 0.0001, 
            gamma = 0.001
        ),
        'rf': sklearn.ensemble.RandomForestRegressor(
            max_depth = 32,
            # max_samples = 0.5,
            n_estimators = 80
        ),
        'en': sklearn.linear_model.ElasticNet(
            alpha = 1e-05, 
            l1_ratio = 0.1
        ),
        'lasso': sklearn.linear_model.Lasso(
            alpha = 1e-06
        ),
        'knn': sklearn.neighbors.KNeighborsRegressor(
            algorithm = 'ball_tree',
            leaf_size = 5,
            n_neighbors = 4,
            weights = 'distance'
        ),
    },
    "butanol": {
        "nn_fine": sklearn.neural_network.MLPRegressor(
            shuffle=True, 
            activation = 'tanh', 
            hidden_layer_sizes = [90, 60, 10, 80],
            max_iter = 5000
        ),
        "nn_coarse": sklearn.neural_network.MLPRegressor(
            shuffle=True, 
            activation = 'tanh', 
            hidden_layer_sizes = [60, 20],
            max_iter = 5000
        ),
        "svm_rbf": sklearn.svm.SVR(
            kernel = 'rbf', 
            C = 1000, 
            epsilon = 0.01, 
            gamma = 0.01
        ),
        'rf': sklearn.ensemble.RandomForestRegressor(
            max_depth = 28,
            n_estimators = 120
        ),
        'en': sklearn.linear_model.ElasticNet(
            alpha =  1e-10,
            l1_ratio = 0.1
        ),
        'lasso': sklearn.linear_model.Lasso(
            alpha = 1e-10
        ),
        'knn': sklearn.neighbors.KNeighborsRegressor(
            algorithm = 'ball_tree',
            leaf_size = 5,
            n_neighbors = 2,
            weights = 'distance'
        ),
    },
    "butyrate": {
        "nn_fine": sklearn.neural_network.MLPRegressor(
            shuffle=True, 
            activation = 'tanh', 
            hidden_layer_sizes = [90, 30, 20],
            max_iter = 5000
        ),
        "nn_coarse": sklearn.neural_network.MLPRegressor(
            shuffle=True, 
            activation = 'tanh', 
            hidden_layer_sizes = [60, 20],
            max_iter = 5000
        ),
        "svm_rbf": sklearn.svm.SVR(
            kernel = 'rbf', 
            C = 10000, 
            epsilon = 0.01, 
            gamma = 0.01
        ),
        'rf': sklearn.ensemble.RandomForestRegressor(
            max_depth = 22,
            # max_samples = 0.5,
            n_estimators = 130
        ),
        'en': sklearn.linear_model.ElasticNet(
            alpha = 0.0001,
            l1_ratio = 0.1
        ),
        'lasso': sklearn.linear_model.Lasso(
            alpha = 1e-10
        ),
        'knn': sklearn.neighbors.KNeighborsRegressor(
            algorithm = 'ball_tree',
            leaf_size = 5,
            n_neighbors = 4,
            weights = 'distance'
        ),
    },
    "ethanol": {
        "nn_fine": sklearn.neural_network.MLPRegressor(
            shuffle=True, 
            activation = 'tanh', 
            hidden_layer_sizes = [80, 50],
            max_iter = 5000
        ),
        "nn_coarse": sklearn.neural_network.MLPRegressor(
            shuffle=True, 
            activation = 'tanh', 
            hidden_layer_sizes = [80, 60],
            max_iter = 5000
        ),
        "svm_rbf": sklearn.svm.SVR(
            kernel = 'rbf', 
            C = 10000, 
            epsilon = 0.0001, 
            gamma = 0.001
        ),
        'rf': sklearn.ensemble.RandomForestRegressor(
            max_depth = 22,
            # max_samples = 0.5,
            n_estimators = 100
        ),
        'en': sklearn.linear_model.ElasticNet(
            alpha =  0.001,
            l1_ratio = 0.1
        ),
        'lasso': sklearn.linear_model.Lasso(
            alpha = 0.01
        ),
        'knn': sklearn.neighbors.KNeighborsRegressor(
            algorithm = 'ball_tree',
            leaf_size = 5,
            n_neighbors = 4,
            weights = 'distance'
        ),
    },
}<jupyter_output><empty_output><jupyter_text>Get the predictions of the models on the training set and test set<jupyter_code># Make a dictionary to save trained models 
trained_models = {}

# Define two dictionaries to hold data for training and test result data
model_train_output_dictionary = {}
model_test_output_dictionary = {}

# Scale both sets of X data
Scaler = sklearn.preprocessing.MinMaxScaler()
X_scaled_train = Scaler.fit_transform(X_train, y_test)
X_scaled_test = Scaler.fit_transform(X_test, y_test)

for index, output in enumerate(['biomass', 'ethanol', 'acetate', 'butanol', 'butyrate']):
    print(output)
    
    # separate out the output of interest
    y_train_output=y_train[:,index]
    y_test_output=y_test[:,index]
    
    # save the measured data in the output dictionaries
    model_train_output_dictionary[f'{output} measured'] = y_train_output
    model_test_output_dictionary[f'{output} measured'] = y_test_output
    
    for algorithm in ['nn_fine', 'nn_coarse', 'svm_rbf', 'rf', 'en', 'lasso', 'knn']:
        print(algorithm)

        # train the model and save it
        trained_model = optimized_models[output][algorithm].fit(X_scaled_train, y_train_output)
        trained_models[f'{output} {algorithm}'] = trained_model
        
        # test the model on training and test data
        y_train_prediction = trained_model.predict(X_scaled_train)
        y_test_prediction = trained_model.predict(X_scaled_test)
        
        # save the outputs in dictionaries
        model_train_output_dictionary[f'{output} {algorithm}'] = y_train_prediction
        model_test_output_dictionary[f'{output} {algorithm}'] = y_test_prediction
        
# convert the dictionaries to dataframes
train_prediction_df = pd.DataFrame(model_train_output_dictionary)
test_prediction_df = pd.DataFrame(model_test_output_dictionary)

# display the data frames
display(train_prediction_df.head())
display(test_prediction_df.head())<jupyter_output>biomass
nn_fine
nn_coarse
svm_rbf
rf
en
lasso
knn
ethanol
nn_fine
nn_coarse
svm_rbf
rf
en
lasso
knn
acetate
nn_fine
nn_coarse
svm_rbf
rf
en
lasso
knn
butanol
nn_fine
nn_coarse
svm_rbf
rf
en
lasso
knn
butyrate
nn_fine<jupyter_text>Save dataframes as .csv for excel based figures<jupyter_code>train_prediction_df.to_csv('../data/train_predictions.csv')
test_prediction_df.to_csv('../data/test_predictions.csv')<jupyter_output><empty_output><jupyter_text>Plot scatter plots of fits<jupyter_code>plt.rcParams.update({'font.size': 11.5})
fig, ax = plt.subplots(5, 7, figsize=(35, 30))
fig.suptitle(f'Rate Predictions')

algorithm_dictionary = {
    'nn_fine': 'Neural Network Fine',
    'nn_coarse': 'Neural Network Coarse',
    'svm_rbf': 'Support Vector', 
    'rf' : 'Random Forest',
    'en': 'Elastic Net', 
    'lasso' : 'Lasso',
    'knn': 'k-Nearest Neighbors'
}

for x_index, output in enumerate(['acetate', 'biomass', 'butanol', 'butyrate', 'ethanol']):
    train_r2_list = []
    test_r2_list = []
    for y_index, algorithm in enumerate(['nn_fine', 'nn_coarse', 'svm_rbf', 'rf', 'en', 'lasso', 'knn']):
 
        # Set title, x label, and y label of the subplot
        ax[x_index, y_index].title.set_text(f'{output.capitalize()} {algorithm_dictionary[algorithm]}')
        if output != 'biomass':
            ax[x_index, y_index].set_xlabel(f'Measured Production Rate (mmol {output}/day)')
            ax[x_index, y_index].set_ylabel(f'Predicted Production Rate (mmol {output}/day)')
        else:
            ax[x_index, y_index].set_xlabel(f'Measured Production Rate (g biomass/day)')
            ax[x_index, y_index].set_ylabel(f'Predicted Production Rate (g biomass/day)')
            
        # Plot training and test data by plotting predicted Production Rates vs measured Production Rates
        ax[x_index, y_index].scatter(train_prediction_df[f'{output} measured'], train_prediction_df[f'{output} {algorithm}'])
        ax[x_index, y_index].scatter(test_prediction_df[f'{output} measured'], test_prediction_df[f'{output} {algorithm}'])
        
        # ensure that all plots are square, so that 45 degree line looks the same in all plots
        minimum = min(pd.concat([
            train_prediction_df[f'{output} measured'],
            train_prediction_df[f'{output} {algorithm}'],
            test_prediction_df[f'{output} measured'],
            test_prediction_df[f'{output} {algorithm}']
        ], axis=0))
        maximum = max(pd.concat([
            train_prediction_df[f'{output} measured'],
            train_prediction_df[f'{output} {algorithm}'],
            test_prediction_df[f'{output} measured'],
            test_prediction_df[f'{output} {algorithm}']
        ], axis=0))
        
        # calculate train data r-squared
        measured_train_data = train_prediction_df[f'{output} measured']
        predicted_train_data = train_prediction_df[f'{output} {algorithm}']
        # could use scikitlearn's r-squared value
        # train_r2 = sklearn.metrics.r2_score(measured_train_data, predicted_train_data)
        _, _, train_r_value, _, _ = scipy.stats.linregress(measured_train_data, predicted_train_data)
        train_r2_list.append(train_r_value**2)

        # calculate test data r-squared
        measured_test_data = test_prediction_df[f'{output} measured']
        predicted_test_data = test_prediction_df[f'{output} {algorithm}']
        # could use scikitlearn's r-squared value
        # test_r2 = sklearn.metrics.r2_score(measured_test_data, predicted_test_data)
        _, _, test_r_value, _, _ = scipy.stats.linregress(measured_test_data, predicted_test_data)
        test_r2_list.append(test_r_value**2)
        
        # add legend to the plot
        ax[x_index, y_index].legend([f'train {train_r_value**2:.3f}', f'test {test_r_value**2:.3f}'])
        
        # plot the subplot
        ax[x_index, y_index].plot([minimum, maximum], [minimum, maximum], 'k--') #row=0, col=0

plt.savefig(f'../figures/scatter_plots.png', dpi=100)
plt.show()<jupyter_output><empty_output><jupyter_text>Plot bar graph of r squared values<jupyter_code>plt.rcParams.update({'font.size': 11.5})
for output in ['acetate', 'biomass', 'butanol', 'butyrate', 'ethanol']:
    train_r2_list = []
    test_r2_list = []
    for algorithm in ['nn_fine', 'nn_coarse', 'svm_rbf', 'rf', 'en', 'lasso', 'knn']:

        # Calculate training data r-squared
        measured_train_data = train_prediction_df[f'{output} measured']
        predicted_train_data = train_prediction_df[f'{output} {algorithm}']
        # could use scikitlearn's r-squared value
        # train_r2 = sklearn.metrics.r2_score(measured_train_data, predicted_train_data)
        _, _, train_r_value, _, _ = scipy.stats.linregress(measured_train_data, predicted_train_data)
        train_r2_list.append(train_r_value**2)

        # Calculate test data r-squared
        measured_test_data = test_prediction_df[f'{output} measured']
        predicted_test_data = test_prediction_df[f'{output} {algorithm}']
        # could use scikitlearn's r-squared value
        # test_r2 = sklearn.metrics.r2_score(measured_test_data, predicted_test_data)
        _, _, test_r_value, _, _ = scipy.stats.linregress(measured_test_data, predicted_test_data)
        test_r2_list.append(test_r_value**2)
    
    # define values to help with plotting
    x = np.arange(6)  # the label locations
    width = 0.25  # the width of the bars

    # plot training and testing r-squared values
    fig, ax = plt.subplots(figsize=(9,5))
    rects1 = ax.bar(x - .5 * width, train_r2_list,      width, label='Train')
    rects3 = ax.bar(x + .5 * width, test_r2_list, width, label='Test')

    # Add labels for title and custom x-axis tick labels, etc.
    ax.set_ylabel('r-squared Value')
    ax.set_title(f'{output.capitalize()} Prediction Quality')
    ax.set_xticks(x)
    ax.set_xticklabels(['Neural Net', 'Support Vector', 'Random Forest', 'Elastic Net', 'Lasso', 'K-nearest'])
    ax.legend()

    fig.tight_layout()
    plt.savefig(f'../figures/bar_graph_{output}.png', dpi=100)
#     plt.show()<jupyter_output><empty_output><jupyter_text>Make data frames to plot feature importance of metabolite and gases<jupyter_code>features = ['acetate', 'biomass', 'butanol', 'butyrate', 'ethanol', 'CO', 'CO2', 'H2', 'flow rate']

metabolites = ['acetate', 'biomass', 'butanol', 'butyrate', 'ethanol']
gases = ['CO', 'CO2', 'H2', 'flow rate']

metabolite_feature_importance_list = []
gases_feature_importance_list = []

for output in ['acetate', 'biomass', 'butanol', 'butyrate', 'ethanol']:
    
    # get trained random forest model for the output
    model = trained_models[f'{output} rf']
    
    # get feature importance
    feature_importance_array = model.feature_importances_
    # print(feature_importance_array)
    
    # label feature importances
    feature_importance_dictionary = {feature: importance for feature, importance in zip(features, feature_importance_array)}
    
    # split dictionary into metabolites and gases and add the split dictionaries to lists
    metabolite_feature_importance_list.append({metabolite: feature_importance_dictionary[metabolite] for metabolite in metabolites})
    gases_feature_importance_list.append({gas: feature_importance_dictionary[gas] for gas in gases})
    
# convert list of dictionaries to dataframes
metabolite_feature_importance_df = pd.DataFrame(metabolite_feature_importance_list, index=metabolites)
gases_feature_importance_df = pd.DataFrame(gases_feature_importance_list, index=metabolites)

# display the dataframes
display(metabolite_feature_importance_df)
display(gases_feature_importance_df)<jupyter_output><empty_output><jupyter_text>Plot metabolite feature importance<jupyter_code>display(metabolite_feature_importance_df)

# This helps with formatting 
x = np.arange(5)  # the label locations
width = 0.15  # the width of the bars

fig, ax = plt.subplots(figsize=(8,5))

for counter, (index, row) in enumerate(metabolite_feature_importance_df.transpose().iterrows()):
    output_feature_importances = list(row)
    
    rects = ax.bar(x - 0.3 + 0.15*counter, output_feature_importances , width, label=f'{index.capitalize()} Concentration')

# # Add some text for labels, title and custom x-axis tick labels, etc.
ax.set_ylabel(f'Importance')

ax.set_title(f'Feature Importance of Metabolites')
ax.set_xticks(x)
labels = ['Acetate Rate', 'Biomass Rate', 'Butanol Rate', 'Butyrate Rate', 'Ethanol Rate']
ax.set_xticklabels(labels)
ax.legend()
# ax.set_yscale('log')

fig.tight_layout()
plt.savefig(f'../figures/feature_importance_metabolites.png', dpi=100)
plt.show()<jupyter_output><empty_output><jupyter_text>Plot gases feature importance<jupyter_code>display(gases_feature_importance_df)

# This helps with formatting 
x = np.arange(5)  # the label locations
width = 0.15  # the width of the bars

fig, ax = plt.subplots(figsize=(8,5))

for counter, (index, row) in enumerate(gases_feature_importance_df.transpose().iterrows()):
    output_feature_importances = list(row)
    if index != 'flow rate':
        index = index + ' Concentration' 
    else:
        index = 'Flow Rate'
    rects = ax.bar(x - 0.3 + 0.15*counter, output_feature_importances , width, label=f'{index}')

# # Add some text for labels, title and custom x-axis tick labels, etc.
ax.set_ylabel(f'Importance')

ax.set_title(f'Feature Importance of Gases')
ax.set_xticks(x)
labels = ['Acetate Rate', 'Biomass Rate', 'Butanol Rate', 'Butyrate Rate', 'Ethanol Rate']
ax.set_xticklabels(labels)
ax.legend()
# ax.set_yscale('log')

fig.tight_layout()
plt.savefig(f'../figures/feature_importance_gases.png', dpi=100)
plt.show()<jupyter_output><empty_output><jupyter_text>Plot feature importance of each gas for each output<jupyter_code>display(gases_feature_importance_df)

plt.rcParams.update({'font.size': 11.5})
fig, ax = plt.subplots(2, 2, figsize=(15, 10))
fig.suptitle(f'Feature importance of each gas for each output')

labels = [label.capitalize() for label in gases_feature_importance_df.index]

ax[0,0].title.set_text(f'Importance of CO on Metabolite Production Rate')
ax[0,0].set_ylabel(f'Importance')
ax[0,0].axis(ymin=0.0,ymax=0.07)
rects = ax[0,0].bar(labels, gases_feature_importance_df['CO'])

ax[0,1].title.set_text(f'Importance of CO2 on Metabolite Production Rate')
ax[0,1].set_ylabel(f'Importance')
ax[0,1].axis(ymin=0.0,ymax=0.07)
rects = ax[0,1].bar(labels, gases_feature_importance_df['CO2'])

ax[1,0].title.set_text(f'Importance of H2 on Metabolite Production Rate')
ax[1,0].set_ylabel(f'Importance')
ax[1,0].axis(ymin=0.0,ymax=0.07)
rects = ax[1,0].bar(labels, gases_feature_importance_df['H2'])

ax[1,1].title.set_text(f'Importance of Flow Rate on Metabolite Production Rate')
ax[1,1].set_ylabel(f'Importance')
ax[1,1].axis(ymin=0.0,ymax=0.07)
rects = ax[1,1].bar(labels, gases_feature_importance_df['flow rate'])

plt.savefig(f'../figures/feature_importance_per_gas.png', dpi=100)

plt.show()<jupyter_output><empty_output><jupyter_text>Save gases_feature_importance as csv<jupyter_code>gases_feature_importance_df.to_csv('../data/gases_importance.csv')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1187135053682426, 'fraction_numerical': 0.021896779054435864, 'mean_word_length': 3.0903312788906008, 'pattern_counts': {'":': 20, '<': 57, '<?xml version=': 0, '>': 57, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '405364', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 7761, 'n_tokens_neox': 7146, 'n_words': 1485, 'path': 'notebooks/C_model_evalutation.ipynb', 'repo_name': 'NKWBTB/SyngasMachineLearning', 'size': '21236'}