diff --git a/doc/DRAGONS/changes.rst b/doc/DRAGONS/changes.rst index 5c387405a..003c60ffb 100644 --- a/doc/DRAGONS/changes.rst +++ b/doc/DRAGONS/changes.rst @@ -8,6 +8,49 @@ Change Logs *********** +4.2.2 +===== + +Bug fixes +--------- +** geminidr.gmos ** + +* Prevent ``addIllumMaskToDQ`` from crashing when the GMOS longslit + bridges cannot be located in the spatial profile, and instead output a + warning. + +Compatibility +------------- +A change to the behavior of Tabs in bokeh v3.9.0 caused the residual/ratio +plot in the interactive tools to be corrupted. This release ensures +compatibility with v3.9.0 as well as earlier versions. + + +4.2.1-goawavecal +================ + +Improvements +------------ + +This internal (used in GOA) release includes the following improvements to +the determineWavelengthSolution primitive: + +* Check the final fit a-posteriori to verify that it is within a factor of two + of the initial fit boundary constraints and rejects it if not +* In Science-Quality mode, if no arc lines are found for all orders being + processed, an exception is raised and processing stops. In Quick-Look mode, + the previous behavior of adding a linear wavelength solution based on the + headers is unchanged +* If only one arc line is found, a linear fit based on the headers is shifted + to match the one arc line, a warning is generated, and that fit is returned as + the wavelength solution + +This release also includes functionality in gempy to embed auxillary files +(intended for eg PDF plot figures) in binary tables in astrodata instances, and +hence FITS files. This is used for example to capture the PDFs of the arc line +identifications into the reduced arc FITS files, and a corresponding script +in gempy to extract them. + 4.2.1 ===== diff --git a/geminidr/core/parameters_spect.py b/geminidr/core/parameters_spect.py index a72897af0..c1d7e8674 100644 --- a/geminidr/core/parameters_spect.py +++ b/geminidr/core/parameters_spect.py @@ -579,6 +579,8 @@ class maskBeyondSlitConfig(config.Config): "Minimum fraction of pixel that must be illuminated to not be masked", float, 0.9, min=0., max=1., inclusiveMax=True) +class monitorWavelengthSolutionConfig(config.Config): + suffix = config.Field("Filename suffix", str, "_wavelengthSolutionMonitoring", optional=True) class normalizeFlatConfig(config.core_1Dfitting_config): suffix = config.Field("Filename suffix", str, "_normalized", optional=True) diff --git a/geminidr/core/primitives_ccd.py b/geminidr/core/primitives_ccd.py index 3570f3614..5ef471d45 100644 --- a/geminidr/core/primitives_ccd.py +++ b/geminidr/core/primitives_ccd.py @@ -225,7 +225,7 @@ def subtractOverscan(self, adinputs=None, **params): # need to match asec location and size pixels = np.arange(asec.y1, asec.y2) data = np.mean(ext.data[asec.y1:asec.y2, x1:x2], axis=axis) - stds = np.std(ext.data[asec.y1:asec.y2, x1:x2], axis=axis) + stds = np.std(ext.data[asec.y1:asec.y2, x1:x2], axis=axis, ddof=1) else: if y1 > asec.y1: # Bias on top y1 += nbiascontam @@ -240,17 +240,29 @@ def subtractOverscan(self, adinputs=None, **params): data = np.mean(ext.data[y1:y2, asec.x1:asec.x2], axis=axis) stds = np.std(ext.data[y1:y2, asec.x1:asec.x2], axis=axis) - # We have 1 sample (of x2-x1 values) from each of N - # populations each of which has a different population mean, - # (ie the bias leve for that row) but we assume they all - # have the same standard deviation (ie the read noise) and - # we want to estimate that population standard deviation. - # We've calculated the mean (data) and standard deviation - # (stds) of each of the samples, and want to estimate the - # population standard deviation. Because the means are - # different, we can't just calculate the std of all the - # samples. - overstd = np.sqrt(np.mean(stds * stds)) + # We have (x2-x1) samples from each of N populations each + # of which has a different population mean, (ie the bias + # level for that row) but we assume they all have the same + # standard deviation (ie the read noise) and we want to + # estimate that population standard deviation. + # We've calculated the mean (data) and estimates of the + # population standard deviation[*] (stds) of each of the set + # of samples, and want to estimate the population standard + # deviation. Because the means are different, we can't just + # calculate the std of all the samples. + # [*] [by which I mean to say we set ddof=1 in the call to + # np.stds to apply the (N/N-1) factor to estimate the + # population standard deviation from the calculaiton of the + # sample standard deviation. + # We use the mean of these per-row population variance + # estimates as our estimate of the global population + # variance. This seems to work well, but I'm not sure if it + # is truly an optimal estimator. + # It turns out that the first and last row of the readout + # often have pixels with massively outlying values, so we + # disregard those samples + trimmedstds = stds[1:-1] + overstd = np.sqrt(np.mean(trimmedstds * trimmedstds)) # Readnoise descriptor always returns value in electrons if ext.is_in_adu(): @@ -309,7 +321,10 @@ def subtractOverscan(self, adinputs=None, **params): ext.hdr.set('OVERSCAN', previous_overscan + bias_level, self.keyword_comments['OVERSCAN']) ext.hdr.set('OVERRMS', sigma, self.keyword_comments['OVERRMS']) - ext.hdr.set('OVERRDNS', overstd, self.keyword_comments['OVERRDNS']) + # By convention, "readnoise" is always in electrons + overrdns = overstd * ext.gain() if ext.is_in_adu() else overstd + ext.hdr.set('OVERRDNS', overrdns, self.keyword_comments['OVERRDNS']) + for desc in ('saturation_level', 'non_linear_level'): with suppress(AttributeError, KeyError): ext.hdr[ad._keyword_for(desc)] -= bias_level diff --git a/geminidr/core/primitives_spect.py b/geminidr/core/primitives_spect.py index 7273fee0f..dea0da32e 100644 --- a/geminidr/core/primitives_spect.py +++ b/geminidr/core/primitives_spect.py @@ -16,6 +16,7 @@ from importlib import import_module import matplotlib +import numpy import numpy as np from astropy import units as u from astropy.io.ascii.core import InconsistentTableError @@ -1218,6 +1219,7 @@ def _get_fit1d_input_data(ext, exptime, spec_table): primitive_name="calculateSensitivity", filename_info=filename_info, help_text=CALCULATE_SENSITIVITY_HELP_TEXT, + plot_ratios=False, ui_params=uiparams) geminidr.interactive.server.interactive_fitter(visualizer) @@ -2909,22 +2911,44 @@ def determineWavelengthSolution(self, adinputs=None, **params): fit1d.image = image wavecal.update_wcs_with_solution(ext, fit1d, other, config) else: + bad_solutions = [] + single_line_solutions = [] for ext in ad: - if len(ad) > 1: - log.stdinfo(f"Determining solution for extension {ext.id}") - - input_data, fit1d, acceptable_fit = wavecal.get_automated_fit( + input_data, fit1d = wavecal.get_automated_fit( ext, uiparams, p=self, linelist=linelist, bad_bits=DQ.not_signal, absorption=absorption) wavecal.update_wcs_with_solution(ext, fit1d, input_data, config) - if not acceptable_fit: - log.warning("No acceptable wavelength solution found") - else: + if fit1d.image.size: figures.append(wavecal.create_pdf_plot( input_data, fit1d.points[~fit1d.mask], fit1d.image[~fit1d.mask], title=f"{ad.filename}:{ext.id}", absorption=absorption)) + if fit1d.image.size == 1: + single_line_solutions.append(ext.id) + else: # no line matches so not an acceptable solution + bad_solutions.append(ext.id) + + if bad_solutions: + msg = f"{ad.filename}: failed to find an acceptable wavelength solution" + if len(ad) > 1: + if len(ad) > len(bad_solutions): + msg += "in extensions " + ", ".join(str(i) for i in bad_solutions) + else: + msg += " in any extensions" + if self.mode == 'sq' and len(ad) == len(bad_solutions): + raise RuntimeError(msg) + log.warning(msg) + + if single_line_solutions: + msg = f"{ad.filename}: a single-line solution has been applied" + if len(ad) > 1: + if len(ad) > len(single_line_solutions): + msg += "in extensions " + ", ".join(str(i) for i in single_line_solutions) + else: + msg += " in all extensions" + log.warning(msg) + ad.update_filename(suffix=sfx, strip=True) if figures: @@ -3995,6 +4019,63 @@ def maskBeyondSlit(self, adinputs=None, **params): return adinputs + def monitorWavelengthSolution(self, adinputs=None, **params): + """ + This captures some info from the wavelength solution for the GOA + instrument monitoring system into various MWS_ headers for the instrument + monitoring system to pick up. + + Values Captured include the polynomial coefficients and RMS from the + .WAVECAL table, the central wavelength from both the WCS and central_wavlength() + descriptor and the difference between them. + :return: + """ + log = self.log + log.debug(gt.log_message("primitive", self.myself(), "starting")) + + tocapture = ('c0', 'c1', 'c2', 'c3', 'rms', 'inv_rms', 'fwidth') + + for ad in adinputs: + for ext in ad: + wavecal = getattr(ext, 'WAVECAL', None) + if wavecal is None: + continue + + for row in wavecal: + # Capture values from WAVEVAL table into monitoring headers + if row['name'] in tocapture: + keyword = 'MWS_' + row['name'].upper()[:4] + if row['name'] == 'inv_rms': + keyword = 'MWS_IRMS' + ext.hdr[keyword] = (row['coefficients'], + f'WAVECAL {row['name']} coefficient') + + # Get the model and evaluate it at the mean of its domain + model = am.get_named_submodel(ext.wcs.forward_transform, "WAVE") + if model: + # Central Wavelength check + center = np.mean(model.domain) + model_central_wlen = model(center) + # We want actual_central_wavelength rather than central_wavelength + header_central_wlen = ext.actual_central_wavelength(asNanometers=True) + + ext.hdr['MWS_DCWL'] = (model_central_wlen - header_central_wlen, + 'Delta between model and header central_wavelength') + + # Dispersion check + model_dispersion = model(center+0.5) - model(center-0.5) + header_dispersion = ext.dispersion(asNanometers=True) + ext.hdr['MWS_DDIS'] = (model_dispersion - header_dispersion, + 'Delta between model and header dispersion') + + + # Count the number of arc lines in the wavecal table + ext.hdr['MWS_NUML'] = (numpy.count_nonzero(wavecal['peaks']), + 'Number of non-zero peaks in the WAVECAL table') + + ad.update_filename(suffix=params["suffix"], strip=True) + + return adinputs def normalizeFlat(self, adinputs=None, **params): """ diff --git a/geminidr/core/tests/test_spect.py b/geminidr/core/tests/test_spect.py index 53e1278a2..295cee3d7 100644 --- a/geminidr/core/tests/test_spect.py +++ b/geminidr/core/tests/test_spect.py @@ -44,11 +44,12 @@ from specutils.utils.wcs_utils import air_to_vac import astrodata, gemini_instruments -from astrodata.testing import ad_compare, assert_most_close +from astrodata.testing import ad_compare, assert_most_close, download_from_archive from gempy.library import astromodels as am from gempy.library.config.config import FieldValidationError from geminidr.core import primitives_spect from geminidr.f2.primitives_f2_longslit import F2Longslit +from geminidr.gmos.primitives_gmos_longslit import GMOSLongslit from geminidr.niri.primitives_niri_image import NIRIImage from geminidr.niri.primitives_niri_longslit import NIRILongslit from geminidr.gnirs import primitives_gnirs_longslit @@ -574,6 +575,26 @@ def test_adjust_wavelength_zero_point_controlled(filename, center, shift, np.testing.assert_allclose(waves, new_waves, atol=0.1*dw) +@pytest.mark.preprocessed_data +def test_determine_wavelength_solution_exits_with_no_solution_in_sq(path_to_inputs): + """Primitive should crash in SQ mode without a solution""" + ad = astrodata.open(os.path.join(path_to_inputs, "S20260331S0149_mosaic.fits")) + p = GMOSLongslit([ad]) + p.mode = "sq" # it's the default, but just to be sure + with pytest.raises(RuntimeError): + p.determineWavelengthSolution() + + +@pytest.mark.preprocessed_data +def test_determine_wavelength_solution_continues_with_no_solution_in_qa(path_to_inputs, caplog): + """Primitive should crash in SQ mode without a solution""" + ad = astrodata.open(os.path.join(path_to_inputs, "S20260331S0149_mosaic.fits")) + p = GMOSLongslit([ad]) + p.mode = "qa" + p.determineWavelengthSolution() + assert any(record.levelname == 'WARNING' for record in caplog.records) + + @pytest.mark.preprocessed_data @pytest.mark.parametrize('in_file,instrument', [# GNIRS 111/mm LongBlue, off right edge of detector diff --git a/geminidr/gemini/lookups/keyword_comments.py b/geminidr/gemini/lookups/keyword_comments.py index 8d6769abf..34a84ada4 100644 --- a/geminidr/gemini/lookups/keyword_comments.py +++ b/geminidr/gemini/lookups/keyword_comments.py @@ -48,7 +48,7 @@ "OBSTYPE": "Observation type", "ORIGNAME": "Original filename prior to processing", "OVERRMS": "RMS of fit to overscan mean values", - "OVERRDNS": "Readnoise estimate derived from overscan", + "OVERRDNS": "Readnoise estimate from overscan [electron]", "OVERSCAN": "Median of Overscan mean values", "OVERSEC": "Section used for overscan calculation", "PIXSCALE": "Pixel scale [arcsec/pixel]", diff --git a/geminidr/gmos/primitives_gmos_longslit.py b/geminidr/gmos/primitives_gmos_longslit.py index 3b22e68f0..a3ff8723d 100644 --- a/geminidr/gmos/primitives_gmos_longslit.py +++ b/geminidr/gmos/primitives_gmos_longslit.py @@ -242,13 +242,20 @@ def addIllumMaskToDQ(self, adinputs=None, suffix=None, illum_mask=None, plt.show() plt.ion() - yshift = mshift - maxima[0] - if len(maxima) > 1 or abs(yshift) > mshift: - log.warning(f"{ad.filename}: cross-correlation peak is" - " untrustworthy so not adding illumination " - "mask. Please re-run with a specified shift.") + if not maxima: + log.warning(f"{ad.filename}: no peaks found in " + "cross-correlation, so not adding " + "illumination mask. Please re-run with a specified shift.") yshift = None log.stdinfo(slit_location_msg) + else: + yshift = mshift - maxima[0] + if len(maxima) > 1 or abs(yshift) > mshift: + log.warning(f"{ad.filename}: cross-correlation peak is" + " untrustworthy so not adding illumination " + "mask. Please re-run with a specified shift.") + yshift = None + log.stdinfo(slit_location_msg) else: yshift = shift diff --git a/geminidr/gmos/recipes/qa/recipes_ARC_LS_SPECT.py b/geminidr/gmos/recipes/qa/recipes_ARC_LS_SPECT.py index 74f2893d9..6493547e0 100644 --- a/geminidr/gmos/recipes/qa/recipes_ARC_LS_SPECT.py +++ b/geminidr/gmos/recipes/qa/recipes_ARC_LS_SPECT.py @@ -19,5 +19,42 @@ def makeProcessedArc(p): p.storeProcessedArc() p.writeOutputs() +def checkProcessedArc(p): + """ + Extracts some values from the .WAVECAL extension for instrument monitoring + :param p: + :return: + """ + + p.monitorWavelengthSolution() + p.writeOutputs() + +def checkArc(p): + """ + Reduce an arc, extract some values from the .WAVECAL extension for + instrument monitoring + :param p: + :return: + """ + # In QA mode, if determineWavelengthSolution cannot find a good arc fit it + # uses a linear model based on the headers, whereas in SQ mode it raises an + # exception. For the purposes of checking arcs, we want everything to + # behave as if it is in SQ mode, so we just set that here + p.mode = 'sq' + + p.prepare(require_wcs=False) + p.addDQ() + p.addVAR(read_noise=True) + p.overscanCorrect() + p.ADUToElectrons() + p.addVAR(poisson_noise=True) + p.mosaicDetectors() + # We force center=None to ensure it always gets the central row. + # This is the default anyway, but this prevents it being overridden. + p.determineWavelengthSolution(center=None) + # We don't need determineDistortion to check Arcs, and it's slow. + # p.determineDistortion() + p.monitorWavelengthSolution() + p.writeOutputs() _default = makeProcessedArc diff --git a/geminidr/gnirs/tests/crossdispersed/test_determine_wavelength_solution.py b/geminidr/gnirs/tests/crossdispersed/test_determine_wavelength_solution.py index 294f449df..590ab1042 100644 --- a/geminidr/gnirs/tests/crossdispersed/test_determine_wavelength_solution.py +++ b/geminidr/gnirs/tests/crossdispersed/test_determine_wavelength_solution.py @@ -328,9 +328,9 @@ def test_regression_determine_wavelength_solution( np.testing.assert_allclose(wavelength[indices], ref_wavelength[indices], atol=tolerance) print(f"Test passed for {ad.filename}, extension {wcalibrated_ext.id}") - except AssertionError: + except AssertionError as e: failed = True - raise + raise AssertionError(f"Test failed for {ad.filename}, extension {wcalibrated_ext.id}\n{repr(model)}\n{repr(ref_model)}") finally: if write_report: do_report(wcalibrated_ext, ref_ext, failed=failed) diff --git a/geminidr/gnirs/tests/longslit/test_determine_wavelength_solution.py b/geminidr/gnirs/tests/longslit/test_determine_wavelength_solution.py index 1d1faaa42..0c6305960 100644 --- a/geminidr/gnirs/tests/longslit/test_determine_wavelength_solution.py +++ b/geminidr/gnirs/tests/longslit/test_determine_wavelength_solution.py @@ -246,6 +246,7 @@ def test_regression_determine_wavelength_solution( pixel_scale = ad[0].pixel_scale() # arcsec / px p.viewer = geminidr.dormantViewer(p, None) + p.mode = 'qa' p.determineWavelengthSolution(**{**determine_wavelength_solution_parameters, **params}) diff --git a/geminidr/interactive/fit/fit1d.py b/geminidr/interactive/fit/fit1d.py index bac54f01b..52e4b3411 100644 --- a/geminidr/interactive/fit/fit1d.py +++ b/geminidr/interactive/fit/fit1d.py @@ -1238,6 +1238,8 @@ def build_figures( ylabel=self.ylabel, model=self.model, enable_user_masking=self.enable_user_masking, + plot_ratios=plot_ratios, + plot_residuals=plot_residuals, ) if self.enable_regions: @@ -2291,6 +2293,10 @@ def fit1d_figure( tabs = bm.Tabs( tabs=[], sizing_mode="stretch_width", stylesheets=dragons_styles() ) + try: # needed for bokeh 3.9 + tabs.link_layouts = True + except AttributeError: # bokeh < 3.9 + pass tabs.tabs.append(bm.TabPanel(child=p_resid, title="Residuals")) tabs.tabs.append(bm.TabPanel(child=p_ratios, title="Ratios")) diff --git a/geminidr/interactive/fit/tracing.py b/geminidr/interactive/fit/tracing.py index 69e45e505..62a65962d 100644 --- a/geminidr/interactive/fit/tracing.py +++ b/geminidr/interactive/fit/tracing.py @@ -75,6 +75,7 @@ def interactive_trace_apertures(ad, tab_labels, fit1d_params, ui_params: UIParam modal_button_label="Trace apertures", modal_message="Tracing apertures...", ui_params=ui_params, + plot_ratios=False, turbo_tabs=True ) diff --git a/gempy/library/fitting.py b/gempy/library/fitting.py index 9e7d59d1d..df08f83fa 100644 --- a/gempy/library/fitting.py +++ b/gempy/library/fitting.py @@ -613,3 +613,32 @@ def translate_params(params): if "grow" in params: new_params["grow"] = params["grow"] return new_params + + @classmethod + def create_with_model(cls, model, image, points): + """ + Create a fit_1D object with an already-defined model and an array + of points. + + Parameters + ---------- + model : callable + A model that + + image : array-like + N-dimensional input array containing the values to be fitted. If + it is a `numpy.ma.MaskedArray`, any masked points are ignored when + fitting. + + points : `~numpy.ndarray`, optional + 1-dimensional input array with the x-values of each 1D slice being + fitted, If not given, the independent variable will be treated as + a sequence of integers starting at zero. + model + """ + fit1d = cls(np.arange(5), function="chebyshev", order=1, niter=0) + fit1d._models = model + fit1d.image = np.asarray(image) + fit1d.points = np.asarray(points) + fit1d.mask = np.zeros_like(fit1d.image, dtype=bool) + return fit1d diff --git a/gempy/library/wavecal.py b/gempy/library/wavecal.py index e02c65bac..2f169a616 100644 --- a/gempy/library/wavecal.py +++ b/gempy/library/wavecal.py @@ -423,7 +423,7 @@ def create_interactive_inputs(ad, ui_params=None, p=None, linelist=None, bad_bits=0, absorption=False): data = {"x": [], "y": [], "meta": []} for ext in ad: - input_data, fit1d, _ = get_automated_fit( + input_data, fit1d = get_automated_fit( ext, ui_params, p=p, linelist=linelist, bad_bits=bad_bits, absorption=absorption) # peak locations and line wavelengths of matched peaks/lines @@ -486,8 +486,6 @@ class holding parameters for the UI, passed from the primitive's Config fit1d : a fit_1D object containing the wavelength solution, plus an "image" attribute that lists the matched arc line wavelengths - acceptable_fit : bool - whether this fit is likely to be good """ input_data = get_all_input_data( ext, p, ui_params.toDict(), linelist=linelist, bad_bits=bad_bits, @@ -499,14 +497,14 @@ class holding parameters for the UI, passed from the primitive's Config dw = np.diff(init_models[0](np.arange(spectrum.size))).mean() kdsigma = fwidth * abs(dw) k = 1 if kdsigma < 3 else 2 - fit1d, acceptable_fit = find_solution( + fit1d = find_solution( init_models, ui_params.toDict(), peaks=peaks, peak_weights=weights[ui_params.weighting], linelist=input_data["linelist"], fwidth=fwidth, kdsigma=kdsigma, k=k, bounds_setter=input_data["bounds_setter"], filename=ext.filename) input_data["fit"] = fit1d - return input_data, fit1d, acceptable_fit + return input_data, fit1d def create_chebyshev(waves, central_wavelength=None, dispersion=None, @@ -773,9 +771,9 @@ def find_solution(init_models, config, peaks=None, peak_weights=None, Returns ------- - length-2 tuple - A tuple of the best-fit model and a boolean denoting whether or not it - is an acceptable fit. + fit_1D: + the best-fit model or the original model, of no acceptable fit is + found. The "image" attribute will be empty if no fit is found. """ log = logutils.get_logger(__name__) min_lines = [int(x) for x in str(config["debug_min_lines"]).split(',')] @@ -805,57 +803,91 @@ def find_solution(init_models, config, peaks=None, peak_weights=None, min_lines_per_fit=min_lines_per_fit, k=k, bounds_setter=bounds_setter) + wavemin, wavemax = sorted(model(model.domain)) + nlines_expected = sum(wavemin <= wline <= wavemax for wline in arc_lines) + # We perform a regular least-squares fit to all the matches # we've made. This allows a high polynomial order to be # used without the risk of it going off the rails matched = np.where(matches > -1)[0] fit_it = fitting.LinearLSQFitter() - if len(matched) > 1: # need at least 2 lines, right? - m_init = models.Chebyshev1D(degree=min(config["order"], len(matched)-1), - domain=domain) + + # We allow this continue if we only matched 1 line but only 1 line is + # expected + if len(matched) > 1 or (len(matched) * nlines_expected == 1 and len(peaks) <= 2): + m_init = models.Chebyshev1D( + degree=min(config["order"], (len(matched) + 1) // 2), + domain=domain) for p, v in zip(model.param_names, model.parameters): if p in m_init.param_names: setattr(m_init, p, v) - #bounds_setter(m_init) + model_bounds = bounds_setter(m_init) + if len(matched) == 1: + m_init.c1.fixed = True #for i in range(len(matched), m_init.degree + 1): # m_init.fixed[f"c{i}"] = True matched_peaks = peaks[matched] matched_arc_lines = arc_lines[matches[matched]] m_final = fit_it(m_init, matched_peaks, matched_arc_lines) + dw = abs(np.diff(m_final(m_final.domain))[0] / np.diff(m_final.domain)[0]) #for p, l in zip(matched_peaks, matched_arc_lines): # print(f"{p:.2f} => {l:.2f}") - # We're close to the correct solution, perform a KDFit - m_init = models.Chebyshev1D(degree=config["order"], domain=domain) - for p, v in zip(m_final.param_names, m_final.parameters): - setattr(m_init, p, v) - dw = abs(np.diff(m_final(m_final.domain))[0] / np.diff(m_final.domain)[0]) - fit_it = matching.KDTreeFitter(sigma=2 * abs(dw), maxsig=5, - k=k, method='Nelder-Mead') - m_final = fit_it(m_init, peaks, arc_lines, in_weights=peak_weights, - ref_weights=arc_weights) - logit(f'{repr(m_final)} {fit_it.statistic}') - - # And then recalculate the matches - match_radius = 4 * fwidth * abs(m_final.c1) / len_data # 2*fwidth pixels - try: - matched = matching.match_sources(m_final(peaks), arc_lines, - radius=match_radius) - incoords, outcoords = zip(*[(peaks[i], arc_lines[m]) - for i, m in enumerate(matched) if m > -1]) - # Probably incoords and outcoords as defined here should go to - # the interactive fitter, but cull to derive the "best" model - fit1d = fit_1D(outcoords, points=incoords, function="chebyshev", - order=min(m_final.degree, len(incoords)-1), - domain=m_final.domain, - niter=config["niter"], sigma_lower=config["lsigma"], - sigma_upper=config["hsigma"]) - fit1d.image = np.asarray(outcoords) - except ValueError: - log.warning("Line-matching failed") - continue - nmatched = np.sum(~fit1d.mask) - logit(f"{filename} {repr(fit1d.model)} {nmatched} {fit1d.rms}") + if len(matched) > 1: + # We're close to the correct solution, perform a KDFit + m_init = models.Chebyshev1D(m_final.degree, domain=domain) + for p, v in zip(m_final.param_names, m_final.parameters): + setattr(m_init, p, v) + fit_it = matching.KDTreeFitter(sigma=2 * abs(dw), maxsig=5, + k=k, method='Nelder-Mead') + m_final = fit_it(m_init, peaks, arc_lines, in_weights=peak_weights, + ref_weights=arc_weights) + logit(f'{repr(m_final)} {fit_it.statistic}') + + # And then recalculate the matches + match_radius = 4 * fwidth * abs(m_final.c1) / len_data # 2*fwidth pixels + try: + matched = matching.match_sources(m_final(peaks), arc_lines, + radius=match_radius) + incoords, outcoords = zip(*[(peaks[i], arc_lines[m]) + for i, m in enumerate(matched) if m > -1]) + # Probably incoords and outcoords as defined here should go to + # the interactive fitter, but cull to derive the "best" model + fit1d = fit_1D(outcoords, points=incoords, function="chebyshev", + order=min(m_final.degree, (len(incoords) + 1) // 2), + domain=m_final.domain, + niter=config["niter"], sigma_lower=config["lsigma"], + sigma_upper=config["hsigma"]) + fit1d.image = np.asarray(outcoords) + except ValueError: + log.warning("Line-matching failed") + continue + nmatched = np.sum(~fit1d.mask) + logit(f"{filename} {repr(fit1d.model)} {nmatched} {fit1d.rms}") + + # A posteriori check that the fit is within the bounds + # (LinearLSQFitter does not allow bounds so we can't force + # the model to be bounded). We allow the tolerances to be twice + # as large as they were during the piecewise fit, as this seems + # to work well empirically. + try: + for p, v in zip(fit1d.model.param_names, fit1d.model.parameters): + bounds = model_bounds.get(p, (-np.inf, np.inf)) + new_bounds = (np.asarray(bounds) + + np.array([-0.5, 0.5]) * np.diff(bounds)[0]) + assert new_bounds[0] <= v <= new_bounds[1] + except AssertionError: + for p, v in zip(fit1d.model.param_names, fit1d.model.parameters): + bounds = model_bounds.get(p, (-np.inf, np.inf)) + new_bounds = (np.asarray(bounds) + + np.array([-0.5, 0.5]) * np.diff(bounds)[0]) + print(p, v, new_bounds) + continue + else: + # Hack a fit1D object with the model and the single match + fit1d = fit_1D.create_with_model( + m_final, image=matched_arc_lines, points=matched_peaks) + nmatched = 1 # Wavelength solution models need to be monotonic. Make that check. waves = fit1d.evaluate(np.arange(len_data)) @@ -870,7 +902,7 @@ def find_solution(init_models, config, peaks=None, peak_weights=None, # Trial and error suggests this criterion works well if fit1d.rms < 0.8 / config["order"] * fwidth * abs(dw) and nmatched >= min_matches_required: #print("RETURNING", fit1d.model.parameters) - return fit1d, True + return fit1d # This seems to be a reasonably ranking for poor models if nmatched > config["order"] + 1: @@ -883,13 +915,8 @@ def find_solution(init_models, config, peaks=None, peak_weights=None, if best_fit1d is None: # Hack a fit1D object that represents the original model with no fitted lines - best_fit1d = fit_1D(np.arange(5), function="chebyshev", order=1, - niter=0) - best_fit1d._models = init_models[0] - best_fit1d.image = np.array([]) - best_fit1d.points = np.array([]) - best_fit1d.mask = np.array([], dtype=bool) - return best_fit1d, True + best_fit1d = fit_1D.create_with_model(init_models[0], [], []) + return best_fit1d def perform_piecewise_fit(model, peaks, arc_lines, pixel_start, kdsigma, @@ -1168,6 +1195,8 @@ def update_wcs_with_solution(ext, fit1d, input_data, config): if len(incoords) > 1: inv_rms = np.std(m_inverse(m_final(incoords)) - incoords) log.stdinfo(f"Inverse model has rms = {inv_rms:.3f} pixels.") + else: + inv_rms = np.nan m_final.name = "WAVE" # always WAVE, never AWAV m_final.inverse = m_inverse @@ -1182,8 +1211,8 @@ def update_wcs_with_solution(ext, fit1d, input_data, config): # while I refactor tests temptable.add_columns([[1], [m_final.degree], [domain[0]], [domain[1]]], names=("ndim", "degree", "domain_start", "domain_end")) - temptable.add_columns([[rms], [input_data["fwidth"]]], - names=("rms", "fwidth")) + temptable.add_columns([[rms], [inv_rms], [input_data["fwidth"]]], + names=("rms", "inv_rms", "fwidth")) if ext.data.ndim > 1: # TODO: Need to update this from the interactive tool's values direction, location = input_data["location"].split() @@ -1276,14 +1305,15 @@ def create_pdf_plot(input_data, peaks, arc_lines, title="", data = input_data["spectrum"] spacing = 0.01 vert_align = "bottom" - xmin, xmax = input_data["init_models"][0].domain + init_model = input_data["init_models"][0] + xmin, xmax = init_model.domain pixels = np.arange(xmin, xmax + 1) data_max = data.max() spacing *= data_max fig, ax = plt.subplots() ax.plot(pixels, data, 'b-') ax.set_ylim(0, data_max * 1.1) - if len(arc_lines) and np.diff(arc_lines)[0] / np.diff(peaks)[0] < 0: + if len(arc_lines) and init_model(xmin) > init_model(xmax): ax.set_xlim(xmax + 1, xmin - 1) else: ax.set_xlim(xmin - 1, xmax + 1) diff --git a/recipe_system/cal_service/file_getter.py b/recipe_system/cal_service/file_getter.py index 039065e2d..9bed303b9 100644 --- a/recipe_system/cal_service/file_getter.py +++ b/recipe_system/cal_service/file_getter.py @@ -125,11 +125,15 @@ def get_request(self, url, filename, calmd5): tuples = [] total_bytes = 0 for dirent in dirents: - statobj = dirent.stat() # To avoid multiple stat calls - tuples.append((statobj.st_atime, - dirent.path, - statobj.st_size)) - total_bytes += statobj.st_size + try: + statobj = dirent.stat() # To avoid multiple stat calls + tuples.append((statobj.st_atime, + dirent.path, + statobj.st_size)) + total_bytes += statobj.st_size + except FileNotFoundError: + # Race condition where another process deleted the file + pass # Sort list by increasing atime (ie oldest first) tuples.sort() while total_bytes > (1E9 * self.cachegbs):