Bulk movement of code from private repo.

"""

Code for processing the Google vs-datasets using keras.

"""

import numpy as np

from keras.models import Graph

from keras.models import Sequential

from keras.layers.core import Dense, Dropout, Activation

from keras.optimizers import SGD

from dataset_arxiv import load_datasets

from dataset_arxiv import multitask_to_singletask

from dataset_arxiv import train_test_random_split

from dataset_arxiv import train_test_scaffold_split

from dataset_arxiv import dataset_to_numpy

from dataset_arxiv import to_one_hot

from dataset_arxiv import eval_model

from dataset_arxiv import compute_r2_scores

from dataset_arxiv import compute_rms_scores

from dataset_arxiv import compute_roc_auc_scores

from dataset_arxiv import load_and_transform_dataset

def process_multitask(paths, task_transforms, desc_transforms, splittype="random",

    seed=None, add_descriptors=False, desc_weight=0.5):

  """Extracts multitask datasets and splits into train/test.

  Returns a tuple of test/train datasets, fingerprints, and labels.

  TODO(rbharath): This function is ugly. Returns way too many arguments. Clean

  it up.

  Parameters

  ----------

  paths: list 

    List of paths to Google vs datasets. 

  task_transforms: dict 

    dict mapping target names to label transform. Each output type must be either

    None or "log". Only for regression outputs.

  desc_transforms: dict

    dict mapping descriptor number to transform. Each transform must be

    either None, "log", "normalize", or "log-normalize"

  splittype: string

    Must be "random" or "scaffold"

  seed: int

    Seed used for random splits.

  """

  dataset = load_and_transform_dataset(paths, task_transforms, desc_transforms,

      add_descriptors=add_descriptors)

  if splittype == "random":

    train, test = train_test_random_split(dataset, seed=seed)

  elif splittype == "scaffold":

    train, test = train_test_scaffold_split(dataset)

  else:

    raise ValueError("Improper splittype. Must be random/scaffold.")

  X_train, y_train, W_train = dataset_to_numpy(train,

      add_descriptors=add_descriptors, desc_weight=desc_weight)

  X_test, y_test, W_test = dataset_to_numpy(test,

      add_descriptors=add_descriptors, desc_weight=desc_weight)

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

def process_singletask(paths, task_transforms, desc_transforms, splittype="random", seed=None,

    add_descriptors=False, desc_weight=0.5):

  """Extracts singletask datasets and splits into train/test.

  Returns a dict that maps target names to tuples.

  Parameters

  ----------

  paths: list 

    List of paths to Google vs datasets. 

  task_transforms: dict 

    dict mapping target names to label transform. Each output type must be either

    None or "log". Only for regression outputs.

  splittype: string

    Must be "random" or "scaffold"

  seed: int

    Seed used for random splits.

  """

  dataset = load_and_transform_dataset(paths, task_transforms, desc_transforms,

      add_descriptors=add_descriptors)

  singletask = multitask_to_singletask(dataset)

  arrays = {}

  for target in singletask:

    data = singletask[target]

    # TODO(rbharath): Remove limitation after debugging.

    if len(data) == 0:

      continue

    if splittype == "random":

      train, test = train_test_random_split(data, seed=seed)

    elif splittype == "scaffold":

      train, test = train_test_scaffold_split(data)

    else:

      raise ValueError("Improper splittype. Must be random/scaffold.")

    X_train, y_train, W_train = dataset_to_numpy(train,

        add_descriptors=add_descriptors, desc_weight=desc_weight)

    X_test, y_test, W_test = dataset_to_numpy(test,

        add_descriptors=add_descriptors, desc_weight=desc_weight)

    arrays[target] = (train, X_train, y_train, W_train, test, X_test, y_test,

        W_test)

  return arrays

