<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 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