How to add specific conversions?

[nstarman]

In galax we have implemented many specific conversions, for example both the gala to galax and the galax to gala conversions, as done here.

The question is what to do with these functions? Should everything be moved here? Or should each galactic dynamics library
host the conversion functions from the other libraries to their own library?
This is partly a question about how best to support long-term compatibility, when one library updates and changes some definitions. Does every other library then have to cut compatibility releases? Or if it's hosted here, then only this package needs to cut a new release.

Thoughts? @jobovy @adrn @eugvas

Working Overview:

• Hosted here: the conversion code is hosted here, outside of the libraries themselves, to prevent each library having take on the tech debt of having to support all the other libraries.
• support policy: committing to some NEP-29-style support policy, to support conversions going back a few versions for each library.
• Reference Implementation: define a standard implementation here for each potential and then have each library convert to and from that.
  • A schema defining each potential
  • A data structure, following the schema, to which each library's potentials can be converted / read from.

[jobovy]

My two cents: I think it would be ideal if all conversion code is hosted outside of the libraries themselves, e.g., here. I think it's easiest for maintaining each library to not take on the tech debt of having to support all the other libraries. But I think we'd have to commit to supporting some sort of NEP-29-style support policy, to support conversions going back a few versions for each library. But again, that would then be effectively done through tests hosted in this repository. We could workflow_dispatch builds here ever time one of the supported libraries pushes to their main to maintain compatibility.

Even with just four libraries, the combinatorics of implementing all library1 <-> library2 conversions separately are painful, so ideally we'd not do this (especially if you think we'd be adding even more libraries). There are a few options:

• Define a standard implementation here for each potential and then have each library convert to and from that, so all conversions can run through it --> lot of work implementing all standard versions
• For each potential, pick one library's version as the standard one and then have every other library convert to and from that. This means that we don't have to implement each potential again here and easily allows libraries not to implement all potentials.

[nstarman]

I think it would be ideal if all conversion code is hosted outside of the libraries themselves, e.g., here. I think it's easiest for maintaining each library to not take on the tech debt of having to support all the other libraries. But I think we'd have to commit to supporting some sort of NEP-29-style support policy, to support conversions going back a few versions for each library. But again, that would then be effectively done through tests hosted in this repository. We could workflow_dispatch builds here ever time one of the supported libraries pushes to their main to maintain compatibility.

💯. I totally agree. I really like the NEP-29 backwards compatibility idea.

Even with just four libraries, the combinatorics of implementing all library1 <-> library2 conversions separately are painful, so ideally we'd not do this (especially if you think we'd be adding even more libraries). There are a few options:

• For each potential, pick one library's version as the standard one and then have every other library convert to and from that. This means that we don't have to implement each potential again here and easily allows libraries not to implement all potentials.

I like this idea in the abstract, but worry that there's no ideal library to choose as the standard. The ideal standard library is

1. pure python, so never has C compilation requirements
2. is in some framework all the libraries support, so jax is kinda iffy

• Define a standard implementation here for each potential and then have each library convert to and from that, so all conversions can run through it --> lot of work implementing all standard versions

I actually think this might be the more feasible! We could define a schema for each potential, with a mathematical definition of the potential (and density) given the parameters in the schema. I would advocate using sympy in some way for the schema...

PotentialSchema(NamedTuple):
    doc: str
    constants:  tuple[Symbol, ...] = ...
    parameters: tuple[Symbol, ...] = ...
    variables: tuple[Symbol, ...] = ...
    potential: Expression = ...

kepler = PotentialSchema(
    """The Kepler Potential.

     Parameters:
         G: Gravitational constant
         M: total mass
         r:  spherical radius
     """
    (<symbols>),
    (<symbols>),
    (<symbols>),
    - G M / r^2
)
...

With sympy libraries could programatically create utilities for conversion...

Each library could then support a conversion function that provides the parameter values.