def fit_multitask_mlp(paths, task_types, task_transforms, desc_transforms,

                      splittype="random", add_descriptors=False, desc_weight=0.5,

                      **training_params):

  """

  Perform stochastic gradient descent optimization for a keras multitask MLP.

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

  Parameters

  ----------

  paths: list 

    List of paths to Google vs datasets. 

  task_types: dict 

    dict mapping target names to output type. Each output type must be either

    "classification" or "regression".

  task_transforms: dict 

    dict mapping target names to label transform. Each output type must be either

    None or "log". Only for regression outputs.

  desc_transforms: dict

    dict mapping descriptor number to transform. Each transform must be

    either None, "log", "normalize", or "log-normalize"

  add_descriptors: bool

    Add descriptor prediction as extra task.

  training_params: dict

    Aggregates keyword parameters to pass to train_multitask_model

  """

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

      process_multitask(paths, task_transforms, desc_transforms,

      splittype=splittype, add_descriptors=add_descriptors, desc_weight=desc_weight))

  print np.shape(y_train)

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

                                desc_transforms, add_descriptors=add_descriptors,

                                **training_params)

  results = eval_model(test, model, task_types, desc_transforms,

      add_descriptors=add_descriptors, modeltype="keras_multitask")

  if add_descriptors:

    local_task_types = task_types.copy()

    for desc in desc_transforms:

      local_task_types[desc] = "regression"

  else:

    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()))

  #rms = compute_rms_scores(results, local_task_types)

  #if rms:

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

  #return (aucs, r2s, rms)

def fit_singletask_mlp(paths, task_types, task_transforms,

                       desc_transforms, splittype="random",

                       add_descriptors=False, desc_weight=0.5,

                       **training_params):

  """

  Perform stochastic gradient descent optimization for a keras MLP.

  paths: list 

    List of paths to Google vs datasets. 

  task_types: dict 

    dict mapping target names to output type. Each output type must be either

    "classification" or "regression".

  task_transforms: dict 

    dict mapping target names to label transform. Each output type must be either

    None or "log". Only for regression outputs.

  desc_transforms: dict

    dict mapping descriptor number to transform. Each transform must be

    either None, "log", "normalize", or "log-normalize"

  training_params: dict

    Aggregates keyword parameters to pass to train_multitask_model

  """

  singletasks = process_singletask(paths, task_transforms, desc_transforms,

    splittype=splittype, add_descriptors=add_descriptors,

    desc_weight=desc_weight)

  ret_vals = {}

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

  for index, target in enumerate(singletasks):

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

        singletasks[target])

    model = train_multitask_model(X_train, y_train, W_train,

        {target: task_types[target]}, desc_transforms, add_descriptors=add_descriptors,

        **training_params)

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

                         desc_transforms,

                         # We run singletask models as special cases of

                         # multitask.

                         modeltype="keras_multitask",

                         add_descriptors=add_descriptors)

    print "Target %s" % 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 "Mean AUC: %f" % np.mean(np.array(aucs.values()))

  if r2s:

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

  if rms:

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

