docs: def fun_residuals(params, parmodel, bounddict,
                    shrinking_factor, numresolution,
                    islitmin, islitmax, debugplot):
      """Function to be minimised.
  
      Parameters
      ----------
      params : :class:`~lmfit.parameter.Parameters`
          Parameters to be employed in the prediction of the distorted
          boundaries.
      parmodel : str
          Model to be assumed. Allowed values are 'longslit' and
          'multislit'.
      bounddict : JSON structure
          Structure employed to store bounddict information.
      shrinking_factor : float
          Fraction of the detected X range (specrtral) to be employed
          in the fit. This must be a number verifying
          0 < shrinking_factor <= 1. The resulting interval will be
          centered within the original one.
      numresolution : int
          Number of points in which the X-range interval is subdivided
          before computing the residuals.
      islitmin : int
          Minimum slitlet number.
      islitmax : int
          Maximum slitlet number.
      debugplot : int
          Debugging level for messages and plots. For details see
          'numina.array.display.pause_debugplot.py'.
  
      Returns
      -------
      global_residual : float
          Squared root of the averaged sum of squared residuals.
  
      """
  
      global FUNCTION_EVALUATIONS
      global_residual = 0.0
      nsummed = 0
  
      read_slitlets = list(bounddict['contents'].keys())
      # read_slitlets.sort()  # this is not really necessary
      for tmp_slitlet in read_slitlets:
          islitlet = int(tmp_slitlet[7:])
          if islitmin <= islitlet <= islitmax:
              read_dateobs = list(bounddict['contents'][tmp_slitlet].keys())
              # read_dateobs.sort()  # this is not really necessary
              for tmp_dateobs in read_dateobs:
                  tmp_dict = bounddict['contents'][tmp_slitlet][tmp_dateobs]
                  csu_bar_slit_center = tmp_dict['csu_bar_slit_center']
                  # expected boundaries using provided parameters
                  list_spectrails = expected_distorted_boundaries(
                          islitlet, csu_bar_slit_center,
                          [0, 1], params, parmodel,
                          numpts=numresolution, deg=5, debugplot=0
                      )
                  poly_lower_expected = list_spectrails[0].poly_funct
                  poly_upper_expected = list_spectrails[1].poly_funct
                  # measured lower boundary
                  poly_lower_measured = np.polynomial.Polynomial(
                      tmp_dict['boundary_coef_lower']
                  )
                  xmin_lower_bound = tmp_dict['boundary_xmin_lower']
                  xmax_lower_bound = tmp_dict['boundary_xmax_lower']
                  dx = (xmax_lower_bound - xmin_lower_bound) * \
                       (1 - shrinking_factor) / 2
                  xdum_lower = np.linspace(xmin_lower_bound + dx,
                                           xmax_lower_bound - dx,
                                           num=numresolution)
                  # distance between expected and measured polynomials
                  poly_diff = poly_lower_expected - poly_lower_measured
                  global_residual += np.sum(poly_diff(xdum_lower)**2)
                  nsummed += numresolution
                  # measured upper boundary
                  poly_upper_measured = np.polynomial.Polynomial(
                      tmp_dict['boundary_coef_upper']
                  )
                  xmin_upper_bound = tmp_dict['boundary_xmin_upper']
                  xmax_upper_bound = tmp_dict['boundary_xmax_upper']
                  dx = (xmax_lower_bound - xmin_lower_bound) * \
                       (1 - shrinking_factor) / 2
                  xdum_upper = np.linspace(xmin_upper_bound + dx,
                                           xmax_upper_bound - dx,
                                           num=numresolution)
                  # distance between expected and measured polynomials
                  poly_diff = poly_upper_expected - poly_upper_measured
                  global_residual += np.sum(poly_diff(xdum_upper)**2)
                  nsummed += numresolution
  
      if nsummed > 0:
          global_residual = np.sqrt(global_residual/nsummed)
      if debugplot >= 10:
          FUNCTION_EVALUATIONS += 1
          print('-' * 79)
          print('>>> Number of function evaluations:', FUNCTION_EVALUATIONS)
          print('>>> global residual...............:', global_residual)
          params.pretty_print()
      return global_residual