Feature Resolution and Residuals Experiment¶

This notebook uses the SageWorks Science Workbench to quickly build an AWS® Machine Learning Pipeline with the AQSolDB public dataset. For this exercise we're going to look at the relationship between feature space and target values, specifically we're going to use SageWorks to help us identify areas where compounds that are close in feature space have significant differences in their target values (solubility in this case).

Data¶

AqSolDB: A curated reference set of aqueous solubility, created by the Autonomous Energy Materials Discovery [AMD] research group, consists of aqueous solubility values of 9,982 unique compounds curated from 9 different publicly available aqueous solubility datasets. AqSolDB also contains some relevant topological and physico-chemical 2D descriptors. Additionally, AqSolDB contains validated molecular representations of each of the compounds. This openly accessible dataset, which is the largest of its kind, and will not only serve as a useful reference source of measured and calculated solubility data, but also as a much improved and generalizable training data source for building data-driven models. (2019-04-10)

Main Reference: https://www.nature.com/articles/s41597-019-0151-1

Data Dowloaded from the Harvard DataVerse: https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/OVHAW8

® Amazon Web Services, AWS, the Powered by AWS logo, are trademarks of Amazon.com, Inc. or its affiliates.

In [1]:

import os
os.environ["SAGEWORKS_CONFIG"] = "/Users/briford/.sageworks/ideaya_sandbox.json"

In [2]:

import sageworks
import logging
logging.getLogger("sageworks").setLevel(logging.WARNING)

In [3]:

# We've already created a FeatureSet so just grab it a sample
from sageworks.api.feature_set import FeatureSet
fs = FeatureSet("test_sol_nightly_log_s")
full_df = fs.pull_dataframe()

In [4]:

full_df.head()

In [5]:

# df["class"].value_counts()

In [6]:

# Sanity check our solubility and solubility_class
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
plt.rcParams['font.size'] = 12.0
plt.rcParams['figure.figsize'] = 14.0, 5.0
sns.set_theme(style='darkgrid')

# Create a box plot
ax = sns.boxplot(x='class', y='log_s', data=full_df, order = [2, 1, 0])
plt.title('Solubility by Solubility Class')
plt.xlabel('Solubility Class')
plt.ylabel('Solubility')
plt.show()

In [7]:

# Let's get the current model metrics
from sageworks.api.model import Model
model = Model("test-sol-nightly-regression")

features = model.features()
target = model.target()
id_column = "udm_mol_bat_id"
model.performance_metrics()

In [8]:

# Okay lets look at our reference model (already deployed to an Endpoint)
from sageworks.api.endpoint import Endpoint
end = Endpoint("test-sol-nightly-regression-end")
pred_df = end.auto_inference()

In [11]:

plot_predictions(pred_df)

In [102]:

show_columns = [id_column, target, "prediction", "residuals", "residuals_abs", "udm_asy_protocol", "smiles"]
high_residuals = pred_df[(pred_df["residuals_abs"] > 2.0) & (pred_df["log_s"] != -16.0)][show_columns]
high_residuals

Feature Spider¶

In [ ]:

# from sageworks.utils.endpoint_utils import fs_training_data
# training_df = fs_training_data(end)

In [ ]:

# training_df.shape

In [17]:

from sageworks.algorithms.dataframe.feature_spider import FeatureSpider

feature_spider = FeatureSpider(full_df, features=features, target_column=target, id_column=id_column)

In [18]:

pred_df["confidence"] = feature_spider.confidence_scores(pred_df, pred_df["prediction"])

In [19]:

plot_predictions(pred_df, color="confidence")

In [20]:

neighbor_info = feature_spider.neighbor_info(pred_df)

In [103]:

pd.set_option("display.width", 1000)
pd.set_option("display.max_colwidth", 50)
pd.set_option("display.max_rows", None)
# idya_high_residuals = high_residuals[high_residuals["udm_asy_protocol"] == "IDYA"]
# idya_high_residuals.sort_values(by='residuals_abs', ascending=False)
high_residuals

In [100]:

pd.set_option("display.width", 1000)
pd.set_option("display.max_colwidth", 150)
pd.set_option("display.max_rows", None)
id = "IDC-5965-1"
index = high_residuals[high_residuals["udm_mol_bat_id"]==id].index[0]
print(index)
print(high_residuals.loc[index])
neighbor_info.iloc[index]

So What?¶

Blah, blah, who cares? The real model works fine...

So lets switch gears and look at the predictions on the real model¶

In [105]:

logging.getLogger("sageworks").setLevel(logging.INFO)

