docs: def entropy(data=None, prob=None, method='nearest-neighbors', bins=None, errorVal=1e-5, units='bits'):
      '''
      given a probability distribution (prob) or an interable of symbols (data) compute and
      return its continuous entropy.
  
      inputs:
      ------
          data:       samples by dimensions ndarray
  
          prob:       iterable with probabilities
  
          method:     'nearest-neighbors', 'gaussian', or 'bin'
  
          bins:       either a list of num_bins, or a list of lists containing
                      the bin edges
  
          errorVal:   if prob is given, 'entropy' checks that the sum is about 1.
                      It raises an error if abs(sum(prob)-1) >= errorVal
  
          units:      either 'bits' or 'nats'
  
          Different Methods:
  
          'nearest-neighbors' computes the binless entropy (bits) of a random vector
          using average nearest neighbors distance (Kozachenko and Leonenko, 1987).
          For a review see Beirlant et al., 2001 or Chandler & Field, 2007.
  
          'gaussian' computes the binless entropy based on estimating the covariance
          matrix and assuming the data is normally distributed.
  
          'bin' discretizes the data and computes the discrete entropy.
  
      '''
  
      if prob is None and data is None:
          raise ValueError("%s.entropy requires either 'prob' or 'data' to be defined" % __name__)
  
      if prob is not None and data is not None:
          raise ValueError("%s.entropy requires only 'prob' or 'data to be given but not both" % __name__)
  
      if prob is not None and not isinstance(prob, np.ndarray):
          raise TypeError("'entropy' in '%s' needs 'prob' to be an ndarray" % __name__)
  
      if prob is not None and abs(prob.sum()-1) > errorVal:
          raise ValueError("parameter 'prob' in '%s.entropy' should sum to 1" % __name__)
  
      if data is not None:
          num_samples = data.shape[0]
          if len(data.shape) == 1:
              num_dimensions = 1
          else:
              num_dimensions = data.shape[1]
  
      if method == 'nearest-neighbors':
          from sklearn.neighbors import NearestNeighbors
          from scipy.special import gamma
  
          if data is None:
              raise ValueError('Nearest neighbors entropy requires original data')
  
          if len(data.shape) > 1:
              k = num_dimensions
          else:
              k = 1
  
          nbrs = NearestNeighbors(n_neighbors=2, algorithm='auto').fit(data)
          distances, indices = nbrs.kneighbors(data)
          rho = distances[:,1] # take nearest-neighbor distance (first column is always zero)
          Ak  = (k*np.pi**(float(k)/float(2)))/gamma(float(k)/float(2)+1)
  
          if units is 'bits':
              # 0.577215... is the Euler-Mascheroni constant (np.euler_gamma)
              return k*np.mean(np.log2(rho)) + np.log2(num_samples*Ak/k) + np.log2(np.exp(1))*np.euler_gamma
          elif units is 'nats':
              # 0.577215... is the Euler-Mascheroni constant (np.euler_gamma)
              return k*np.mean(np.log(rho)) + np.log(num_samples*Ak/k) + np.log(np.exp(1))*np.euler_gamma
          else:
              print('Units not recognized: {}'.format(units))
  
  
      elif method == 'gaussian':
          from numpy.linalg import det
  
          if data is None:
              raise ValueError('Nearest neighbors entropy requires original data')
  
          detCov = det(np.dot(data.transpose(), data)/num_samples)
          normalization = (2*np.pi*np.exp(1))**num_dimensions
  
          if detCov == 0:
              return -np.inf
          else:
              if units is 'bits':
                  return 0.5*np.log2(normalization*detCov)
              elif units is 'nats':
                  return 0.5*np.log(normalization*detCov)
              else:
                  print('Units not recognized: {}'.format(units))
  
      elif method == 'bin':
          if prob is None and bins is None:
              raise ValueError('Either prob or bins must be specified.')
  
          if data is not None:
              prob = symbols_to_prob(data, bins=bins)
  
          if units is 'bits':
              # compute the log2 of the probability and change any -inf by 0s
              logProb = np.log2(prob)
              logProb[logProb == -np.inf] = 0
          elif units is 'nats':
              # compute the log2 of the probability and change any -inf by 0s
              logProb = np.log(prob)
              logProb[logProb == -np.inf] = 0
          else:
              print('Units not recognized: {}'.format(units))
  
          # return sum of product of logProb and prob
          # (not using np.dot here because prob, logprob are nd arrays)
          return -float(np.sum(prob * logProb))