From 9ef02ce90879bb6bbd672ee740393ec70630a1be Mon Sep 17 00:00:00 2001 From: Michael Palumbo Date: Thu, 3 Feb 2022 13:29:23 -0500 Subject: [PATCH 1/5] quick fix for broken build bc of dead url --- deps/download_soap_example.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/deps/download_soap_example.jl b/deps/download_soap_example.jl index bfe282b..97a48b3 100644 --- a/deps/download_soap_example.jl +++ b/deps/download_soap_example.jl @@ -3,7 +3,7 @@ include("download.jl") -download_url = "https://psu.box.com/shared/static/nikva44r7xai9h0ulfj2pxewzv8ez4fo.h5" +download_url = "https://pennstateoffice365-my.sharepoint.com/:u:/g/personal/ebf11_psu_edu/EYyYhOHEjM5RqvZn1D-4I58BDbpbcYg6QZKhd0KQ3N7YoQ?e=Is4PcG" download_filename = joinpath("..","data","spectra","soap_demo.h5") download_md5 = "18a6291039d7684ef9187762e4b2dfd5" From 12898bec9a0b63bd0baad0526c06d32e56d2b481 Mon Sep 17 00:00:00 2001 From: Michael Palumbo Date: Tue, 8 Feb 2022 11:48:27 -0500 Subject: [PATCH 2/5] new quad regression --- src/calc_rv/fit_quadratic_to_ccf.jl | 46 ++++++++++++----------------- 1 file changed, 19 insertions(+), 27 deletions(-) diff --git a/src/calc_rv/fit_quadratic_to_ccf.jl b/src/calc_rv/fit_quadratic_to_ccf.jl index b5f34c2..a16001f 100644 --- a/src/calc_rv/fit_quadratic_to_ccf.jl +++ b/src/calc_rv/fit_quadratic_to_ccf.jl @@ -6,7 +6,7 @@ Based on code by Alex Wise (aw@psu.edu) Refactors & optimized by Eric Ford """ - +import PyPlot; plt = PyPlot; mpl = plt.matplotlib; plt.ioff() """ Functor to estimate RV based on fitting quadratic near minimum of CCF. TODO: Revist the logic here and see if need to perform transformation first. @@ -25,7 +25,7 @@ Optional Arguments: function MeasureRvFromCCFQuadratic(; frac_of_width_to_fit::Real = default_frac_of_width_to_fit, measure_width_at_frac_depth::Real = default_measure_width_at_frac_depth ) @assert 0.25 <= measure_width_at_frac_depth <= 0.75 - @assert 0.1 <= frac_of_width_to_fit <= 2.0 + @assert 0.0 <= frac_of_width_to_fit <= 2.0 MeasureRvFromCCFQuadratic(frac_of_width_to_fit,measure_width_at_frac_depth) end @@ -33,24 +33,24 @@ using Statistics function (mrv::MeasureRvFromCCFQuadratic)(vels::A1, ccf::A2 ) where {T1<:Real, A1<:AbstractArray{T1,1}, T2<:Real, A2<:AbstractArray{T2,1} } # find the min and fit only the part near the minimum of the CCF amin, inds = find_idx_at_and_around_minimum(vels, ccf, frac_of_width_to_fit=mrv.frac_of_width_to_fit, measure_width_at_frac_depth=mrv.measure_width_at_frac_depth) - if isnan(amin) return ( rv=NaN, σ_rv=NaN ) end + if isnan(amin) return (rv=NaN, σ_rv=NaN) end - mean_v = mean(view(vels,inds)) - X = ones(length(inds),3) - X[:,2] .= view(vels,inds) .-mean_v - X[:,3] .= (view(vels,inds) .-mean_v).^2 - (c, b, a) = (X'*X) \ (X'*view(ccf,inds)) - #= - # do the polyfit - pfit = Polynomials.fit(vels[inds], ccf[inds], 2) - @assert length(Polynomials.coeffs(pfit)) >= 3 # just in case fails to fit a quadratic + # set up linear system to solve + A = ones(3,3) + A[[3,5,7]] .= sum(view(vels, inds).^2) # set anti-diagonal indices + A[1,1] = length(view(vels, inds)) + A[1,2] = A[2,1] = sum(view(vels, inds)) + A[2,3] = A[3,2] = sum(view(vels, inds).^3) + A[3,3] = sum(view(vels, inds).^4) + + rhs = ones(3) + rhs[1] = sum(view(ccf,inds)) + rhs[2] = sum(view(ccf,inds) .* view(vels,inds)) + rhs[3] = sum(view(ccf,inds) .* view(vels,inds).^2) + (c, b, a) = A \ rhs - # get center from coeffs - c, b, a = Polynomials.coeffs(pfit) - @assert a>0 - =# - v_at_min_of_quadratic = -b/(2*a) + mean_v - return ( rv=v_at_min_of_quadratic, σ_rv=NaN ) + v_at_min_of_quadratic = -b/(2*a) + return (rv=v_at_min_of_quadratic, σ_rv=NaN) end function (mrv::MeasureRvFromCCFQuadratic)(vels::A1, ccf::A2, ccf_var::A2 ) where {T1<:Real, A1<:AbstractArray{T1,1}, T2<:Real, A2<:AbstractArray{T2,1}, T3<:Real, A3<:AbstractArray{T3,1} } @@ -68,13 +68,5 @@ function (mrv::MeasureRvFromCCFQuadratic)(vels::A1, ccf::A2, ccf_var::A2 ) where v_at_min_of_quadratic = -b/(2*a) sigma_rv = abs(v_at_min_of_quadratic) * sqrt(covar_beta[2,2]/b^2+covar_beta[3,3]/a^2) v_at_min_of_quadratic += mean_v - #= - # do the polyfit - pfit = Polynomials.fit(vels[inds], ccf[inds], 2) - @assert length(Polynomials.coeffs(pfit)) >= 3 # just in case fails to fit a quadratic - - # get center from coeffs - c, b, a = Polynomials.coeffs(pfit) - =# - return ( rv=v_at_min_of_quadratic, σ_rv=sigma_rv ) + return (rv=v_at_min_of_quadratic, σ_rv=sigma_rv) end From 41178bcb28b6edaa67224f2a93aeecc67129dc45 Mon Sep 17 00:00:00 2001 From: Michael Palumbo Date: Thu, 10 Feb 2022 08:50:53 -0500 Subject: [PATCH 3/5] fix download, added rv injection/extraction test --- deps/download_soap_example.jl | 3 +- src/calc_rv/fit_quadratic_to_ccf.jl | 4 +-- test/Project.toml | 1 + test/calc_rv.jl | 50 ++++++++++++++++++++++++++++- 4 files changed, 53 insertions(+), 5 deletions(-) diff --git a/deps/download_soap_example.jl b/deps/download_soap_example.jl index 97a48b3..646176d 100644 --- a/deps/download_soap_example.jl +++ b/deps/download_soap_example.jl @@ -3,7 +3,8 @@ include("download.jl") -download_url = "https://pennstateoffice365-my.sharepoint.com/:u:/g/personal/ebf11_psu_edu/EYyYhOHEjM5RqvZn1D-4I58BDbpbcYg6QZKhd0KQ3N7YoQ?e=Is4PcG" +# download_url = "https://pennstateoffice365-my.sharepoint.com/:u:/g/personal/ebf11_psu_edu/EYyYhOHEjM5RqvZn1D-4I58BDbpbcYg6QZKhd0KQ3N7YoQ?e=Is4PcG" +download_url = "https://pennstateoffice365-my.sharepoint.com/:u:/g/personal/ebf11_psu_edu/EYyYhOHEjM5RqvZn1D-4I58BDbpbcYg6QZKhd0KQ3N7YoQ?download=1" download_filename = joinpath("..","data","spectra","soap_demo.h5") download_md5 = "18a6291039d7684ef9187762e4b2dfd5" diff --git a/src/calc_rv/fit_quadratic_to_ccf.jl b/src/calc_rv/fit_quadratic_to_ccf.jl index a16001f..c661611 100644 --- a/src/calc_rv/fit_quadratic_to_ccf.jl +++ b/src/calc_rv/fit_quadratic_to_ccf.jl @@ -6,8 +6,6 @@ Based on code by Alex Wise (aw@psu.edu) Refactors & optimized by Eric Ford """ -import PyPlot; plt = PyPlot; mpl = plt.matplotlib; plt.ioff() - """ Functor to estimate RV based on fitting quadratic near minimum of CCF. TODO: Revist the logic here and see if need to perform transformation first. """ @@ -25,7 +23,7 @@ Optional Arguments: function MeasureRvFromCCFQuadratic(; frac_of_width_to_fit::Real = default_frac_of_width_to_fit, measure_width_at_frac_depth::Real = default_measure_width_at_frac_depth ) @assert 0.25 <= measure_width_at_frac_depth <= 0.75 - @assert 0.0 <= frac_of_width_to_fit <= 2.0 + @assert 0.1 <= frac_of_width_to_fit <= 2.0 MeasureRvFromCCFQuadratic(frac_of_width_to_fit,measure_width_at_frac_depth) end diff --git a/test/Project.toml b/test/Project.toml index 61fef04..e629e19 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -2,4 +2,5 @@ DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45" RvSpectMLBase = "c48404b2-35ea-40e7-ac7f-06a53de703d6" +Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/test/calc_rv.jl b/test/calc_rv.jl index fbab2ac..e5e801e 100644 --- a/test/calc_rv.jl +++ b/test/calc_rv.jl @@ -1,6 +1,6 @@ using EchelleCCFs using Test - +import Statistics: mean @testset "Check CCF accuracy" begin # TODO repeat for multiple mask shapes @testset "Tophat mask w/ 1 line" begin @@ -135,4 +135,52 @@ using Test @test vfit.rv ≈ 0 atol = 500 # TODO: tune tolerance better =# end + @testset "RV injection/extraction" begin + # velocities + amp = 1.0 # meters per second + vel = amp .* sin.(range(0, 2π, length=100)) + mid = 5434.5 + dep = 0.75 + wid = 0.1 + mid_shift = mid .* EchelleCCFs.calc_doppler_factor.(vel) + + # compute gaussian line spectrum + buffer = 0.75 + resolution = 7e5 + Δlnλ = (1.0 / resolution) + wavs = exp.(range(log(mid - buffer), log(mid + buffer), step=Δlnλ)) + + flux = zeros(length(wavs), length(vel)) + for i in eachindex(vel) + flux[:,i] .= @. 1.0 - dep * exp(-((wavs-mid_shift[i])/wid)^2 / 2.0) + end + + # measure ccf + line_list = EchelleCCFs.BasicLineList([mid], [dep]) + speed_of_light = EchelleCCFs.speed_of_light_mps + mask_width = speed_of_light/resolution + mask_shape = TopHatCCFMask(mask_width) + ccf_plan = BasicCCFPlan(line_list=line_list, mask_shape=mask_shape) + v_grid = EchelleCCFs.calc_ccf_v_grid(ccf_plan) + ccf = EchelleCCFs.ccf_1D(wavs, flux, ccf_plan) + + # measure velocities from CCF + mrv_gauss = MeasureRvFromCCFGaussian(frac_of_width_to_fit=0.5) + mrv_quad = MeasureRvFromCCFQuadratic(frac_of_width_to_fit=0.1) + + rvs_gauss = zeros(length(vel)) + sig_gauss = zeros(length(vel)) + rvs_quad = zeros(length(vel)) + sig_quad = zeros(length(vel)) + for i in eachindex(vel) + rvs_gauss[i], sig_gauss[i] = mrv_gauss(v_grid, view(ccf, :, i)) + rvs_quad[i], sig_quad[i] = mrv_quad(v_grid, view(ccf, :, i)) + end + err_gauss = abs.(rvs_gauss .- mean(rvs_gauss) .- vel) + err_quad = abs.(rvs_quad .- mean(rvs_quad) .- vel) + + # test maximum errors are less than a cm/s + @assert maximum(err_gauss) < 1e-2 + @assert maximum(err_quad) < 1e-2 + end end From ada684d0e89ea10c760adc99ccb9209691df9539 Mon Sep 17 00:00:00 2001 From: Michael Palumbo Date: Tue, 15 Feb 2022 11:14:56 -0500 Subject: [PATCH 4/5] zenodo url --- deps/download_soap_example.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/deps/download_soap_example.jl b/deps/download_soap_example.jl index 646176d..b49e024 100644 --- a/deps/download_soap_example.jl +++ b/deps/download_soap_example.jl @@ -3,8 +3,7 @@ include("download.jl") -# download_url = "https://pennstateoffice365-my.sharepoint.com/:u:/g/personal/ebf11_psu_edu/EYyYhOHEjM5RqvZn1D-4I58BDbpbcYg6QZKhd0KQ3N7YoQ?e=Is4PcG" -download_url = "https://pennstateoffice365-my.sharepoint.com/:u:/g/personal/ebf11_psu_edu/EYyYhOHEjM5RqvZn1D-4I58BDbpbcYg6QZKhd0KQ3N7YoQ?download=1" +download_url = "https://zenodo.org/record/6077993/files/soap_demo.h5?download=1" download_filename = joinpath("..","data","spectra","soap_demo.h5") download_md5 = "18a6291039d7684ef9187762e4b2dfd5" From 97f5289e95a7977337e1cce949295b3d2ba14c41 Mon Sep 17 00:00:00 2001 From: Michael Palumbo Date: Tue, 15 Feb 2022 11:36:50 -0500 Subject: [PATCH 5/5] revert to eric's code --- src/calc_rv/fit_quadratic_to_ccf.jl | 22 +++++++--------------- 1 file changed, 7 insertions(+), 15 deletions(-) diff --git a/src/calc_rv/fit_quadratic_to_ccf.jl b/src/calc_rv/fit_quadratic_to_ccf.jl index c661611..eea2a8b 100644 --- a/src/calc_rv/fit_quadratic_to_ccf.jl +++ b/src/calc_rv/fit_quadratic_to_ccf.jl @@ -31,23 +31,15 @@ using Statistics function (mrv::MeasureRvFromCCFQuadratic)(vels::A1, ccf::A2 ) where {T1<:Real, A1<:AbstractArray{T1,1}, T2<:Real, A2<:AbstractArray{T2,1} } # find the min and fit only the part near the minimum of the CCF amin, inds = find_idx_at_and_around_minimum(vels, ccf, frac_of_width_to_fit=mrv.frac_of_width_to_fit, measure_width_at_frac_depth=mrv.measure_width_at_frac_depth) - if isnan(amin) return (rv=NaN, σ_rv=NaN) end - - # set up linear system to solve - A = ones(3,3) - A[[3,5,7]] .= sum(view(vels, inds).^2) # set anti-diagonal indices - A[1,1] = length(view(vels, inds)) - A[1,2] = A[2,1] = sum(view(vels, inds)) - A[2,3] = A[3,2] = sum(view(vels, inds).^3) - A[3,3] = sum(view(vels, inds).^4) + if isnan(amin) return ( rv=NaN, σ_rv=NaN ) end - rhs = ones(3) - rhs[1] = sum(view(ccf,inds)) - rhs[2] = sum(view(ccf,inds) .* view(vels,inds)) - rhs[3] = sum(view(ccf,inds) .* view(vels,inds).^2) - (c, b, a) = A \ rhs + mean_v = mean(view(vels,inds)) + X = ones(length(inds),3) + X[:,2] .= view(vels,inds) .-mean_v + X[:,3] .= (view(vels,inds) .-mean_v).^2 + (c, b, a) = (X'*X) \ (X'*view(ccf,inds)) - v_at_min_of_quadratic = -b/(2*a) + v_at_min_of_quadratic = -b/(2*a) + mean_v return (rv=v_at_min_of_quadratic, σ_rv=NaN) end