Dataset loading merged, best results included

from deepchem.models.tensorflow_models.fcnet import TensorflowMultiTaskRegressor

from deepchem.models.tensorflow_models.fcnet import TensorflowMultiTaskClassifier

from deepchem.models.tensorflow_models.robust_multitask import RobustMultitaskRegressor

from deepchem.models.tensorflow_models.lr import TensorflowLogisticRegression

 No newline at end of file

# -*- coding: utf-8 -*-

"""

Created on Tue Nov 08 14:10:02 2016

@author: Zhenqin Wu

"""

import tensorflow as tf

import numpy as np

import os

import time

from deepchem.metrics import from_one_hot

from deepchem.models.tensorflow_models import TensorflowGraph

from deepchem.models.tensorflow_models import TensorflowGraphModel

from deepchem.models.tensorflow_models import model_ops

from deepchem.utils.save import log

from deepchem.datasets import pad_features

from deepchem.metrics import to_one_hot

from deepchem.data import pad_features

def weight_decay(penalty_type, penalty):

  variables = []

  # exclude bias variables

  for v in tf.trainable_variables():

    if v.get_shape().as_list()[0] > 1:

      variables.append(v)

  with tf.name_scope('weight_decay'):

    if penalty_type == 'l1':

      cost = tf.add_n([tf.reduce_sum(tf.Abs(v)) for v in variables])

    elif penalty_type == 'l2':

      cost = tf.add_n([tf.nn.l2_loss(v) for v in variables])

    else:

      raise NotImplementedError('Unsupported penalty_type %s' % penalty_type)

    cost *= penalty

    tf.scalar_summary('Weight Decay Cost', cost)

  return cost    

