docs: def estimate_hessian(objective_func, parameters,
                       lower_bounds=None, upper_bounds=None,
                       step_ratio=2, nmr_steps=5,
                       max_step_sizes=None,
                       data=None, cl_runtime_info=None):
      """Estimate and return the upper triangular elements of the Hessian of the given function at the given parameters.
  
      This calculates the Hessian using central difference (using a 2nd order Taylor expansion) with a Richardson
      extrapolation over the proposed sequence of steps. If enough steps are given, we apply a Wynn epsilon extrapolation
      on top of the Richardson extrapolated results. If more steps are left, we return the estimate with the lowest error,
      taking into account outliers using a median filter.
  
      The Hessian is evaluated at the steps:
  
      .. math::
          \quad  ((f(x + d_j e_j + d_k e_k) - f(x + d_j e_j - d_k e_k)) -
                  (f(x - d_j e_j + d_k e_k) - f(x - d_j e_j - d_k e_k)) /
                  (4 d_j d_k)
  
      where :math:`e_j` is a vector where element :math:`j` is one and the rest are zero
      and :math:`d_j` is a scalar spacing :math:`steps_j`.
  
      Steps are generated according to an exponentially diminishing ratio, defined as:
  
          steps = max_step * step_ratio**-i, i = 0,1,..,nmr_steps-1.
  
      Where the maximum step can be provided. For example, a maximum step of 2 with a step ratio of 2, computed for
      4 steps gives: [2.0, 1.0, 0.5, 0.25]. If lower and upper bounds are given, we use as maximum step size the largest
      step size that fits between the Hessian point and the boundaries.
  
      The steps define the order of the estimation, with 2 steps resulting in a O(h^2) estimate, 3 steps resulting in a
      O(h^4) estimate and 4 or more steps resulting in a O(h^6) derivative estimate.
  
      Args:
          objective_func (mot.lib.cl_function.CLFunction): The function we want to differentiate.
              A CL function with the signature:
  
              .. code-block:: c
  
                  double <func_name>(local const mot_float_type* const x, void* data);
  
              The objective function has the same signature as the minimization function in MOT. For the numerical
              hessian, the ``objective_list`` parameter is ignored.
  
          parameters (ndarray): The parameters at which to evaluate the gradient. A (d, p) matrix with d problems,
              and p parameters
          lower_bounds (list or None): a list of length (p,) for p parameters with the lower bounds.
              Each element of the list can be a scalar or a vector (of the same length as the number
              of problem instances). To disable bounds for this parameter use -np.inf.
          upper_bounds (list or None): a list of length (p,) for p parameters with the upper bounds.
              Each element of the list can be a scalar or a vector (of the same length as the number
              of problem instances). To disable bounds for this parameter use np.inf.
          step_ratio (float): the ratio at which the steps diminish.
          nmr_steps (int): the number of steps we will generate. We will calculate the derivative for each of these
              step sizes and extrapolate the best step size from among them. The minimum number of steps is 1.
          max_step_sizes (float or ndarray or None): the maximum step size, or the maximum step size per parameter.
              If None is given, we use 0.1 for all parameters. If a float is given, we use that for all parameters.
              If a list is given, it should be of the same length as the number of parameters.
          data (mot.lib.kernel_data.KernelData): the user provided data for the ``void* data`` pointer.
          cl_runtime_info (mot.configuration.CLRuntimeInfo): the runtime information
  
      Returns:
          ndarray: per problem instance a vector with the upper triangular elements of the Hessian matrix.
              This array can hold NaN's, for elements where the Hessian failed to approximate.
      """
      if len(parameters.shape) == 1:
          parameters = parameters[None, :]
  
      nmr_voxels = parameters.shape[0]
      nmr_params = parameters.shape[1]
      nmr_derivatives = nmr_params * (nmr_params + 1) // 2
  
      initial_step = _get_initial_step(parameters, lower_bounds, upper_bounds, max_step_sizes)
  
      kernel_data = {
          'parameters': Array(parameters, ctype='mot_float_type'),
          'initial_step': Array(initial_step, ctype='float'),
          'derivatives': Zeros((nmr_voxels, nmr_derivatives), 'double'),
          'errors': Zeros((nmr_voxels, nmr_derivatives), 'double'),
          'x_tmp': LocalMemory('mot_float_type', nmr_params),
          'data': data,
          'scratch': LocalMemory('double', nmr_steps + (nmr_steps - 1) + nmr_steps)
      }
  
      hessian_kernel = SimpleCLFunction.from_string('''
          void _numdiff_hessian(
                  global mot_float_type* parameters,
                  global float* initial_step,
                  global double* derivatives,
                  global double* errors,
                  local mot_float_type* x_tmp,
                  void* data,
                  local double* scratch){
  
              if(get_local_id(0) == 0){
                  for(uint i = 0; i < ''' + str(nmr_params) + '''; i++){
                      x_tmp[i] = parameters[i];
                  }
              }
              barrier(CLK_LOCAL_MEM_FENCE);
  
              double f_x_input = ''' + objective_func.get_cl_function_name() + '''(x_tmp, data);
  
              // upper triangle loop
              uint coord_ind = 0;
              for(int i = 0; i < ''' + str(nmr_params) + '''; i++){
                  for(int j = i; j < ''' + str(nmr_params) + '''; j++){
                      _numdiff_hessian_element(
                          data, x_tmp, f_x_input, i, j, initial_step, 
                          derivatives + coord_ind, errors + coord_ind, scratch);
  
                      coord_ind++;
                  }
              }
          }
      ''', dependencies=[objective_func,
                         _get_numdiff_hessian_element_func(objective_func, nmr_steps, step_ratio)])
  
      hessian_kernel.evaluate(kernel_data, nmr_voxels, use_local_reduction=True, cl_runtime_info=cl_runtime_info)
  
      return kernel_data['derivatives'].get_data()