text:
from builtins import range
from builtins import object
import numpy as np

from past.builtins import xrange

from layers import *
from layer_utils import *

class TwoLayerNet(object):
    """
    A two-layer fully-connected neural network with ReLU nonlinearity and
    softmax loss that uses a modular layer design. We assume an input dimension
    of D, a hidden dimension of H, and perform classification over C classes.

    The architecure should be affine - relu - affine - softmax.

    Note that this class does not implement gradient descent; instead, it
    will interact with a separate Solver object that is responsible for running
    optimization.

    The learnable parameters of the model are stored in the dictionary
    self.params that maps parameter names to numpy arrays.
    """

    def __init__(self, input_dim=3 * 32 * 32, hidden_dim=100, num_classes=10,
                 weight_scale=1e-3, reg=0.0):
        """
        Initialize a new network.

        Inputs:
        - input_dim: An integer giving the size of the input
        - hidden_dim: An integer giving the size of the hidden layer
        - num_classes: An integer giving the number of classes to classify
        - weight_scale: Scalar giving the standard deviation for random
          initialization of the weights.
        - reg: Scalar giving L2 regularization strength.
        """
        self.params = {}
        self.reg = reg

        ############################################################################
        # TODO: Initialize the weights and biases of the two-layer net. Weights    #
        # should be initialized from a Gaussian with standard deviation equal to   #
        # weight_scale, and biases should be initialized to zero. All weights and  #
        # biases should be stored in the dictionary self.params, with first layer  #
        # weights and biases using the keys 'W1' and 'b1' and second layer weights #
        # and biases using the keys 'W2' and 'b2'.                                 #
        # See also: http://cs231n.github.io/neural-networks-2/#init                #
        ############################################################################
        # Initialize the weights in the self.params dictionary, with random numbers
        # multiplied by the weight scales. Dimensions should follow network layer
        # dimensions, i.e. first layer dim [input x hidden], second layer [hiddenxoutput]
        self.params["W1"] = weight_scale * np.random.randn(input_dim, hidden_dim)
        self.params["W2"] = weight_scale * np.random.randn(hidden_dim, num_classes)
        # Initialize the biases with 0s corresponding to the right vector dimensions
        self.params["b1"] = np.zeros(hidden_dim)
        self.params["b2"] = np.zeros(num_classes)
        ############################################################################
        #                             END OF YOUR CODE                             #
        ############################################################################

    def loss(self, X, y=None):
        """
        Compute loss and gradient for a minibatch of data.

        Inputs:
        - X: Array of input data of shape (N, d_1, ..., d_k)
        - y: Array of labels, of shape (N,). y[i] gives the label for X[i].

        Returns:
        If y is None, then run a test-time forward pass of the model and return:
        - scores: Array of shape (N, C) giving classification scores, where
          scores[i, c] is the classification score for X[i] and class c.

        If y is not None, then run a training-time forward and backward pass and
        return a tuple of:
        - loss: Scalar value giving the loss
        - grads: Dictionary with the same keys as self.params, mapping parameter
          names to gradients of the loss with respect to those parameters.
        """
        scores = None
        ############################################################################
        # TODO: Implement the forward pass for the two-layer net, computing the    #
        # class scores for X and storing them in the scores variable.              #
        ############################################################################
        # Initialize Weights and Biases for use in the functions below
        W1 = self.params["W1"]
        W2 = self.params["W2"]
        b1 = self.params["b1"]
        b2 = self.params["b2"]

        # Perform forward pass with Relu Activation for the hidden layer
        h, cachel1 = affine_relu_forward(X, W1, b1)
        # Use output of hidden layer to complete the forward pass
        scores, cachel2 = affine_forward(h, W2, b2)
        ############################################################################
        #                             END OF YOUR CODE                             #
        ############################################################################

        # If y is None then we are in test mode so just return scores
        if y is None:
            return scores

        loss, grads = 0, {}
        ############################################################################
        # TODO: Implement the backward pass for the two-layer net. Store the loss  #
        # in the loss variable and gradients in the grads dictionary. Compute data #
        # loss using softmax, and make sure that grads[k] holds the gradients for  #
        # self.params[k]. Don't forget to add L2 regularization!                   #
        #                                                                          #
        # NOTE: To ensure that your implementation matches ours and you pass the   #
        # automated tests, make sure that your L2 regularization includes a factor #
        # of 0.5 to simplify the expression for the gradient.                      #
        ############################################################################
        # Evaluate the Softmax loss and derivative of scores
        loss, dscores = softmax_loss(scores, y)
        # Apply Regularization to the loss
        loss += 0.5 * self.reg * (np.sum(W1 ** 2) + np.sum(W2 ** 2))
        # Use the Affine_backward module to compute backward pass of affine layer
        dh, grads["W2"], grads["b2"] = affine_backward(dscores, cachel2)
        # Use the Affine_relu_backward to compute the backward pass for a layer of
        # rectified linear units (ReLUs).
        dx, grads["W1"], grads["b1"] = affine_relu_backward(dh, cachel1)
        # Apply Regularization to the Weight gradients
        grads["W2"] += W2 * self.reg
        grads["W1"] += W1 * self.reg
        ############################################################################
        #                             END OF YOUR CODE                             #
        ############################################################################

        return loss, grads

