Skip to content

Commit 8f1a664

Browse files
authored
PropagationBuffer (#68)
* Add OhMyThreads compat entry * Multithreaded evaleph * PropagationBuffer * DynamicalParameters * Update jetcoeffs.jl * Remove taylorize macro from newtonian * Adapt propagate_lyap to PropagationBuffer * In-place propres * Fix rvelea * _gausstriplets1/2 * Jet Transport Least Squares * Fix tests * Fix residuals * Rename orbit_determination -> orbitdetermination * Update orbitdetermination functions * Orbit Determination script * Names script * Minor fixes * Incremental radec in jtls * _adam! helps gauss jtls * epoch(::NEOSolution) * Do not use jplcompare in OD tests * Optional ADAM helps Gauss * Update scripts/orbitdetermination.jl * Update scripts/orbitdetermination.jl * Bump patch version * Minor fix * Suggestions by @lbenet
1 parent 6c9df93 commit 8f1a664

21 files changed

+2914
-2145
lines changed

Project.toml

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
name = "NEOs"
22
uuid = "b41c07a2-2abb-11e9-070a-c3c1b239e7df"
33
authors = ["Jorge A. Pérez Hernández", "Luis Benet", "Luis Eduardo Ramírez Montoya"]
4-
version = "0.8.5"
4+
version = "0.8.6"
55

66
[deps]
77
Artifacts = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
@@ -20,6 +20,7 @@ JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
2020
LazyArtifacts = "4af54fe1-eca0-43a8-85a7-787d91b784e3"
2121
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
2222
LsqFit = "2fda8390-95c7-5789-9bda-21331edee243"
23+
OhMyThreads = "67456a42-1dca-4109-a031-0a68de7e3ad5"
2324
PlanetaryEphemeris = "d83715d0-7e5f-11e9-1a59-4137b20d8363"
2425
Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
2526
Quadmath = "be4d8f0f-7fa4-5f49-b795-2f01399ab2dd"
@@ -55,6 +56,7 @@ JSON = "0.21"
5556
LazyArtifacts = "1"
5657
LinearAlgebra = "1"
5758
LsqFit = "0.15"
59+
OhMyThreads = "0.5"
5860
PlanetaryEphemeris = "0.8"
5961
Printf = "1"
6062
Quadmath = "0.5"

scripts/names.jl

Lines changed: 106 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,106 @@
1+
using ArgParse, Downloads, Random, HORIZONS, NEOs
2+
using NEOs: numberofdays
3+
4+
# Potentially Hazardous Asteroids MPC list
5+
const PHAS_URL = "https://cgi.minorplanetcenter.net/cgi-bin/textversion.cgi?f=lists/PHAs.html"
6+
# Amors MPC list
7+
const AMORS_URL = "https://cgi.minorplanetcenter.net/cgi-bin/textversion.cgi?f=lists/Amors.html"
8+
# Apollos MPC list
9+
const APOLLOS_URL = "https://cgi.minorplanetcenter.net/cgi-bin/textversion.cgi?f=lists/Apollos.html"
10+
# Atens MPC list
11+
const ATENS_URL = "https://cgi.minorplanetcenter.net/cgi-bin/textversion.cgi?f=lists/Atens.html"
12+
13+
function parse_commandline()
14+
s = ArgParseSettings()
15+
16+
# Program name (for usage & help screen)
17+
s.prog = "names.jl"
18+
# Desciption (for help screen)
19+
s.description = "Select a list of NEOs for orbitdetermination.jl"
20+
21+
@add_arg_table! s begin
22+
"--N", "-N"
23+
help = "Number of NEOs"
24+
arg_type = Int
25+
default = 100
26+
"--output", "-o"
27+
help = "Output file"
28+
arg_type = String
29+
default = "names.txt"
30+
end
31+
32+
s.epilog = """
33+
examples:\n
34+
\n
35+
julia --project names.jl -N 250 -o names.txt\n
36+
\n
37+
"""
38+
39+
return parse_args(s)
40+
end
41+
42+
function main()
43+
44+
# Parse arguments from commandline
45+
parsed_args = parse_commandline()
46+
47+
println("NEOs selector for orbitdetermination.jl")
48+
49+
# Number of NEOs
50+
N::Int = parsed_args["N"]
51+
println("• Number of NEOs: ", N)
52+
# Output file
53+
outfile::String = parsed_args["output"]
54+
55+
# Parse PHAs, Atens, Apollos and Amors
56+
provdesig = Vector{SubString{String}}(undef, 0)
57+
58+
for url in [PHAS_URL, ATENS_URL, APOLLOS_URL, AMORS_URL]
59+
# Download MPC file
60+
Downloads.download(url, "tmp.txt")
61+
# Parse lines
62+
lines = readlines("tmp.txt")[3:end]
63+
filter!(!isempty, lines)
64+
filter!(l -> isempty(strip(l[1:27])), lines)
65+
# Parse provisional designations
66+
provdesig = vcat(provdesig, map(l -> strip(l[28:40]), lines))
67+
# Delete MPC file
68+
rm("tmp.txt")
69+
end
70+
# Delete repeated NEOs
71+
unique!(provdesig)
72+
# We can only process asteroids discovered between 2000 and 2024
73+
filter!(desig -> "2000" <= desig[1:4] <= "2024", provdesig)
74+
# Shuffle the provisional designations list
75+
shuffle!(provdesig)
76+
77+
# Assemble a list with N NEOs
78+
names = Vector{String}(undef, N)
79+
80+
for i in eachindex(names)
81+
while !isempty(provdesig)
82+
neo = String(pop!(provdesig))
83+
radec = fetch_radec_mpc("designation" => neo)
84+
jplorbit = sbdb("des" => neo)["orbit"]
85+
if (numberofdays(radec) <= 30.0) && (jplorbit["n_obs_used"] == length(radec)) &&
86+
isnothing(jplorbit["n_del_obs_used"]) && isnothing(jplorbit["n_dop_obs_used"])
87+
names[i] = neo
88+
break
89+
end
90+
end
91+
end
92+
sort!(names)
93+
94+
# Save names list
95+
open(outfile, "w") do file
96+
for i in eachindex(names)
97+
write(file, names[i], i == N ? "" : "\n")
98+
end
99+
end
100+
101+
println("• Saved ", N, " names to ", outfile)
102+
103+
nothing
104+
end
105+
106+
main()

scripts/orbitdetermination.jl

Lines changed: 208 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,208 @@
1+
using Distributed, ArgParse
2+
3+
function parse_commandline()
4+
s = ArgParseSettings()
5+
6+
# Program name (for usage & help screen)
7+
s.prog = "orbitdetermination.jl"
8+
# Desciption (for help screen)
9+
s.description = """
10+
Determines the orbit of a list of NEOs via jet transport
11+
using Julia's interface for distributed computing."""
12+
13+
@add_arg_table! s begin
14+
"--names", "-n"
15+
help = "File containing the names of the NEOs to be processed"
16+
arg_type = String
17+
default = "names.txt"
18+
"--output", "-o"
19+
help = "Output directory"
20+
arg_type = String
21+
default = pwd()
22+
end
23+
24+
s.epilog = """
25+
example:\n
26+
\n
27+
# Run orbitdetermination.jl with 10 workers and 5 threads each\n
28+
julia -p 10 -t 5 --project scripts/orbitdetermination.jl -n names.txt -o orbits/\n
29+
\n
30+
"""
31+
32+
return parse_args(s)
33+
end
34+
35+
@everywhere begin
36+
using NEOs, Dates, JLD2
37+
using NEOs: Tracklet, reduce_tracklets, isgauss
38+
39+
# Initial orbit determination routines
40+
41+
# First filter
42+
function iod1(neo::String, outdir::String)
43+
# Output file
44+
filename = joinpath(outdir, replace(neo, " " => "_") * ".jld2")
45+
# Download optical astrometry
46+
radec = fetch_radec_mpc("designation" => neo)
47+
if length(radec) < 3
48+
jldsave(filename; tracklets = Vector{Tracklet}(undef, 0))
49+
return false
50+
end
51+
# Parameters
52+
params = NEOParameters(coeffstol = 10.0, bwdoffset = 0.007,
53+
fwdoffset = 0.007, jtlsiter = 10, adamhelp = false)
54+
# Select at most three tracklets
55+
tracklets = reduce_tracklets(radec)
56+
if length(tracklets) > 3
57+
tracklets = tracklets[1:3]
58+
radec = reduce(vcat, getfield.(tracklets, :radec))
59+
sort!(radec)
60+
end
61+
62+
try
63+
# Start of computation
64+
init_time = now()
65+
# Orbit determination
66+
sol = orbitdetermination(radec, params)
67+
# Time of computation
68+
Δ = (now() - init_time).value
69+
# Save orbit
70+
if length(sol.res) != length(radec)
71+
jldsave(filename; tracklets = tracklets, Δ = Δ)
72+
return false
73+
else
74+
jldsave(filename; sol = sol, Δ = Δ)
75+
return true
76+
end
77+
catch
78+
# An error ocurred
79+
jldsave(filename; tracklets = tracklets)
80+
return false
81+
end
82+
end
83+
84+
# Second filter
85+
function iod2(neo::String, outdir::String)
86+
# Output from last filter
87+
filename = joinpath(outdir, replace(neo, " " => "_") * ".jld2")
88+
dict = JLD2.load(filename)
89+
# Previous filter already computed an orbit
90+
"sol" in keys(dict) && return true
91+
# Check tracklets are non empty
92+
tracklets = dict["tracklets"]
93+
isempty(tracklets) && return false
94+
# Computation time so far
95+
if "Δ" in keys(dict)
96+
Δ = dict["Δ"]
97+
else
98+
Δ = 0
99+
end
100+
# Optical astrometry
101+
radec = reduce(vcat, getfield.(tracklets, :radec))
102+
sort!(radec)
103+
# Parameters
104+
params = NEOParameters(coeffstol = Inf, bwdoffset = 0.5,
105+
fwdoffset = 0.5, jtlsiter = 10, adamhelp = true)
106+
107+
try
108+
# Start of computation
109+
init_time = now()
110+
# Orbit determination
111+
sol = orbitdetermination(radec, params)
112+
# Time of computation
113+
Δ += (now() - init_time).value
114+
# Save orbit
115+
if length(sol.res) != length(radec)
116+
jldsave(filename; tracklets = tracklets, Δ = Δ)
117+
return false
118+
else
119+
jldsave(filename; sol = sol, Δ = Δ)
120+
return true
121+
end
122+
catch
123+
# An error ocurred
124+
jldsave(filename; tracklets = tracklets)
125+
return false
126+
end
127+
end
128+
129+
# Third filter
130+
function iod3(neo::String, outdir::String)
131+
# Output from last filter
132+
filename = joinpath(outdir, replace(neo, " " => "_") * ".jld2")
133+
dict = JLD2.load(filename)
134+
# Previous filter already computed an orbit
135+
"sol" in keys(dict) && return true
136+
# Check tracklets are non empty
137+
tracklets = dict["tracklets"]
138+
isempty(tracklets) && return false
139+
# Computation time so far
140+
if "Δ" in keys(dict)
141+
Δ = dict["Δ"]
142+
else
143+
Δ = 0
144+
end
145+
# Optical astrometry
146+
radec = reduce(vcat, getfield.(tracklets, :radec))
147+
sort!(radec)
148+
# Parameters
149+
params = NEOParameters(coeffstol = Inf, bwdoffset = 0.5,
150+
fwdoffset = 0.5, jtlsiter = 10, adamhelp = true)
151+
152+
try
153+
# Start of computation
154+
init_time = now()
155+
# Orbit determination
156+
if isgauss(tracklets)
157+
sol = tooshortarc(radec, tracklets, params)
158+
else
159+
sol = gaussinitcond(radec, tracklets, params)
160+
end
161+
# Time of computation
162+
Δ += (now() - init_time).value
163+
# Save orbit
164+
if length(sol.res) != length(radec)
165+
jldsave(filename; tracklets = tracklets, Δ = Δ)
166+
return false
167+
else
168+
jldsave(filename; sol = sol, Δ = Δ)
169+
return true
170+
end
171+
catch
172+
# An error ocurred
173+
jldsave(filename; tracklets = tracklets)
174+
return false
175+
end
176+
end
177+
end
178+
179+
function main()
180+
# Parse arguments from commandline
181+
parsed_args = parse_commandline()
182+
183+
println("Orbit determination for NEOs via jet transport")
184+
println()
185+
186+
# NEOs to be processed
187+
namesfile = parsed_args["names"]
188+
neos = readlines(namesfile)
189+
println(length(neos), " NEOs to be processed with ", nworkers(), " workers (",
190+
Threads.nthreads(), " threads each)")
191+
println()
192+
# Output directory
193+
outdir = parsed_args["output"]
194+
195+
# Distributed orbit determination
196+
197+
# First filter
198+
mask = pmap(neo -> iod1(neo, outdir), neos)
199+
all(mask) && return nothing
200+
# Second filter
201+
mask = pmap(neo -> iod2(neo, outdir), neos)
202+
all(mask) && return nothing
203+
# Third filter
204+
mask = pmap(neo -> iod3(neo, outdir), neos)
205+
return nothing
206+
end
207+
208+
main()

0 commit comments

Comments
 (0)