docs: def grad_numerical(self, x, func, epsilon=None):
          """symmetric gradient"""
          eps = 1e-8 * (1 + abs(x)) if epsilon is None else epsilon
          grad = np.zeros(len(x))
          ei = np.zeros(len(x))  # float is 1.6 times faster than int
          for i in rglen(x):
              ei[i] = eps[i]
              grad[i] = (func(x + ei) - func(x - ei)) / (2*eps[i])
              ei[i] = 0
          return grad