def nearest_neighbors(query_pts, target_pts=None, metric='euclidean',
                      k=None, epsilon=None, return_dists=False,
                      precomputed=False):
  '''Find nearest neighbors of query points from a matrix of target points.

  Returns a list of indices of neighboring points, one list per query.
  If no target_pts are specified, distances are calculated within query_pts.
  When return_dists is True, returns two lists: (distances, indices)
  '''
  if k is None and epsilon is None:
    raise ValueError('Must provide `k` or `epsilon`.')

  # TODO: deprecate the precomputed kwarg
  precomputed = precomputed or (metric == 'precomputed')

  if precomputed and target_pts is not None:
    raise ValueError('`target_pts` cannot be used with precomputed distances')

  query_pts = np.array(query_pts)
  if len(query_pts.shape) == 1:
    query_pts = query_pts.reshape((1,-1))  # ensure that the query is a 1xD row

  if precomputed:
    dists = query_pts.copy()
  else:
    dists = pairwise_distances(query_pts, Y=target_pts, metric=metric)

  if epsilon is not None:
    if k is not None:
      # kNN filtering
      _, not_nn = _min_k_indices(dists, k, inv_ind=True)
      dists[np.arange(dists.shape[0]), not_nn.T] = np.inf
    # epsilon-ball
    is_close = dists <= epsilon
    if return_dists:
      nnis,nnds = [],[]
      for i,row in enumerate(is_close):
        nns = np.nonzero(row)[0]
        nnis.append(nns)
        nnds.append(dists[i,nns])
      return nnds, nnis
    return np.array([np.nonzero(row)[0] for row in is_close])

  # knn
  nns = _min_k_indices(dists,k)
  if return_dists:
    # index each row of dists by each row of nns
    row_inds = np.arange(len(nns))[:,np.newaxis]
    nn_dists = dists[row_inds, nns]
    return nn_dists, nns
  return nns