A thin C++ shim for Starlink AST

Data Management
#1

As per RFC-211 we are proposing to use Starlink AST package to represent WCS and other coordinate transforms in the LSST, thus replacing LSST's afw::image::Wcs, XYTransform and similar classes. Starlink AST is written in C, in an object-oriented style. This makes it somewhat difficult to extend (e.g. add new mappings) and somewhat clumsy to use from C++, so we hope to assist David Berry, AST's creator, to write a version in C++. Note that AST's API is somewhat outdated, but given how large AST is, its extensive documentation and its significant user base, we propose not making substantive changes.

Meanwhile, I propose to write a thin, generic C++ and Python shim for those parts of Starlink AST that we know we want to use now, a wrapper that represents an obvious translation of the existing API and closely matches PyAST, an existing Python wrapper for Starlink AST. This shim can easily be extended to additional parts of AST, as needed, and used as the API (possibly with some refinement) for a C++ rewrite of AST. Thus if the rewrite is a success we can easily adopt it. This document gives the details.

Basic Features of a C++ Shim for Starlink AST

The following features are definitely wanted in a C++ shim for AST. These basically already exist in PyAST:

  • Create C++ classes for fundamental AST classes such astObject, astMapping and astFrameSet.
  • Change free functions that act like methods, such as astTranP and astSimplify, into methods (member functions) of the appropriate fundamental classes.
  • Put the code in namespace ast, eliminate the ast prefix on function and class names, and use a lowercase letter to start function and method names. Thus astMapping (a class) becomes ast::Mapping in the shim, and astTranP (a function that acts like a method a.k.a. member function) becomes ast::Mapping::tranP.
  • Raise exceptions for errors.
  • Use std::vector and ndarray::Array where appropriate, instead of C style arrays.

Additional Changes Wanted Now

The following are desirable and fairly easy, so I feel they are worth doing now:

  • Simplify I/O with some convenience functions or methods
  • Make forward vs. inverse transforms more obvious than a flag provided to astTranP and similar functions, e.g. with separately named functions. This will make it easier to understand the intent of calling code.
  • For transforming collections of points, support data in C order, e.g. x0, y0, x1, y1, .... This is more natural for C/C++ and Python users, and allows LSST to naturally represent such data as ndarray::Array<afw::geom::Point2D>.
  • Minor cleanups, such as
    • Make it easy to chain more than two transforms. This can be done as an additional constructor for ast::CmpMap that accepts a vector of mappings.
    • Frame set method astAddFrame does two very different things, depending on the value of its first argument; I propose replacing it with two methods: addFrame and addAxes, where the former adds one frame and the latter adds axes to all frames in the frame set.

Features to Postpone

I would like to have the following changes, but feel they are out of scope for a shim and should be postponed for a possible rewrite of AST in C++:

  • Improve the API for specifying and getting attributes, to allow direct access and numbers instead of strings.
  • Simplify how inverse transforms are handled. I consider the current system too complicated, especially when trying to make a compound mapping with an inverse transform, and internally, when trying to simplify mappings. My preference is to eliminate the invert flag, and instead be able to obtain an inverse transform from any other transform. In many cases the inverse is trivial, for example the inverse of an affine transform is another affine transform with appropriate coefficients. In other cases we can provide an inverse mapping that contains the original transform. However, I doubt it is practical for the shim to do this.
  • Replace AST__BAD with IEEE nan. This matches modern user expectations and would simplify and speed up AST, because it contains many explicit tests for AST__BAD which could be eliminated.

Prototype

A prototype is available as astshim.

Details

The remainder of this document describes the shim in more detail. The I/O has not yet been fleshed out.

Wrap Now

This section describes those aspects of AST that we know we need now, and so will be part of the shim.