def train_multitask_model(X, y, W, task_types, desc_transforms, add_descriptors=False,

                      learning_rate=0.01, decay=1e-6,

                      momentum=0.9, nesterov=True, activation="relu",

                      dropout=0.5, nb_epoch=20, batch_size=50, n_hidden=500,

                      n_input=1024, validation_split=0.1):

  """

  Perform stochastic gradient descent optimization for a keras multitask MLP.

  Returns a trained model.

  TODO(rbharath): The handling of add_descriptors for semi-supervised learning

  is horrible. Refactor.

  Parameters

  ----------

  X: np.ndarray

    Feature matrix

  y: np.ndarray

    Label matrix

  W: np.ndarray

    Weight matrix

  task_types: dict 

    dict mapping target names to output type. Each output type must be either

    "classification" or "regression".

  desc_transforms: dict

    dict mapping descriptor number to transform. Each transform must be

    either None, "log", "normalize", or "log-normalize"

  add_descriptors: bool

    Add descriptor prediction as extra task.

  learning_rate: float

    Learning rate used.

  decay: float

    Learning rate decay.

  momentum: float

    Momentum used in SGD.

  nesterov: bool

    Use Nesterov acceleration

  n_epochs: int

    maximal number of epochs to run the optimizer

  """

  eps = .001

  num_tasks = len(task_types)

  sorted_targets = sorted(task_types.keys())

  if add_descriptors:

    sorted_descriptors = sorted(desc_transforms.keys())

    endpoints = sorted_targets + sorted_descriptors

    local_task_types = task_types.copy()

    for desc in desc_transforms:

      local_task_types[desc] = "regression"

  else:

    local_task_types = task_types.copy()

    endpoints = sorted_targets

  print "endpoints: " + str(endpoints)

  # Add eps weight to avoid minibatches with zero weight (causes theano to crash).

  W = W + eps * np.ones(np.shape(W))

  model = Graph()

  model.add_input(name="input", ndim=n_input)

  model.add_node(

      Dense(n_input, n_hidden, init='uniform', activation=activation),

      name="dense", input="input")

  model.add_node(Dropout(dropout), name="dropout", input="dense")

  top_layer = "dropout"

  for task, target in enumerate(endpoints):

    task_type = local_task_types[target]

    if task_type == "classification":

      model.add_node(

          Dense(n_hidden, 2, init='uniform', activation="softmax"),

          name="dense_head%d" % task, input=top_layer)

    elif task_type == "regression":

      model.add_node(

          Dense(n_hidden, 1, init='uniform'),

          name="dense_head%d" % task, input=top_layer)

    model.add_output(name="task%d" % task, input="dense_head%d" % task)

  data_dict, loss_dict, sample_weights = {}, {}, {}

  data_dict["input"] = X

  for task, target in enumerate(endpoints):

    task_type = local_task_types[target]

    taskname = "task%d" % task

    sample_weights[taskname] = W[:, task]

    if task_type == "classification":

      loss_dict[taskname] = "binary_crossentropy"

      data_dict[taskname] = to_one_hot(y[:,task])

    elif task_type == "regression":

      loss_dict[taskname] = "mean_squared_error"

      data_dict[taskname] = y[:,task]

  sgd = SGD(lr=learning_rate, decay=decay, momentum=momentum, nesterov=nesterov)

  print "About to compile model!"

  model.compile(optimizer=sgd, loss=loss_dict)

  print "Done compiling. About to fit model!"

  print "validation_split: " + str(validation_split)

  print "decay: " + str(decay)

  model.fit(data_dict, nb_epoch=nb_epoch, batch_size=batch_size, validation_split=validation_split,

            sample_weight=sample_weights)

  return model

"""

Code for processing the Google vs-datasets using scikit-learn.

"""

import numpy as np

from dataset_arxiv import load_and_transform_dataset

from dataset_arxiv import multitask_to_singletask

from dataset_arxiv import train_test_random_split

from dataset_arxiv import train_test_scaffold_split

from dataset_arxiv import dataset_to_numpy

from dataset_arxiv import eval_model

from dataset_arxiv import compute_r2_scores

from dataset_arxiv import compute_rms_scores

from dataset_arxiv import compute_roc_auc_scores

from sklearn.ensemble import RandomForestClassifier

from sklearn.ensemble import RandomForestRegressor

from sklearn.linear_model import MultiTaskLasso 

from sklearn.linear_model import LogisticRegression

from sklearn.linear_model import LinearRegression

from sklearn.linear_model import RidgeCV

from sklearn.linear_model import LassoCV

from sklearn.linear_model import ElasticNetCV

from sklearn.linear_model import LassoLarsCV

from sklearn.svm import SVR

