Run the discrete-Fourier-transform photometry algorithm.

    See the module-level documentation and the output of ``casatask dftphotom
    --help`` for help. All of the algorithm configuration is specified in the
    *cfg* argument, which is an instance of :class:`Config`.

def dftphotom(cfg):
    """Run the discrete-Fourier-transform photometry algorithm.

    See the module-level documentation and the output of ``casatask dftphotom
    --help`` for help. All of the algorithm configuration is specified in the
    *cfg* argument, which is an instance of :class:`Config`.

    """
    tb = util.tools.table()
    ms = util.tools.ms()
    me = util.tools.measures()

    # Read stuff in. Even if the weight values don't have their
    # absolute scale set correctly, we can still use them to set the
    # relative weighting of the data points.
    #
    # datacol is (ncorr, nchan, nchunk)
    # flag is (ncorr, nchan, nchunk)
    # weight is (ncorr, nchunk)
    # uvw is (3, nchunk)
    # time is (nchunk)
    # axis_info.corr_axis is (ncorr)
    # axis_info.freq_axis.chan_freq is (nchan, 1) [for now?]
    #
    # Note that we apply msselect() again when reading the data because
    # selectinit() is broken, but the invocation here is good because it
    # affects the results from ms.range() and friends.

    if ':' in (cfg.spw or ''):
        warn('it looks like you are attempting to select channels within one or more spws')
        warn('this is NOT IMPLEMENTED; I will average over the whole spw instead')

    ms.open(b(cfg.vis))
    totrows = ms.nrow()
    ms_sels = dict((n, cfg.get(n)) for n in util.msselect_keys
                   if cfg.get(n) is not None)
    ms.msselect(b(ms_sels))

    rangeinfo = ms.range(b'data_desc_id field_id'.split())
    ddids = rangeinfo['data_desc_id']
    fields = rangeinfo['field_id']
    colnames = [cfg.datacol] + 'flag weight time axis_info'.split()
    rephase = (cfg.rephase is not None)

    if fields.size != 1:
        # I feel comfortable making this a fatal error, even if we're
        # not rephasing.
        die('selected data should contain precisely one field; got %d', fields.size)

    if rephase:
        fieldid = fields[0]
        tb.open(b(os.path.join(cfg.vis, 'FIELD')))
        phdirinfo = tb.getcell(b'PHASE_DIR', fieldid)
        tb.close()

        if phdirinfo.shape[1] != 1:
            die('trying to rephase but target field (#%d) has a '
                'time-variable phase center, which I can\'t handle', fieldid)
        ra0, dec0 = phdirinfo[:,0] # in radians.

        # based on intflib/pwflux.py, which was copied from
        # hex/hex-lib-calcgainerr:

        dra = cfg.rephase[0] - ra0
        dec = cfg.rephase[1]
        l = np.sin(dra) * np.cos(dec)
        m = np.sin(dec) * np.cos(dec0) - np.cos(dra) * np.cos(dec) * np.sin(dec0)
        n = np.sin(dec) * np.sin(dec0) + np.cos(dra) * np.cos(dec) * np.cos(dec0)
        n -= 1 # makes the work below easier
        lmn = np.asarray([l, m, n])
        colnames.append('uvw')

        # Also need this although 99% of the time `ddid` and `spwid` are the same
        tb.open(b(os.path.join(cfg.vis, 'DATA_DESCRIPTION')))
        ddspws = np.asarray(tb.getcol(b'SPECTRAL_WINDOW_ID'))
        tb.close()

    tbins = {}
    colnames = b(colnames)

    for ddindex, ddid in enumerate(ddids):
        # Starting in CASA 4.6, selectinit(ddid) stopped actually filtering
        # your data to match the specified DDID! What garbage. Work around
        # with our own filtering.
        ms_sels['taql'] = 'DATA_DESC_ID == %d' % ddid
        ms.msselect(b(ms_sels))

        ms.selectinit(ddid)
        if cfg.polarization is not None:
            ms.selectpolarization(b(cfg.polarization.split(',')))
        ms.iterinit(maxrows=4096)
        ms.iterorigin()

        while True:
            cols = ms.getdata(items=colnames)

            if rephase:
                # With appropriate spw/DDID selection, `freqs` has shape
                # (nchan, 1). Convert to m^-1 so we can multiply against UVW
                # directly.
                freqs = cols['axis_info']['freq_axis']['chan_freq']
                assert freqs.shape[1] == 1, 'internal inconsistency, chan_freq??'
                freqs = freqs[:,0] * util.INVERSE_C_MS

            for i in range(cols['time'].size): # all records
                time = cols['time'][i]
                # get out of UTC as fast as we can! For some reason
                # giving 'unit=s' below doesn't do what one might hope it would.
                # CASA can convert to a variety of timescales; TAI is probably
                # the safest conversion in terms of being helpful while remaining
                # close to the fundamental data, but TT is possible and should
                # be perfectly precise for standard applications.
                mq = me.epoch(b'utc', b({'value': time / 86400., 'unit': 'd'}))
                mjdtt = me.measure(b(mq), b'tt')['m0']['value']

                tdata = tbins.get(mjdtt, None)
                if tdata is None:
                    tdata = tbins[mjdtt] = [0., 0., 0., 0., 0]

                if rephase:
                    uvw = cols['uvw'][:,i]
                    ph = np.exp((0-2j) * np.pi * np.dot(lmn, uvw) * freqs)

                for j in range(cols['flag'].shape[0]): # all polns
                    # We just average together all polarizations right now!
                    # (Not actively, but passively by just iterating over them.)
                    data = cols[cfg.datacol][j,:,i]
                    flags = cols['flag'][j,:,i]

                    # XXXXX casacore is currently (ca. 2012) broken and
                    # returns the raw weights from the dataset rather than
                    # applying the polarization selection. Fortunately all of
                    # our weights are the same, and you can never fetch more
                    # pol types than the dataset has, so this bit works
                    # despite the bug.

                    w = np.where(~flags)[0]
                    if not w.size:
                        continue # all flagged

                    if rephase:
                        data *= ph

                    d = data[w].mean()
                    # account for flagged parts. 90% sure this is the
                    # right thing to do:
                    wt = cols['weight'][j,i] * w.size / data.size
                    wd = wt * d
                    # note a little bit of a hack here to encode real^2 and
                    # imag^2 separately:
                    wd2 = wt * (d.real**2 + (1j) * d.imag**2)

                    tdata[0] += wd
                    tdata[1] += wd2
                    tdata[2] += wt
                    tdata[3] += wt**2
                    tdata[4] += 1

            if not ms.iternext():
                break

        ms.reset() # reset selection filter so we can get next DDID

    ms.close()

    # Could gain some efficiency by using a better data structure than a dict().
    smjd = sorted(six.iterkeys(tbins))
    cfg.format.header(cfg)

    for mjd in smjd:
        wd, wd2, wt, wt2, n = tbins[mjd]
        if n < 3: # not enough data for meaningful statistics
            continue

        dtmin = 1440 * (mjd - smjd[0])
        r_sc = wd.real / wt * cfg.datascale
        i_sc = wd.imag / wt * cfg.datascale
        r2_sc = wd2.real / wt * cfg.datascale**2
        i2_sc = wd2.imag / wt * cfg.datascale**2

        if cfg.believeweights:
            ru_sc = wt**-0.5 * cfg.datascale
            iu_sc = wt**-0.5 * cfg.datascale
        else:
            rv_sc = r2_sc - r_sc**2 # variance among real/imag msmts
            iv_sc = i2_sc - i_sc**2
            ru_sc = np.sqrt(rv_sc * wt2) / wt # uncert in mean real/img values
            iu_sc = np.sqrt(iv_sc * wt2) / wt

        mag = np.sqrt(r_sc**2 + i_sc**2)
        umag = np.sqrt(r_sc**2 * ru_sc**2 + i_sc**2 * iu_sc**2) / mag
        cfg.format.row(cfg, mjd, dtmin, r_sc, ru_sc, i_sc, iu_sc, mag, umag, n)