The following classes and their methods need wrapping

  • Object: base object for AST
    • Object(astObject * object)
    • clear(string attrib)
    • clone() -> shared_ptr<Object> # do we really need this?
    • fromString(string string) -> shared_ptr<Object> # this should probably be a constructor or class method on each class. That would be more convenient than trying to figure out what class was reconstructed.
    • getT(string attrib) ->, T where T is one of C, F, D, I, L
    • hasAttribute(string attrib) -> bool
    • same(Object obj) -> bool
    • set(string setting)
    • setT(string attrib, T value), T where T is one of C, F, D, I, L
    • test(string attrib) -> bool
  • Mapping(Object):
    • Mapping(astMapping * mapping)
    • decompose() -> struct containing map1, map2, series, invert1, invert2
    • invert()
    • linearApprox(vector<double> lbnd, vector<double> ubnd, double tol) -> vector<double>
    • mapBox(vector<double> lbnd_in, vector<double> ubnd_in, int forward, vector<double> &lbnd_out, vector<double> &ubnd_out, vector<double> *xl, vector<double> *xu)
    • quadApprox(vector<double> lbnd[2], vector<double> ubnd[2], int nx, int ny) -> struct containing fit and rms vectors of double
    • rate(vector<double> at, int ax1, int ax2)
    • removeRegions() -> Mapping
    • simplify()
    • tran(ndarray data) -> ndarray # data as an ndarray with dimensions #axes x #points
    • tranGrid(vector<int> lbnd, vector<int> ubnd, double tol, int maxpix, int forward) -> ndarray
    • report() -> string
    • property getters:
      • isSimple()
      • isLinear()
      • getNin()
      • getNout()
  • Wrap the following concrete mappings. If they have no additional methods (as in most cases) they can be wrapped as a simple factory function, which reduces boilerplate:
    • CmpMap has two constructors/factory functions, but no additional methods:
      • CmpMap(Mapping map1, Mapping map2, bool series, std::string options="")
      • CmpMap(std::vector<Mapping>, bool series, std::string options="")
    • PolyMap has one constructor and one additional method:
      • PolyMap(int nin, int nout, vector<double> coeff_f, vector<double> coeff_i, std::string options="")
      • polyTran(int forward, double acc, double maxacc, int maxorder, vector<int> lbnd, vector<int> ubnd) -> shared_ptr<PolyMap>
    • SlaMap has one trivial constructor and one addition method:
      • SlaMap(int flags=0, std::string options="")
      • slaAdd(string cvt, vector<double> args)
    • LutMap(vector<lut>, double start, double inc, std::string options="")
    • MathMap(int nin, int nout, vector<string> fwd, vector<string> inv, std::string options="")
    • MatrixMap has two constructors (AST offers a third with no matrix which results in a unit map; why bother?):
      • MatrixMap(ndarray::Array<double, 2, 2> matrix, std::string options="") full matrix
      • MatrixMap(ndarray::Array<double, 1, 1 matrix, std::string options="") diagonal elements as a 1-D array; nin=nout=length
    • NormMap(Frame frame, std::string options="")
    • PcdMap(double disco, vector<double> pcdcen, std::string options="")
    • PermMap(int nin, vector<int> inperm, int nout, vector<int> outperm, vector<double> constant, std::string options="")
    • RateMap(Mapping map, int ax1, int ax2, std::string options="")
    • ShiftMap(int ncoord, vector<double> shift, std::string options="")
    • SkyOffsetMap(Frame frame)
    • SphMap(std::string options="")
    • TranMap(Mapping map1, Mapping map2, std::string options="")
    • UnitMap(int ncoord, std::string options="")
    • UnitNormMap(int ncoord, std::string options="")
    • WcsMap(int ncoord, int type, double lonax, double latex, std::string options="")
    • WinMap(int ncoord, vector<double> in, vector<double> ina, vector<double> outa, vector<double> outb, std::string options="")
    • ZoomMap(int ncoord, double zoom, std::string options="")
  • Frame(Mapping)
    • Frame(astFrame * frame)
    • angle(vector<double> a, vector<double> b, vector<double> c) -> double
    • axAngle(vector<double> a, vector<double> b, int axis) -> double
    • axDistance(int axis, double v1, double v2) -> double
    • axOffset(int axis, double v1, double dist)
    • convert(Frame to, std::string domainlist) -> Frameset
    • distanc(vector<double> point1, vector<double> point2) -> double
    • matchAxes(Frame frm) -> vector<int>
    • norm(vector<double> &value)
    • offset(vector<double> point1, vector<double> point2, double offset, vector<double> &point_out)
    • offset2(vector<double> point1[2], double angle, double offset, vector<double> &point2[2]) -gt; double
    • permAxes(vector<int> perm)
  • CmpFrame(Frame)
  • SkyFrame(Frame)
  • FrameSet(Mapping), including methods:
    • FrameSet(astFrameSet * frameSet)
    • FrameSet.readFits(std::string path)
    • FrameSet(Frame frame)
    • FrameSet(Frame fromFrame, Mapping mapping, Frame toFrame)
    • addFrame(int iframe, Mapping map, Frame frame)
    • getMapping(ind1, ind2) -> Mapping
    • getFrame(ind) -> Frame
    • removeFrame(ind)
    • remapFrame(int iframe, Mapping map)

The following free functions need wrapping:

  • tune(string name, int value)
  • tuneC(string name, string value) -> string oldvalue
  • version() -> int

We need the following classes for I/O. Probably start by wrapping them directly, but I'd like syntactic sugar to simplify I/O.

  • Channel, including putChannelData, read, write
  • channelData
  • FitsChan, including getFits, purgeWcs, putCards, putFits, readFits, removeTables, setFits..., showFits, testFits, writeFits
  • XmlChan

Perhaps Wrap Later