class FullyConnectedNet(object):
    """
    A fully-connected neural network with an arbitrary number of hidden layers,
    ReLU nonlinearities, and a softmax loss function. This will also implement
    dropout and batch normalization as options. For a network with L layers,
    the architecture will be

    {affine - [batch norm] - relu - [dropout]} x (L - 1) - affine - softmax

    where batch normalization and dropout are optional, and the {...} block is
    repeated L - 1 times.

    Similar to the TwoLayerNet above, learnable parameters are stored in the
    self.params dictionary and will be learned using the Solver class.
    """

    def __init__(self, hidden_dims, input_dim=3 * 32 * 32, num_classes=10,
                 dropout=0, use_batchnorm=False, reg=0.0,
                 weight_scale=1e-2, dtype=np.float32, seed=None):
        """
        Initialize a new FullyConnectedNet.

        Inputs:
        - hidden_dims: A list of integers giving the size of each hidden layer.
        - input_dim: An integer giving the size of the input.
        - num_classes: An integer giving the number of classes to classify.
        - dropout: Scalar between 0 and 1 giving dropout strength. If dropout=0 then
          the network should not use dropout at all.
        - use_batchnorm: Whether or not the network should use batch normalization.
        - reg: Scalar giving L2 regularization strength.
        - weight_scale: Scalar giving the standard deviation for random
          initialization of the weights.
        - dtype: A numpy datatype object; all computations will be performed using
          this datatype. float32 is faster but less accurate, so you should use
          float64 for numeric gradient checking.
        - seed: If not None, then pass this random seed to the dropout layers. This
          will make the dropout layers deteriminstic so we can gradient check the
          model.
        """
        self.use_batchnorm = use_batchnorm
        self.use_dropout = dropout > 0
        self.reg = reg
        self.num_layers = 1 + len(hidden_dims)
        self.dtype = dtype
        self.params = {}

        ############################################################################
        # TODO: Initialize the parameters of the network, storing all values in    #
        # the self.params dictionary. Store weights and biases for the first layer #
        # in W1 and b1; for the second layer use W2 and b2, etc. Weights should be #
        # initialized from a normal distribution with standard deviation equal to  #
        # weight_scale and biases should be initialized to zero.                   #
        #                                                                          #
        # When using batch normalization, store scale and shift parameters for the #
        # first layer in gamma1 and beta1; for the second layer use gamma2 and     #
        # beta2, etc. Scale parameters should be initialized to one and shift      #
        # parameters should be initialized to zero.                                #
        ############################################################################
        # Find random weights (with mean 0) for the first layer based on the scale
        # defined while maintianing the right dimensions
        W1_layer1 = np.random.normal(0, weight_scale, input_dim * hidden_dims[0])
        # Store first layer weights in teh params dictionary, reshaping just to be sure
        self.params["W1"] = W1_layer1.reshape((input_dim, hidden_dims[0]))
        # Initialize first layer biases with zeros
        self.params['b1'] = np.zeros((hidden_dims[0]))

        # Loop through initialization for intermediary layers
        for i in xrange(1, self.num_layers - 1):
            # Initialize the indexed weights in the dictionary with random numbers with mean 0 and defined weight scale,
            # trying to maintain dimensionality
            self.params['W' + str(i + 1)] = np.random.normal(0, weight_scale,
                                                             hidden_dims[i - 1] * hidden_dims[i]).reshape(
                (hidden_dims[i - 1], hidden_dims[i]))
            # Initialize the indexed biases in the dictionary with zeros
            self.params['b' + str(i + 1)] = np.zeros((hidden_dims[i]))

            # Initialize the final layer weights with random numbers with mean 0 and defined weight scale.
        self.params['W' + str(self.num_layers)] = np.random.normal(0, weight_scale,
                                                                   hidden_dims[-1] * num_classes).reshape(
            (hidden_dims[-1], num_classes))
        # Initialize the final layer biases in the dictionary with zeros
        self.params['b' + str(self.num_layers)] = np.zeros((num_classes))

        # If batchnorm is being used, initialize the relevant parameters (gamma and beta).
        # Maintain the right dimensionality
        if self.use_batchnorm:
            for i in xrange(1, self.num_layers):
                self.params['beta' + str(i)] = np.zeros(hidden_dims[i - 1])
                self.params['gamma' + str(i)] = np.ones(hidden_dims[i - 1])

        ############################################################################
        #                             END OF YOUR CODE                             #
        ############################################################################

        # When using dropout we need to pass a dropout_param dictionary to each
        # dropout layer so that the layer knows the dropout probability and the mode
        # (train / test). You can pass the same dropout_param to each dropout layer.
        self.dropout_param = {}
        if self.use_dropout:
            self.dropout_param = {'mode': 'train', 'p': dropout}
            if seed is not None:
                self.dropout_param['seed'] = seed

        # With batch normalization we need to keep track of running means and
        # variances, so we need to pass a special bn_param object to each batch
        # normalization layer. You should pass self.bn_params[0] to the forward pass
        # of the first batch normalization layer, self.bn_params[1] to the forward
        # pass of the second batch normalization layer, etc.
        self.bn_params = []
        if self.use_batchnorm:
            self.bn_params = [{'mode': 'train'} for i in range(self.num_layers - 1)]

        # Cast all parameters to the correct datatype
        for k, v in self.params.items():
            self.params[k] = v.astype(dtype)

    def loss(self, X, y=None):
        """
        Compute loss and gradient for the fully-connected net.

        Input / output: Same as TwoLayerNet above.
        """
        X = X.astype(self.dtype)
        mode = 'test' if y is None else 'train'
        # Set train/test mode for batchnorm params and dropout param since they
        # behave differently during training and testing.
        if self.use_dropout:
            self.dropout_param['mode'] = mode
        if self.use_batchnorm:
            for bn_param in self.bn_params:
                bn_param['mode'] = mode
        scores = None
        ############################################################################
        # TODO: Implement the forward pass for the fully-connected net, computing  #
        # the class scores for X and storing them in the scores variable.          #
        #                                                                          #
        # When using dropout, you'll need to pass self.dropout_param to each       #
        # dropout forward pass.                                                    #
        #                                                                          #
        # When using batch normalization, you'll need to pass self.bn_params[0] to #
        # the forward pass for the first batch normalization layer, pass           #
        # self.bn_params[1] to the forward pass for the second batch normalization #
        # layer, etc.                                                              #
        ############################################################################
        # Container to store the cache data from each layer
        list_cache = []
        # Make a copy of the input data to make it easier to handle during iterations
        input_data = X.copy()
        # If dropout is being used
        if self.use_dropout:
            # Loop through the layers
            for i in xrange(1, self.num_layers):
                # Perfrom Affine Forward Pass
                a, fc_cache = affine_forward(input_data, self.params['W' + str(i)], self.params['b' + str(i)])
                # Perfrom Relu activation for forward pass
                r, relu_cache = relu_forward(a)
                # Use the ouptut from Relu activation to perform dropout forward pass
                dropout_output, dropout_cache = dropout_forward(r, self.dropout_param)
                # Store the cached values
                cache = (fc_cache, relu_cache, dropout_cache)
                # Store current cache to the Cache list
                list_cache.append(cache)
                # Set input data as current output for next layer
                input_data = dropout_output

            # Perfrom the affine forward pass for the last layer
            scores, dr_cache_ln = affine_forward(input_data, self.params['W' + str(self.num_layers)],
                                                 self.params['b' + str(self.num_layers)])
            # Store current cache to Cache list
            list_cache.append(dr_cache_ln)

        # If use of batchnorm is desired
        elif self.use_batchnorm:
            # Loop throught the layers
            for i in xrange(1, self.num_layers):
                # Perform Affine forward
                affine_out, fc_cache = affine_forward(input_data, self.params['W' + str(i)], self.params['b' + str(i)])
                # Perform Batchnorm first then Relu (Combo 1)
                bn_out, bn_cache = batchnorm_forward(affine_out, self.params['gamma' + str(i)],
                                                     self.params['beta' + str(i)], self.bn_params[i - 1])
                relu_out, relu_cache = relu_forward(bn_out)
                # Perform Relu activation & then Batchnorm (Combo 2) for Question 2
                # relu_out, relu_cache = relu_forward(affine_out)
                # bn_out, bn_cache = batchnorm_forward(relu_out, self.params['gamma'+str(i)], self.params['beta'+str(i)], self.bn_params[i-1])
                # Store all the cache
                batchnorm_cache = (fc_cache, bn_cache, relu_cache)
                list_cache.append(batchnorm_cache)
                input_data = relu_out

        # Otherwise run normal mode
        else:
            # Loop through the layers
            for i in xrange(1, self.num_layers):
                # Perform forward pass with Relu Activation at given indexed layer
                layer_output, layer_cache = affine_relu_forward(input_data, self.params['W' + str(i)],
                                                                self.params['b' + str(i)])
                # Store current cache to the cache list
                list_cache.append(layer_cache)
                # Set input data as current output for next layer
                input_data = layer_output

        # Perfrom the affine forward pass for the last layer
        scores, cache_ln = affine_forward(input_data, self.params['W' + str(self.num_layers)],
                                          self.params['b' + str(self.num_layers)])
        # Store current cache to Cache list
        list_cache.append(cache_ln)

        ############################################################################
        #                             END OF YOUR CODE                             #
        ############################################################################

        # If test mode return early
        if mode == 'test':
            return scores

        loss, grads = 0.0, {}
        ############################################################################
        # TODO: Implement the backward pass for the fully-connected net. Store the #
        # loss in the loss variable and gradients in the grads dictionary. Compute #
        # data loss using softmax, and make sure that grads[k] holds the gradients #
        # for self.params[k]. Don't forget to add L2 regularization!               #
        #                                                                          #
        # When using batch normalization, you don't need to regularize the scale   #
        # and shift parameters.                                                    #
        #                                                                          #
        # NOTE: To ensure that your implementation matches ours and you pass the   #
        # automated tests, make sure that your L2 regularization includes a factor #
        # of 0.5 to simplify the expression for the gradient.                      #
        ############################################################################
        # Evaluate the Softmax loss and derivative of scores
        loss, dout = softmax_loss(scores, y)
        # Initialize containers to store Derivatives of weights and biases
        list_dw = []
        list_db = []
        list_dgamma = []
        list_dbeta = []

        # Extract and remove the last layer entry in the Cache List
        cache = list_cache.pop()
        # Perform affine backward pass on the cached data for the last layer
        dx, dw, db = affine_backward(dout, cache)
        # Store the derivative of weights and biases at index 0 of respective lists
        list_dw.insert(0, dw)
        list_db.insert(0, db)
        dout = dx

        # If Dropout is being used
        if self.use_dropout:

            # Loop through the cached entries for all the intermediary layers
            for i in xrange(len(list_cache)):
                # Extract and remove the last entry in the cache list
                cache = list_cache.pop()
                # Extract the data from the cache file
                fc_cache, relu_cache, dropout_cache = cache
                # Perform dropout backward pass
                dd = dropout_backward(dout, dropout_cache)
                # Perfrom Relu activation for backward pass
                dr = relu_backward(dd, relu_cache)
                # Perfrom Affine backward Pass
                dx, dw, db = affine_backward(dr, fc_cache)
                # Update list of derivatives of weights and biases
                list_dw.insert(0, dw)
                list_db.insert(0, db)
                # Set derivative of output as derivative of x
                dout = dx

            para_loss = 0

            # Loop through the values in list of derivatives of weights
            for i in xrange(len(list_dw)):
                # Apply regularization to the weights
                W = self.params['W' + str(i + 1)]
                list_dw[i] += self.reg * W
                # Use para_loss variable to store the iterative penalty terms for the regularization
                para_loss += np.sum(W ** 2)
            # Regularize the loss
            loss += 0.5 * self.reg * para_loss

            # Loop through and update the grads dictionary entries for derivatives of weights and biases
            for i in xrange(len(list_dw)):
                grads['W' + str(i + 1)] = list_dw[i]
                grads['b' + str(i + 1)] = list_db[i]

        # If use of batchnorm is desired
        elif self.use_batchnorm:
            # Loop through the cached entries for all the intermediary layers
            for i in xrange(len(list_cache)):
                # Get the last entry from the cache list and store relevant layer caches
                cache = list_cache.pop()
                fc_cache, bn_cache, relu_cache = cache
                # Perform Relu Backward Pass then Batchnorm Backward (Combo1)
                drelu_out = relu_backward(dout, relu_cache)
                dbn_out, dgamma, dbeta = batchnorm_backward(drelu_out, bn_cache)
                dx, dw, db = affine_backward(dbn_out, fc_cache)
                # Perform Batchnorm Backward then Relu Backward Pass  (Combo2) for Question 2
                # dbn_out, dgamma, dbeta = batchnorm_backward(dout, bn_cache)
                # drelu_out = relu_backward(dbn_out, relu_cache)
                # dx, dw, db = affine_backward(drelu_out, fc_cache)
                # Store relevant gradients in the containers and update computed dx for next layer
                list_dw.insert(0, dw)
                list_db.insert(0, db)
                list_dgamma.insert(0, dgamma)
                list_dbeta.insert(0, dbeta)
                dout = dx

                # Regularize the loss and compute the parametric loss for all layers
            para_loss = 0
            for i in xrange(len(list_dw)):
                list_dw[i] += self.reg * self.params['W' + str(i + 1)]
                para_loss += np.sum((self.params['W' + str(i + 1)]) ** 2)
                # Store the gradients of Weights and biases for all the layers
                grads['W' + str(i + 1)] = list_dw[i]
                grads['b' + str(i + 1)] = list_db[i]

            # Calculate the loss and store sore respective gradients of Gamma and beta
            for i in xrange(len(list_dgamma)):
                list_dgamma[i] += self.reg * self.params['gamma' + str(i + 1)]
                list_dbeta[i] += self.reg * self.params['beta' + str(i + 1)]
                para_loss += np.sum((self.params['gamma' + str(i + 1)]) ** 2) + np.sum(
                    (self.params['beta' + str(i + 1)]) ** 2)
                grads['gamma' + str(i + 1)] = list_dgamma[i]
                grads['beta' + str(i + 1)] = list_dbeta[i]

            # Put it all together to find the total loss.
            loss += 0.5 * self.reg * para_loss

        # If dropout is not specified, run normal mode
        else:
            # Loop through the cached entries for all the intermediary layers
            for i in xrange(len(list_cache)):
                # Extract and remove the last entry in the cache list
                cache = list_cache.pop()
                # Perform Backward pass with Relu activation
                dx, dw, db = affine_relu_backward(dout, cache)
                # Update list of derivatives of weights and biases
                list_dw.insert(0, dw)
                list_db.insert(0, db)
                # Set derivative of output as derivative of x
                dout = dx
            para_loss = 0

            # Loop through the values in list of derivatives of weights
            for i in xrange(len(list_dw)):
                # Apply regularization to the weights
                W = self.params['W' + str(i + 1)]
                list_dw[i] += self.reg * W
                # Use para_loss variable to store the iterative penalty terms for the regularization
                para_loss += np.sum(W ** 2)
            # Regularize the loss
            loss += 0.5 * self.reg * para_loss

            # Loop through and update the grads dictionary entries for derivatives of weights and biases
            for i in xrange(len(list_dw)):
                grads['W' + str(i + 1)] = list_dw[i]
                grads['b' + str(i + 1)] = list_db[i]
        ############################################################################
        #                             END OF YOUR CODE                             #
        ############################################################################

        return loss, grads


text:
import torch
from torch import nn, optim
from copy import deepcopy
import numpy as np

