|
1 | | -#function mock_spectrum(n::Int = Int(1e3); use_units::Bool = false) |
2 | | -# wave = range(1e4, 3e4, length = n) |
3 | | -# sigma = 0.1 .* sin.(wave) |
4 | | -# T = 6700 |
5 | | -# flux = @. 1e14 / (wave^5 * (exp(1 / (wave * T)) - 1)) ± sigma |
6 | | -# if use_units |
7 | | -# wave *= u"angstrom" |
8 | | -# flux *= u"erg/s/cm^2/angstrom" |
9 | | -# end |
10 | | -# spectrum(wave, flux, name="Test Spectrum") |
11 | | -#end |
12 | | -# |
13 | | -#@testset "resample" begin |
14 | | -# @testset "Resampling" begin |
15 | | -# spec = mock_spectrum() |
16 | | -# new_wave = range(minimum(spec.wave), maximum(spec.wave), length=Integer(length(spec.wave) ÷ 2.4)) |
17 | | -# res_spec = resample(spec, new_wave) |
18 | | -# |
19 | | -# @test res_spec.wave == new_wave |
20 | | -# @test length(res_spec.flux) == length(new_wave) |
21 | | -# |
22 | | -# resample!(spec, new_wave) |
23 | | -# @test res_spec.wave == spec.wave |
24 | | -# @test res_spec.flux == spec.flux |
25 | | -# |
26 | | -# # Unitful |
27 | | -# spec = mock_spectrum(use_units=true) |
28 | | -# new_wave = range(minimum(spec.wave), maximum(spec.wave), length=Integer(length(spec.wave) ÷ 2.4)) |
29 | | -# @assert unit(eltype(new_wave)) == unit(eltype(spec.wave)) |
30 | | -# res_spec = resample(spec, new_wave) |
31 | | -# |
32 | | -# @test res_spec.wave == new_wave |
33 | | -# @test length(res_spec.flux) == length(new_wave) |
34 | | -# |
35 | | -# resample!(spec, new_wave) |
36 | | -# @test res_spec.wave == spec.wave |
37 | | -# @test res_spec.flux == spec.flux |
38 | | -# |
39 | | -# # Test resampling to another Spectrum |
40 | | -# spec1 = mock_spectrum(Integer(1e4)) |
41 | | -# spec2 = mock_spectrum(Integer(1e3)) |
42 | | -# res_spec = resample(spec1, spec2) |
43 | | -# @test res_spec.wave == spec2.wave |
44 | | -# |
45 | | -# resample!(spec1, spec2) |
46 | | -# @test spec1.wave == spec2.wave |
47 | | -# @test spec1.flux == res_spec.flux |
48 | | -# |
49 | | -# # Unitful |
50 | | -# spec1 = mock_spectrum(Integer(1e4), use_units=true) |
51 | | -# spec2 = mock_spectrum(Integer(1e3), use_units=true) |
52 | | -# res_spec = resample(spec1, spec2) |
53 | | -# |
54 | | -# @test res_spec.wave == spec2.wave |
55 | | -# @test unit(eltype(spec1.flux)) == unit(eltype(spec2.flux)) == unit(eltype(res_spec.flux)) |
56 | | -# |
57 | | -# # Test when spectra2 has different units |
58 | | -# spec1 = mock_spectrum(Integer(1e4), use_units=true) |
59 | | -# spec2 = mock_spectrum(Integer(1e3), use_units=true) |
60 | | -# spec1.wave = uconvert.(u"cm", spec1.wave) |
61 | | -# resample!(spec1, spec2) |
62 | | -# @test ustrip.(spec1.wave) ≈ ustrip.(unit(eltype(spec1.wave)), spec2.wave) |
63 | | -# |
64 | | -# # address bug when doing the same thing but affecting spec2 |
65 | | -# spec1 = mock_spectrum(Integer(1e4), use_units=true) |
66 | | -# spec2 = mock_spectrum(Integer(1e3), use_units=true) |
67 | | -# spec2.wave = uconvert.(u"cm", spec2.wave) |
68 | | -# resample!(spec1, spec2) |
69 | | -# @test ustrip.(spec1.wave) ≈ ustrip.(unit(eltype(spec1.wave)), spec2.wave) |
70 | | -# |
71 | | -# end |
72 | | -#end |
| 1 | +using Spectra: Spectra, AbstractSpectrum, SpectrumResampler, spectrum, spectral_axis, flux_axis |
| 2 | +using DataInterpolations: LinearInterpolation, ExtrapolationType |
| 3 | +using Unitful: @u_str, uconvert |
| 4 | +using UnitfulAstro |
| 5 | +using Measurements: ± |
| 6 | + |
| 7 | +# TODO: See if it makes sense to have an exported API for resample/resample! even though we are using external interpolators. |
| 8 | + |
| 9 | +function resample(spec, new_wave) |
| 10 | + interp = LinearInterpolation(flux_axis(spec), spectral_axis(spec); extrapolation = ExtrapolationType.Constant) |
| 11 | + resampler = SpectrumResampler(spec, interp) |
| 12 | + return resampler(new_wave) |
| 13 | +end |
| 14 | + |
| 15 | +function resample(spec1::AbstractSpectrum, spec2::AbstractSpectrum) |
| 16 | + new_wave = spectral_axis(spec2) |
| 17 | + return resample(spec1, new_wave) |
| 18 | +end |
| 19 | + |
| 20 | +function mock_spectrum(n::Int = Int(1e3); use_units::Bool = false) |
| 21 | + wave = range(1e4, 3e4, length = n) |
| 22 | + sigma = 0.1 .* sin.(wave) |
| 23 | + T = 6700 |
| 24 | + flux = @. 1e14 / (wave^5 * (exp(1 / (wave * T)) - 1)) ± sigma |
| 25 | + if use_units |
| 26 | + wave *= u"angstrom" |
| 27 | + flux *= u"erg/s/cm^2/angstrom" |
| 28 | + end |
| 29 | + spectrum(wave, flux; name = "Test Spectrum") |
| 30 | +end |
| 31 | + |
| 32 | +@testset "Resampler" begin |
| 33 | + spec = mock_spectrum() |
| 34 | + s, f = spectral_axis(spec), flux_axis(spec) |
| 35 | + interp = LinearInterpolation(f, s; extrapolation = ExtrapolationType.Constant) |
| 36 | + resampler = SpectrumResampler(spec, interp) |
| 37 | + expected = """ |
| 38 | + SpectrumResampler(Float64, Measurements.Measurement{Float64}) |
| 39 | + spec: Spectra.SingleSpectrum{Float64, Measurements.Measurement{Float64}} |
| 40 | + 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}}""" |
| 41 | + |
| 42 | + @test sprint(show, resampler) == expected |
| 43 | + @test spectral_axis(resampler) == s |
| 44 | + @test flux_axis(resampler) == f |
| 45 | +end |
| 46 | + |
| 47 | +@testset "Resampling" begin |
| 48 | + spec = mock_spectrum() |
| 49 | + s, f = spectral_axis(spec), flux_axis(spec) |
| 50 | + new_wave = range(minimum(s), maximum(s); length = Integer(length(s) ÷ 2.4)) |
| 51 | + res_spec = resample(spec, new_wave) |
| 52 | + expected = """ |
| 53 | + SingleSpectrum(Float64, Measurements.Measurement{Float64}) |
| 54 | + spectral axis (416,): 10000.0 .. 30000.0 |
| 55 | + flux axis (416,): 67.0 ± 0.031 .. 0.827 ± 0.08 |
| 56 | + meta: Dict{Symbol, Any}(:name => "Test Spectrum")""" |
| 57 | + |
| 58 | + @test sprint(show, res_spec) == expected |
| 59 | + @test spectral_axis(res_spec) == new_wave |
| 60 | + @test length(flux_axis(res_spec)) == length(new_wave) |
| 61 | + |
| 62 | + # Unitful |
| 63 | + spec = mock_spectrum(; use_units = true) |
| 64 | + s, f = spectral_axis(spec), flux_axis(spec) |
| 65 | + new_wave = range(minimum(s), maximum(s); length = Integer(length(s) ÷ 2.4)) |
| 66 | + @assert unit(eltype(new_wave)) == unit(eltype(s)) |
| 67 | + res_spec = resample(spec, new_wave) |
| 68 | + |
| 69 | + @test spectral_axis(res_spec) == new_wave |
| 70 | + @test length(flux_axis(res_spec)) == length(new_wave) |
| 71 | + |
| 72 | + # Test resampling to another Spectrum |
| 73 | + spec1 = mock_spectrum(Integer(1e4)) |
| 74 | + spec2 = mock_spectrum(Integer(1e3)) |
| 75 | + res_spec = resample(spec1, spec2) |
| 76 | + @test spectral_axis(res_spec) == spectral_axis(spec2) |
| 77 | + |
| 78 | + # Unitful |
| 79 | + spec1 = mock_spectrum(Integer(1e4); use_units = true) |
| 80 | + spec2 = mock_spectrum(Integer(1e3); use_units = true) |
| 81 | + res_spec = resample(spec1, spec2) |
| 82 | + |
| 83 | + @test spectral_axis(res_spec) == spectral_axis(spec2) |
| 84 | + @test unit(eltype(flux_axis(spec1))) == unit(eltype(flux_axis(spec2))) == unit(eltype(flux_axis(res_spec))) |
| 85 | + |
| 86 | + # Test when spectra2 has different units |
| 87 | + #spec1 = mock_spectrum(Integer(1e4), use_units=true) |
| 88 | + #spec2 = mock_spectrum(Integer(1e3), use_units=true) |
| 89 | + #spec1.wave = uconvert.(u"cm", spec1.wave) |
| 90 | + #resample!(spec1, spec2) |
| 91 | + #@test ustrip.(spectral_axis(spec1)) ≈ ustrip.(unit(eltype(spectral_axis(spec1))), spectral_axis(spec2)) |
| 92 | + |
| 93 | + # address bug when doing the same thing but affecting spec2 |
| 94 | + #spec1 = mock_spectrum(Integer(1e4), use_units=true) |
| 95 | + #spec2 = mock_spectrum(Integer(1e3), use_units=true) |
| 96 | + #spec2.wave = uconvert.(u"cm", spec2.wave) |
| 97 | + #resample!(spec1, spec2) |
| 98 | + #@test ustrip.(spectral_axis(spec1)) ≈ ustrip.(unit(eltype(spectral_axis(spec1))), spectral_axis(spec2)) |
| 99 | +end |
0 commit comments