From fcea42b732351b682e5b36974ac9cac0a87e3ace Mon Sep 17 00:00:00 2001 From: icweaver Date: Sat, 1 Nov 2025 03:31:45 -0700 Subject: [PATCH 01/25] started resampler ext --- Project.toml | 7 +++++++ ext/DataInterpolationsExt.jl | 11 +++++++++++ src/Spectra.jl | 32 ++++++++++++++++++++++++++++++++ 3 files changed, 50 insertions(+) create mode 100644 ext/DataInterpolationsExt.jl diff --git a/Project.toml b/Project.toml index cfbb101..4174e9d 100644 --- a/Project.toml +++ b/Project.toml @@ -15,7 +15,14 @@ RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" UnitfulAstro = "6112ee07-acf9-5e0f-b108-d242c714bf9f" +[weakdeps] +DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0" + +[extensions] +DataInterpolationsExt = "DataInterpolations" + [compat] +DataInterpolations = "8.6.1" DustExtinction = ">=0.6.0" LinearAlgebra = "1.5" Measurements = "2" diff --git a/ext/DataInterpolationsExt.jl b/ext/DataInterpolationsExt.jl new file mode 100644 index 0000000..78cc850 --- /dev/null +++ b/ext/DataInterpolationsExt.jl @@ -0,0 +1,11 @@ +module DataInterpolationsExt + using Spectra: Spectra, Spectrum, spectrum + using DataInterpolations: AbstractInterpolation + + function Spectra.resample(spec::Spectrum, wave_sampled, interp::AbstractInterpolation) + p = sortperm(spec.wave) # x-scale must be monitonically increasing + flux_sampled = interp(wave_sampled) + return spectrum(wave_sampled, flux_sampled; meta = spec.meta) + end + +end # module diff --git a/src/Spectra.jl b/src/Spectra.jl index f804c52..ccea2e4 100644 --- a/src/Spectra.jl +++ b/src/Spectra.jl @@ -6,6 +6,8 @@ export AbstractSpectrum, spectrum export blackbody # transforms/redden.jl export redden, redden!, deredden, deredden! +# resampling: ../ext/DataInterpolationsExt.jl +export resample using RecipesBase: @recipe using Measurements: Measurements, Measurement @@ -91,6 +93,36 @@ function spectrum(wave::AbstractMatrix{<:Quantity}, flux::AbstractMatrix{<:Quant EchelleSpectrum(wave, flux, Dict{Symbol,Any}(kwds)) end +# Stub +""" + function resample(spec, wave_sampled, interp) + +Resample a spectrum onto the given wavelength grid `wave_sampled` using the supplied interpolator `interp`. Currently supports interpolators from the [DataInterpolations.jl](https://github.com/SciML/DataInterpolations.jl) package. + +# Examples + +```julia-repl +julia> using DataInterpolations + +julia> spec = spectrum([2, 4, 12, 16, 20], [1, 3, 7, 6, 20]); + +julia> wave_sampled = [1, 5, 9, 13, 14, 17, 21, 22, 23]; + + # BYO interpolator +julia> interp = LinearInterpolation(spec.flux, spec.wave; extrapolation = ExtrapolationType.Constant); + +julia> spec_sampled = resample(spec, wave_sampled, interp); +``` +""" +function resample(spec, wave_sampled, interp) + error("""Supported interpolation package not loaded. Please try adding one of the supported interpolation packages below: + + - DataInterpolations.jl + + PRs to add additional interpolation packages are welcome! + """) +end + # tools include("utils.jl") include("transforms/transforms.jl") From e40a7246c837fb6c7baf7f9412706f0496393078 Mon Sep 17 00:00:00 2001 From: cgarling Date: Sat, 1 Nov 2025 11:42:17 -0400 Subject: [PATCH 02/25] Another resampling idea --- Project.toml | 2 +- ext/DataInterpolationsExt.jl | 41 +++++++++++++++++++++++++++++------- 2 files changed, 34 insertions(+), 9 deletions(-) diff --git a/Project.toml b/Project.toml index 4174e9d..70b1183 100644 --- a/Project.toml +++ b/Project.toml @@ -22,7 +22,7 @@ DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0" DataInterpolationsExt = "DataInterpolations" [compat] -DataInterpolations = "8.6.1" +DataInterpolations = "8" DustExtinction = ">=0.6.0" LinearAlgebra = "1.5" Measurements = "2" diff --git a/ext/DataInterpolationsExt.jl b/ext/DataInterpolationsExt.jl index 78cc850..cb290be 100644 --- a/ext/DataInterpolationsExt.jl +++ b/ext/DataInterpolationsExt.jl @@ -1,11 +1,36 @@ module DataInterpolationsExt - using Spectra: Spectra, Spectrum, spectrum - using DataInterpolations: AbstractInterpolation - - function Spectra.resample(spec::Spectrum, wave_sampled, interp::AbstractInterpolation) - p = sortperm(spec.wave) # x-scale must be monitonically increasing - flux_sampled = interp(wave_sampled) - return spectrum(wave_sampled, flux_sampled; meta = spec.meta) - end + +using Spectra: Spectra, Spectrum, spectrum, wave, flux +import Spectra: resample +using DataInterpolations: AbstractInterpolation + +function Spectra.resample(spec::Spectrum, wave_sampled, interp::AbstractInterpolation) + # p = sortperm(spec.wave) # x-scale must be monitonically increasing; not being used + flux_sampled = interp(wave_sampled) + return spectrum(wave_sampled, flux_sampled; meta = spec.meta) +end + +""" + resample(spec::Spectrum, wave_sampled, interp::Type{<:DataInterpolations.AbstractInterpolation}; kws...) +Constructs an interpolator of type `interp` from the provided spectrum `spec` and evaluates it at `wave_sampled`, returning a `Spectrum`. The `kws...` control extrapolation and are passed through to the interpolator. + +```jldoctest +julia> using Spectra, DataInterpolations + +julia> spec = spectrum([2, 4, 12, 16, 20], [1, 3, 7, 6, 20]); + +julia> wave_sampled = [1, 5, 9, 13, 14, 17, 21, 22, 23]; + +julia> resample(spec, wave_sampled, LinearInterpolation; extrapolation = ExtrapolationType.Constant) +Spectrum(Int64, Float64) +``` +""" +function Spectra.resample(spec::Spectrum, wave_sampled, interp::Type{<:AbstractInterpolation}; kws...) + # Not sure if we have to check wave is monotonically increasing or if that is guaranteed by the constructor? + p = sortperm(wave(spec)) + itp = interp(flux(spec)[p], wave(spec)[p]; kws...) + flux_sampled = itp.(wave_sampled) + return spectrum(wave_sampled, flux_sampled; meta = spec.meta) +end end # module From fdc9c9ceeab6acf198cfc9c521a5fa64ef71d450 Mon Sep 17 00:00:00 2001 From: Ian Weaver Date: Sat, 1 Nov 2025 22:56:08 -0700 Subject: [PATCH 03/25] Update src/Spectra.jl Co-authored-by: Chris Garling --- src/Spectra.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Spectra.jl b/src/Spectra.jl index ccea2e4..968bf2f 100644 --- a/src/Spectra.jl +++ b/src/Spectra.jl @@ -109,7 +109,7 @@ julia> spec = spectrum([2, 4, 12, 16, 20], [1, 3, 7, 6, 20]); julia> wave_sampled = [1, 5, 9, 13, 14, 17, 21, 22, 23]; # BYO interpolator -julia> interp = LinearInterpolation(spec.flux, spec.wave; extrapolation = ExtrapolationType.Constant); +julia> interp = LinearInterpolation(flux(spec), wave(spec); extrapolation = ExtrapolationType.Constant); julia> spec_sampled = resample(spec, wave_sampled, interp); ``` From 333c95748349812f70b7378d2b0c421c4cdef637 Mon Sep 17 00:00:00 2001 From: Ian Weaver Date: Sat, 1 Nov 2025 22:56:33 -0700 Subject: [PATCH 04/25] Update src/Spectra.jl Co-authored-by: Chris Garling --- src/Spectra.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Spectra.jl b/src/Spectra.jl index 968bf2f..e160597 100644 --- a/src/Spectra.jl +++ b/src/Spectra.jl @@ -115,7 +115,7 @@ julia> spec_sampled = resample(spec, wave_sampled, interp); ``` """ function resample(spec, wave_sampled, interp) - error("""Supported interpolation package not loaded. Please try adding one of the supported interpolation packages below: + error("""No supported interpolation package loaded. To enable resampling, load one of the supported interpolation packages below: - DataInterpolations.jl From 6426edba311e9df005b58b9c27937353fec991eb Mon Sep 17 00:00:00 2001 From: icweaver Date: Sun, 2 Nov 2025 01:33:41 -0800 Subject: [PATCH 05/25] added Interpolations ext, cleaned up interface a bit --- Project.toml | 3 ++ ext/DataInterpolationsExt.jl | 33 ++++++++++++--- ext/InterpolationsExt.jl | 80 ++++++++++++++++++++++++++++++++++++ src/Spectra.jl | 15 ++++--- 4 files changed, 120 insertions(+), 11 deletions(-) create mode 100644 ext/InterpolationsExt.jl diff --git a/Project.toml b/Project.toml index 70b1183..73f6a88 100644 --- a/Project.toml +++ b/Project.toml @@ -17,13 +17,16 @@ UnitfulAstro = "6112ee07-acf9-5e0f-b108-d242c714bf9f" [weakdeps] DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0" +Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" [extensions] DataInterpolationsExt = "DataInterpolations" +InterpolationsExt = "Interpolations" [compat] DataInterpolations = "8" DustExtinction = ">=0.6.0" +Interpolations = "0.16.2" LinearAlgebra = "1.5" Measurements = "2" PhysicalConstants = "0.2" diff --git a/ext/DataInterpolationsExt.jl b/ext/DataInterpolationsExt.jl index cb290be..9d8ac32 100644 --- a/ext/DataInterpolationsExt.jl +++ b/ext/DataInterpolationsExt.jl @@ -2,10 +2,31 @@ module DataInterpolationsExt using Spectra: Spectra, Spectrum, spectrum, wave, flux import Spectra: resample -using DataInterpolations: AbstractInterpolation +using DataInterpolations: DataInterpolations, AbstractInterpolation +""" + function Spectra.resample(spec::Spectrum, wave_sampled, interp::AbstractInterpolation) + +Resample `spec` using the given interpolator `interp`. + +```jldoctest +julia> using Spectra: spectrum, resample, wave, flux + +julia> using DataInterpolations + +julia> spec = spectrum([2, 4, 12, 16, 20], [1, 3, 7, 6, 20]); + +julia> wave_sampled = [1, 5, 9, 13, 14, 17, 21, 22, 23]; + +julia> interp = LinearInterpolation(flux(spec), wave(spec); extrapolation = ExtrapolationType.Constant); + +julia> resample(spec, wave_sampled, interp) +Spectrum(Int64, Float64) + +See [resample(spec::Spectrum, wave_sampled, interp::Type{<:DataInterpolations.AbstractInterpolation} ; kws...)](@ref) for a convenience function version. +``` +""" function Spectra.resample(spec::Spectrum, wave_sampled, interp::AbstractInterpolation) - # p = sortperm(spec.wave) # x-scale must be monitonically increasing; not being used flux_sampled = interp(wave_sampled) return spectrum(wave_sampled, flux_sampled; meta = spec.meta) end @@ -23,13 +44,13 @@ julia> wave_sampled = [1, 5, 9, 13, 14, 17, 21, 22, 23]; julia> resample(spec, wave_sampled, LinearInterpolation; extrapolation = ExtrapolationType.Constant) Spectrum(Int64, Float64) + +For better performance, see [Spectra.resample(spec::Spectrum, wave_sampled, interp::DataInterpolations.AbstractInterpolation)](@ref). ``` """ function Spectra.resample(spec::Spectrum, wave_sampled, interp::Type{<:AbstractInterpolation}; kws...) - # Not sure if we have to check wave is monotonically increasing or if that is guaranteed by the constructor? - p = sortperm(wave(spec)) - itp = interp(flux(spec)[p], wave(spec)[p]; kws...) - flux_sampled = itp.(wave_sampled) + itp = interp(flux(spec), wave(spec); kws...) + flux_sampled = itp(wave_sampled) return spectrum(wave_sampled, flux_sampled; meta = spec.meta) end diff --git a/ext/InterpolationsExt.jl b/ext/InterpolationsExt.jl new file mode 100644 index 0000000..b43ec9a --- /dev/null +++ b/ext/InterpolationsExt.jl @@ -0,0 +1,80 @@ +module InterpolationsExt + +using Spectra: Spectra, Spectrum, spectrum, wave, flux +import Spectra: resample +using Interpolations: + Interpolations, + AbstractInterpolation, + Degree, + Constant, + Linear, + Cubic, + constant_interpolation, + linear_interpolation, + cubic_interpolation + +""" + Spectra.resample(spec::Spectrum, wave_sampled, interp::AbstractInterpolation) + +Resample `spec` using the given interpolator `interp`. + +```jldoctest +julia> using Spectra: spectrum, resample, wave, flux + +julia> using Interpolations + +julia> spec = spectrum([2, 4, 12, 16, 20], [1, 3, 7, 6, 20]); + +julia> wave_sampled = [1, 5, 9, 13, 14, 17, 21, 22, 23]; + +julia> interp = extrapolate(interpolate((wave(spec), ), flux(spec), Gridded(Linear())), Flat()); + +julia> resample(spec, wave_sampled, interp) +Spectrum(Int64, Float64) +``` + +See [resample(spec::Spectrum, wave_sampled, deg::Interpolations.Degree; kws...)](@ref) for a convenience function version with limited flexibility. +""" +function Spectra.resample(spec::Spectrum, wave_sampled, interp::AbstractInterpolation) + flux_sampled = interp(wave_sampled) + return spectrum(wave_sampled, flux_sampled; meta = spec.meta) +end + +""" + resample(spec::Spectrum, wave_sampled, deg::Interpolations.Degree; kws...) + +Constructs an interpolator of type `interp` from the provided spectrum `spec` and evaluates it at `wave_sampled`, returning a `Spectrum`. The `kws...` control extrapolation and are passed through to the interpolator. + +```jldoctest +julia> using Spectra: spectrum, resample + +julia> using Interpolations: Flat + +julia> spec = spectrum([2, 4, 12, 16, 20], [1, 3, 7, 6, 20]); + +julia> wave_sampled = [1, 5, 9, 13, 14, 17, 21, 22, 23]; + +julia> resample(spec, wave_sampled, Linear(); extrapolation_bc = Flat()) +Spectrum(Int64, Float64) +``` + +This convenience function currently only supports `Constant`, `Linear`, and `Cubic` splines from Interpolations.jl. For more interpolators and better performance, see [Spectra.resample(spec::Spectrum, wave_sampled, interp::Interpolations.AbstractInterpolation)](@ref). +""" +function Spectra.resample(spec::Spectrum, wave_sampled, deg::Degree; kws...) + + deg_func = if deg isa Constant + constant_interpolation + elseif deg isa Linear + linear_interpolation + elseif deg isa Cubic + cubic_interpolation + else + throw(ArgumentError(deg, "Invalid degree passed. Must be one of `Constant`, `Linear`, `Quadratic`, or `Cubic`.")) + end + + itp = deg_func(wave(spec), flux(spec); kws...) + flux_sampled = itp(wave_sampled) + return spectrum(wave_sampled, flux_sampled; meta = spec.meta) +end + +end # module diff --git a/src/Spectra.jl b/src/Spectra.jl index e160597..0ba81c6 100644 --- a/src/Spectra.jl +++ b/src/Spectra.jl @@ -6,7 +6,7 @@ export AbstractSpectrum, spectrum export blackbody # transforms/redden.jl export redden, redden!, deredden, deredden! -# resampling: ../ext/DataInterpolationsExt.jl +# resampling: ../ext/DataInterpolationsExt.jl ../ext/Interpolations.jl export resample using RecipesBase: @recipe @@ -97,7 +97,10 @@ end """ function resample(spec, wave_sampled, interp) -Resample a spectrum onto the given wavelength grid `wave_sampled` using the supplied interpolator `interp`. Currently supports interpolators from the [DataInterpolations.jl](https://github.com/SciML/DataInterpolations.jl) package. +Resample a spectrum onto the given wavelength grid `wave_sampled` using the supplied interpolator `interp`. Currently supports interpolators from the following packages: + + * [DataInterpolations.jl](https://github.com/SciML/DataInterpolations.jl) + * [Interpolations.jl](https://github.com/JuliaMath/Interpolations.jl) # Examples @@ -109,15 +112,17 @@ julia> spec = spectrum([2, 4, 12, 16, 20], [1, 3, 7, 6, 20]); julia> wave_sampled = [1, 5, 9, 13, 14, 17, 21, 22, 23]; # BYO interpolator -julia> interp = LinearInterpolation(flux(spec), wave(spec); extrapolation = ExtrapolationType.Constant); +julia> interp = ... -julia> spec_sampled = resample(spec, wave_sampled, interp); +julia> spec_sampled = resample(spec, wave_sampled, interp) +Spectrum(Int64, Float64) ``` """ function resample(spec, wave_sampled, interp) error("""No supported interpolation package loaded. To enable resampling, load one of the supported interpolation packages below: - - DataInterpolations.jl + * [DataInterpolations.jl](https://github.com/SciML/DataInterpolations.jl) + * [Interpolations.jl](https://github.com/JuliaMath/Interpolations.jl) PRs to add additional interpolation packages are welcome! """) From 248b9fe568c30862127959fea717edf73d8ad50c Mon Sep 17 00:00:00 2001 From: cgarling Date: Sun, 2 Nov 2025 13:25:06 -0500 Subject: [PATCH 06/25] Try `SpectrumResampler` type --- src/Spectra.jl | 42 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 42 insertions(+) diff --git a/src/Spectra.jl b/src/Spectra.jl index 0ba81c6..f5cd3f9 100644 --- a/src/Spectra.jl +++ b/src/Spectra.jl @@ -128,6 +128,48 @@ function resample(spec, wave_sampled, interp) """) end +""" + SpectrumResampler(spec::Spectrum, interp) +Type representing the spectrum `spec` with interpolator `interp`. + +```jldoctest +julia> using Spectra: spectrum, flux, wave, SpectrumResampler + +julia> using DataInterpolations: LinearInterpolation, ExtrapolationType + +julia> spec = spectrum([2, 4, 12, 16, 20], [1, 3, 7, 6, 20]); + +julia> wave_sampled = [1, 5, 9, 13, 14, 17, 21, 22, 23]; + +julia> interp = LinearInterpolation(flux(spec), wave(spec); extrapolation = ExtrapolationType.Constant); + +julia> resampler = SpectrumResampler(spec, interp); + +julia> result = resampler(wave_sampled); + +julia> result isa SpectrumResampler +true + +julia> wave(result) == wave_sampled +true + +julia> flux(result) == interp.(wave_sampled) +true +``` +""" +struct SpectrumResampler{A <: Spectrum, B} + spectrum::A + interp::B +end +wave(s::SpectrumResampler) = wave(s.spectrum) +flux(s::SpectrumResampler) = flux(s.spectrum) +Spectrum(s::SpectrumResampler) = s.spectrum +function (spec::SpectrumResampler)(wave_sampled) + newspec = spectrum(wave_sampled, spec.interp.(wave_sampled); meta = Spectrum(spec).meta) + return SpectrumResampler(newspec, spec.interp) +end +resample(spec::SpectrumResampler, wave_sampled) = spec(wave_sampled) + # tools include("utils.jl") include("transforms/transforms.jl") From 623b4469f4fdcfdd21d7eb15f56efee5afa9e3c8 Mon Sep 17 00:00:00 2001 From: cgarling Date: Sun, 2 Nov 2025 13:43:00 -0500 Subject: [PATCH 07/25] Use the `spectrum` constructor not `Spectrum` --- src/Spectra.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Spectra.jl b/src/Spectra.jl index f5cd3f9..30ff569 100644 --- a/src/Spectra.jl +++ b/src/Spectra.jl @@ -163,9 +163,9 @@ struct SpectrumResampler{A <: Spectrum, B} end wave(s::SpectrumResampler) = wave(s.spectrum) flux(s::SpectrumResampler) = flux(s.spectrum) -Spectrum(s::SpectrumResampler) = s.spectrum +spectrum(s::SpectrumResampler) = s.spectrum function (spec::SpectrumResampler)(wave_sampled) - newspec = spectrum(wave_sampled, spec.interp.(wave_sampled); meta = Spectrum(spec).meta) + newspec = spectrum(wave_sampled, spec.interp.(wave_sampled); meta = spectrum(spec).meta) return SpectrumResampler(newspec, spec.interp) end resample(spec::SpectrumResampler, wave_sampled) = spec(wave_sampled) From 24714e6f52ba2a873a1fff61d4853d3b91c292b4 Mon Sep 17 00:00:00 2001 From: icweaver Date: Sun, 2 Nov 2025 20:39:40 -0800 Subject: [PATCH 08/25] switched to explicit interpolation + docs --- Project.toml | 10 ---- docs/Project.toml | 2 + docs/src/transforms.md | 68 +++++++++++++++++++++++++ ext/DataInterpolationsExt.jl | 57 --------------------- src/Spectra.jl | 96 +++++------------------------------- src/transforms/resampler.jl | 51 +++++++++++++++++++ 6 files changed, 132 insertions(+), 152 deletions(-) delete mode 100644 ext/DataInterpolationsExt.jl create mode 100644 src/transforms/resampler.jl diff --git a/Project.toml b/Project.toml index 73f6a88..cfbb101 100644 --- a/Project.toml +++ b/Project.toml @@ -15,18 +15,8 @@ RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" UnitfulAstro = "6112ee07-acf9-5e0f-b108-d242c714bf9f" -[weakdeps] -DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0" -Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" - -[extensions] -DataInterpolationsExt = "DataInterpolations" -InterpolationsExt = "Interpolations" - [compat] -DataInterpolations = "8" DustExtinction = ">=0.6.0" -Interpolations = "0.16.2" LinearAlgebra = "1.5" Measurements = "2" PhysicalConstants = "0.2" diff --git a/docs/Project.toml b/docs/Project.toml index 4d4d0b2..837fa96 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,6 +1,8 @@ [deps] +DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" FITSIO = "525bcba6-941b-5504-bd06-fd0dc1a4d2eb" +Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Spectra = "391af1a9-06f1-59d3-8d21-0be089654739" diff --git a/docs/src/transforms.md b/docs/src/transforms.md index a5a05c5..5fa3a63 100644 --- a/docs/src/transforms.md +++ b/docs/src/transforms.md @@ -56,6 +56,73 @@ julia> red.flux โ‰ˆ spec.flux true ``` +## Resampling + +External interpolators, e.g., from [DataInterpolations.jl](https://github.com/SciML/DataInterpolations.jl), [Interpolations.jl](https://github.com/JuliaMath/Interpolations.jl), can be used to resample spectra onto a given wavelength grid. Starting with a sample spectrum `spec`, we first create a [`SpectrumResampler`](@ref) object `resampler` which stores the initial spectrum and interpolator `interp` together. We then apply this object to the wavelength grid our our choice to produce our desired resampled spectrum. Below, we show how this may be accomplished using either package. + +### DataInterpolations.jl + +```@repl resample +using Spectra: SpectrumResampler, spectrum, wave, flux + +using DataInterpolations: LinearInterpolation, ExtrapolationType + +spec = spectrum([20, 40, 120, 160, 200], [1, 3, 7, 6, 20]) + +interp = LinearInterpolation(flux(spec), wave(spec); extrapolation = ExtrapolationType.Constant) + +resampler = SpectrumResampler(spec, interp) +``` + +We can now apply `resampler` to our desired wavelength grid `wave_sampled` to produce our resampled spectrum: + +```@repl resample +wave_sampled = [10, 50, 90, 130, 140, 170, 210, 220, 230]; + +spec_resampled = resampler(wave_sampled) +``` + +Calling this on different wavelength grids uses the same stored interpolator, allowing for efficient resampling. For convenience, the resampled spectrum is returned as another `SpectrumResampler` object, which keeps information about the resampled spectrum and applied interpolator coupled. The resulting interpolated flux and wavelength grid can be retrieved with the `flux` and `wave` getters exported by Spectra.jl, respectively: + +```@repl resample +wave(spec_resampled) + +flux(spec_resampled) + +# Check that this is equivalent to calling the interpolator directly +flux(spec_resampled) == interp(wave_sampled) +``` + +### Interpolations.jl + +For completeness, here is how a similar resampling scheme could be accomplished with Interpolations.jl: + +```@repl +using Spectra: SpectrumResampler, spectrum, wave, flux + +using Interpolations: linear_interpolation, Flat + +spec = spectrum([20, 40, 120, 160, 200], [1, 3, 7, 6, 20]) + +# Note that `wave` and `flux` are flipped relative to the DataInterpolations.jl example +interp = linear_interpolation(wave(spec), flux(spec); extrapolation_bc = Flat()) + +resampler = SpectrumResampler(spec, interp) + +wave_sampled = [10, 50, 90, 130, 140, 170, 210, 220, 230]; + +spec_resampled = resampler(wave_sampled) + +wave(spec_resampled) + +flux(spec_resampled) + +# Check that this is equivalent to calling the interpolator directly +flux(spec_resampled) == interp(wave_sampled) +``` + +See the relevant documentation for the external interpolation package used for more information on what configurations can be passed for your desired `interp` object. + ### API/Reference ```@docs @@ -63,6 +130,7 @@ redden redden! deredden deredden! +SpectrumResampler ``` ```@meta diff --git a/ext/DataInterpolationsExt.jl b/ext/DataInterpolationsExt.jl deleted file mode 100644 index 9d8ac32..0000000 --- a/ext/DataInterpolationsExt.jl +++ /dev/null @@ -1,57 +0,0 @@ -module DataInterpolationsExt - -using Spectra: Spectra, Spectrum, spectrum, wave, flux -import Spectra: resample -using DataInterpolations: DataInterpolations, AbstractInterpolation - -""" - function Spectra.resample(spec::Spectrum, wave_sampled, interp::AbstractInterpolation) - -Resample `spec` using the given interpolator `interp`. - -```jldoctest -julia> using Spectra: spectrum, resample, wave, flux - -julia> using DataInterpolations - -julia> spec = spectrum([2, 4, 12, 16, 20], [1, 3, 7, 6, 20]); - -julia> wave_sampled = [1, 5, 9, 13, 14, 17, 21, 22, 23]; - -julia> interp = LinearInterpolation(flux(spec), wave(spec); extrapolation = ExtrapolationType.Constant); - -julia> resample(spec, wave_sampled, interp) -Spectrum(Int64, Float64) - -See [resample(spec::Spectrum, wave_sampled, interp::Type{<:DataInterpolations.AbstractInterpolation} ; kws...)](@ref) for a convenience function version. -``` -""" -function Spectra.resample(spec::Spectrum, wave_sampled, interp::AbstractInterpolation) - flux_sampled = interp(wave_sampled) - return spectrum(wave_sampled, flux_sampled; meta = spec.meta) -end - -""" - resample(spec::Spectrum, wave_sampled, interp::Type{<:DataInterpolations.AbstractInterpolation}; kws...) -Constructs an interpolator of type `interp` from the provided spectrum `spec` and evaluates it at `wave_sampled`, returning a `Spectrum`. The `kws...` control extrapolation and are passed through to the interpolator. - -```jldoctest -julia> using Spectra, DataInterpolations - -julia> spec = spectrum([2, 4, 12, 16, 20], [1, 3, 7, 6, 20]); - -julia> wave_sampled = [1, 5, 9, 13, 14, 17, 21, 22, 23]; - -julia> resample(spec, wave_sampled, LinearInterpolation; extrapolation = ExtrapolationType.Constant) -Spectrum(Int64, Float64) - -For better performance, see [Spectra.resample(spec::Spectrum, wave_sampled, interp::DataInterpolations.AbstractInterpolation)](@ref). -``` -""" -function Spectra.resample(spec::Spectrum, wave_sampled, interp::Type{<:AbstractInterpolation}; kws...) - itp = interp(flux(spec), wave(spec); kws...) - flux_sampled = itp(wave_sampled) - return spectrum(wave_sampled, flux_sampled; meta = spec.meta) -end - -end # module diff --git a/src/Spectra.jl b/src/Spectra.jl index 30ff569..90e9857 100644 --- a/src/Spectra.jl +++ b/src/Spectra.jl @@ -6,8 +6,8 @@ export AbstractSpectrum, spectrum export blackbody # transforms/redden.jl export redden, redden!, deredden, deredden! -# resampling: ../ext/DataInterpolationsExt.jl ../ext/Interpolations.jl -export resample +# transforms/resampler +export SpectrumResampler using RecipesBase: @recipe using Measurements: Measurements, Measurement @@ -21,6 +21,15 @@ include("common.jl") include("spectrum.jl") include("EchelleSpectrum.jl") +# Resampler type +include("transforms/resampler.jl") + +# tools +include("utils.jl") +include("transforms/transforms.jl") +include("plotting.jl") +include("fitting/fitting.jl") + """ spectrum(wave, flux; kwds...) @@ -93,87 +102,4 @@ function spectrum(wave::AbstractMatrix{<:Quantity}, flux::AbstractMatrix{<:Quant EchelleSpectrum(wave, flux, Dict{Symbol,Any}(kwds)) end -# Stub -""" - function resample(spec, wave_sampled, interp) - -Resample a spectrum onto the given wavelength grid `wave_sampled` using the supplied interpolator `interp`. Currently supports interpolators from the following packages: - - * [DataInterpolations.jl](https://github.com/SciML/DataInterpolations.jl) - * [Interpolations.jl](https://github.com/JuliaMath/Interpolations.jl) - -# Examples - -```julia-repl -julia> using DataInterpolations - -julia> spec = spectrum([2, 4, 12, 16, 20], [1, 3, 7, 6, 20]); - -julia> wave_sampled = [1, 5, 9, 13, 14, 17, 21, 22, 23]; - - # BYO interpolator -julia> interp = ... - -julia> spec_sampled = resample(spec, wave_sampled, interp) -Spectrum(Int64, Float64) -``` -""" -function resample(spec, wave_sampled, interp) - error("""No supported interpolation package loaded. To enable resampling, load one of the supported interpolation packages below: - - * [DataInterpolations.jl](https://github.com/SciML/DataInterpolations.jl) - * [Interpolations.jl](https://github.com/JuliaMath/Interpolations.jl) - - PRs to add additional interpolation packages are welcome! - """) -end - -""" - SpectrumResampler(spec::Spectrum, interp) -Type representing the spectrum `spec` with interpolator `interp`. - -```jldoctest -julia> using Spectra: spectrum, flux, wave, SpectrumResampler - -julia> using DataInterpolations: LinearInterpolation, ExtrapolationType - -julia> spec = spectrum([2, 4, 12, 16, 20], [1, 3, 7, 6, 20]); - -julia> wave_sampled = [1, 5, 9, 13, 14, 17, 21, 22, 23]; - -julia> interp = LinearInterpolation(flux(spec), wave(spec); extrapolation = ExtrapolationType.Constant); - -julia> resampler = SpectrumResampler(spec, interp); - -julia> result = resampler(wave_sampled); - -julia> result isa SpectrumResampler -true - -julia> wave(result) == wave_sampled -true - -julia> flux(result) == interp.(wave_sampled) -true -``` -""" -struct SpectrumResampler{A <: Spectrum, B} - spectrum::A - interp::B -end -wave(s::SpectrumResampler) = wave(s.spectrum) -flux(s::SpectrumResampler) = flux(s.spectrum) -spectrum(s::SpectrumResampler) = s.spectrum -function (spec::SpectrumResampler)(wave_sampled) - newspec = spectrum(wave_sampled, spec.interp.(wave_sampled); meta = spectrum(spec).meta) - return SpectrumResampler(newspec, spec.interp) -end -resample(spec::SpectrumResampler, wave_sampled) = spec(wave_sampled) - -# tools -include("utils.jl") -include("transforms/transforms.jl") -include("plotting.jl") -include("fitting/fitting.jl") - end # module diff --git a/src/transforms/resampler.jl b/src/transforms/resampler.jl new file mode 100644 index 0000000..453d534 --- /dev/null +++ b/src/transforms/resampler.jl @@ -0,0 +1,51 @@ +""" + SpectrumResampler(spec::Spectrum, interp) + +Type representing the spectrum `spec` with interpolator `interp`. + +```jldoctest +julia> using Spectra: spectrum, flux, wave, SpectrumResampler + +julia> using DataInterpolations: LinearInterpolation, ExtrapolationType + +julia> spec = spectrum([20, 40, 120, 160, 200], [1, 3, 7, 6, 20]); + +julia> interp = LinearInterpolation(flux(spec), wave(spec); extrapolation = ExtrapolationType.Constant); + +julia> resampler = SpectrumResampler(spec, interp); + +julia> wave_sampled = [10, 50, 90, 130, 140, 170, 210, 220, 230]; + +julia> result = resampler(wave_sampled); + +julia> result isa SpectrumResampler +true + +julia> wave(result) == wave_sampled +true + +julia> flux(result) == interp.(wave_sampled) +true +``` +""" +struct SpectrumResampler{A <: Spectrum, B} + spectrum::A + interp::B +end + +wave(s::SpectrumResampler) = wave(s.spectrum) +flux(s::SpectrumResampler) = flux(s.spectrum) +spectrum(s::SpectrumResampler) = s.spectrum + +function (s::SpectrumResampler)(wave_sampled) + interp = s.interp + spec_resampled = interp(wave_sampled) + s_new = spectrum(wave_sampled, spec_resampled; meta = spectrum(s).meta) + return SpectrumResampler(s_new, interp) +end + +function Base.show(io::IO, s::SpectrumResampler) + println(io, "SpectrumResampler($(eltype(wave(s))), $(eltype(flux(s))))") + println(io, " spec: ", spectrum(s)) + print(io, " interpolator: ", s.interp) +end From 921cbaa883cee259ef187b6cb28737b956f217d6 Mon Sep 17 00:00:00 2001 From: icweaver Date: Sun, 2 Nov 2025 20:48:47 -0800 Subject: [PATCH 09/25] update make.jl --- docs/make.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/docs/make.jl b/docs/make.jl index 950496e..dfee206 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -27,5 +27,7 @@ makedocs(sitename = "Spectra.jl", deploydocs(; repo = "github.com/JuliaAstro/Spectra.jl.git", + devbranch = "main", + push_preview = true, versions = ["stable" => "v^", "v#.#"] # Restrict to minor releases ) From 7d596265ecfe3b244364ec46ece581a80c2b3801 Mon Sep 17 00:00:00 2001 From: cgarling Date: Mon, 3 Nov 2025 08:06:57 -0500 Subject: [PATCH 10/25] Fix doctest on local --- docs/src/index.md | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/docs/src/index.md b/docs/src/index.md index 92b6a7f..6765f79 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -23,9 +23,9 @@ Here is a quick demo of some of our features ```jldoctest guide julia> using Spectra, FITSIO, Unitful, UnitfulAstro, Plots -julia> fitsurl = "https://dr14.sdss.org/optical/spectrum/view/data/format=fits/spec=lite?plateid=1323&mjd=52797&fiberid=12"; +julia> # fitsurl = "https://dr14.sdss.org/optical/spectrum/view/data/format=fits/spec=lite?plateid=1323&mjd=52797&fiberid=12"; -julia> # download(fitsurl, "sdss.fits"); +julia> # f = FITS(HTTP.get(fitsurl).body) julia> f = FITS("sdss.fits") File: sdss.fits @@ -49,10 +49,13 @@ julia> plot(spec); ![](assets/sdss.svg) ```jldoctest guide -julia> cont_fit = continuum(spec) -Spectrum(Quantity{Float32, ๐‹, Unitful.FreeUnits{(ร…,), ๐‹, nothing}}, Quantity{Float64, ๐Œ ๐‹^-1 ๐“^-3, Unitful.FreeUnits{(ร…^-1, erg, cm^-2, s^-1), ๐Œ ๐‹^-1 ๐“^-3, nothing}}) - coeffs: Quantity{Float64, ๐Œ ๐‹^-1 ๐“^-3, Unitful.FreeUnits{(ร…^-1, erg, cm^-2, s^-1), ๐Œ ๐‹^-1 ๐“^-3, nothing}}[1.983152216046405e-15 erg ร…^-1 cm^-2 s^-1, -1.8822245369267038e-16 erg ร…^-1 cm^-2 s^-1, -1.0422750370065006e-16 erg ร…^-1 cm^-2 s^-1, 4.8112282273206135e-17 erg ร…^-1 cm^-2 s^-1] - normalized: true +julia> cont_fit = continuum(spec); + +julia> cont_fit isa Spectra.Spectrum +true + +julia> cont_fit.normalized +true julia> plot(cont_fit, xlims=(6545, 6600)); ``` From 8a2dbc0185e36093432080d1bc26e5bb77051747 Mon Sep 17 00:00:00 2001 From: cgarling Date: Mon, 3 Nov 2025 08:38:55 -0500 Subject: [PATCH 11/25] Expand `SpectrumResampler` docstring --- src/transforms/resampler.jl | 43 ++++++++++++++++++++++++++++++++++--- 1 file changed, 40 insertions(+), 3 deletions(-) diff --git a/src/transforms/resampler.jl b/src/transforms/resampler.jl index 453d534..d96de8b 100644 --- a/src/transforms/resampler.jl +++ b/src/transforms/resampler.jl @@ -3,18 +3,55 @@ Type representing the spectrum `spec` with interpolator `interp`. -```jldoctest +Interpolation methods from many packages can be used without issue. Below we show example usage of [DataInterpolations.jl](https://github.com/SciML/DataInterpolations.jl) and [Interpolations.jl](https://github.com/JuliaMath/Interpolations.jl). + +First, we set up a an arbitrary spectrum and a linear interpolator from DataInterpolations.jl +```jldoctest resampling julia> using Spectra: spectrum, flux, wave, SpectrumResampler julia> using DataInterpolations: LinearInterpolation, ExtrapolationType julia> spec = spectrum([20, 40, 120, 160, 200], [1, 3, 7, 6, 20]); -julia> interp = LinearInterpolation(flux(spec), wave(spec); extrapolation = ExtrapolationType.Constant); +julia> interp = LinearInterpolation(flux(spec), wave(spec); + extrapolation = ExtrapolationType.Constant); +``` + +Now, we construct the `SpectrumResampler` and define the new wavelength grid that we want to resample the original spectrum to. +```jldoctest resampling julia> resampler = SpectrumResampler(spec, interp); julia> wave_sampled = [10, 50, 90, 130, 140, 170, 210, 220, 230]; +``` + +To perform the resampling, you call the resampler with the desired wavelength grid. + +```jldoctest resampling +julia> result = resampler(wave_sampled); +``` + +The resampled wavelength and flux can be obtained with the `wave` and `flux` methods. + +```jldoctest resampling +julia> result isa SpectrumResampler +true + +julia> wave(result) == wave_sampled +true + +julia> flux(result) == interp(wave_sampled) +true +``` + +Use of [Interpolations.jl](https://github.com/JuliaMath/Interpolations.jl) follows the same general procedure but using a different `interp`. + +```jldoctest resampling +julia> using Interpolations: linear_interpolation, Flat + +julia> interp = linear_interpolation(wave(spec), flux(spec); extrapolation_bc = Flat()); + +julia> resampler = SpectrumResampler(spec, interp); julia> result = resampler(wave_sampled); @@ -24,7 +61,7 @@ true julia> wave(result) == wave_sampled true -julia> flux(result) == interp.(wave_sampled) +julia> flux(result) == interp(wave_sampled) true ``` """ From bf1373a4c281114f3d42008212ffda0f69e20f81 Mon Sep 17 00:00:00 2001 From: cgarling Date: Mon, 3 Nov 2025 09:03:56 -0500 Subject: [PATCH 12/25] Use `Spectrum` constructor in resampler Saves allocations over using `spectrum` probably due how the `meta` dict is handled internally. --- src/transforms/resampler.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/transforms/resampler.jl b/src/transforms/resampler.jl index d96de8b..a7e3ade 100644 --- a/src/transforms/resampler.jl +++ b/src/transforms/resampler.jl @@ -77,7 +77,8 @@ spectrum(s::SpectrumResampler) = s.spectrum function (s::SpectrumResampler)(wave_sampled) interp = s.interp spec_resampled = interp(wave_sampled) - s_new = spectrum(wave_sampled, spec_resampled; meta = spectrum(s).meta) + s_new = Spectrum(wave_sampled, spec_resampled, spectrum(s).meta) + # s_new = spectrum(wave_sampled, spec_resampled; meta = spectrum(s).meta) return SpectrumResampler(s_new, interp) end From 32ea04300a21fffb7240acaf2a37041b936ab42c Mon Sep 17 00:00:00 2001 From: icweaver Date: Mon, 3 Nov 2025 09:49:33 -0800 Subject: [PATCH 13/25] point docs to docstring --- docs/src/transforms.md | 74 ++++--------------------------------- src/transforms/resampler.jl | 12 +++--- 2 files changed, 14 insertions(+), 72 deletions(-) diff --git a/docs/src/transforms.md b/docs/src/transforms.md index 5fa3a63..e745d54 100644 --- a/docs/src/transforms.md +++ b/docs/src/transforms.md @@ -56,73 +56,6 @@ julia> red.flux โ‰ˆ spec.flux true ``` -## Resampling - -External interpolators, e.g., from [DataInterpolations.jl](https://github.com/SciML/DataInterpolations.jl), [Interpolations.jl](https://github.com/JuliaMath/Interpolations.jl), can be used to resample spectra onto a given wavelength grid. Starting with a sample spectrum `spec`, we first create a [`SpectrumResampler`](@ref) object `resampler` which stores the initial spectrum and interpolator `interp` together. We then apply this object to the wavelength grid our our choice to produce our desired resampled spectrum. Below, we show how this may be accomplished using either package. - -### DataInterpolations.jl - -```@repl resample -using Spectra: SpectrumResampler, spectrum, wave, flux - -using DataInterpolations: LinearInterpolation, ExtrapolationType - -spec = spectrum([20, 40, 120, 160, 200], [1, 3, 7, 6, 20]) - -interp = LinearInterpolation(flux(spec), wave(spec); extrapolation = ExtrapolationType.Constant) - -resampler = SpectrumResampler(spec, interp) -``` - -We can now apply `resampler` to our desired wavelength grid `wave_sampled` to produce our resampled spectrum: - -```@repl resample -wave_sampled = [10, 50, 90, 130, 140, 170, 210, 220, 230]; - -spec_resampled = resampler(wave_sampled) -``` - -Calling this on different wavelength grids uses the same stored interpolator, allowing for efficient resampling. For convenience, the resampled spectrum is returned as another `SpectrumResampler` object, which keeps information about the resampled spectrum and applied interpolator coupled. The resulting interpolated flux and wavelength grid can be retrieved with the `flux` and `wave` getters exported by Spectra.jl, respectively: - -```@repl resample -wave(spec_resampled) - -flux(spec_resampled) - -# Check that this is equivalent to calling the interpolator directly -flux(spec_resampled) == interp(wave_sampled) -``` - -### Interpolations.jl - -For completeness, here is how a similar resampling scheme could be accomplished with Interpolations.jl: - -```@repl -using Spectra: SpectrumResampler, spectrum, wave, flux - -using Interpolations: linear_interpolation, Flat - -spec = spectrum([20, 40, 120, 160, 200], [1, 3, 7, 6, 20]) - -# Note that `wave` and `flux` are flipped relative to the DataInterpolations.jl example -interp = linear_interpolation(wave(spec), flux(spec); extrapolation_bc = Flat()) - -resampler = SpectrumResampler(spec, interp) - -wave_sampled = [10, 50, 90, 130, 140, 170, 210, 220, 230]; - -spec_resampled = resampler(wave_sampled) - -wave(spec_resampled) - -flux(spec_resampled) - -# Check that this is equivalent to calling the interpolator directly -flux(spec_resampled) == interp(wave_sampled) -``` - -See the relevant documentation for the external interpolation package used for more information on what configurations can be passed for your desired `interp` object. - ### API/Reference ```@docs @@ -130,6 +63,13 @@ redden redden! deredden deredden! +``` + +## Resampling + +External interpolators, e.g., from [DataInterpolations.jl](https://github.com/SciML/DataInterpolations.jl) or [Interpolations.jl](https://github.com/JuliaMath/Interpolations.jl), can be used to resample spectra onto a given wavelength grid. Starting with a sample spectrum `spec`, we first create a [`SpectrumResampler`](@ref) object `resampler` which stores the initial spectrum and interpolator `interp` together. We then apply this object to the wavelength grid of our choice to produce the resampled spectrum. We show example usage in the docstring below: + +```@docs SpectrumResampler ``` diff --git a/src/transforms/resampler.jl b/src/transforms/resampler.jl index a7e3ade..39eb094 100644 --- a/src/transforms/resampler.jl +++ b/src/transforms/resampler.jl @@ -5,7 +5,8 @@ Type representing the spectrum `spec` with interpolator `interp`. Interpolation methods from many packages can be used without issue. Below we show example usage of [DataInterpolations.jl](https://github.com/SciML/DataInterpolations.jl) and [Interpolations.jl](https://github.com/JuliaMath/Interpolations.jl). -First, we set up a an arbitrary spectrum and a linear interpolator from DataInterpolations.jl +First, we set up an arbitrary spectrum and a linear interpolator from DataInterpolations.jl: + ```jldoctest resampling julia> using Spectra: spectrum, flux, wave, SpectrumResampler @@ -13,11 +14,12 @@ julia> using DataInterpolations: LinearInterpolation, ExtrapolationType julia> spec = spectrum([20, 40, 120, 160, 200], [1, 3, 7, 6, 20]); -julia> interp = LinearInterpolation(flux(spec), wave(spec); - extrapolation = ExtrapolationType.Constant); +julia> interp = LinearInterpolation(flux(spec), wave(spec); + extrapolation = ExtrapolationType.Constant + ); ``` -Now, we construct the `SpectrumResampler` and define the new wavelength grid that we want to resample the original spectrum to. +Now, we construct the `SpectrumResampler` and define the new wavelength grid that we want to resample the original spectrum to: ```jldoctest resampling julia> resampler = SpectrumResampler(spec, interp); @@ -44,7 +46,7 @@ julia> flux(result) == interp(wave_sampled) true ``` -Use of [Interpolations.jl](https://github.com/JuliaMath/Interpolations.jl) follows the same general procedure but using a different `interp`. +Use of [Interpolations.jl](https://github.com/JuliaMath/Interpolations.jl) follows the same general procedure, but using a different `interp`: ```jldoctest resampling julia> using Interpolations: linear_interpolation, Flat From 2768c0f8c667a69aefb26fdacf8f8b4f8fffee57 Mon Sep 17 00:00:00 2001 From: icweaver Date: Mon, 3 Nov 2025 11:57:11 -0800 Subject: [PATCH 14/25] tests up --- test/Project.toml | 1 + test/transforms/resample.jl | 174 +++++++++++++++++++++--------------- 2 files changed, 103 insertions(+), 72 deletions(-) diff --git a/test/Project.toml b/test/Project.toml index 99cdd30..512392c 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,5 +1,6 @@ [deps] Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" +DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0" DustExtinction = "fb44c06c-c62f-5397-83f5-69249e0a3c8e" Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7" ParallelTestRunner = "d3525ed8-44d0-4b2c-a655-542cee43accc" diff --git a/test/transforms/resample.jl b/test/transforms/resample.jl index b5c219f..f20571b 100644 --- a/test/transforms/resample.jl +++ b/test/transforms/resample.jl @@ -1,72 +1,102 @@ -#function mock_spectrum(n::Int = Int(1e3); use_units::Bool = false) -# wave = range(1e4, 3e4, length = n) -# sigma = 0.1 .* sin.(wave) -# T = 6700 -# flux = @. 1e14 / (wave^5 * (exp(1 / (wave * T)) - 1)) ยฑ sigma -# if use_units -# wave *= u"angstrom" -# flux *= u"erg/s/cm^2/angstrom" -# end -# spectrum(wave, flux, name="Test Spectrum") -#end -# -#@testset "resample" begin -# @testset "Resampling" begin -# spec = mock_spectrum() -# new_wave = range(minimum(spec.wave), maximum(spec.wave), length=Integer(length(spec.wave) รท 2.4)) -# res_spec = resample(spec, new_wave) -# -# @test res_spec.wave == new_wave -# @test length(res_spec.flux) == length(new_wave) -# -# resample!(spec, new_wave) -# @test res_spec.wave == spec.wave -# @test res_spec.flux == spec.flux -# -# # Unitful -# spec = mock_spectrum(use_units=true) -# new_wave = range(minimum(spec.wave), maximum(spec.wave), length=Integer(length(spec.wave) รท 2.4)) -# @assert unit(eltype(new_wave)) == unit(eltype(spec.wave)) -# res_spec = resample(spec, new_wave) -# -# @test res_spec.wave == new_wave -# @test length(res_spec.flux) == length(new_wave) -# -# resample!(spec, new_wave) -# @test res_spec.wave == spec.wave -# @test res_spec.flux == spec.flux -# -# # Test resampling to another Spectrum -# spec1 = mock_spectrum(Integer(1e4)) -# spec2 = mock_spectrum(Integer(1e3)) -# res_spec = resample(spec1, spec2) -# @test res_spec.wave == spec2.wave -# -# resample!(spec1, spec2) -# @test spec1.wave == spec2.wave -# @test spec1.flux == res_spec.flux -# -# # Unitful -# spec1 = mock_spectrum(Integer(1e4), use_units=true) -# spec2 = mock_spectrum(Integer(1e3), use_units=true) -# res_spec = resample(spec1, spec2) -# -# @test res_spec.wave == spec2.wave -# @test unit(eltype(spec1.flux)) == unit(eltype(spec2.flux)) == unit(eltype(res_spec.flux)) -# -# # Test when spectra2 has different units -# spec1 = mock_spectrum(Integer(1e4), use_units=true) -# spec2 = mock_spectrum(Integer(1e3), use_units=true) -# spec1.wave = uconvert.(u"cm", spec1.wave) -# resample!(spec1, spec2) -# @test ustrip.(spec1.wave) โ‰ˆ ustrip.(unit(eltype(spec1.wave)), spec2.wave) -# -# # address bug when doing the same thing but affecting spec2 -# spec1 = mock_spectrum(Integer(1e4), use_units=true) -# spec2 = mock_spectrum(Integer(1e3), use_units=true) -# spec2.wave = uconvert.(u"cm", spec2.wave) -# resample!(spec1, spec2) -# @test ustrip.(spec1.wave) โ‰ˆ ustrip.(unit(eltype(spec1.wave)), spec2.wave) -# -# end -#end +using Spectra: Spectra, AbstractSpectrum, SpectrumResampler, spectrum +using DataInterpolations: LinearInterpolation, ExtrapolationType +using Unitful: @u_str, uconvert +using UnitfulAstro +using Measurements: ยฑ + +# TODO: See if it makes sense to have an exported API for resample/resample! even though we are using external interpolators. + +function resample(spec, new_wave) + interp = LinearInterpolation(Spectra.flux(spec), Spectra.wave(spec); extrapolation = ExtrapolationType.Constant) + resampler = SpectrumResampler(spec, interp) + return resampler(new_wave) +end + +function resample(spec1::AbstractSpectrum, spec2::AbstractSpectrum) + new_wave = Spectra.wave(spec2) + return resample(spec1, new_wave) +end + +function resample!(spec::AbstractSpectrum, new_wave) + spec_new = resample(spec, new_wave) + spec.flux = Spectra.flux(spec_new) + spec.wave = Spectra.wave(spec_new) + return spec +end + +function resample!(spec1::AbstractSpectrum, spec2::AbstractSpectrum) + new_wave = Spectra.wave(spec2) + resample!(spec1, new_wave) +end + +function mock_spectrum(n::Int = Int(1e3); use_units::Bool = false) + wave = range(1e4, 3e4, length = n) + sigma = 0.1 .* sin.(wave) + T = 6700 + flux = @. 1e14 / (wave^5 * (exp(1 / (wave * T)) - 1)) ยฑ sigma + if use_units + wave *= u"angstrom" + flux *= u"erg/s/cm^2/angstrom" + end + spectrum(wave, flux; name = "Test Spectrum") +end + +@testset "resample" begin + @testset "Resampling" begin + spec = mock_spectrum() + new_wave = range(minimum(spec.wave), maximum(spec.wave); length = Integer(length(spec.wave) รท 2.4)) + res_spec = resample(spec, new_wave) + + @test Spectra.wave(res_spec) == new_wave + @test length(Spectra.flux(res_spec)) == length(new_wave) + + resample!(spec, new_wave) + @test Spectra.wave(res_spec) == spec.wave + @test Spectra.flux(res_spec) == spec.flux + + # Unitful + spec = mock_spectrum(use_units=true) + new_wave = range(minimum(spec.wave), maximum(spec.wave), length=Integer(length(spec.wave) รท 2.4)) + @assert unit(eltype(new_wave)) == unit(eltype(spec.wave)) + res_spec = resample(spec, new_wave) + + @test Spectra.wave(res_spec) == new_wave + @test length(Spectra.flux(res_spec)) == length(new_wave) + + resample!(spec, new_wave) + @test Spectra.wave(res_spec) == spec.wave + @test Spectra.flux(res_spec) == spec.flux + + # Test resampling to another Spectrum + spec1 = mock_spectrum(Integer(1e4)) + spec2 = mock_spectrum(Integer(1e3)) + res_spec = resample(spec1, spec2) + @test Spectra.wave(res_spec) == Spectra.wave(spec2) + + resample!(spec1, spec2) + @test Spectra.wave(spec1) == spec2.wave + @test Spectra.flux(spec1) == Spectra.flux(res_spec) + + # Unitful + spec1 = mock_spectrum(Integer(1e4), use_units=true) + spec2 = mock_spectrum(Integer(1e3), use_units=true) + res_spec = resample(spec1, spec2) + + @test Spectra.wave(res_spec) == spec2.wave + @test unit(eltype(Spectra.flux(spec1))) == unit(eltype(Spectra.flux(spec2))) == unit(eltype(Spectra.flux(res_spec))) + + # Test when spectra2 has different units + spec1 = mock_spectrum(Integer(1e4), use_units=true) + spec2 = mock_spectrum(Integer(1e3), use_units=true) + spec1.wave = uconvert.(u"cm", spec1.wave) + resample!(spec1, spec2) + @test ustrip.(Spectra.wave(spec1)) โ‰ˆ ustrip.(unit(eltype(Spectra.wave(spec1))), Spectra.wave(spec2)) + + # address bug when doing the same thing but affecting spec2 + spec1 = mock_spectrum(Integer(1e4), use_units=true) + spec2 = mock_spectrum(Integer(1e3), use_units=true) + spec2.wave = uconvert.(u"cm", spec2.wave) + resample!(spec1, spec2) + @test ustrip.(Spectra.wave(spec1)) โ‰ˆ ustrip.(unit(eltype(Spectra.wave(spec1))), Spectra.wave(spec2)) + end +end From 17bb315db770b64cfe5dd21cf64b30b12def60d6 Mon Sep 17 00:00:00 2001 From: icweaver Date: Mon, 3 Nov 2025 12:43:54 -0800 Subject: [PATCH 15/25] added ext by mistake --- ext/InterpolationsExt.jl | 80 ---------------------------------------- 1 file changed, 80 deletions(-) delete mode 100644 ext/InterpolationsExt.jl diff --git a/ext/InterpolationsExt.jl b/ext/InterpolationsExt.jl deleted file mode 100644 index b43ec9a..0000000 --- a/ext/InterpolationsExt.jl +++ /dev/null @@ -1,80 +0,0 @@ -module InterpolationsExt - -using Spectra: Spectra, Spectrum, spectrum, wave, flux -import Spectra: resample -using Interpolations: - Interpolations, - AbstractInterpolation, - Degree, - Constant, - Linear, - Cubic, - constant_interpolation, - linear_interpolation, - cubic_interpolation - -""" - Spectra.resample(spec::Spectrum, wave_sampled, interp::AbstractInterpolation) - -Resample `spec` using the given interpolator `interp`. - -```jldoctest -julia> using Spectra: spectrum, resample, wave, flux - -julia> using Interpolations - -julia> spec = spectrum([2, 4, 12, 16, 20], [1, 3, 7, 6, 20]); - -julia> wave_sampled = [1, 5, 9, 13, 14, 17, 21, 22, 23]; - -julia> interp = extrapolate(interpolate((wave(spec), ), flux(spec), Gridded(Linear())), Flat()); - -julia> resample(spec, wave_sampled, interp) -Spectrum(Int64, Float64) -``` - -See [resample(spec::Spectrum, wave_sampled, deg::Interpolations.Degree; kws...)](@ref) for a convenience function version with limited flexibility. -""" -function Spectra.resample(spec::Spectrum, wave_sampled, interp::AbstractInterpolation) - flux_sampled = interp(wave_sampled) - return spectrum(wave_sampled, flux_sampled; meta = spec.meta) -end - -""" - resample(spec::Spectrum, wave_sampled, deg::Interpolations.Degree; kws...) - -Constructs an interpolator of type `interp` from the provided spectrum `spec` and evaluates it at `wave_sampled`, returning a `Spectrum`. The `kws...` control extrapolation and are passed through to the interpolator. - -```jldoctest -julia> using Spectra: spectrum, resample - -julia> using Interpolations: Flat - -julia> spec = spectrum([2, 4, 12, 16, 20], [1, 3, 7, 6, 20]); - -julia> wave_sampled = [1, 5, 9, 13, 14, 17, 21, 22, 23]; - -julia> resample(spec, wave_sampled, Linear(); extrapolation_bc = Flat()) -Spectrum(Int64, Float64) -``` - -This convenience function currently only supports `Constant`, `Linear`, and `Cubic` splines from Interpolations.jl. For more interpolators and better performance, see [Spectra.resample(spec::Spectrum, wave_sampled, interp::Interpolations.AbstractInterpolation)](@ref). -""" -function Spectra.resample(spec::Spectrum, wave_sampled, deg::Degree; kws...) - - deg_func = if deg isa Constant - constant_interpolation - elseif deg isa Linear - linear_interpolation - elseif deg isa Cubic - cubic_interpolation - else - throw(ArgumentError(deg, "Invalid degree passed. Must be one of `Constant`, `Linear`, `Quadratic`, or `Cubic`.")) - end - - itp = deg_func(wave(spec), flux(spec); kws...) - flux_sampled = itp(wave_sampled) - return spectrum(wave_sampled, flux_sampled; meta = spec.meta) -end - -end # module From 46dfee47702c5b896dd05648e8e6fffd1b55e769 Mon Sep 17 00:00:00 2001 From: icweaver Date: Tue, 4 Nov 2025 10:46:16 -0800 Subject: [PATCH 16/25] updated show method + test --- src/transforms/resampler.jl | 2 +- test/transforms/resample.jl | 6 ++++++ 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/src/transforms/resampler.jl b/src/transforms/resampler.jl index 39eb094..2886957 100644 --- a/src/transforms/resampler.jl +++ b/src/transforms/resampler.jl @@ -87,5 +87,5 @@ end function Base.show(io::IO, s::SpectrumResampler) println(io, "SpectrumResampler($(eltype(wave(s))), $(eltype(flux(s))))") println(io, " spec: ", spectrum(s)) - print(io, " interpolator: ", s.interp) + print(io, " interpolator: ", typeof(s.interp)) end diff --git a/test/transforms/resample.jl b/test/transforms/resample.jl index f20571b..28478f1 100644 --- a/test/transforms/resample.jl +++ b/test/transforms/resample.jl @@ -46,7 +46,13 @@ end spec = mock_spectrum() new_wave = range(minimum(spec.wave), maximum(spec.wave); length = Integer(length(spec.wave) รท 2.4)) res_spec = resample(spec, new_wave) + expected = """ + SpectrumResampler(Float64, Measurements.Measurement{Float64}) + spec: Spectrum(Float64, Measurements.Measurement{Float64}) + name: Test Spectrum + interpolator: DataInterpolations.LinearInterpolation{Vector{Measurements.Measurement{Float64}}, Vector{Float64}, Vector{Measurements.Measurement{Float64}}, Vector{Measurements.Measurement{Float64}}, Measurements.Measurement{Float64}}""" + @test sprint(show, res_spec) == expected @test Spectra.wave(res_spec) == new_wave @test length(Spectra.flux(res_spec)) == length(new_wave) From e2ceff8dc50f785ada7799e4d58a1b61f190d53d Mon Sep 17 00:00:00 2001 From: icweaver Date: Tue, 4 Nov 2025 18:24:42 -0800 Subject: [PATCH 17/25] drop mutation tests, using immutable structs in #35 --- test/transforms/resample.jl | 60 ++++++++++++++++++------------------- 1 file changed, 30 insertions(+), 30 deletions(-) diff --git a/test/transforms/resample.jl b/test/transforms/resample.jl index 28478f1..2d3545d 100644 --- a/test/transforms/resample.jl +++ b/test/transforms/resample.jl @@ -17,17 +17,17 @@ function resample(spec1::AbstractSpectrum, spec2::AbstractSpectrum) return resample(spec1, new_wave) end -function resample!(spec::AbstractSpectrum, new_wave) - spec_new = resample(spec, new_wave) - spec.flux = Spectra.flux(spec_new) - spec.wave = Spectra.wave(spec_new) - return spec -end - -function resample!(spec1::AbstractSpectrum, spec2::AbstractSpectrum) - new_wave = Spectra.wave(spec2) - resample!(spec1, new_wave) -end +#function resample!(spec::AbstractSpectrum, new_wave) +# spec_new = resample(spec, new_wave) +# spec.flux = Spectra.flux(spec_new) +# spec.wave = Spectra.wave(spec_new) +# return spec +#end + +#function resample!(spec1::AbstractSpectrum, spec2::AbstractSpectrum) +# new_wave = Spectra.wave(spec2) +# resample!(spec1, new_wave) +#end function mock_spectrum(n::Int = Int(1e3); use_units::Bool = false) wave = range(1e4, 3e4, length = n) @@ -56,9 +56,9 @@ end @test Spectra.wave(res_spec) == new_wave @test length(Spectra.flux(res_spec)) == length(new_wave) - resample!(spec, new_wave) - @test Spectra.wave(res_spec) == spec.wave - @test Spectra.flux(res_spec) == spec.flux + #resample!(spec, new_wave) + #@test Spectra.wave(res_spec) == spec.wave + #@test Spectra.flux(res_spec) == spec.flux # Unitful spec = mock_spectrum(use_units=true) @@ -69,9 +69,9 @@ end @test Spectra.wave(res_spec) == new_wave @test length(Spectra.flux(res_spec)) == length(new_wave) - resample!(spec, new_wave) - @test Spectra.wave(res_spec) == spec.wave - @test Spectra.flux(res_spec) == spec.flux + #resample!(spec, new_wave) + #@test Spectra.wave(res_spec) == spec.wave + #@test Spectra.flux(res_spec) == spec.flux # Test resampling to another Spectrum spec1 = mock_spectrum(Integer(1e4)) @@ -79,9 +79,9 @@ end res_spec = resample(spec1, spec2) @test Spectra.wave(res_spec) == Spectra.wave(spec2) - resample!(spec1, spec2) - @test Spectra.wave(spec1) == spec2.wave - @test Spectra.flux(spec1) == Spectra.flux(res_spec) + #resample!(spec1, spec2) + #@test Spectra.wave(spec1) == spec2.wave + #@test Spectra.flux(spec1) == Spectra.flux(res_spec) # Unitful spec1 = mock_spectrum(Integer(1e4), use_units=true) @@ -92,17 +92,17 @@ end @test unit(eltype(Spectra.flux(spec1))) == unit(eltype(Spectra.flux(spec2))) == unit(eltype(Spectra.flux(res_spec))) # Test when spectra2 has different units - spec1 = mock_spectrum(Integer(1e4), use_units=true) - spec2 = mock_spectrum(Integer(1e3), use_units=true) - spec1.wave = uconvert.(u"cm", spec1.wave) - resample!(spec1, spec2) - @test ustrip.(Spectra.wave(spec1)) โ‰ˆ ustrip.(unit(eltype(Spectra.wave(spec1))), Spectra.wave(spec2)) + #spec1 = mock_spectrum(Integer(1e4), use_units=true) + #spec2 = mock_spectrum(Integer(1e3), use_units=true) + #spec1.wave = uconvert.(u"cm", spec1.wave) + #resample!(spec1, spec2) + #@test ustrip.(Spectra.wave(spec1)) โ‰ˆ ustrip.(unit(eltype(Spectra.wave(spec1))), Spectra.wave(spec2)) # address bug when doing the same thing but affecting spec2 - spec1 = mock_spectrum(Integer(1e4), use_units=true) - spec2 = mock_spectrum(Integer(1e3), use_units=true) - spec2.wave = uconvert.(u"cm", spec2.wave) - resample!(spec1, spec2) - @test ustrip.(Spectra.wave(spec1)) โ‰ˆ ustrip.(unit(eltype(Spectra.wave(spec1))), Spectra.wave(spec2)) + #spec1 = mock_spectrum(Integer(1e4), use_units=true) + #spec2 = mock_spectrum(Integer(1e3), use_units=true) + #spec2.wave = uconvert.(u"cm", spec2.wave) + #resample!(spec1, spec2) + #@test ustrip.(Spectra.wave(spec1)) โ‰ˆ ustrip.(unit(eltype(Spectra.wave(spec1))), Spectra.wave(spec2)) end end From 2de9ca3178dd79dec3dad52a641a2081f85338b3 Mon Sep 17 00:00:00 2001 From: cgarling Date: Tue, 4 Nov 2025 21:56:15 -0500 Subject: [PATCH 18/25] Simplify docs --- docs/src/transforms.md | 15 ++------------- src/transforms/resampler.jl | 3 +-- 2 files changed, 3 insertions(+), 15 deletions(-) diff --git a/docs/src/transforms.md b/docs/src/transforms.md index e745d54..a8dc048 100644 --- a/docs/src/transforms.md +++ b/docs/src/transforms.md @@ -1,18 +1,11 @@ -```@meta -DocTestSetup = quote - using Spectra, Random - Random.seed!(11894) -end -``` - # Transformations ## Extinction -By levaraging [DustExtinction.jl](https://github.com/juliaastro/dustextinction.jl) we can apply common reddening laws to our spectra. +By leveraging [DustExtinction.jl](https://github.com/juliaastro/dustextinction.jl) we can apply common reddening laws to our spectra. ```jldoctest -julia> using Unitful, Measurements, Random +julia> using Spectra, Unitful, Measurements, Random julia> rng = Random.seed!(0); @@ -72,7 +65,3 @@ External interpolators, e.g., from [DataInterpolations.jl](https://github.com/Sc ```@docs SpectrumResampler ``` - -```@meta -DocTestSetup = nothing -``` diff --git a/src/transforms/resampler.jl b/src/transforms/resampler.jl index 2886957..ad08c6e 100644 --- a/src/transforms/resampler.jl +++ b/src/transforms/resampler.jl @@ -15,8 +15,7 @@ julia> using DataInterpolations: LinearInterpolation, ExtrapolationType julia> spec = spectrum([20, 40, 120, 160, 200], [1, 3, 7, 6, 20]); julia> interp = LinearInterpolation(flux(spec), wave(spec); - extrapolation = ExtrapolationType.Constant - ); + extrapolation = ExtrapolationType.Constant); ``` Now, we construct the `SpectrumResampler` and define the new wavelength grid that we want to resample the original spectrum to: From 44e236a942fc68d20b2535d651097bc08faa59e1 Mon Sep 17 00:00:00 2001 From: cgarling Date: Tue, 4 Nov 2025 22:22:57 -0500 Subject: [PATCH 19/25] use getters --- src/transforms/resampler.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/transforms/resampler.jl b/src/transforms/resampler.jl index ad08c6e..ec92d0d 100644 --- a/src/transforms/resampler.jl +++ b/src/transforms/resampler.jl @@ -71,9 +71,9 @@ struct SpectrumResampler{A <: Spectrum, B} interp::B end -wave(s::SpectrumResampler) = wave(s.spectrum) -flux(s::SpectrumResampler) = flux(s.spectrum) spectrum(s::SpectrumResampler) = s.spectrum +wave(s::SpectrumResampler) = wave(spectrum(s)) +flux(s::SpectrumResampler) = flux(spectrum(s)) function (s::SpectrumResampler)(wave_sampled) interp = s.interp From 5b15e7bdb3a95d8305bc5213cd398d5b5e1ce823 Mon Sep 17 00:00:00 2001 From: Ian Weaver Date: Thu, 6 Nov 2025 18:09:10 -0800 Subject: [PATCH 20/25] just return a regular Spectrum --- src/transforms/resampler.jl | 7 ++----- test/transforms/resample.jl | 6 ++---- 2 files changed, 4 insertions(+), 9 deletions(-) diff --git a/src/transforms/resampler.jl b/src/transforms/resampler.jl index ec92d0d..852a69c 100644 --- a/src/transforms/resampler.jl +++ b/src/transforms/resampler.jl @@ -76,11 +76,8 @@ wave(s::SpectrumResampler) = wave(spectrum(s)) flux(s::SpectrumResampler) = flux(spectrum(s)) function (s::SpectrumResampler)(wave_sampled) - interp = s.interp - spec_resampled = interp(wave_sampled) - s_new = Spectrum(wave_sampled, spec_resampled, spectrum(s).meta) - # s_new = spectrum(wave_sampled, spec_resampled; meta = spectrum(s).meta) - return SpectrumResampler(s_new, interp) + flux_resampled = (s.interp)(wave_sampled) + return Spectrum(wave_sampled, flux_resampled, spectrum(s).meta) end function Base.show(io::IO, s::SpectrumResampler) diff --git a/test/transforms/resample.jl b/test/transforms/resample.jl index 2d3545d..43c09bf 100644 --- a/test/transforms/resample.jl +++ b/test/transforms/resample.jl @@ -47,10 +47,8 @@ end new_wave = range(minimum(spec.wave), maximum(spec.wave); length = Integer(length(spec.wave) รท 2.4)) res_spec = resample(spec, new_wave) expected = """ - SpectrumResampler(Float64, Measurements.Measurement{Float64}) - spec: Spectrum(Float64, Measurements.Measurement{Float64}) - name: Test Spectrum - interpolator: DataInterpolations.LinearInterpolation{Vector{Measurements.Measurement{Float64}}, Vector{Float64}, Vector{Measurements.Measurement{Float64}}, Vector{Measurements.Measurement{Float64}}, Measurements.Measurement{Float64}}""" + Spectrum(Float64, Measurements.Measurement{Float64}) + name: Test Spectrum""" @test sprint(show, res_spec) == expected @test Spectra.wave(res_spec) == new_wave From 54fbc98f1839e0fee75ecbe9bc33f110eb01f47d Mon Sep 17 00:00:00 2001 From: Ian Weaver Date: Fri, 7 Nov 2025 04:23:13 -0800 Subject: [PATCH 21/25] more test cov --- test/transforms/resample.jl | 149 ++++++++++++++++++++---------------- 1 file changed, 81 insertions(+), 68 deletions(-) diff --git a/test/transforms/resample.jl b/test/transforms/resample.jl index 43c09bf..5b83134 100644 --- a/test/transforms/resample.jl +++ b/test/transforms/resample.jl @@ -1,4 +1,4 @@ -using Spectra: Spectra, AbstractSpectrum, SpectrumResampler, spectrum +using Spectra: Spectra, AbstractSpectrum, SpectrumResampler, spectrum, wave, flux using DataInterpolations: LinearInterpolation, ExtrapolationType using Unitful: @u_str, uconvert using UnitfulAstro @@ -7,25 +7,25 @@ using Measurements: ยฑ # TODO: See if it makes sense to have an exported API for resample/resample! even though we are using external interpolators. function resample(spec, new_wave) - interp = LinearInterpolation(Spectra.flux(spec), Spectra.wave(spec); extrapolation = ExtrapolationType.Constant) + interp = LinearInterpolation(flux(spec), wave(spec); extrapolation = ExtrapolationType.Constant) resampler = SpectrumResampler(spec, interp) return resampler(new_wave) end function resample(spec1::AbstractSpectrum, spec2::AbstractSpectrum) - new_wave = Spectra.wave(spec2) + new_wave = wave(spec2) return resample(spec1, new_wave) end #function resample!(spec::AbstractSpectrum, new_wave) # spec_new = resample(spec, new_wave) -# spec.flux = Spectra.flux(spec_new) -# spec.wave = Spectra.wave(spec_new) +# spec.flux = flux(spec_new) +# spec.wave = wave(spec_new) # return spec #end #function resample!(spec1::AbstractSpectrum, spec2::AbstractSpectrum) -# new_wave = Spectra.wave(spec2) +# new_wave = wave(spec2) # resample!(spec1, new_wave) #end @@ -41,66 +41,79 @@ function mock_spectrum(n::Int = Int(1e3); use_units::Bool = false) spectrum(wave, flux; name = "Test Spectrum") end -@testset "resample" begin - @testset "Resampling" begin - spec = mock_spectrum() - new_wave = range(minimum(spec.wave), maximum(spec.wave); length = Integer(length(spec.wave) รท 2.4)) - res_spec = resample(spec, new_wave) - expected = """ - Spectrum(Float64, Measurements.Measurement{Float64}) - name: Test Spectrum""" - - @test sprint(show, res_spec) == expected - @test Spectra.wave(res_spec) == new_wave - @test length(Spectra.flux(res_spec)) == length(new_wave) - - #resample!(spec, new_wave) - #@test Spectra.wave(res_spec) == spec.wave - #@test Spectra.flux(res_spec) == spec.flux - - # Unitful - spec = mock_spectrum(use_units=true) - new_wave = range(minimum(spec.wave), maximum(spec.wave), length=Integer(length(spec.wave) รท 2.4)) - @assert unit(eltype(new_wave)) == unit(eltype(spec.wave)) - res_spec = resample(spec, new_wave) - - @test Spectra.wave(res_spec) == new_wave - @test length(Spectra.flux(res_spec)) == length(new_wave) - - #resample!(spec, new_wave) - #@test Spectra.wave(res_spec) == spec.wave - #@test Spectra.flux(res_spec) == spec.flux - - # Test resampling to another Spectrum - spec1 = mock_spectrum(Integer(1e4)) - spec2 = mock_spectrum(Integer(1e3)) - res_spec = resample(spec1, spec2) - @test Spectra.wave(res_spec) == Spectra.wave(spec2) - - #resample!(spec1, spec2) - #@test Spectra.wave(spec1) == spec2.wave - #@test Spectra.flux(spec1) == Spectra.flux(res_spec) - - # Unitful - spec1 = mock_spectrum(Integer(1e4), use_units=true) - spec2 = mock_spectrum(Integer(1e3), use_units=true) - res_spec = resample(spec1, spec2) - - @test Spectra.wave(res_spec) == spec2.wave - @test unit(eltype(Spectra.flux(spec1))) == unit(eltype(Spectra.flux(spec2))) == unit(eltype(Spectra.flux(res_spec))) - - # Test when spectra2 has different units - #spec1 = mock_spectrum(Integer(1e4), use_units=true) - #spec2 = mock_spectrum(Integer(1e3), use_units=true) - #spec1.wave = uconvert.(u"cm", spec1.wave) - #resample!(spec1, spec2) - #@test ustrip.(Spectra.wave(spec1)) โ‰ˆ ustrip.(unit(eltype(Spectra.wave(spec1))), Spectra.wave(spec2)) - - # address bug when doing the same thing but affecting spec2 - #spec1 = mock_spectrum(Integer(1e4), use_units=true) - #spec2 = mock_spectrum(Integer(1e3), use_units=true) - #spec2.wave = uconvert.(u"cm", spec2.wave) - #resample!(spec1, spec2) - #@test ustrip.(Spectra.wave(spec1)) โ‰ˆ ustrip.(unit(eltype(Spectra.wave(spec1))), Spectra.wave(spec2)) - end +@testset "Resampler" begin + spec = mock_spectrum() + interp = LinearInterpolation(flux(spec), wave(spec); extrapolation = ExtrapolationType.Constant) + resampler = SpectrumResampler(spec, interp) + expected = """ + SpectrumResampler(Float64, Measurements.Measurement{Float64}) + spec: Spectrum(Float64, Measurements.Measurement{Float64}) + name: Test Spectrum + interpolator: LinearInterpolation{Vector{Measurements.Measurement{Float64}}, Vector{Float64}, Vector{Measurements.Measurement{Float64}}, Vector{Measurements.Measurement{Float64}}, Measurements.Measurement{Float64}}""" + + @test sprint(show, resampler) == expected + @test wave(resampler) == wave(spec) + @test flux(resampler) == flux(spec) +end + +@testset "Resampling" begin + spec = mock_spectrum() + new_wave = range(minimum(spec.wave), maximum(spec.wave); length = Integer(length(spec.wave) รท 2.4)) + res_spec = resample(spec, new_wave) + expected = """ + Spectrum(Float64, Measurements.Measurement{Float64}) + name: Test Spectrum""" + + @test sprint(show, res_spec) == expected + @test wave(res_spec) == new_wave + @test length(flux(res_spec)) == length(new_wave) + + #resample!(spec, new_wave) + #@test wave(res_spec) == spec.wave + #@test flux(res_spec) == spec.flux + + # Unitful + spec = mock_spectrum(use_units=true) + new_wave = range(minimum(spec.wave), maximum(spec.wave), length=Integer(length(spec.wave) รท 2.4)) + @assert unit(eltype(new_wave)) == unit(eltype(spec.wave)) + res_spec = resample(spec, new_wave) + + @test wave(res_spec) == new_wave + @test length(flux(res_spec)) == length(new_wave) + + #resample!(spec, new_wave) + #@test wave(res_spec) == spec.wave + #@test flux(res_spec) == spec.flux + + # Test resampling to another Spectrum + spec1 = mock_spectrum(Integer(1e4)) + spec2 = mock_spectrum(Integer(1e3)) + res_spec = resample(spec1, spec2) + @test wave(res_spec) == wave(spec2) + + #resample!(spec1, spec2) + #@test wave(spec1) == spec2.wave + #@test flux(spec1) == flux(res_spec) + + # Unitful + spec1 = mock_spectrum(Integer(1e4), use_units=true) + spec2 = mock_spectrum(Integer(1e3), use_units=true) + res_spec = resample(spec1, spec2) + + @test wave(res_spec) == spec2.wave + @test unit(eltype(flux(spec1))) == unit(eltype(flux(spec2))) == unit(eltype(flux(res_spec))) + + # Test when spectra2 has different units + #spec1 = mock_spectrum(Integer(1e4), use_units=true) + #spec2 = mock_spectrum(Integer(1e3), use_units=true) + #spec1.wave = uconvert.(u"cm", spec1.wave) + #resample!(spec1, spec2) + #@test ustrip.(wave(spec1)) โ‰ˆ ustrip.(unit(eltype(wave(spec1))), wave(spec2)) + + # address bug when doing the same thing but affecting spec2 + #spec1 = mock_spectrum(Integer(1e4), use_units=true) + #spec2 = mock_spectrum(Integer(1e3), use_units=true) + #spec2.wave = uconvert.(u"cm", spec2.wave) + #resample!(spec1, spec2) + #@test ustrip.(wave(spec1)) โ‰ˆ ustrip.(unit(eltype(wave(spec1))), wave(spec2)) end From c9a03fa0b8331c1d49c004f0d8d946781e1e9ae8 Mon Sep 17 00:00:00 2001 From: Ian Weaver Date: Fri, 7 Nov 2025 04:28:03 -0800 Subject: [PATCH 22/25] update tests --- src/Spectra.jl | 2 +- src/transforms/resampler.jl | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/Spectra.jl b/src/Spectra.jl index 90e9857..9a21918 100644 --- a/src/Spectra.jl +++ b/src/Spectra.jl @@ -1,7 +1,7 @@ module Spectra # common.jl -export AbstractSpectrum, spectrum +export AbstractSpectrum, Spectrum, spectrum # utils.jl export blackbody # transforms/redden.jl diff --git a/src/transforms/resampler.jl b/src/transforms/resampler.jl index 852a69c..07f1d2f 100644 --- a/src/transforms/resampler.jl +++ b/src/transforms/resampler.jl @@ -35,7 +35,7 @@ julia> result = resampler(wave_sampled); The resampled wavelength and flux can be obtained with the `wave` and `flux` methods. ```jldoctest resampling -julia> result isa SpectrumResampler +julia> result isa Spectrum true julia> wave(result) == wave_sampled @@ -56,7 +56,7 @@ julia> resampler = SpectrumResampler(spec, interp); julia> result = resampler(wave_sampled); -julia> result isa SpectrumResampler +julia> result isa Spectrum true julia> wave(result) == wave_sampled From 6f022f5582bd27b1aa0557d0d6d05e1d5fccf273 Mon Sep 17 00:00:00 2001 From: Ian Weaver Date: Fri, 7 Nov 2025 04:36:29 -0800 Subject: [PATCH 23/25] qualify interp --- test/transforms/resample.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/transforms/resample.jl b/test/transforms/resample.jl index 5b83134..6f1846b 100644 --- a/test/transforms/resample.jl +++ b/test/transforms/resample.jl @@ -49,7 +49,7 @@ end SpectrumResampler(Float64, Measurements.Measurement{Float64}) spec: Spectrum(Float64, Measurements.Measurement{Float64}) name: Test Spectrum - interpolator: LinearInterpolation{Vector{Measurements.Measurement{Float64}}, Vector{Float64}, Vector{Measurements.Measurement{Float64}}, Vector{Measurements.Measurement{Float64}}, Measurements.Measurement{Float64}}""" + interpolator: DataInterpolations.LinearInterpolation{Vector{Measurements.Measurement{Float64}}, Vector{Float64}, Vector{Measurements.Measurement{Float64}}, Vector{Measurements.Measurement{Float64}}, Measurements.Measurement{Float64}}""" @test sprint(show, resampler) == expected @test wave(resampler) == wave(spec) From 68683ff207f14c3ceaa6a51d01814bfa7d918e19 Mon Sep 17 00:00:00 2001 From: Ian Weaver Date: Thu, 18 Dec 2025 15:49:56 -0800 Subject: [PATCH 24/25] cleanup --- src/Spectra.jl | 31 ++++++-------- src/transforms/resampler.jl | 11 ++--- test/transforms/resample.jl | 80 ++++++++++++++----------------------- 3 files changed, 48 insertions(+), 74 deletions(-) diff --git a/src/Spectra.jl b/src/Spectra.jl index 209277d..7dbab5a 100644 --- a/src/Spectra.jl +++ b/src/Spectra.jl @@ -2,16 +2,18 @@ module Spectra # Uniform API export AbstractSpectrum, Spectrum, spectrum, spectral_axis, flux_axis -# spectra_single.jl, spectra_ifu.jl, spectra_echelle.jl, spectra_binned.jl -export SingleSpectrum, IFUSpectrum, EchelleSpectrum, spectral_axis, flux_axis -# utils.jl + +# AbstractSpectrum types +export SingleSpectrum, IFUSpectrum, EchelleSpectrum + +# Transforms +export SpectrumResampler, redden, redden!, deredden, deredden! + +# Utilities export blackbody #, line_flux, equivalent_width -# fitting/fitting.jl + +# Fitting #export continuum, continuum! -# transforms/redden.jl -export redden, redden!, deredden, deredden! -# transforms/resampler -export SpectrumResampler using RecipesBase: @recipe using Measurements: Measurements, Measurement @@ -208,15 +210,6 @@ include("spectrum_echelle.jl") include("spectrum_ifu.jl") #include("spectrum_binned.jl") -# Resampler type -include("transforms/resampler.jl") - -# tools -include("utils.jl") -include("transforms/transforms.jl") -include("plotting.jl") -include("fitting/fitting.jl") - """ spectrum(spectral_axis, flux_axis, [meta]) @@ -305,10 +298,10 @@ function spectrum(spectral_axis::AbstractMatrix{<:Quantity}, flux_axis::Abstract Spectrum(spectral_axis, flux_axis, Dict{Symbol,Any}(kwds)) end -# tools include("utils.jl") +include("transforms/resampler.jl") include("transforms/transforms.jl") include("plotting.jl") #include("fitting/fitting.jl") -end # module \ No newline at end of file +end # module diff --git a/src/transforms/resampler.jl b/src/transforms/resampler.jl index 07f1d2f..9f7be6f 100644 --- a/src/transforms/resampler.jl +++ b/src/transforms/resampler.jl @@ -72,16 +72,17 @@ struct SpectrumResampler{A <: Spectrum, B} end spectrum(s::SpectrumResampler) = s.spectrum -wave(s::SpectrumResampler) = wave(spectrum(s)) -flux(s::SpectrumResampler) = flux(spectrum(s)) +spectral_axis(s::SpectrumResampler) = (spectral_axis โˆ˜ spectrum)(s) +flux_axis(s::SpectrumResampler) = (flux_axis โˆ˜ spectrum)(s) +meta(s::SpectrumResampler) = (meta โˆ˜ spectrum)(s) function (s::SpectrumResampler)(wave_sampled) flux_resampled = (s.interp)(wave_sampled) - return Spectrum(wave_sampled, flux_resampled, spectrum(s).meta) + return Spectrum(wave_sampled, flux_resampled, meta(s)) end function Base.show(io::IO, s::SpectrumResampler) - println(io, "SpectrumResampler($(eltype(wave(s))), $(eltype(flux(s))))") - println(io, " spec: ", spectrum(s)) + println(io, "SpectrumResampler(", (eltype โˆ˜ spectral_axis)(s), ", ", (eltype โˆ˜ flux_axis)(s), ")") + println(io, " spec: ", (typeof โˆ˜ spectrum)(s)) print(io, " interpolator: ", typeof(s.interp)) end diff --git a/test/transforms/resample.jl b/test/transforms/resample.jl index 6f1846b..1a599f1 100644 --- a/test/transforms/resample.jl +++ b/test/transforms/resample.jl @@ -1,4 +1,4 @@ -using Spectra: Spectra, AbstractSpectrum, SpectrumResampler, spectrum, wave, flux +using Spectra: Spectra, AbstractSpectrum, SpectrumResampler, spectrum, spectral_axis, flux_axis using DataInterpolations: LinearInterpolation, ExtrapolationType using Unitful: @u_str, uconvert using UnitfulAstro @@ -7,28 +7,16 @@ using Measurements: ยฑ # TODO: See if it makes sense to have an exported API for resample/resample! even though we are using external interpolators. function resample(spec, new_wave) - interp = LinearInterpolation(flux(spec), wave(spec); extrapolation = ExtrapolationType.Constant) + interp = LinearInterpolation(flux_axis(spec), spectral_axis(spec); extrapolation = ExtrapolationType.Constant) resampler = SpectrumResampler(spec, interp) return resampler(new_wave) end function resample(spec1::AbstractSpectrum, spec2::AbstractSpectrum) - new_wave = wave(spec2) + new_wave = spectral_axis(spec2) return resample(spec1, new_wave) end -#function resample!(spec::AbstractSpectrum, new_wave) -# spec_new = resample(spec, new_wave) -# spec.flux = flux(spec_new) -# spec.wave = wave(spec_new) -# return spec -#end - -#function resample!(spec1::AbstractSpectrum, spec2::AbstractSpectrum) -# new_wave = wave(spec2) -# resample!(spec1, new_wave) -#end - function mock_spectrum(n::Int = Int(1e3); use_units::Bool = false) wave = range(1e4, 3e4, length = n) sigma = 0.1 .* sin.(wave) @@ -43,77 +31,69 @@ end @testset "Resampler" begin spec = mock_spectrum() - interp = LinearInterpolation(flux(spec), wave(spec); extrapolation = ExtrapolationType.Constant) + s, f = spectral_axis(spec), flux_axis(spec) + interp = LinearInterpolation(f, s; extrapolation = ExtrapolationType.Constant) resampler = SpectrumResampler(spec, interp) expected = """ SpectrumResampler(Float64, Measurements.Measurement{Float64}) - spec: Spectrum(Float64, Measurements.Measurement{Float64}) - name: Test Spectrum - interpolator: DataInterpolations.LinearInterpolation{Vector{Measurements.Measurement{Float64}}, Vector{Float64}, Vector{Measurements.Measurement{Float64}}, Vector{Measurements.Measurement{Float64}}, Measurements.Measurement{Float64}}""" + spec: Spectra.SingleSpectrum{Float64, Measurements.Measurement{Float64}} + interpolator: DataInterpolations.LinearInterpolation{Vector{Measurements.Measurement{Float64}}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Vector{Measurements.Measurement{Float64}}, Vector{Measurements.Measurement{Float64}}, Measurements.Measurement{Float64}}""" @test sprint(show, resampler) == expected - @test wave(resampler) == wave(spec) - @test flux(resampler) == flux(spec) + @test spectral_axis(resampler) == s + @test flux_axis(resampler) == f end @testset "Resampling" begin spec = mock_spectrum() - new_wave = range(minimum(spec.wave), maximum(spec.wave); length = Integer(length(spec.wave) รท 2.4)) + s, f = spectral_axis(spec), flux_axis(spec) + new_wave = range(minimum(s), maximum(s); length = Integer(length(s) รท 2.4)) res_spec = resample(spec, new_wave) expected = """ - Spectrum(Float64, Measurements.Measurement{Float64}) - name: Test Spectrum""" + SingleSpectrum(Float64, Measurements.Measurement{Float64}) + spectral axis (416,): 10000.0 .. 30000.0 + flux axis (416,): 67.0 ยฑ 0.031 .. 0.827 ยฑ 0.08 + meta: Dict{Symbol, Any}(:name => "Test Spectrum")""" @test sprint(show, res_spec) == expected - @test wave(res_spec) == new_wave - @test length(flux(res_spec)) == length(new_wave) - - #resample!(spec, new_wave) - #@test wave(res_spec) == spec.wave - #@test flux(res_spec) == spec.flux + @test spectral_axis(res_spec) == new_wave + @test length(flux_axis(res_spec)) == length(new_wave) # Unitful - spec = mock_spectrum(use_units=true) - new_wave = range(minimum(spec.wave), maximum(spec.wave), length=Integer(length(spec.wave) รท 2.4)) - @assert unit(eltype(new_wave)) == unit(eltype(spec.wave)) + spec = mock_spectrum(; use_units = true) + s, f = spectral_axis(spec), flux_axis(spec) + new_wave = range(minimum(s), maximum(s); length = Integer(length(s) รท 2.4)) + @assert unit(eltype(new_wave)) == unit(eltype(s)) res_spec = resample(spec, new_wave) - @test wave(res_spec) == new_wave - @test length(flux(res_spec)) == length(new_wave) - - #resample!(spec, new_wave) - #@test wave(res_spec) == spec.wave - #@test flux(res_spec) == spec.flux + @test spectral_axis(res_spec) == new_wave + @test length(flux_axis(res_spec)) == length(new_wave) # Test resampling to another Spectrum spec1 = mock_spectrum(Integer(1e4)) spec2 = mock_spectrum(Integer(1e3)) res_spec = resample(spec1, spec2) - @test wave(res_spec) == wave(spec2) - - #resample!(spec1, spec2) - #@test wave(spec1) == spec2.wave - #@test flux(spec1) == flux(res_spec) + @test spectral_axis(res_spec) == spectral_axis(spec2) # Unitful - spec1 = mock_spectrum(Integer(1e4), use_units=true) - spec2 = mock_spectrum(Integer(1e3), use_units=true) + spec1 = mock_spectrum(Integer(1e4); use_units = true) + spec2 = mock_spectrum(Integer(1e3); use_units = true) res_spec = resample(spec1, spec2) - @test wave(res_spec) == spec2.wave - @test unit(eltype(flux(spec1))) == unit(eltype(flux(spec2))) == unit(eltype(flux(res_spec))) + @test spectral_axis(res_spec) == spectral_axis(spec2) + @test unit(eltype(flux_axis(spec1))) == unit(eltype(flux_axis(spec2))) == unit(eltype(flux_axis(res_spec))) # Test when spectra2 has different units #spec1 = mock_spectrum(Integer(1e4), use_units=true) #spec2 = mock_spectrum(Integer(1e3), use_units=true) #spec1.wave = uconvert.(u"cm", spec1.wave) #resample!(spec1, spec2) - #@test ustrip.(wave(spec1)) โ‰ˆ ustrip.(unit(eltype(wave(spec1))), wave(spec2)) + #@test ustrip.(spectral_axis(spec1)) โ‰ˆ ustrip.(unit(eltype(spectral_axis(spec1))), spectral_axis(spec2)) # address bug when doing the same thing but affecting spec2 #spec1 = mock_spectrum(Integer(1e4), use_units=true) #spec2 = mock_spectrum(Integer(1e3), use_units=true) #spec2.wave = uconvert.(u"cm", spec2.wave) #resample!(spec1, spec2) - #@test ustrip.(wave(spec1)) โ‰ˆ ustrip.(unit(eltype(wave(spec1))), wave(spec2)) + #@test ustrip.(spectral_axis(spec1)) โ‰ˆ ustrip.(unit(eltype(spectral_axis(spec1))), spectral_axis(spec2)) end From fab8c54a4a7707ae90e2a5b53c34ccb65227e9df Mon Sep 17 00:00:00 2001 From: Ian Weaver Date: Thu, 18 Dec 2025 19:43:26 -0800 Subject: [PATCH 25/25] cleanup --- src/Spectra.jl | 2 +- src/transforms/resampler.jl | 14 +++++++------- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/Spectra.jl b/src/Spectra.jl index 7dbab5a..4c61508 100644 --- a/src/Spectra.jl +++ b/src/Spectra.jl @@ -86,7 +86,7 @@ Return the spectral axis of `spec`. spectral_axis(spec::AbstractSpectrum) = spec.spectral_axis """ - flux(spec::AbstractSpectrum) + flux_axis(spec::AbstractSpectrum) Return the flux axis of `spec`. """ diff --git a/src/transforms/resampler.jl b/src/transforms/resampler.jl index 9f7be6f..1de68ad 100644 --- a/src/transforms/resampler.jl +++ b/src/transforms/resampler.jl @@ -8,13 +8,13 @@ Interpolation methods from many packages can be used without issue. Below we sho First, we set up an arbitrary spectrum and a linear interpolator from DataInterpolations.jl: ```jldoctest resampling -julia> using Spectra: spectrum, flux, wave, SpectrumResampler +julia> using Spectra: SpectrumResampler, spectrum, spectral_axis, flux_axis julia> using DataInterpolations: LinearInterpolation, ExtrapolationType julia> spec = spectrum([20, 40, 120, 160, 200], [1, 3, 7, 6, 20]); -julia> interp = LinearInterpolation(flux(spec), wave(spec); +julia> interp = LinearInterpolation(flux_axis(spec), spectral_axis(spec); extrapolation = ExtrapolationType.Constant); ``` @@ -38,10 +38,10 @@ The resampled wavelength and flux can be obtained with the `wave` and `flux` met julia> result isa Spectrum true -julia> wave(result) == wave_sampled +julia> spectral_axis(result) == wave_sampled true -julia> flux(result) == interp(wave_sampled) +julia> flux_axis(result) == interp(wave_sampled) true ``` @@ -50,7 +50,7 @@ Use of [Interpolations.jl](https://github.com/JuliaMath/Interpolations.jl) follo ```jldoctest resampling julia> using Interpolations: linear_interpolation, Flat -julia> interp = linear_interpolation(wave(spec), flux(spec); extrapolation_bc = Flat()); +julia> interp = linear_interpolation(spectral_axis(spec), flux_axis(spec); extrapolation_bc = Flat()); julia> resampler = SpectrumResampler(spec, interp); @@ -59,10 +59,10 @@ julia> result = resampler(wave_sampled); julia> result isa Spectrum true -julia> wave(result) == wave_sampled +julia> spectral_axis(result) == wave_sampled true -julia> flux(result) == interp(wave_sampled) +julia> flux_axis(result) == interp(wave_sampled) true ``` """