From f08c44577efbf214fd526cbd544b8cc7a7a7ec14 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Luis=20Eduardo=20Ram=C3=ADrez=20Montoya?= Date: Mon, 7 Aug 2023 15:28:05 -0600 Subject: [PATCH] Relaxation factor --- ext/DataFramesExt.jl | 86 +++++++++++++------------ pha/apophis.jl | 24 +++---- src/NEOs.jl | 2 +- src/orbit_determination/gauss_method.jl | 3 +- src/postprocessing/least_squares.jl | 60 ++++++++++++----- 5 files changed, 103 insertions(+), 72 deletions(-) diff --git a/ext/DataFramesExt.jl b/ext/DataFramesExt.jl index 8bc7fc49..d1f07cd2 100644 --- a/ext/DataFramesExt.jl +++ b/ext/DataFramesExt.jl @@ -4,10 +4,10 @@ using Dates: Date, Period, Hour, Day, datetime2julian, julian2datetime using TaylorSeries: get_numvars using PlanetaryEphemeris: J2000, selecteph, su, ea, yr, daysec, auday2kmsec using NEOs: RadecMPC, date, gauss_triplets, propagate, RNp1BP_pN_A_J23E_J2S_eph_threads!, order, abstol, sseph, - scaled_variables, gauss_method, residuals, bwdfwdeph, newtonls, diffcorr, nrms, hascoord + scaled_variables, gauss_method, residuals, bwdfwdeph, newtonls, diffcorr, nrms, hascoord, tryls import Base: convert -import NEOs: AbstractAstrometry, extrapolation, reduce_nights, gaussinitcond +import NEOs: AbstractAstrometry, extrapolation, reduce_nights, gaussinitcond, relax_factor if isdefined(Base, :get_extension) using DataFrames: AbstractDataFrame, DataFrame, nrow, eachcol, eachrow, groupby, combine @@ -72,8 +72,6 @@ via polynomial interpolation. function reduce_nights(radec::Vector{RadecMPC{T}}) where {T <: AbstractFloat} # Convert to DataFrame df = DataFrame(radec) - # Eliminate observatories without coordinates - filter!(:observatory => hascoord, df) # Group by observatory and Date df.Date = Date.(df.date) gdf = groupby(df, [:observatory, :Date]) @@ -88,6 +86,20 @@ function reduce_nights(radec::Vector{RadecMPC{T}}) where {T <: AbstractFloat} return cdf.observatory, cdf.date, cdf.α, cdf.δ end +function relax_factor(radec::Vector{RadecMPC{T}}) where {T <: AbstractFloat} + # Convert to DataFrame + df = DataFrame(radec) + # Group by observatory and Date + df.Date = Date.(df.date) + gdf = groupby(df, [:observatory, :Date]) + # Interpolate observation nights + cdf = combine(gdf, nrow) + # Count observations in each group + Nv = cdf[gdf.groups, :nrow] + # Relaxation factor + return map(x -> x > 4.0 ? x/4.0 : 1.0, Nv) +end + @doc raw""" gaussinitcond(radec::Vector{RadecMPC{T}}; Δ::Period = Day(1), Q_max::T = 0.75, niter::Int = 5, maxsteps::Int = 100, varorder::Int = 5, order::Int = order, abstol::T = abstol, parse_eqs::Bool = true) where {T <: AbstractFloat} @@ -110,7 +122,7 @@ See also [`gauss_method`](@ref). !!! warning This function will set the (global) `TaylorSeries` variables to `δα₁ δα₂ δα₃ δδ₁ δδ₂ δδ₃`. """ -function gaussinitcond(radec::Vector{RadecMPC{T}}; max_triplets::Int = 10, Q_max::T = 0.5, niter::Int = 5, maxsteps::Int = 100, +function gaussinitcond(radec::Vector{RadecMPC{T}}; max_triplets::Int = 10, Q_max::T = 100., niter::Int = 5, maxsteps::Int = 100, varorder::Int = 5, order::Int = order, abstol::T = abstol, parse_eqs::Bool = true) where {T <: AbstractFloat} # Sun's ephemeris @@ -118,12 +130,16 @@ function gaussinitcond(radec::Vector{RadecMPC{T}}; max_triplets::Int = 10, Q_max # Earth's ephemeris eph_ea = selecteph(sseph, ea) + # Eliminate observatories without coordinates + filter!(x -> hascoord(x.observatory), radec) # Reduce nights by interpolation observatories, dates, α, δ = reduce_nights(radec) # Observations triplets triplets = gauss_triplets(dates, max_triplets) - - # Initial date of integration [julian days] + + # Julian day of first (last) observation + t0, tf = datetime2julian(radec[1].date), datetime2julian(radec[end].date) + # Julian day when to start propagation jd0 = zero(T) # Jet transport perturbation (ra/dec) @@ -133,9 +149,6 @@ function gaussinitcond(radec::Vector{RadecMPC{T}}; max_triplets::Int = 10, Q_max best_Q = T(Inf) # Initial conditions best_Q0 = zeros(T, 6) - - # Global counter - k = 1 # Break flag flag = false @@ -143,20 +156,17 @@ function gaussinitcond(radec::Vector{RadecMPC{T}}; max_triplets::Int = 10, Q_max for j in eachindex(triplets) # Current triplet - idxs = triplets[j] + triplet = triplets[j] - # [julian days] - t0, jd0, tf = datetime2julian.(dates[idxs]) + # Julian day when to start propagation + jd0 = datetime2julian(dates[triplet[2]]) # Number of years in forward integration nyears_fwd = (tf - jd0 + 2) / yr # Number of years in backward integration nyears_bwd = -(jd0 - t0 + 2) / yr - # Subset of radec for residuals - sub_radec = filter(x -> dates[idxs[1]] <= date(x) <= dates[idxs[3]], radec) - # Gauss method solution - sol = gauss_method(observatories[idxs], dates[idxs], α[idxs] .+ dq[1:3], δ[idxs] .+ dq[4:6]; niter = niter) + sol = gauss_method(observatories[triplet], dates[triplet], α[triplet] .+ dq[1:3], δ[triplet] .+ dq[4:6]; niter = niter) # Iterate over Gauss solutions for i in eachindex(sol) @@ -169,45 +179,37 @@ function gaussinitcond(radec::Vector{RadecMPC{T}}; max_triplets::Int = 10, Q_max order = order, abstol = abstol, parse_eqs = parse_eqs) # O-C residuals - res, w = residuals(sub_radec; xvs = et -> auday2kmsec(eph_su(et/daysec)), xve = et -> auday2kmsec(eph_ea(et/daysec)), - xva = et -> bwdfwdeph(et, bwd, fwd)) + res, w = residuals(radec; xvs = et -> auday2kmsec(eph_su(et/daysec)), xve = et -> auday2kmsec(eph_ea(et/daysec)), + xva = et -> bwdfwdeph(et, bwd, fwd)) + + # Subset of radec for orbit fit + i_O = findall(x -> dates[triplet[1]] <= date(x) <= dates[triplet[3]], radec) + i_ξ = vcat(i_O, i_O .+ length(radec)) + + # Relaxation factor + rex = relax_factor(radec) + rex = vcat(rex, rex) + w = w ./ rex + + # Orbit fit + fit = tryls(res[i_ξ], w[i_ξ], zeros(get_numvars()), niter) - # Orbit fit (Newton) - success, x_new, Γ = newtonls(res, w, zeros(get_numvars()), niter) # NRMS - Q = nrms(res(x_new), w) + Q = nrms(res(fit.x), w) # TO DO: check cases where newton converges but diffcorr no - if success + if fit.success # Update NRMS and initial conditions if Q < best_Q best_Q = Q - best_Q0 .= bwd(bwd.t0)(x_new) + best_Q0 .= bwd(bwd.t0)(fit.x) end # Break condition if Q <= Q_max flag = true end - else - # Orbit fit (differential corrections) - success, x_new, Γ = diffcorr(res, w, zeros(get_numvars()), niter) - # NRMS - Q = nrms(res(x_new), w) - - if success - # Update NRMS and initial conditions - if Q < best_Q - best_Q = Q - best_Q0 .= bwd(bwd.t0)(x_new) - end - # Break condition - if Q <= Q_max - flag = true - end - end end - k += 1 # Break condition if flag break diff --git a/pha/apophis.jl b/pha/apophis.jl index fdca5417..e7442acb 100644 --- a/pha/apophis.jl +++ b/pha/apophis.jl @@ -231,9 +231,9 @@ function main(dynamics::D, maxsteps::Int, jd0_datetime::DateTime, nyears_bwd::T, res = vcat(res_radec, res_del, res_dop) w = vcat(w_radec, w_del, w_dop) - success, δx_OR8, Γ_OR8 = newtonls(res, w, zeros(get_numvars()), 10) - x_OR8 = sol_fwd(sol_fwd.t0)(δx_OR8) - σ_OR8 = sqrt.(diag(Γ_OR8)).*scalings + fit_OR8 = newtonls(res, w, zeros(get_numvars()), 10) + x_OR8 = sol_fwd(sol_fwd.t0)(fit_OR8.x) + σ_OR8 = sqrt.(diag(fit_OR8.Γ)).*scalings nradec = length(res_radec) res_ra = view(res_radec, 1:nradec÷2) @@ -246,19 +246,19 @@ function main(dynamics::D, maxsteps::Int, jd0_datetime::DateTime, nyears_bwd::T, print_header("Orbital fit (8-DOF) and post-fit statistics", 2) # orbital fit - println("Success flag : ", success, "\n") + println("Success flag : ", fit_OR8.success, "\n") println("Nominal solution [au,au,au,au/d,au/d,au/d,au/d²,au/d²]: ", x_OR8, "\n") println("1-sigma formal uncertainties [au,au,au,au/d,au/d,au/d,au/d²,au/d²]: ", σ_OR8, "\n") # post-fit statistics - println("Normalized RMS (optical-only) [adimensional] : ", nrms(res_radec(δx_OR8),w_radec)) - println("Normalized RMS (radar-only) [adimensional] : ", nrms(vcat(res_del,res_dop)(δx_OR8),vcat(w_del,w_dop))) - println("Normalized RMS (combined optical and radar) [adimensional] : ", nrms(res(δx_OR8),w), "\n") - println("Mean weighted right-ascension residual [arcseconds] : ", mean(res_ra(δx_OR8), weights(w_ra))) - println("Mean weighted declination residual [arcseconds] : ", mean(res_dec(δx_OR8), weights(w_dec))) - println("Mean weighted time-delay residual [micro-seconds]: ", mean(res_del(δx_OR8), weights(w_del))) - println("Mean weighted Doppler-shift residual [Hz] : ", mean(res_dop(δx_OR8), weights(w_dop)), "\n") - println("Chi-squared statistic (χ²): [adimensional] : ", chi2(res(δx_OR8),w)) + println("Normalized RMS (optical-only) [adimensional] : ", nrms(res_radec(fit_OR8.x),w_radec)) + println("Normalized RMS (radar-only) [adimensional] : ", nrms(vcat(res_del,res_dop)(fit_OR8.x),vcat(w_del,w_dop))) + println("Normalized RMS (combined optical and radar) [adimensional] : ", nrms(res(fit_OR8.x),w), "\n") + println("Mean weighted right-ascension residual [arcseconds] : ", mean(res_ra(fit_OR8.x), weights(w_ra))) + println("Mean weighted declination residual [arcseconds] : ", mean(res_dec(fit_OR8.x), weights(w_dec))) + println("Mean weighted time-delay residual [micro-seconds]: ", mean(res_del(fit_OR8.x), weights(w_del))) + println("Mean weighted Doppler-shift residual [Hz] : ", mean(res_dop(fit_OR8.x), weights(w_dop)), "\n") + println("Chi-squared statistic (χ²): [adimensional] : ", chi2(res(fit_OR8.x),w)) return sol_bwd, sol_fwd, res_radec, res_del, res_dop, w_radec, w_del, w_dop end diff --git a/src/NEOs.jl b/src/NEOs.jl index 720a0e5a..fe0d7fa2 100644 --- a/src/NEOs.jl +++ b/src/NEOs.jl @@ -64,7 +64,7 @@ export RNp1BP_pN_A_J23E_J2S_ng_eph_threads!, RNp1BP_pN_A_J23E_J2S_eph_threads! # Propagate export propagate, propagate_lyap, propagate_root -export valsecchi_circle, nrms, chi2, newtonls, newtonls_6v, diffcorr, newtonls_Q, bopik +export valsecchi_circle, nrms, chi2, newtonls, newtonls_6v, diffcorr, newtonls_Q, bopik, tryls, project include("constants.jl") include("observations/process_radar.jl") diff --git a/src/orbit_determination/gauss_method.jl b/src/orbit_determination/gauss_method.jl index c853d631..0be7286b 100644 --- a/src/orbit_determination/gauss_method.jl +++ b/src/orbit_determination/gauss_method.jl @@ -399,4 +399,5 @@ end # Empty methods to be overloaded by DataFramesExt function reduce_nights end -function gaussinitcond end \ No newline at end of file +function gaussinitcond end +function relax_factor end \ No newline at end of file diff --git a/src/postprocessing/least_squares.jl b/src/postprocessing/least_squares.jl index 3bb2025d..5e70f170 100644 --- a/src/postprocessing/least_squares.jl +++ b/src/postprocessing/least_squares.jl @@ -1,5 +1,33 @@ include("b_plane.jl") +struct OrbitFit{T <: Real} + success::Bool + x::Vector{T} + Γ::Matrix{T} + routine::Symbol + OrbitFit{T}(success::Bool, x::Vector{T}, Γ::Matrix{T}, routine::Symbol) where {T <: Real} = new{T}(success, x, Γ, routine) +end + +OrbitFit(success::Bool, x::Vector{T}, Γ::Matrix{T}, routine::Symbol) where {T <: Real} = OrbitFit{T}(success, x, Γ, routine) + +# Print method for OrbitFit +# Examples: +# N00hp15 α: 608995.65 δ: -25653.3 t: 2020-12-04T10:41:43.209 obs: 703 +# 99942 α: 422475.3 δ: 97289.49 t: 2021-05-12T06:28:35.904 obs: F51 +function show(io::IO, fit::OrbitFit{T}) where {T <: Real} + success_s = fit.success ? "Succesful" : "Unsuccesful" + routine_s = fit.routine == :newton ? "Newton" : "differential corrections" + print(io, success_s, " ", routine_s) +end + +function project(y::Vector{TaylorN{T}}, fit::OrbitFit{T}) where {T <: Real} + J = Matrix{T}(undef, get_numvars(), length(y)) + for i in eachindex(y) + J[:, i] = TS.gradient(y[i])(fit.x) + end + return (J') * fit.Γ * J +end + @doc raw""" nrms(res, w) @@ -194,7 +222,7 @@ function diffcorr(res::Vector{TaylorN{T}}, w::Vector{T}, x0::Vector{T}, niters:: error[i+1] = sqrt(error2) # The method do not converge else - return false, x[:, i+1], inv(C) + return OrbitFit(false, x[:, i+1], inv(C), :diffcorr) end end # Index with the lowest error @@ -206,7 +234,7 @@ function diffcorr(res::Vector{TaylorN{T}}, w::Vector{T}, x0::Vector{T}, niters:: # Covariance matrix Γ = inv(C) - return true, x_new, Γ + return OrbitFit(true, x_new, Γ, :diffcorr) end @doc raw""" @@ -274,7 +302,7 @@ function newtonls(res::Vector{TaylorN{T}}, w::Vector{T}, x0::Vector{T}, niters:: error[i+1] = sqrt(error2) # The method do not converge else - return false, x[:, i+1], inv(C) + return OrbitFit(false, x[:, i+1], inv(C), :newton) end end # TO DO: study Gauss method solution dependence on jt order @@ -290,28 +318,28 @@ function newtonls(res::Vector{TaylorN{T}}, w::Vector{T}, x0::Vector{T}, niters:: # Covariance matrix Γ = inv(C) - return true, x_new, Γ + return OrbitFit(true, x_new, Γ, :newton) end function tryls(res::Vector{TaylorN{T}}, w::Vector{T}, x0::Vector{T}, niters::Int = 5) where {T <: Real} # Newton's method - success_1, x_1, Γ_1 = newtonls(res, w, x0, niters) + fit_1 = newtonls(res, w, x0, niters) # Differential corrections - success_2, x_2, Γ_2 = diffcorr(res, w, x0, niters) - if success_1 && success_2 - Q_1 = nrms(res(x_1), w) - Q_2 = nrms(res(x_2), w) + fit_2 = diffcorr(res, w, x0, niters) + if fit_1.success && fit_2.success + Q_1 = nrms(res(fit_1.x), w) + Q_2 = nrms(res(fit_2.x), w) if Q_1 <= Q_2 - return success_1, x_1, Γ_1 + return fit_1 else - return success_2, x_2, Γ_2 + return fit_2 end - elseif success_1 - return success_1, x_1, Γ_1 - elseif success_2 - return success_2, x_2, Γ_2 + elseif fit_1.success + return fit_1 + elseif fit_2.success + return fit_2 else - return false, x_1, Γ_1 + return fit_1 end end