Commit 8623479b authored by Yutong Zhao's avatar Yutong Zhao
Browse files

Add unittest for gradient and ANIRegression

parent 46602deb
Loading
Loading
Loading
Loading
+102 −11
Original line number Diff line number Diff line
@@ -136,7 +136,26 @@ class ANIRegression(TensorGraph):
  def build_grad(self):
    self.grad = tf.gradients(self.outputs, self.atom_feats)

  def compute_grad(self, dataset, batch_size=1):
  def compute_grad(self, dataset, upper_lim=1):
    """
    Computes a batched gradients given an input dataset.

    Parameters
    ----------
    dataset: dc.Dataset
      dataset-like object whose X values will be used to compute
      gradients from
    upper_lim: int
      subset of dataset used.

    Returns
    -------
    np.array
      Gradients of the input of shape (max_atoms, 4). Note that it is up to
      the end user to slice this matrix into the correct shape, since it's very
      likely the derivatives with respect to the atomic numbers are zero.

    """
    with self._get_tf("Graph").as_default():
      if not self.built:
        self.build()
@@ -145,41 +164,110 @@ class ANIRegression(TensorGraph):

      feed_dict = dict()
      X = dataset.X
      flags = np.sign(np.array(X[:batch_size, :, 0]))
      flags = np.sign(np.array(X[:upper_lim, :, 0]))
      feed_dict[self.atom_flags] = np.stack([flags]*self.max_atoms, axis=2)*\
          np.stack([flags]*self.max_atoms, axis=1)
      feed_dict[self.atom_numbers] = np.array(X[:batch_size, :, 0], dtype=int)
      feed_dict[self.atom_feats] = np.array(X[:batch_size, :, :], dtype=float)
      feed_dict[self.atom_numbers] = np.array(X[:upper_lim, :, 0], dtype=int)
      feed_dict[self.atom_feats] = np.array(X[:upper_lim, :, :], dtype=float)
      return self.session.run([self.grad], feed_dict=feed_dict)

  def pred_one(self, X, atomic_nums):
    """
    Makes an energy prediction for a set of atomic coordinates.

    Parameters
    ----------
    X: np.array
      numpy array of shape (a, 3) where a <= max_atoms and
      dtype is float-like
    atomic_nums: np.array
      numpy array of shape (a,) where a is the same as that of X. 

    Returns
    -------
    float
      Predicted energy. Note that the meaning of the returned value is
      dependent on the training y-values both in semantics (relative vs absolute)
      and units (kcal/mol vs Hartrees)

    """
    num_atoms = atomic_nums.shape[0]
    X = X.reshape((num_atoms, 3))
    A = atomic_nums.reshape((atomic_nums.shape[0], 1))
    Z = np.zeros((23, 4))
    Z = np.zeros((self.max_atoms, 4))
    Z[:X.shape[0], 1:X.shape[1]+1] = X
    Z[:A.shape[0], :A.shape[1]] = A
    X = Z
    dd = dc.data.NumpyDataset(np.array(X).reshape((1, 23, 4)), np.array(0), np.array(1))
    dd = dc.data.NumpyDataset(np.array(X).reshape((1, self.max_atoms, 4)), np.array(0), np.array(1))
    return self.predict(dd)[0]

  def grad_one(self, X, atomic_nums):
    y = self.pred_one(X, atomic_nums)
    """
    Computes gradients for that of a single structure.

    Parameters
    ----------
    X: np.array
      numpy array of shape (a, 3) where a <= max_atoms and
      dtype is float-like
    atomic_nums: np.array
      numpy array of shape (a,) where a is the same as that of X. 

    Returns
    -------
    np.array
      derivatives of the same shape and type as input parameter X.

    """
    num_atoms = atomic_nums.shape[0]
    X = X.reshape((num_atoms, 3))
    A = atomic_nums.reshape((atomic_nums.shape[0], 1))
    Z = np.zeros((23, 4))
    Z = np.zeros((self.max_atoms, 4))
    Z[:X.shape[0], 1:X.shape[1]+1] = X
    Z[:A.shape[0], :A.shape[1]] = A
    X = Z
    inp = np.array(X).reshape((1, 23, 4))
    dd = dc.data.NumpyDataset(inp, np.array([y]), np.array([1]))
    inp = np.array(X).reshape((1, self.max_atoms, 4))
    dd = dc.data.NumpyDataset(inp, np.array([1]), np.array([1]))
    res = self.compute_grad(dd)[0][0][0]
    res = res[:num_atoms, 1:]
    return res.reshape((num_atoms*3,))

  def minimize_structure(self, X, atomic_nums):
    """
    Minimizes a structure, as defined by a set of coordinates and their atomic
    numbers.

    Parameters
    ----------
    X: np.array
      numpy array of shape (a, 3) where a <= max_atoms and
      dtype is float-like
    atomic_nums: np.array
      numpy array of shape (a,) where a is the same as that of X. 

    Returns
    -------
    np.array
      minimized coordinates of the same shape and type as input parameter X.

    """
    num_atoms = atomic_nums.shape[0]

    # minimizer_kwargs = {
    #   "jac": self.grad_one,
    #   "args": (atomic_nums,),
    #   "method": "CG",
    #   "options": {"disp": True},
    # }

    # res = scipy.optimize.basinhopping(
    #   self.pred_one,
    #   X,
    #   niter=75,
    #   T=0.015, # ~10kcal/mol
    #   stepsize=0.17,
    #   minimizer_kwargs=minimizer_kwargs)

    res = scipy.optimize.minimize(
      self.pred_one,
      X,
@@ -188,6 +276,7 @@ class ANIRegression(TensorGraph):
      method="BFGS",
      tol=1e-6,
      options={'disp': True})

    return res.x.reshape((num_atoms, 3))

  def build_graph(self):
