diff --git a/geminidr/gmos/parameters_gmos_longslit.py b/geminidr/gmos/parameters_gmos_longslit.py index 4592dff86..4a958c901 100644 --- a/geminidr/gmos/parameters_gmos_longslit.py +++ b/geminidr/gmos/parameters_gmos_longslit.py @@ -31,21 +31,37 @@ def setDefaults(self): class makeSlitIllumConfig(config.Config): - bins = config.Field("Total number of bins across the dispersion axis.", - int, None, optional=True) - border = config.Field("Size of the border added to the reconstructed slit illumination image", - int, 0, optional=True) - debug_plot = config.Field("Create diagnosis plots?", - bool, False, optional=True) - smooth_order = config.Field("Spline order to smooth binned data", - int, 3, optional=True) + bins = config.Field("Either total number of bins across the dispersion axis, " + "or a comma-separated list \n " + "of pixel coordinate pairs defining then dispersion bins, e.g. 1:300,301:500", + (int, str), None, optional=True) + regions = config.Field("Sample regions along the slit", str, None, optional=True) suffix = config.Field("Filename suffix", str, "_slitIllum", optional=True) - x_order = config.Field("Order of the x-component of the Chebyshev2D model used to reconstruct data", - int, 4, optional=True) - y_order = config.Field("Order of the y-component of the Chebyshev2D model used to reconstruct data", - int, 4, optional=True) - + function = config.ChoiceField("Fitting function to use for bin fitting (spatial direction)", str, + allowed={"spline3": "Cubic spline", + "chebyshev": "Chebyshev polynomial", + "legendre": "Legendre polynomial", + "spline1": "Linear spline"}, + default="spline3", optional=False) + order = config.Field("Order of the bin fitting function", + int, 20, optional=True) + hsigma = config.RangeField("High rejection threshold (sigma) of the bin fit", + float, 3., min=0) + lsigma = config.RangeField("Low rejection threshold (sigma) of the bin fit", + float, 3., min=0) + niter = config.RangeField("Maximum number of iterations", + int, 3, min = 0, optional=True) + grow = config.RangeField("Growth radius for rejected pixels of the bin fit", + float, 0, min=0, optional=True) + interp_order = config.RangeField("Order of the spline interpolator", + int, 3, min=1, max=5, optional=True) + debug_boundary_ext = config.Field("Extrapolate outside interpolation interval? If 'False'" + "boundary values are used", bool, False, optional=True) + interactive = config.Field("Set to activate an interactive preview to fine tune the input parameters", + bool, False, optional=True) + border = config.Field("Size of the border added to the reconstructed slit illumination image", + int, 2, optional=True) class normalizeFlatConfig(config.core_1Dfitting_config): suffix = config.Field("Filename suffix", str, "_normalized", optional=True) diff --git a/geminidr/gmos/primitives_gmos_longslit.py b/geminidr/gmos/primitives_gmos_longslit.py index 96ed82d26..513c43e17 100644 --- a/geminidr/gmos/primitives_gmos_longslit.py +++ b/geminidr/gmos/primitives_gmos_longslit.py @@ -4,35 +4,33 @@ # primtives_gmos_longslit.py # ------------------------------------------------------------------------------ -from copy import copy, deepcopy +from copy import deepcopy from importlib import import_module from itertools import groupby +import re + from gempy.library.config import RangeField import astrodata import geminidr import numpy as np from scipy.signal import correlate +from scipy.interpolate.fitpack2 import InterpolatedUnivariateSpline from astrodata.provenance import add_provenance -from astropy import visualization as vis -from astropy.modeling import models, fitting +from astropy.modeling import models from geminidr.gemini.lookups import DQ_definitions as DQ from gempy.gemini import gemini_tools as gt from gempy.library.fitting import fit_1D -from gempy.library import astromodels, transform +from gempy.library import transform from gempy.library import astrotools as at from gwcs import coordinate_frames from gwcs.wcs import WCS as gWCS -from matplotlib import gridspec -from matplotlib import pyplot as plt -from mpl_toolkits.axes_grid1 import make_axes_locatable - from recipe_system.utils.decorators import parameter_override, capture_provenance from recipe_system.utils.md5 import md5sum @@ -43,6 +41,7 @@ # ------------------------------------------------------------------------------ from ..interactive.fit import fit1d +from ..interactive.fit import bineditor from ..interactive.fit.help import NORMALIZE_FLAT_HELP_TEXT from ..interactive.interactive import UIParameters @@ -63,7 +62,6 @@ def __new__(cls, adinputs, **kwargs): " specifying 'adinputs'. Please instantiate either" "'GMOSClassicLongslit' or 'GMOSNSLongslit' instead.") - @parameter_override @capture_provenance class GMOSClassicLongslit(GMOSSpect): @@ -246,11 +244,12 @@ def addIllumMaskToDQ(self, adinputs=None, suffix=None, illum_mask=None, def makeSlitIllum(self, adinputs=None, **params): """ Makes the processed Slit Illumination Function by binning a 2D - spectrum along the dispersion direction, fitting a smooth function - for each bin, fitting a smooth 2D model, and reconstructing the 2D - array using this last model. + spectrum in wavelength, averaging data points withing each bin along dispersion + axis, then fitting a smooth function to each bin in spatial direction. + Each fitted function is then normalized to slit center, and a 2D slit illumination + function is created by interpolating between dispersion points for each row. - Its implementation based on the IRAF's `noao.twodspec.longslit.illumination` + The implementation is based on the IRAF's `noao.twodspec.longslit.illumination` task following the algorithm described in [Valdes, 1968]. It expects an input calibration image to be an a dispersed image of the @@ -264,30 +263,49 @@ def makeSlitIllum(self, adinputs=None, **params): adinputs : list List of AstroData objects containing the dispersed image of the slit of a source free of illumination problems. The data needs to - have been overscan and bias corrected and is expected to have a + have been overscan and bias corrected and it is expected to have a Data Quality mask. - bins : {None, int}, optional - Total number of bins across the dispersion axis. If None, - the number of bins will match the number of extensions on each - input AstroData object. It it is an int, it will create N bins - with the same size. - border : int, optional + bins : int/str/None + Either an integer number of equally-spaced bins covering the whole dispersion + range, or a comma-separated list of pixel coordinate pairs defining the + dispersion bins to be used for the slit profile fitting. If None, the number + of bins is max of 12 or the number of extensions. + regions : str/None + Sample region(s) selected along the spatial axis to use for fitting each dispersion + bin, as a comma-separated list of pixel ranges. Any pixels outside these + ranges will be ignored when fitting each dispersion bin. + function : str/None + Type of function to fit along the the spatial axis (for dispersion bin fitting). + None => "spline3". + order : int/None + Order of the bin fitting function (None => 20) + lsigma : float/None + Lower rejection limit in standard deviations of the bin fit (None => 3). + hsigma : float/None + Upper rejection limit in standard deviations of the bin fit (None => 3). + niter : int/None + Maximum number of rejection iterations of the bin fit (None => 3). + grow : float/None + Growth radius for rejected pixels of the bin fit (None => 0). + interp_order: int/None + Order of the spline interpolator (1 <= interp_order <= 5). None => 3. + debug_boundary_ext: bool/None + Controls the extrapolation mode for the elements outside the interpolation interval + (before the first bin center and after the last bin center). If set to "False" (default + value), the row ends are set to constant values equal to the interpolation boundary values. + If set to "True", the row ends are extrapolated. + interactive : bool/None + Activate the interactive preview of bin selection and bin fitting steps to fine tune the + input parameters (None => "False") + border : int/None Border size that is added on every edge of the slit illumination - image before cutting it down to the input AstroData frame. - smooth_order : int, optional - Order of the spline that is used in each bin fitting to smooth - the data (Default: 3) - x_order : int, optional - Order of the x-component in the Chebyshev2D model used to - reconstruct the 2D data from the binned data. - y_order : int, optional - Order of the y-component in the Chebyshev2D model used to - reconstruct the 2D data from the binned data. + image before cutting it down to the input AstroData frame (None => 2) Return ------ - List of AstroData : containing an AstroData with the Slit Illumination - Response Function for each of the input object. + List of AstroData : + Contains an AstroData with the Slit Illumination Response Function for each + of the input object. References ---------- @@ -302,14 +320,13 @@ def makeSlitIllum(self, adinputs=None, **params): suffix = params["suffix"] bins = params["bins"] border = params["border"] - debug_plot = params["debug_plot"] - smooth_order = params["smooth_order"] - cheb2d_x_order = params["x_order"] - cheb2d_y_order = params["y_order"] - + interactive_reduce = params["interactive"] + spat_params = fit_1D.translate_params(params) + interp_order = params["interp_order"] + boundary_ext = params["debug_boundary_ext"] ad_outputs = [] - for ad in adinputs: + for ad in adinputs: if len(ad) > 1 and "mosaic" not in ad[0].wcs.available_frames: log.info('Add "mosaic" gWCS frame to input data') @@ -321,9 +338,7 @@ def makeSlitIllum(self, adinputs=None, **params): log.info("Temporarily mosaicking multi-extension file") mosaicked_ad = transform.resample_from_wcs( ad, "mosaic", attributes=None, order=1, process_objcat=False) - else: - log.info('Input data already has one extension and has a ' '"mosaic" frame.') @@ -344,67 +359,156 @@ def makeSlitIllum(self, adinputs=None, **params): std = np.sqrt(variance) # Easier to work with log.info("Creating bins for data and variance") + height = data.shape[0] width = data.shape[1] - if bins is None: nbins = max(len(ad), 12) bin_limits = np.linspace(0, height, nbins + 1, dtype=int) + bin_list = list(zip(bin_limits[:-1], bin_limits[1:])) elif isinstance(bins, int): nbins = bins bin_limits = np.linspace(0, height, nbins + 1, dtype=int) + bin_list = list(zip(bin_limits[:-1], bin_limits[1:])) else: - # ToDo: Handle input bins as array - raise TypeError("Expected None or Int for `bins`. " - "Found: {}".format(type(bins))) + bin_list = _parse_user_bins(bins, height) + nbins = len(bin_list) - bin_top = bin_limits[1:] - bin_bot = bin_limits[:-1] - binned_data = np.zeros_like(data) - binned_std = np.zeros_like(std) + if ad.filename: + filename_info = ad.filename + else: + filename_info = '' + + bin_parameters = { + 'nbins': nbins, + 'bin_list': bin_list, + 'height': height + } + # Interactive interface for bin inspection and editing + if interactive_reduce: + model = { + 'x': np.arange(height), + 'y': np.ma.mean(data, axis=1), + 'regions': [(left, right) for (left, right) in + bin_list] + } + visualizer = bineditor.BinVisualizer(model, domain=[0, height], + title='Set Dispersion Bins', + primitive_name='makeSlitIllum', + central_plot=False, + bin_parameters=bin_parameters, + filename_info=filename_info) + geminidr.interactive.server.interactive_fitter(visualizer) + bin_list = _parse_user_bins(''.join(visualizer.results().split())) + nbins = len(bin_list) + + cols_val = np.arange(-border, height+border) + rows_val = np.arange(-border, width+border) + binned_shape = (nbins, len(rows_val)) + bin_data_avg = np.ma.empty((nbins, width)) + bin_std_avg = np.ma.empty((nbins, width)) + bin_data_fits = np.ma.zeros(binned_shape) + bin_std_fits = np.ma.zeros(binned_shape) + for i, bin in enumerate(bin_list): + bin_data_avg[i] = np.ma.mean(data[bin[0]:bin[1]], axis=0) + bin_std_avg[i] = np.ma.mean(std[bin[0]:bin[1]], axis=0) + spat_fitting_pars = [] + for _ in range(nbins): + spat_fitting_pars.append(spat_params) + spat_fit_points = np.arange(width) log.info("Smooth binned data and variance, and normalize them by " "smoothed central value") - for bin_idx, (b0, b1) in enumerate(zip(bin_bot, bin_top)): - rows = np.arange(width) + config = self.params[self.myself()] - avg_data = np.ma.mean(data[b0:b1], axis=0) - model_1d_data = astromodels.UnivariateSplineWithOutlierRemoval( - rows, avg_data, order=smooth_order) + # Interactive interface for fitting dispersion bins + if interactive_reduce: + all_pixels = [] + all_domains = [] + for _ in range(nbins): + all_pixels.append(spat_fit_points) + all_domains.append([-border, width + border]) - avg_std = np.ma.mean(std[b0:b1], axis=0) - model_1d_std = astromodels.UnivariateSplineWithOutlierRemoval( - rows, avg_std, order=smooth_order) + data_with_weights = {"x": [], "y": [], "weights": []} + for rppixels, avg_data, avg_std in zip(all_pixels, bin_data_avg, bin_std_avg): + data_with_weights["x"].append(rppixels) + data_with_weights["y"].append(avg_data) + data_with_weights["weights"].append(at.divide0(1., avg_std)) - slit_central_value = model_1d_data(rows)[width // 2] - binned_data[b0:b1] = model_1d_data(rows) / slit_central_value - binned_std[b0:b1] = model_1d_std(rows) / slit_central_value + config.update(**params) + uiparams = UIParameters(config) - log.info("Reconstruct 2D mosaicked data") - bin_center = np.array(0.5 * (bin_bot + bin_top), dtype=int) - cols_fit, rows_fit = np.meshgrid(np.arange(width), bin_center) + x_label = "Rows" if dispaxis == 1 else "Columns" + first_label = 'cols' if dispaxis == 1 else 'rows' + tab_names = [f'[{start+1}:{end}]' for (start, end) in bin_list] + tab_names[0] = f"Mean of {first_label} {tab_names[0]}" - fitter = fitting.SLSQPLSQFitter() - model_2d_init = models.Chebyshev2D( - x_degree=cheb2d_x_order, x_domain=(0, width), - y_degree=cheb2d_y_order, y_domain=(0, height)) + visualizer = fit1d.Fit1DVisualizer(data_with_weights, spat_fitting_pars, + tab_names=tab_names, + xlabel=x_label, ylabel='Counts', + domains=all_domains, + title="Make Slit Illumination Function", + primitive_name="makeSlitIllum", + filename_info=filename_info, + enable_user_masking=True, + enable_regions=True, + # help_text=NORMALIZE_FLAT_HELP_TEXT, + recalc_inputs_above=True, + modal_message="Recalculating", + pad_buttons=True, + ui_params=uiparams) + geminidr.interactive.server.interactive_fitter(visualizer) + fits = visualizer.results() + spat_fitting_pars = [fit.extract_params() for fit in fits] + for bin_idx, bin_fit in enumerate(fits): + bin_data_fit = bin_fit.evaluate(rows_val) + bin_data_fits[bin_idx,:] = bin_data_fit - model_2d_data = fitter(model_2d_init, cols_fit, rows_fit, - binned_data[rows_fit, cols_fit]) + else: + for bin_idx, (fit_params, avg_data, avg_std) in \ + enumerate(zip(spat_fitting_pars, bin_data_avg, bin_std_avg)): + bin_data_fit = fit_1D(avg_data, + points=spat_fit_points, + weights=at.divide0(1., avg_std), + domain=(-border, width+border), + **fit_params, axis=0).evaluate(rows_val) + bin_data_fits[bin_idx,:] = bin_data_fit + + for bin_idx, (fit_params, avg_std, bin_data_fit) in \ + enumerate(zip(spat_fitting_pars, bin_std_avg, bin_data_fits)): + bin_std_fit = fit_1D(avg_std, + points=spat_fit_points, + **fit_params, + axis=0).evaluate(rows_val) + slit_central_value = bin_data_fit[(width+2*border) // 2] + bin_data_fits[bin_idx,:] = bin_data_fit / slit_central_value + bin_std_fits[bin_idx,:] = bin_std_fit / slit_central_value - model_2d_std = fitter(model_2d_init, cols_fit, rows_fit, - binned_std[rows_fit, cols_fit]) + log.info("Reconstruct 2D mosaicked data") - rows_val, cols_val = \ - np.mgrid[-border:height+border, -border:width+border] + bin_center = np.array([0.5 * (bin_start + bin_end) for (bin_start, bin_end) in bin_list], + dtype=int) + slit_response_data = np.zeros((len(cols_val), len(rows_val))) + slit_response_std = np.zeros((len(cols_val), len(rows_val))) + + # Interpolation between dispersion points along the rows + if nbins > 1: + for k, (data_row, std_row) in enumerate(zip(bin_data_fits.T, bin_std_fits.T)): + # Set extrapolated row ends to interpolation boundary value, or extrapolate + ext_val = 0 if boundary_ext else 3 + f1 = InterpolatedUnivariateSpline(bin_center, data_row, k=interp_order, ext=ext_val) + f2 = InterpolatedUnivariateSpline(bin_center, std_row, k=interp_order, ext=ext_val) + slit_response_data[:,k] = f1(cols_val) + slit_response_std[:,k] = f2(cols_val) + + # If there is only one bin, copy the slit profile to each column + elif nbins == 1: + slit_response_data[:] = bin_data_fits + slit_response_std[:] = bin_std_fits - slit_response_data = model_2d_data(cols_val, rows_val) - slit_response_mask = np.pad(mask, border, mode='edge') # ToDo: any update to the mask? - slit_response_std = model_2d_std(cols_val, rows_val) slit_response_var = slit_response_std ** 2 - - del cols_fit, cols_val, rows_fit, rows_val + slit_response_mask = np.pad(mask, border, mode='edge') _data, _mask, _variance = at.transpose_if_needed( slit_response_data, slit_response_mask, slit_response_var, @@ -418,7 +522,7 @@ def makeSlitIllum(self, adinputs=None, **params): if "mosaic" in ad[0].wcs.available_frames: - log.info("Map coordinates between slit function and mosaicked data") # ToDo: Improve message? + log.info("Map coordinates between slit function and mosaicked data") # ToDo: Improve message? slit_response_ad = _split_mosaic_into_extensions( ad, slit_response_ad, border_size=border) @@ -440,186 +544,6 @@ def makeSlitIllum(self, adinputs=None, **params): slit_response_ad.update_filename(suffix=suffix, strip=True) ad_outputs.append(slit_response_ad) - # Plotting ------ - if debug_plot: - - log.info("Creating plots") - palette = copy(plt.cm.cividis) - palette.set_bad('r', 0.75) - - norm = vis.ImageNormalize(data[~data.mask], - stretch=vis.LinearStretch(), - interval=vis.PercentileInterval(97)) - - fig = plt.figure( - num="Slit Response from MEF - {}".format(ad.filename), - figsize=(12, 9), dpi=110) - - gs = gridspec.GridSpec(nrows=2, ncols=3, figure=fig) - - # Display raw mosaicked data and its bins --- - ax1 = fig.add_subplot(gs[0, 0]) - im1 = ax1.imshow(data, cmap=palette, origin='lower', - vmin=norm.vmin, vmax=norm.vmax) - - ax1.set_title("Mosaicked Data\n and Spectral Bins", fontsize=10) - ax1.set_xlim(-1, data.shape[1]) - ax1.set_xticks([]) - ax1.set_ylim(-1, data.shape[0]) - ax1.set_yticks(bin_center) - ax1.tick_params(axis=u'both', which=u'both', length=0) - - ax1.set_yticklabels( - ["Bin {}".format(i) for i in range(len(bin_center))], - fontsize=6) - - _ = [ax1.spines[s].set_visible(False) for s in ax1.spines] - _ = [ax1.axhline(b, c='w', lw=0.5) for b in bin_limits] - - divider = make_axes_locatable(ax1) - cax1 = divider.append_axes("right", size="5%", pad=0.05) - plt.colorbar(im1, cax=cax1) - - # Display non-smoothed bins --- - ax2 = fig.add_subplot(gs[0, 1]) - im2 = ax2.imshow(binned_data, cmap=palette, origin='lower') - - ax2.set_title("Binned, smoothed\n and normalized data ", fontsize=10) - ax2.set_xlim(0, data.shape[1]) - ax2.set_xticks([]) - ax2.set_ylim(0, data.shape[0]) - ax2.set_yticks(bin_center) - ax2.tick_params(axis=u'both', which=u'both', length=0) - - ax2.set_yticklabels( - ["Bin {}".format(i) for i in range(len(bin_center))], - fontsize=6) - - _ = [ax2.spines[s].set_visible(False) for s in ax2.spines] - _ = [ax2.axhline(b, c='w', lw=0.5) for b in bin_limits] - - divider = make_axes_locatable(ax2) - cax2 = divider.append_axes("right", size="5%", pad=0.05) - plt.colorbar(im2, cax=cax2) - - # Display reconstructed slit response --- - vmin = slit_response_data.min() - vmax = slit_response_data.max() - - ax3 = fig.add_subplot(gs[1, 0]) - im3 = ax3.imshow(slit_response_data, cmap=palette, - origin='lower', vmin=vmin, vmax=vmax) - - ax3.set_title("Reconstructed\n Slit response", fontsize=10) - ax3.set_xlim(0, data.shape[1]) - ax3.set_xticks([]) - ax3.set_ylim(0, data.shape[0]) - ax3.set_yticks([]) - ax3.tick_params(axis=u'both', which=u'both', length=0) - _ = [ax3.spines[s].set_visible(False) for s in ax3.spines] - - divider = make_axes_locatable(ax3) - cax3 = divider.append_axes("right", size="5%", pad=0.05) - plt.colorbar(im3, cax=cax3) - - # Display extensions --- - ax4 = fig.add_subplot(gs[1, 1]) - ax4.set_xticks([]) - ax4.set_yticks([]) - _ = [ax4.spines[s].set_visible(False) for s in ax4.spines] - - sub_gs4 = gridspec.GridSpecFromSubplotSpec( - nrows=len(ad), ncols=1, subplot_spec=gs[1, 1], hspace=0.03) - - # The [::-1] is needed to put the fist extension in the bottom - for i, ext in enumerate(slit_response_ad[::-1]): - - ext_data, ext_mask, ext_variance = at.transpose_if_needed( - ext.data, ext.mask, ext.variance, transpose=dispaxis == 1) - - ext_data = np.ma.masked_array(ext_data, mask=ext_mask) - - sub_ax = fig.add_subplot(sub_gs4[i]) - - im4 = sub_ax.imshow(ext_data, origin="lower", vmin=vmin, - vmax=vmax, cmap=palette) - - sub_ax.set_xlim(0, ext_data.shape[1]) - sub_ax.set_xticks([]) - sub_ax.set_ylim(0, ext_data.shape[0]) - sub_ax.set_yticks([ext_data.shape[0] // 2]) - - sub_ax.set_yticklabels( - ["Ext {}".format(len(slit_response_ad) - i - 1)], - fontsize=6) - - _ = [sub_ax.spines[s].set_visible(False) for s in sub_ax.spines] - - if i == 0: - sub_ax.set_title("Multi-extension\n Slit Response Function") - - divider = make_axes_locatable(ax4) - cax4 = divider.append_axes("right", size="5%", pad=0.05) - plt.colorbar(im4, cax=cax4) - - # Display Signal-To-Noise Ratio --- - snr = data / np.sqrt(variance) - - norm = vis.ImageNormalize(snr[~snr.mask], - stretch=vis.LinearStretch(), - interval=vis.PercentileInterval(97)) - - ax5 = fig.add_subplot(gs[0, 2]) - - im5 = ax5.imshow(snr, cmap=palette, origin='lower', - vmin=norm.vmin, vmax=norm.vmax) - - ax5.set_title("Mosaicked Data SNR", fontsize=10) - ax5.set_xlim(-1, data.shape[1]) - ax5.set_xticks([]) - ax5.set_ylim(-1, data.shape[0]) - ax5.set_yticks(bin_center) - ax5.tick_params(axis=u'both', which=u'both', length=0) - - ax5.set_yticklabels( - ["Bin {}".format(i) for i in range(len(bin_center))], - fontsize=6) - - _ = [ax5.spines[s].set_visible(False) for s in ax5.spines] - _ = [ax5.axhline(b, c='w', lw=0.5) for b in bin_limits] - - divider = make_axes_locatable(ax5) - cax5 = divider.append_axes("right", size="5%", pad=0.05) - plt.colorbar(im5, cax=cax5) - - # Display Signal-To-Noise Ratio of Slit Illumination --- - slit_response_snr = np.ma.masked_array( - slit_response_data / np.sqrt(slit_response_var), - mask=slit_response_mask) - - ax6 = fig.add_subplot(gs[1, 2]) - - im6 = ax6.imshow(slit_response_snr, origin="lower", - vmin=norm.vmin, vmax=norm.vmax, cmap=palette) - - ax6.set_xlim(0, slit_response_snr.shape[1]) - ax6.set_xticks([]) - ax6.set_ylim(0, slit_response_snr.shape[0]) - ax6.set_yticks([]) - ax6.set_title("Reconstructed\n Slit Response SNR") - - _ = [ax6.spines[s] .set_visible(False) for s in ax6.spines] - - divider = make_axes_locatable(ax6) - cax6 = divider.append_axes("right", size="5%", pad=0.05) - plt.colorbar(im6, cax=cax6) - - # Save plots --- - fig.tight_layout(rect=[0, 0, 0.95, 1], pad=0.5) - fname = slit_response_ad.filename.replace(".fits", ".png") - log.info("Saving plots to {}".format(fname)) - plt.savefig(fname) - return ad_outputs def normalizeFlat(self, adinputs=None, **params): @@ -895,6 +819,11 @@ def slitIllumCorrect(self, adinputs=None, slit_illum=None, "slit illumination file: \n{}".format(ad.filename, slit_illum_ad.filename)) ad_out = deepcopy(ad) + + + log.stdinfo(f"{ad.filename}: dividing by the slit illumination function " + f"{slit_illum_ad.filename}") + ad_out.divide(slit_illum_ad) # Update the header and filename, copying QECORR keyword from flat @@ -1027,6 +956,69 @@ def _split_mosaic_into_extensions(ref_ad, mos_ad, border_size=0): return ad_out +def _parse_user_bins(bins, frame_size:int=None): + """ + Parse a string of bin ranges containing a comma-separated list of colon- or hyphen-separated + pixel sections into a list of Python tuples. Arrange the bins in ascending order, remove the ones outside + the frame pixel range, merge the overlapping bins. + + Parameters + ---------- + bins: str + Comma-separated list of colon- or hyphen-separated pixel ranges + frame_size: int + The size of the frame dimension along which the bins are selected. + + Returns + ------- + A sorted list of 2-value tuples lying within the specified frame range, with no overlapping. + """ + bin_list = [] + for bin in re.split(",|;| ", bins.strip("[]()'")): + bin = bin.strip("()[]' ") + bin_limits = re.split(":|-", bin) + if not len(bin_limits) == 2: + raise TypeError("Bin limits must be specified as comma-separated list " + "of colon- or hyphen-separated pixel sections, e.g. 1:300,301:500") + int_bin_limits = [] + for bin_limit in bin_limits: + try: + int_bin_limit = int(bin_limit) + except ValueError: + raise TypeError("Bin ranges must be integer") + if frame_size is not None: + int_bin_limits.append(min(int_bin_limit, frame_size)) + else: + int_bin_limits.append(int_bin_limit) + int_bin_limits.sort() + + # trim off the bins that are outside the frame pixel range + if frame_size is not None and int_bin_limits[0] == frame_size: + break + + bin_list.append(tuple(int_bin_limits)) + bin_list = sorted(bin_list, key=lambda tup: tup[0]) + + # merge overlapping bins + merged_list = [bin_list[0]] + for bin in bin_list[1:]: + last = merged_list[-1] + if last[1] > bin[0]: + print("Merging the overlapping bins") + merged_list[-1] = (last[0], max(last[1], bin[1])) + elif last[1] == bin[0]: + merged_list.append((bin[0]+1, bin[1])) + else: + merged_list.append(bin) + + adjusted_bin_list = [] + for i, bin in enumerate(merged_list): + if i == 0 and bin[0] == 0: + adjusted_bin_list.append(bin) + else: + adjusted_bin_list.append((bin[0]-1, bin[1])) + return adjusted_bin_list + @parameter_override @capture_provenance diff --git a/geminidr/gmos/tests/longslit/test_make_slit_illum.py b/geminidr/gmos/tests/longslit/test_make_slit_illum.py index 3fdb94003..e1d26f771 100755 --- a/geminidr/gmos/tests/longslit/test_make_slit_illum.py +++ b/geminidr/gmos/tests/longslit/test_make_slit_illum.py @@ -1,8 +1,6 @@ #!/usr/bin/env python """ -Tests for the `makeSlitIllum` primitive. The primitive itself is -defined in :mod:`~geminidr.core.primitives_spect` but these tests use GMOS Spect -data. +Tests for the `makeSlitIllum` primitive. """ import os import pytest @@ -49,7 +47,7 @@ def test_create_slit_illumination_with_mosaicked_data(ad, change_working_dir, re constant. There are several ways of doing this but, given the noise levels, we bin - the data, fit a polynomium, and check that the fitted polynomium has its 1st + the data, fit a polynomial, and check that the fitted polynomial has its 1st and 2nd coefficients almost zero. """ plot = request.config.getoption("--do-plots") @@ -63,7 +61,7 @@ def test_create_slit_illumination_with_mosaicked_data(ad, change_working_dir, re assert hasattr(ad[0], "wcs") p = primitives_gmos_longslit.GMOSLongslit([ad]) - p.makeSlitIllum(bins=25, border=10, debug_plot=plot) + p.makeSlitIllum(bins=25, border=10) slit_illum_ad = p.writeOutputs( suffix="_mosaickedSlitIllum", strip=True)[0] @@ -122,7 +120,7 @@ def test_create_slit_illumination_with_multi_extension_data(ad, change_working_d cwd = os.getcwd() print("Running tests inside folder:\n {}".format(cwd)) p = primitives_gmos_longslit.GMOSLongslit([ad]) - p.makeSlitIllum(bins=25, border=10, debug_plot=plot) + p.makeSlitIllum(bins=25, border=10) slit_illum_ad = p.writeOutputs()[0] for ext, slit_ext in zip(ad, slit_illum_ad): diff --git a/geminidr/interactive/fit/bineditor.py b/geminidr/interactive/fit/bineditor.py new file mode 100644 index 000000000..6df8a1fb4 --- /dev/null +++ b/geminidr/interactive/fit/bineditor.py @@ -0,0 +1,586 @@ +import numpy as np + +from bokeh.layouts import column, row +from bokeh.models import (Button, Div, Spacer, NumericInput) +from bokeh import models as bm +from bokeh.plotting import figure + +from geminidr.interactive.controls import Controller +from geminidr.interactive.fit.fit1d import Fit1DRegionListener, InteractiveModel1D +from geminidr.interactive.interactive import (PrimitiveVisualizer, RegionEditor, + GIRegionModel, connect_region_model) +from gempy.utils import logutils + +log = logutils.get_logger(__name__) + +DETAILED_HELP = """ + +
Interface to inspect and edit dispersion bins.
+Allows the user to generate a number of equally-spaced bins, and to modify +bin limits, either manually entering them, or using a point-and-click interface.
+""" + +def tuples_to_slices(tuples): + """ + Helping function to translate into slices a list of bin limit expressed + as tuples (start, end) + """ + return [slice(start, end) for (start, end) in tuples] + +def bin_figure(width=None, height=None, xpoint='x', ypoint='y', + xlabel=None, ylabel=None, model=None): + """ + Function to produce bokeh objects for the bin editing plot. + Listeners are not added here. + + Parameters + ---------- + width : int + width of the plots + height : int + height of the main plot (ratios/residuals are half-height) + xpoint, ypoint : str + column names in model.data containing x and y data for points + xlabel, ylabel : str + label for axes of main plot + model : InteractiveModel1D + object containing the fit information + + Returns + ------- + Figure + The plotting figure + """ + + tools = "pan,wheel_zoom,box_zoom,reset" + + p_main = figure(plot_width=width, plot_height=height, min_width=400, + title='Illumination bins', x_axis_label=xlabel, y_axis_label=ylabel, + tools=tools, + output_backend="webgl", x_range=None, y_range=None, + min_border_left=80) + p_main.height_policy = 'fixed' + p_main.width_policy = 'fit' + p_main.scatter(x=xpoint, y=ypoint, source=model.data, + size=5, legend_field='mask', + **model.mask_rendering_kwargs()) + + return p_main + +class BinEditor(RegionEditor): + """ + Specialized RegionEditor. Just for cosmetic changes (changes the title) + """ + def __init__(self, *args, **kw): + super().__init__(*args, **kw) + self.text_input.title = "Bin limits (i.e. 101:500,511:900,951: Press 'Enter' to apply):" + +class BinModel1D(InteractiveModel1D): + """ + Specialized InteractiveModle1D that removes the fitting functionality, as + it's not needed for the bin editor. + """ + def perform_fit(self, *args): + """ + Dummy function. Needs to be here to be compliant with the InteractiveModel1D + interface, but no calculation will be performed. + """ + ... + + def evaluate(self, x): + """ + Returns the `x` parameter itself. This model does not perform fitting evaluation. + """ + return x + +class NumericInputControl: + def __init__(self, title, value, mode='int'): + self.orig_value = value + self.title = title + self.component = NumericInput(width=64, value=value, mode=mode) + + def build(self): + return row([Div(text=self.title, align='center'), + Spacer(width_policy='max'), + self.component]) + + def update(self, value): + self.component.value = value + + def reset(self): + self.update(self.orig_value) + + @property + def value(self): + return self.component.value + +class BinResettingUI: + def __init__(self, vis, model, bin_parameters): + """ + Class to manage the set of UI controls to generate equally-sized bins over the whole dispersion + range, and to reset to original values. + + Parameters + ---------- + vis : :class:`~geminidr.interactive.fit.bineditor.BinVisualizer` + The visualizer related to these inputs + bin_parameters: dict + Initial number of bins and generated regions + fit : :class:`~geminidr.interactive.fit.fit1d.InteractiveModel1D` + The model information for doing the 1-D fit + """ + self.vis = vis + self.model = model + self.original_parameters = bin_parameters + + # self.num_input = NumericInput(width=64, value=self.original_parameters['nbins'], mode='int') + # self.number_bins = row([Div(text="Number of bins", align='center'), + # Spacer(width_policy='max'), + # self.num_input]) + + self.number_bins = NumericInputControl(title="Number of bins", + value=self.original_parameters['nbins']) + + def _generate_handler(result): + if result: + generate_button.disabled = True + def fn(): + self.generate_model_regions(self.number_bins.value) + generate_button.disabled = False + vis.do_later(fn) + + + generate_button = Button(label="Generate bins", button_type='primary', + default_size=200) + + vis.make_ok_cancel_dialog(generate_button, + 'All bin limits will be recomputed and changes will be lost. Proceed?', + _generate_handler) + reset_button = bm.Button(label="Reset", align='center', + button_type='warning', width_policy='min') + self.reset_dialog = self.vis.make_ok_cancel_dialog( + reset_button, 'Reset will change all inputs for this tab back ' + 'to their original values. Proceed?', self.reset_dialog_handler) + + self.controls_column = ( + self.number_bins.build(), + generate_button, + reset_button, + ) + + def get_bokeh_components(self): + """ + Return the bokeh components to be added with all the input widgets. + + Returns + ------- + list : :class:`~bokeh.models.layout.LayoutDOM` + List of bokeh components to add to the UI + """ + return self.controls_column + + def generate_model_regions(self, nbins): + """ + Handle the 'Generate bins' button being clicked. + + This will generate a new set of bin limits and update the model, which + in turn will update the interface. + """ + bin_limits = np.linspace(0, self.original_parameters['height'], nbins + 1, dtype=int) + bin_list = list(zip(bin_limits[:-1], bin_limits[1:])) + self.model.load_from_tuples(tuples_to_slices(bin_list)) + + def reset_model_regions(self): + """ + Handle the 'Reset' button being clicked. + + This will update the model with the initial bin limit list, which in + turn will udpate the interface. + """ + # self.num_input.value = self.original_parameters['nbins'] + self.number_bins.reset() + self.model.load_from_tuples(tuples_to_slices(self.original_parameters['bin_list'])) + + def reset_dialog_handler(self, result): + """ + Reset bin limits values. + + Parameters + ---------- + result : bool + This is the user response to an ok/cancel confirmation dialog. If False, do not reset. + """ + if result: + self.reset_model_regions() + +class BinPanel: + def __init__(self, visualizer, regions, bin_parameters, domain=None, + x=None, y=None, weights=None, xlabel='x', ylabel='y', + plot_width=600, plot_height=400, central_plot=True): + """ + Panel for visualizing a 1-D fit, perhaps in a tab + + Parameters + ---------- + visualizer : :class:`~geminidr.interactive.fit.fit1d.Fit1DVisualizer` + visualizer to associate with + regions : list of regions + ... + bin_parameters: dict + Initial number of bins and generated regions + domain : list of pixel coordinates + Used for new fit_1D fitter + x : :class:`~numpy.ndarray` + X coordinate values + y : :class:`~numpy.ndarray` + Y coordinate values + weights : None or :class:`~numpy.ndarray` + weights of individual points + xlabel : str + label for X axis + ylabel : str + label for Y axis + plot_width : int + width of plot area in pixels + plot_height : int + height of plot area in pixels + central_plot : bool + If True, the main plot will be on the left and the control column + on the right. If False, the opposite. + """ + # Just to get the doc later + self.visualizer = visualizer + + # self.title = "" + + self.width = plot_width + self.height = plot_height + self.xlabel = xlabel + self.ylabel = ylabel + self.xpoint = 'x' + self.ypoint = 'y' + self.p_main = None + + # Avoids having to check whether this is None all the time + band_model = GIRegionModel(domain=domain, support_adjacent=True) + self.model = BinModel1D({}, domain, x, y, weights, band_model=band_model) + self.model.add_listener(self.model_change_handler) + + self.bin_resetting_ui = BinResettingUI(visualizer, band_model, bin_parameters) + controls_column = self.bin_resetting_ui.get_bokeh_components() + # reset_button = bm.Button(label="Reset", align='center', + # button_type='warning', width_policy='min') + + # self.reset_dialog = self.visualizer.make_ok_cancel_dialog( + # reset_button, 'Reset will change all inputs for this tab back ' + # 'to their original values. Proceed?', self.reset_dialog_handler) + + controller_div = Div(margin=(20, 0, 0, 0), width=220, + style={"color": "gray", "padding": "5px"}) + controls = column(*controls_column, controller_div, + width=220) + + fig_column = self.build_figures(domain=domain, controller_div=controller_div) + + # Initializing regions here ensures the listeners are notified of the region(s) + band_model.load_from_tuples(tuples_to_slices(regions)) + + region_editor = BinEditor(band_model) + fig_column.append(region_editor.get_widget()) + col = column(*fig_column) + col.sizing_mode = 'scale_width' + + col_order = [col, controls] if central_plot else [controls, col] + self.component = row(*col_order, css_classes=["tab-content"], + spacing=10) + + def build_figures(self, domain=None, controller_div=None): + """ + Construct the figure containing the plot needed for this + Visualizer. + + Parameters + ---------- + domain : 2-tuple/None + the domain over which the model is defined + controller_div : Div + Div object accessible by Controller for updating help text + + Returns + ------- + fig_column : list + list of bokeh objects with attached listeners + """ + + p_main = bin_figure(width=self.width, height=self.height, + xpoint=self.xpoint, ypoint=self.ypoint, + xlabel=self.xlabel, ylabel=self.ylabel, model=self.model) + self.model.band_model.add_listener(Fit1DRegionListener(self.update_bin_limits)) + connect_region_model(p_main, self.model.band_model) + + Controller(p_main, None, self.model.band_model, controller_div, + mask_handlers=None, domain=domain, helpintrotext= + "While the mouse is over the upper plot, " + "choose from the following commands:") + + fig_column = [p_main] + + self.p_main = p_main + + # Do a custom padding for the ranges + self.reset_view() + + return fig_column + + def reset_view(self): + """ + This calculates the x and y ranges for the figure with some custom padding. + + This is used when initially building the figure, but also as a listener for + whenever the data changes. + """ + if not hasattr(self, 'p_main') or self.p_main is None: + # This may be a subclass, p_main is not being stored so nothing to reset + return + + x_range = None + y_range = None + try: + xdata = self.model.data.data[self.xpoint] + ydata = self.model.data.data[self.ypoint] + except (AttributeError, KeyError): + pass + else: + x_min, x_max = min(xdata), max(xdata) + if x_min != x_max: + x_pad = (x_max - x_min) * 0.1 + self.p_main.x_range.update(start=x_min - x_pad, end=x_max + x_pad * 2) + y_min, y_max = min(ydata), max(ydata) + if y_min != y_max: + y_pad = (y_max - y_min) * 0.1 + self.p_main.y_range.update(start=y_min - y_pad, end=y_max + y_pad) + if x_range is not None: + self.p_main.x_range = x_range + if y_range is not None: + self.p_main.y_range = y_range + + def update_bin_limits(self): + """ Update bin limits """ + self.model.regions = self.model.band_model.build_regions() + + def model_change_handler(self, model): + """ + If the model changes, this gets called to save the results. + + Parameters + ---------- + model : :class:`~geminidr.interactive.fit.fit1d.InteractiveModel1D` + The model that changed. + """ + # We're not evaluating fits, but we're reusing existing models so + # we'll follow the interface used elsewhere. + model.evaluation.data['model'] = model.evaluate(model.evaluation.data['xlinspace']) + +class BinVisualizer(PrimitiveVisualizer): + """ + Specialized visualizer for displaying and editing bin limits. + + Attributes: + tabs: layout containing all the stuff required for an interactive 1D fit + submit_button: the button signifying a successful end to the interactive session + + config: the Config object describing the parameters are their constraints + widgets: a dict of (param_name, widget) elements that allow the properties + of the widgets to be set/accessed by the calling primitive. So far + I'm only including the widgets in the reinit_panel + """ + def __init__(self, data_source, xlabel='Column', ylabel='Signal', + domain=None, title=None, primitive_name=None, filename_info=None, + template="fit1d.html", help_text=None, + ui_params=None, pad_buttons=False, + **kwargs): + """ + Parameters + ---------- + data_source : dict or dict-returning function + input data or the function to calculate the input data. The dict + must have keys of "x" and "y" indicating the input data values, + with each being an array or a list of arrays. There are also + optional "weights" (array of weights), "*_mask" (additional input masks), + and "meta" (any additional data) + xlabel : str + String label for X axis + ylabel : str + String label for Y axis + domain : list + Domains for the input + title : str + Title for UI (Interactive