Add adammode to NEOParameters

scripts/orbitdetermination.jl

@@ -87,8 +87,8 @@ end
         params[1] = NEOParameters(
             coeffstol = Inf, bwdoffset = 0.042, fwdoffset = 0.042, # Propagation
             gaussorder = 6, gaussQmax = 1.0,                       # Gauss method
-            adamiter = 500, adamQtol = 1e-5, tsaQmax = 1.0,        # ADAM
-            jtlsiter = 20, lsiter = 10,                            # Least squares
+            adamiter = 500, adammode = false, adamQtol = 1e-5,     # ADAM
+            tsaQmax = 1.0, jtlsiter = 20, lsiter = 10,             # Least squares
             outrej = true, χ2_rec = 7.0, χ2_rej = 8.0,             # Outlier rejection
             fudge = 400.0, max_per = 20.0
         )

src/observations/tracklet.jl

@@ -48,6 +48,8 @@ vdec(x::Tracklet) = x.v_δ
 mag(x::Tracklet) = x.mag
 nobs(x::Tracklet) = x.nobs
 indices(x::Tracklet) = x.indices
+indices(x::AbstractVector{Tracklet{T}}) where {T <: AbstractFloat} =
+    sort!(reduce(vcat, indices.(x)))
 
 # Print method for Tracklet{T}
 # Examples:

src/orbitdetermination/orbitdetermination.jl

@@ -101,7 +101,7 @@ function _initialtracklets(trksa::AbstractVector{Tracklet{T}},
     dts = @. abs(dtutc2et(date(tout)) - et)
     permute!(tout, sortperm(dts))
     # Starting observations
-    rin = sort!(reduce(vcat, indices.(tin)))
+    rin = indices(tin)
     # jtls needs at least three observations
     while length(rin) < 3 && !isempty(tout)
         tracklet = popfirst!(tout)

src/orbitdetermination/tooshortarc.jl

@@ -23,8 +23,8 @@ function adam(od::ODProblem{D, T}, i::Int, A::AdmissibleRegion{T}, ρ::T, v_ρ::
     μ::T = 0.75, ν::T = 0.9, ϵ::T = 1e-8, adamorder::Int = 2) where {D, T <: Real}
     # Initial time of integration [julian days TDB]
     jd0 = dtutc2jdtdb(A.date)
-    # Unfold maximum number of iterations and relative tolerance
-    maxiter, Qtol = params.adamiter, params.adamQtol
+    # Unfold parameters
+    maxiter, mode, Qtol = params.adamiter, params.adammode, params.adamQtol
     # Allocate memory
     aes = Matrix{T}(undef, 6, maxiter+1)
     Qs = fill(T(Inf), maxiter+1)
@@ -43,7 +43,7 @@ function adam(od::ODProblem{D, T}, i::Int, A::AdmissibleRegion{T}, ρ::T, v_ρ::
     dae = [scalings[i] * TaylorN(i, order = adamorder) for i in 1:6]
     AE = aes[:, 1] .+ dae
     # Subset of radec
-    idxs = indices(od.tracklets[i])
+    idxs = mode ? indices(od.tracklets) : indices(od.tracklets[i])
     # Propagation buffer
     buffer = PropagationBuffer(od, jd0, idxs[1], idxs[end], AE, params)
     # Vector of O-C residuals
@@ -132,10 +132,12 @@ function tsaiod(od::ODProblem{D, T}, params::NEOParameters{T};
     initcond::I = iodinitcond) where {D, I, T <: Real}
     # Allocate memory for orbit
     sol = zero(NEOSolution{T, T})
+    # Unfold parameters
+    varorder, mode, Qmax = params.tsaorder, params.adammode, params.tsaQmax
     # Iterate tracklets
     for i in eachindex(od.tracklets)
         # ADAM requires a minimum of 2 observations
-        od.tracklets[i].nobs < 2 && continue
+        nobs(od.tracklets[i]) < 2 && continue
         # Admissible region
         A = AdmissibleRegion(od.tracklets[i], params)
         # List of naive initial conditions
@@ -155,13 +157,14 @@ function tsaiod(od::ODProblem{D, T}, params::NEOParameters{T};
             # Scaling factors
             scalings = abs.(q0) ./ 10^5
             # Jet Transport initial condition
-            q = [q0[k] + scalings[k] * TaylorN(k, order = params.tsaorder) for k in 1:6]
+            q = [q0[k] + scalings[k] * TaylorN(k, order = varorder) for k in 1:6]
             # Jet Transport Least Squares
-            _sol_ = jtls(od, jd0, q, od.tracklets[i:i], params, false)
+            trks = mode ? od.tracklets[:] : od.tracklets[i:i]
+            _sol_ = jtls(od, jd0, q, trks, params, mode)
             # Update solution
             sol = updatesol(sol, _sol_, od.radec)
             # Termination condition
-            nrms(sol) <= params.tsaQmax && return sol
+            nrms(sol) <= Qmax && return sol
         end
     end
 

src/propagation/parameters.jl

@@ -31,7 +31,9 @@ Parameters for all orbit determination functions.
 - `H_max::T`: maximum absolute magnitude (default: `34.5`).
 - `a_max::T`: maximum semimajor axis (default: `100.0`).
 - `adamiter::Int`: maximum number of iterations for `ADAM` optimizer (default: `200`).
-- `Qtol::T`: target function relative tolerance (default: `0.001`).
+- `adammode::Bool`: whether to perform ADAM iterations with all the observations
+    (default: `false`).
+- `adamQtol::T`: target function relative tolerance (default: `0.001`).
 - `tsaorder::Int`: order of the jet transport perturbation (default: `6`).
 - `tsaQmax::T`: nrms threshold (default: `1.5`).
 
@@ -72,6 +74,7 @@ Parameters for all orbit determination functions.
     H_max::T = 34.5
     a_max::T = 100.0
     adamiter::Int = 200
+    adammode::Bool = false
     adamQtol::T = 0.001
     tsaorder::Int = 6
     tsaQmax::T = 1.5

test/orbitdetermination.jl

@@ -11,7 +11,7 @@ using NEOs: NEOSolution, RadecMPC, reduce_tracklets,
 
 function iodsubradec(radec::Vector{RadecMPC{T}}, N::Int = 3) where {T <: Real}
     tracklets = reduce_tracklets(radec)
-    idxs = sort!(reduce(vcat, indices.(tracklets[1:N])))
+    idxs = indices(tracklets[1:N])
     subradec = radec[idxs]
     return subradec
 end