This section describes aspects of AST that we may want someday. We will not wrap them now. If the C++ rewrite of AST goes through, then these will be updated as part of that.

The following classes will not be wrapped until needed:

  • Circle, including circlePars
  • DsbSpecFrame
  • Ellipse, including ellipsePars
  • FitsTable, including getTableHeader, putColumnData,
  • FluxFrame
  • GrismMap
  • IntraMap
  • KeyMap, including mapCopy, mapDefined, mapGet0<X>, mapGet1<X>, mapGetElem<X>, mapHasKey, mapKey, mapLenC, mapLength, mapPut0<X>, mapPut1<X>, mapPutElem<X>, mapPutU, mapRemove, mapRename, mapSize, mapType,
  • NullRegion
  • Plot, including curve, eBuf, genCurve, grfPop, grfPush, grfSet, grid, gridLine, mark, polyCurve, regionOutline, text
  • Plot3D
  • PointList
  • Polygon, including downsize
  • Prism
  • Region, including getRegionBounds, getRegionFrame, getRegionFrameSet, getRegionMesh, getRegionPoints, getUnc, mapRegion, mask<X>, negate, overlap, setUnc, showMesh
  • SelectorMap
  • SpecFluxFrame
  • SpecFrame, including getRefPos, setRefPos
  • SpecFluxFrame
  • SpecMap, including specAdd
  • Stc, including getStcCoord, getStcNCoord, getStcRegion
  • StcCatalogEntryLocation
  • StcObsDataLocation
  • StcResourceProfile
  • StcSearchLocation
  • StcsChan
  • SwitchMap
  • Table, including columnName, columnNull, columnShape, columnSize, getColumnData, getGrfContext, getTables, hasColumn, hasParameter, parameterName, purgeRows, removeColumn, removeParameter, removeRow
  • TimeFrame, including currentTime
  • TimeMap, including timeAdd

The following methods will not be wrapped until needed:

  • Channel:
    • warnings() -> KeyMap # can't wrap unless we wrap KeyMap or return a different type
  • FitsChan:
    • getTables, putTable, putTableHeader, putTables, retainFits, tableSource
  • Frame:
    • circle
    • findFrame(prototypeFrame) # important, but will modify the FrameSet; applies to Frame and FrameSet
    • format(int axis, double value) -> std::string
    • getActiveUnit() -> int
    • intersect
    • interval
    • resolve
    • setActiveUnit(int value)
    • pickAxes(vector<int> axes) -> shared_ptr<Mapping>
    • unformat(int axis, std::string name) -> double # perhaps take char * instead of std::string
  • FrameSet
    • addVariant(Mapping map, string name)
    • mirrorVariants(int iframe)
  • Mapping:
    • mapSplit
    • resample
    • rebin
    • rebinSeq
  • Object:
    • lock
    • thread
    • unlock

The following free functions will not be wrapped until needed:

  • convex
  • delete
  • escapes
  • exempt
  • export
  • import
  • intraReg
  • outline # this should probably be a variant constructor for Polygon
  • stripEscapes

The following function prototypes will not be wrapped:

  • uinterp
  • ukern1

Never Wrap

The following code will not be wrapped because the wrapper code eliminates the need for users to call it. The wrapper code will call these functions internally.

Our wrapper code will perform all status checking and raise exceptions as appropriate, so we will not wrap the following:

  • clearStatus
  • setStatus
  • status
  • watch
  • ok

Our wrapper will automatically set up AST and tear it down, so we will not wrap the following:

  • begin
  • end
DM Highlights 2016-08-03 to 2016-08-09
#4

I can’t quite see what API you are proposing to support from this description. Can you document the complete API that we will be using to replace the current afw::image::Wcs class?

I am not trying to reopen our decision to implement this using AST under the covers.

#7

Excellent question. That is described in RFD-211. I was trying to post two things at once that refer to each other, and this one ended up visible first. I’ll be adding a link to the text above, now that I have one.

#8

What is the proposal for dealing with the reinterpret_cast problem?

#9

As specified in the proposal, we will have a factory function for each concrete mapping class (e.g. astZoomMap, etc.), unless the mapping has additional methods, in which case it will have a class instead.

If the pointers had been compatible then it might have been acceptable to not wrap the mapping classes and instead use this shortcut:

auto zoomMap = ast::Mapping(astZoomMap(...));

However, using such a shortcut is quite unpleasant for ast mappings that take arrays (in C++ we can use ndarray or std::vector and know how big the array is, but in C we need to pass a pointer and shape information separately), so I’m happier with the proposal in any case.

If you need to construct some AST object that does not have a shim, you can do the following:

auto foo = ast::Object(std::reinterpret_cast<astObject *>astFoo(...));

where Object can be replaced by Mapping, Frame or whatever other main class seems appropriate. But that is clearly not the way we should expect users to normally use the ast shim.