class DynArr():
  def __init__(self, d=None, resizemult=2, initsize=0, initbuffersize=10, device='cpu', **kw):
    self.d = d
    self.device = device
    self.resizemult = resizemult
    if d is not None:
      self.arr = torch.zeros([initbuffersize, d], device=device)
    else:
      self.arr = torch.zeros([initbuffersize], device=device)
    self.size = initsize
  def isempty(self):
    return self.size == 0
  def _cat(self, arr):
    size = arr.shape[0]
    if self.size + size > len(self.arr):
      if self.d is not None:
        assert arr.shape[1] == self.d
        new_arr = torch.zeros([int((self.resizemult-1) * (self.size + size)), self.d], device=self.arr.device)
      else:
        assert len(arr.shape) == 1
        new_arr = torch.zeros([int((self.resizemult-1) * (self.size + size))], device=self.arr.device)
      self.arr = torch.cat([self.arr, new_arr])
    self.arr[self.size:self.size+size] = arr
    self.size += size
  def cat(self, *arrs):
    for arr in arrs:
      self._cat(arr)
  def checkpoint(self):
    self._checkpoint = self.size
  def restore(self):
    self.size = self._checkpoint
  @property
  def a(self):
    return self.arr[:self.size]
  def cuda(self):
    self.arr = self.arr.cuda()
    return self
  def cpu(self):
    self.arr = self.arr.cpu()
    return self
  def half(self):
    self.arr = self.arr.half()
    return self
  def float(self):
    self.arr = self.arr.float()
    return self

class CycArr():
  def __init__(self, d=None, maxsize=10000, initsize=0, device='cpu', **kw):
    assert initsize <= maxsize
    self.size = initsize
    if initsize == maxsize:
      self.end = 0
    else:
      self.end = None
    self.d = d
    self.maxsize = maxsize
    self.device = device
    if d is not None:
      self.arr = torch.zeros([maxsize, d], device=device)
    else:
      self.arr = torch.zeros([maxsize], device=device)

  @property
  def a(self):
    if self.size == self.maxsize:
      return self.arr
    return self.arr[:self.size]

  def cuda(self):
    self.arr = self.arr.cuda()
    return self

  def half(self):
    self.arr = self.arr.half()
    return self

  def isempty(self):
    return self.size == 0

  def _cat(self, arr):
    size = arr.shape[0]
    if self.size == self.maxsize:
      # cyclic writing
      if self.end + size < self.maxsize:
        self.arr[self.end:self.end+size] = arr
        self.end += size
      else:
        p1size = self.maxsize - self.end
        p2size = size - p1size
        self.arr[self.end:] = arr[:p1size]
        self.arr[:p2size] = arr[p1size:]
        self.end = p2size
    elif self.size + size >= self.maxsize:
      assert size < self.maxsize
      # writing at the end and spill over to beginning
      p1size = self.maxsize - self.size
      p2size = size - p1size
      self.arr[self.size:] = arr[:p1size]
      self.arr[:p2size] = arr[p1size:]
      self.end = p2size
      self.size = self.maxsize
    else:
      # noncyclic writing
      self.arr[self.size:self.size+size] = arr
      self.size += size

  def cat(self, *arrs):
    for arr in arrs:
      self._cat(arr)

ext: py

reason: import numpy


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

# @Time    : 19-6-20 下午4:18
# @Author  : zj

import numpy as np

def load_xor():
    xor_data = np.array([[0, 0], [0, 1], [1, 0], [1, 1]])
    xor_labels = np.array([0, 1, 1, 0])

    return xor_data, xor_labels

max_stars_repo_name: deep-learning-algorithm/PyNet

max_issues_repo_name: zjZSTU/PyNet

max_forks_repo_name: zjZSTU/PyNet


text:
from obspy import read
import numpy as np

def loadData(uploadedFile):
  if uploadedFile is not None:
    try:
      dataSGY = read(uploadedFile)
      nTrace = len(dataSGY)
      nLength = len(dataSGY[0].data)

      rawData = np.zeros((nLength, nTrace))
      axisYData = np.zeros((nLength, nTrace))

      maxTime = dataSGY[0].stats.npts
      time = np.cumsum([dataSGY[0].stats.delta for i in range(maxTime)])

      for i in range(len(dataSGY)):
        rawData[:,i] = dataSGY[i].data
        axisYData[:,i] = time

    except UnicodeDecodeError as e:
      statError = (f'error loading sgy data: {e}')
      return statError
  else:
    rawData = None
    axisYData = None

  return rawData, axisYData

max_stars_repo_name: putupradnya/streamlit-app

max_issues_repo_name: putupradnya/streamlit-app

max_forks_repo_name: putupradnya/streamlit-app


text:
"""Operations related to evaluating predictions.
"""

import numpy as np

def mse(y, y_pred):
    """Returns the mean squared error between
    ground truths and predictions.
    """
    return np.mean((y - y_pred) ** 2)

def rmse(y, y_pred):
    """Returns the root mean squared error between
    ground truths and predictions.
    """
    return np.sqrt(mse(y, y_pred))

def mae(y, y_pred):
    """Returns mean absolute error between
    ground truths and predictions.
    """
    return np.mean(np.fabs(y - y_pred))

max_stars_repo_name: zhrlove/tf-recsys

max_issues_repo_name: zhrlove/tf-recsys

max_forks_repo_name: zhrlove/tf-recsys


text:
import os
from utilities.fileOperations import FileOperations
import numpy as np

if __name__ == '__main__':
    fo = FileOperations()
    path = os.path.dirname(__file__)

    photos_path = os.path.join(path, '../yelpData/') + 'photos.json'
    photos_df = fo.getfile(photos_path).set_index('photo_id')

    # Images count
    print('Total images: ' + str(len(photos_df)))

    # Business count
    print('Total Businesses: ' + str(photos_df['business_id'].unique().size))

    # Duplicate images
    print('Duplicated images: ' + str(photos_df.reset_index().duplicated(subset=['photo_id']).any()))

    # Sanity check for missing yelpData
    print('Any images without labels: ' + str(photos_df['label'].isnull().any()))
    print('Any missing business id: ' + str(photos_df['business_id'].isnull().any()))

    print('-----------------------------------------')
    print('Distribution of images among 5 labels:')
    # Display image labels count
    print(photos_df['label'].value_counts().reset_index().
          rename(columns={'index': 'label', 'label': 'count'}))

    print('-----------------------------------------')
    photos_df['caption'].replace('', np.nan, inplace=True)
    photos_df.dropna(subset=['caption'], inplace=True)
    print('Number of photos with caption: ' + str(len(photos_df)))

max_stars_repo_name: naveenrc/YelpChallenge

max_issues_repo_name: naveenrc/YelpChallenge

max_forks_repo_name: naveenrc/YelpChallenge


text:
# Lint as: python3
# Copyright 2019 DeepMind Technologies Limited. All Rights Reserved.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
# ==============================================================================
"""Basic Haiku modules and functions."""

import functools
from typing import Any, Callable, Iterable, Optional, Text

from haiku._src import base
from haiku._src import initializers
from haiku._src import module
import jax
import jax.nn
import jax.numpy as jnp
import tree

# Utility and activation functions.
def one_hot(indices, num_classes, dtype=jnp.float32):
  """Returns a one-hot version of indices.

  Args:
    indices: A tensor of indices.
    num_classes: Number of classes in the one-hot dimension.
    dtype: The dtype.

  Returns:
    The one-hot tensor. If indices' shape is [A, B, ...], shape is
      [A, B, ... num_classes].
  """

  return jnp.array(
      indices[..., jnp.newaxis] == jnp.arange(num_classes), dtype=dtype)

def multinomial(rng, logits, num_samples):
  """Draws samples from a multinomial distribution.

  Args:
    rng: A JAX PRNGKey.
    logits: Unnormalized log-probabilities, of shape [batch_size, categories] or
      [categories].
    num_samples: Number of samples to draw.

  Returns:
    Chosen categories, of shape [batch_size, num_samples] or [num_samples].
  """
  # NOTE(tycai): Currently, tf.multinomial uses CDF for non-XLA CPU only.
  # We may want to switch to the Gumbel trick as used in XLA.
  if len(logits.shape) > 2 or not logits.shape:
    raise ValueError("Logits must be rank-1 or rank-2.")
  probs = jax.nn.softmax(logits)
  probs = jnp.cumsum(probs, axis=-1)
  # Special-case num_samples == 1 due to TPU padding, as in TF2XLA.
  # https://github.com/tensorflow/tensorflow/blob/b1608511d5a50d05825c4025b0c347e8689a241f/tensorflow/compiler/tf2xla/kernels/categorical_op.cc#L79
  if num_samples == 1:
    a = jax.random.uniform(rng, logits.shape[:-1] + (1,))
    out = jnp.argmin(a > probs, axis=-1)
    return out[..., None]
  else:
    a = jax.random.uniform(rng, (num_samples,) + logits.shape[:-1] + (1,))
    out = jnp.argmin(a > probs, axis=-1)
    return jnp.transpose(out)

# Common modules.
class Sequential(module.Module):
  """Sequentially calls the given list of layers.

  Note that `Sequential` is limited in the range of possible architectures
  it can handle. This is a deliberate design decision; `Sequential` is only
  meant to be used for the simple case of fusing together modules/ops where
  the input of a particular module/op is the output of the previous one.

  Another restriction is that it is not possible to have extra arguments in the
  `__call__` method that are passed to the constituents of the module - for
  example, if there is a `BatchNorm` module in `Sequential` and the user wishes
  to switch the `is_training` flag. If this is the desired use case, the
  recommended solution is to subclass `snt.Module` and implement `__call__`:

      >>> class CustomModule(hk.Module):
      ...   def __init__(self, name=None):
      ...     super().__init__(name=name)
      ...     self.conv2d = hk.Conv2D(32, 4, 2)
      ...     self.bn = hk.BatchNorm()
      ...
      ...   def __call__(self, inputs, is_training):
      ...     outputs = self.conv2d(inputs)
      ...     outputs = self.bn(outputs, is_training=is_training)
      ...     outputs = jax.nn.relu(outputs)
      ...     return outputs
  """

  def __init__(self,
               layers: Iterable[Callable[..., Any]] = None,
               name: Optional[Text] = None):
    super(Sequential, self).__init__(name=name)
    self.layers = layers

  def __call__(self, inputs, *args, **kwargs):
    """Connects all layers. *args and **kwargs are passed to the first layer."""
    outputs = inputs
    for i, layer in enumerate(self.layers):
      if i == 0:
        outputs = layer(outputs, *args, **kwargs)
      else:
        outputs = layer(outputs)
    return outputs

