docs: def _lm_numdiff_jacobian(eval_func, nmr_params, nmr_observations):
      """Get a numerical differentiated Jacobian function.
  
      This computes the Jacobian of the observations (function vector) with respect to the parameters.
  
      Args:
          eval_func (mot.lib.cl_function.CLFunction): the evaluation function
          nmr_params (int): the number of parameters
          nmr_observations (int): the number of observations (the length of the function vector).
  
      Returns:
          mot.lib.cl_function.CLFunction: CL function for numerically estimating the Jacobian.
      """
      return SimpleCLFunction.from_string(r'''
          /**
           * Compute the Jacobian for use in the Levenberg-Marquardt optimization routine.
           *
           * This should place the output in the ``fjac`` matrix.
           *
           * Parameters:
           *
           *   model_parameters: (nmr_params,) the current point around which we want to know the Jacobian
           *   data: the current modeling data, used by the objective function
           *   fvec: (nmr_observations,), the function values corresponding to the current model parameters
           *   fjac: (nmr_parameters, nmr_observations), the memory location for the Jacobian
           */
          void _lm_numdiff_jacobian(local mot_float_type* model_parameters,
                                    void* data,
                                    local mot_float_type* fvec,
                                    local mot_float_type* const fjac){
              
              const uint nmr_params = ''' + str(nmr_params) + ''';
              const uint nmr_observations = ''' + str(nmr_observations) + ''';
              
              local mot_float_type* lower_bounds = ((_lm_eval_func_data*)data)->lower_bounds;
              local mot_float_type* upper_bounds = ((_lm_eval_func_data*)data)->upper_bounds;
              local mot_float_type* jacobian_x_tmp = ((_lm_eval_func_data*)data)->jacobian_x_tmp;
              
              mot_float_type step_size = 30 * MOT_EPSILON;
              
              for (int i = 0; i < nmr_params; i++) {
                  if(model_parameters[i] + step_size > upper_bounds[i]){
                      _lm_numdiff_jacobian_backwards(model_parameters, i, step_size, data, fvec, 
                                                     fjac + i*nmr_observations, jacobian_x_tmp);
                  }
                  else if(model_parameters[i] - step_size < lower_bounds[i]){
                      _lm_numdiff_jacobian_forwards(model_parameters, i, step_size, data, fvec, 
                                                    fjac + i*nmr_observations, jacobian_x_tmp);
                  }
                  else{
                      _lm_numdiff_jacobian_centered(model_parameters, i, step_size, data, fvec,  
                                                    fjac + i*nmr_observations, jacobian_x_tmp);
                  }
              }
          }
      ''', dependencies=[eval_func, SimpleCLFunction.from_string('''
          void _lm_numdiff_jacobian_centered(
                  local mot_float_type* model_parameters,
                  uint px,
                  float step_size,
                  void* data,
                  local mot_float_type* fvec,
                  local mot_float_type* const fjac,
                  local mot_float_type* fjac_tmp){
              
              const uint nmr_observations = ''' + str(nmr_observations) + ''';
              
              uint observation_ind;
              mot_float_type temp;
              
              uint local_id = get_local_id(0);
              uint workgroup_size = get_local_size(0);
              
              // F(x + h)     
              if(get_local_id(0) == 0){
                  temp = model_parameters[px];
                  model_parameters[px] = temp + step_size;
              }
              barrier(CLK_LOCAL_MEM_FENCE);
              ''' + eval_func.get_cl_function_name() + '''(model_parameters, data, fjac);
              
              
              // F(x - h)
              if(get_local_id(0) == 0){
                  model_parameters[px] = temp - step_size;
              }
              barrier(CLK_LOCAL_MEM_FENCE);
              ''' + eval_func.get_cl_function_name() + '''(model_parameters, data, fjac_tmp);
              
              
              // combine
              for(int i = 0; i < (nmr_observations + workgroup_size - 1) / workgroup_size; i++){
                  observation_ind = i * workgroup_size + local_id;
                  if(observation_ind < nmr_observations){
                      fjac[observation_ind] = (fjac[observation_ind] - fjac_tmp[observation_ind]) / (2 * step_size);
                  }
              }
              
              // restore parameter vector
              if(get_local_id(0) == 0){
                  model_parameters[px] = temp;
              }
              barrier(CLK_LOCAL_MEM_FENCE);
          }
      '''), SimpleCLFunction.from_string('''
          void _lm_numdiff_jacobian_backwards(
                 local mot_float_type* model_parameters,
                 uint px,
                 float step_size,
                 void* data,
                 local mot_float_type* fvec,
                 local mot_float_type* const fjac,
                 local mot_float_type* fjac_tmp){
  
              const uint nmr_observations = ''' + str(nmr_observations) + ''';
              
              uint observation_ind;
              mot_float_type temp;
              
              uint local_id = get_local_id(0);
              uint workgroup_size = get_local_size(0);
              
              // F(x - h)     
              if(get_local_id(0) == 0){
                 temp = model_parameters[px];
                 model_parameters[px] = temp - step_size;
              }
              barrier(CLK_LOCAL_MEM_FENCE);
              ''' + eval_func.get_cl_function_name() + '''(model_parameters, data, fjac);
                          
              
              // F(x - 2*h)
              if(get_local_id(0) == 0){
                 model_parameters[px] = temp - 2 * step_size;
              }
              barrier(CLK_LOCAL_MEM_FENCE);
              ''' + eval_func.get_cl_function_name() + '''(model_parameters, data, fjac_tmp);
              
              // combine
              for(int i = 0; i < (nmr_observations + workgroup_size - 1) / workgroup_size; i++){
                 observation_ind = i * workgroup_size + local_id;
                 if(observation_ind < nmr_observations){
                     fjac[observation_ind] = (  3 * fvec[observation_ind] 
                                              - 4 * fjac[observation_ind] 
                                              + fjac_tmp[observation_ind]) / (2 * step_size);
                 }
              }
              
              // restore parameter vector
              if(get_local_id(0) == 0){
                 model_parameters[px] = temp;
              }
              barrier(CLK_LOCAL_MEM_FENCE);
          }
      '''), SimpleCLFunction.from_string('''
          void _lm_numdiff_jacobian_forwards(
                 local mot_float_type* model_parameters,
                 uint px,
                 float step_size,
                 void* data,
                 local mot_float_type* fvec,
                 local mot_float_type* const fjac,
                 local mot_float_type* fjac_tmp){
  
              const uint nmr_observations = ''' + str(nmr_observations) + ''';
              
              uint observation_ind;
              mot_float_type temp;
              
              uint local_id = get_local_id(0);
              uint workgroup_size = get_local_size(0);
              
              // F(x + h)     
              if(get_local_id(0) == 0){
                 temp = model_parameters[px];
                 model_parameters[px] = temp + step_size;
              }
              barrier(CLK_LOCAL_MEM_FENCE);
              ''' + eval_func.get_cl_function_name() + '''(model_parameters, data, fjac);
              
            
              // F(x + 2*h)
              if(get_local_id(0) == 0){
                 model_parameters[px] = temp + 2 * step_size;
              }
              barrier(CLK_LOCAL_MEM_FENCE);
              ''' + eval_func.get_cl_function_name() + '''(model_parameters, data, fjac);
              
              // combine
              for(int i = 0; i < (nmr_observations + workgroup_size - 1) / workgroup_size; i++){
                 observation_ind = i * workgroup_size + local_id;
                 if(observation_ind < nmr_observations){
                     fjac[observation_ind] = (- 3 * fvec[observation_ind] 
                                              + 4 * fjac[observation_ind] 
                                              - fjac_tmp[observation_ind]) / (2 * step_size);
                 }
              }
              
              // restore parameter vector
              if(get_local_id(0) == 0){
                 model_parameters[px] = temp;
              }
              barrier(CLK_LOCAL_MEM_FENCE);
          }
      ''')])