@@ -196,7 +285,9 @@ class ANIRegression(TensorGraph):
    self.atom_flags = Feature(shape=(None, self.max_atoms, self.max_atoms))
    self.atom_feats = Feature(shape=(None, self.max_atoms, 4))

    previous_layer = ANIFeat(self.atom_feats)
    previous_layer = ANIFeat(
      in_layers=self.atom_feats,
      max_atoms=self.max_atoms)

    self.featurized = previous_layer

+66 −0
Original line number Diff line number Diff line
import scipy
import numpy as np
import unittest

from deepchem.models import ANIRegression
import deepchem as dc

class TestANIRegression(unittest.TestCase):

  def test_gradients(self):

    max_atoms = 3

    X = np.array([
      [1, 5.0, 3.2, 1.1],
      [6, 1.0, 3.4, -1.1],
      [1, 2.3, 3.4, 2.2]
    ])

    X = X.reshape((1, X.shape[0], X.shape[1]))

    y = np.array([2.0])
    y = y.reshape((1,1))

    layer_structures = [128, 128, 64]
    atom_number_cases = [1, 6, 7, 8]

    model = ANIRegression(
      1,
      max_atoms,
      layer_structures=layer_structures,
      atom_number_cases=atom_number_cases,
      batch_size=1,
      learning_rate=0.001,
      use_queue=False,
      mode="regression")

    print(X.shape, y.shape)

    train_dataset = dc.data.NumpyDataset(X, y, n_tasks=1)

    model.fit(train_dataset, nb_epoch=2, checkpoint_interval=100)

    new_x = np.array([
      -2.0, 1.2, 2.1,
      1.3, -6.4, 3.1,
      -2.5, 2.4, 5.6,
    ])

    new_atomic_nums = np.array([1,1,6])

    grad_approx = scipy.optimize.approx_fprime(
      new_x,
      model.pred_one,
      1e-4,
      new_atomic_nums)

    grad_exact = model.grad_one(new_x, new_atomic_nums)

    print(grad_approx)
    print(grad_exact)

    np.testing.assert_array_almost_equal(grad_approx, grad_exact)

