From 868a99f129a5f218cd71abffe35f0cefc36719fe Mon Sep 17 00:00:00 2001 From: astrofle Date: Thu, 16 Jan 2025 14:50:29 -0500 Subject: [PATCH 1/8] Fix: multiple baseline exclusion regions can use units Leverages SpectralRegion with some tweaks to account for decreasing spectral axis --- src/dysh/spectra/core.py | 106 ++++++++++++++------------------------- 1 file changed, 39 insertions(+), 67 deletions(-) diff --git a/src/dysh/spectra/core.py b/src/dysh/spectra/core.py index 6cef9a86..25ee7cec 100644 --- a/src/dysh/spectra/core.py +++ b/src/dysh/spectra/core.py @@ -158,7 +158,7 @@ def find_nonblank_ints(cycle1, cycle2, cycle3=None, cycle4=None): return goodrows -def exclude_to_region(exclude, refspec, fix_exclude=False): +def exclude_to_region(exclude, refspec, fix_exclude=True): """Convert an exclude list to a list of ~specutuls.SpectralRegion. Parameters @@ -188,81 +188,53 @@ def exclude_to_region(exclude, refspec, fix_exclude=False): fix_exclude: bool If True, fix exclude regions that are out of bounds of the specctral axis to be within the spectral axis. Default:False - Returns - ------- - regionlist : list of `~specutil.SpectralRegion` - A list of `~specutil.SpectralRegion` corresponding to `exclude` with units of the `refspec.spectral_axis`. + Returns + ------- + regionlist : list of `~specutil.SpectralRegion` + A list of `~specutil.SpectralRegion` corresponding to `exclude` with units of the `refspec.spectral_axis`. """ + regionlist = [] p = refspec sa = refspec.spectral_axis + lastchan = len(sa) - 1 + + o = 1 + # If the spectral axis is inverted, flip the order of the elements. + if refspec.spectral_axis_direction == "decreasing": + o = -1 + if exclude is not None: - regionlist = [] - # a single SpectralRegion was given + # A single SpectralRegion was given. if isinstance(exclude, SpectralRegion): b = exclude.bounds - if b[0] < sa[0] or b[1] > sa[1]: - msg = f"Exclude limits {b} are not fully within the spectral axis {sa}" - raise Exception(msg) regionlist.append(exclude) - # list of int or Quantity or SpectralRegion was given + # list of int or Quantity or SpectralRegion was given. else: - # if user provided a single list, we have to - # add another set of brackets so we an iterate. - # If SpectralRegion took a list argument, we wouldn't - # have to do this. - if len(np.shape(exclude[0])) == 0: - exclude = [exclude] - # NB: we are assuming that a SpectralAxis is always [lower...upper]. Is this true??? - for pair in exclude: - if type(pair[0]) == int: - # convert channel to spectral axis units - lastchan = len(sa) - 1 - msg = f"Exclude limits {pair} are not fully within the spectral axis [0,{lastchan}]." - if pair[0] < 0 or pair[1] > lastchan: - if fix_exclude: - msg += f" Setting upper limit to {lastchan}." - pair[1] = lastchan - warnings.warn(msg) - else: - raise Exception(msg) - pair = [sa[pair[0]], sa[pair[1]]] - # if it is already a spectral region no additional - # work is needed - # @todo we should test that the SpectralRegion is not out of bounds - if isinstance(pair[0], SpectralRegion): - b = pair[0].bounds - if b[0] < sa[0] or b[1] > sa[1]: - msg = f"Exclude limits {pair} are not fully within the spectral axis {p.spectral_axis}" - raise Exception(msg) - regionlist.append(pair) - else: # it is a Quantity that may need conversion to spectral_axis units - q = [pair[0].value, pair[1].value] * pair[0].unit - if q.unit.is_equivalent("km/s"): - veldef = p.meta.get("VELDEF", None) - if veldef is None: - raise KeyError("Input spectrum has no VELDEF in header, can't convert to frequency units.") - pair = veltofreq(q, p.rest_value, veldef) - # offset = p.rest_value - p.radial_velocity.to(sa.unit, equivalencies=p.equivalencies) - else: - pair[0] = pair[0].to(sa.unit, equivalencies=p.equivalencies) - pair[1] = pair[1].to(sa.unit, equivalencies=p.equivalencies) - # Ensure test is with sorted [lower,upper] - pair = sorted(pair) - salimits = sorted([sa[0], sa[-1]]) - if pair[0] < salimits[0] or pair[1] > salimits[-1]: - msg = ( - f"Exclude limits {pair} are not fully within the spectral axis" - f" {[salimits[0],salimits[-1]]}." - ) - if fix_exclude: - msg += f" Setting upper limit to {p.spectral_axis[-1]}." - pair[1] = sa[-1] - warnings.warn(msg) - else: - raise Exception(msg) - sr = SpectralRegion(pair[0], pair[1]) - regionlist.append(sr) + # If user provided a single list, we have to + # turn it into a list of tuples. If SpectralRegion + # took a list argument, we wouldn't have to do this. + if type(exclude[0]) is not tuple: + it = iter(exclude) + exclude = list(zip(it, it)) + try: + sr = SpectralRegion(exclude) + # The above will error if the elements are not quantities. + # In that case use the spectral axis to define the exclusion regions. + except ValueError: + # If the user requested to fix the exclude range. + if fix_exclude: + exclude = np.array(exclude) + mask = exclude > len(sa) + if mask.sum() > 0: + msg = f"Setting upper limit to {lastchan}." + exclude[exclude > len(sa)] = lastchan + warnings.warn(msg) + # If the spectral_axis is decreasing, flip it. + sr = SpectralRegion(sa[exclude][:, ::o]) + # The continuum fitting routines use a list of SpectralRegions as input. + for r in sr.subregions: + regionlist.append(SpectralRegion([r])) return regionlist From 0de7a0c9f6f10e354eada781bb4b00ed923fcd7b Mon Sep 17 00:00:00 2001 From: astrofle Date: Thu, 16 Jan 2025 15:02:33 -0500 Subject: [PATCH 2/8] Add: test for baseline subtraction including inclusion regions --- src/dysh/spectra/tests/test_spectrum.py | 30 +++++++++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/src/dysh/spectra/tests/test_spectrum.py b/src/dysh/spectra/tests/test_spectrum.py index dc1f5b9b..20295fa1 100644 --- a/src/dysh/spectra/tests/test_spectrum.py +++ b/src/dysh/spectra/tests/test_spectrum.py @@ -631,3 +631,33 @@ def test_single_baseline(sdf, ex_reg, order, model, gbtidl_bmodel, exclude_regio ex_reg = None test_single_baseline(sdf, ex_reg, order, "chebyshev", gbtidl_no_reg) + + def test_baseline_include(self): + """ + Test for comparing GBTIDL baseline to dysh baselines + using inclusion regions. + """ + + def test_single_baseline(sdf, in_reg, order, model, gbtidl_bmodel, exclude_region_upper_bounds=True): + """For use with TestSpectrum.test_baseline()""" + dysh_spec = sdf.getspec(0) + temp_bmodel = np.copy(dysh_spec.data) + dysh_spec.baseline( + order, include=in_reg, remove=True, model=model, exclude_region_upper_bounds=exclude_region_upper_bounds + ) + dysh_bmodel = temp_bmodel - np.copy(dysh_spec.data) + diff = np.sum(np.abs(dysh_bmodel - gbtidl_bmodel)) + assert diff < 1.5e-6 + + data_dir = get_project_testdata() / "AGBT17A_404_01" + sdf_file = data_dir / "AGBT17A_404_01_scan_19_prebaseline.fits" + sdf = GBTFITSLoad(sdf_file) + gbtidl_two_reg = loadfits(data_dir / "AGBT17A_404_01_scan_19_bmodel.fits") + gbtidl_no_reg = loadfits(data_dir / "AGBT17A_404_01_scan_19_noregion_bmodel.fits") + + order = 3 + in_reg = [(99, 381), (449, 721)] + + test_single_baseline(sdf, in_reg, order, "chebyshev", gbtidl_two_reg) + test_single_baseline(sdf, in_reg, order, "legendre", gbtidl_two_reg) + test_single_baseline(sdf, in_reg, order, "hermite", gbtidl_two_reg) From 59f80210dd2592ef9ca285212ab3078b7d074dca Mon Sep 17 00:00:00 2001 From: astrofle Date: Thu, 16 Jan 2025 16:13:52 -0500 Subject: [PATCH 3/8] Fix: return support for other units --- src/dysh/spectra/core.py | 21 ++++++++++++++------- 1 file changed, 14 insertions(+), 7 deletions(-) diff --git a/src/dysh/spectra/core.py b/src/dysh/spectra/core.py index 25ee7cec..47d0e532 100644 --- a/src/dysh/spectra/core.py +++ b/src/dysh/spectra/core.py @@ -207,13 +207,13 @@ def exclude_to_region(exclude, refspec, fix_exclude=True): if exclude is not None: # A single SpectralRegion was given. if isinstance(exclude, SpectralRegion): - b = exclude.bounds - regionlist.append(exclude) - # list of int or Quantity or SpectralRegion was given. + sr = exclude + # `list` of `int` or `Quantity` or `SpectralRegion` was given. else: # If user provided a single list, we have to # turn it into a list of tuples. If SpectralRegion # took a list argument, we wouldn't have to do this. + print(exclude) if type(exclude[0]) is not tuple: it = iter(exclude) exclude = list(zip(it, it)) @@ -232,11 +232,18 @@ def exclude_to_region(exclude, refspec, fix_exclude=True): warnings.warn(msg) # If the spectral_axis is decreasing, flip it. sr = SpectralRegion(sa[exclude][:, ::o]) - # The continuum fitting routines use a list of SpectralRegions as input. - for r in sr.subregions: - regionlist.append(SpectralRegion([r])) - return regionlist + # Change units to match those of the `Spectrum`. + qt = sr.as_table() + lb = qt["lower_bound"].to(p._spectral_axis.unit, equivalencies=p.equivalencies) + ub = qt["upper_bound"].to(p._spectral_axis.unit, equivalencies=p.equivalencies) + lsr = list(zip(lb, ub)) + + # The continuum fitting routines use a list of `SpectralRegions` as input. + for r in lsr: + regionlist.append(SpectralRegion([r])) + + return regionlist def region_to_axis_indices(region, refspec): From cf36d1bb58a9a437bbaa7f4b7d63ea4a301fa627 Mon Sep 17 00:00:00 2001 From: astrofle Date: Fri, 17 Jan 2025 13:55:14 -0500 Subject: [PATCH 4/8] Add: support for units in include regions This expands include, for baseline fitting, to use units. --- src/dysh/spectra/core.py | 236 ++++++++++++++++++++++++++++++----- src/dysh/spectra/spectrum.py | 3 +- 2 files changed, 207 insertions(+), 32 deletions(-) diff --git a/src/dysh/spectra/core.py b/src/dysh/spectra/core.py index 47d0e532..48e63297 100644 --- a/src/dysh/spectra/core.py +++ b/src/dysh/spectra/core.py @@ -158,29 +158,108 @@ def find_nonblank_ints(cycle1, cycle2, cycle3=None, cycle4=None): return goodrows -def exclude_to_region(exclude, refspec, fix_exclude=True): +def sort_spectral_region(spectral_region): + """ + Sort the elements of a `~specutils.SpectralRegion`. + + Parameters + ---------- + spectral_region : `~specutils.SpectralRegion` + `~specutils.SpectralRegion` to be sorted. + + Returns + ------- + sorted_spectral_region : `~specutils.SpectralRegion` + Sorted `~specutils.SpectralRegion`. + """ + + unit = spectral_region.lower.unit + bound_list = np.sort([srb.value for sr in spectral_region.subregions for srb in sr]) * unit + it = iter(bound_list) + sorted_spectral_region = SpectralRegion(list(zip(it, it))) + + return sorted_spectral_region + + +def include_to_exclude_spectral_region(include, refspec): + """ + Convert an inclusion region to an exclude region. + + Parameters + ---------- + include : list of 2-tuples of int or `~astropy.units.Quantity`, or `~specutils.SpectralRegion` + List of region(s) to exclude from the fit. The tuple(s) represent a range in the form [lower,upper], inclusive. + Examples: + + One channel-based region: + + >>> [11,51] + + Two channel-based regions: + + >>> [(11,51),(99,123)] + + One `~astropy.units.Quantity` region: + + >>> [110.198*u.GHz,110.204*u.GHz]. + + One compound `~specutils.SpectralRegion`: + + >>> SpectralRegion([(110.198*u.GHz,110.204*u.GHz),(110.196*u.GHz,110.197*u.GHz)]). + + refspec: `~spectra.spectrum.Spectrum` + The reference spectrum whose spectral axis will be used + when converting between include and axis units (e.g. channels to GHz). + + Returns + ------- + exclude_region : `~specutils.SpectralRegion` + `include` as a region to be excluded. + """ + + unit = "cm" # This has to be wavelength units. cm are familiar to observers? + + spectral_region = exclude_to_spectral_region(include, refspec) + spectral_region = spectral_region_to_unit(spectral_region, refspec, unit=unit) + # Convert unit of reference spectrum to wavelength. + # For some reason `SpectralRegion.invert_from_spectrum` + # only works in wavelength units. + org_unit = refspec._spectral_axis.unit + refspec._spectral_axis = refspec._spectral_axis.to(unit, equivalencies=refspec.equivalencies) + exclude_region = spectral_region.invert_from_spectrum(refspec) + refspec._spectral_axis = refspec._spectral_axis.to(org_unit, equivalencies=refspec.equivalencies) + exclude_region = spectral_region_to_unit(exclude_region, refspec) + # Sort the output. + # For some reason it is scrambled. + exclude_region = sort_spectral_region(exclude_region) + + return exclude_region + + +def exclude_to_spectral_region(exclude, refspec, fix_exclude=True): """Convert an exclude list to a list of ~specutuls.SpectralRegion. Parameters ---------- exclude : list of 2-tuples of int or `~astropy.units.Quantity`, or `~specutils.SpectralRegion` - List of region(s) to exclude from the fit. The tuple(s) represent a range in the form [lower,upper], inclusive. Examples: + List of region(s) to exclude from the fit. The tuple(s) represent a range in the form [lower,upper], inclusive. + Examples: - One channel-based region: + One channel-based region: - >>> [11,51] + >>> [11,51] - Two channel-based regions: + Two channel-based regions: - >>> [(11,51),(99,123)] + >>> [(11,51),(99,123)] - One `~astropy.units.Quantity` region: + One `~astropy.units.Quantity` region: - >>> [110.198*u.GHz,110.204*u.GHz]. + >>> [110.198*u.GHz,110.204*u.GHz]. - One compound `~specutils.SpectralRegion`: + One compound `~specutils.SpectralRegion`: - >>> SpectralRegion([(110.198*u.GHz,110.204*u.GHz),(110.196*u.GHz,110.197*u.GHz)]). + >>> SpectralRegion([(110.198*u.GHz,110.204*u.GHz),(110.196*u.GHz,110.197*u.GHz)]). refspec: `~spectra.spectrum.Spectrum` The reference spectrum whose spectral axis will be used @@ -190,11 +269,10 @@ def exclude_to_region(exclude, refspec, fix_exclude=True): Returns ------- - regionlist : list of `~specutil.SpectralRegion` - A list of `~specutil.SpectralRegion` corresponding to `exclude` with units of the `refspec.spectral_axis`. + sr : `~specutil.SpectralRegion` + A `~specutil.SpectralRegion` corresponding to `exclude`. """ - regionlist = [] p = refspec sa = refspec.spectral_axis lastchan = len(sa) - 1 @@ -213,7 +291,6 @@ def exclude_to_region(exclude, refspec, fix_exclude=True): # If user provided a single list, we have to # turn it into a list of tuples. If SpectralRegion # took a list argument, we wouldn't have to do this. - print(exclude) if type(exclude[0]) is not tuple: it = iter(exclude) exclude = list(zip(it, it)) @@ -233,17 +310,64 @@ def exclude_to_region(exclude, refspec, fix_exclude=True): # If the spectral_axis is decreasing, flip it. sr = SpectralRegion(sa[exclude][:, ::o]) - # Change units to match those of the `Spectrum`. - qt = sr.as_table() - lb = qt["lower_bound"].to(p._spectral_axis.unit, equivalencies=p.equivalencies) - ub = qt["upper_bound"].to(p._spectral_axis.unit, equivalencies=p.equivalencies) - lsr = list(zip(lb, ub)) + return sr + + +def spectral_region_to_unit(spectral_region, refspec, unit=None): + """ + Change the unit of `spectral_region` to `unit` using the equivalencies of `refspec`. + If no `unit` is provided, it will change to the units of `refspec._spectral_axis`. + + Parameters + ---------- + spectral_region : `~specutil.SpectralRegion` + `~specutil.SpectralRegion` whose units will be converted. + refspec : `~spectra.spectrum.Spectrum` + The reference spectrum whose spectral axis will be used + when converting to `unit` (e.g. channels to GHz). + unit : str or `~astropy.units.Quantity` + The target units for `spectral_region`. + + Returns + ------- + spectral_region : `~specutil.SpectralRegion` + SpectralRegion with units of `unit`. + """ + + qt = spectral_region.as_table() + + if unit is None: + unit = refspec._spectral_axis.unit + + lb = qt["lower_bound"].to(unit, equivalencies=refspec.equivalencies) + ub = qt["upper_bound"].to(unit, equivalencies=refspec.equivalencies) + + return SpectralRegion(list(zip(lb, ub))) + + +def spectral_region_to_list(spectral_region): + """ + Turn `spectral_region` into a list of `~specutil.SpectralRegion`. + Each subregion in `spectral_region` will be a list element. - # The continuum fitting routines use a list of `SpectralRegions` as input. - for r in lsr: - regionlist.append(SpectralRegion([r])) + Parameters + ---------- + spectral_region : `~specutil.SpectralRegion` + `~specutil.SpectralRegion` to convert into a list of `~specutil.SpectralRegion`. + + Returns + ------- + region_list : list of `~specutil.SpectralRegion` + Subregions of `spectral_region` in a list of `~specutil.SpectralRegion`. + """ + + region_list = [] - return regionlist + # The continuum fitting routines use a list of `SpectralRegion` as input. + for r in spectral_region.subregions: + region_list.append(SpectralRegion([r])) + + return region_list def region_to_axis_indices(region, refspec): @@ -277,18 +401,65 @@ def exclude_to_mask(exclude, refspec): pass +def exclude_to_region_list(exclude, spectrum, fix_exclude=True): + """ + Convert an exclusion region, `exclude`, to a list of `~SpectralRegion`. + This is used for baseline fitting. + + Parameters + ---------- + exclude : list of 2-tuples of int or `~astropy.units.Quantity`, or `~specutils.SpectralRegion` + List of region(s) to exclude from the fit. The tuple(s) represent a range in the form [lower,upper], inclusive. + Examples: + + One channel-based region: + + >>> [11,51] + + Two channel-based regions: + + >>> [(11,51),(99,123)] + + One `~astropy.units.Quantity` region: + + >>> [110.198*u.GHz,110.204*u.GHz]. + + One compound `~specutils.SpectralRegion`: + + >>> SpectralRegion([(110.198*u.GHz,110.204*u.GHz),(110.196*u.GHz,110.197*u.GHz)]). + + spectrum : `~spectra.spectrum.Spectrum` + The reference spectrum whose spectral axis will be used + when converting between exclude and axis units (e.g. channels to GHz). + fix_exclude : bool + See `~spectra.core.exclude_to_spectral_region` for details. Default: True + + Returns + ------- + region_list : list of `~specutil.SpectralRegion` + Regions defined in `exclude` as a list of `~specutil.SpectralRegion`. + """ + + spectral_region = exclude_to_spectral_region(exclude, spectrum, fix_exclude=fix_exclude) + spectral_region = spectral_region_to_unit(spectral_region, spectrum) + region_list = spectral_region_to_list(spectral_region) + + return region_list + + @log_function_call() def baseline(spectrum, order, exclude=None, exclude_region_upper_bounds=True, **kwargs): - """Fit a baseline for a spectrum + """Fit a baseline to `spectrum`. Parameters ---------- spectrum : `~spectra.spectrum.Spectrum` - The input spectrum + The input spectrum. order : int - The order of the polynomial series, a.k.a. baseline order + The order of the polynomial series, a.k.a. baseline order. exclude : list of 2-tuples of int or `~astropy.units.Quantity`, or `~specutils.SpectralRegion` - List of region(s) to exclude from the fit. The tuple(s) represent a range in the form [lower,upper], inclusive. Examples: + List of region(s) to exclude from the fit. The tuple(s) represent a range in the form [lower,upper], inclusive. + Examples: One channel-based region: @@ -305,14 +476,17 @@ def baseline(spectrum, order, exclude=None, exclude_region_upper_bounds=True, ** One compound `~specutils.SpectralRegion`: >>> SpectralRegion([(110.198*u.GHz,110.204*u.GHz),(110.196*u.GHz,110.197*u.GHz)]). + Default: no exclude region model : str One of 'polynomial' or 'chebyshev', Default: 'polynomial' fitter : `~astropy.fitting._FitterMeta` - The fitter to use. Default: `~astropy.fitter.LinearLSQFitter` (with `calc_uncertaintes=True`). Be care when choosing a different fitter to be sure it is optimized for this problem. + The fitter to use. Default: `~astropy.fitter.LinearLSQFitter` (with `calc_uncertaintes=True`). + Be careful when choosing a different fitter to be sure it is optimized for this problem. exclude_region_upper_bounds : bool - Makes the upper bound of any excision region(s) inclusive. Allows excising channel 0 for lower-sideband data, and the last channel for upper-sideband data. + Makes the upper bound of any excision region(s) inclusive. + Allows excising channel 0 for lower-sideband data, and the last channel for upper-sideband data. Returns ------- @@ -325,7 +499,7 @@ def baseline(spectrum, order, exclude=None, exclude_region_upper_bounds=True, ** #'show': False, "model": "chebyshev", "fitter": LinearLSQFitter(calc_uncertainties=True), - "fix_exclude": False, + "fix_exclude": True, "exclude_action": "replace", # {'replace','append', None} } kwargs_opts.update(kwargs) @@ -353,7 +527,7 @@ def baseline(spectrum, order, exclude=None, exclude_region_upper_bounds=True, ** # @todo handle masks return None # or raise exception if exclude is not None: - regionlist = exclude_to_region(exclude, spectrum, fix_exclude=kwargs_opts["fix_exclude"]) + regionlist = exclude_to_region_list(exclude, spectrum, fix_exclude=kwargs_opts["fix_exclude"]) if kwargs_opts["exclude_action"] == "replace": p._exclude_regions = regionlist elif kwargs_opts["exclude_action"] == "append": diff --git a/src/dysh/spectra/spectrum.py b/src/dysh/spectra/spectrum.py index 6c573d99..ad350d06 100644 --- a/src/dysh/spectra/spectrum.py +++ b/src/dysh/spectra/spectrum.py @@ -256,7 +256,8 @@ def baseline(self, degree, exclude=None, include=None, **kwargs): if exclude != None: print(f"Warning: ignoring exclude={exclude}") nchan = len(self._spectral_axis) - exclude = self._toggle_sections(nchan, include) + exclude = core.include_to_exclude_spectral_region(include, self) + # exclude = self._toggle_sections(nchan, include) self._baseline_model = baseline(self, degree, exclude, **kwargs) From 692e77ad1be2031880ff80e5b65dd4afbb88ff8d Mon Sep 17 00:00:00 2001 From: astrofle Date: Fri, 17 Jan 2025 14:22:59 -0500 Subject: [PATCH 5/8] Cleanup: update baseline docstring and remove deprecated _toggle_selection method --- src/dysh/spectra/spectrum.py | 56 ++++++++++++++++-------------------- 1 file changed, 24 insertions(+), 32 deletions(-) diff --git a/src/dysh/spectra/spectrum.py b/src/dysh/spectra/spectrum.py index ad350d06..b91c7e79 100644 --- a/src/dysh/spectra/spectrum.py +++ b/src/dysh/spectra/spectrum.py @@ -144,28 +144,6 @@ def exclude_regions(self): """The baseline exclusion region(s) of this spectrum""" return self._exclude_regions - def _toggle_sections(self, nchan, s): - """helper routine to toggle between an include= and exclude= - only works in channel (0..nchan-1) units - sections s need to be a list of (start_chan,end_chan) tuples, - for example [(100,200),(500,600)] would be an include= - An exclude= needs to start with 0 - channels need to be ordered low to high, but there is no check - for this yet! - """ - ns = len(s) - s1 = [] - e = 0 # set this to 1 if you want to be exact complementary - if s[0][0] == 0: - for i in range(ns - 1): - s1.append((s[i][1] + e, s[i + 1][0] - e)) - else: - s1.append((0, s[0][0])) - for i in range(ns - 1): - s1.append((s[i][1], s[i + 1][0])) - s1.append((s[ns - 1][1], nchan - 1)) - return s1 - ##@todo # def exclude_region(self,region): # where region is SpectralRegion, channels, velocity, etc. See core.py baseline method. @@ -182,23 +160,39 @@ def baseline_model(self): def baseline(self, degree, exclude=None, include=None, **kwargs): # fmt: off """ - Compute and optionally remove a baseline. The model for the + Compute and optionally remove a baseline. The model for the baseline can be either a `1D polynomial model `_ or a `1D Chebyshev polynomial of the first kind `_. - The code uses `astropy.modeling` - and `astropy.fitter` to compute the baseline. See the documentation for those modules. + The code uses `~astropy.modeling` + and `~astropy.fitter` to compute the baseline. See the documentation for those modules. This method will set the `baseline_model` attribute to the fitted model function which can be evaluated over a domain. - Note that include= and exclude= are mutually exclusive. + Note that `include` and `exclude` are mutually exclusive. If both are present, only `include` will be used. Parameters ---------- degree : int The degree of the polynomial series, a.k.a. baseline order - exclude : list of 2-tuples of int or ~astropy.units.Quantity, or ~specutils.SpectralRegion + exclude : list of 2-tuples of int or `~astropy.units.Quantity`, or `~specutils.SpectralRegion` List of region(s) to exclude from the fit. The tuple(s) represent a range in the form [lower,upper], inclusive. - In channel units. + If both `exclude` and `include` are given, only `include` is used. + + Examples: + + One channel-based region: [11,51] + + Two channel-based regions: [(11,51),(99,123)]. + + One `~astropy.units.Quantity` region: [110.198*u.GHz,110.204*u.GHz]. + + One compound `~specutils.SpectralRegion`: SpectralRegion([(110.198*u.GHz,110.204*u.GHz),(110.196*u.GHz,110.197*u.GHz)]). + + Default: no `exclude` region + + include: list of 2-tuples of int or `~astropy.units.Quantity`, or `~specutils.SpectralRegion` + List of region(s) to include in the fit. The tuple(s) represent a range in the form [lower,upper], inclusive. + If both `exclude` and `include` are given, only `include` is used. Examples: @@ -206,13 +200,12 @@ def baseline(self, degree, exclude=None, include=None, **kwargs): Two channel-based regions: [(11,51),(99,123)]. - One ~astropy.units.Quantity region: [110.198*u.GHz,110.204*u.GHz]. + One `~astropy.units.Quantity` region: [110.198*u.GHz,110.204*u.GHz]. One compound `~specutils.SpectralRegion`: SpectralRegion([(110.198*u.GHz,110.204*u.GHz),(110.196*u.GHz,110.197*u.GHz)]). - Default: no exclude region + Default: no `include` region - include: list of 2-tuples of int (currently units not supported yet, pending issue 251/260) model : str One of 'polynomial' 'chebyshev', 'legendre', or 'hermite' @@ -257,7 +250,6 @@ def baseline(self, degree, exclude=None, include=None, **kwargs): print(f"Warning: ignoring exclude={exclude}") nchan = len(self._spectral_axis) exclude = core.include_to_exclude_spectral_region(include, self) - # exclude = self._toggle_sections(nchan, include) self._baseline_model = baseline(self, degree, exclude, **kwargs) From 7ec23e3102b311f7f71330a62be1639c30cececf Mon Sep 17 00:00:00 2001 From: astrofle Date: Fri, 17 Jan 2025 14:30:50 -0500 Subject: [PATCH 6/8] Fix: docstrings --- src/dysh/spectra/core.py | 48 +++++++++++++++++++++------------------- 1 file changed, 25 insertions(+), 23 deletions(-) diff --git a/src/dysh/spectra/core.py b/src/dysh/spectra/core.py index 48e63297..050113b9 100644 --- a/src/dysh/spectra/core.py +++ b/src/dysh/spectra/core.py @@ -209,7 +209,7 @@ def include_to_exclude_spectral_region(include, refspec): refspec: `~spectra.spectrum.Spectrum` The reference spectrum whose spectral axis will be used - when converting between include and axis units (e.g. channels to GHz). + when converting between `include` and axis units (e.g. channels to GHz). Returns ------- @@ -237,12 +237,12 @@ def include_to_exclude_spectral_region(include, refspec): def exclude_to_spectral_region(exclude, refspec, fix_exclude=True): - """Convert an exclude list to a list of ~specutuls.SpectralRegion. + """Convert `exclude` to a `~specutils.SpectralRegion`. Parameters ---------- exclude : list of 2-tuples of int or `~astropy.units.Quantity`, or `~specutils.SpectralRegion` - List of region(s) to exclude from the fit. The tuple(s) represent a range in the form [lower,upper], inclusive. + List of region(s) to exclude. The tuple(s) represent a range in the form [lower,upper], inclusive. Examples: One channel-based region: @@ -263,14 +263,15 @@ def exclude_to_spectral_region(exclude, refspec, fix_exclude=True): refspec: `~spectra.spectrum.Spectrum` The reference spectrum whose spectral axis will be used - when converting between exclude and axis units (e.g. channels to GHz). + when converting between `exclude` and axis units (e.g. channels to GHz). fix_exclude: bool - If True, fix exclude regions that are out of bounds of the specctral axis to be within the spectral axis. Default:False + If True, fix exclude regions that are out of bounds of the spectral axis to be within the spectral axis. + Default: True Returns ------- - sr : `~specutil.SpectralRegion` - A `~specutil.SpectralRegion` corresponding to `exclude`. + sr : `~specutils.SpectralRegion` + A `~specutils.SpectralRegion` corresponding to `exclude`. """ p = refspec @@ -320,8 +321,8 @@ def spectral_region_to_unit(spectral_region, refspec, unit=None): Parameters ---------- - spectral_region : `~specutil.SpectralRegion` - `~specutil.SpectralRegion` whose units will be converted. + spectral_region : `~specutils.SpectralRegion` + `~specutils.SpectralRegion` whose units will be converted. refspec : `~spectra.spectrum.Spectrum` The reference spectrum whose spectral axis will be used when converting to `unit` (e.g. channels to GHz). @@ -330,7 +331,7 @@ def spectral_region_to_unit(spectral_region, refspec, unit=None): Returns ------- - spectral_region : `~specutil.SpectralRegion` + spectral_region : `~specutils.SpectralRegion` SpectralRegion with units of `unit`. """ @@ -347,18 +348,18 @@ def spectral_region_to_unit(spectral_region, refspec, unit=None): def spectral_region_to_list(spectral_region): """ - Turn `spectral_region` into a list of `~specutil.SpectralRegion`. + Turn `spectral_region` into a list of `~specutils.SpectralRegion`. Each subregion in `spectral_region` will be a list element. Parameters ---------- - spectral_region : `~specutil.SpectralRegion` - `~specutil.SpectralRegion` to convert into a list of `~specutil.SpectralRegion`. + spectral_region : `~specutils.SpectralRegion` + `~specutils.SpectralRegion` to convert into a list of `~specutils.SpectralRegion`. Returns ------- - region_list : list of `~specutil.SpectralRegion` - Subregions of `spectral_region` in a list of `~specutil.SpectralRegion`. + region_list : list of `~specutils.SpectralRegion` + Subregions of `spectral_region` in a list of `~specutils.SpectralRegion`. """ region_list = [] @@ -403,13 +404,13 @@ def exclude_to_mask(exclude, refspec): def exclude_to_region_list(exclude, spectrum, fix_exclude=True): """ - Convert an exclusion region, `exclude`, to a list of `~SpectralRegion`. + Convert an exclusion region, `exclude`, to a list of `~specutils.SpectralRegion`. This is used for baseline fitting. Parameters ---------- exclude : list of 2-tuples of int or `~astropy.units.Quantity`, or `~specutils.SpectralRegion` - List of region(s) to exclude from the fit. The tuple(s) represent a range in the form [lower,upper], inclusive. + List of region(s) to exclude. The tuple(s) represent a range in the form [lower,upper], inclusive. Examples: One channel-based region: @@ -430,14 +431,15 @@ def exclude_to_region_list(exclude, spectrum, fix_exclude=True): spectrum : `~spectra.spectrum.Spectrum` The reference spectrum whose spectral axis will be used - when converting between exclude and axis units (e.g. channels to GHz). + when converting between `exclude` and axis units (e.g. channels to GHz). fix_exclude : bool - See `~spectra.core.exclude_to_spectral_region` for details. Default: True + See `~spectra.core.exclude_to_spectral_region` for details. + Default: True Returns ------- - region_list : list of `~specutil.SpectralRegion` - Regions defined in `exclude` as a list of `~specutil.SpectralRegion`. + region_list : list of `~specutils.SpectralRegion` + Regions defined in `exclude` as a list of `~specutils.SpectralRegion`. """ spectral_region = exclude_to_spectral_region(exclude, spectrum, fix_exclude=fix_exclude) @@ -458,7 +460,7 @@ def baseline(spectrum, order, exclude=None, exclude_region_upper_bounds=True, ** order : int The order of the polynomial series, a.k.a. baseline order. exclude : list of 2-tuples of int or `~astropy.units.Quantity`, or `~specutils.SpectralRegion` - List of region(s) to exclude from the fit. The tuple(s) represent a range in the form [lower,upper], inclusive. + List of region(s) to exclude from the fit. The tuple(s) represent a range in the form [lower,upper], inclusive. Examples: One channel-based region: @@ -492,7 +494,7 @@ def baseline(spectrum, order, exclude=None, exclude_region_upper_bounds=True, ** ------- models : list of `~astropy.modeling.Model` The list of models that contain the fitted model parameters. - See `~specutuls.fitting.fit_continuum`. + See `~specutils.fitting.fit_continuum`. """ kwargs_opts = { From c1538c7c2bc515ff3127b73430e7b7f2c72ce58a Mon Sep 17 00:00:00 2001 From: astrofle Date: Fri, 17 Jan 2025 14:31:20 -0500 Subject: [PATCH 7/8] Add: test for baseline fitting using regions with units --- src/dysh/spectra/tests/test_spectrum.py | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/src/dysh/spectra/tests/test_spectrum.py b/src/dysh/spectra/tests/test_spectrum.py index 20295fa1..e83cf2df 100644 --- a/src/dysh/spectra/tests/test_spectrum.py +++ b/src/dysh/spectra/tests/test_spectrum.py @@ -632,6 +632,16 @@ def test_single_baseline(sdf, ex_reg, order, model, gbtidl_bmodel, exclude_regio ex_reg = None test_single_baseline(sdf, ex_reg, order, "chebyshev", gbtidl_no_reg) + # Test with units. Nothing to compare to though. + dysh_spec = sdf.getspec(0) + kms = u.km / u.s + ex_reg = [(0 * kms, 4200 * kms), (6000 * kms, 7000 * kms), (8800 * kms, 90000 * kms)] + dysh_spec.baseline(order, ex_reg, model="chebyshev") + ex_reg = [(1 * u.GHz, 1.38 * u.GHz), (1.388 * u.GHz, 1.392 * u.GHz), (1.4 * u.GHz, 2 * u.GHz)] + dysh_spec.baseline(order, ex_reg, model="chebyshev") + ex_reg = [21 * u.cm, 21.5 * u.cm] + dysh_spec.baseline(order, ex_reg, model="chebyshev") + def test_baseline_include(self): """ Test for comparing GBTIDL baseline to dysh baselines @@ -661,3 +671,13 @@ def test_single_baseline(sdf, in_reg, order, model, gbtidl_bmodel, exclude_regio test_single_baseline(sdf, in_reg, order, "chebyshev", gbtidl_two_reg) test_single_baseline(sdf, in_reg, order, "legendre", gbtidl_two_reg) test_single_baseline(sdf, in_reg, order, "hermite", gbtidl_two_reg) + + # Test with units. Nothing to compare to though. + dysh_spec = sdf.getspec(0) + kms = u.km / u.s + in_reg = [(4200 * kms, 6000 * kms), (7000 * kms, 8800 * kms)] + dysh_spec.baseline(order, include=in_reg, model="chebyshev") + in_reg = [(1.38 * u.GHz, 1.388 * u.GHz), (1.392 * u.GHz, 1.4 * u.GHz)] + dysh_spec.baseline(order, include=in_reg, model="chebyshev") + in_reg = [21 * u.cm, 21.5 * u.cm] + dysh_spec.baseline(order, include=in_reg, model="chebyshev") From 3d4d2fbd35707d19ceece2814f27d2b5ed942224 Mon Sep 17 00:00:00 2001 From: Pedro Salas Date: Wed, 22 Jan 2025 13:58:53 -0500 Subject: [PATCH 8/8] Fix: change print to logger.info --- src/dysh/spectra/core.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/dysh/spectra/core.py b/src/dysh/spectra/core.py index 050113b9..f61eb46b 100644 --- a/src/dysh/spectra/core.py +++ b/src/dysh/spectra/core.py @@ -539,7 +539,7 @@ def baseline(spectrum, order, exclude=None, exclude_region_upper_bounds=True, ** # use the spectrum's preset exclude regions if they # exist (they will be a list of SpectralRegions or None) regionlist = p._exclude_regions - print(f"EXCLUDING {regionlist}") + logger.info(f"EXCLUDING {regionlist}") return fit_continuum( spectrum=p, model=selected_model,