class TensorflowLogisticRegression(TensorflowGraphModel):

  """ A simple tensorflow based logistic regression model. """

  def build(self, graph, name_scopes, training):

    """Constructs the graph architecture of model: n_tasks * sigmoid nodes.

    This method creates the following Placeholders:

      mol_features: Molecule descriptor (e.g. fingerprint) tensor with shape

        batch_size x n_features.

    """

    placeholder_scope = TensorflowGraph.get_placeholder_scope(

        graph, name_scopes)

    n_features = self.n_features

    with graph.as_default():

      with placeholder_scope:

        self.mol_features = tf.placeholder(

            tf.float32,

            shape=[None, n_features],

            name='mol_features')

      weight_init_stddevs = self.weight_init_stddevs

      bias_init_consts = self.bias_init_consts

      lg_list = []

      for task in range(self.n_tasks):

        lg = model_ops.fully_connected_layer(

            tensor=self.mol_features,

            size = 1,

            weight_init=tf.truncated_normal(

                shape=[self.n_features, 1],

                stddev=weight_init_stddevs[0]),

            bias_init=tf.constant(value=bias_init_consts[0],

                                  shape=[1]))

        lg_list.append(lg)

    return lg_list

  def add_label_placeholders(self, graph, name_scopes):

    labels = []

    placeholder_scope = TensorflowGraph.get_placeholder_scope(graph, name_scopes)

    with placeholder_scope:

      for task in range(self.n_tasks):

        labels.append(tf.identity(

            tf.placeholder(tf.float32, shape=[None,1],

                           name='labels_%d' % task)))

    return labels

  def add_training_cost(self, graph, name_scopes, output, labels, weights):

    with graph.as_default():

      epsilon = 1e-3  # small float to avoid dividing by zero

      weighted_costs = []  # weighted costs for each example

      gradient_costs = []  # costs used for gradient calculation

      with TensorflowGraph.shared_name_scope('costs', graph, name_scopes):

        for task in range(self.n_tasks):

          task_str = str(task).zfill(len(str(self.n_tasks)))

          with TensorflowGraph.shared_name_scope(

              'cost_{}'.format(task_str), graph, name_scopes):

            with tf.name_scope('weighted'):

              weighted_cost = self.cost(output[task], labels[task],

                                        weights[task])

              weighted_costs.append(weighted_cost)

            with tf.name_scope('gradient'):

              # Note that we divide by the batch size and not the number of

              # non-zero weight examples in the batch.  Also, instead of using

              # tf.reduce_mean (which can put ops on the CPU) we explicitly

              # calculate with div/sum so it stays on the GPU.

              gradient_cost = tf.div(tf.reduce_sum(weighted_cost),

                                     self.batch_size)

              gradient_costs.append(gradient_cost)

        # aggregated costs

        with TensorflowGraph.shared_name_scope('aggregated', graph, name_scopes):

          with tf.name_scope('gradient'):

            loss = tf.add_n(gradient_costs)

          # weight decay

          if self.penalty != 0.0:

            penalty = weight_decay(self.penalty_type, self.penalty)

            loss += penalty

      return loss 

  def cost(self, logits, labels, weights):

    return tf.mul(tf.nn.sigmoid_cross_entropy_with_logits(logits, labels),

                  weights)

  def add_output_ops(self, graph, output):

    with graph.as_default():

      sigmoid = []

      with tf.name_scope('inference'):

        for i, logits in enumerate(output):

          sigmoid.append(tf.nn.sigmoid(logits, name='sigmoid_%d' % i))

      output = sigmoid

    return output 

  def construct_feed_dict(self, X_b, y_b=None, w_b=None, ids_b=None):

    orig_dict = {}

    orig_dict["mol_features"] = X_b

    for task in range(self.n_tasks):

      if y_b is not None:

        y_2column = to_one_hot(y_b[:, task])

        orig_dict["labels_%d" % task] = y_2column[:,1:2]

      else:

        # Dummy placeholders

        orig_dict["labels_%d" % task] = np.zeros((self.batch_size,1))

      if w_b is not None:

        orig_dict["weights_%d" % task] = w_b[:, task]

      else:

        # Dummy placeholders

        orig_dict["weights_%d" % task] = np.ones(

            (self.batch_size,)) 

    return TensorflowGraph.get_feed_dict(orig_dict)

  def predict_proba_on_batch(self, X, pad_batch=False):

    if pad_batch:

      X = pad_features(self.batch_size, X)

    if not self._restored_model:

      self.restore()

    with self.eval_graph.graph.as_default():

      # run eval data through the model

      n_tasks = self.n_tasks

      with self._get_shared_session(train=False).as_default():

        feed_dict = self.construct_feed_dict(X)

        data = self._get_shared_session(train=False).run(

            self.eval_graph.output, feed_dict=feed_dict)

        batch_outputs = np.asarray(data[:n_tasks], dtype=float)

        # reshape to batch_size x n_tasks x ...

        complimentary = np.ones(np.shape(batch_outputs))

        complimentary = complimentary - batch_outputs

        batch_outputs = np.squeeze(np.stack(arrays = [complimentary,

						      batch_outputs],

                                            axis = 2))

        if batch_outputs.ndim == 3:

          batch_outputs = batch_outputs.transpose((1, 0, 2))

        elif batch_outputs.ndim == 2:

          batch_outputs = batch_outputs.transpose((1, 0))

        else:

          raise ValueError(

              'Unrecognized rank combination for output: %s ' %

              (batch_outputs.shape,))

      outputs = batch_outputs

    return np.copy(outputs)

  def predict_on_batch(self, X, pad_batch=False):

    if pad_batch:

      X = pad_features(self.batch_size, X)

    if not self._restored_model:

      self.restore()

    with self.eval_graph.graph.as_default():

      # run eval data through the model

      n_tasks = self.n_tasks

      output = []

      start = time.time()

      with self._get_shared_session(train=False).as_default():

        feed_dict = self.construct_feed_dict(X)

        data = self._get_shared_session(train=False).run(

            self.eval_graph.output, feed_dict=feed_dict)

        batch_output = np.asarray(data[:n_tasks], dtype=float)

        # reshape to batch_size x n_tasks x ...

        complimentary = np.ones(np.shape(batch_output))

        complimentary = complimentary - batch_output

        batch_output = np.squeeze(np.stack(arrays = [complimentary,

                                                     batch_output],

                                            axis = 2))

        if batch_output.ndim == 3:

          batch_output = batch_output.transpose((1, 0, 2))

        elif batch_output.ndim == 2:

          batch_output = batch_output.transpose((1, 0))

        else:

          raise ValueError(

              'Unrecognized rank combination for output: %s' %

              (batch_output.shape,))

        output.append(batch_output)

        outputs = np.array(from_one_hot(

            np.squeeze(np.concatenate(output)), axis=-1))

    return np.copy(outputs)

