<jupyter_start><jupyter_text>Attempting to use molecular data with fastai/pytorch<jupyter_code># must use new version of pandas to correctly work with pkl
!pip install pandas==0.24.1
# after restarting the runtime
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
import torch
import fastai
from fastai.tabular import *
# Doing work so we can access data from Google drive
# Note that a brief authorization step is necessary 
!pip install PyDrive
from pydrive.auth import GoogleAuth
from pydrive.drive import GoogleDrive
from google.colab import auth
from oauth2client.client import GoogleCredentials

auth.authenticate_user()
gauth = GoogleAuth()
gauth.credentials = GoogleCredentials.get_application_default()
drive = GoogleDrive(gauth)
# Now download the csv, which I have made shareable via google drive
# note that the id is the ending part of the shareable link url (after open?id=)
download = drive.CreateFile({'id': '1uc0dI2z0bPYuW2ptOM5M1cFQqXM2eXto'})
download.GetContentFile('acidic_pka_df.csv')
pka_df = pd.read_csv('acidic_pka_df.csv')
pka_df.head()
# do a 80, 10, 10 split
train_df, valid_df, test_df = np.split(
  pka_df.sample(frac=1, random_state=42),
  [int(0.80*len(pka_df)), int(0.90*len(pka_df))]  
)
X_train = train_df.values[:,1:]
y_train = train_df.values[:,0]

X_valid = valid_df.values[:,1:]
y_valid = valid_df.values[:,0]

X_test = test_df.values[:,1:]
y_test = test_df.values[:,0]<jupyter_output><empty_output><jupyter_text>Test first on simple GBRT to make sure we are still getting decent results<jupyter_code>from sklearn.ensemble import GradientBoostingRegressor
from sklearn.model_selection import cross_validate

gbrt=GradientBoostingRegressor(n_estimators=1200)
cv_scores = cross_validate(gbrt,
                            X_train,
                            y_train,
                            cv=5,
                            n_jobs=-1,
                            verbose=1,
                            scoring=('r2',
                                     'neg_mean_squared_error',
                                     'neg_mean_absolute_error'
                                    )
                           )

RMSE_scores = [np.sqrt(-1*x) for x in cv_scores['test_neg_mean_squared_error']]
print('RMSE SCORES: ')
print("scores : ", RMSE_scores)
print("mean score : ", np.mean(RMSE_scores))
print("std of score : ", np.std(RMSE_scores))
print("*** \n")

MAE_scores = [-1*x for x in cv_scores['test_neg_mean_absolute_error']]
print('MAE SCORES: ')
print("scores : ", MAE_scores)
print("mean score : ", np.mean(MAE_scores))
print("std of score : ", np.std(MAE_scores))
print("*** \n")

R2_scores = [x for x in cv_scores['test_r2']]
print('R2 SCORES: ')
print("scores : ", R2_scores)
print("mean score : ", np.mean(R2_scores))
print("std of score : ", np.std(R2_scores))
print("*** \n")<jupyter_output>[Parallel(n_jobs=-1)]: Using backend LokyBackend with 2 concurrent workers.<jupyter_text>Good, the simple GBRT model with no hyperameter opimization does alright Now let's try some neural network approaches with pytorch and fastai:<jupyter_code>X_train, y_train, X_valid, y_valid, X_test, y_test = \
map(
  partial(torch.tensor, dtype=torch.float),
  (X_train, y_train, X_valid, y_valid, X_test, y_test)
)
bs = 64
train_ds = TensorDataset(X_train, y_train)
valid_ds = TensorDataset(X_valid, y_valid)
data = DataBunch.create(train_ds, valid_ds, bs=bs)<jupyter_output><empty_output><jupyter_text>First create a really simple architecture with 2 hidden layers (256 nodes each). I'll also add the capability for some dropout because that is what seemed to work with Deepchem. I was able to get pretty good results with this neural net and Deepchem, also another paper mentioned good results (slightly better than the GBRT) with a neural net in Matlab of similar architecture<jupyter_code>import torch.nn as nn
import torch.nn.functional as F

class Simple_DNN(nn.Module):
  def __init__(self, num_features):
    super(Simple_DNN, self).__init__()
    
    # DEFINE CONSTANTS
    HIDDEN_DIM = 256
    DROPOUT_PCT = 0.0
    
    # implement two hidden layers
    self.fc1 = nn.Linear(num_features, HIDDEN_DIM)
    self.fc2 = nn.Linear(HIDDEN_DIM, HIDDEN_DIM)
    self.fc3 = nn.Linear(HIDDEN_DIM, HIDDEN_DIM)
    self.fc4 = nn.Linear(HIDDEN_DIM, 1)
    self.dropout = nn.Dropout(p=DROPOUT_PCT)
    
  def forward(self, xb):
    xb = F.relu(self.fc1(xb))
    xb = self.dropout(xb)
    xb = F.relu(self.fc2(xb))
    xb = self.dropout(xb)
    xb = F.relu(self.fc3(xb))
    xb = self.fc4(xb)
    return xb
