Factored out evaluation code into separate function.

from deep_chem.utils.evaluate import compute_r2_scores

from deep_chem.utils.evaluate import compute_rms_scores

from deep_chem.utils.evaluate import compute_roc_auc_scores

from deep_chem.utils.load import load_and_transform_dataset

def fit_multitask_mlp(train_data, test_data, task_types, **training_params):

def fit_multitask_mlp(per_task_data, task_types, **training_params):

  """

  Perform stochastic gradient descent optimization for a keras multitask MLP.

  Returns AUCs, R^2 scores, and RMS values.

  training_params: dict

    Aggregates keyword parameters to pass to train_multitask_model

  """

  models = {}

  # Follows convention from process_datasets that the data for multitask models

  # is grouped under key "all"

  (train, X_train, y_train, W_train), (test, X_test, y_test, W_test) = (

      train_data, test_data)

  model = train_multitask_model(X_train, y_train, W_train, task_types,

      per_task_data["all"])

  models["all"] = train_multitask_model(X_train, y_train, W_train, task_types,

                                **training_params)

  results = eval_model(test, model, task_types,

      modeltype="keras_multitask")

  local_task_types = task_types.copy()

  aucs = compute_roc_auc_scores(results, local_task_types)

  if aucs:

    print "Mean AUC: %f" % np.mean(np.array(aucs.values()))

  r2s = compute_r2_scores(results, local_task_types)

  if r2s:

    print "Mean R^2: %f" % np.mean(np.array(r2s.values()))

  return results

  return models

def fit_singletask_mlp(per_task_data, task_types, num_to_train=None, **training_params):

