Repository: StingraySoftware/stingray
Pull request #214 (2017-07-19 21:25:24)

Functions to fit the epoch folding periodograms
These functions can be used to fit the periodograms obtained by epoch folding, using a Sinc or a Gaussian function.

TODO:

- [x] Gaussian fit for non-perfect peaks

- [x] Physically meaningful constraints on the parameters


Patch:
diff --git a/stingray/pulse/modeling.py b/stingray/pulse/modeling.py
new file mode 100644
index 00000000..c58e83fa
--- /dev/null
+++ b/stingray/pulse/modeling.py
@@ -0,0 +1,192 @@
+import numpy as np
+from astropy.modeling import models, fitting
+
+
+__all__ = ["sinc_square_model", "sinc_square_deriv", "fit_sinc",
+           "fit_gaussian", "SincSquareModel"]
+
+
+def sinc(x):
+    """
+    Calculate a sinc function.
+
+    sinc(x)=sin(x)/x
+
+    Parameters
+    ----------
+    x : array-like
+
+    Returns
+    -------
+    values : array-like
+    """
+    values = np.sinc(x/np.pi)
+    return values
+
+
+def sinc_square_model(x, amplitude=1., mean=0., width=1.):
+    """
+    Calculate a sinc-squared function.
+
+    (sin(x)/x)**2
+
+    Parameters
+    ----------
+    x: array-like
+
+    Other Parameters
+    ----------
+    amplitude : float
+        the value for x=mean
+    mean : float
+        mean of the sinc function
+    width : float
+        width of the sinc function
+
+    Returns
+    -------
+    sqvalues : array-like
+         Return square of sinc function
+
+    Examples
+    --------
+    >>> sinc_square_model(0, amplitude=2.)
+    2.0
+    """
+    sqvalues = amplitude * sinc((x-mean)/width) ** 2
+    return sqvalues
+
+
+def sinc_square_deriv(x, amplitude=1., mean=0., width=1.):
+    """
+    Calculate partial derivatives of sinc-squared.
+
+    Parameters
+    ----------
+    x: array-like
+
+    Other Parameters
+    ----------
+    amplitude : float
+        the value for x=mean
+    mean : float
+        mean of the sinc function
+    width : float
+        width of the sinc function
+
+    Returns
+    -------
+    d_amplitude : array-like
+         partial derivative of sinc-squared function
+         with respect to the amplitude
+    d_mean : array-like
+         partial derivative of sinc-squared function
+         with respect to the mean
+    d_width : array-like
+         partial derivative of sinc-squared function
+         with respect to the width
+
+    Examples
+    --------
+    >>> np.all(sinc_square_deriv(0, amplitude=2.) == [1., 0., 0.])
+    True
+    """
+    x_is_zero = x == mean
+
+    d_x = 2 * amplitude * sinc((x-mean)/width) * (x * np.cos((x-mean)/width) \
+        - np.sin((x-mean)/width)) / ((x-mean)/width)**2
+    d_x = np.asarray(d_x)
+    d_amplitude = sinc((x-mean)/width)**2
+    d_x[x_is_zero] = 0
+
+    d_mean = d_x*(-1/width)
+    d_width = d_x*(-(x-mean)/(width)**2)
+
+    return [d_amplitude, d_mean, d_width]
+
+
+SincSquareModel = models.custom_model(sinc_square_model,
+                                      fit_deriv=sinc_square_deriv)
+
+
+def fit_sinc(x, y, amp=1.5, mean=0., width=1., tied={}, fixed={}, bounds={},
+             obs_length=None):
+    """
+    Fit a sinc function to x,y values.
+
+    Parameters
+    ----------
+    x : array-like
+    y : array-like
+
+    Other Parameters
+    ----------------
+    amp : float
+        The initial value for the amplitude
+    mean : float
+        The initial value for the mean of the sinc
+    obs_length : float
+        The length of the observation. Default None. If it's defined, it
+        fixes width to 1/(pi*obs_length), as expected from epoch folding
+        periodograms
+    width : float
+        The initial value for the width of the sinc. Only valid if
+        obs_length is 0
+    tied : dict
+    fixed : dict
+    bounds : dict
+        Parameters to be passed to the astropy models
+        http://docs.astropy.org/en/stable/api/astropy.modeling.functional_models.Gaussian1D.html
+
+    Returns
+    -------
+    sincfit : function
+        The best-fit function, accepting x as input
+        and returning the best-fit model as output
+    """
+    if obs_length is not None:
+        width = 1 / (np.pi * obs_length)
+        fixed["width"] = True
+
+    sinc_in = SincSquareModel(amplitude=amp, mean=mean,width=width, tied=tied,
+                              fixed=fixed, bounds=bounds)
+    fit_s = fitting.LevMarLSQFitter()
+    sincfit = fit_s(sinc_in, x, y)
+    return sincfit
+
+
+def fit_gaussian(x, y, amplitude=1.5,mean=0.,stddev=2., tied={}, fixed={},
+                 bounds={}):
+    """
+    Fit a gaussian function to x,y values.
+
+    Parameters
+    ----------
+    x : array-like
+    y : array-like
+
+    Other Parameters
+    ----------------
+    amplitude : float
+        The initial value for the amplitude
+    mean : float
+        The initial value for the mean of the gaussian function
+    stddev : float
+        The initial value for the standard deviation of the gaussian function
+    tied : dict
+    fixed : dict
+    bounds : dict
+        Parameters to be passed to the astropy models
+        http://docs.astropy.org/en/stable/api/astropy.modeling.functional_models.Gaussian1D.html
+
+    Returns
+    -------
+    g : function
+        The best-fit function, accepting x as input
+        and returning the best-fit model as output
+    """
+    g_in = models.Gaussian1D(amplitude=amplitude,mean=mean,stddev=stddev,
+                             tied=tied, fixed=fixed, bounds=bounds)
+    fit_g = fitting.LevMarLSQFitter()
+    g = fit_g(g_in, x, y)
+    return g