# use fastai for finding the learning rate and the training 
# start with no dropout
learn = Learner(data, Simple_DNN(num_features=2048),
                loss_func=nn.MSELoss(),
                metrics=[root_mean_squared_error, r2_score])
learn.lr_find()
learn.recorder.plot()
learn.fit_one_cycle(50, 1e-02)<jupyter_output><empty_output><jupyter_text>It's not even overfitting (and there is no dropout). Not getting above R2 score of 0 is really bad. Let's see if we can at least get it to overfit with a small number of samples:<jupyter_code>small_train_ds = TensorDataset(X_train[0:5], y_train[0:5])
small_valid_ds = TensorDataset(X_valid[0:5], y_valid[0:5])
small_data = DataBunch.create(small_train_ds, small_valid_ds)
# use fastai for finding the learning rate and the training 
# start with no dropout
learn_small_data_test = Learner(small_data, Simple_DNN(num_features=2048),
                loss_func=nn.MSELoss(),
                metrics=root_mean_squared_error)
learn_small_data_test.lr_find()
learn_small_data_test.recorder.plot()<jupyter_output><empty_output><jupyter_text>I can't get it to work with any batch size on this small number of samples! Ok, I will try another way based on Jeremy's nice pytorch tutorial https://pytorch.org/tutorials/beginner/nn_tutorial.html<jupyter_code>from torch import optim

# creat function to fetch the model and optimizer
def get_model(num_features, lr, GPU=False, GPU_dev=None):
  if GPU:
    model = Simple_DNN(num_features)
    return model.to(GPU_dev), optim.SGD(model.parameters(), lr=lr)
  else:
    model = Simple_DNN(num_features)
    return model, optim.SGD(model.parameters(), lr=lr)

# function to calculate the loss, and do opt step in training
def loss_batch(model, loss_func, xb, yb, opt=None):
  loss = loss_func(model(xb), yb)
  
  if opt is not None:
    loss.backward()
    opt.step()
    opt.zero_grad()
    
  return loss.item(), len(xb)

# function to fit the model for number of epochs given loss func
def fit(epochs, model, loss_func, opt, train_dl, valid_dl):
  for epoch in range(epochs):
    model.train()
    for xb, yb in train_dl:
      loss_batch(model, loss_func, xb, yb, opt)
      
    model.eval()
    with torch.no_grad():
      train_losses, train_nums  = zip(
        *[loss_batch(model, loss_func, xb, yb) for xb, yb in train_dl]
      )
      trn_loss = ( 
          np.sum(np.multiply(train_losses, train_nums)) / np.sum(train_nums)
      )
      
      valid_losses, valid_nums  = zip(
        *[loss_batch(model, loss_func, xb, yb) for xb, yb in valid_dl]
      )
      val_loss = (
          np.sum(np.multiply(valid_losses, valid_nums)) / np.sum(valid_nums)
      )
      
      if epoch % 5 == 0:
        print('epoch: ', epoch, 'train_loss: ', trn_loss, 'valid_loss: ',val_loss)
# creating custom dataloader to use with GPU
class GPUDataLoader:
  def __init__(self, dl, dev):
    self.dl = dl
    self.dev = dev
    
  def __len__(self):
    return len(self.dl)
  
  def __iter__(self):
    batches = iter(self.dl)
    for b in batches:
      yield(self._move_to_GPU(*b))
      
  def _move_to_GPU(self, X, y):
    return X.to(self.dev), y.to(self.dev)

# function that fetches data either on/off GPU
def get_data(train_ds, valid_ds, bs, GPU=False, GPU_dev=None):
  if GPU and GPU_dev is not None:
    # train_dl = DataLoader(train_ds, batch_size=bs, shuffle=True),
    # valid_dl = DataLoader(valid_ds, batch_size=bs*2)
    return (
      GPUDataLoader(DataLoader(train_ds, batch_size=bs, shuffle=True), GPU_dev),
      GPUDataLoader(DataLoader(valid_ds, batch_size=bs*2), GPU_dev)
    )
  else:
    return (
      DataLoader(train_ds, batch_size=bs, shuffle=True),
      DataLoader(valid_ds, batch_size=bs*2)
    )
small_train_ds = TensorDataset(X_train[0:5], y_train[0:5])
small_valid_ds = TensorDataset(X_valid[0:5], y_valid[0:5])
# Finally, train the model!

dev = torch.device('cuda') if torch.cuda.is_available() else torch.device('cpu')

epochs = 50
bs = 64
lr = 0.1
loss_func = nn.MSELoss()

