Refine wavelength calibration using an initial polynomial.

    Parameters
    ----------
    sp : numpy array
        1D array of length NAXIS1 containing the input spectrum.
    poly_initial : Polynomial instance
        Initial wavelength calibration polynomial, providing the
        wavelength as a function of pixel number (running from 1 to
        NAXIS1).
    wv_master : numpy array
        Array containing the master list of arc line wavelengths.
    poldeg : int
        Polynomial degree of refined wavelength calibration. Note
        that this degree can be different from the polynomial degree
        of poly_initial.
    nrepeat : int
        Number of times lines are iteratively included in the initial
        fit.
    ntimes_match_wv : int
        Number of pixels around each line peak where the expected
        wavelength must match the tabulated wavelength in the master
        list.
    nwinwidth_initial : int
        Initial window width to search for line peaks in spectrum.
    nwinwidth_refined : int
        Window width to refine line peak location.
    times_sigma_reject : float
        Times sigma to reject points in the fit.
    interactive : bool
        If True, the function allows the user to modify the fit
        interactively.
    threshold : float
        Minimum signal in the peaks.
    plottitle : string or None
        Plot title.
    decimal_places : int
        Number of decimal places to be employed when displaying relevant
        fitted parameters.
    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
    -------
    poly_refined : Polynomial instance
        Refined wavelength calibration polynomial.
    yres_summary : dictionary
        Statistical summary of the residuals.