def fit_singletask_models(paths, modeltype, task_types, task_transforms,

    add_descriptors=False, desc_transforms={}, splittype="random",

    seed=None):

  """Fits singletask linear regression models to potency.

  Parameters

  ----------

  paths: list 

    List of paths to Google vs datasets. 

  modeltype: String

    A string describing the model to be trained. Options are RandomForest,

  splittype: string

    Type of split for train/test. Either random or scaffold.

  seed: int (optional)

    Seed to initialize np.random.

  task_types: dict 

    dict mapping target names to output type. Each output type must be either

    "classification" or "regression".

  task_transforms: dict 

    dict mapping target names to label transform. Each output type must be either

    None or "log". Only for regression outputs.

  desc_transforms: dict

    dict mapping descriptor number to transform. Each transform must be

    either None, "log", "normalize", or "log-normalize"

  """

  dataset = load_and_transform_dataset(paths, task_transforms, desc_transforms,

      add_descriptors=add_descriptors)

  singletask = multitask_to_singletask(dataset)

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

  for target in singletask:

    data = singletask[target]

    if splittype == "random":

      train, test = train_test_random_split(data, seed=seed)

    elif splittype == "scaffold":

      train, test = train_test_scaffold_split(data)

    else:

      raise ValueError("Improper splittype. Must be random/scaffold.")

    X_train, y_train, W_train = dataset_to_numpy(train)

    X_test, y_test, W_test = dataset_to_numpy(test)

    if modeltype == "RandomForestRegressor":

      model = RandomForestRegressor(n_estimators=500, n_jobs=-1,

          warm_start=True, max_features="sqrt")

    elif modeltype == "RandomForestClassifier":

      model = RandomForestClassifier(n_estimators=500, n_jobs=-1,

          warm_start=True, max_features="sqrt")

    elif modeltype == "LogisticRegression":

      model = LogisticRegression(class_weight="auto")

    elif modeltype == "LinearRegression":

      model = LinearRegression(normalize=True)

    elif modeltype == "RidgeRegression":

      model = RidgeCV(alphas=[0.01, 0.1, 1.0, 10.0], normalize=True) 

    elif modeltype == "Lasso":

      model = LassoCV(max_iter=2000, n_jobs=-1) 

    elif modeltype == "LassoLars":

      model = LassoLarsCV(max_iter=2000, n_jobs=-1) 

    elif modeltype == "ElasticNet":

      model = ElasticNetCV(max_iter=2000, n_jobs=-1) 

    elif modeltype == "SVR-rbf":

      model = SVR(kernel="rbf") 

    elif modeltype == "SVR-poly":

      model = SVR(kernel="poly") 

    elif modeltype == "SVR-linear":

      model = SVR(kernel="linear") 

    else:

      raise ValueError("Invalid model type provided.")

    model.fit(X_train, y_train.ravel())

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

        desc_transforms, modeltype="sklearn", add_descriptors=add_descriptors)

    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 "Mean AUC: %f" % np.mean(np.array(aucs.values()))

  if r2s:

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

  if rms:

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

def fit_multitask_rf(dataset, splittype="random"):

  """Fits a multitask RF model to provided dataset.

  Performs a random 80-20 train/test split.

  Parameters

  ----------

  dataset: dict 

    A dictionary of type produced by load_datasets. 

  splittype: string

    Type of split for train/test. Either random or scaffold.

  """

  if splittype == "random":

    train, test = train_test_random_split(data, seed=0)

  elif splittype == "scaffold":

    train, test = train_test_scaffold_split(data)

  else:

    raise ValueError("Improper splittype. Must be random/scaffold.")

  X_train, y_train, W_train = dataset_to_numpy(train)

  classifier = RandomForestClassifier(n_estimators=100, n_jobs=-1,

      class_weight="auto")

  classifier.fit(X_train, y_train)

  results = eval_model(test, classifier)

  scores = compute_roc_auc_scores(results)

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

#!/bin/bash

# Download PCBA JSON files for all AIDs into the current directory

# usage: download_pcba_json.sh [n_jobs=8]

if [ -z $1 ]

then

  n_jobs=8

else

  n_jobs=$1

fi

url=ftp://ftp.ncbi.nlm.nih.gov/pubchem/Bioassay/JSON/

echo "Downloading PCBA JSON files (${n_jobs} jobs)..."

curl -sl ${url} | awk "{print \"${url}\"\$1}" | xargs -n 1 -P ${n_jobs} wget -q

"""

Extract JSON assay descriptions (leave assay data behind).

"""

import argparse

import glob

import gzip

import json

import os

from vs_utils.utils.public_data import PcbaJsonParser

def parse_args(input_args=None):

  """

  Parse command-line arguments.

  Parameters

  ----------

  input_args : list, optional

    Input arguments. If not provided, defaults to sys.argv[1:].

  """

  parser = argparse.ArgumentParser()

  parser.add_argument('files', nargs='+',

                      help='Directories containing PCBA JSON files.')

  return parser.parse_args(input_args)

def main(dirs):

  for this_dir in dirs:

    print this_dir

    for filename in glob.glob(os.path.join(this_dir, '*.json.gz')):

      parser = PcbaJsonParser(filename)

      tree = parser.tree

      aid = parser.get_aid()

      try:

        del tree['PC_AssaySubmit']['data']

      except KeyError as e:

        print 'JSON is not properly formatted. Please follow NCBI FTP format.'

        raise e

      with gzip.open(os.path.join(

              this_dir, '{}-desc.json.gz'.format(aid)), 'wb') as f:

          json.dump(tree, f, indent=2)

if __name__ == '__main__':

    args = parse_args()

    main(args.files)