rv error + GPU updates

scripts/NEID_RV_ccf.jl

@@ -359,28 +359,35 @@ function line_rvs_ccf(line_names, vacwav, orders, timestamps, path)
     # iterate through lines and determine line RV for eclipse EM curve
     Threads.@threads for i in 1:length(line_names)
         order_index = orders[i]
-        min_wav = vacwav[i] - 2
-        max_wav = vacwav[i] + 2
+        if i == 11
+            min_wav = vacwav[i] - 0.2
+            max_wav = vacwav[i] + 0.2
+        else
+            min_wav = vacwav[i] - 0.35
+            max_wav = vacwav[i] + 0.35
+        end
 
         lines = [vacwav[i]]
 
         RV_list = Vector{Float64}(undef,length(timestamps)...)
         RV_error_list = Vector{Float64}(undef,length(timestamps)...)
         Threads.@threads for j in 1:length(timestamps)
             f = FITS(joinpath(path, timestamps[j]))
-            # header = read_header(f[1])
-            # bc_ms = (header["SSBRV0$order_index"]) * 1000
+            header = read_header(f[1])
+            bc_ms = (header["SSBRV0$order_index"]) * 1000
 
             spectrum = NEID.read_data(joinpath(path, timestamps[j]); normalization = :blaze)
             # find where line is - NEID flux
             chunk = NEID.ChunkOfSpectrum(spectrum, order_index, full_pixels) 
             pixels = NEID.find_pixels_for_line_in_chunk(chunk, min_wav, max_wav)
             chunk_flux_full = chunk.flux[pixels]
             chunk_flux_full ./= maximum(chunk_flux_full)
-            chunck_vac_wav = chunk.λ[pixels]
+            chunck_vac_wav = chunk.λ[pixels] 
             chunck_var = chunk.var[pixels] ./ maximum(chunk_flux_full)^2
 
-            v_grid_cpu, ccf_cpu, ccf_var_out = GRASS.calc_ccf(chunck_vac_wav, chunk_flux_full, chunck_var, lines, [maximum(chunk_flux_full) - minimum(chunk_flux_full)], resolution)
+            v_grid_cpu, ccf_cpu, ccf_var_out = GRASS.calc_ccf(chunck_vac_wav, chunk_flux_full, chunck_var, lines, 
+                                                [maximum(chunk_flux_full) - minimum(chunk_flux_full)], resolution,
+                                                Δv_max=8000.0)         
             rvs_cpu, sigs_cpu = GRASS.calc_rvs_from_ccf(v_grid_cpu, ccf_cpu, ccf_var_out)  
 
             RV_list[j] = rvs_cpu
@@ -456,21 +463,6 @@ function neid_nxt_day()
     end
 end
 
-function neid_nxt_day()
-    path_october = "/storage/group/ebf11/default/pipeline/neid_solar/data/v1.3/L2/2023/10/15/"
-    files_and_dirs = readdir(path_october)
-    files = filter(f -> isfile(joinpath(path_october, f)), files_and_dirs)
-    timestamps_october = files[4:length(files)-145]
-
-    RV_all_lines, RV_error_all_lines = line_rvs_ccf(line_names, vacwav, orders, timestamps_october, path_october)
-    @save "neid_RVlinebyline_nxt_day.jld2"
-    jldopen("neid_RVlinebyline_nxt_day.jld2", "a+") do file
-        file["name"] = line_names 
-        file["rv"] = RV_all_lines 
-        file["rv_error"] = RV_error_all_lines 
-    end
-end
-
 function compare_spec_sources() # out of transit line comparsion between IAS, GRASS, and NEID
     neid_timestamps = ["2023-10-14T19:03:06.500000"] # last timestamp (out of transit)
     path_october = "/storage/group/ebf11/default/pipeline/neid_solar/data/v1.3/L2/2023/10/14/"
@@ -484,6 +476,6 @@ function compare_spec_sources() # out of transit line comparsion between IAS, GR
                             timestamps, path_october, "neid_october_N_50", "SSD")
 end
 
-# neid_eclipse()
+neid_eclipse()
 # neid_nxt_day()
-compare_spec_sources()
+# compare_spec_sources()

scripts/debug_gpu_eclipse.jl

