docs: def _get_numdiff_hessian_steps_func(objective_func, nmr_steps, step_ratio):
      """Get a function to compute the multiple step sizes for a single element of the Hessian."""
      return SimpleCLFunction.from_string('''
          /**
           * Compute one element of the Hessian for a number of steps.
           * 
           * This uses the initial steps in the data structure, indexed by the parameters to change (px, py).
           *
           * Args:
           *  data: the data container
           *  x_tmp: the array with the input parameters, needs to be writable, although it will return
           *         the same values.
           *  f_x_input: the objective function value at the original set of parameters  
           *  px: the index of the first parameter to perturbate
           *  py: the index of the second parameter to perturbate
           *  steps: storage location for the output steps
           *  initial_step: the initial steps, array of same length as x_temp
           */
          void _numdiff_hessian_steps(void* data, local mot_float_type* x_tmp, 
                                      mot_float_type f_x_input,
                                      uint px, uint py, 
                                      local double* steps, 
                                      global float* initial_step){
  
              double step_x;
              double step_y;
              double tmp;
              bool is_first_workitem = get_local_id(0) == 0;
  
              if(px == py){
                  for(uint step_ind = 0; step_ind < ''' + str(nmr_steps) + '''; step_ind++){
                      step_x = initial_step[px] / pown(''' + str(float(step_ratio)) + ''', step_ind);
  
                      tmp = (
                            _numdiff_hessian_eval_step_mono(data, x_tmp, px, 2 * step_x)
                          + _numdiff_hessian_eval_step_mono(data, x_tmp, px, -2 * step_x)
                          - 2 * f_x_input
                      ) / (4 * step_x * step_x);    
  
                      if(is_first_workitem){
                          steps[step_ind] = tmp;
                      }
                  }
              }
              else{
                  for(uint step_ind = 0; step_ind < ''' + str(nmr_steps) + '''; step_ind++){
                      step_x = initial_step[px] / pown(''' + str(float(step_ratio)) + ''', step_ind);
                      step_y = initial_step[py] / pown(''' + str(float(step_ratio)) + ''', step_ind);
  
                      tmp = (
                            _numdiff_hessian_eval_step_bi(data, x_tmp, px, step_x, py, step_y)
                          - _numdiff_hessian_eval_step_bi(data, x_tmp, px, step_x, py, -step_y)
                          - _numdiff_hessian_eval_step_bi(data, x_tmp, px, -step_x, py, step_y)
                          + _numdiff_hessian_eval_step_bi(data, x_tmp, px, -step_x, py, -step_y)
                      ) / (4 * step_x * step_y);
  
                      if(is_first_workitem){
                          steps[step_ind] = tmp;
                      }                       
                  }
              }
          }
      ''', dependencies=[SimpleCLFunction.from_string('''
          /**
           * Evaluate the model with a perturbation in one dimensions.
           *
           * Args:
           *  data: the data container
           *  x_tmp: the array with the input parameters, needs to be writable, although it will return
           *         the same values.
           *  perturb_dim0: the index (into the x_tmp parameters) of the parameter to perturbate
           *  perturb_0: the added perturbation of the index corresponding to ``perturb_dim_0``
           *
           * Returns:
           *  the function evaluated at the parameters plus their perturbation.
           */
          double _numdiff_hessian_eval_step_mono(
                  void* data, local mot_float_type* x_tmp, 
                  uint perturb_dim_0, mot_float_type perturb_0){
  
              mot_float_type old_0;
              double return_val;
  
              if(get_local_id(0) == 0){
                  old_0 = x_tmp[perturb_dim_0];
                  x_tmp[perturb_dim_0] += perturb_0;
              }
              barrier(CLK_LOCAL_MEM_FENCE);
  
              return_val = ''' + objective_func.get_cl_function_name() + '''(x_tmp, data);
              barrier(CLK_LOCAL_MEM_FENCE);
  
              if(get_local_id(0) == 0){
                  x_tmp[perturb_dim_0] = old_0;
              }
              barrier(CLK_LOCAL_MEM_FENCE);
  
              return return_val;
          }
      '''), SimpleCLFunction.from_string('''
          /**
           * Evaluate the model with a perturbation in two dimensions.
           *
           * Args:
           *  data: the data container
           *  x_tmp: the array with the input parameters, needs to be writable, although it will return
           *         the same values.
           *  perturb_dim_0: the index (into the x_tmp parameters) of the first parameter to perturbate
           *  perturb_0: the added perturbation of the index corresponding to ``perturb_dim_0``
           *  perturb_dim_1: the index (into the x_tmp parameters) of the second parameter to perturbate
           *  perturb_1: the added perturbation of the index corresponding to ``perturb_dim_1``
           *
           * Returns:
           *  the function evaluated at the parameters plus their perturbation.
           */
          double _numdiff_hessian_eval_step_bi(
                  void* data, local mot_float_type* x_tmp, 
                  uint perturb_dim_0, mot_float_type perturb_0,
                  uint perturb_dim_1, mot_float_type perturb_1){
  
              mot_float_type old_0;
              mot_float_type old_1;
              double return_val;
  
              if(get_local_id(0) == 0){
                  old_0 = x_tmp[perturb_dim_0];
                  old_1 = x_tmp[perturb_dim_1];
  
                  x_tmp[perturb_dim_0] += perturb_0;
                  x_tmp[perturb_dim_1] += perturb_1;
              }
              barrier(CLK_LOCAL_MEM_FENCE);
  
              return_val = ''' + objective_func.get_cl_function_name() + '''(x_tmp, data);
              barrier(CLK_LOCAL_MEM_FENCE);
  
              if(get_local_id(0) == 0){
                  x_tmp[perturb_dim_0] = old_0;
                  x_tmp[perturb_dim_1] = old_1;
              }
              barrier(CLK_LOCAL_MEM_FENCE);
  
              return return_val;
          }
      ''')])