Introduce ODProblem (#87)

* Add errormodel.jl

* Add odproblem.jl

* Update process_radec.jl

* Update parameters.jl

* Update orbit determination functions

* Update orbit determination tests

* Update propagation tests

* Update bplane tests

* Update scripts/orbitdetermination.jl

* Minor additions

* New PropagationBuffer constructor

* Update scripts/mcmov.jl

* Update examples/recipes.jl

* Bump minor version

* Update scripts/adam.jl

* Minor fix

Project.toml

@@ -1,7 +1,7 @@
 name = "NEOs"
 uuid = "b41c07a2-2abb-11e9-070a-c3c1b239e7df"
 authors = ["Jorge A. Pérez Hernández", "Luis Benet", "Luis Eduardo Ramírez Montoya"]
-version = "0.9.1"
+version = "0.10.0"
 
 [deps]
 Artifacts = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

examples/recipes.jl

@@ -2,23 +2,24 @@
 # ext/NEOsRecipesBaseExt.jl
 using NEOs, Plots
 
-# Download optical astrometry of asteroid 2023 DW
-radec = fetch_radec_mpc("2023 DW")
+# Download optical astrometry of asteroid 2016 TU93
+radec = fetch_radec_mpc("2016 TU93")
 
 # Orbit determination
 params = NEOParameters(bwdoffset = 0.007, fwdoffset = 0.007)
-sol = orbitdetermination(radec, params)
+od = ODProblem(newtonian!, radec)
+sol = orbitdetermination(od, params)
 
 # Plot optical residuals
 plot(sol.res, label = "", xlabel = "αcos(δ) [arcsec]", ylabel = "δ [arcsec]",
-     title = "2023 DW optical O-C residuals")
+     title = "2016 TU93 optical O-C residuals")
 
 # Plot asteroid's orbit (in barycentric cartesian coordinates)
 # N: number of points between t0 and tf (100 by default)
 # projection: which coordinates to plot, options are: :x, :y, :z, :xy,
 #             :xz, :yz and :xyz (default)
 t0, tf = sol.bwd.t0 + sol.bwd.t[end], sol.fwd.t0 + sol.fwd.t[end]