Test patch:
diff --git a/stingray/pulse/tests/test_modeling.py b/stingray/pulse/tests/test_modeling.py
new file mode 100644
index 00000000..5660da95
--- /dev/null
+++ b/stingray/pulse/tests/test_modeling.py
@@ -0,0 +1,84 @@
+import numpy as np
+from stingray.pulse.modeling import fit_sinc, fit_gaussian
+
+np.random.seed(0)
+
+
+def test_sinc_function():
+    x = np.linspace(-5., 5., 200)
+    y = 2 * (np.sin(x)/x)**2
+    y += np.random.normal(0., 0.1, x.shape)
+
+    s = fit_sinc(x, y)
+
+    assert np.abs(s.mean) < 0.1
+    assert np.abs(s.amplitude - 2) < 0.1
+    assert np.abs(s.width - 1) < 0.1
+
+
+def test_sinc_fixed():
+    x = np.linspace(-5., 5., 200)
+    y = 2 * (np.sin(x)/x)**2
+    y += np.random.normal(0., 0.1, x.shape)
+
+    sf = fit_sinc(x, y, mean=1., fixed={"mean": True, "amplitude": False})
+    assert sf.mean.fixed
+    assert not sf.amplitude.fixed
+
+
+def test_sinc_obs():
+    obs_length = 0.32
+    x = np.linspace(-5., 5., 200)
+    w = 1 / (np.pi*obs_length)
+    y = 2 * (np.sin(x / w) / (x / w))**2
+    y += np.random.normal(0., 0.1, x.shape)
+
+    s = fit_sinc(x, y, obs_length=obs_length )
+
+    assert np.abs(1 / (np.pi*obs_length) - s.width) < 0.1
+    assert s.width.fixed
+
+
+def test_gaussian_function():
+    x = np.linspace(-5., 5., 200)
+    y = 2 * np.exp(-0.5 * (x - 1.3)**2 / 0.7**2)
+    y += np.random.normal(0., 0.1, x.shape)
+
+    gs = fit_gaussian(x, y)
+
+    assert np.abs(gs.mean - 1.3) < 0.1
+    assert np.abs(gs.amplitude - 2) < 0.1
+    assert np.abs(gs.stddev - 0.7) < 0.1
+
+
+def test_gaussian_bounds():
+    x = np.linspace(-5., 5., 200)
+    y = 2 * np.exp(-0.5 * (x - 1.3)**2 / 0.7**2)
+    y += np.random.normal(0., 0.1, x.shape)
+
+    gs = fit_gaussian(x, y,
+                      bounds={"mean": [1., 1.6], "amplitude":[1.7, 2.3]})
+
+
+def test_gaussian_fixed():
+    x = np.linspace(-5., 5., 200)
+    y = 2 * np.exp(-0.5 * (x - 1.3)**2 / 0.7**2)
+    y += np.random.normal(0., 0.1, x.shape)
+
+    gs = fit_gaussian(x, y, mean=1.3, fixed={"mean": True, "amplitude": False})
+    assert gs.mean.fixed
+    assert not gs.amplitude.fixed
+
+
+def test_gaussian_tied():
+    x = np.linspace(-5., 5., 200)
+    y = 2 * np.exp(-0.5 * (x - 1.3)**2 / 0.7**2)
+    y += np.random.normal(0., 0.1, x.shape)
+
+    def tiedgaussian(model):
+        mean = model.amplitude / 2
+        return mean
+
+    gs = fit_gaussian(x, y, tied={"mean": tiedgaussian})
+
+    assert np.abs(gs.mean/gs.amplitude - 0.5) < 0.1