SPICE only implementation #110
Conversation
|
Regarding the machine dependent results you observe, this depends on the compiler. You might need to turn off floating point contractions. I don't know what compiler you're using. If it's clang, use the |
Thanks @hannorein . I'll investigate that. |
akoumjian
force-pushed
the
ak/spk_only
branch
from
2187131
to
dbc4ed7
Compare
akoumjian
force-pushed
the
ak/spk_only
branch
from
871da37
to
ce03e94
Compare
|
|
||
| # Iterate through each directory in the unit_tests directory and compile the test files | ||
| # then run them | ||
| test: |
There was a problem hiding this comment.
I wanted an easy way to run all the c unit tests locally.
| To use use ASSIST, you also need to download ephemeris data files. One file for planet ephemeris and another suplementary file for asteroid ephemeris. You can do this with Python packages or by downloading files directly. Note that these are large files, almost 1GB in size. | ||
|
|
||
| pip install naif-de440 | ||
| pip install jpl-small-bodies-de441-n16 | ||
|
|
||
|
|
||
| You can initialize the ephemeris data from the packages like so: | ||
|
|
||
| python3 | ||
|
|
||
| >>> import assist | ||
| >>> from naif_de440 import de440 | ||
| >>> from jpl_small_bodies_de441_n16 import de441_n16 | ||
| >>> ephem = assist.Ephem(de440, de441_n16) |
There was a problem hiding this comment.
I hope it's okay that I demonstrate how the python packages can be used in lieu of downloading the files from JPL.
| printf("Difference: %.2f m\n",diff); | ||
| assert(diff < 250); | ||
| fprintf(stderr,"diff to JPL Small Body code %fm\n",diff); | ||
| assert(diff < 200); |
There was a problem hiding this comment.
This particular moment in time gets closer to JPL with SPICE ephem.
| @@ -51,7 +51,7 @@ int main(int argc, char* argv[]){ | |||
|
|
|||
| // Compare results | |||
| for (int i=0; i<Noutput; i++){ | |||
| assert(fabs(x_direct[i]-x_interpolate[i]) < 1e-13); // Check that results agree | |||
| assert(fabs(x_direct[i]-x_interpolate[i]) < 2e-13); // Check that results agree | |||
There was a problem hiding this comment.
I believe the switch to SPICE is making the interpolate function agree ever so slightly less. After doing many comparisons with JPL Horizons and seeing equal or better residuals (compared to the .440 version of ASSIST), I'm inclined to believe this has to do more with interpolation than it does less accuracy with the integration.
|
@matthewholman @cylon359 I think this is ready for your review. If it does get merged, it seems like a major version increment is warranted given the lack of backwards compatibility with the other data files. |
|
What was the reason to remove support for the other data files? |
|
Rob and Matt suggested that maintaining two distinct paths for loading the coefficients wasn't necessary, which is why I came back with this simpler implementation (prior PR supported both implementations). There was also indication that the ASCII-derived files are missing a couple bits of data, likely due to the serialization process, though there's no knowing if those bits have meaningful precision or not. |
|
I see. Oh well. I would have voted to keep backwards compatibility. 😉 |
|
I'm happy to accommodate whatever consensus is reached. If I do go back to the multiple file implementation, I will need to add back a couple small improvements and I would want to more aggressively re-organize how the code runs. Initially I wanted to modify as little as possible being a new contributor. |
|
IMO, the old implementation could just live in a separate "legacy" file with the intend to not update it anymore or provide much support. |
|
It's not quite that straight forward, as there are two slightly different interfaces for how the masses are loaded and the coefficients generated. I would probably update both to have a similar or identical interface so that all the places that differ in |
|
Yeah. I see there is a lot change in the backend. What was the idea behind using a getter function for the constants rather than just populating constants in a struct like it used to be? It probably doesn't matter too much for the run time, but now every time a constant is used in the code, it might require up to 8 string comparisons just to read the value. |
|
Yes, that was left over from the two implementations version. I was wondering if it might impact performance and was surprised that it benchmarked faster. Perhaps I could get it even faster removing those comparisons. |
|
It probably doesn't matter much for the runtime, just seems a bit unnecessary. In fact, I would remove Also, is the data in |
|
I'm definitely open to that. If we maintain the file type readers, I would consolidate both implementations into the assist_ephem struct. |
|
FYI: I'm also happy to put in some clean-up work. (I just don't want to mess up your pull request too much!) It looks like the actual new "datafile reading part" is great and doesn't need any changes, so the precision tests should not be affected by any of the more or less stylistic code changes. |
|
If you want to branch off this, or merge it in, then clean it up, all are fine with me. I would like to get a version published to pypi that supports this format soon, as it will make a lot of our pipelines cleaner and we won't need to download multiple versions of the coefficient files. I won't have much time this next week to put towards it, but if you'd like me to make those changes I can do it later. |
|
I'm also pretty busy this and (most of) next week. If you want to use it in your pipeline I wouldn't worry too much about putting it on pypi, just install it from your copy. Then you're also safe from any sudden changes which would break the pipeline and can update manually when a new feature is added to the upstream version that you need. But I have no objections to merging this PR and then later changing things (it's just more likely to make other users unhappy if it breaks something for them). |
Update ---
I’ve finally had a chance to do the analysis to compare this version of ASSIST (SPICE) with the current implementation (ASCII derived binary files). I wanted to do this analysis to make sure that my changes did not negatively impact the accuracy or performance.
I have run two similar types of tests. The first compares against future state vectors provided by JPL Horizons, the second compares against real historical observations in the MPC (using our ephemeris generation on top of ASSIST).
In addition to testing spk vs ascii, I also wanted to include comparisons for ephemeris version (440 vs 441), FFP contraction (on vs off), and coefficient truncation (50 bits vs the full 52). The latter is only tested for the SPICE implementation and only because there were some indications that the truncated values were aligning closer to JPL.
Since we are analyzing multiple objects, and each object has a different number of dates we are comparing, I’m providing two numbers to represent the performance.
Mean residual – For a given environment, object result pair, we take the mean of all residuals. We then take the mean across all objects for each environment. Smaller is better.
Inverse, normalized residual – Take the worst mean residual for each environment, object pair. Normalize all other environment, object pair mean residuals against it (e.g. the worst performer gets a 1, an environment with half the residuals compared to the worst would be 0.5). Then take the 1 - normalized residual. Lastly, take the mean of those inverse, normalized residuals across all objects for each distinct environment. Since this is an inverse value, larger is better.
Summary of Results
Based on how small the differences are, how main and spk both perform slightly better in different scenarios, I am comfortable making this transition. To me the open question is whether or not to truncate the values. I am inclined not to, as it uses extra computation and we don’t know if it’s lost in the noise.
Minor Planet Center
The main branch performs better, on order of 0.05 milliarcseconds against the spk branch. I am not sure about the statistical significance. Interesting that spk with truncation also aligns closer to MPC and that the main branch is also using truncated coefficients via the ASCII file.
JPL Horizons
The SPICE branch with truncation performs best, on the order of 5-10 meters. Once again, I don’t see any statistical significance. Interesting again that the spice implementation performs better with coefficient truncation.
I’m using a branch of our adam-assist tool to do the analysis. You can find the relevant files here:
https://github.com/B612-Asteroid-Institute/adam-assist/blob/ak/test-assist-versions/tests/test_complete_residuals.py
https://github.com/B612-Asteroid-Institute/adam-assist/blob/ak/test-assist-versions/test-assist-versions.sh
Performance
I ran the benchmark test but updated the sample size to 1000, here are the results on my Macbook. I wasn't expecting it to be faster, but it is slightly.
Bug in Loading DE441
This PR fixes a bug in the SPICE kernel loading. SPICE files do not require that summary records for the same target are adjacent to each other. The previous
assist_spk_initwould miss the second string of records in DE441.As per meeting with @matthewholman and @cylon359 ...
This PR moves to support only the SPICE (.bsp) files for planet ephem. Most of the significant changes are the same as #107 . References to the
.440ASCII derived file have been replaced.I have removed the truncation as we no longer need to find agreement with the ASCII derived files. I still intend to determine if truncation provides better agreement with Horizons on a wide scale, as Davide seemed rightly skeptical.
This PR remains as Work in Progress as I am attempting to determine why the
onthefly_interpolation,onthefly_backwards_interpolationunit tests are failing on Github CI. I did notice that these tests produce different results when running on Apple M machines vs the Github X86_64 VMs. The difference is enough to cause one test to pass on Apple but fail on x86_64.I have been able to narrow the difference down to the step size being used by
rebound. I'm wondering if the higher precision from the SPICE ephem is triggering different behavior on some maths happening in rebound (but only for certain architectures / libraries).