>>> from galactic_dynamics_interop import to_schema_dict
>>> import galax.potential as gp
>>> to_schema_dict(gp.KeplerPotential(...))
{"type": "kepler", "constants: {"G": ...}, "parameters": {"M": ...}}   # note: no variables.

[eugvas]

can you more clearly define the scope of this package? I could think of two rather different possibilities:

1. have a universal backend for computing potentials and their derivatives, which would be used by all dynamics codes; or
2. have a standardized interface for defining the parameters of the potential (and/or density, for that matter), which all codes understand.

Now, (1) seems to be out of question for the reasons you listed above – impossible to find the common language base. Pure Python is slow, C/C++ is not JAXable, and vice versa JAX is not easy to port into the C++ workflow (I haven't tried though, but in any case it's a separate language, adding an unwelcome dependency).
(2) seems more realistic, but also with some caveats. It's straightforward to have a common ground for simple models (e.g. Plummer – just two parameters), not so much for complicated cases such as a multicomponent evolving MW + moving LMC expressed in some time-dependent basis sets.
What shall we aim for when designing such a partially supported interface?

[adrn]

I agree with @jobovy and what has been said already! (@nstarman and I talked about this a bit earlier as well) I think we're talking about a small enough number of potential models that it's not worth the complexity of routing everything through sympy, so I would prefer to keep the translations more explicit for now.

[nstarman]

Ok. If we drop the sympy then the schema idea becomes:

LatexStr = NewType(str)

PotentialSchema(NamedTuple):
    description: str
    constants:  tuple[str, ...] = ...
    parameters: tuple[str, ...] = ...
    variables: tuple[str, ...] = ...
    potential: LatexStr = ...

kepler = PotentialSchema(
    """The Kepler Potential.

     Constants:
     - G: Gravitational constant

     Parameters:
     - M: total mass

     Variables:
     - r:  spherical radius
     """
    ("G",),
    ("M",),
    ("r",),
    "-\frac{G M}{r^2}"
)

Then libraries convert instances of the potential to some form, e.g. a dict, with all the schema fields
and the other library grabs that and turns it back into a potential.
For now users never see the schema, that's just the inter-library agreement for what are the constants, parameters, input variables, and mathematical expression which define each potential.

Maybe instead of just a dict we have an official object

class Potential(NamedTuple):
    kind: str
    constants: dict[str, Any]
    parameters: dict[str, Any]
    metadata: dict[str, Any] | None = None  # e.g. for supporting out-of-schema info like "positive_density" in gala's MN3Potential.

Then full-blown conversion would approximately look like

>>> from galactic_dynamics_interop import to_schema, from_schema, GalaxLibrary, GalpyLibrary
>>> from galax.potential import KeplerPotential as GalaxKeplerPotential
>>> from galpy.potential import KeplerPotential as GalpyKeplerPotentia

>>> galax_potential = GalaxKeplerPotential(...)
>>> schematized = to_schema(galax_potential)
>>> schematized
Potential(kind="kepler", "constants": {"G": 4.498e-6}, "parameters": {"M": 1e12}, metadata=None)

>>> galpy_potential = from_schema(GalpyLibrary, schematized)

>>> round_trip_galax_potential = from_schema(GalaxLibrary, to_schema(galpy_potential))

Of course libraries are free to obscure the existence of galactic_dynamics_interop, e.g.

>>> import galax.potential as gp
>>> pot = gp.KeplerPotential(...)
>>> galpy_pot = pot.as_interop("galpy")  # calls `from_schema(GalaxLibrary, to_schema(pot))

[jobovy]

That looks good! So the idea is that every library would implement their own version of to_schema(MyLibrary and from_schema(MyLibrary,...) here rather than in the library itself, right? Or we might want to add internal methods that basically do the to_schema at least and this code could just call those? I guess each library could make a different choice here.

To not just get bogged down in discussion, maybe you should implement a very basic version for just the KeplerPotential so we can play around with it, maybe in a branch, to see how everything comes together.

[eugvas]

worth trying on a trivial example, indeed.
I'm not sure I understand the difference between "constants" and "parameters" (and in particular, the need for "G" at all – my stance is that the gravitational framework should operate within the convention that G=1, and it's the job of the calling code to do any unit conversions, if necessary. Or at the very least, "G" can be absorbed into "GM" or "G rho0" etc. in all cases). And what are the "metadata"?
Speaking for Agama, the above schema would work well in simple cases, when the potential can be fully specified by parameters passed to the constructor. But that's not the most general way to create a potential; for example, a BasisSet potential can be created from parameters + input density or potential profile, or from parameters + input N-body snapshot, or read from a text file with coefficients. Technically, the input density profile or snapshot could be called a "parameter", but it needs to be specified in the way compatible with all other codes, while the coefficients file name can also be a "parameter", but its contents/format are certainly specific to each code.

[nstarman]

I'm not sure I understand the difference between "constants" and "parameters" (and in particular, the need for "G" at all – my stance is that the gravitational framework should operate within the convention that G=1, and it's the job of the calling code to do any unit conversions, if necessary

I think that's a good suggestion that eliminates the need to worry about units within this package.

And what are the "metadata"?

Additional information about a potential that aren't within the schema, but that another library might want to support. E.g. gala's MN3Potential has the flag "positive_density".

But that's not the most general way to create a potential;

Agreed. Let's certainly iterate on the schema!

---

I'm going to push a few refactor PRs before making an example branch for the schema. I need to 1. uncouple this from galax and 2. clean up the code structure to support having all the libraries be dependencies (following agreement on #5 (comment)).