import numpy as np
import pandas as pd
import scipy.linalg as linalg
import rdkit.Chem as Chem
from rdkit.Chem import AllChem
import torch
import matplotlib.pyplot as plt
def generate_graph(smile): #block copolymer code
def get_atom_hash(atomic_number):
""" A helper function to quickly encode atomic number into a one-hot vector """
atomic_number_list = [6, 7, 8, 9, 15, 16, 17, 35, 53]
if atomic_number in atomic_number_list:
return atomic_number_list.index(atomic_number)
else:
return len(atomic_number_list)
def encode_atom(atom):
""" Generates the vector of features for an atom in a molecule"""
atom_matrix = np.zeros(10)
atom_matrix[get_atom_hash(atom.GetAtomicNum())] = 1 #only information we are getting is atomic # in atom feature matrix-- 1 hot encoding
return atom_matrix
m = AllChem.MolFromSmiles(smile)
m = AllChem.AddHs(m)
AllChem.ComputeGasteigerCharges(m)
atom_matrix = np.array(list(map(encode_atom, (atom for atom in m.GetAtoms()))))
adj_matrix = AllChem.GetAdjacencyMatrix(m) + np.identity(m.GetNumAtoms()) # Augmented Adjacency Matrix #can generate adjacency matrix easily from RdKit
return(atom_matrix, adj_matrix, m.GetNumAtoms())
esol_data = pd.read_csv('delaney-processed.csv', sep=',')
amigdalin = esol_data['smiles'][0]
print(amigdalin)
print(generate_graph(amigdalin)[2])
OCC3OC(OCC2OC(OC(C#N)c1ccccc1)C(O)C(O)C2O)C(O)C(O)C3O
59
n_samples = 100 #change to 100
feature_matrices = []
adj_matrices = [] #block adj matrix
molecule_membership = [] #keeps track of which atoms respond to molecule 1, 2, etc.
for i in range(n_samples):
feature_matrix, adj_matrix, n_atoms = generate_graph(esol_data['smiles'][i])
feature_matrices.append(feature_matrix)
adj_matrices.append(adj_matrix)
molecule_membership = np.concatenate((molecule_membership, [i]*n_atoms))
big_feature_matrix = np.concatenate(feature_matrices, axis = 0)
print(np.shape(big_feature_matrix))
big_adj_matrix = linalg.block_diag(*adj_matrices)
print(np.shape(big_adj_matrix))
# need to aggregate atom-level features into something for molecule: (num atoms by num features) * mask --> transform into matrix (num_mols x num_features)
# sum only atoms that belong to that particular molecule
masks = torch.tensor(np.array([np.where(molecule_membership == i, 1, 0) for i in range(n_samples)])).double()
print(np.shape(masks))
(2616, 10)
(2616, 2616)
torch.Size([100, 2616])
### Calculate normalized graph laplacian D^-1/2AD^-1/2
degree_vector = np.sum(big_adj_matrix, axis=1)
degree_vector = np.sqrt(1/degree_vector)
D = torch.tensor(np.diag(degree_vector))
A = torch.tensor(big_adj_matrix)
norm_graph_laplacian = D @ A @ D
print(np.shape(norm_graph_laplacian))
#input size is # feature
class SolNet(torch.nn.Module):
def __init__(self, laplacian, molecule_masks, input_size = 10, internal_size = 10, output_size = 1, num_conv_layers = 3, activation_function=torch.nn.functional.relu):
super().__init__()
### Your Code Here ###
#GCN layer h = nonlinearity(D^-1/2AD^-1/2 *h*Weight)
self.W0 = torch.nn.Parameter(torch.randn((input_size, internal_size)))
self.W1 = torch.nn.Parameter(torch.randn((internal_size, internal_size)))
self.W2 = torch.nn.Parameter(torch.randn((internal_size, internal_size)))
self.W3 = torch.nn.Parameter(torch.randn((internal_size, output_size)))
self.laplacian = laplacian
self.molecule_masks = molecule_masks
self.activation_function = activation_function
def forward(self, x):
"""Computes the prediction of the neural network"""
### Your Code Here ###
x = self.activation_function(self.laplacian @ x @ self.W0)
x = self.activation_function(self.laplacian @ x @ self.W1)
x = self.activation_function(self.laplacian @ x @ self.W2)
x = self.molecule_masks @ x #Aggregation of atom level features in each molecule
y = x @ self.W3
return y
torch.Size([2616, 2616])
## Create Targets
targets = torch.tensor(esol_data['measured log solubility in mols per litre'][:n_samples])
# Only train on 2/3 of points
labelled_indices = np.random.choice(np.arange(n_samples), int(2*n_samples/3), replace = False)
unlabelled_indices = np.setdiff1d(np.arange(n_samples), labelled_indices)
### Train your model here (You can use a very similar approach to the one used in previous practicals) ###
model = SolNet(norm_graph_laplacian, masks)
model.double()
num_epochs=500
learning_rate=1E-2
opt = torch.optim.Adam(model.parameters(), lr=learning_rate)
all_train_loss = []
all_test_loss = []
for epoch in range(num_epochs):
pred = model(big_feature_matrix).squeeze()
train_loss = torch.mean((pred[labelled_indices]-targets[labelled_indices])**2) # sum over last two dimensions and take average over batch
opt.zero_grad()
train_loss.backward()
opt.step()
all_train_loss.append(train_loss)
test_loss = torch.mean((pred[unlabelled_indices]-targets[unlabelled_indices])**2)
all_test_loss.append(test_loss)
plt.plot(all_train_loss, label="Train")
plt.plot(all_test_loss, label="Test")
plt.title('Epoch Loss')
plt.xlabel('Number of Epochs')
plt.ylabel('Loss')
plt.legend()
print(all_train_loss[-1], all_test_loss[-1])
tensor(1.8123, dtype=torch.float64, grad_fn=<MeanBackward0>) tensor(1.6626, dtype=torch.float64, grad_fn=<MeanBackward0>)
# relies on local information (graph kernel vs whole aggregation graph)
# >3 layers --> >3rd order NN degrades performance
class MPNet(torch.nn.Module):
def __init__(self, n_features, activation_function=torch.nn.functional.relu):
super().__init__()
self.W0 = torch.nn.Parameter(0.1*torch.randn((n_features*2, n_features)))
self.W1 = torch.nn.Parameter(0.1*torch.randn((n_features*2, n_features))) #current vector and the message
self.W2 = torch.nn.Parameter(0.1*torch.randn((n_features, 1))) #linear layer down to output size of 1 to predict solubilities
self.n_features = n_features
self.activation_function = activation_function
def forward(self, feats, adj):
"""Computes the prediction of the neural network"""
messages = torch.stack([sum(self.activation_function(torch.hstack((feats[i], feats[j])) @ self.W0) for j in np.where(adj[i]==1)[0]) for i in range(len(feats))])
feats = self.activation_function(torch.hstack((feats, messages)) @ self.W1)
y = torch.sum(feats, axis=0) @ self.W2
# for i in range(len(feats)):
# message = 0
# for j in np.where(adj[i]==1)[0]:
# message += self.activation_function(torch.hstack((feats[i], feats[i])) @ self.W0)
# feats[i] = self.activation_function(torch.hstack((feats[i]), message)) @ self.W1) / len(np.where(adj[i]==1)[0])
#y = torch.sum(feats, axis=0) @ self.W2
return y
### Train your model here (You can use a very similar approach to the one used in previous practicals) ###
# Note that training here can be slow if your module is not sufficiently vectorized. It might
# be helpful to reduce n_samples to make sure your architecture is working properly before trying to apply
# to larger data.
#loop through labeled indices, at end of epoch
model = MPNet(10)
model.double()
num_epochs=100
learning_rate=1E-2
opt = torch.optim.Adam(model.parameters(), lr=learning_rate)
all_train_loss = []
all_test_loss = []
for epoch in range(num_epochs):
epoch_loss = 0
for i in labelled_indices:
pred = model(torch.tensor(generate_graph(esol_data['smiles'][i])[0]), torch.tensor(generate_graph(esol_data['smiles'][i])[1]))
loss = torch.mean((pred-targets[i])**2)
epoch_loss += loss
opt.zero_grad()
epoch_loss.backward()
opt.step()
all_train_loss.append(epoch_loss/len(labelled_indices))
print(epoch_loss/len(labelled_indices))
epoch_test_loss=0
for j in unlabelled_indices:
with torch.no_grad():
pred = model(torch.tensor(generate_graph(esol_data['smiles'][j])[0]), torch.tensor(generate_graph(esol_data['smiles'][j])[1]))
loss = torch.mean((pred-targets[j])**2)
epoch_test_loss += loss
all_test_loss.append(epoch_test_loss/len(unlabelled_indices))
print(epoch_test_loss/len(unlabelled_indices))
plt.plot(all_train_loss, label="Train")
plt.plot(all_test_loss, label="Test")
plt.title('Epoch Loss')
plt.xlabel('Number of Epochs')
plt.ylabel('Loss')
plt.legend()
print(all_train_loss[-1], all_test_loss[-1])
tensor(15.7882, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(17.0197, dtype=torch.float64)
tensor(13.4477, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(14.2287, dtype=torch.float64)
tensor(11.2771, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(11.4754, dtype=torch.float64)
tensor(9.1945, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(8.6294, dtype=torch.float64)
tensor(7.1567, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(5.8892, dtype=torch.float64)
tensor(5.4452, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(3.9241, dtype=torch.float64)
tensor(4.7026, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(3.7421, dtype=torch.float64)
tensor(5.6001, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(4.0735, dtype=torch.float64)
tensor(6.2468, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(3.7523, dtype=torch.float64)
tensor(5.7438, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(3.2917, dtype=torch.float64)
tensor(4.8739, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(3.1157, dtype=torch.float64)
tensor(4.1880, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(3.2648, dtype=torch.float64)
tensor(3.8541, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(3.5449, dtype=torch.float64)
tensor(3.7635, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(3.8088, dtype=torch.float64)
tensor(3.7573, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(3.9666, dtype=torch.float64)
tensor(3.7369, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(3.9715, dtype=torch.float64)
tensor(3.6500, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(3.8175, dtype=torch.float64)
tensor(3.4820, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(3.5284, dtype=torch.float64)
tensor(3.2456, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(3.1477, dtype=torch.float64)
tensor(2.9699, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(2.7229, dtype=torch.float64)
tensor(2.6879, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(2.3205, dtype=torch.float64)
tensor(2.4508, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.9991, dtype=torch.float64)
tensor(2.3017, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.7970, dtype=torch.float64)
tensor(2.2530, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.7025, dtype=torch.float64)
tensor(2.2478, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.6380, dtype=torch.float64)
tensor(2.1805, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.5708, dtype=torch.float64)
tensor(2.0313, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.5197, dtype=torch.float64)
tensor(1.8543, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.5172, dtype=torch.float64)
tensor(1.6903, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.5719, dtype=torch.float64)
tensor(1.5808, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.6677, dtype=torch.float64)
tensor(1.5209, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.7709, dtype=torch.float64)
tensor(1.4931, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.8461, dtype=torch.float64)
tensor(1.4745, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.8785, dtype=torch.float64)
tensor(1.4555, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.8335, dtype=torch.float64)
tensor(1.4242, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.7401, dtype=torch.float64)
tensor(1.3939, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.6310, dtype=torch.float64)
tensor(1.3815, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.5380, dtype=torch.float64)
tensor(1.3942, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.4920, dtype=torch.float64)
tensor(1.4170, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.4914, dtype=torch.float64)
tensor(1.4225, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.5201, dtype=torch.float64)
tensor(1.4051, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.5767, dtype=torch.float64)
tensor(1.3771, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.6483, dtype=torch.float64)
tensor(1.3521, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.7202, dtype=torch.float64)
tensor(1.3346, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.7702, dtype=torch.float64)
tensor(1.3211, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.7827, dtype=torch.float64)
tensor(1.3056, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.7547, dtype=torch.float64)
tensor(1.2858, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.6969, dtype=torch.float64)
tensor(1.2646, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.6274, dtype=torch.float64)
tensor(1.2482, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.5644, dtype=torch.float64)
tensor(1.2403, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.5187, dtype=torch.float64)
tensor(1.2389, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.4947, dtype=torch.float64)
tensor(1.2384, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.4940, dtype=torch.float64)
tensor(1.2345, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.5142, dtype=torch.float64)
tensor(1.2270, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.5514, dtype=torch.float64)
tensor(1.2193, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.5981, dtype=torch.float64)
tensor(1.2148, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.6439, dtype=torch.float64)
tensor(1.2136, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.6773, dtype=torch.float64)
tensor(1.2133, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.6908, dtype=torch.float64)
tensor(1.2110, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.6829, dtype=torch.float64)
tensor(1.2059, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.6590, dtype=torch.float64)
tensor(1.1995, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.6283, dtype=torch.float64)
tensor(1.1940, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.6009, dtype=torch.float64)
tensor(1.1904, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.5841, dtype=torch.float64)
tensor(1.1878, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.5817, dtype=torch.float64)
tensor(1.1845, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.5941, dtype=torch.float64)
tensor(1.1797, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.6194, dtype=torch.float64)
tensor(1.1743, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.6521, dtype=torch.float64)
tensor(1.1697, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.6843, dtype=torch.float64)
tensor(1.1667, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.7090, dtype=torch.float64)
tensor(1.1644, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.7209, dtype=torch.float64)
tensor(1.1618, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.7188, dtype=torch.float64)
tensor(1.1586, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.7063, dtype=torch.float64)
tensor(1.1552, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.6896, dtype=torch.float64)
tensor(1.1524, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.6757, dtype=torch.float64)
tensor(1.1504, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.6695, dtype=torch.float64)
tensor(1.1486, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.6734, dtype=torch.float64)
tensor(1.1463, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.6866, dtype=torch.float64)
tensor(1.1435, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.7059, dtype=torch.float64)
tensor(1.1407, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.7259, dtype=torch.float64)
tensor(1.1383, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.7412, dtype=torch.float64)
tensor(1.1361, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.7478, dtype=torch.float64)
tensor(1.1337, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.7450, dtype=torch.float64)
tensor(1.1310, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.7354, dtype=torch.float64)
tensor(1.1281, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.7235, dtype=torch.float64)
tensor(1.1254, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.7142, dtype=torch.float64)
tensor(1.1230, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.7107, dtype=torch.float64)
tensor(1.1206, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.7143, dtype=torch.float64)
tensor(1.1181, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.7237, dtype=torch.float64)
tensor(1.1155, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.7361, dtype=torch.float64)
tensor(1.1130, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.7474, dtype=torch.float64)
tensor(1.1108, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.7530, dtype=torch.float64)
tensor(1.1087, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.7534, dtype=torch.float64)
tensor(1.1066, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.7497, dtype=torch.float64)
tensor(1.1044, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.7445, dtype=torch.float64)
tensor(1.1025, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.7406, dtype=torch.float64)
tensor(1.1009, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.7400, dtype=torch.float64)
tensor(1.0993, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.7396, dtype=torch.float64)
tensor(1.0980, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.7436, dtype=torch.float64)
tensor(1.0966, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.7509, dtype=torch.float64)
tensor(1.0952, dtype=torch.float64, grad_fn=<DivBackward0>)
tensor(1.7593, dtype=torch.float64)
tensor(1.0952, dtype=torch.float64, grad_fn=<DivBackward0>) tensor(1.7593, dtype=torch.float64)