Check wavelength calibration of the provided spectrum.

    Parameters
    ----------
    sp : numpy array
        Wavelength calibrated spectrum.
    crpix1: float
        CRPIX1 keyword.
    crval1: float
        CRVAL1 keyword.
    cdelt1: float
        CDELT1 keyword.
    wv_master: numpy array
        Array with the detailed list of expected arc lines.
    coeff_ini : array like
        Coefficients initially employed to obtain the wavelength
        calibration of the provided spectrum. When this coefficients
        are provided, this function computes a refined version of
        them, incorporating the corrections derived from the fit to
        the residuals.
    naxis1_ini : int
        NAXIS1 in original spectrum employed to fit the initial
        wavelength calibration.
    min_nlines_to_refine : int
        Minimum number of identified lines necessary to perform the
        wavelength calibration refinement. If zero, no minimum number
        is required.
    interactive : bool
        If True, the function allows the user to modify the residuals
        fit.
    threshold : float
        Minimum signal in the peaks.
    nwinwidth_initial : int
        Width of the window where each peak must be initially found.
    nwinwidth_refined : int
        Width of the window where each peak must be refined.
    ntimes_match_wv : float
        Times CDELT1 to match measured and expected wavelengths.
    poldeg_residuals : int
        Polynomial degree for fit to residuals.
    times_sigma_reject : float or None
        Number of times the standard deviation to reject points
        iteratively. If None, the fit does not reject any point.
    use_r : bool
        If True, additional statistical analysis is performed using R.
    title : string
        Plot title.
    remove_null_borders : bool
        If True, remove leading and trailing zeros in spectrum.
    ylogscale : bool
        If True, the spectrum is displayed in logarithmic units. Note
        that this is only employed for display purposes. The line peaks
        are found in the original spectrum.
    geometry : tuple (4 integers) or None
        x, y, dx, dy values employed to set the window geometry.
    pdf : PdfFile object or None
        If not None, output is sent to PDF file.
    debugplot : int
        Debugging level for messages and plots. For details see
        'numina.array.display.pause_debugplot.py'.

    Returns
    -------
    coeff_refined : numpy array
        Refined version of the initial wavelength calibration
        coefficients. These coefficients are computed only when
        the input parameter 'coeff_ini' is not None.

