Remplace PRESTO FFTFIT with native

src/pint/profile/__init__.py

@@ -1,26 +1,57 @@
 """Tools for working with pulse profiles.
 
-The key tool here is `fftfit`, which allows one to find the phase shift that
-optimally aligns a template with a profile, but there are also tools here for
-doing those shifts and generating useful profiles.
+The key tool here is FFTFIT (:func:`pint.profile.fftfit_full`), which allows
+one to find the phase shift that optimally aligns a template with a profile,
+but there are also tools here for doing those shifts and generating useful
+profiles.
 """
+
 import numpy as np
 import scipy.stats
+from numpy.fft import rfft, irfft
+
 import pint.profile.fftfit_aarchiba
 import pint.profile.fftfit_nustar
 import pint.profile.fftfit_presto
 
-
 __all__ = [
+    "fftfit_full",
+    "fftfit_basic",
+    "FFTFITResult",
+    "fftfit_cprof",
+    "fftfit_classic",
     "wrap",
     "vonmises_profile",
     "upsample",
     "shift",
-    "fftfit_full",
-    "fftfit_basic",
 ]
 
 
+class FFTFITResult:
+    """Summary of the results of an FFTFit operation.
+
+    Not all of these attributes may be present in every object; which are
+    returned depends on the algorithm used and the options it is passed.
+
+    If these quantities are available, then
+    r.scale*shift(template, r.shift) + r.offset
+    should be as close as possible to the profile used in the fitting.
+
+    Attributes
+    ----------
+    shift : float
+        The shift required to make the template match the profile. Between 0 and 1.
+    scale : float
+        The amount the template must be scaled by to match the profile.
+    offset : float
+        The amount to add to the scaled template to match the profile.
+    uncertainty : float
+        The estimated one-sigma uncertainty in the shift attribute.
+    """
+
+    pass
+
+
 def wrap(a):
     """Wrap a floating-point number or array to the range -0.5 to 0.5."""
     return (a + 0.5) % 1 - 0.5
@@ -112,10 +143,27 @@ def irfft_value(c, phase, n=None):
 def fftfit_full(template, profile, code="aarchiba"):
     """Match template to profile and return match properties.
 
-    The returned object has a `.shift` attribute indicating the optimal shift,
+    The returned object, a :class:`pint.profile.FFTFITResult`, has a
+    `.shift` attribute indicating the optimal shift,
     a `.uncertainty` attribute containting an estimate of the uncertainty, and
     possibly certain other attributes depending on which version of the code is
     run.
+
+    The ``.shift`` attribute is computed so that ``shift(template, r.shift)`` is
+    as closely aligned with ``profile`` as possible.
+
+    Parameters
+    ----------
+    template : array
+        The template representing the ideal pulse profile.
+    profile : array
+        The observed profile the template should be aligned with.
+    code : "aarchiba", "nustar", "presto"
+        Which underlying algorithm and code should be used to carry out the
+        operation. Generally the "aarchiba" code base is the best tested
+        under idealized circumstances, and "presto" is the compiled FORTRAN
+        code FFTFIT that has been in use for a very long time (but is not
+        available unless PINT has access to the compiled code base of PRESTO).
     """
     if code == "aarchiba":
         return pint.profile.fftfit_aarchiba.fftfit_full(template, profile)
@@ -130,8 +178,104 @@ def fftfit_full(template, profile, code="aarchiba"):
 
 
 def fftfit_basic(template, profile, code="aarchiba"):
-    """Return the optimal phase shift to match template to profile."""
+    """Return the optimal phase shift to match template to profile.
+
+    This calls :func:`pint.profile.fftfit_full` and extracts the ``.shift`` attribute.
+
+    Parameters
+    ----------
+    template : array
+        The template representing the ideal pulse profile.
+    profile : array
+        The observed profile the template should be aligned with.
+    code : "aarchiba", "nustar", "presto"
+        Which code to use. See :func:`pint.profile.fftfit_full` for details.
+    """
     if code == "aarchiba":
         return pint.profile.fftfit_aarchiba.fftfit_basic(template, profile)
     else:
         return fftfit_full(template, profile, code=code).shift
+
+
+def fftfit_cprof(template):
+    """Transform a template for use with fftfit_classic.
+
+    Emulate the version of fftfit.cprof in PRESTO.
+    Returns results suitable for :func:`pint.profile.fftfit_classic`.
+
+    Parameter
+    ---------
+    template : array
+
+    Returns
+    -------
+    c : float
+        The constant term? This may not agree with PRESTO.
+    amp : array
+        Real values representing the Fourier amplitudes of the template, not
+        including the constant term.
+    pha : array
+        Real values indicating the angles of the Fourier coefficients of the
+        template, not including the constant term; the sign is the negative
+        of that returned by ``np.fft.rfft``.
+    """
+    tc = rfft(template)
+    return tc[0], np.abs(tc)[1:], -np.angle(tc)[1:]
+
+
+def fftfit_classic(profile, template_amplitudes, template_angles, code="aarchiba"):
+    """Emulate the version of fftfit in PRESTO.
+
+    This has a different calling and return convention.
+    The template can be transformed appropriately with
+    :func:`pint.profile.fftfit_cprof`.
+
+    Parameters
+    ----------
+    profile : array
+        The observed profile.
+    template_amplitudes : array
+        Real values representing the Fourier amplitudes of the template, not
+        including the constant term.
+    template_angles : array
+        Real values indicating the angles of the Fourier coefficients of the
+        template, not including the constant term; the sign is the negative
+        of that returned by ``np.fft.rfft``.
+
+    Returns
+    -------
+    shift : float
+        The shift, in bins, plus one.
+    eshift : float
+        The uncertainty in the shift, in bins.
+    snr : float
+        Some kind of signal-to-noise ratio; not implemented.
+    esnr : float
+        Uncertainty in the above; not implemented.
+    b : float
+        Unknown; not implemented.
+    errb : float
+        Uncertainty in the above; not implemented.
+    ngood
+        Unknown, maybe the number of harmonics used? Not implemented.
+    """
+    if code == "presto":
+        import presto.fftfit
+
+        return presto.fftfit.fftfit(profile, template_amplitudes, template_angles)
+    if len(profile) % 2:
+        raise ValueError("fftfit_classic only works on even-length profiles")
+    template_f = np.zeros(len(template_amplitudes) + 1, dtype=complex)
+    template_f[1:] = template_amplitudes * np.exp(-1.0j * template_angles)
+    template = irfft(template_f)
+    r = fftfit_full(template, profile, code=code)
+
+    shift = (r.shift % 1) * len(profile) + 1
+    eshift = r.uncertainty * len(profile)
+    snr = np.nan
+    esnr = np.nan
+    b = np.nan
+    errb = np.nan
+    ngood = np.nan
+
+    return shift, eshift, snr, esnr, b, errb, ngood

