matthewholman / assist Goto Github PK

Our conda build of sorcha is getting errors and it looks like its associated with the assist libraries.

OSError: /home/conda/feedstock_root/build_artifacts/sorcha_1702463297198/_test_env_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold/lib/python3.11/site-packages/assist/../libassist.cpython-311-x86_64-linux-gnu.so: undefined symbol: reb_error

Would you be able to take a look?

I tried to install the C version of ASSIST and I get the following errors:

https://github.com/matthewholman/assist/blob/main/README.md

Next, go to one of the example directories and compile the problem file. This will also trigger the installation of the REBOUND and ASSIST shared libraries.

cd assist/examples/plain_interface
make

The directory assist/examples/plain_interface does not exist. Instead, I find:
/Volumes/dione_raid2/Research/ASSIST/assist/examples

asteroid/ ephemeris/ interpolation/

Then, when I cd to ephemeris I get the following errors:

cd ephemeris/
ls
Makefile problem.c
make
Makefile:14: ../../..//src/Makefile.defs: No such file or directory
make: *** No rule to make target `../../..//src/Makefile.defs'. Stop.

Do you have an example (from the paper) that we can use as a unit test for the non-gravitational forces?

Note that I haven't put much thought into how to handle the parameters for the non-gravitational forces (A1 A2 A3). I think just letting the user setup an array an passing that to assist would be the easiest. So basically what's already there - with a few modifications for using it with the plain interface.

I think I see a problem with the rotation equations in the J2,J3,J4 subroutine -- lines 492-495 in
forces.c.

For starters, assume the rotations were exactly zero in both ra and dec (ie no rotation at all).
I would think that you would expect that the rotation equations would yield:
dxp = dx
dyp = dy
dzp = dz

But if you set a=0 and d=0 with cosa=1, sina=0, cosd=1, sind=0 the equations give:

dxp = dy
dyp = dz
dzp = dx

which CAN'T be right and implies that the rotation equations are wrong.

Second, you're taking the value of a and d to be the actual RA and Dec of the Earth's pole (at J2000), ie 0
and 90 respectively and not the rotation angle of the pole from the J2000 position, both of which should both
be very small angles and are actually both zero if you want to fix things at J2000. But, you're effectively saying you want to
rotate the pole through 90 degrees in Dec and 0 degrees in RA.

And a very minor point, you calculate costheta2 = dz*dz/r2. But it's actually the sine of an angle, not
a cosine. And I think it's actually the sine of the "dec" of the object (but calculated from the non-light time corrected
positions).

Larry Wasserman
Lowell Observatory
[email protected]

I'm linking to ASSIST as a library. I have my own ways I'd like to handle errors and provide warnings to the user. It'd be nice if I could disable the fprintf(stderr, ...) calls in ASSIST, maybe with a -D flag during compilation or something.

I don't really feel that I need to be able to recover those messages in any way; just suppressing them is fine, especially if there are well-defined error codes, which is mostly the case throughout ASSIST.

I'm trying to add a few more unit tests to make sure I'm not messing things up. For that, I'm converting the 5303 encounter from your notebook to a unit test. Two questions:

• Why do you turn off the various harmonics for this case? I
• It looks like the initial and final position in the notebook are coming from Horizons and our results differ by ~500m. Is that right? The paper mentions that we get to within 1m when comparing our results to the small body code. Can you give me the initial/final positions for that? That would make the test more stringent.

Right now, the version number is hardcoded in setup.py; every time one makes a release one needs to remember to manually update it.

It'd be better to use something like setuptools_scm to automatically infer the version number from git.

This issue is to make sure we don't forget about it.

You might want to add a link to the documentation, a short description, and a few keywords to the "about" section of the repository on GitHub (right side bar on the main assist page).

I get a couple of compiler warnings about the pointer arithmetic in spk.c. For examples:

spk.c:260:24: warning: arithmetic on a pointer to void is a GNU extension [-Wpointer-arith]
        val = (void *)pl->map + sizeof(double) * (pl->one[m][n] - 1)

I've fixed this in 860c72c but you've changed it back. I'm wondering why. Did that not work?

Would ASSIST be open to using de440.bsp instead of linux_p1550p2650.440 ?

I ask because Asteroid Institute has a package for DE440 (naif-de440) and is also in the process of publishing sb441-n16.bsp to a package jpl-small-bodies-de441-n16. This package is still awaiting the file size limit increase.

If ASSIST moved to using the .bsp format for both files, it would mean users could install all the ephemeris files as simple Python packages, rather than having to download them separately.

If you're open to it, I'd be happy to put in a PR with tests.

I suggest not using environment variables. I think this can be confusing to users who might not be familiar with the concept (one rarely comes across environment variables when using standard python packages).

Right now environment variables are only used in two locations to specify the path for data files. I suggest having an explicit init function that has to be called by the user before using assist. Something like

#include "assist.h"
assist_init("path/to/planet/ephem", NULL);
// or
assist_init("path/to/planet/ephem", "path/to/smallbody/ephem");

and similarly in python

import assist
assist.init("path/to/planet/ephem")
# or 
assist.init("path/to/planet/ephem", "path/to/smallbody/ephem")

If this function call results in an error for the user, then it's immediately clear what the issue is and how to fix it. It removes the stumbling block of learning what an environment variable is and how to set it.

Somewhat related to this, the metadata and the pointer to the data structures get stored in a static variable, which is effectively a global variable. If you only every read from that variable after it's initialized (which I think is the case), then the reading part is thread safe. However, the initialization is not. Imagine two threads both initializing this structure at the same time. It might work in practice, because it's a simple data structure, but there really should be a mutex around the init function if we want to claim it's thread safe (and to prevent we're accidentally loading in the data and allocating memory more than once).

Tagging the code:

git tag -a -m "Version 1.1.7" v1.1.7
git push --tags

Conda-forge "feedstock" (== the recipe on how to build the conda package) is located at https://github.com/conda-forge/assist-feedstock

Going through the constants in the JPL files, I just noticed that we're not taking account the sun's mass loss. The size of the effect is ~6e-14 per year, so it might be important for some calculations. Any thoughts?

The geocentric flag currently doesn't do anything at the moment. Do you plan to add some functionality in the future or can I remove it?

Unless I'm missing something, the user has no real control over the output cadence right now. Would it make sense to allow the user to request an output at regular intervals (similar to what NASA Horizons does). This could just be an additional option. But I think we might even want to make this the default. So instead of the user defining the number of substeps as a fraction of the timestep (which is changing because IAS15 is adaptive), the user would provide a time interval they want between outputs. This would make use of the same interpolation that is currently being used in the store_function, just evaluating it at different points. One advantage would be that the user knows exactly how much memory needs to be allocated for the output variables.

I'm happy to implement that but I want to check if that's something you think is worth having and/or if you prefer the current substeps idea.

I am seeing significantly higher roundtrip integration errors if I use the new adaptive stepping scheme (mode 2, see attached figures). The figures are generated from the Figure 2 section of the publication_figures.ipynb notebook by just setting the sim.ri_ias15.adaptive_mode to either 1 or 2. I have also implemented the new timestepping scheme in my small-body OD/propagation library and am seeing worse behavior with long-term integration of (65803) Didymos (I suspect this is just because (3666) Holman is a main-belt asteroid vs. Didymos being an NEA).

Is this an expected consequence of the new adaptive timestep algorithm? If so, @hannorein do you have recommendations for when to use each scheme? Because I am seeing great stability for integrations with close approaches (e.g., Apophis) when I slightly modify the minimum timestep with the new algorithm (the previous one would oscillate significantly in those scenarios). So, my current intuition dictates that if I am expecting to propagate deep encounters/impacts, I should use the new timestep algorithm, but not in other cases such as long integrations (~100 years). Although I'm not sure what to prefer if I don't have this information beforehand...

When writing a shared library, people tend to add a prefix to all non-static variables and functions to prevent collisions between symbols with the same name from different libraries. If for example rebound and assist were both implementing a function with the name heartbeat, then there would be a conflict. At the moment a lot of symbols in assist do not have the assist_ prefix. You can check which ones by running:

nm -g --defined-only libassist.so

All the symbols that show up when running this command should start with assist_.

../../../../../../../../mambaforge/envs/sorcha/lib/python3.10/site-packages/assist/init.py:44
/Users/mschwamb/mambaforge/envs/sorcha/lib/python3.10/site-packages/assist/init.py:44: DeprecationWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html
import pkg_resources

We should make the PyPi (and conda-forge) publishing be automatically triggered from a CI step in github, once the repo is tagged with something like vX.Y.Z.

It would ensure pip+conda releases are always consistent with the code on github (and make the whole release process much easier -- just tag, push & let the CI do the rest).

Originally posted by @mjuric in #94 (comment)

Since the rebound update to 4.0.3, we noticed some of our orbits were taking a long time to run, or in some cases hanging indefinitely. Pinning the version to 4.0.2 resolves the issue. The following code recreates the behavior:

import rebound
import assist

this_breaks = True

print(f"rebound version: {rebound.__version__}")

ephem = assist.Ephem("../data/linux_p1550p2650.440", "../data/sb441-n16.bsp")

sim = rebound.Simulation()

sim.t = 14658.412205805536

if this_breaks:
    sim.add(
        x=-0.9716337587843651,
        y=0.5557887162674847,
        z=0.15913299510349357,
        vx=-0.004071385625678586,
        vy=-0.016941109480633685,
        vz=-0.004809720968068568,
    )
else:
    sim.add(
        x=-0.971634,
        y=0.555789,
        z=0.159133,
        vx=-0.004071,
        vy=-0.016941,
        vz=-0.00481,
    )

ax = assist.Extras(sim, ephem)

sim.integrate(14718.412205805536)

I don't know if it's helpful or not, but the truncated values do not produce the same error, which is why I was having trouble recreating the issue.

I'm looking into the jd_ref thing. But I'm confused. I'm either missing something (most likely) or there is a bug.

In the examples, both jd_ref and tstart are full (not relative) Julian dates. The rebound variable r->t is being set to tstart, so also a full Julian date. But then this line doesn't make sense. Following through the various function calls, I believe the code is just adding jd_ref to r->t here and then comparing it to pl->beg and pl->end. All of these are in full Julian dates so it doesn't work. It would only work if r->t is a relative Julian date.

Maybe this specific line didn't trigger an error because the code was never tested with dates before 11 Jan 17000 CE. But I'm also not sure that the lines below are correct where the record number blk is calculated. But that would have surely shown up on a test, so I assume that part must be correct. So I must be missing something.

The documentation is still a bit sparse. I don't think we need to add much more. The examples show most of the functionality. However, would you have any objections to me converting the documentation from sphinx to mkdocs? Mkdocs has the nice feature which lets you put C and python code next to each other, so we don't need two different pages explaining the same feature in two different languages?

I'm getting errors when trying to install assist with pip and Python 3.9.18 | packaged by conda-forge | (main, Aug 30 2023, 03:53:08)

mschwamb@cider:~$ pip install assist
Collecting assist
Downloading assist-1.1.4.tar.gz (47 kB)
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 47.0/47.0 kB 769.9 kB/s eta 0:00:00
Installing build dependencies ... done
Getting requirements to build wheel ... done
Installing backend dependencies ... done
Preparing metadata (pyproject.toml) ... done
Collecting rebound>=3.10.0 (from assist)
Using cached rebound-4.0.1-cp39-cp39-macosx_11_0_arm64.whl.metadata (9.7 kB)
Collecting numpy (from assist)
Using cached numpy-1.26.2-cp39-cp39-macosx_11_0_arm64.whl.metadata (61 kB)
Using cached rebound-4.0.1-cp39-cp39-macosx_11_0_arm64.whl (233 kB)
Using cached numpy-1.26.2-cp39-cp39-macosx_11_0_arm64.whl (14.0 MB)
Building wheels for collected packages: assist
Building wheel for assist (pyproject.toml) ... error
error: subprocess-exited-with-error

× Building wheel for assist (pyproject.toml) did not run successfully.
│ exit code: 1
╰─> [58 lines of output]
fatal: not a git repository (or any of the parent directories): .git
running bdist_wheel
running build
running build_py
creating build
creating build/lib.macosx-11.1-arm64-cpython-39
creating build/lib.macosx-11.1-arm64-cpython-39/assist
copying assist/extras.py -> build/lib.macosx-11.1-arm64-cpython-39/assist
copying assist/init.py -> build/lib.macosx-11.1-arm64-cpython-39/assist
copying assist/ephem.py -> build/lib.macosx-11.1-arm64-cpython-39/assist
copying assist/py.typed -> build/lib.macosx-11.1-arm64-cpython-39/assist
running build_ext
building 'libassist' extension
creating build/temp.macosx-11.1-arm64-cpython-39
creating build/temp.macosx-11.1-arm64-cpython-39/src
clang -Wno-unused-result -Wsign-compare -Wunreachable-code -DNDEBUG -fwrapv -O2 -Wall -fPIC -O2 -isystem /Users/mschwamb/mambaforge/envs/assist/include -arch arm64 -I/Users/mschwamb/mambaforge/envs/assist/include -fPIC -O2 -isystem /Users/mschwamb/mambaforge/envs/assist/include -arch arm64 -DLIBASSIST -Isrc -I/Users/mschwamb/mambaforge/envs/assist/include/python3.9 -I/private/var/folders/ss/8zqns7xj6jd1kj4q3y81kg500000gn/T/pip-build-env-5suqi78x/normal/lib/python3.9/site-packages/rebound -c src/assist.c -o build/temp.macosx-11.1-arm64-cpython-39/src/assist.o -fstrict-aliasing -O3 -std=c99 -fPIC -D_GNU_SOURCE -Wpointer-arith -DASSISTGITHASH=76ae72ccd4dc2b706ce638a30bb59f597476f717
src/assist.c:195:32: warning: comparison of integers of different signs: 'int' and 'unsigned long' [-Wsign-compare]
for (int m=0; m<sizeof(massmap); m+=2){
~^~~~~~~~~~~~~~~~
src/assist.c:363:9: error: call to undeclared function 'reb_error'; ISO C99 and later do not support implicit function declarations [-Wimplicit-function-declaration]
reb_error(assist->sim, msg);
^
src/assist.c:452:36: error: variable has incomplete type 'enum reb_input_binary_messages'
enum reb_input_binary_messages warnings = REB_INPUT_BINARY_WARNING_NONE;
^
src/assist.c:452:10: note: forward declaration of 'enum reb_input_binary_messages'
enum reb_input_binary_messages warnings = REB_INPUT_BINARY_WARNING_NONE;
^
src/assist.c:453:33: error: call to undeclared function 'reb_create_simulation'; ISO C99 and later do not support implicit function declarations [-Wimplicit-function-declaration]
struct reb_simulation* r2 = reb_create_simulation();
^
src/assist.c:453:28: warning: incompatible integer to pointer conversion initializing 'struct reb_simulation ' with an expression of type 'int' [-Wint-conversion]
struct reb_simulation r2 = reb_create_simulation();
^ ~~~~~~~~~~~~~~~~~~~~~~~
src/assist.c:454:5: error: call to undeclared function 'reb_create_simulation_from_simulationarchive_with_messages'; ISO C99 and later do not support implicit function declarations [-Wimplicit-function-declaration]
reb_create_simulation_from_simulationarchive_with_messages(r2, sa, blob-1, &warnings);
^
src/assist.c:454:5: note: did you mean 'reb_simulation_create_from_simulationarchive_with_messages'?
/private/var/folders/ss/8zqns7xj6jd1kj4q3y81kg500000gn/T/pip-build-env-5suqi78x/normal/lib/python3.9/site-packages/rebound/rebound.h:1231:16: note: 'reb_simulation_create_from_simulationarchive_with_messages' declared here
DLLEXPORT void reb_simulation_create_from_simulationarchive_with_messages(struct reb_simulation* r, struct reb_simulationarchive* sa, int64_t snapshot, enum reb_simulation_binary_error_codes* warnings);
^
src/assist.c:457:28: warning: incompatible integer to pointer conversion initializing 'struct reb_simulation ' with an expression of type 'int' [-Wint-conversion]
struct reb_simulation r3 = reb_create_simulation();
^ ~~~~~~~~~~~~~~~~~~~~~~~
src/assist.c:463:5: error: call to undeclared function 'reb_free_simulation'; ISO C99 and later do not support implicit function declarations [-Wimplicit-function-declaration]
reb_free_simulation(r3);
^
src/assist.c:490:9: error: call to undeclared function 'reb_integrate'; ISO C99 and later do not support implicit function declarations [-Wimplicit-function-declaration]
reb_integrate(sim, t);
^
src/assist.c:579:5: error: call to undeclared function 'reb_update_acceleration'; ISO C99 and later do not support implicit function declarations [-Wimplicit-function-declaration]
reb_update_acceleration(sim); // This will later be recalculated. Could be optimized.
^
src/assist.c:452:47: error: use of undeclared identifier 'REB_INPUT_BINARY_WARNING_NONE'
enum reb_input_binary_messages warnings = REB_INPUT_BINARY_WARNING_NONE;
^
3 warnings and 8 errors generated.
error: command '/usr/bin/clang' failed with exit code 1
[end of output]

note: This error originates from a subprocess, and is likely not a problem with pip.
ERROR: Failed building wheel for assist
Failed to build assist
ERROR: Could not build wheels for assist, which is required to install pyproject.toml-based projects

In some cases, it's not necessary for my purposes to include asteroid perturbers in the force model. But assist_ephem_init still requires a full small bodies ephemeris file, even if I don't intend to use it.

That's a consequential requirement because the SB441-N16 file is over 600MB, which is enough that (for example) it significantly inflates application runtime container sizes. That makes application launch time slower and causes unnecessary network traffic.

So it would be nice if I could use assist without any small bodies ephemeris file. I'm assuming that if "ASTEROIDS" is not in my set of forces, that file is not necessary, which could be wrong, though.