# use the small data
small_train_dl, small_valid_dl = get_data(small_train_ds, small_valid_ds,
                              bs, GPU=True, GPU_dev=dev)

model, opt = get_model(num_features=2048, lr=lr, GPU=True, GPU_dev=dev)

fit(epochs, model, loss_func, opt, small_train_dl, small_valid_dl)<jupyter_output>epoch:  0 train_loss:  28.1943416595459 valid_loss:  25.76837921142578
epoch:  5 train_loss:  14.899280548095703 valid_loss:  8.270212173461914
epoch:  10 train_loss:  13.446062088012695 valid_loss:  4.915328502655029
epoch:  15 train_loss:  13.446062088012695 valid_loss:  4.915293216705322
epoch:  20 train_loss:  13.446062088012695 valid_loss:  4.915318012237549
epoch:  25 train_loss:  13.446061134338379 valid_loss:  4.915342807769775
epoch:  30 train_loss:  13.446061134338379 valid_loss:  4.9153666496276855
epoch:  35 train_loss:  13.446061134338379 valid_loss:  4.9153900146484375
epoch:  40 train_loss:  13.446061134338379 valid_loss:  4.915412425994873
epoch:  45 train_loss:  13.446061134338379 valid_loss:  4.91543436050415<jupyter_text>Why does that train loss converge, but not go to 0!!! Let's try with one sample!<jupyter_code>one_sample_train_ds = TensorDataset(X_train[0,:].reshape(1,-1), y_train[[0]])
one_sample_valid_ds = TensorDataset(X_valid[0,:].reshape(1,-1), y_valid[[0]])
# Train the model with just one sample!

dev = torch.device('cuda') if torch.cuda.is_available() else torch.device('cpu')

epochs = 50
bs = 64
lr = 0.1
loss_func = nn.MSELoss()

# use the small data
small_train_dl, small_valid_dl = get_data(one_sample_train_ds, one_sample_valid_ds,
                              bs, GPU=True, GPU_dev=dev)

model, opt = get_model(num_features=2048, lr=lr, GPU=True, GPU_dev=dev)

fit(epochs, model, loss_func, opt, small_train_dl, small_valid_dl)<jupyter_output>epoch:  0 train_loss:  38.313716888427734 valid_loss:  5.290582180023193
epoch:  5 train_loss:  108.30901336669922 valid_loss:  41.19266128540039
epoch:  10 train_loss:  5091.599609375 valid_loss:  4538.103515625
epoch:  15 train_loss:  546.7064819335938 valid_loss:  376.0401306152344
epoch:  20 train_loss:  58.7021484375 valid_loss:  13.481624603271484
epoch:  25 train_loss:  6.303093433380127 valid_loss:  2.1886346340179443
epoch:  30 train_loss:  0.6767898201942444 valid_loss:  10.031963348388672
epoch:  35 train_loss:  0.07266971468925476 valid_loss:  13.841573715209961
epoch:  40 train_loss:  0.007802910171449184 valid_loss:  15.222994804382324
epoch:  45 train_loss:  0.0008378112688660622 valid_loss:  15.68995475769043<jupyter_text>Well at least that overfit correctly Let's see If I can get anywhere better using the fastai tabular learner<jupyter_code># Note that this isn't quite a fair comparison because I am not using 
# the same train/valid/test split! But I am just trying to get something working
data = (TabularList.from_df(pka_df, path='.', 
                            cont_names=list(pka_df.columns)[1:])
        .split_by_idx(valid_idx=range(3000,pka_df.shape[0]))
        .label_from_df(cols='pka')
        .databunch()
       )
learn = tabular_learner(data,
                        layers=[256,256],
                        loss_func=nn.MSELoss(),
                        metrics=[root_mean_squared_error, r2_score])
learn.lr_find()
learn.recorder.plot()
learn.fit_one_cycle(5, 1e-02)<jupyter_output><empty_output><jupyter_text>Wow that looks like clear overfitting since the validation loss shoots way up! But the training error is still not even good... There is only one main difference I can think of; in the Deepchem models I had standardized the pka values before training the neural net. Let's see if that helps at all.<jupyter_code># cheating because using all of data for standardization,
# but just want something to work!
pka_df['std_pka'] = (pka_df['pka'] - pka_df['pka'].mean()) / pka_df['pka'].std()
pka_df.head()
# Note that this isn't quite a fair comparison because I am not using 
# the same train/valid/test split! But I am just trying to get something working
data = (TabularList.from_df(pka_df, path='.',
                            cont_names=list(pka_df.columns)[1:-1])
        .split_by_idx(valid_idx=range(3000,pka_df.shape[0]))
        .label_from_df(cols='std_pka')
        .databunch()
       )