def fit_singletask_mlp(per_task_data, task_types, **training_params):

  """

  Perform stochastic gradient descent optimization for a keras MLP.

  training_params: dict

    Aggregates keyword parameters to pass to train_multitask_model

  """

  ret_vals = {}

  aucs, r2s, rms = {}, {}, {}

  sorted_targets = sorted(per_task_data.keys())

  if num_to_train:

    sorted_targets = sorted_targets[:num_to_train]

  all_results = {}

  for index, target in enumerate(sorted_targets):

  models = {}

  for index, target in enumerate(sorted(per_task_data.keys())):

    print "Training model %d" % index

    print "Target %s" % target

    (train, X_train, y_train, W_train), (test, X_test, y_test, W_test) = (

        per_task_data[target])

    print "len(train)"

    print len(train)

    print "len(test)"

    print len(test)

    model = train_multitask_model(X_train, y_train, W_train,

    print "%d compounds in Train" % len(train)

    print "%d compounds in Test" % len(test)

    models[target] = train_multitask_model(X_train, y_train, W_train,

        {target: task_types[target]}, **training_params)

    results = eval_model(test, model, {target: task_types[target]}, 

                         # We run singletask models as special cases of

                         # multitask.

                         modeltype="keras_multitask")

    all_results[target] = results[target]

    target_aucs = compute_roc_auc_scores(results, task_types)

    target_r2s = compute_r2_scores(results, task_types)

    target_rms = compute_rms_scores(results, task_types)

    aucs.update(target_aucs)

    r2s.update(target_r2s)

    rms.update(target_rms)

  if aucs:

    print aucs

    print "Mean AUC: %f" % np.mean(np.array(aucs.values()))

  if r2s:

    print r2s

    print "Mean R^2: %f" % np.mean(np.array(r2s.values()))

  return all_results

  return models

def train_multitask_model(X, y, W, task_types,

  learning_rate=0.01, decay=1e-6, momentum=0.9, nesterov=True, activation="relu",

from keras.layers.convolutional import Convolution3D, MaxPooling3D

from keras.utils import np_utils

from deep_chem.utils.preprocess import split_dataset

from deep_chem.utils.load import load_and_transform_dataset

from deep_chem.utils.preprocess import tensor_dataset_to_numpy

from deep_chem.utils.evaluate import eval_model

from deep_chem.utils.evaluate import compute_r2_scores

def fit_3D_convolution(train_data, test_data, task_types, axis_length=32, **training_params):

def fit_3D_convolution(per_task_data, task_types, axis_length=32, **training_params):

  """

  Perform stochastic gradient descent for a 3D CNN.

  """

  models = {}

  (X_train, y_train, W_train, train), (X_test, y_test, W_test, test) = (

      train_data, test_data)

      per_task_data["all"] 

  nb_classes = 2

  model = train_3D_convolution(X_train, y_train, axis_length, **training_params)

  results = eval_model(test, model, task_types,

      modeltype="keras", mode="tensor")

  local_task_types = task_types.copy()

  r2s = compute_r2_scores(results, local_task_types)

  print "Mean R^2: %f" % np.mean(np.array(r2s.values()))

  return results

  models["all"] = train_3D_convolution(X_train, y_train, axis_length, **training_params)

  return models

def train_3D_convolution(X, y, axis_length=32, batch_size=50, nb_epoch=1):

  """

                       help="Name of measured endpoint to predict.")

  parser.add_argument("--split-endpoint", type=str, default=None,

                       help="Name of endpoint specifying train/test split.")

  # SUBCOMMAND Training params

  parser.add_argument("--n-hidden", type=int, default=500,

                      help="Number of hidden neurons for NN models.")

  parser.add_argument("--learning-rate", type=float, default=0.01,

                  help="Percent of training data to use for validation.")

  parser.add_argument("--weight-positives", type=bool, default=False,

                  help="Weight positive examples to have same total weight as negatives.")

  # TODO(rbharath): Remove this once debugging is complete.

  parser.add_argument("--num-to-train", type=int, default=None,

                  help="Number of datasets to train on. Only for debug.")

  parser.add_argument("--axis-length", type=int, default=32,

                  help="Size of a grid axis for 3D CNN input.")

  # SUBCOMMAND Evaluation types

  parser.add_argument("--compute-aucs", type=bool, default=False,

                      help="Compute AUC for trained models on test set.")

  parser.add_argument("--compute-r2s", type=bool, default=False,

                      help="Compute R^2 for trained models on test set.")

  parser.add_argument("--compute-rms", type=bool, default=False,

                      help="Compute RMS for trained models on test set.")

  return parser.parse_args(input_args)

def main():

  input_transforms = args.input_transforms 

  output_transforms = {target: args.output_transforms for target in targets}

  # TODO(rbharath): The datatype (vector vs. tensor) should be automatically

  # detected in dataset_to_numpy

  datatype = "tensor" if args.model == "3D_cnn" else "vector"

  processed = process_datasets(paths,

  per_task_data = process_datasets(paths,

      input_transforms, output_transforms, feature_types=args.feature_types, 

      prediction_endpoint=args.prediction_endpoint,

      split_endpoint=args.split_endpoint,

      splittype=args.splittype, weight_positives=args.weight_positives,

      datatype=datatype, mode=args.mode)

  if args.mode == "multitask":

    train_data, test_data = processed

  else:

    per_task_data = processed

  # TODO(rbharath): Bundle training params into a training_param dict that's passed

  # down to these functions.

  if args.model == "singletask_deep_network":

    results = fit_singletask_mlp(per_task_data, task_types, n_hidden=args.n_hidden,

    models = fit_singletask_mlp(per_task_data, task_types, n_hidden=args.n_hidden,

      learning_rate=args.learning_rate, dropout=args.dropout,

      nb_epoch=args.n_epochs, decay=args.decay, batch_size=args.batch_size,

      validation_split=args.validation_split,

      num_to_train=args.num_to_train)

  elif args.model == "multitask_deep_network":

    results = fit_multitask_mlp(train_data, test_data, task_types,

    models = fit_multitask_mlp(per_task_data, task_types,

      n_hidden=args.n_hidden, learning_rate = args.learning_rate,

      dropout = args.dropout, batch_size=args.batch_size,

      nb_epoch=args.n_epochs, decay=args.decay,

      validation_split=args.validation_split)

  elif args.model == "3D_cnn":

    results = fit_3D_convolution(train_data, test_data, task_types,

    models = fit_3D_convolution(train_data, test_data, task_types,

        axis_length=args.axis_length, nb_epoch=args.n_epochs,

        batch_size=args.batch_size)

  else:

    results = fit_singletask_models(per_task_data, args.model, task_types,

    models = fit_singletask_models(per_task_data, args.model, task_types,

                                    num_to_train=args.num_to_train)

  results, aucs, r2s, rms = compute_model_performance(per_task_data, models,

    args.compute_aucs, args.compute_r2s, args.compute_rms) 

  if args.csv_out is not None:

    results_to_csv(results, args.csv_out, task_type=args.task_type)

from rdkit import Chem

from rdkit.Chem.Descriptors import ExactMolWt

def compute_model_performance(per_task_data, models, aucs=True, r2s=False, rms=False):

  """Computes statistics for model performance on test set."""

  all_results, auc_vals, r2s, rms = {}, {}, {}, {}

  for index, target in enumerate(sorted(per_task_data.keys())):

    print "Evaluating model %d" % index

    print "Target %s" % target

    (train, _, _, _), (test, _, _, _) = per_task_data[target]

    model = models[target]

    results = eval_model(test, model, {target: task_types[target]}, 

                         # We run singletask models as special cases of

                         # multitask.

                         modeltype="keras_multitask")

    all_results[target] = results[target]

    if aucs:

      auc_vals.update(compute_roc_auc_scores(results, task_types))

    if r2s: 

      r2_vals.update(compute_r2_scores(results, task_types))

    if rms:

      rms_vals.update(compute_rms_scores(results, task_types))

  if aucs:

    print "Mean AUC: %f" % np.mean(np.array(auc_vals.values()))

  if r2s:

    print "Mean R^2: %f" % np.mean(np.array(r2_vals.values()))

  if rms:

    print "Mean RMS: %f" % np.mean(np.array(rms_vals.values()))

  return all_results, aucs, r2s, rms

def model_predictions(test_set, model, n_targets, task_types,

    modeltype="sklearn", mode="regular"):

    modeltype="sklearn", datatype="vector"):

  """Obtains predictions of provided model on test_set.

  Returns a list of per-task predictions.

    Either sklearn, keras, or keras_multitask

  """

  # Extract features for test set and make preds

  if mode == "regular":

  if datatype == "vector":

    X, _, _ = dataset_to_numpy(test_set)

  elif mode == "tensor":

  elif datatype == "tensor":

    X, _, _ = tensor_dataset_to_numpy(test_set)

    (n_samples, axis_length, _, _, n_channels) = np.shape(X)

    # TODO(rbharath): Modify the featurization so that it matches desired shaped. 

    X = np.reshape(X, (n_samples, axis_length, n_channels, axis_length, axis_length))

  else:

    raise ValueError("Improper mode: " + str(mode))

    raise ValueError("Datatype must be vector or tensor.")

  if modeltype == "keras_multitask":

    predictions = model.predict({"input": X})

    ypreds = []

    ypreds = [ypreds]

  return ypreds

def size_eval_model(test_set, model, task_types, modeltype="sklearn"):

  """Split test set based on size of molecule."""

  weights = {}

  for smiles in test_set:

    weights[smiles] = ExactMolWt(Chem.MolFromSmiles(smiles))

  #print weights

  weight_arr = np.array(weights.values())

  print "mean: " + str(np.mean(weight_arr))

  print "std: " + str(np.std(weight_arr))

  print "max: " + str(np.amax(weight_arr))

  print "min: " + str(np.amin(weight_arr))

  buckets = {250: {}, 500: {}, 750: {}, 1000: {}, 1250: {}, 1500: {}, 1750: {}, 2000: {}, 2250: {}, 2500: {}, 2750: {}}

  buckets_to_labels = {250: "0-250", 500: "250-500", 750: "500-750", 1000: "750-1000", 1250: "1000-1250", 1500: "1250-1500", 1750: "1500-1750", 2000: "1750-2000", 2250: "2000-2250", 2500: "2250-2500", 2750: "2500-2750"}

  for smiles in test_set:

    weight = weights[smiles]

    if weight < 250:

      buckets[250][smiles] = test_set[smiles]

    elif weight < 500:

      buckets[500][smiles] = test_set[smiles]

    elif weight < 750:

      buckets[750][smiles] = test_set[smiles]

    elif weight < 1000:

      buckets[1000][smiles] = test_set[smiles]

    elif weight < 1250:

      buckets[1250][smiles] = test_set[smiles]

    elif weight < 1500:

      buckets[1500][smiles] = test_set[smiles]

    elif weight < 1750:

      buckets[1750][smiles] = test_set[smiles]

    elif weight < 2000:

      buckets[2000][smiles] = test_set[smiles]

    elif weight < 2250:

      buckets[2250][smiles] = test_set[smiles]

    elif weight < 2500:

      buckets[2500][smiles] = test_set[smiles]

    elif weight < 2750:

      buckets[2750][smiles] = test_set[smiles]

    else:

      raise ValueError("High Weight: " + str(weight))

  for weight_class in sorted(buckets.keys()):

    test_bucket = buckets[weight_class]

    if len(test_bucket) == 0:

      continue

    print "Evaluating model for %s dalton molecules" % buckets_to_labels[weight_class]

    print "%d compounds in bucket" % len(test_bucket)

    results = eval_model(test_bucket, model, task_types, modeltype=modeltype)

    target_r2s = compute_r2_scores(results, task_types)

    target_rms = compute_rms_scores(results, task_types)

    print "R^2: " + str(target_r2s)

    print "RMS: " + str(target_rms)

  print "Performing Global Evaluation"

  results = eval_model(test_set, model, task_types, modeltype=modeltype)

  target_r2s = compute_r2_scores(results, task_types)

  target_rms = compute_rms_scores(results, task_types)

  print "R^2: " + str(target_r2s)

  print "RMS: " + str(target_rms)

def eval_model(test_set, model, task_types, modeltype="sklearn", mode="regular"):

def eval_model(test_set, model, task_types, modeltype="sklearn", datatype="vector"):

  """Evaluates the provided model on the test-set.

  Returns a dict which maps target-names to pairs of np.ndarrays (ytrue,

  local_task_types = task_types.copy()

  endpoints = sorted_targets

  ypreds = model_predictions(test_set, model, len(sorted_targets),

      local_task_types, modeltype=modeltype, mode=mode)

      local_task_types, modeltype=modeltype, datatype=datatype)

  results = {}

  for target in endpoints:

    results[target] = ([], [], [])  # (smiles, ytrue, yscore)

  dataset = load_and_transform_dataset(paths, input_transforms, output_transforms,

      prediction_endpoint, split_endpoint=split_endpoint,

      feature_types=feature_types, weight_positives=weight_positives)

  arrays = {}

  if mode == "singletask":

    singletask = multitask_to_singletask(dataset)

    arrays = {}

    for target in singletask:

      data = singletask[target]

      if len(data) == 0:

        continue

      train, test = split_dataset(dataset, splittype)

      train_data, test_data = to_arrays(train, test, datatype)

      arrays[target] = train_data, test_data 

    return arrays

      arrays[target] = (train_data, test_data)

  elif mode == "multitask":

    sorted_targets = sorted(dataset.keys())

    train, test = split_dataset(dataset, splittype)

    train_data, test_data = to_arrays(train, test, datatype)

    return train_data, test_data

    arrays["all"] = (train_data, test_data)

  else:

    raise ValueError("Unsupported mode for process_datasets.")

  return arrays

def load_molecules(paths, feature_types=["fingerprints"]):

          labels[smiles][target_name] = measurement 

  return labels, splits

def load_datasets(paths, prediction_endpoint, split_endpoint, datatype="vs",

    **load_args):

  """Dispatches to correct loader depending on type of data."""

  if datatype == "vs":

    return load_vs_datasets(paths, prediction_endpoint,

                            split_endpoint, **load_args)

  elif datatype == "pdbbind":

    return load_pdbbind_datasets(paths, prediction_endpoint, **load_args)

  else:

    raise ValueError("Unsupported datatype.")

def load_vs_datasets(paths, prediction_endpoint, split_endpoint, target_dir_name="targets",

def load_datasets(paths, prediction_endpoint, split_endpoint, target_dir_name="targets",

    feature_types=["fingerprints"]):

  """Load both labels and fingerprints.

def load_and_transform_dataset(paths, input_transforms, output_transforms,

    prediction_endpoint, split_endpoint=None, labels_endpoint="labels", weight_positives=True,

    datatype="vs", feature_types=["fingerprints"]):

    datatype="tensor", feature_types=["fingerprints"]):

  """Transform data labels as specified

  Parameters

    are performed in the order specified. An empty list corresponds to no

    transformations. Only for regression outputs.

  """

  dataset = load_datasets(paths, prediction_endpoint, split_endpoint, datatype=datatype,

  dataset = load_datasets(paths, prediction_endpoint, split_endpoint,

      feature_types=feature_types)

  if datatype == "vs":

  if datatype == "vector":

    X, y, W = dataset_to_numpy(dataset, weight_positives=weight_positives)

  elif datatype == "pdbbind":

  elif datatype == "tensor":

    X, y, W = tensor_dataset_to_numpy(dataset)

  y = transform_outputs(y, W, output_transforms,

      weight_positives=weight_positives)