def refine_arccalibration(sp, poly_initial, wv_master, poldeg,
                          nrepeat=3,
                          ntimes_match_wv=2,
                          nwinwidth_initial=7,
                          nwinwidth_refined=5,
                          times_sigma_reject=5,
                          interactive=False,
                          threshold=0,
                          plottitle=None,
                          decimal_places=4,
                          ylogscale=False,
                          geometry=None,
                          pdf=None,
                          debugplot=0):
    """Refine wavelength calibration using an initial polynomial.

    Parameters
    ----------
    sp : numpy array
        1D array of length NAXIS1 containing the input spectrum.
    poly_initial : Polynomial instance
        Initial wavelength calibration polynomial, providing the
        wavelength as a function of pixel number (running from 1 to
        NAXIS1).
    wv_master : numpy array
        Array containing the master list of arc line wavelengths.
    poldeg : int
        Polynomial degree of refined wavelength calibration. Note
        that this degree can be different from the polynomial degree
        of poly_initial.
    nrepeat : int
        Number of times lines are iteratively included in the initial
        fit.
    ntimes_match_wv : int
        Number of pixels around each line peak where the expected
        wavelength must match the tabulated wavelength in the master
        list.
    nwinwidth_initial : int
        Initial window width to search for line peaks in spectrum.
    nwinwidth_refined : int
        Window width to refine line peak location.
    times_sigma_reject : float
        Times sigma to reject points in the fit.
    interactive : bool
        If True, the function allows the user to modify the fit
        interactively.
    threshold : float
        Minimum signal in the peaks.
    plottitle : string or None
        Plot title.
    decimal_places : int
        Number of decimal places to be employed when displaying relevant
        fitted parameters.
    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
    -------
    poly_refined : Polynomial instance
        Refined wavelength calibration polynomial.
    yres_summary : dictionary
        Statistical summary of the residuals.

    """

    # check that nrepeat is larger than zero
    if nrepeat <= 0:
        raise ValueError("Unexpected nrepeat=", str(nrepeat))

    # check that the requested polynomial degree is equal or larger than
    # the degree of the initial polynomial
    if poldeg < len(poly_initial.coef) - 1:
        raise ValueError("Polynomial degree of refined polynomial must be "
                         "equal or larger than that of the initial polynomial")

    # check that interactive use takes place when plotting
    if interactive:
        if abs(debugplot) % 10 != 0:
            local_debugplot = debugplot
        else:
            local_debugplot = -12
    else:
        local_debugplot = debugplot

    local_ylogscale = ylogscale

    # latest limits
    global xmin_previous
    global xmax_previous
    global ymin_previous
    global ymax_previous
    xmin_previous = None
    xmax_previous = None
    ymin_previous = None
    ymax_previous = None

    # spectrum length
    naxis1 = sp.shape[0]

    # define default values in case no useful lines are identified
    fxpeaks = np.array([])
    poly_refined = np.polynomial.Polynomial([0.0])
    poldeg_effective = len(poly_refined.coef) - 1
    yres_summary = summary(np.array([]))

    # compute linear values from initial polynomial
    crmin1_linear = poly_initial(1)
    crmax1_linear = poly_initial(naxis1)
    cdelt1_linear = (crmax1_linear - crmin1_linear) / (naxis1 - 1)

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

    if npeaks > 0:

        # refine line peak locations
        fxpeaks, sxpeaks = refine_peaks_spectrum(
            sp, ixpeaks,
            nwinwidth=nwinwidth_refined,
            method="gaussian"
        )

        # expected wavelength of all identified peaks
        wv_expected_all_peaks = poly_initial(fxpeaks + 1.0)

        # assign individual arc lines from master list to spectrum
        # line peaks when the expected wavelength is within the maximum
        # allowed range (+/- ntimes_match_wv * CDELT1 around the peak)
        delta_wv_max = ntimes_match_wv * cdelt1_linear
        wv_verified_all_peaks = match_wv_arrays(
            wv_master,
            wv_expected_all_peaks,
            delta_wv_max=delta_wv_max
        )

    loop = True

    nrepeat_eff = nrepeat
    while loop:

        nlines_ok = 0
        xdum = np.array([])
        ydum = np.array([])
        reject = np.array([])

        if npeaks > 0:

            for irepeat in range(nrepeat_eff):
                # fit with sigma rejection
                lines_ok = np.where(wv_verified_all_peaks > 0)
                nlines_ok = len(lines_ok[0])

                # there are matched lines
                if nlines_ok > 0:
                    # select points to be fitted
                    xdum = (fxpeaks + 1.0)[lines_ok]
                    ydum = wv_verified_all_peaks[lines_ok]

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

                    # fit polynomial
                    poly_refined, yres, reject = \
                        polfit_residuals_with_sigma_rejection(
                            x=xdum,
                            y=ydum,
                            deg=poldeg_effective,
                            times_sigma_reject=times_sigma_reject,
                            debugplot=0
                        )

                    # effective number of points
                    yres_summary = summary(yres[np.logical_not(reject)])

                else:
                    poly_refined = np.polynomial.Polynomial([0.0])
                    yres_summary = summary(np.array([]))

                if irepeat < nrepeat_eff - 1:
                    delta_wv_max = ntimes_match_wv * cdelt1_linear
                    wv_verified_all_peaks = match_wv_arrays(
                        wv_master,
                        poly_refined(fxpeaks + 1.0),
                        delta_wv_max=delta_wv_max
                    )

        # update poldeg_effective
        poldeg_effective = len(poly_refined.coef) - 1

        # update linear values
        crpix1_linear = 1.0
        crval1_linear = poly_refined(crpix1_linear)
        crmin1_linear = poly_refined(1)
        crmax1_linear = poly_refined(naxis1)
        cdelt1_linear = (crmax1_linear - crmin1_linear) / (naxis1 - 1)

        if abs(local_debugplot) >= 10:
            print(79 * '=')
            print(">>> poldeg (requested, effective)..:",
                  poldeg, poldeg_effective)
            print(">>> Fitted coefficients............:\n", poly_refined.coef)
            print(">>> NAXIS1.........................:", naxis1)
            print(">>> CRVAL1 linear scale............:", crval1_linear)
            print(">>> CDELT1 linear scale............:", cdelt1_linear)
            print(79 * '.')
            print(">>> Number of peaks................:", npeaks)
            print(">>> nlines identified (total, used): ",
                  nlines_ok, yres_summary['npoints'])
            print(">>> robust_std.....................:",
                  yres_summary['robust_std'])
            print(79 * '-')

        if (abs(local_debugplot) % 10 != 0) or (pdf is not None):
            from numina.array.display.matplotlib_qt import plt

            def handle_close(evt):
                global xmin_previous
                global xmax_previous
                global ymin_previous
                global ymax_previous
                xmin_previous, xmax_previous = ax2.get_xlim()
                ymin_previous, ymax_previous = ax2.get_ylim()

            if pdf is not None:
                fig = plt.figure(figsize=(11.69, 8.27), dpi=100)
            else:
                fig = plt.figure()
                fig.canvas.mpl_connect('close_event', handle_close)
            set_window_geometry(geometry)

            grid = plt.GridSpec(2, 1)
            grid.update(left=0.10, right=0.98,
                        bottom=0.10, top=0.90, hspace=0.01)

            # differences between linear fit and fitted polynomial
            # polynomial fit
            xpol = np.arange(1, naxis1 + 1)
            ylinear = crval1_linear + (xpol - crpix1_linear) * cdelt1_linear
            ypol = poly_refined(xpol) - ylinear
            # identified lines
            yp = ydum - (crval1_linear + (xdum - crpix1_linear) *
                         cdelt1_linear)

            # upper plot
            ax1 = fig.add_subplot(grid[0, 0])
            if nlines_ok > 0:
                ymin = min(ypol)
                ymax = max(ypol)
                dy = ymax - ymin
                if dy > 0:
                    ymin -= dy/50
                    ymax += dy/50
                else:
                    ymin -= 0.5
                    ymax += 0.5
            else:
                ymin = -1.0
                ymax = 1.0
            ax1.set_xlim(1 - 0.02 * naxis1, naxis1 * 1.02)
            ax1.set_ylim(ymin, ymax)
            if nlines_ok > 0:
                ax1.plot(xdum, yp, 'mo', label='identified')
                if sum(reject) > 0:
                    ax1.plot(xdum[reject], yp[reject], 'o',
                             color='tab:gray', label='ignored')
                ax1.plot(xpol, ypol, 'c-', label='fit')
            ax1.set_ylabel('polynomial - linear fit ' + r'($\AA$)')
            ax1.text(0.01, 0.99, 'CRVAL1 (' + r'$\AA$' + '):' +
                     str(round(crval1_linear, decimal_places)),
                     horizontalalignment='left',
                     verticalalignment='top',
                     transform=ax1.transAxes)
            ax1.text(0.01, 0.91, 'CDELT1 (' + r'$\AA$' + '/pixel):' +
                     str(round(cdelt1_linear, decimal_places)),
                     horizontalalignment='left',
                     verticalalignment='top',
                     transform=ax1.transAxes)
            ax1.text(0.99, 0.99, 'robust std (' + r'$\AA$' + '):' +
                     str(round(yres_summary['robust_std'], decimal_places)),
                     horizontalalignment='right',
                     verticalalignment='top',
                     transform=ax1.transAxes)
            ax1.text(0.5, 0.55, 'No. points (total / used): ' +
                     str(nlines_ok) + ' / ' +
                     str(yres_summary['npoints']),
                     horizontalalignment='center',
                     verticalalignment='top',
                     transform=ax1.transAxes)
            ax1.text(0.5, 0.4, 'Polynomial degree: ' + str(poldeg_effective),
                     horizontalalignment='center',
                     verticalalignment='bottom',
                     transform=ax1.transAxes)
            if plottitle is None:
                ax1.set_title('Refined wavelength calibration')
            else:
                ax1.set_title(plottitle)
            ax1.legend(numpoints=1, ncol=3, fancybox=True,
                       loc='lower center')
            #           bbox_to_anchor=(0.5, 1.00))

            # lower plot
            if local_ylogscale:
                spectrum = sp - sp.min() + 1.0
                spectrum = np.log10(spectrum)
            else:
                spectrum = sp.copy()
            ax2 = fig.add_subplot(grid[1, 0], sharex=ax1)
            ax2.set_xlim(1 - 0.02 * naxis1, naxis1 * 1.02)
            if local_ylogscale:
                ymin = spectrum[ixpeaks].min()
            else:
                ymin = spectrum.min()
            ymax = spectrum.max()
            dy = ymax - ymin
            ymin -= dy / 40.
            ymax += dy / 40.
            ax2.set_ylim(ymin, ymax)
            if xmin_previous is not None:
                ax2.set_xlim(xmin_previous, xmax_previous)
                ax2.set_ylim(ymin_previous, ymax_previous)
            ax2.plot(xpol, spectrum, '-')
            ax2.set_xlabel('pixel position (from 1 to NAXIS1)')
            if local_ylogscale:
                ax2.set_ylabel('~ log10(number of counts)')
            else:
                ax2.set_ylabel('number of counts')
            # mark peak location
            # ax2.plot(ixpeaks + 1, spectrum[ixpeaks], 'co',
            #          label="initial location")
            # ax2.plot(fxpeaks + 1, spectrum[ixpeaks], 'go',
            #          label="refined location")
            ax2.plot((fxpeaks + 1)[lines_ok], spectrum[ixpeaks][lines_ok],
                     'mo', label="identified lines")
            for i in range(len(ixpeaks)):
                if wv_verified_all_peaks[i] > 0:
                    ax2.text(fxpeaks[i] + 1.0, spectrum[ixpeaks[i]],
                             str(wv_verified_all_peaks[i]) +
                             '(' + str(i + 1) + ')',
                             fontsize=8,
                             horizontalalignment='center')
                else:
                    ax2.text(fxpeaks[i] + 1.0, spectrum[ixpeaks[i]],
                             '(' + str(i + 1) + ')',
                             fontsize=8,
                             horizontalalignment='center')
            # display expected location of lines in master file
            for i in range(len(wv_master)):
                tempol = poly_refined.copy()
                tempol.coef[0] -= wv_master[i]
                # compute roots
                tmproots = tempol.roots()
                # select real solutions
                tmproots = tmproots.real[abs(tmproots.imag) < 1e-5]
                # choose values within valid channel range
                tmproots = tmproots[(tmproots >= 1) * (tmproots <= naxis1)]
                if len(tmproots) > 0:
                    ax2.plot([tmproots[0], tmproots[0]], [ymin, ymax],
                             color='grey', linestyle='dotted')
            # legend
            ax2.legend(numpoints=1)

            if pdf is not None:
                pdf.savefig()
            else:
                if local_debugplot in [-22, -12, 12, 22]:
                    pause_debugplot(
                        debugplot=local_debugplot,
                        optional_prompt='Zoom/Unzoom or ' +
                                        'press RETURN to continue...',
                        tight_layout=False,
                        pltshow=True
                    )
                else:
                    pause_debugplot(debugplot=local_debugplot,
                                    tight_layout=False, pltshow=True)

            # request next action in interactive session
            if interactive:
                nrepeat_eff = 1
                print('Recalibration menu')
                print('------------------')
                print('[i] (i)nsert new peak and restart')
                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('[e] (e)valuate current polynomial at a given pixel')
                print('[w] replot (w)hole spectrum')
                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='adeilrwx')
                if ioption == 'd':
                    wv_verified_all_peaks = np.zeros(npeaks)
                elif ioption == 'r':
                    delta_wv_max = ntimes_match_wv * cdelt1_linear
                    wv_expected_all_peaks = poly_initial(fxpeaks + 1.0)
                    wv_verified_all_peaks = match_wv_arrays(
                        wv_master,
                        wv_expected_all_peaks,
                        delta_wv_max=delta_wv_max
                    )
                elif ioption == 'a':
                    delta_wv_max = ntimes_match_wv * cdelt1_linear
                    wv_verified_all_peaks = match_wv_arrays(
                        wv_master,
                        poly_refined(fxpeaks + 1.0),
                        delta_wv_max=delta_wv_max
                    )
                elif ioption == 'l':
                    xmin_previous = None
                    xmax_previous = None
                    ymin_previous = None
                    ymax_previous = None
                    if local_ylogscale:
                        local_ylogscale = False
                    else:
                        local_ylogscale = True
                elif ioption == 'e':
                    pixel = 1
                    while pixel != 0:
                        pixel = readf("Pixel coordinate (0=exit)",
                                      default=0)
                        print("--> Wavelength:", poly_refined(pixel))
                elif ioption == 'i':
                    ipixel = 1
                    # include new peaks
                    while ipixel != 0:
                        ipixel = readi("Closest pixel coordinate (integer) "
                                       "to insert peak (0=exit)",
                                       default=0, minval=0, maxval=naxis1)
                        if ipixel > 0:
                            ixpeaks = np.concatenate((ixpeaks,
                                                      np.array([ipixel-1])))
                    # sort updated array
                    ixpeaks.sort()
                    npeaks = len(ixpeaks)
                    # refine line peak locations
                    fxpeaks, sxpeaks = refine_peaks_spectrum(
                        sp, ixpeaks,
                        nwinwidth=nwinwidth_refined,
                        method="gaussian"
                    )
                    # expected wavelength of all identified peaks
                    delta_wv_max = ntimes_match_wv * cdelt1_linear
                    wv_verified_all_peaks = match_wv_arrays(
                        wv_master,
                        poly_refined(fxpeaks + 1.0),
                        delta_wv_max=delta_wv_max
                    )
                elif ioption == 'w':
                    xmin_previous = None
                    xmax_previous = None
                    ymin_previous = None
                    ymax_previous = None
                elif ioption == 'x':
                    loop = False
                else:
                    print(wv_master)
                    expected_value = \
                        poly_refined(fxpeaks[ioption - 1] + 1.0)
                    print('>>> Current expected wavelength for line #' +
                          str(ioption) + ": ", expected_value)
                    delta_wv_max = ntimes_match_wv * cdelt1_linear
                    close_value = match_wv_arrays(
                        wv_master,
                        np.array([expected_value]),
                        delta_wv_max=delta_wv_max)
                    newvalue = readf('New value for line #' + str(ioption) +
                                     ' (0 to delete line)',
                                     default=close_value[0])
                    wv_verified_all_peaks[ioption - 1] = newvalue

            else:
                loop = False

        else:
            loop = False

    # if effective degree of poly_refined < poldeg, add zeros
    if poldeg_effective < poldeg:
        numzeros = poldeg - poldeg_effective
        final_coefficients = np.concatenate((poly_refined.coef,
                                             np.zeros(numzeros)))
        poly_refined = np.polynomial.Polynomial(final_coefficients)

    if abs(local_debugplot) >= 10:
        print(">>> Initial coefficients:\n", poly_initial.coef)
        print(">>> Refined coefficients:")
        for cdum in poly_refined.coef:
            print(cdum)
        print(">>> Final CRVAL1 linear scale............:", crval1_linear)
        print(">>> Final CDELT1 linear scale............:", cdelt1_linear)

    return poly_refined, yres_summary