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