class Linear(module.Module):
  """Linear module."""

  def __init__(self,
               output_size: int,
               with_bias: bool = True,
               w_init: Optional[base.Initializer] = None,
               b_init: Optional[base.Initializer] = None,
               name: Optional[Text] = None):
    """Constructs the Linear module.

    Args:
      output_size: Output dimensionality.
      with_bias: Whether to add a bias to the output.
      w_init: Optional initializer for weights. By default, uses random values
        from truncated normal, with stddev `1 / sqrt(fan_in)`. See
        https://arxiv.org/abs/1502.03167v3.
      b_init: Optional initializer for bias. By default, zero.
      name: Name of the module.
    """
    super(Linear, self).__init__(name=name)
    self.input_size = None
    self.output_size = output_size
    self.with_bias = with_bias
    self.w_init = w_init
    self.b_init = b_init or jnp.zeros

  def __call__(self, inputs):
    if not inputs.shape:
      raise ValueError("Input must not be scalar.")

    self.input_size = inputs.shape[-1]
    default_stddev = 1. / jnp.sqrt(self.input_size)
    w_init = self.w_init or initializers.TruncatedNormal(stddev=default_stddev)

    w = base.get_parameter("w", [self.input_size, self.output_size],
                           inputs.dtype, init=w_init)
    out = jnp.dot(inputs, w)
    if self.with_bias:
      out += base.get_parameter("b", [self.output_size], inputs.dtype,
                                init=self.b_init)
    return out

def ndim_at_least(x, num_dims):
  if x is None:
    return False
  x = jnp.asarray(x)
  return len(x.shape) >= num_dims

def arbitrary_mergeable_leaf(min_num_dims, args, kwargs):
  for a in tree.flatten(args):
    if ndim_at_least(a, min_num_dims):
      return a
  for k in tree.flatten(kwargs):
    if ndim_at_least(k, min_num_dims):
      return k
  # Couldn't find a satisfactory leaf.
  return None

def merge_leading_dims(x, num_dims):
  """Merge leading dimensions."""
  # Don't merge if there aren't dimensions to merge.
  if not ndim_at_least(x, num_dims):
    return x

  new_shape = (jnp.prod(x.shape[:num_dims]),) + x.shape[num_dims:]
  return jnp.reshape(x, new_shape)

def split_leading_dim(x, to_dim):
  if x is None:
    return None
  new_shape = to_dim + x.shape[1:]
  return jnp.reshape(x, new_shape)

class BatchApply(object):
  """Temporarily merges leading dimensions of input tensors.

  Merges the leading dimensions of a tensor into a single dimension, runs the
  given callable, then splits the leading dimension of the result to match the
  input.

  Input arrays whose rank is smaller than the number of dimensions to collapse
  are passed unmodified.

  This may be useful for applying a module to each timestep of e.g. a
  [Time, Batch, ...] array.

  For some `f`s and platforms, this may be more efficient than `jax.vmap`,
  especially when combined with other transformations like `jax.grad`.
  """

  def __init__(self, f, num_dims=2):
    """Constructs a BatchApply module.

    Args:
      f: The callable to be applied to the reshaped array.
      num_dims: The number of dimensions to merge.
    """
    self._f = f
    self.num_dims = num_dims

  def __call__(self, *args, **kwargs):
    example = arbitrary_mergeable_leaf(self.num_dims, args, kwargs)
    if example is None:
      msg = "BatchApply requires at least one input with ndim >= {}."
      raise ValueError(msg.format(self.num_dims))
    merge = lambda x: merge_leading_dims(x, self.num_dims)
    split = lambda x: split_leading_dim(x, example.shape[:self.num_dims])
    args = tree.map_structure(merge, args)
    kwargs = tree.map_structure(merge, kwargs)
    outputs = self._f(*args, **kwargs)
    return tree.map_structure(split, outputs)

def expand_apply(f, axis=0):
  """Wraps f to temporarily add a size-1 axis to its inputs.

  Syntactic sugar for:
  ```
      in = tree.map_structure(lambda t: np.expand_dims(t, axis=axis), in)
      out = f(in)
      out = tree.map_structure(lambda t: np.squeeze(t, axis=axis), out)
  ```

  This may be useful for applying a function built for [Time, Batch, ...] arrays
  to a single timestep.

  Args:
    f: The callable to be applied to the expanded inputs.
    axis: Where to add the extra axis.

  Returns:
    f, wrapped as described above.
  """
  if axis not in [0, -1]:
    raise ValueError("expand_apply currently only supports axis=0 or axis=-1.")

  @functools.wraps(f)
  def wrapper(*args, **kwargs):
    expand = lambda t: jnp.expand_dims(t, axis=axis)
    args = tree.map_structure(expand, args)
    kwargs = tree.map_structure(expand, kwargs)
    outputs = f(*args, **kwargs)
    return tree.map_structure(lambda t: jnp.squeeze(t, axis=axis), outputs)

  return wrapper

def dropout(rng, rate, x):
  """Randomly drop units in the input at a given rate.

  See: http://www.cs.toronto.edu/~hinton/absps/dropout.pdf

  Args:
    rng: The RNGKey.
    rate: Probability that each element of x is discarded. Must be a scalar in
      the range `[0, 1)`.
    x: The value to be dropped out.
  Returns:
    x, but dropped out and scaled by 1 / (1 - rate).
  """
  if rate < 0 or rate >= 1:
    raise ValueError("rate must be in [0, 1).")
  keep_rate = 1.0 - rate
  keep = jax.random.bernoulli(rng, keep_rate, shape=x.shape)
  return keep * x / keep_rate

def to_module(f):
  """Converts a function into a module.

  Sample usage:
  >>> def add_bias(x):
  ...   b = hk.get_parameter("b", [], init=hk.initializers.RandomNormal())
  ...   return x + b
  >>> Bias = hk.to_module(add_bias)
  >>> def net(x, y):
  ...   b = Bias()
  ...   # Bias x and y by the same amount.
  ...   return b(x) * b(y)

  Args:
    f: The function to convert.

  Returns:
    A module class which runs runs `f` when called.
  """

  class Wrapper(module.Module):

    def __init__(self, name=None):
      if name is None:
        name = f.__name__
      super(Wrapper, self).__init__(name=name)

    def __call__(self, *a, **k):
      return f(*a, **k)

  return Wrapper

max_stars_repo_name: marload/dm-haiku

max_issues_repo_name: marload/dm-haiku

max_forks_repo_name: marload/dm-haiku

reason: import jax


text:
import os
import subprocess
import unittest

import pandas as pd
import numpy as np

import tempfile

from .mixins import FileMixins
from . import config

class TestTargeted(FileMixins, unittest.TestCase):
    def setUp(self):
        super(TestTargeted, self).setUp()
        self.output = os.path.join(self.out_dir, "targeted")

    def test_pyquant_trypsin(self):
        # This searches for the y1 ions of arginine. It also is a check that the label-scheme parameter works.
        f = tempfile.NamedTemporaryFile("w")
        f.write("\t".join(["0", "R", "10.008269", "R10"]) + "\n")
        f.write("\t".join(["1", "R", "0", "Light"]) + "\n")
        f.seek(0)
        com = [
            self.executable,
            "--scan-file",
            self.ecoli_mzml,
            "-p",
            str(config.CORES),
            "-o",
            self.output,
            "--msn-ion",
            "175",
            "--html",
            "--label-scheme",
            f.name,
        ]
        subprocess.call(com)
        pyquant = pd.read_table(self.output)
        label = "R10"
        pq_sel = "{}/Light".format(label)
        pyquant[pq_sel] = np.log2(pyquant[pq_sel] + 0.000001)
        r10_med = pyquant[pq_sel].median()
        self.assertNotAlmostEqual(r10_med, -2.0, places=2)

if __name__ == "__main__":
    unittest.main()


text:
import argparse
import networkx as nx
import numpy as np
import os
import pandas as pd
import multiprocessing as mp
import pickle
import re
from functools import partial
from itertools import combinations
import logging

if not __name__=="__main__":
    from . import get_HPO_similarity_score as gs

coding_variants = ["splice_acceptor_variant",
                   "splice_donor_variant",
                   "stop_gained",
                   "frameshift_variant",
                   "stop_lost",
                   "start_lost",
                   "inframe_insertion",
                   "inframe_deletion",
                   "missense_variant",
                   "protein_altering_variant",
                   "splice_region_variant",
                   "incomplete_terminal_codon_variant",
                   "start_retained_variant",
                   "stop_retained_variant",
                   "synonymous_variant",
                   "coding_sequence_variant"]
                   #"5_prime_UTR_variant",
                   #"3_prime_UTR_variant"]

logger = logging.getLogger(__name__)

cadd_identifier = "CADD_PHRED"
duplication_identifier = "segmentDuplication"
repeat_identifier = "simpleRepeat"
low_conf_region_identifier = "lowConfRegion"

