code: def compute_gradient(self, *args, **kwargs):
          
          _EPS = 1.254e-5
          vector = self.get_parameter_vector()
          value0 = self.get_value(*args, **kwargs)
          grad = np.empty([len(vector)] + list(value0.shape), dtype=np.float64)
          for i, v in enumerate(vector):
              vector[i] = v + _EPS
              self.set_parameter_vector(vector)
              value = self.get_value(*args, **kwargs)
              vector[i] = v
              self.set_parameter_vector(vector)
              grad[i] = (value - value0) / _EPS
          return grad
docs: Compute the "gradient" of the model for the current parameters
  
          The default implementation computes the gradients numerically using
          a first order forward scheme. For better performance, this method
          should be overloaded by subclasses. The output of this function
          should be an array where the first dimension is ``full_size``.