@@ -0,0 +1,79 @@
+using Revise
+using GRASS
+using SPICE
+using CSV
+using DataFrames
+import PyPlot
+plt = PyPlot
+
+#NEID location
+obs_lat = 31.9583 
+obs_long = -111.5967  
+alt = 2.097938
+
+#convert from utc to et as needed by SPICE
+time_stamps = utc2et.(["2023-10-14T14:29:23.500000"])
+
+N = 50
+Nt = length(time_stamps)
+Nsubgrid = 4
+disk = GRASS.DiskParamsEclipse(N=N, Nt=Nt, Nsubgrid=Nsubgrid)
+
+lines = [5434.5232] # array of line centers 
+depths = [0.6]   # array of line depths
+templates = ["FeI_5434"] # template data to use
+variability = trues(length(lines))  # whether or not the bisectors should "dance"
+blueshifts = zeros(length(lines))   # set convective blueshift value
+resolution = 7e5                    # spectral resolution
+
+# make the disk and spec composite type instances
+spec = GRASS.SpecParams(lines=lines, depths=depths, variability=variability,
+                        blueshifts=blueshifts, templates=templates, resolution=resolution) 
+
+"""
+CPU GEOMETRY CODE CALLS BELOW
+"""
+
+wsp = GRASS.SynthWorkspaceEclipse(disk, 1, Nt, verbose=true)
+mem = GRASS.GeoWorkspaceEclipse(disk, 1, Nt)
+LD_type = "SSD"
+
+mu_grid_matrix = Matrix{Matrix{Float64}}(undef, length(disk.ϕc), maximum(disk.Nθ))    
+dA_total_proj_matrix = Matrix{Matrix{Float64}}(undef, length(disk.ϕc), maximum(disk.Nθ))
+idx1_matrix = Matrix{Matrix{Int}}(undef, length(disk.ϕc), maximum(disk.Nθ))
+idx3_matrix = Matrix{Matrix{Int}}(undef, length(disk.ϕc), maximum(disk.Nθ))
+z_rot_matrix = Matrix{Matrix{Float64}}(undef, length(disk.ϕc), maximum(disk.Nθ))
+zenith_matrix = Matrix{Matrix{Float64}}(undef, length(disk.ϕc), maximum(disk.Nθ))
+
+# compute geometry for timestamp
+t = 1
+
+GRASS.eclipse_compute_quantities!(disk, time_stamps[t], t, obs_long, obs_lat, alt, wsp.ϕc, wsp.θc, 
+                                  wsp.μs, wsp.dA, wsp.xyz, wsp.ax_codes, 
+                                  mem.dA_total_proj_mean, mem.mean_weight_v_no_cb,
+                                  mem.mean_weight_v_earth_orb, mem.pole_vector_grid,
+                                  mem.SP_sun_pos, mem.SP_sun_vel, mem.SP_bary, mem.SP_bary_pos,
+                                  mem.SP_bary_vel, mem.OP_bary, mem.mu_grid, mem.projected_velocities_no_cb, 
+                                  mem.distance, mem.v_scalar_grid, mem.v_earth_orb_proj, mu_grid_matrix,
+                                  dA_total_proj_matrix, idx1_matrix, idx3_matrix, z_rot_matrix, zenith_matrix)    
+
+lambdas_cpu, outspec_cpu = GRASS.eclipse_compute_intensity(disk, lines, 1.0, LD_type, idx1_matrix, idx3_matrix,
+                                mu_grid_matrix, z_rot_matrix, dA_total_proj_matrix, wsp.ld, wsp.z_rot, zenith_matrix, 
+                                wsp.μs, wsp.ax_codes, wsp.dA, wsp.μs, wsp.ax_codes, wsp.dA, true, t, wsp.ext)
+
+"""
+GPU GEOMETRY CODE CALLS BELOW
+"""
+gpu_allocs = GRASS.GPUAllocsEclipse(spec, disk, 1)
+
+GRASS.calc_eclipse_quantities_gpu!(time_stamps[t], obs_long, obs_lat, alt, lines, LD_type, 1.0, 1.0, disk, gpu_allocs)
+
+# make a plot
+plt.imshow(wsp.dA .- Array(gpu_allocs.dA))
+# plt.imshow(wsp.z_rot .* GRASS.c_ms, vmin=-2000, vmax=2000, cmap="seismic")
+# plt.imshow(Array(gpu_allocs.z_rot) .* GRASS.c_ms .- wsp.z_rot .* GRASS.c_ms)
+plt.title("dA")
+plt.colorbar()
+plt.savefig("gpu_test.png")
+plt.clf()
+#z_rot & dA