def prioritize_variants(variant_data, hpo_resources_folder, num_cores, family_file=None, family_type="SINGLE", hpo_list=None, gene_exclusion_list=None):
    # load HPO resources
    gene_2_HPO_f = hpo_resources_folder + "gene2hpo.pkl"
    hgnc_2_gene_f = hpo_resources_folder + "hgnc2gene.pkl"
    gene_2_interacting_f = hpo_resources_folder + "gene2interacting.pkl"
    HPO_graph_file = hpo_resources_folder + "hpo_graph.pkl"
    hpo_list_file = hpo_list
    gene_exclusion_file = gene_exclusion_list

    hgnc_2_gene = pickle.load(open(hgnc_2_gene_f, "rb"))
    gene_2_interacting = pickle.load(open(gene_2_interacting_f, "rb"))
    gene_2_HPO = pickle.load(open(gene_2_HPO_f, "rb"))
    hpo_nodes, hpo_edges = pickle.load(open(HPO_graph_file, "rb"))

    HPO_graph = nx.Graph()
    HPO_graph.add_nodes_from(hpo_nodes)
    HPO_graph.add_edges_from(hpo_edges)

    genes2exclude = set()
    if gene_exclusion_file is not None:
        if os.path.isfile(gene_exclusion_file):
            with open(gene_exclusion_file, "r") as exclusion_file:
                for line in exclusion_file:
                    if line.startswith("#"):
                        continue
                    if line == "\n":
                        continue
                    gene = line.rstrip().split("\t")[0].upper()
                    genes2exclude.add(gene)
        else:
            logger.error("The specified gene exclusion list %s is not a valid file" % (gene_exclusion_file))
            logger.warn("No genes will be excluded during filtering!")

    HPO_query = set()
    HPO_query_distances = 0
    if hpo_list_file is not None:
        if os.path.isfile(hpo_list_file):
            with open(hpo_list_file, "r") as hpo_file:
                for line in hpo_file:
                    HPO_term = line.rstrip()
                    HPO_query.add(HPO_term)
            HPO_query = list(HPO_query) # removes duplicate entries in the list
            HPO_query.sort() # makes sure that the gene symbols are ordered (could lead to problems otherwise)
            HPO_query_distances = gs.precompute_query_distances(HPO_graph, HPO_query)
        else:
            logger.error("The specified HPO list %s is not a valid file" % (hpo_list_file))
            logger.warn("HPO score finalization will be skipped!")

    # read family relationship (from PED file)
    family = dict()
    if family_file is not None:
        if os.path.isfile(family_file):
            with open(family_file, "r") as fam_file:
                for line in fam_file:
                    line = line.rstrip()
                    splitline = line.split("\t")

                    if splitline[5] == "2":
                        family[splitline[1]] = 1
                    elif splitline[5] == "1":
                        family[splitline[1]] = 0
                    else:
                        logger.error("There was a problem with the given PED file describing the family relations.")
        else:
            logger.error("The specified family file %s is not a valid file" % (family_file))
            logger.warn("Inheritance assessment will be skipped!")

    variant_data = parallelize_dataframe_processing(variant_data, partial(parallelized_variant_processing, family, family_type, genes2exclude, gene_2_HPO, hgnc_2_gene, gene_2_interacting, HPO_graph, HPO_query, HPO_query_distances), num_cores)
    variant_data = variant_data.sort_values(["FINAL_AIDIVA_SCORE"], ascending=[False])
    variant_data = variant_data.reset_index(drop=True)

    return variant_data

def parallelize_dataframe_processing(variant_data, function, num_cores):
    num_partitions = num_cores * 2

    if len(variant_data) <= num_partitions:
        dataframe_splitted = np.array_split(variant_data, 1)
    else:
        dataframe_splitted = np.array_split(variant_data, num_partitions)

    try:
        pool = mp.Pool(num_cores)
        variant_data = pd.concat(pool.map(function, dataframe_splitted))
    finally:
        pool.close()
        pool.join()

    return variant_data

def parallelized_variant_processing(family, family_type, genes2exclude, gene_2_HPO, hgnc_2_gene, gene_2_interacting, HPO_graph, HPO_query, HPO_query_distances, variant_data):
    variant_data = check_inheritance(variant_data, family_type, family)
    variant_data[["HPO_RELATEDNESS", "HPO_RELATEDNESS_INTERACTING", "FINAL_AIDIVA_SCORE"]] = variant_data.apply(lambda variant: pd.Series(compute_hpo_relatedness_and_final_score(variant, genes2exclude, gene_2_HPO, hgnc_2_gene, gene_2_interacting, HPO_graph, HPO_query, HPO_query_distances)), axis=1)
    variant_data[["FILTER_PASSED", "FILTER_COMMENT"]] = variant_data.apply(lambda variant: pd.Series(check_filters(variant, genes2exclude, HPO_query)), axis=1)

    return variant_data

def compute_hpo_relatedness_and_final_score(variant, genes2exclude, gene_2_HPO, hgnc_2_gene, gene_2_interacting, HPO_graph, HPO_query, HPO_query_distances):
    if HPO_query:
        if np.isnan(variant["AIDIVA_SCORE"]) and ((str(variant["rf_score"]) == "nan") or (str(variant["rf_score"]) == "")) and ((str(variant["ada_score"]) == "nan") or (str(variant["ada_score"]) == "")):
            final_score = np.nan
            hpo_relatedness = np.nan
            hpo_relatedness_interacting = np.nan

        else:
            if np.isnan(variant["AIDIVA_SCORE"]):
                # if both scores are present use maximum to benefit of the strengths of each model
                if ((str(variant["rf_score"]) != "nan") and (str(variant["rf_score"]) != "")) and ((str(variant["ada_score"]) != "nan") and (str(variant["ada_score"]) != "")):
                    pathogenictiy_prediction = max([float(variant["rf_score"]), float(variant["ada_score"])], default=np.nan)
                elif ((str(variant["rf_score"]) != "") and (str(variant["rf_score"]) != "nan")) and ((str(variant["ada_score"]) == "") or (str(variant["ada_score"]) == "nan")):
                    pathogenictiy_prediction = float(variant["rf_score"])
                elif ((str(variant["ada_score"]) != "") and (str(variant["ada_score"]) != "nan")) and ((str(variant["rf_score"]) == "") or (str(variant["rf_score"]) == "nan")):
                    pathogenictiy_prediction = float(variant["ada_score"])
                else:
                    pathogenictiy_prediction = np.nan

                logger.info("Using dbscSNV prediction instead of AIdiva prediction for splicing variant!")

            else:
                pathogenictiy_prediction = float(variant["AIDIVA_SCORE"])

            variant_gene = str(variant["SYMBOL"]).upper()
            hgnc_id = str(variant["HGNC_ID"])
            gene_distances = []
            gene_distances_interacting = []
            processed_HPO_genes = dict()

            if variant_gene in gene_2_interacting.keys():
                interacting_genes = [gene_interaction["interacting_gene"] for gene_interaction in gene_2_interacting[variant_gene]]
            else:
                interacting_genes = []

            # we use the hgnc ID to prevent problems if a given gene symbol isn't used anymore
            # let variants with a high AIDIVA_SCORE (>=0.8) pass to prevent missed hits
            if (variant_gene not in genes2exclude) or (pathogenictiy_prediction >= 0.8):
                if variant_gene in gene_2_HPO.keys():
                    gene_HPO_list = gene_2_HPO.get(variant_gene, [])

                else:
                    if (str(hgnc_id) != "nan") and (hgnc_id in hgnc_2_gene.keys()):
                        gene_symbol = hgnc_2_gene[hgnc_id]
                        gene_HPO_list = gene_2_HPO.get(gene_symbol, [])

                    else:
                        logger.warn("The processed gene %s (HGNC:%s) is not found in the current HPO resources" % (variant_gene, hgnc_id))
                        gene_HPO_list = []

                g_dist = gs.list_distance(HPO_graph, HPO_query, gene_HPO_list, HPO_query_distances)
                gene_distances.append(g_dist)

                # TODO: check with HGNC ID to prevent use of previous gene symbols
                for interacting_gene in interacting_genes:
                    if (interacting_gene in genes2exclude) and (pathogenictiy_prediction < 0.8):
                        continue

                    if interacting_gene in gene_2_HPO.keys():
                        gene_HPO_list = gene_2_HPO.get(interacting_gene, [])

                    else:
                        logger.warn("The processed interacting gene %s is not found in the current HPO resources" % (interacting_gene))
                        gene_HPO_list = []

                    g_dist = gs.list_distance(HPO_graph, HPO_query, gene_HPO_list, HPO_query_distances)
                    gene_distances_interacting.append(g_dist)

                if gene_distances or gene_distances_interacting:
                    hpo_relatedness = max(gene_distances, default=0.0)
                    hpo_relatedness_interacting = max(gene_distances_interacting, default=0.0)

                    ## TODO: try different weighting of AIDIVA_SCORE and HPO_RELATEDNESS and HPO_RELATEDNESS_INTERACTING (eg 0.6 and 0.3 and 0.1)
                    # predicted pathogenicity has a higher weight than the HPO relatedness
                    final_score = (pathogenictiy_prediction * 0.6 + float(hpo_relatedness) * 0.3 + float(hpo_relatedness_interacting) * 0.1)

            else:
                final_score = np.nan
                hpo_relatedness = np.nan
                hpo_relatedness_interacting = np.nan

    else:
        if np.isnan(variant["AIDIVA_SCORE"]):
            if ((str(variant["rf_score"]) == "") or (str(variant["rf_score"]) == "nan")) and ((str(variant["ada_score"]) == "") or (str(variant["ada_score"]) == "nan")):
                pathogenictiy_prediction = np.nan

            else:
                if (str(variant["rf_score"]) != "") and (str(variant["rf_score"]) != "nan"):
                    rf_score = float(variant["rf_score"])

                else:
                    rf_score = np.nan

                if (str(variant["ada_score"]) != "") and (str(variant["ada_score"]) != "nan"):
                    ada_score = float(variant["ada_score"])

                else:
                    ada_score = np.nan

                pathogenictiy_prediction = max([rf_score, ada_score], default=np.nan)
                logger.info("Using dbscSNV prediction instead of AIdiva prediction for splicing variant!")

        else:
            pathogenictiy_prediction = float(variant["AIDIVA_SCORE"])

        final_score = pathogenictiy_prediction
        hpo_relatedness = np.nan
        hpo_relatedness_interacting = np.nan

    return [hpo_relatedness, hpo_relatedness_interacting, final_score]

