docs: def _get_penalty_function(nmr_parameters, constraints_func=None):
      """Get a function to compute the penalty term for the boundary conditions.
  
      This is meant to be used in the evaluation function of the optimization routines.
  
      Args:
          nmr_parameters (int): the number of parameters in the model
          constraints_func (mot.optimize.base.ConstraintFunction): function to compute (inequality) constraints.
              Should hold a CL function with the signature:
  
              .. code-block:: c
  
                  void <func_name>(local const mot_float_type* const x,
                                   void* data,
                                   local mot_float_type* constraint_values);
  
              Where ``constraints_values`` is filled as:
  
              .. code-block:: c
  
                  constraint_values[i] = g_i(x)
  
              That is, for each constraint function :math:`g_i`, formulated as :math:`g_i(x) <= 0`, we should return
              the function value of :math:`g_i`.
  
      Returns:
          tuple: Struct and SimpleCLFunction, the required data for the penalty function and the penalty function itself.
      """
      dependencies = []
      data_requirements = {'scratch': LocalMemory('double', 1)}
      constraints_code = ''
  
      if constraints_func and constraints_func.get_nmr_constraints() > 0:
          nmr_constraints = constraints_func.get_nmr_constraints()
          dependencies.append(constraints_func)
          data_requirements['constraints'] = LocalMemory('mot_float_type', nmr_constraints)
  
          constraints_code = '''
              local mot_float_type* constraints = ((_mle_penalty_data*)scratch_data)->constraints;
              
              ''' + constraints_func.get_cl_function_name() + '''(x, data, constraints);
              
              for(int i = 0; i < ''' + str(nmr_constraints) + '''; i++){
                  *penalty_sum += pown(max((mot_float_type)0, constraints[i]), 2); 
              }
          '''
  
      data = Struct(data_requirements, '_mle_penalty_data')
      func = SimpleCLFunction.from_string('''
          double _mle_penalty(
                  local mot_float_type* x,
                  void* data, 
                  local mot_float_type* lower_bounds, 
                  local mot_float_type* upper_bounds,
                  float penalty_weight,
                  void* scratch_data){
              
              local double* penalty_sum = ((_mle_penalty_data*)scratch_data)->scratch;
              
              if(get_local_id(0) == 0){
                  *penalty_sum = 0;
                  
                  // boundary conditions
                  for(int i = 0; i < ''' + str(nmr_parameters) + '''; i++){
                      if(isfinite(upper_bounds[i])){
                          *penalty_sum += pown(max((mot_float_type)0, x[i] - upper_bounds[i]), 2);    
                      }
                      if(isfinite(lower_bounds[i])){
                          *penalty_sum += pown(max((mot_float_type)0, lower_bounds[i] - x[i]), 2);
                      }
                  }
              }
              barrier(CLK_LOCAL_MEM_FENCE);
              
              // constraints
              ''' + constraints_code + '''
              
              return penalty_weight * *penalty_sum;
          }
      ''', dependencies=dependencies)
      return data, func