from sider.sider_datasets import load_sider

def benchmark_loading_datasets(base_dir_o, hyper_parameters, n_features = 1024, 

                               dataset_name='all',model='all',reload = True,

                               dataset_name='all',model='tf',reload = True,

                               verbosity='high',out_path='/tmp'):

  """

  Loading dataset for benchmark test

  dataset_name : string, optional (default='all')

      choice of which dataset to use, 'all' = computing all the datasets

  model : string,  optional (default='all')

      choice of which model to use, 'all' = running all models on the dataset

  model : string,  optional (default='tf')

      choice of which model to use, should be: tf, logreg, graphconv

  out_path : string, optional(default='/tmp')

      path of result file

  """

  assert dataset_name in ['all', 'muv', 'nci', 'pcba', 'tox21','sider',

                          'toxcast']

  if not dataset_name in ['all','muv','nci','pcba','tox21','sider','toxcast']:

    raise ValueError('Dataset not supported')

  if model in ['graphconv']:

    method = 'GraphConv'

  elif model in ['tf','logreg','rf']:

    method = 'ECFP'

  else:

    raise ValueError('Model not supported')

  if dataset_name == 'all':

    #currently not including the nci dataset

    time_start = time.time()

    #loading datasets     

    tasks,datasets,transformers = loading_functions[dname]()

    tasks,datasets,transformers = loading_functions[dname](method = method)

    train_dataset, valid_dataset, test_dataset = datasets

    time_finish_loading = time.time()

    #time_finish_loading-time_start is the time(s) used for dataset loading

    #running model

    train_score,valid_score = benchmark_train_and_valid(base_dir,train_dataset,

                                                        valid_dataset, tasks,

                                                        transformers,

                                                        hyper_parameters,

                                                        n_features=n_features,

                                                        model = model,

                                                        verbosity = verbosity)

    time_finish_running = time.time()

    #time_finish_running-time_finish_loading is the time(s) used for fitting and evaluating

    for i, hp in enumerate(hyper_parameters[model]):

      time_start_fitting = time.time()

      train_score,valid_score = benchmark_train_and_valid(base_dir,

                                    train_dataset,valid_dataset,tasks,

                                    transformers,hp,n_features=n_features,

                                    model = model,verbosity = verbosity)      

      time_finish_fitting = time.time()

      with open(os.path.join(out_path,'results.csv'),'a') as f:

        f.write('\n\n'+str(i))

        f.write('\n'+dname+',train')

        for i in train_score:

          f.write(','+i+','+str(train_score[i]['mean-roc_auc_score']))

        f.write('\n'+dname+',valid')

        for i in valid_score:

          f.write(','+i+','+str(valid_score[i]['mean-roc_auc_score'])) 

      #output timing data: running time include all the model

      f.write('\n'+dname+',time_for_loading,,'+

              str(time_finish_loading-time_start)+'seconds')

      f.write('\n'+dname+',time_for_running,,'+

              str(time_finish_running-time_finish_loading)+'seconds')

        f.write('\n'+dname+',time_for_running,'+

              str(time_finish_fitting-time_start_fitting))

    #clear workspace         

    del tasks,datasets,transformers

    del train_dataset,valid_dataset, test_dataset

    del time_start,time_finish_loading,time_finish_running

  return None

                                            verbosity=verbosity,

                                            mode="classification")

  assert model in ['all', 'tf', 'rf']

  assert model in ['graphconv', 'tf', 'rf','logreg']

  if model == 'all' or model == 'tf':

    # Initialize model folder

    model_dir_tf = os.path.join(base_dir, "model_tf")

    dropouts = hyper_parameters['tf'][0]

    learning_rate = hyper_parameters['tf'][1]

    weight_init_stddevs = hyper_parameters['tf'][2]

    batch_size = hyper_parameters['tf'][3]

      # Building tensorflow MultiTaskDNN model

    tensorflow_model = dc.models.TensorflowMultiTaskClassifier(

        len(tasks), n_features, dropouts=[dropouts],

        learning_rate=learning_rate, weight_init_stddevs=[weight_init_stddevs],

        batch_size=batch_size, verbosity=verbosity)

    dropouts = hyper_parameters['dropouts']

    learning_rate = hyper_parameters['learning_rate']

    layer_sizes = hyper_parameters['layer_sizes']

    penalty = hyper_parameters['penalty']

    batch_size = hyper_parameters['batch_size']

    nb_epoch = hyper_parameters['nb_epoch']

    tensorflow_model = dc.models.TensorflowMultiTaskClassifier(len(tasks),

          n_features,  learning_rate=learning_rate, layer_sizes=layer_sizes, 

          dropouts=dropouts, penalty=penalty, batch_size=batch_size, 

          verbosity=verbosity)

    model_tf = dc.models.TensorflowModel(tensorflow_model)

    print('-------------------------------------')

    print('Start fitting by tensorflow')

    model_tf.fit(train_dataset)

    model_tf.fit(train_dataset,nb_epoch = nb_epoch)

    train_scores['tensorflow'] = model_tf.evaluate(train_dataset,

                                        [classification_metric],transformers)

    shutil.rmtree(base_dir_o)

  os.makedirs(base_dir_o)

  #Datasets and models used in the benchmark test, all=all the datasets(models)

  dataset_name = sys.argv[1]

  model = sys.argv[2]

  #Datasets and models used in the benchmark test, all=all the datasets

  dataset_name = 'tox21'

  model = 'tf'

  #input hyperparameters

  #tf: dropouts, learning rate, weight initial stddev, batch_size

  hyper_parameters = {'tf':[0.25, 0.0003, 0.1, 50]}

  #tf: dropouts, learning rate, layer_sizes, weight initial stddev,penalty,

  #    batch_size

  hps = {}

  hps['tf'] = [{'dropouts':[0.25],'learning_rate':0.001,'layer_sizes':[1000],

                'penalty':0.0, 'batch_size':50, 'nb_epoch':10}]

  benchmark_loading_datasets(base_dir_o,hyper_parameters,n_features = 1024,

  benchmark_loading_datasets(base_dir_o,hps,n_features = 1024,

                             dataset_name = dataset_name, model = model,

                             reload = reload, verbosity = verbosity)

                             reload = reload, verbosity = 'high')

