Relaxation factor

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,20 +122,24 @@ 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
     eph_su = selecteph(sseph, su)
     # 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,30 +149,24 @@ 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
 
     # Iterate over triplets
     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

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

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")

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
+function gaussinitcond end
+function relax_factor end

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