def check_inheritance(variant_data, family_type="SINGLE", family=None):
    variant_columns = variant_data.columns
    if family_type == "SINGLE":
        variant_data["COMPOUND"] = 0
        variant_data["RECESSIVE"] = variant_data.apply(lambda variant: check_recessive_single(variant, variant_columns), axis=1)

        variant_data_grouped = [group for key, group in variant_data.groupby("SYMBOL")]

        for group in variant_data_grouped:
            check_compound_single(group, variant_columns)
        variant_data = pd.concat(variant_data_grouped)

    elif family_type == "TRIO":
        if not family is None:
            variant_data["COMPOUND"] = 0
            variant_data["DOMINANT_DENOVO"] = variant_data.apply(lambda variant: check_denovo(variant, family), axis=1)

            # TODO: do we need a check for affected family members?
            variant_data["DOMINANT"] = variant_data.apply(lambda variant: check_dominant(variant, family), axis=1)
            variant_data["XLINKED"] = variant_data.apply(lambda variant: check_xlinked(variant, family), axis=1)
            variant_data["RECESSIVE"] = variant_data.apply(lambda variant: check_recessive(variant, family, family_type), axis=1)

            variant_data_grouped = [group for key, group in variant_data.groupby("SYMBOL")]

            affected_child = ""
            parent_1 = ""
            parent_2 = ""

            # TODO: handle the case if affected status is given as unknown
            # TODO: change to be less strict and to not rely on unaffected family members
            for name in family.keys():
                if family[name] == 1:
                    affected_child = name
                elif family[name] == 0:
                    if not parent_1:
                        parent_1 = name
                        continue
                    elif not parent_2:
                        parent_2 = name
                        continue
                    else:
                        logger.error("There was a problem processing the loaded PED file")

            if  affected_child and parent_1 and parent_2:
                for group in variant_data_grouped:
                    check_compound(group, affected_child, parent_1, parent_2)

            variant_data = pd.concat(variant_data_grouped)

        else:
            logger.error("If family type (TRIO) is used a proper PED file defining the family relations is required!")
            logger.warn("Inheritance assessment will be skipped!")

    elif (family_type == "MULTI") and (family is not None):
        if family is not None:
            variant_data["DOMINANT"] = variant_data.apply(lambda variant: check_dominant(variant, family), axis=1)
            variant_data["XLINKED"] = variant_data.apply(lambda variant: check_xlinked(variant, family), axis=1)
            variant_data["RECESSIVE"] = variant_data.apply(lambda variant: check_recessive(variant, family, family_type), axis=1)

        else:
            logger.error("If family type (MULTI) is used a proper PED file defining the family relations (if known) is required!")
            logger.warn("Inheritance assessment will be skipped!")

    else:
        logger.error("Unsupported family type (%s) was given!" % family_type)
        logger.warn("Inheritance assessment will be skipped!")

    variant_columns = variant_data.columns
    variant_data["INHERITANCE"] = variant_data.apply(lambda variant: add_inheritance_mode(variant, variant_columns), axis=1)

    return variant_data

def add_inheritance_mode(variant, variant_columns):
    inheritance_list = []

    if "DOMINANT" in variant_columns:
        if variant["DOMINANT"] == 1:
            inheritance_list.append("DOMINANT")

    if "DOMINANT_DENOVO" in variant_columns:
        if variant["DOMINANT_DENOVO"] == 1:
            inheritance_list.append("DOMINANT_DENOVO")

    if "RECESSIVE" in variant_columns:
        if variant["RECESSIVE"] == 1:
            inheritance_list.append("RECESSIVE")

    if "COMPOUND" in variant_columns:
        if variant["COMPOUND"] == 1:
            inheritance_list.append("COMPOUND")

    if "XLINKED" in variant_columns:
        if variant["XLINKED"] == 1:
            inheritance_list.append("XLINKED")

    inheritance_mode = "&".join(inheritance_list)

    return inheritance_mode

def check_filters(variant, genes2exclude, HPO_query):
    variant_genes = re.sub("\(.*?\)", "", str(variant["SYMBOL"]))
    genenames = set(variant_genes.split(";"))

    consequences = str(variant["Consequence"])
    found_consequences = consequences.split("&")
    seg_dup = float(variant[duplication_identifier])
    repeat = str(variant[repeat_identifier])
    cadd = float(variant[cadd_identifier])
    filter_comment = ""

    try:
        low_conf_region = int(variant[low_conf_region_identifier])
    except Exception as e:
        low_conf_region = 0

    try:
        maf = float(variant["MAX_AF"])
    except Exception as e:
        logger.warn("Allele frequency entry could not be identified, use 0.0 instead")
        maf = 0.0

    # filter for low confidence regions (these are regions where often false positives were observed)
    if "low_conf_region" in variant["FILTER"]:
        filter_passed = 0 # low confidence region
        filter_comment = "low confidence region"

        return filter_passed, filter_comment

    elif "off-target" in variant["FILTER"]:
        filter_passed = 0 # off target
        filter_comment = "off target"

        return filter_passed, filter_comment

    # high confidence: 0 - 0.15
    # mid confidence: 0.15 - 0.75
    # see https://onlinelibrary.wiley.com/doi/full/10.1002/humu.23674 for more detailed information about the abb score
    # if high confidence (> 0.15) filtering is to strict mid confidence (> 0.75) filtering could be applied
    if float(variant["ABB_SCORE"]) > 0.75:
        filter_passed = 0 # abb score to high -> crap region (low confidence variant)
        filter_comment = "ABB_SCORE"

        return filter_passed, filter_comment

    # exclude gene, if it is in the exclusion list
    if len(genes2exclude & genenames) > 0:
        for gene in genenames:
            # let variants with a high AIDIVA_SCORE (>=0.8) pass to be more sensitive
            if (gene.upper() in genes2exclude) and (float(variant["AIDIVA_SCORE"]) < 0.8):
                filter_passed = 0 # gene in exclusion list
                filter_comment = "gene exclusion"

                return filter_passed, filter_comment

    # let variants with a high AIDIVA_SCORE (>=0.8) pass to be more sensitive
    if (repeat != "NA") and (repeat != "") and (repeat != "nan") and (float(variant["AIDIVA_SCORE"]) < 0.8):
        filter_passed = 0 # tandem repeat
        filter_comment = "tandem repeat"

        return filter_passed, filter_comment

    # let variants with a high AIDIVA_SCORE (>=0.8) pass to be more sensitive
    if (seg_dup > 0) and (float(variant["AIDIVA_SCORE"]) < 0.8):
        filter_passed = 0 # segmental duplication
        filter_comment = "segmental duplication"

        return filter_passed, filter_comment

    ## TODO: change frequency based on inheritance mode (hom/het)
    if maf <= 0.01:
        if any(term for term in coding_variants if term in found_consequences):
            if not np.isnan(variant["FINAL_AIDIVA_SCORE"]):
                if len(HPO_query) >= 1:
                    if float(variant["HPO_RELATEDNESS"]) > 0.0:
                        filter_passed = 1
                        filter_comment = "passed all"

                    elif float(variant["HPO_RELATEDNESS_INTERACTING"]) > 0.0:
                        filter_passed = 1
                        filter_comment = "HPO related to interacting genes"

                    else:
                        filter_passed = 0 # no relation to reported HPO terms
                        filter_comment = "no HPO relation"

                else:
                    filter_passed = 1 # skip hpo filter if no terms are present
                    filter_comment = "no HPO terms given"

            else:
                filter_passed = 0 # no prediction present (eg. variant type not covered by the used ML models)
                filter_comment = "missing AIDIVA_SCORE"

        else:
            filter_passed = 0 # non coding
            filter_comment = "non coding"

    else:
        filter_passed = 0 # allele frequency to high
        filter_comment = "high MAF"

    return filter_passed, filter_comment

def check_compound(gene_variants, affected_child, parent_1, parent_2):
    num_variant_candidates = gene_variants.shape[0]

    if num_variant_candidates >= 2:
        candidate_indices = [x for x in combinations(gene_variants.index.tolist(), 2)]
        for index_pair in candidate_indices:
            affected_child_zygosity_a = gene_variants.loc[index_pair[0], "GT." + affected_child].replace("|", "/")
            affected_child_zygosity_b = gene_variants.loc[index_pair[1], "GT." + affected_child].replace("|", "/")

            parent_1_zygosity_a = gene_variants.loc[index_pair[0], "GT." + parent_1].replace("|", "/")
            parent_1_zygosity_b = gene_variants.loc[index_pair[1], "GT." + parent_1].replace("|", "/")

            parent_2_zygosity_a = gene_variants.loc[index_pair[0], "GT." + parent_2].replace("|", "/")
            parent_2_zygosity_b = gene_variants.loc[index_pair[1], "GT." + parent_2].replace("|", "/")

            if (affected_child_zygosity_a == "0/1") and (affected_child_zygosity_b == "0/1"):
                if ((parent_1_zygosity_a == "0/0") and (parent_2_zygosity_a == "0/1")) and ((parent_1_zygosity_b == "0/1") and (parent_2_zygosity_b == "0/0")):
                    gene_variants.loc[index_pair[0], "COMPOUND"] = 1
                    gene_variants.loc[index_pair[1], "COMPOUND"] = 1

                elif ((parent_1_zygosity_a == "0/1") and (parent_2_zygosity_a == "0/0")) and ((parent_1_zygosity_b == "0/0") and (parent_2_zygosity_b == "0/1")):
                    gene_variants.loc[index_pair[0], "COMPOUND"] = 1
                    gene_variants.loc[index_pair[1], "COMPOUND"] = 1

            if ("." in affected_child_zygosity_a) or ("." in affected_child_zygosity_b):
                logger.warn("Skip variant pair, uncalled genotype in affected sample!")

def check_compound_single(gene_variants, variant_columns):
    genotype_column = [column for column in variant_columns if column.startswith("GT.")][0]
    num_variant_candidates = gene_variants.shape[0]

    if num_variant_candidates >= 2:
        candidate_indices = [x for x in combinations(gene_variants.index.tolist(), 2)]
        for index_pair in candidate_indices:
            genotype_variant_a = gene_variants.loc[index_pair[0], genotype_column].replace("|", "/")
            genotype_variant_b = gene_variants.loc[index_pair[1], genotype_column].replace("|", "/")

            if (((genotype_variant_a == "0/1") and (genotype_variant_b == "0/1"))):
                gene_variants.loc[index_pair[0], "COMPOUND"] = 1
                gene_variants.loc[index_pair[1], "COMPOUND"] = 1

def check_denovo(variant, family):
    judgement = 0
    check_samples = dict()

    # create data structure for completeness check
    for name in family.keys():
        check_samples[name] = 0

        zygosity = variant["GT." + name].replace("|", "/")

        # check if sample is found in pedigree
        # sample info complete?
        ## TODO: do we need this check???
        if name in check_samples:
            check_samples[name] = 1

        # heterozygous in affected individual - good
        if zygosity == "0/1" and family[name] == 1:
            judgement = 1
            continue

        # hom ref, not affected - good
        elif zygosity == "0/0" and family[name] == 0 :
            judgement = 1
            continue

        # heterozygous in non-affected - bad
        elif zygosity == "0/1" and family[name] == 0:
            judgement = 0
            break

        # hom ref in affected - bad
        elif zygosity == "0/0" and family[name] == 1:
            judgement = 0
            break

        # homozygous can"t be denovo
        elif zygosity == "1/1":
            judgement = 0
            break

        else:
            # reject missing genotype here or beforehand
            pass

    for vals in check_samples.values():
       if vals == 0:
            judgement = 0
            break

    return judgement