In [107]:

nightly_end = Endpoint("solubility-class-0-nightly")
nightly_pred = nightly_end.auto_inference()

In [130]:

# ReInitialize the feature spider with "class" as the target column
feature_spider = FeatureSpider(full_df, features=features, target_column="class", id_column=id_column)

In [132]:

nightly_pred["confidence"] = feature_spider.confidence_scores(nightly_pred, nightly_pred["prediction"])
nightly_pred["confidence"].hist()

In [133]:

show_columns = ["udm_mol_bat_id", "udm_asy_protocol", "udm_asy_cnd_format", "udm_asy_res_value", "class", "prediction"]

# Filter rows where 'class' and 'prediction' differ by at least 2
filtered_df = nightly_pred[abs(nightly_pred["class"] - nightly_pred["prediction"]) >= 2][show_columns]
filtered_df

In [136]:

pd.set_option("display.max_colwidth", 150)
id = "IDC-26481-1"
index = nightly_pred[nightly_pred["udm_mol_bat_id"]==id].index[0]
print(index)
feature_spider.neighbor_info(nightly_pred).iloc[index]

In [ ]:

high_residuals = pred_df[pred_df["residuals_abs"] > 4.0][show_columns + ["confidence"]]
high_residuals

Feature Resolution¶

In [ ]:

from sageworks.algorithms.dataframe.feature_resolution import FeatureResolution

# Grab Model and Get the target and feature columns
m = Model("aqsol-mol-regression-full")
target_column = m.target()
feature_columns = m.features()

# Create the class and run the report
resolution = FeatureResolution(df, features=feature_columns, target_column=target_column, id_column="id")
resolve_df = resolution.compute(within_distance=0.05, min_target_difference=0.5)
print(resolve_df)

Hypothesis¶

There should be some overlap of the predictions with high residuals and the features that had resolution issues

In [ ]:

show_columns = ["id", "solubility", "prediction", "residuals", "smiles"]
high_residuals = pred_df[pred_df["residuals_abs"] > 2.0][show_columns]
high_residuals

In [ ]:

query_ids = ["B-1027", "A-3099"]
show_columns = ["id", "solubility", "prediction", "residuals"]
pred_df[pred_df["id"].isin(query_ids)][show_columns]

In [ ]:

high_error_ids = set(high_residuals["id"])
feature_issue_ids = set(list(resolve_df["id"]) + list(resolve_df["n_id"]))

In [ ]:

# Output the overlap
print(f"High Model Error Ids: {len(high_error_ids)}")
print(f"Feature Resolution Ids: {len(feature_issue_ids)}")
print(f"Overlap: {len(high_error_ids.intersection(feature_issue_ids))}")

In [ ]:

feature_issue_ids - high_error_ids

In [ ]:

pred_df[pred_df["id"]=='F-465'][show_columns]

In [ ]:

feature_issue_id

In [ ]:

# Lets construct another model
import xgboost as xgb
from sklearn.model_selection import train_test_split

# Grab all the data and do a Train/Test split
full_df = fs.pull_dataframe()
train_df, test_df = train_test_split(full_df, test_size=0.2, random_state=42)


X = train_df[features]
y = train_df[target]

# Train the main XGBRegressor model
model = xgb.XGBRegressor()
model.fit(X, y)

# Run Predictions on the the hold out
test_df["prediction"] = model.predict(test_df[features])
plot_predictions(test_df)

In [ ]:

len(full_df.columns)

In [ ]:

from sageworks.algorithms.dataframe.feature_spider import FeatureSpider

# Create the FeatureSpider class and run the various methods
f_spider = FeatureSpider(full_df, features, id_column="id", target_column=target, neighbors=5)
coincident = f_spider.coincident(target_diff=1)
print("COINCIDENT")
print(coincident)
high_gradients = f_spider.high_gradients(within_distance=0.25, target_diff=2)
print("\nHIGH GRADIENTS")
print(high_gradients)

In [ ]:

len([6913, 644, 1799, 9739, 2834, 5014, 8729, 8741, 5292, 1713, 320, 7360, 5698, 9923, 1355, 4305, 5713, 4820, 6874, 7004, 8165, 3174, 364, 1781, 6140, 4735])

In [ ]:

from rdkit import Chem
from rdkit.Chem import Draw
from rdkit.Chem.Draw.rdMolDraw2D import SetDarkMode

In [ ]:

test_df = FeatureSet("solubility_features_test

In [ ]:

def show(id):
    smile = df[df["id"]==id]["smiles"].values[0]
    print(smile)
    _features = df[df["id"]==id][features].values[0]
    # print(_features)
    _target = df[df["id"]==id][target].values[0]
    print(_target)
    return Chem.MolFromSmiles(smile)

In [ ]:

show("IDC-10845-1")

In [ ]:

show("B-3121")

In [ ]:

show("B-962")

In [ ]:

show("C-846")

In [ ]:

show("A-2377")

In [ ]:

show("B-3122")

In [ ]:

from sageworks.algorithms.dataframe.target_gradients import TargetGradients
gradients = TargetGradients()
grad_df = gradients.compute(df, features, target, min_target_distance=1.0)
grad_df.head()

In [ ]:

grad_df[grad_df["target_gradient"] == float("inf")]

In [ ]:

import numpy as np
finite_gradients = grad_df[grad_df["target_gradient"].apply(np.isfinite)]
finite_gradients.plot.hist(y="target_gradient", bins=20, alpha=0.5, logy=True)

In [ ]:

# Remove Observations with Coincident Features and large target differences
print(f"Before Removal: {len(full_df)}")
indexes_to_remove = list(set(coincident).union(set(high_gradients)))
clean_df = full_df.drop(indexes_to_remove).copy().reset_index(drop=True)
print(f"After Removal: {len(clean_df)}")

In [ ]:

f_spider = FeatureSpider(clean_df, features, id_column="id", target_column=target)
preds = f_spider.predict(clean_df)
print(preds)
coincident = f_spider.coincident(target_diff=1)
print("COINCIDENT")
print(coincident)
high_gradients = f_spider.high_gradients(within_distance=0.25, target_diff=2)
print("\nHIGH GRADIENTS")
print(high_gradients)

In [ ]:

print(f"Before Removal: {len(clean_df)}")
indexes_to_remove = list(set(coincident).union(set(high_gradients)))
clean_df = clean_df.drop(indexes_to_remove).copy().reset_index(drop=True)
print(f"After Removal: {len(clean_df)}")

In [ ]:

train_df, test_df = train_test_split(clean_df, test_size=0.2, random_state=42)


X = clean_df[features]
y = clean_df[target]

# Train the main XGBRegressor model
model = xgb.XGBRegressor()
model.fit(X, y)

# Run Predictions on the the hold out
test_df["prediction"] = model.predict(test_df[features])
plot_predictions(test_df)

In [ ]:

f_spider = FeatureSpider(clean_df, features, id_column="id", target_column=target)
preds = f_spider.predict(clean_df)
print(preds)
coincident = f_spider.coincident(target_diff=0.1)
print("COINCIDENT")
print(coincident)
high_gradients = f_spider.high_gradients(within_distance=0.25, target_diff=0.1)
print("\nHIGH GRADIENTS")
print(high_gradients)

In [ ]:

print(f"Before Removal: {len(clean_df)}")
indexes_to_remove = list(set(coincident).union(set(high_gradients)))
clean_df = clean_df.drop(indexes_to_remove).copy().reset_index(drop=True)
print(f"After Removal: {len(clean_df)}")

In [ ]:

clean_df

In [ ]:

# Spin up a test model so get some info
target = model.target()
print(target)
features = model.features()
print(features)

In [ ]:

import xgboost as xgb
xgb_model = xgb.XGBRegressor()

# Grab our Features, Target and Train the Model
y = clean_df[target]
X = clean_df[features]
xgb_model.fit(X, y)
predictions = xgb_model.predict(clean_df[features])
clean_df["prediction"] = predictions

In [ ]:

plot_predictions(clean_df)

In [ ]:

clean_train_df, clean_test_df = train_test_split(clean_df, test_size=0.2, random_state=42)
# Grab our Features, Target and Train the Model
y = clean_train_df[target]
X = clean_train_df[features]
xgb_model.fit(X, y)
predictions = xgb_model.predict(clean_test_df[features])
clean_test_df["prediction"] = predictions
plot_predictions(clean_test_df)

In [ ]:

list(df[df["id"].isin(["C-2581", "B-13"])]["smiles"])

In [ ]:

Chem.MolFromSmiles("CO[C@H]1[C@@H](C[C@@H]2CN3CCc4c([nH]c5cc(OC)ccc45)[C@H]3C[C@@H]2[C@@H]1C(=O)OC)OC(=O)c6cc(OC)c(OC)c(OC)c6")

In [ ]:

Chem.MolFromSmiles("COC1C(CC2CN3CCC4=C([NH]C5=C4C=CC(=C5)OC)C3CC2C1C(=O)OC)OC(=O)C6=CC(=C(OC)C(=C6)OC)OC")

In [ ]:

df[df["id"].isin(["C-2581", "B-13"])][features]

In [ ]:

df[df["id"].isin(["C-2581", "B-13"])][target]

In [ ]:

df[df["id"]=="B-13"]

In [ ]:

df[df["id"].isin(["C-846", "B-962"])]

In [ ]:

show("C-846")

In [ ]:

show("B-962")

In [ ]:

Chem.MolFromSmiles("CCC(=O)OC(CC1=CC=CC=C1)(C(C)CN(C)C)C2=CC=CC=C2")

In [ ]:

Chem.MolFromSmiles("CCC(=O)O[C@@](Cc1ccccc1)([C@H](C)CN(C)C)c2ccccc2")

In [ ]:

def show(id):
    smile = df[df["id"]==id]["smiles"].values[0]
    print(smile)
    _features = df[df["id"]==id][features].values[0]
    print(_features)
    _target = df[df["id"]==id][target].values[0]
    print(_target)
    return Chem.MolFromSmiles(smile)

In [ ]:

close_ids = ["E-1200", "A-3473", "B-1665"]
show("E-1200")

In [ ]:

show("A-3473")

In [ ]:

show("B-1665")

In [ ]:

from rdkit import Chem

# Create an RDKit molecule from a SMILES string
smiles = "CCC(=O)OC(CC1=CC=CC=C1)(C(C)CN(C)C)C2=CC=CC=C2"
mol = Chem.MolFromSmiles(smiles)

# Assign stereochemistry using RDKit
Chem.AssignStereochemistry(mol, cleanIt=True, force=True)

# Find chiral centers and their configurations
chiral_centers = Chem.FindMolChiralCenters(mol, includeUnassigned=True)

# Print the results
for center in chiral_centers:
    index, configuration = center
    print(f"Atom index: {index}, Configuration: {configuration}")

In [ ]:

# Create an RDKit molecule from a SMILES string
smiles = "CCC(=O)O[C@@](Cc1ccccc1)([C@H](C)CN(C)C)c2ccccc2"
mol = Chem.MolFromSmiles(smiles)

# Assign stereochemistry using RDKit
Chem.AssignStereochemistry(mol, cleanIt=True, force=True)

# Find chiral centers and their configurations
chiral_centers = Chem.FindMolChiralCenters(mol, includeUnassigned=True)

# Print the results
for center in chiral_centers:
    index, configuration = center
    print(f"Atom index: {index}, Configuration: {configuration}")

In [ ]:

from sageworks.algorithms.dataframe.row_tagger import RowTagger
row_tagger = RowTagger(
    df,
    features=features,
    id_column="id",
    target_column=target,
    min_dist=0.0,
    min_target_diff=1.0,
)
data_df = row_tagger.tag_rows()
print(data_df)

In [ ]:

data_df["tags"].value_counts()

In [ ]:

data_df[data_df["tags"]=={"htg", "coincident"}]

Helper Methods¶

In [10]:

# Helper to look at predictions vs target
from math import sqrt
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
plt.rcParams['font.size'] = 12.0
plt.rcParams['figure.figsize'] = 14.0, 5.0
sns.set_theme(style='darkgrid')
def plot_predictions(df, line=True, color="PredError"):
    
    # Dataframe of the targets and predictions
    target = 'Actual Solubility'
    pred = 'Predicted Solubility'
    df_plot = pd.DataFrame({target: df['log_s'], pred: df['prediction']})
    
    # Compute Error per prediction
    if color == "PredError":
        df_plot["PredError"] = df_plot.apply(lambda x: abs(x[pred] - x[target]), axis=1)
    else:
        df_plot[color] = df[color]

    #df_plot['error'] = df_plot.apply(lambda x: abs(x[pred] - x[target]), axis=1)
    ax = df_plot.plot.scatter(x=target, y=pred, c=color, cmap='coolwarm', sharex=False)
    
    # Just a diagonal line
    if line:
        ax.axline((1, 1), slope=1, linewidth=2, c='black')
        x_pad = (df_plot[target].max() - df_plot[target].min())/10.0 
        y_pad = (df_plot[pred].max() - df_plot[pred].min())/10.0
        plt.xlim(df_plot[target].min()-x_pad, df_plot[target].max()+x_pad)
        plt.ylim(df_plot[pred].min()-y_pad, df_plot[pred].max()+y_pad)

    

In [ ]: