MNGuenther / allesfitter

allesfitter is a convenient wrapper around the packages ellc (light curve and RV models), dynesty (static and dynamic nested sampling) emcee (Markov Chain Monte Carlo sampling) and celerite (Gaussian Process models).
MIT License
62 stars 36 forks source link

Quasi Periodic Kernal? #55

Open aspolanski opened 2 years ago

aspolanski commented 2 years ago

Just wondering if there's any plans to implement a QP kernal for stellar variability? It's become a fairly popular option for RV analysis

MNGuenther commented 2 years ago

Hi @aspolanski,

Thanks for this interesting feature request! (And thanks for your patience, you caught me right on my first day of holidays.)

Allesfitter's GP kernels are drawn from celerite's default selection for now (Real, Complex, SHO, Matern 3/2). Out of genuine curiosity: Is the QP (or QPC) kernel behaving significantly different from the SHO kernel in your experience?

Would you have a formulation of your desired QP (or QPC) kernel that is conform with celerite (i.e. a coded up line to create a new celerite kernel)? Then I can plug it into allesfiter for you!

Cheers Max

szunigaf commented 1 year ago

Hi @MNGuenther

I found in celerite2 a term that can be used for stellar rotation, and indeed is a combination of two SHO terms that maybe can be integrated in allesfitter.

https://celerite2.readthedocs.io/en/latest/api/python/#celerite2.terms.RotationTerm

class RotationTerm(TermSum):
    r"""A mixture of two SHO terms that can be used to model stellar rotation

    This term has two modes in Fourier space: one at ``period`` and one at
    ``0.5 * period``. This can be a good descriptive model for a wide range of
    stochastic variability in stellar time series from rotation to pulsations.

    More precisely, the parameters of the two :class:`SHOTerm` terms are
    defined as

    .. math::

        Q_1 = 1/2 + Q_0 + \delta Q \\
        \omega_1 = \frac{4\,\pi\,Q_1}{P\,\sqrt{4\,Q_1^2 - 1}} \\
        S_1 = \frac{\sigma^2}{(1 + f)\,\omega_1\,Q_1}

    for the primary term, and

    .. math::

        Q_2 = 1/2 + Q_0 \\
        \omega_2 = \frac{8\,\pi\,Q_1}{P\,\sqrt{4\,Q_1^2 - 1}} \\
        S_2 = \frac{f\,\sigma^2}{(1 + f)\,\omega_2\,Q_2}

    for the secondary term.

    Args:
        sigma: The standard deviation of the process.
        period: The primary period of variability.
        Q0: The quality factor (or really the quality factor minus one half;
            this keeps the system underdamped) for the secondary oscillation.
        dQ: The difference between the quality factors of the first and the
            second modes. This parameterization (if ``dQ > 0``) ensures that
            the primary mode alway has higher quality.
        f: The fractional amplitude of the secondary mode compared to the
            primary. This should probably always be ``0 < f < 1``, but that
            is not enforced.
    """

    @staticmethod
    def get_test_parameters():
        return dict(sigma=1.5, period=3.45, Q0=1.3, dQ=1.05, f=0.5)

    def __init__(self, *, sigma, period, Q0, dQ, f):
        self.sigma = float(sigma)
        self.period = float(period)
        self.Q0 = float(Q0)
        self.dQ = float(dQ)
        self.f = float(f)

        self.amp = self.sigma**2 / (1 + self.f)

        # One term with a period of period
        Q1 = 0.5 + self.Q0 + self.dQ
        w1 = 4 * np.pi * Q1 / (self.period * np.sqrt(4 * Q1**2 - 1))
        S1 = self.amp / (w1 * Q1)

        # Another term at half the period
        Q2 = 0.5 + self.Q0
        w2 = 8 * np.pi * Q2 / (self.period * np.sqrt(4 * Q2**2 - 1))
        S2 = self.f * self.amp / (w2 * Q2)

        super().__init__(
            SHOTerm(S0=S1, w0=w1, Q=Q1), SHOTerm(S0=S2, w0=w2, Q=Q2)
        )