src/pint/profile/fftfit_aarchiba.py

@@ -13,10 +13,6 @@
 import pint.profile
 
 
-class FFTFITResult:
-    pass
-
-
 def fftfit_full(
     template, profile, compute_scale=True, compute_uncertainty=True, std=None
 ):
@@ -54,8 +50,8 @@ def gof(x):
     if not res.success:
         raise ValueError("FFTFIT failed: %s" % res.message)
     # assert gof(res.x) <= gof(x)
-    r = FFTFITResult()
-    r.shift = res.x % 1
+    r = pint.profile.FFTFITResult()
+    r.shift = pint.profile.wrap(res.x)
 
     if compute_scale or compute_uncertainty:
         # shifted template corefficients

src/pint/profile/fftfit_nustar.py

@@ -1,8 +1,11 @@
 # by mateobachetti
 # from https://github.com/NuSTAR/nustar-clock-utils/blob/master/nuclockutils/diagnostics/fftfit.py
 from collections import namedtuple
+
 import numpy as np
-from scipy.optimize import minimize, brentq
+from scipy.optimize import brentq, minimize
+
+import pint.profile
 
 
 def find_delay_with_ccf(amp, pha):
@@ -195,13 +198,9 @@ def fftfit_basic(template, profile):
     return shift
 
 
-class FullResult:
-    pass
-
-
 def fftfit_full(template, profile):
     r = fftfit(profile, template)
-    ro = FullResult()
+    ro = pint.profile.FFTFITResult()
     ro.shift = r.mean_phase
     ro.uncertainty = r.std_phase
     return ro

src/pint/profile/fftfit_presto.py

@@ -1,5 +1,6 @@
 import numpy as np
 from numpy.fft import rfft
+
 import pint.profile
 
 try:
@@ -8,10 +9,6 @@
     presto = None
 
 
-class FFTFITResult:
-    pass
-
-
 def fftfit_full(template, profile):
     if len(template) != len(profile):
         raise ValueError(
@@ -23,11 +20,11 @@ def fftfit_full(template, profile):
         raise ValueError(
             "template has length {} which is too long".format(len(template))
         )
-    tc = rfft(template)
+    _, amp, pha = pint.profile.fftfit_cprof(template)
     shift, eshift, snr, esnr, b, errb, ngood = presto.fftfit.fftfit(
-        profile, np.abs(tc)[1:], -np.angle(tc)[1:]
+        profile,
     )
-    r = FFTFITResult()
+    r = pint.profile.FFTFITResult()
     # Need to add 1 to the shift for some reason
     r.shift = pint.profile.wrap((shift + 1) / len(template))
     r.uncertainty = eshift / len(template)

src/pint/scripts/event_optimize.py

@@ -27,6 +27,7 @@
     UniformUnboundedRV,
 )
 from pint.observatory.satellite_obs import get_satellite_observatory
+from pint.profile import fftfit_cprof, fftfit_classic
 
 __all__ = ["read_gaussfitfile", "marginalize_over_phase", "main"]
 # log.setLevel('DEBUG')
@@ -141,13 +142,11 @@ def measure_phase(profile, template, rotate_prof=True):
             (returned as a tuple).  These are defined as in Taylor's
             talk at the Royal Society.
     """
-    import fftfit
-
-    c, amp, pha = fftfit.cprof(template)
+    c, amp, pha = fftfit_cprof(template)
     pha1 = pha[0]
     if rotate_prof:
         pha = np.fmod(pha - np.arange(1, len(pha) + 1) * pha1, 2.0 * np.pi)
-    shift, eshift, snr, esnr, b, errb, ngood = fftfit.fftfit(profile, amp, pha)
+    shift, eshift, snr, esnr, b, errb, ngood = fftfit_classic(profile, amp, pha)
     return shift, eshift, snr, esnr, b, errb, ngood
 
 
@@ -238,8 +237,7 @@ def marginalize_over_phase(
 
 
 def get_fit_keyvals(model, phs=0.0, phserr=0.1):
-    """Read the model to determine fitted keys and their values and errors from the par file
-    """
+    """Read the model to determine fitted keys and their values and errors from the par file"""
     fitkeys = [p for p in model.params if not getattr(model, p).frozen]
     fitvals = []
     fiterrs = []

tests/test_event_optimize.py

@@ -1,6 +1,4 @@
 #!/usr/bin/env python
-# This test is DISABLED because event_optimize requires PRESTO to be installed
-# to get the fftfit module.  It can be run manually by people who have PRESTO
 from __future__ import division, print_function
 
 import os