code: def approximate_gradient(decision_function, sample, num_evals,
                           delta, constraint, shape, clip_min, clip_max):
    
    # Generate random vectors.
    noise_shape = [num_evals] + list(shape)
    if constraint == 'l2':
      rv = np.random.randn(*noise_shape)
    elif constraint == 'linf':
      rv = np.random.uniform(low=-1, high=1, size=noise_shape)
  
    axis = tuple(range(1, 1 + len(shape)))
    rv = rv / np.sqrt(np.sum(rv ** 2, axis=axis, keepdims=True))
    perturbed = sample + delta * rv
    perturbed = clip_image(perturbed, clip_min, clip_max)
    rv = (perturbed - sample) / delta
  
    # query the model.
    decisions = decision_function(perturbed)
    decision_shape = [len(decisions)] + [1] * len(shape)
    fval = 2 * decisions.astype(np_dtype).reshape(decision_shape) - 1.0
  
    # Baseline subtraction (when fval differs)
    if np.mean(fval) == 1.0:  # label changes.
      gradf = np.mean(rv, axis=0)
    elif np.mean(fval) == -1.0:  # label not change.
      gradf = - np.mean(rv, axis=0)
    else:
      fval = fval - np.mean(fval)
      gradf = np.mean(fval * rv, axis=0)
  
    # Get the gradient direction.
    gradf = gradf / np.linalg.norm(gradf)
  
    return gradf
docs: Gradient direction estimation