data.train_ds.y
learn = tabular_learner(data,
                        layers=[256,256],
                        loss_func=nn.MSELoss(),
                        metrics=[root_mean_squared_error, r2_score])
learn.lr_find()
learn.recorder.plot()
learn.fit_one_cycle(20, 1e-03)<jupyter_output><empty_output><jupyter_text>See if dropout helps<jupyter_code>learn = tabular_learner(data,
                        layers=[256,256],
                        loss_func=nn.MSELoss(),
                        ps=0.5,
                        metrics=[root_mean_squared_error, r2_score])
learn.lr_find()
learn.recorder.plot()
learn.fit_one_cycle(50, 1e-03)<jupyter_output><empty_output><jupyter_text>Try more layers:<jupyter_code>learn = tabular_learner(data,
                        layers=[1000,500],
                        loss_func=nn.MSELoss(),
                        ps=0.5,
                        metrics=[root_mean_squared_error, r2_score])
learn.lr_find()
learn.recorder.plot()
learn.fit_one_cycle(10, 1e-02)
learn = tabular_learner(data,
                        layers=[256,256,256],
                        loss_func=nn.MSELoss(),
                        ps=0.5,
                        metrics=[root_mean_squared_error, r2_score])
learn.lr_find()
learn.recorder.plot()
learn.fit_one_cycle(5, 1e-02)<jupyter_output><empty_output><jupyter_text>I can't even get the coefficient of determination above 0... that is bad. Time to go back to no standardization and try different architectures:<jupyter_code># Note that this isn't quite a fair comparison because I am not using 
# the same train/valid/test split! But I am just trying to get something working
data = (TabularList.from_df(pka_df, path='.',
                            cont_names=list(pka_df.columns)[1:-1])
        .split_by_idx(valid_idx=range(3000,pka_df.shape[0]))
        .label_from_df(cols='pka')
        .databunch()
       ) # no standardization anymore
learn = tabular_learner(data,
                        layers=[256,256],
                        loss_func=nn.MSELoss(),
                        ps=0.5,
                        metrics=[root_mean_squared_error, r2_score])
learn.lr_find()
learn.recorder.plot()
learn.fit_one_cycle(5, 1e-02)<jupyter_output><empty_output><jupyter_text>Again seems we are overfitting, but the train loss is not even low<jupyter_code>learn = tabular_learner(data,
                        layers=[256],
                        loss_func=nn.MSELoss(),
                        ps=[0.5],
                        metrics=[root_mean_squared_error, r2_score])
learn.lr_find()
learn.recorder.plot()
learn.fit_one_cycle(5, 1e-02)<jupyter_output><empty_output><jupyter_text>Since nothing is working right now..., let's try one last (kind of crazy) approach<jupyter_code>fps = [np.array(x) for x in pka_df.values[:,1:-1]]
fps[0]
# pad the pictures a bit so we can make them square
fps = [np.hstack((x,np.zeros(68))) for x in fps]
fps[0]
fps[0].shape
plt.imshow(fps[0].reshape((46,46)), cmap='gray')
!mkdir ./img_train/
import scipy
import imageio
from fastai.basics import *

for idx, fp in enumerate(fps):
  imageio.imwrite('./img_train/'+ str(idx)+'.png',255*fps[idx].reshape((46,46)).astype(np.uint8))
labels_df = pka_df[['pka']]
labels_df.rename(columns={'pka':'label'},inplace=True)
labels_df['name'] = ''
for idx, fp in enumerate(fps):
  labels_df.loc[idx,'name'] = './img_train/'+ str(idx)+'.png'
labels_df.head()
labels_df.to_csv('img_labels.csv')
from fastai.vision import *
from fastai.basics import *
data = (ImageItemList.from_csv('.', 'img_labels.csv', folder='img_train', suffix='.png')
        .random_split_by_pct()
        .label_from_df()
        .databunch())
data.show_batch(rows=2)
learn = create_cnn(data,
                   models.resnet34,
                   loss_func=nn.MSELoss(),
                   metrics=[root_mean_squared_error, r2_score] )
learn.lr_find()
learn.recorder.plot()
lr = 1e-2
learn.fit_one_cycle(5,slice(lr))
learn.unfreeze()
learn.lr_find()
learn.recorder.plot()
learn.fit_one_cycle(5, slice(1e-5, lr/5))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11043380481742321, 'fraction_numerical': 0.057584186541989706, 'mean_word_length': 3.5893536121673004, 'pattern_counts': {'":': 0, '<': 54, '<?xml version=': 0, '>': 54, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '74894', 'license': 'MIT', 'n_tokens_mistral': 6600, 'n_tokens_neox': 5974, 'n_words': 1382, 'path': 'pytorch_and_fastai_molecular_tabular_learner.ipynb', 'repo_name': 'mc-robinson/random_tutorials', 'size': '16897'}