if __name__ == '__main__':
  unittest.main()
 No newline at end of file
+1 −0
Original line number Diff line number Diff line
import numpy as np
from flask import request, abort, Flask

import flask
+41 −38
Original line number Diff line number Diff line
@@ -27,7 +27,7 @@ def load_roiterberg_ANI(relative=True):
  hdf5files = [
      'ani_gdb_s01.h5',
      'ani_gdb_s02.h5',
      # 'ani_gdb_s03.h5',
      'ani_gdb_s03.h5',
      # 'ani_gdb_s04.h5',
      # 'ani_gdb_s05.h5',
      # 'ani_gdb_s06.h5',
@@ -138,14 +138,31 @@ def broadcast(dataset, metadata):

if __name__ == "__main__":
  

  max_atoms = 23
  batch_size = 64  # CHANGED FROM 16
  layer_structures = [128, 128, 64]
  atom_number_cases = [1, 6, 7, 8]

  metric = [
      dc.metrics.Metric(dc.metrics.mean_absolute_error, mode="regression"),
      dc.metrics.Metric(dc.metrics.pearson_r2_score, mode="regression")
  ]

  model_dir = "/tmp/ani.pkl"

  if os.path.exists(model_dir):
    print("Restoring existing model...")
    model = dc.models.ANIRegression.load_from_dir(model_dir=model_dir)
  else:
    print("Fitting new model...")

    train_valid_dataset, test_dataset, all_groups = load_roiterberg_ANI(relative=True)

    splitter = dc.splits.RandomGroupSplitter(broadcast(train_valid_dataset, all_groups))

    print("performing 1-fold split...")

  # n_folds = 5

    # for fold in range(n_folds):
    print("Folding once....")
    train_dataset, valid_dataset = splitter.train_test_split(train_valid_dataset)
@@ -160,23 +177,6 @@ if __name__ == "__main__":
      valid_dataset = transformer.transform(valid_dataset)
      test_dataset = transformer.transform(test_dataset)

  max_atoms = 23
  batch_size = 64  # CHANGED FROM 16
  layer_structures = [128, 128, 64]
  atom_number_cases = [1, 6, 7, 8]

  metric = [
      dc.metrics.Metric(dc.metrics.mean_absolute_error, mode="regression"),
      dc.metrics.Metric(dc.metrics.pearson_r2_score, mode="regression")
  ]

  model_dir = "/tmp/ani.pkl"

  if os.path.exists(model_dir):
    print("Restoring existing model...")
    model = dc.models.ANIRegression.load_from_dir(model_dir=model_dir)
  else:
    print("Fitting new model...")
    model = dc.models.ANIRegression(
        1,
        max_atoms,
@@ -189,10 +189,12 @@ if __name__ == "__main__":
        mode="regression")

    # For production, set nb_epoch to 100+
    model.fit(train_dataset, nb_epoch=1, checkpoint_interval=100)
    for i in range(10):
      model.fit(train_dataset, nb_epoch=10, checkpoint_interval=100)

      print("Saving model...")
      model.save()
      print("Done.")

    print("Evaluating model")
    train_scores = model.evaluate(train_dataset, metric, transformers)
@@ -216,13 +218,14 @@ if __name__ == "__main__":

  atomic_nums = np.array([1, 8, 1])

  print("Prediction of a single test set element:")
  print("Prediction of a single test set structure:")
  print(model.pred_one(coords, atomic_nums))

  print("Gradient of a single test set element:")
  print("Gradient of a single test set structure:")
  print(model.grad_one(coords, atomic_nums))


  # print("Minimization of a single test set structure:")
  # print(model.minimize_structure(coords, atomic_nums))

  app.webapp.model = model
  app.webapp.run(host='0.0.0.0', debug=False)