-plot(sol, t0, tf, N = 10_000, projection = :xyz, label = "2023 DW",
+plot(sol, t0, tf, N = 10_000, projection = :xyz, label = "2016 TU93",
      xlabel = "x [au]", ylabel = "y [au]", zlabel = "z [au]", color = :blue)
 
 # We can also visualize the orbits of the Sun and the Earth
@@ -33,6 +34,5 @@ plot!(params.eph_ea, t0, tf, N = 10_000, projection = :xyz,
 using NEOs: AdmissibleRegion
 A = AdmissibleRegion(sol.tracklets[1], params)
 plot(A, N = 10_000, ρscale = :log, label = "", xlabel = "log10(ρ) [au]",
-     ylabel = "v_ρ [au/day]", title = "2023 DW first tracklet AR",
-     c = :red, linewidth = 2)
-
+     ylabel = "v_ρ [au/day]", title = "2016 TU93 first tracklet AR",
+     c = :red, linewidth = 2)

scripts/adam.jl

@@ -43,21 +43,21 @@ end
 @everywhere begin
     using NEOs, Dates, TaylorSeries, PlanetaryEphemeris, JLD2
     using NEOs: RadecMPC, AdmissibleRegion, PropagationBuffer, OpticalResidual,
-        attr2bary, propres!, boundary_projection, reduce_tracklets, arboundary
+        attr2bary, propres!, boundary_projection, reduce_tracklets, arboundary,
+        indices
 
-    function adam(radec::Vector{RadecMPC{T}}, A::AdmissibleRegion{T}, ρ::T, v_ρ::T,
+    function adam(od::ODProblem{D, T}, i::Int, A::AdmissibleRegion{T}, ρ::T, v_ρ::T,
         params::NEOParameters{T}; scale::Symbol = :linear, η::T = 25.0,
-        μ::T = 0.75, ν::T = 0.9, ϵ::T = 1e-8, Qtol::T = 0.001, adamorder::Int = 2,
-        dynamics::D = newtonian!) where {T <: Real, D}
-        # Initial time of integration [julian days]
+        μ::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)
-        # Maximum number of iterations
-        maxiter = params.adamiter
+        # Unfold maximum number of iterations and relative tolerance
+        maxiter, Qtol = params.adamiter, params.adamQtol
         # Allocate memory
         aes = Matrix{T}(undef, 6, maxiter+1)
         Qs = fill(T(Inf), maxiter+1)
         # Initial attributable elements
-        aes[:, 1] .= [A.α, A.δ, A.v_α, A.v_δ, ρ, v_ρ]
+        aes[:, 1] .= A.α, A.δ, A.v_α, A.v_δ, ρ, v_ρ
         # Scaling factors
         scalings = Vector{T}(undef, 6)
         scalings[1:4] .= abs.(aes[1:4, 1]) ./ 1e6
@@ -67,46 +67,44 @@ end
             scalings[5] = (log10(A.ρ_domain[2]) - log10(A.ρ_domain[1])) / 1_000
         end
         scalings[6] = (A.v_ρ_domain[2] - A.v_ρ_domain[1]) / 1_000
-        # Jet transport variables
-        set_variables(Float64, "dx"; order = adamorder, numvars = 6)
+        # Jet transport variables and initial condition
         dae = [scalings[i] * TaylorN(i, order = adamorder) for i in 1:6]
+        AE = aes[:, 1] .+ dae
+        # Subset of radec
+        idxs = indices(od.tracklets[i])
         # Propagation buffer
-        t0, tf = dtutc2days(date(radec[1])), dtutc2days(date(radec[end]))
-        tlim = (t0 - params.bwdoffset, tf + params.fwdoffset)
-        buffer = PropagationBuffer(dynamics, jd0, tlim, aes[:, 1] .+ dae, params)
+        buffer = PropagationBuffer(od, jd0, idxs[1], idxs[end], AE, params)
         # Vector of O-C residuals
-        res = Vector{OpticalResidual{T, TaylorN{T}}}(undef, length(radec))
+        res = [zero(OpticalResidual{T, TaylorN{T}}) for _ in eachindex(idxs)]
         # Origin
-        x0 = zeros(T, 6)
-        x1 = zeros(T, 6)
+        x0, x1 = zeros(T, 6), zeros(T, 6)
         # Gradient of objective function wrt (ρ, v_ρ)
         g_t = Vector{T}(undef, 2)
         # First momentum
-        m = zeros(T, 2)
-        _m_ = zeros(T, 2)
+        m, _m_ = zeros(T, 2), zeros(T, 2)
         # Second momentum
-        n = zeros(T, 2)
-        _n_ = zeros(T, 2)
+        n, _n_ = zeros(T, 2), zeros(T, 2)
         # Gradient descent
         for t in 1:maxiter
             # Current attributable elements (plain)
             ae = aes[:, t]
             # Attributable elements (JT)
             if scale == :linear
-                AE = ae + dae
+                AE .= ae + dae
             elseif scale == :log
-                AE = [ae[1] + dae[1], ae[2] + dae[2], ae[3] + dae[3],
-                    ae[4] + dae[4], 10^(log10(ae[5]) + dae[5]), ae[6] + dae[6]]
+                AE .= [ae[1] + dae[1], ae[2] + dae[2], ae[3] + dae[3],
+                      ae[4] + dae[4], 10^(log10(ae[5]) + dae[5]), ae[6] + dae[6]]
             end
             # Barycentric state vector
             q = attr2bary(A, AE, params)
             # Propagation and residuals
             # TO DO: `ρ::TaylorN` is too slow for `adam` due to evaluations
             # within the dynamical model
-            propres!(res, radec, jd0 - ae[5]/c_au_per_day, q, params; buffer, dynamics)
+            propres!(res, od, jd0 - ae[5]/c_au_per_day, q, params; buffer, idxs)
             iszero(length(res)) && break
             # Least squares fit
             fit = tryls(res, x0, 5, 1:4)
+            !fit.success && break
             x1 .= fit.x
             # Current Q
             Q = nms(res)
@@ -177,11 +175,15 @@ function main()
     points_per_worker = [points[i] for i in idxs]
     # Evaluate `NEOs.adam` in each point
     result = pmap(points_per_worker) do points
+        # Orbit determination problem
+        od = ODProblem(newtonian!, radec)
+        # Pre-allocate output
         aes = Vector{Matrix{Float64}}(undef, length(points))
         Qs = Vector{Vector{Float64}}(undef, length(points))
+        # Main loop
         for i in eachindex(points)
             x, v_ρ = points[i]
-            aes[i], Qs[i] = adam(radec, A, 10^x, v_ρ, params; scale = :log,
+            aes[i], Qs[i] = adam(od, 1, A, 10^x, v_ρ, params; scale = :log,
                 adamorder = varorder)
             j = findlast(!isinf, Qs[i])
             if isnothing(j)

scripts/mcmov.jl

@@ -42,31 +42,35 @@ end
 @everywhere begin
     using NEOs, Dates, TaylorSeries, PlanetaryEphemeris, JLD2
     using NEOs: AdmissibleRegion, OpticalResidual, RadecMPC, PropagationBuffer,
-                reduce_tracklets, argoldensearch, scaled_variables, attr2bary, propres!
+                reduce_tracklets, argoldensearch, scaled_variables, attr2bary,
+                propres!, nobs
 
     function mcmov(points::Vector{Tuple{T, T}}, A::AdmissibleRegion{T},
         radec::Vector{RadecMPC{T}}, params::NEOParameters{T},
         varorder::Int = 2) where {T <: Real}
+        # Orbit determination problem
+        od = ODProblem(newtonian!, radec)
         # JT variables
         dae = scaled_variables("dx", ones(T, 6), order = varorder)
-        # Attributable elements
+        # Attributable elements (plain)
         ae = Vector{T}(undef, 6)
-        AE = Vector{TaylorN{T}}(undef, 6)
         ae[1:4] .= A.α, A.δ, A.v_α, A.v_δ
+        ae[5:6] .= points[1]
         # Scaling factors
         scalings = Vector{T}(undef, 6)
         scalings[1:4] .= abs.(ae[1:4]) ./ 1e6
         x_min, x_max = log10.(A.ρ_domain)
         _, y_min = argoldensearch(A, 10^x_min, 10^x_max, :min, :outer, 1e-20)
         _, y_max = argoldensearch(A, 10^x_min, 10^x_max, :max, :outer, 1e-20)
         scalings[5:6] .= (x_max - x_min) / 100, (y_max - y_min) / 100
+        # Attributable elements (jet transport)
+        AE = ae .+ scalings .* dae
+        # TDB epoch of admissible region
+        _jd0_ = dtutc2jdtdb(A.date)
         # Propagation buffer
-        t0, tf = dtutc2days(date(radec[1])), dtutc2days(date(radec[end]))
-        tlim = (t0 - params.bwdoffset, tf + params.fwdoffset)
-        buffer = PropagationBuffer(newtonian!, dtutc2jdtdb(A.date), tlim,
-                                   ae .+ scalings .* dae, params)
+        buffer = PropagationBuffer(od, _jd0_, 1, nobs(od), AE, params)
         # Vector of residuals
-        res = Vector{OpticalResidual{T, TaylorN{T}}}(undef, length(radec))
+        res = [zero(OpticalResidual{T, TaylorN{T}}) for _ in eachindex(radec)]
         # Origin
         x0 = zeros(T, 6)
         # Manifold of variations
@@ -82,9 +86,9 @@ end
             # Barycentric initial conditions (JT)
             q = attr2bary(A, AE, params)
             # Reference epoch in julian days (corrrected for light-time)
-            jd0 = dtutc2jdtdb(A.date) - constant_term(AE[5]) / c_au_per_day
+            jd0 = _jd0_ - constant_term(AE[5]) / c_au_per_day
             # Propagation and residuals
-            propres!(res, radec, jd0, q, params; buffer)
+            propres!(res, od, jd0, q, params; buffer)
             # Least squares fit
             fit = newtonls(res, x0, 10, 1:4)
             # Objective function

scripts/orbitdetermination.jl

@@ -78,8 +78,8 @@ end
         # Get at most 3 tracklets for orbit determination
         flag, radec = radecfilter(radec)
         !flag && return false
-        # Dynamical function
-        dynamics = newtonian!
+        # Orbit determination problem
+        od = ODProblem(newtonian!, radec)
         # Parameters
         params = NEOParameters(
             coeffstol = Inf, bwdoffset = 0.007, fwdoffset = 0.007,
@@ -90,7 +90,7 @@ end
         # Start of computation
         init_time = now()
         # Initial orbit determination
-        sol = orbitdetermination(radec, params; dynamics, initcond)
+        sol = orbitdetermination(od, params; initcond)
         # Termination condition
         if length(sol.res) == length(radec) && nrms(sol) < 1.5
             # Time of computation
@@ -108,7 +108,7 @@ end
             jtlsiter = 20, newtoniter = 5
         )
         # Initial orbit determination
-        _sol_ = orbitdetermination(radec, params; dynamics, initcond)
+        _sol_ = orbitdetermination(od, params; initcond)
         # Time of computation
         Δ = (now() - init_time).value
         # Choose best orbit

src/NEOs.jl

@@ -12,6 +12,7 @@ using Dates, InteractiveUtils, LazyArtifacts, LinearAlgebra, Printf, JSON,
       TaylorSeries, SatelliteToolboxTransformations, TaylorIntegration,
       SPICE, JLD2, Scratch
 using AutoHashEquals.Compat
+using Base: RefValue
 using Dates: epochms2datetime
 using Downloads: download
 using DelimitedFiles: readdlm
@@ -58,14 +59,14 @@ export loadpeeph, bwdfwdeph
 # Observer position in ECEF and ECI frames
 export obsposECEF, obsposvelECI
 # Optical astrometry processing
-export compute_radec, select_debiasing_table, debiasing, w8sveres17, residuals, unfold,
-       relax_factor, outlier
+export UniformWeights, Veres17, Farnocchia15, Eggl20, compute_radec, unfold,
+       residuals, outlier
 # Radar astrometry processing
 export compute_delay, radar_astrometry
 # Asteroid dynamical models
 export RNp1BP_pN_A_J23E_J2S_ng_eph_threads!, RNp1BP_pN_A_J23E_J2S_eph_threads!, newtonian!
 # Propagation
-export NEOParameters, propagate, propagate_lyap, propagate_root
+export propagate, propagate_lyap, propagate_root
 # Osculating
 export pv2kep, yarkp2adot
 # Least squares
@@ -77,7 +78,7 @@ export tsaiod
 # Gauss method
 export gauss_method, gaussiod
 # Orbit determination
-export curvature, issinglearc, orbitdetermination
+export NEOParameters, ODProblem, curvature, issinglearc, orbitdetermination
 # B plane
 export valsecchi_circle, bopik, mtp
 # Outlier rejection