sider,tf,train,0.931,valid,0.647,learning_rate,0.0003~0.003,dropouts,0.2~0.3,layer_sizes,1000~1200

tox21,tf,train,0.996,valid,0.763,learning_rate,0.001~0.003,dropouts,0.4~0.5,layer_sizes,1000~1200

toxcast,tf,train,0.944,valid,0.699,learning_rate,0.0004,dropouts,0.4,layer_sizes,500

muv,tf,train,0.980,valid,0.710,learning_rate,0.0008~0.0013,dropouts,0.25~0.4,layer_sizes,1200

pcba,graphconv,train,0.866,valid,0.836,learning_rate,0.0001,learning_rate_decay,1500,layer_structure,2layers&64filters

muv,graphconv,train,0.881,valid,0.832,learning_rate,0.0004~0.0008,learning_rate_decay,1500,layer_structure,2layers&64~128filters

tox21,graphconv,train,0.930,valid,0.819,learning_rate,0.0003~0.002,learning_rate_decay,1000~2000,layer_structure,2layers&128filters

import shutil

import deepchem as dc

def load_muv():

def load_muv(method = 'ECFP'):

  """Load MUV datasets. Does not do train/test split"""

  # Load MUV dataset

  print("About to load MUV dataset.")

      current_dir, "../../datasets/muv.csv.gz")

  # Featurize MUV dataset

  print("About to featurize MUV dataset.")

  if method == 'ECFP':

      featurizer = dc.feat.CircularFingerprint(size=1024)

  elif method == 'GraphConv':

      featurizer = dc.feat.ConvMolFeaturizer()

  MUV_tasks = sorted(['MUV-692', 'MUV-689', 'MUV-846', 'MUV-859', 'MUV-644',

                      'MUV-548', 'MUV-852', 'MUV-600', 'MUV-810', 'MUV-712',

                      'MUV-737', 'MUV-858', 'MUV-713', 'MUV-733', 'MUV-652',

    dataset = transformer.transform(dataset)

  splitter = dc.splits.IndexSplitter()

  train, valid, test = splitter.train_valid_test_split(dataset)

  train, valid, test = splitter.train_valid_test_split(

	dataset, compute_feature_statistics=False)

  return MUV_tasks, (train, valid, test), transformers