astshim

A C++ and Python shim around a subset of Starlink AST a library for handling world coordinate systems in astronomy.

For detailed documentation of AST see http://starlink.eao.hawaii.edu/devdocs/sun211.htx/sun211.html.

The focus of astshim is on support for spatial mappings for use by LSST. Thus few of AST's functions that support time, spectra and tables have yet been wrapped. The python wrapper uses pybind11.

Differences between this and Starlink AST

• Most functions in AST are classes or methods in this shim.
• New class FrameDict which is a FrameSet that can reference frames by domain name.
• Mapping::applyForward and Mapping::applyInverse methods transform single points or lists of points. These replace AST's astTran<X> functions and no invert flag is supported. There are three versions of each method:
  • Accept a std::vector<double> and return a new std::vector<double>. This is the simplest version to use if you wish to transorm a single point.
  • Accept a 2-dimensional ndarray and return a newly allocated 2-dimensional ndarray.
  • Accept a 2-dimensional ndarray and fill a pre-allocated 2-dimensional ndarray.
• Mappings may not be inverted in place. Instead call Mapping::inverted to get an inverse mapping, and Mapping::isInverted to find out if a mapping is inverted.
• Compound mappings and frames have a few minor changes:
  • New classes SeriesMap and ParallelMap are subclasses of CmpMap and any compound mapping will always be an instance of one or the other of these subclasses. Thus Object.getClassName returns "SeriesMap" or "ParallelMap" for all compound maps.
  • Added Mapping::then as a convenient way to make a SeriesMap and Mapping::under as a convenient way to make a ParallelMap.
  • Added Frame.under as a convenient way to make a CmpFrame.
  • Instead of wrapping astDecompose:
    • Call CmpMap.operator[] to retrieve a mapping from a compound mapping.
    • Call CmpMap::getSeries to find out if a compound mapping is in series.
    • Call CmpFrame.operator[] to retrieve a Frame from a CmpFrame
• Exceptions are raised on errors, leaving AST in a normal (non-error) state:
  • Where practical, the wrapper checks arguments and throws std::invalid_argument before calling AST code.
  • After calling AST code the wrapper checks AST's error state and if invalid, the wrapper resets the AST status code and throws std::runtime_error.
• The AST functions set, set<X> and get<X> are hidden; instead each class has explicit accessors for its attributes, such as Object::getID. Mappings are mostly immutable, so they have getters, but no setters beyond a few generic setters from Object. SlaMap and TimeMap both violate immutability by having add methods; if this is a problem we can replace the add methods with constructor arguments. Frames are mutable, so all frame-specific attributes have setters as well as getters.
• Channels are constructed with a Stream; subclasses are available for files and strings, and it is easy (in C++, but not Python) to construct a Stream to use standard in and/or out.
• astshim manages memory using C++ smart pointers. Thus the following AST functions are not wrapped: astAnnul, astBegin, astClone, astDelete, astEnd, and astExport.
• Methods that output floating point data have AST__BAD replaced with nan.

Smaller differences (not a complete list):

• Methods such as Object::set and the attributes argument of constructors do not support printf formatting and extra arguments (and this may cause warnings in some compilers).
• Object::show prints its output to a provided output stream or returns a string, rather than printing to stdout.
• Get the class name using Object::getClassName() instead of getClass() because the latter sounds like a class, not a string, and Python doesn't allow class as a property name.
• FrameSet::addAxes and FrameSet::addFrame implement AST's astAddFrame(AST__ALLFRAMES, map, frame) function, because the AST function does two very different things.
• FitsChan does not offer methods for getting or setting complex integers (astGetFitsCI and astSetFitsCI), because that data type is not supported by standard C++.
• KeyMap has several differences:
  • Methods get<X> and put<X> work with both scalars and vectors (implementing astMapGet0<X>, astMapGet1<X>, astMapPut0<X> and astMapPut1<X>).
  • Methods append and replace are used to alter values in existing entries (implementing astMapPutElem...<X>).
• PolyMap has two separate constructors, one that takes foward and inverse coefficients, one that takes only forward coefficients. Neither provides an iterative inverse by default.

Supported AST versions

astshim works against Starlink AST 9.3.x and 9.5.0 or later. AST 9.4.x will not pass the packages tests.** This is because there is an edge case in one of the tests that causes a FITS solution to fail. This problem is fixed in 9.5.0. Where appropriate the tests do change behavior based on AST version since the role of this library is to assert that the result agrees with the C library even if that result depended on a bug.

Thread safety

astshim requires a build of AST that does not have POSIX thread support, which is how the starlink-ast conda-forge package is built (its recipe passes --without-pthreads to AST's configure). Against an AST built with thread support, astshim is unsafe to use from more than one Python thread, and the failure is a process-level crash rather than an exception.

In a thread-safe AST every object is locked to the thread that created it, and moving one to another thread requires the creating thread to unlock it and the receiving thread to lock it. Object::lock and Object::unlock wrap those calls, but nothing inside astshim uses them, so any object that is created on one thread and used on another is rejected by AST.

This is easy to hit without writing threaded code. Pickling an Object serializes it to text and rebuilds it with Channel::read, so the object that arrives is a genuinely separate AST object rather than a shared pointer - but it is created by whichever thread performs the unpickling. multiprocessing.Pool unpickles results on its internal result-handler thread while the main thread consumes them, so pool.map returning astshim objects is enough to trigger it. Against a thread-safe AST that reports "Invalid Object pointer ... owned by another thread" and then terminates the process, because the error escapes from a destructor.

What makes the current arrangement safe is CPython's global interpreter lock. The pybind11 bindings never release it, so two threads are never inside AST at the same time, and no astshim method leaves an AST operation open across calls: AST's global error status is checked, and the static buffers that functions such as astGetC return are copied into a std::string, before each wrapper function returns. AST's remaining global state is therefore only ever read and written while one thread holds the GIL.

Two changes would remove that protection and require a real locking discipline: running under a free-threaded CPython build (PEP 703), or releasing the GIL around AST calls in the bindings.

Missing Functionality

Many portions of AST have not yet been wrapped. Here are some highlights:

• Rebinning and resampling (astRebin<X>, astRebinSeq<X> and astResample<X>
• astDBSPecFrame
• astFluxFrame
• astRegion, astRemoveRegions and other region support
• astPlot and other plotting support
• astSpecFluxFrame
• astStcsChan
• astTable and other table support
• "astIsA<Class>". I don't think we need these (one can always use rtti or call Object::getClassName) but if they are needed then it probably makes sense to wrap these as isInstance(Object const & obj) static methods on all classes.
• and astTuneC.

The Python interface could present a more dict-like view of KeyMap and FitsChan, as pyast does.

To build

• Install a minimal LSST DM stack that includes sconsUtils, pybind11, ndarray and doxygen.
• Install starlink-ast as any LSST package.
• Install this package in the same way as any LSST package.

To build only the documentation:

• Install doxygen
• Execute this command: doxygen doc/doxygen.conf
• See file doc/html/index.html

Examples

See the unit tests in tests and the examples in examples

License

This product includes software developed by the LSST Project (http://www.lsst.org/).

This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

See https://www.lsstcorp.org/LegalNotices/ for the LSST License Statement and the GNU General Public License.