def check_dominant(variant, family):
    judgement = 0
    check_samples = dict()

    # create data structure for completeness check
    for name in family.keys():
        check_samples[name] = 0

        zygosity = variant["GT." + name].replace("|", "/")

        if name in check_samples:
            check_samples[name] = 1

        # affected family members should have the mutation (hom ref not allowed)
        if zygosity == "0/0" and family[name] == 1:
            judgement = 0
            break

        # affected family members might be het
        elif zygosity == "0/1" and family[name] == 1:
            judgement = 1
            continue

        # affected family members might be hom alt
        # that"s the major difference to de novo...
        elif zygosity == "1/1" and family[name] == 1:
            judgement = 1
            continue

        # non-affected family members must not have the mutation - hom ref is OK
        elif zygosity == "0/0" and family[name] == 0:
            judgement = 1
            continue

        # non-affected family members must not have the mutation - het is bad
        elif zygosity == "0/1" and family[name] == 0:
            judgement = 0
            break

        # non-affected family members must not have the mutation - hom alt is worst
        elif zygosity == "1/1" and family[name] == 0:
            judgement = 0
            break

        else:
            # reject missing genotype here or beforehand
            pass

    for vals in check_samples.values():
       if vals == 0:
            judgement = 0
            break

    return judgement

def check_dominant_single(variant, variant_columns):
    genotype_column = [column for column in variant_columns if column.startswith("GT.")][0]
    judgement = 0

    if (variant[genotype_column] == "0/1"):
                judgement = 1

    return judgement

def check_recessive(variant, family, family_type):
    judgement = 0
    check_samples = dict()

    # create data structure for completeness check
    for name in family.keys():
        check_samples[name] = 0

        zygosity = variant["GT." + name].replace("|", "/")

        if name in check_samples:
            check_samples[name] = 1

        # affected individuals have to be homozygous
        if zygosity == "1/1" and family[name] == 1:
            judgement = 1
            continue

        # affected individuals should not be hom ref or het
        elif ( zygosity == "0/0" or zygosity == "0/1" ) and family[name] == 1:
            judgement = 0
            break

        # non-affected individuals might be het
        elif zygosity == "0/1" and family[name] == 0:
            judgement = 1
            continue

        # non-affected individuals might be hom ref, if a family is interrogated
        elif zygosity == "0/0" and family[name] == 0 and family_type == "MULTI":
            judgement = 1
            continue

        # non-affected individuals in a trio are the parents and have to be het
        elif zygosity == "0/0" and family[name] == 0 and family_type == "TRIO":
            judgement = 0
            break

        # non-affected individuals must not be hom alt
        elif zygosity == "1/1" and family[name] == 0:
            judgement = 0
            break

        else:
            # reject missing genotype here or beforehand
            pass

    for vals in check_samples.values():
        if vals == 0:
            judgement = 0
            break

    return judgement

def check_recessive_single(variant, variant_columns):
    genotype_column = [column for column in variant_columns if column.startswith("GT.")][0]
    variant_genotype = variant[genotype_column].replace("|", "/")
    is_recessive = 0

    # recessive variants have to be homozygous in affected samples
    if (variant_genotype == "1/1"):
        is_recessive = 1

    return is_recessive

def check_xlinked(variant, family):
    judgement = 0
    check_samples = dict()
    inheritance_logic = dict()

    if not ((variant["CHROM"] == "X") or (variant["CHROM"] == "x") or (variant["CHROM"] == "chrX") or (variant["CHROM"] == "chrx") or (variant["CHROM"] == "23")):
        return 0

    # create data structure for completeness check
    for name in family.keys():
        check_samples[name] = 0

        zygosity = variant["GT." + name].replace("|", "/")

        if name in check_samples:
            check_samples[name] = 1

        if family[name] == 0:
            inheritance_logic[name] = zygosity

        # affected individuals have to be homozygous
        if zygosity == "1/1" and family[name] == 1:
            judgement = 1
            continue

        # affected individuals should not be hom ref or het
        elif ( zygosity == "0/0" or zygosity == "0/1" ) and family[name] == 1:
            judgement = 0
            break

        # non-affected individuals might be het
        elif zygosity == "0/1" and family[name] == 0:
            judgement = 1
            continue

        # non-affected individuals might be hom ref
        elif zygosity == "0/0" and family[name] == 0:
            judgement = 1
            continue

        # non-affected individuals must not be hom alt
        elif zygosity == "1/1" and family[name] == 0:
            judgement = 0
            break

        else:
            # reject missing genotype here or beforehand
            pass

    # sanity check
    het_checker = 0
    hom_checker = 0

    for values in inheritance_logic.values():
        # mother
        if values == "0/1":
            het_checker = 1

        # father
        if values == "0/0":
            hom_checker = 1

    if het_checker == 1 and hom_checker == 1:
        judgement = 1
    else:
        judgement = 0

    # another sanity check
    for vals in check_samples.values():
        if vals == 0:
            judgement = 0
            break

    return judgement

if __name__=="__main__":
    import get_HPO_similarity_score as gs

    parser = argparse.ArgumentParser(description = "Filter variants and finalize the AIDIVA_SCORE based on the given HPO terms (if this information is present)")
    parser.add_argument("--in_file", type=str, dest="in_file", required=True, help="Tab separated input annotated and scored file [required]")
    parser.add_argument("--out_file", type=str, dest="out_filename", required=True, help="Name to save the results [required]")
    parser.add_argument("--family", type=str, dest="family", required=False, help="Tab separated list of samples annotated with affection status. [required]")
    parser.add_argument("--family_type", type=str, choices=["TRIO", "FAMILY", "SINGLE"], dest="family_type", required=False, help="Choose if the data you provide is a trio or a larger family [required]")
    parser.add_argument("--gene_exclusion", type=str, dest="gene_exclusion_list", required=False, help="List of genes that should be excluded in the prioritization")
    parser.add_argument("--hpo_list", type=str, dest="hpo_list", default=None, required=False, help="List of HPO terms that are observed in the patient. These terms are used to adjust the AIDIVA_SCORE\n")
    parser.add_argument("--hpo_resources", type=str, dest="hpo_resources", default="../../data/", required=True, help="Folder where the HPO resources (HPO_graph,...) are found\n")
    args = parser.parse_args()

    input_data = pd.read_csv(args.in_file, sep="\t", low_memory=False)

    prioritized_variants = prioritize_variants(input_data, args.hpo_resources, args.family, args.family_type, args.hpo_list, args.gene_exclusion_list, args.low_conf_region)
    prioritized_variants.to_csv(args.out_filename, sep="\t", index=False)

max_stars_repo_name: imgag/AIdiva

max_issues_repo_name: imgag/AIdiva

max_forks_repo_name: imgag/AIdiva

reason: import numpy,import networkx


text:
# Author: xiaotaw@qq.com (Any bug report is welcome)
# Time Created: Nov 2016
# Time Last Updated: Dec 2016
# Addr: Shenzhen, China
# Description:

import os
import getpass
import numpy as np
import pandas as pd
from scipy import sparse
from collections import defaultdict

data_dir = "/home/%s/Documents/chembl/data_files/" % getpass.getuser()
fp_dir = os.path.join(data_dir, "fp_files")
mask_dir = os.path.join(data_dir, "mask_files")
structure_dir = os.path.join(data_dir, "structure_files")

# the newly picked out 15 targets, include 9 targets from 5 big group, and 6 targets from others.
target_list = ["CHEMBL279", "CHEMBL203", # Protein Kinases
               "CHEMBL217", "CHEMBL253", # GPCRs (Family A)
               "CHEMBL235", "CHEMBL206", # Nuclear Hormone Receptors
               "CHEMBL240", "CHEMBL4296", # Voltage Gated Ion Channels
               "CHEMBL4805", # Ligand Gated Ion Channels
               "CHEMBL204", "CHEMBL244", "CHEMBL4822", "CHEMBL340", "CHEMBL205", "CHEMBL4005" # Others
              ]

def dense_to_one_hot(labels_dense, num_classes=2, dtype=np.int):
  """Convert class labels from scalars to one-hot vectors.
  Args:
    labels_dense: <type 'numpy.ndarray'> dense label
    num_classes: <type 'int'> the number of classes in one hot label
    dtype: <type 'type'> data type
  Return:
    labels_ont_hot: <type 'numpy.ndarray'> one hot label
  """
  num_labels = labels_dense.shape[0]
  index_offset = np.arange(num_labels) * num_classes
  labels_one_hot = np.zeros((num_labels, num_classes))
  labels_one_hot.flat[index_offset + labels_dense.ravel().astype(dtype)] = 1
  return labels_one_hot

def sparse_features(fps_list, target_columns_dict, num_features, is_log=True):
  """construct a sparse matrix(csr_matrix) for features according to target_columns_dict.
  Args:
    fps_list: <type 'list'> a list of apfps for the molecules
    is_log: <type 'bool'> flag whether apply np.log to data, default is True
  Return:
    features: the sparse matrix of features
  """
  data = []
  indices = []
  indptr = [0]
  for fps_str in fps_list:
    n = indptr[-1]
    for fp in fps_str[1:-1].split(","):
      if ":" in fp:
        k, v = fp.split(":")
        indices.append(target_columns_dict[int(k)])
        data.append(int(v))
        n += 1
    indptr.append(n)
  data = np.array(data)
  if is_log:
    data = np.log(data).astype(np.float32)
  # here we add one to num_features, because any apfp not founded in target_apfp_picked will be mapped
  # to the last column of the features matrix, though the last column will not be used ultimately.
  features = sparse.csr_matrix((data, indices, indptr), shape=(len(fps_list), num_features + 1))
  return features

class DatasetBase(object):
  def __init__(self, target):
    # read count and the apfps that were picked out
    counts = np.genfromtxt(mask_dir + "/%s_apfp.count" % target, delimiter="\t", dtype=int)
    self.target_apfp_picked = counts[counts[:, 1] > 10][:, 0]
    self.target_apfp_picked.sort()
    self.num_features = len(self.target_apfp_picked)
    # columns and sparse features
    # here we use a defaultdict, where any apfp not founded in target_apfp_picked will be mapped
    # to the last column of the features matrix, though the last column will not be used ultimately.
    self.target_columns_dict = defaultdict(lambda : self.num_features)
    for i, apfp in enumerate(self.target_apfp_picked):
      self.target_columns_dict[apfp] = i

  def batch_generator_base(self, size, batch_size):
    begin = 0
    end = 0
    while True:
      begin = end
      if begin >= size:
        raise StopIteration()
      end += batch_size
      if end > size:
        end = size
      yield begin, end