src/gpu/gpu_precomps_eclipse.jl

@@ -249,15 +249,15 @@ function calc_eclipse_quantities_gpu!(wavelength, μs, z_rot, ax_codes,
                 dA_sub *= μ_sub
                 dA_sum += dA_sub
 
-                # iterate counter 
-                count += 1
-
                 # calculate distance
                 n2 = CUDA.sqrt(OM_bary[1]^2.0 + OM_bary[2]^2.0 + OM_bary[3]^2.0)  
                 d2 = acos((OM_bary[1] * OP_bary_x + OM_bary[2] * OP_bary_y + OM_bary[3] * OP_bary_z) / (n2 * n1))
                 if (d2 < atan(moon_radius/n2))
                     continue
                 end
+
+                # iterate counter 
+                count += 1
 
                 # zenith
                 n3 = CUDA.sqrt(EO_bary[1]^2.0 + EO_bary[2]^2.0 + EO_bary[3]^2.0) 

src/velocities.jl

@@ -6,6 +6,7 @@ import EchelleCCFs: MeasureRvFromCCFGaussian as GaussianFit
 import EchelleCCFs: MeasureRvFromCCFQuadratic as QuadraticFit
 import EchelleCCFs: AbstractCCFMaskShape as MaskShape
 import EchelleCCFs: AbstractMeasureRvFromCCF as FitType
+import RvSpectMLBase: searchsortednearest
 
 """
     calc_ccf(lambdas, flux, lines, depths, resolution; normalize=true)
@@ -238,41 +239,104 @@ Calculate apparent radial velocity from a CCF and velocity grid.
 - `v_grid::AbstractArray{Float64,1}`: List of velocities returned by calc_ccf.
 - `ccf::AbstractArray{Float64,1}`: CCF values returned by calc_ccf.
 """
-function calc_rvs_from_ccf(v_grid::AA{Float64,1}, ccf::AA{Float64,1}, var::AA{Float64,1};
-                           frac_of_width_to_fit::Float64=0.75,
-                           fit_type::Type{T}=GaussianFit) where {T<:FitType}
+function calc_rvs_from_ccf(v_grid::AA{Float64,1}, ccf::AA{Float64,1};
+                            frac_of_width_to_fit::Float64=0.75,
+                            fit_type::Type{T}=GaussianFit) where {T<:FitType}
     mrv = T(frac_of_width_to_fit=frac_of_width_to_fit)
-    return mrv(v_grid, ccf, var)
+    return mrv(v_grid, ccf)
 end
 
-function calc_rvs_from_ccf(v_grid::AA{T,1}, ccf::AA{T,2}, var::AA{T,2}; kwargs...) where {T<:Float64}
-    func = x -> calc_rvs_from_ccf(v_grid, x, var; kwargs...)
+function calc_rvs_from_ccf(v_grid::AA{T,1}, ccf::AA{T,2}; kwargs...) where {T<:Float64}
+    func = x -> calc_rvs_from_ccf(v_grid, x; kwargs...)
     out = mapslices(func, ccf, dims=1)
     return vec.(unzip(out))
-end  
+end   
 
-function calc_rvs_from_ccf(v_grid::AA{Float64,1}, ccf::AA{Float64,1};
+function calc_rvs_from_ccf(v_grid::AbstractArray{Float64,1}, ccf::AbstractArray{Float64,1}, var::AbstractArray{Float64,1};
                            frac_of_width_to_fit::Float64=0.75,
                            fit_type::Type{T}=GaussianFit) where {T<:FitType}
     mrv = T(frac_of_width_to_fit=frac_of_width_to_fit)
-    return mrv(v_grid, ccf)
+    return MeasureRvFromCCFGaussian_New(v_grid, ccf, var, mrv)
 end
 
-function calc_rvs_from_ccf(v_grid::AA{T,1}, ccf::AA{T,2}; kwargs...) where {T<:Float64}
-    func = x -> calc_rvs_from_ccf(v_grid, x; kwargs...)
-    out = mapslices(func, ccf, dims=1)
-    return vec.(unzip(out))
-end  
+@. gaussian_line_helper(x, p) = p[4] + p[3] * exp(-0.5*((x-p[1])/p[2])^2)
+
+function MeasureRvFromCCFGaussian_New(vels::A1, ccf::A2, ccf_var::A3, mrv) where {T1<:Real, A1<:AbstractArray{T1,1}, T2<:Real, A2<:AbstractArray{T2,1}, T3<:Real, A3<:AbstractArray{T3,1} }
+    if all(ccf .== zero(eltype(ccf))) return (rv=NaN, σ_rv=NaN)     end
+    # 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
+
+    # make initial guess parameters
+    μ = vels[amin]
+    σ = mrv.init_guess_ccf_σ                   
+    minccf, maxccf = extrema(ccf)
+    amp = minccf - maxccf
+    y0 = maxccf
+    p0 = [μ, σ, amp, y0]
+
+    local rvfit
+   # fit and return the mean of the distribution
+    result = curve_fit(gaussian_line_helper, view(vels,inds), view(ccf,inds), (1.0 ./ view(ccf_var,inds)),  p0)
+
+    if result.converged
+          rv = coef(result)[1]
+          sigma_rv = stderror_new(result)[1]
+          rvfit = (rv=rv, σ_rv=sigma_rv)
+    end
+    return rvfit
+end
 
-function calc_rvs_from_ccf(v_grid::AA{Float64,1}, ccf::AA{Float64,1};
-                            frac_of_width_to_fit::Float64=0.75,
-                            fit_type::Type{T}=GaussianFit) where {T<:FitType}
-    mrv = T(frac_of_width_to_fit=frac_of_width_to_fit)
-    return mrv(v_grid, ccf)
+function vcov_new(fit)
+    # computes covariance matrix of fit parameters
+    J = fit.jacobian
+
+    covar = inv(J' * J) * mse(fit)
+    return covar
 end
 
-function calc_rvs_from_ccf(v_grid::AA{T,1}, ccf::AA{T,2}; kwargs...) where {T<:Float64}
-    func = x -> calc_rvs_from_ccf(v_grid, x; kwargs...)
-    out = mapslices(func, ccf, dims=1)
-    return vec.(unzip(out))
-end   
+function stderror_new(fit; rtol::Real=NaN, atol::Real=0)
+    covar = vcov_new(fit)
+    # then the standard errors are given by the sqrt of the diagonal
+    vars = diag(covar)
+    return sqrt.(abs.(vars))
+end
+
+function est_full_width(vels::A1, ccf::A2; measure_width_at_frac_depth::Real = default_measure_width_at_frac_depth ) where  {T1<:Real, A1<:AbstractArray{T1,1}, T2<:Real, A2<:AbstractArray{T2,1} }
+    minccf, maxccf = extrema(ccf)
+    depth = maxccf - minccf
+    target_val = minccf + (1-measure_width_at_frac_depth) * depth
+    ind1 = findfirst(ccf .<= target_val)
+    ind2 = findlast(ccf .<= target_val)
+    if isnothing(ind1) || isnothing(ind2)
+        return NaN
+        println("ccf = ",ccf)
+        println("minccf= ", minccf, " maxccf= ", maxccf, " depth= ", depth, " measure_width_at_frac_depth= ", measure_width_at_frac_depth, " targetval= ",target_val, " ind1= ", ind1, " ind2= ", ind2)
+        @error "est_full_width failed."
+    end
+    return vels[ind2] - vels[ind1]
+end
+
+function find_idx_at_and_around_minimum(vels::A1, ccf::A2; frac_of_width_to_fit::Real = default_frac_of_width_to_fit, measure_width_at_frac_depth::Real = default_measure_width_at_frac_depth) where {T1<:Real, A1<:AbstractArray{T1,1}, T2<:Real, A2<:AbstractArray{T2,1} }
+    # do a prelim fit to get the width
+    full_width = est_full_width(vels, ccf, measure_width_at_frac_depth=measure_width_at_frac_depth)
+    if isnan(full_width)
+       return (NaN, 1:length(vels))
+    end
+
+    # find the min and fit only that
+    amin = argmin(ccf)
+    if amin == 1 || amin==length(vels)
+        offset = max(1,floor(Int64,length(vels)//4))
+        amin = argmin(view(ccf,offset:(length(vels)-offset)))
+        amin += offset-1
+    end
+    lend = vels[amin] - frac_of_width_to_fit * full_width
+    rend = vels[amin] + frac_of_width_to_fit * full_width
+    # get the indices
+    lind = searchsortednearest(view(vels,1:amin), lend)
+    rind = amin + searchsortednearest(view(vels,(amin+1):length(vels)), rend)
+    inds = lind:rind
+
+    return (amin, inds)
+end