Add class and methods to get centroid residuals

[jeanconn]

Add class and methods to get centroid residuals from variety of possible inputs.

[jeanconn]

This is a first draft to support #18. Could use refactoring, but probably good to make a PR now to see if the API is somewhat working out as expected.

[jeanconn]

Regarding "Add this convenience method as a class method on CentroidResiduals. ", assuming you mean something like get_obs_slot_residuals, I'm not sure if that just complicates CentroidResiduals, which I now have as taking a start/stop on initialization. What is your thought now about how the class method would work?

[taldcroft]

@classmethod
def for_slot(cls, obsid=None, start=None, stop=None,
             slot=None, att_source='ground', centroid_source='ground'):
    if obsid is not None:
        if start is not None or stop is not None:
            raise ValueError('cannot specify both obsid and start / stop')
        ds = events.dwells.filter(obsid=obsid)
        start = ds[0].start
        stop = ds[len(ds) - 1].stop

    if start is None or stop is None:
        raise ValueError('must specify obsid or start / stop')

    cr = cls(start, stop)
    cr.set_atts(att_source)
    cr.set_centroids(centroid_source, slot)
    cr.set_star(obsid=obsid, slot=slot, date=start)
    cr.calc_residuals()  # instead of get_residuals
    return cr

cr = CentroidResiduals.for_slot(obsid=19387, slot=6)

I know there was a good reason why this for_slot class method could not be simply the class __init__ in terms of use cases / max flexibility. But I don't quite remember the detailed reason...

[taldcroft]

I think the full interface (using the raw class methods instead of for_slot) could be useful for evaluating simulated or non-flight data.

[jeanconn]

Regarding time offsets for the centroids, I used the centroid residual code in loops over a set of bright stars with small V&V residuals, with a range of time offsets and minimized the sum of the squared residuals. Here are some plots of those "solved" values of the time offsets (dt). On the X axis there are just individual obsids. 'dt' fit with obc centroids have star markers; ground centroids are circles. 'dt' fit with ground aspect solution have red markers; obc attitudes are blue.

I processed several ERs to do something similar with ER 8x8 data:

[taldcroft]

That's interesting. It seems like the OBC attitude solution wanders in phase. We might want to come back to this later because I have a nagging feeling there is some detail we are not understanding. For now, just take the median in each case and call it a day.

Looks like the class now needs to know about the fetch data source fetch_source = 'cxc' (default), 'maude' which can be used in the internal fetching when needed. Probably the ER/OR distinction can be automatic based on the obsid. (Where obsid is not provided one can always try to be clever and find obsid by querying MAUDE for COBSRQID over a small time range centered on the start/stop interval).

[jeanconn]

Regarding time offsets, my first action item out of our meeting was to quickly check to see if the obc and ground aspect solutions are offset in time. I just looked at dwells in 2014, and fitting an offset for the obc solution relative to the ground solution using 1ks of data in each available science dwell gives this: Median offset of -0.32 seconds. Do you want to do anything with this? It explains some of the offsets of the centers of the blue symbols compared to respective red symbols, but there's again, a huge amount of spread and I don't understand which is the "right" answer.

[taldcroft]

Can you send by email the obsid for the worst-case outlier (below -5 secs) and a pointer to the code you used to make that plot?

[jeanconn]

I was totally brute-forcing those attitude fits, so really no idea if they are reasonable. With no use of sherpa/etc to give me errors, I was just hoping that if I ran enough observations that I'd see if there was a systematic that made some sense.

import chandra_aca.centroid_resid
from kadi import events
import numpy as np
from Ska.Numpy import interpolate
from mica.quaternion import normalize
from Quaternion import Quat

def rss_atts(atts1, atts1_times, atts2, atts2_times):
    # I suppose this doesn't take a root, so it is just a ss_atts
    # interpolate att2 atts at att1 times
    x = interpolate(atts2[:, 0], atts2_times, atts1_times)
    y = interpolate(atts2[:, 1], atts2_times, atts1_times)
    z = interpolate(atts2[:, 2], atts2_times, atts1_times)
    w = interpolate(atts2[:, 3], atts2_times, atts1_times)
    atts = np.vstack([x, y, z, w]).transpose()
    atts = normalize(atts)
    errs = []
    for att1, att2 in zip(atts1, atts):
        dq = Quat(att1).dq(Quat(att2))
        errs.append(dq.pitch**2 + dq.yaw**2)
    return np.sum(errs)

ds = events.dwells.filter('2014:001', '2015:001', dur__gt=3000)
best = []
for dwell in ds:
    print(dwell.get_obsid())
    try:
        cr_obc = chandra_aca.centroid_resid.CentroidResiduals(dwell.tstart + 1000, dwell.tstart + 2000)
        cr_obc.set_atts('obc')
        cr_ground = chandra_aca.centroid_resid.CentroidResiduals(dwell.tstart + 1000, dwell.tstart + 2000)
        cr_ground.set_atts('ground')
    except:
        continue
    dts = np.arange(-5, 3, .1)
    dt_errsum = []
    for dt in dts:
        dt_errsum.append(rss_atts(cr_obc.atts, cr_obc.att_times + dt,
                                  cr_ground.atts, cr_ground.att_times))
    p = np.polyfit(dts, np.array(dt_errsum), 2)
    best.append(-1 * p[1] / (2.0 * p[0]))
    print(best[-1], dwell.get_obsid())

[taldcroft]

See: http://nbviewer.jupyter.org/url/asc.harvard.edu/mta/ASPECT/ipynb/calc_obc_ground_att_time_offset.ipynb

There are a few outliers but overall the distribution is fairly tight around 0.08 sec. Takeaways:

• Use the cleanest possible data from early 2008
• Use at least 3 dither cycles
• Never interpolate quaternion values directly
• I took out any mean offset in pitch/yaw as being "uninteresting"

[jeanconn]

Regarding "Never interpolate quaternion values directly", is that really the source of the problem (I will review)? I thought for relatively small changes, this does work out OK even without slerp (and I thought that we are doing some quaternion by-element interpolation in the ground aspect pipeline). If you just interpolate pitches and yaws, doesn't that throw away roll noise?

I think when I was just looking at OBC vs ground solutions recently I went with "nearest" ground value to avoid the by-element interpolation, but something made me think it was OK for this....

[taldcroft]

Looks like the COBSRQID ..

That seems like a bug. Ask Dave about it.

[jeanconn]

New offsets and plots of the source data are in the notebooks referenced from the code: http://nbviewer.jupyter.org/url/cxc.harvard.edu/mta/ASPECT/ipynb/centroid_time_offsets/OR.ipynb and http://nbviewer.jupyter.org/url/cxc.harvard.edu/mta/ASPECT/ipynb/centroid_time_offsets/ER.ipynb

[jeanconn]

I probably should have bumped the version after merging. Now I think I need to rebase or merge.