class DatasetTarget(DatasetBase):
  def __init__(self, target, year_split=2014):
    DatasetBase.__init__(self, target)
    # read chembl id and apfp
    self.chembl_id = []
    self.chembl_apfp = {}
    f = open(fp_dir + "/chembl.apfp", "r")
    for line in f:
      id_, fps_str = line.split("\t")
      id_ = id_.strip()
      fps_str = fps_str.strip()
      self.chembl_id.append(id_)
      self.chembl_apfp[id_] = fps_str
    f.close()
    # read top 79 targets' label data, and get the specific target's label data
    clf_label_79 = pd.read_csv(structure_dir + "/chembl_top79.label", usecols=[0, 2, 3, 4], delimiter="\t")
    self.target_clf_label = clf_label_79[clf_label_79["TARGET_CHEMBLID"] == target]
    # remove compounds whose apfp cannot be caculated
    m = self.target_clf_label["CMPD_CHEMBLID"].isin(self.chembl_id)
    self.target_clf_label = self.target_clf_label[m.values]
    # time split
    time_mask = self.target_clf_label["YEAR"] > year_split
    time_split_train = self.target_clf_label[~time_mask]
    time_split_test = self.target_clf_label[time_mask]
    # ids
    self.target_ids_train = time_split_train["CMPD_CHEMBLID"].values
    self.target_ids_test = time_split_test["CMPD_CHEMBLID"].values
    # features
    self.target_features_train = sparse_features([self.chembl_apfp[k] for k in self.target_ids_train], self.target_columns_dict, self.num_features)[:, :-1]
    self.target_features_test = sparse_features([self.chembl_apfp[k] for k in self.target_ids_test], self.target_columns_dict, self.num_features)[:, :-1]
    # labels
    self.target_labels_train = (time_split_train["CLF_LABEL"] > 0).astype(int).values
    self.target_labels_test = (time_split_test["CLF_LABEL"] > 0).astype(int).values

class DatasetCNS(DatasetTarget):
  def __init__(self, target, year_split=2014):
    DatasetTarget.__init__(self, target, year_split=year_split)
    # read mask
    self.cns_mask = pd.Series.from_csv(mask_dir + "/%s_cns.mask" % target, header=None, sep="\t")
    # features
    self.cns_features = sparse_features([self.chembl_apfp[k] for k in self.chembl_id], self.target_columns_dict, self.num_features)[:, :-1]
    #
    m = self.cns_mask.index.isin(self.target_ids_test)
    self.cns_features_train = self.cns_features[~m]
    self.cns_mask_train = self.cns_mask[~m]

  def batch_generator_cns(self, batch_size):
    for begin, end in self.batch_generator_base(self.cns_features.shape[0], batch_size):
      ids = self.chembl_id[begin: end]
      features = self.cns_features[begin: end].toarray()
      mask = self.cns_mask[begin: end].values
      yield ids, features, mask

class DatasetPNS(DatasetBase):
  def __init__(self, target):
    DatasetBase.__init__(self, target)
    # read pns apfp
    self.pns_id = []
    self.pns_apfp = {}
    f = open(fp_dir + "/pubchem_neg_sample.apfp", "r")
    for line in f:
      id_, fps_str = line.split("\t")
      id_ = id_.strip()
      fps_str = fps_str.strip()
      self.pns_id.append(id_)
      self.pns_apfp[id_] = fps_str
    f.close()
    # read mask
    self.pns_mask = pd.Series.from_csv(mask_dir + "/%s_pns.mask" % target, header=None, sep="\t")
    # features
    self.pns_features = sparse_features([self.pns_apfp[k] for k in self.pns_id], self.target_columns_dict, self.num_features)[:, :-1]

  def batch_generator_pns(self, batch_size):
    for begin, end in self.batch_generator_base(self.pns_features.shape[0], batch_size):
      ids = self.pns_id[begin: end]
      features = self.pns_features[begin: end].toarray()
      mask = self.pns_mask[begin: end].values
      yield ids, features, mask

class Dataset(DatasetCNS, DatasetPNS):
  """Base dataset class for chembl inhibitors
  """
  def __init__(self, target, one_hot=True, is_shuffle_train=True,  train_pos_multiply=0):
    """Constructor, create a dataset container.
    Args:
      target: <type 'str'> the chemblid of the target, e.g. "CHEMBL203".
      one_hot: <type 'bool'> flag whether create one_hot label, default is True.
      is_shuffle: <type 'bool'> flag whether shuffle samples when the dataset created.
      year_split: <type 'int'> time split year,
        if a molecule's year > year_split, it will be split into test data,
        otherwise, if a molecule's year <= year_split, it will be split into train data.
    Return:
      None
    """
    DatasetCNS.__init__(self, target, year_split=2014)
    DatasetPNS.__init__(self, target)
    # cns train pos
    self.cns_features_train_pos = self.cns_features_train[self.cns_mask_train.values]
    self.cns_mask_train_pos = self.cns_mask_train[self.cns_mask_train.values]
    # train, if train_pos_multiply > 0, cns_train_pos will be extra added for train_pos_multiply times .
    tf_list = [self.cns_features_train, self.pns_features]
    tl_list = [self.cns_mask_train, self.pns_mask]
    for _ in range(train_pos_multiply):
      tf_list.append(self.cns_features_train_pos)
      tl_list.append(self.cns_mask_train_pos)
    self.train_features = sparse.vstack(tf_list)
    self.train_labels = np.hstack(tl_list).astype(int)
    # test
    self.test_features = self.target_features_test
    self.test_labels = self.target_labels_test
    # one_hot
    if one_hot:
      self.train_labels_one_hot = dense_to_one_hot(self.train_labels)
      self.test_labels_one_hot = dense_to_one_hot(self.test_labels)
    # batch related
    self.train_size = self.train_features.shape[0] # (954049, 9412)
    self.train_perm = np.arange(self.train_size)
    if is_shuffle_train:
      np.random.shuffle(self.train_perm)
    self.train_begin = 0
    self.train_end = 0

  def generate_perm_for_train_batch(self, batch_size):
    """Create the permutation for a batch of train samples
    Args:
      batch_size: <type 'int'> the number of samples in the batch
    Return:
      perm: <type 'numpy.ndarray'> the permutation of samples which form a batch
    """
    self.train_begin = self.train_end
    self.train_end += batch_size
    if self.train_end > self.train_size:
      np.random.shuffle(self.train_perm)
      self.train_begin = 0
      self.train_end = batch_size
    perm = self.train_perm[self.train_begin: self.train_end]
    return perm

  def generate_train_batch(self, batch_size):
    perm = self.generate_perm_for_train_batch(batch_size)
    return self.train_features[perm].toarray().astype(np.float32), self.train_labels_one_hot[perm]

  def reset_begin_end(self):
    self.train_begin = 0
    self.train_end = 0

  def generate_train_batch_once(self, batch_size):
    self.train_begin = self.train_end
    self.train_end += batch_size
    if self.train_end > self.train_size:
      self.train_end = self.train_size
    perm = self.train_perm[self.train_begin: self.train_end]
    return self.train_features[perm].toarray().astype(np.float32), self.train_labels_one_hot[perm]

# dataset for virtual screening(vs)
class DatasetVS(DatasetBase):
  def __init__(self, target):
    DatasetBase.__init__(self, target)

  def reset(self, fp_fn):
    # read chembl id and apfp
    self.pubchem_id = []
    self.pubchem_apfp = {}
    f = open(fp_fn, "r")
    for line in f:
      id_, fps_str = line.split("\t")
      id_ = id_.strip()
      fps_str = fps_str.strip()
      self.pubchem_id.append(id_)
      self.pubchem_apfp[id_] = fps_str
    f.close()
    # generate features
    self.features = sparse_features([self.pubchem_apfp[k] for k in self.pubchem_id], self.target_columns_dict, self.num_features)[:, :-1]
    self.features_dense = self.features.toarray()

class DatasetChemDiv(DatasetBase):
  def __init__(self, target):
    DatasetBase.__init__(self, target)
    # read ids and apfps
    ChemDiv_dir = "/raid/xiaotaw/ChemDiv"
    fn_list = ["DC01_350000.apfp", "DC02_350000.apfp", "DC03_222773.apfp", "DC_saltdata_not-available_124145.apfp", "IC_non-excl_82693.apfp", "NC_340320.apfp"]
    self.chemdiv_ids = []
    self.chemdiv_apfps = {}
    for fn in fn_list:
      f = open(ChemDiv_dir + "/" + fn, "r")
      for line in f:
        id_, fps_str = line.split("\t")
        id_ = id_.strip()
        fps_str = fps_str.strip()
        self.chemdiv_ids.append(id_)
        self.chemdiv_apfps[id_] = fps_str
      f.close()
    # batch related
    self.begin = 0
    self.end = 0
    self.size = len(self.chemdiv_ids)

  def generate_batch(self, batch_size):
    self.begin = self.end
    if self.begin >= self.size:
      raise StopIteration()
    self.end += batch_size
    if self.end > self.size:
      self.end = self.size
    ids = self.chemdiv_ids[self.begin: self.end]
    apfp_list = [self.chemdiv_apfps[k] for k in ids]
    features = sparse_features(apfp_list, self.target_columns_dict, self.num_features)[:, :-1]
    return ids, features

  def batch_generator_chemdiv(self, batch_size):
    for begin, end in self.batch_generator_base(self.size, batch_size):
      ids = self.chemdiv_ids[begin: end]
      apfp_list = [self.chemdiv_apfps[k] for k in ids]
      features = sparse_features(apfp_list, self.target_columns_dict, self.num_features)[:, :-1].toarray()
      yield ids, features

def compute_performance(label, prediction):
  """sensitivity(SEN), specificity(SPE), accuracy(ACC), matthews correlation coefficient(MCC)
  """
  assert label.shape[0] == prediction.shape[0], "label number should be equal to prediction number"
  N = label.shape[0]
  APP = sum(prediction)
  ATP = sum(label)
  TP = sum(prediction * label)
  FP = APP - TP
  FN = ATP - TP
  TN = N - TP - FP - FN
  SEN = float(TP) / (ATP) if ATP != 0 else np.nan
  SPE = float(TN) / (N - ATP)
  ACC = float(TP + TN) / N
  MCC = (TP * TN - FP * FN) / (np.sqrt(long(N - APP) * long(N - ATP) * APP * ATP)) if not (N - APP) * (N - ATP) * APP * ATP == 0 else 0.0
  return TP, TN, FP, FN, SEN, SPE, ACC, MCC

reason: import numpy,from scipy