def check_wlcalib_sp(sp, crpix1, crval1, cdelt1, wv_master,
                     coeff_ini=None, naxis1_ini=None,
                     min_nlines_to_refine=0,
                     interactive=False,
                     threshold=0,
                     nwinwidth_initial=7,
                     nwinwidth_refined=5,
                     ntimes_match_wv=2,
                     poldeg_residuals=1,
                     times_sigma_reject=5,
                     use_r=False,
                     title=None,
                     remove_null_borders=True,
                     ylogscale=False,
                     geometry=None,
                     pdf=None,
                     debugplot=0):
    """Check wavelength calibration of the provided spectrum.

    Parameters
    ----------
    sp : numpy array
        Wavelength calibrated spectrum.
    crpix1: float
        CRPIX1 keyword.
    crval1: float
        CRVAL1 keyword.
    cdelt1: float
        CDELT1 keyword.
    wv_master: numpy array
        Array with the detailed list of expected arc lines.
    coeff_ini : array like
        Coefficients initially employed to obtain the wavelength
        calibration of the provided spectrum. When this coefficients
        are provided, this function computes a refined version of
        them, incorporating the corrections derived from the fit to
        the residuals.
    naxis1_ini : int
        NAXIS1 in original spectrum employed to fit the initial
        wavelength calibration.
    min_nlines_to_refine : int
        Minimum number of identified lines necessary to perform the
        wavelength calibration refinement. If zero, no minimum number
        is required.
    interactive : bool
        If True, the function allows the user to modify the residuals
        fit.
    threshold : float
        Minimum signal in the peaks.
    nwinwidth_initial : int
        Width of the window where each peak must be initially found.
    nwinwidth_refined : int
        Width of the window where each peak must be refined.
    ntimes_match_wv : float
        Times CDELT1 to match measured and expected wavelengths.
    poldeg_residuals : int
        Polynomial degree for fit to residuals.
    times_sigma_reject : float or None
        Number of times the standard deviation to reject points
        iteratively. If None, the fit does not reject any point.
    use_r : bool
        If True, additional statistical analysis is performed using R.
    title : string
        Plot title.
    remove_null_borders : bool
        If True, remove leading and trailing zeros in spectrum.
    ylogscale : bool
        If True, the spectrum is displayed in logarithmic units. Note
        that this is only employed for display purposes. The line peaks
        are found in the original spectrum.
    geometry : tuple (4 integers) or None
        x, y, dx, dy values employed to set the window geometry.
    pdf : PdfFile object or None
        If not None, output is sent to PDF file.
    debugplot : int
        Debugging level for messages and plots. For details see
        'numina.array.display.pause_debugplot.py'.

    Returns
    -------
    coeff_refined : numpy array
        Refined version of the initial wavelength calibration
        coefficients. These coefficients are computed only when
        the input parameter 'coeff_ini' is not None.

    """

    # protections
    if type(sp) is not np.ndarray:
        raise ValueError("sp must be a numpy.ndarray")
    elif sp.ndim != 1:
        raise ValueError("sp.ndim is not 1")

    if coeff_ini is None and naxis1_ini is None:
        pass
    elif coeff_ini is not None and naxis1_ini is not None:
        pass
    else:
        raise ValueError("coeff_ini and naxis1_ini must be simultaneously "
                         "None of both different from None")

    # check that interactive use takes place when plotting
    if interactive:
        if abs(debugplot) % 10 == 0:
            raise ValueError("ERROR: interative use of this function is not "
                             "possible when debugplot=", debugplot)

   # interactive and pdf are incompatible
    if interactive:
        if pdf is not None:
            raise ValueError("ERROR: interactive use of this function is not "
                             "possible when pdf is not None")

    # display list of expected arc lines
    if abs(debugplot) in (21, 22):
        print('wv_master:', wv_master)

    # determine spectrum length
    naxis1 = sp.shape[0]

    # define default values in case no useful lines are identified
    fxpeaks = np.array([])
    ixpeaks_wv = np.array([])
    fxpeaks_wv = np.array([])
    wv_verified_all_peaks = np.array([])
    nlines_ok = 0
    xresid = np.array([], dtype=float)
    yresid = np.array([], dtype=float)
    reject = np.array([], dtype=bool)
    polyres = np.polynomial.Polynomial([0])
    poldeg_effective = 0
    ysummary = summary(np.array([]))
    local_ylogscale = ylogscale

    # find initial line peaks
    ixpeaks = find_peaks_spectrum(sp,
                                  nwinwidth=nwinwidth_initial,
                                  threshold=threshold)
    npeaks = len(ixpeaks)

    if npeaks > 0:
        # refine location of line peaks
        fxpeaks, sxpeaks = refine_peaks_spectrum(
            sp, ixpeaks,
            nwinwidth=nwinwidth_refined,
            method="gaussian"
        )
        ixpeaks_wv = fun_wv(ixpeaks + 1, crpix1, crval1, cdelt1)
        fxpeaks_wv = fun_wv(fxpeaks + 1, crpix1, crval1, cdelt1)

        # match peaks with expected arc lines
        delta_wv_max = ntimes_match_wv * cdelt1
        wv_verified_all_peaks = match_wv_arrays(
            wv_master,
            fxpeaks_wv,
            delta_wv_max=delta_wv_max
        )

    loop = True

    while loop:

        if npeaks > 0:

            lines_ok = np.where(wv_verified_all_peaks > 0)
            nlines_ok = len(lines_ok[0])

            # there are matched lines
            if nlines_ok > 0:
                # compute residuals
                xresid = fxpeaks_wv[lines_ok]
                yresid = wv_verified_all_peaks[lines_ok] - fxpeaks_wv[lines_ok]

                # determine effective polynomial degree
                if nlines_ok > poldeg_residuals:
                    poldeg_effective = poldeg_residuals
                else:
                    poldeg_effective = nlines_ok  - 1

                # fit polynomial to residuals
                polyres, yresres, reject = \
                    polfit_residuals_with_sigma_rejection(
                        x=xresid,
                        y=yresid,
                        deg=poldeg_effective,
                        times_sigma_reject=times_sigma_reject,
                        use_r=use_r,
                        debugplot=0
                    )
                ysummary = summary(yresres)
            else:
                polyres = np.polynomial.Polynomial([0.0])

        list_wv_found = [str(round(wv, 4))
                         for wv in wv_verified_all_peaks if wv != 0]
        list_wv_master = [str(round(wv, 4)) for wv in wv_master]
        set1 = set(list_wv_master)
        set2 = set(list_wv_found)
        missing_wv = list(set1.symmetric_difference(set2))
        missing_wv.sort()
        if abs(debugplot) >= 10:
            print('-' * 79)
            print(">>> Number of arc lines in master file:", len(wv_master))
        if abs(debugplot) in [21, 22]:
            print(">>> Unmatched lines...................:", missing_wv)
        elif abs(debugplot) >= 10:
            print(">>> Number of unmatched lines.........:", len(missing_wv))
        if abs(debugplot) >= 10:
            print(">>> Number of line peaks found........:", npeaks)
            print(">>> Number of identified lines........:", nlines_ok)
            print(">>> Number of unmatched lines.........:", len(missing_wv))
            print(">>> Polynomial degree in residuals fit:", poldeg_effective)
            print(">>> Polynomial fit to residuals.......:\n", polyres)

        # display results
        if (abs(debugplot) % 10 != 0) or (pdf is not None):
            from numina.array.display.matplotlib_qt import plt
            if pdf is not None:
                fig = plt.figure(figsize=(11.69, 8.27), dpi=100)
            else:
                fig = plt.figure()
            set_window_geometry(geometry)

            # residuals
            ax2 = fig.add_subplot(2, 1, 1)
            if nlines_ok > 0:
                ymin = min(yresid)
                ymax = max(yresid)
                dy = ymax - ymin
                if dy > 0:
                    ymin -= dy/20
                    ymax += dy/20
                else:
                    ymin -= 0.5
                    ymax += 0.5
            else:
                ymin = -1.0
                ymax = 1.0
            ax2.set_ylim(ymin, ymax)
            if nlines_ok > 0:
                ax2.plot(xresid, yresid, 'o')
                ax2.plot(xresid[reject], yresid[reject], 'o', color='tab:gray')
            ax2.set_ylabel('Offset ' + r'($\AA$)')
            ax2.yaxis.label.set_size(10)
            if title is not None:
                ax2.set_title(title, **{'size': 12})
            xwv = fun_wv(np.arange(naxis1) + 1.0, crpix1, crval1, cdelt1)
            ax2.plot(xwv, polyres(xwv), '-')
            ax2.text(1, 0, 'CDELT1 (' + r'$\AA$' + '/pixel)=' + str(cdelt1),
                     horizontalalignment='right',
                     verticalalignment='bottom',
                     transform=ax2.transAxes)
            ax2.text(0, 0, 'Wavelength ' + r'($\AA$) --->',
                     horizontalalignment='left',
                     verticalalignment='bottom',
                     transform=ax2.transAxes)
            ax2.text(0, 1, 'median=' +
                     str(round(ysummary['median'], 4)) + r' $\AA$',
                     horizontalalignment='left',
                     verticalalignment='top',
                     transform=ax2.transAxes)
            ax2.text(0.5, 1, 'npoints (total / used / removed)',
                     horizontalalignment='center',
                     verticalalignment='top',
                     transform=ax2.transAxes)
            ax2.text(0.5, 0.92,
                     str(ysummary['npoints']) + ' / ' +
                     str(ysummary['npoints'] - sum(reject)) + ' / ' +
                     str(sum(reject)),
                     horizontalalignment='center',
                     verticalalignment='top',
                     transform=ax2.transAxes)
            ax2.text(1, 1, 'robust_std=' +
                     str(round(ysummary['robust_std'], 4)) + r' $\AA$',
                     horizontalalignment='right',
                     verticalalignment='top',
                     transform=ax2.transAxes)

            # median spectrum and peaks
            # remove leading and trailing zeros in spectrum when requested
            if remove_null_borders:
                nonzero = np.nonzero(sp)[0]
                j1 = nonzero[0]
                j2 = nonzero[-1]
                xmin = xwv[j1]
                xmax = xwv[j2]
            else:
                xmin = min(xwv)
                xmax = max(xwv)
            dx = xmax - xmin
            if dx > 0:
                xmin -= dx / 80
                xmax += dx / 80
            else:
                xmin -= 0.5
                xmax += 0.5

            if local_ylogscale:
                spectrum = sp - sp.min() + 1.0
                spectrum = np.log10(spectrum)
                ymin = spectrum[ixpeaks].min()
            else:
                spectrum = sp.copy()
                ymin = min(spectrum)
            ymax = max(spectrum)
            dy = ymax - ymin
            if dy > 0:
                ymin -= dy/20
                ymax += dy/20
            else:
                ymin -= 0.5
                ymax += 0.5
            ax1 = fig.add_subplot(2, 1, 2, sharex=ax2)
            ax1.set_xlim(xmin, xmax)
            ax1.set_ylim(ymin, ymax)
            ax1.plot(xwv, spectrum)
            if npeaks > 0:
                ax1.plot(ixpeaks_wv, spectrum[ixpeaks], 'o',
                         fillstyle='none', label="initial location")
                ax1.plot(fxpeaks_wv, spectrum[ixpeaks], 'o',
                         fillstyle='none', label="refined location")
                lok = wv_verified_all_peaks > 0
                ax1.plot(fxpeaks_wv[lok], spectrum[ixpeaks][lok], 'go',
                         label="valid line")
            if local_ylogscale:
                ax1.set_ylabel('~ log10(number of counts)')
            else:
                ax1.set_ylabel('number of counts')
            ax1.yaxis.label.set_size(10)
            ax1.xaxis.tick_top()
            ax1.xaxis.set_label_position('top')
            for i in range(len(ixpeaks)):
                # identified lines
                if wv_verified_all_peaks[i] > 0:
                    ax1.text(fxpeaks_wv[i], spectrum[ixpeaks[i]],
                             str(wv_verified_all_peaks[i]) +
                             '(' + str(i + 1) + ')',
                             fontsize=8,
                             horizontalalignment='center')
                else:
                    ax1.text(fxpeaks_wv[i], spectrum[ixpeaks[i]],
                             '(' + str(i + 1) + ')',
                             fontsize=8,
                             horizontalalignment='center')
                # estimated wavelength from initial calibration
                if npeaks > 0:
                    estimated_wv = fun_wv(fxpeaks[i] + 1,
                                          crpix1, crval1, cdelt1)
                    estimated_wv = str(round(estimated_wv, 4))
                    ax1.text(fxpeaks_wv[i], ymin,  # spmedian[ixpeaks[i]],
                             estimated_wv, fontsize=8, color='grey',
                             rotation='vertical',
                             horizontalalignment='center',
                             verticalalignment='top')
            if len(missing_wv) > 0:
                tmp = [float(wv) for wv in missing_wv]
                ax1.vlines(tmp, ymin=ymin, ymax=ymax,
                           colors='grey', linestyles='dotted',
                           label='missing lines')
            ax1.legend()
            if pdf is not None:
                pdf.savefig()
            else:
                if debugplot in [-22, -12, 12, 22]:
                    pause_debugplot(
                        debugplot=debugplot,
                        optional_prompt='Zoom/Unzoom or ' +
                                        'press RETURN to continue...',
                        pltshow=True
                    )
                else:
                    pause_debugplot(debugplot=debugplot, pltshow=True)

            # display results and request next action
            if interactive:
                print('Recalibration menu')
                print('------------------')
                print('[d] (d)elete all the identified lines')
                print('[r] (r)estart from begining')
                print('[a] (a)utomatic line inclusion')
                print('[l] toggle (l)ogarithmic scale on/off')
                print('[p] modify (p)olynomial degree')
                print('[o] (o)utput data with identified line peaks')
                print('[x] e(x)it without additional changes')
                print('[#] from 1 to ' + str(len(ixpeaks)) +
                      ' --> modify line #')
                ioption = readi('Option', default='x',
                                minval=1, maxval=len(ixpeaks),
                                allowed_single_chars='adloprx')
                if ioption == 'd':
                    wv_verified_all_peaks = np.zeros(npeaks)
                elif ioption == 'r':
                    delta_wv_max = ntimes_match_wv * cdelt1
                    wv_verified_all_peaks = match_wv_arrays(
                        wv_master,
                        fxpeaks_wv,
                        delta_wv_max=delta_wv_max
                    )
                elif ioption == 'a':
                    fxpeaks_wv_corrected = np.zeros_like(fxpeaks_wv)
                    for i in range(npeaks):
                        fxpeaks_wv_corrected[i] = fxpeaks_wv[i] + \
                            polyres(fxpeaks_wv[i])
                    delta_wv_max = ntimes_match_wv * cdelt1
                    wv_verified_all_peaks = match_wv_arrays(
                        wv_master,
                        fxpeaks_wv_corrected,
                        delta_wv_max=delta_wv_max
                    )
                elif ioption == 'l':
                    if local_ylogscale:
                        local_ylogscale = False
                    else:
                        local_ylogscale = True
                elif ioption == 'p':
                    poldeg_residuals = readi('New polynomial degree',
                                             minval=0)
                elif ioption == 'o':
                    for i in range(len(ixpeaks)):
                        # identified lines
                        if wv_verified_all_peaks[i] > 0:
                            print(wv_verified_all_peaks[i],
                                  spectrum[ixpeaks[i]])
                elif ioption == 'x':
                    loop = False
                else:
                    print(wv_master)
                    expected_value = fxpeaks_wv[ioption - 1] + \
                                     polyres(fxpeaks_wv[ioption - 1])
                    print(">>> Current expected wavelength: ", expected_value)
                    delta_wv_max = ntimes_match_wv * cdelt1
                    close_value = match_wv_arrays(
                        wv_master,
                        np.array([expected_value]),
                        delta_wv_max=delta_wv_max)
                    newvalue = readf('New value (0 to delete line)',
                                     default=close_value[0])
                    wv_verified_all_peaks[ioption - 1] = newvalue
            else:
                loop = False
        else:
            loop = False

    # refined wavelength calibration coefficients
    if coeff_ini is not None:
        npoints_total = len(xresid)
        npoints_removed = sum(reject)
        npoints_used = npoints_total - npoints_removed
        if abs(debugplot) >= 10:
            print('>>> Npoints (total / used / removed)..:',
                npoints_total, npoints_used, npoints_removed)
        if npoints_used < min_nlines_to_refine:
            print('Warning: number of lines insuficient to refine '
                  'wavelength calibration!')
            copc = 'n'
        else:
            if interactive:
                copc = readc('Refine wavelength calibration coefficients: '
                             '(y)es, (n)o', default='y', valid='yn')
            else:
                copc = 'y'
        if copc == 'y':
            coeff_refined = update_poly_wlcalib(
                    coeff_ini=coeff_ini,
                    coeff_residuals=polyres.coef,
                    naxis1_ini=naxis1_ini,
                    debugplot=0
                )
        else:
            coeff_refined = np.array(coeff_ini)
    else:
        coeff_refined = None

    if abs(debugplot) % 10 != 0:
        if coeff_refined is not None:
            for idum, fdum in \
                    enumerate(zip(coeff_ini, coeff_refined)):
                print(">>> coef#" + str(idum) + ':  ', end='')
                print("%+.8E  -->  %+.8E" % (decimal.Decimal(fdum[0]),
                                         decimal.Decimal(fdum[1])))

    return coeff_refined