Merge bd7ce9e into a86327a

exoplanet/gp/celerite_test.py

@@ -153,7 +153,4 @@ def test_integrated(seed=1234):
     kernel = terms.IntegratedTerm(kernel, dt)
     a = kernel.get_celerite_matrices(x, diag)[0].eval()
     k0 = kernel.value(tt.zeros(1)).eval()
-    print(a, k0 + diag)
     assert np.allclose(a, k0 + diag)
-
-    assert 0

exoplanet/gp/terms.py

@@ -356,12 +356,12 @@ def get_coefficients(self):
         dd = d * self.delta
         c2 = c ** 2
         d2 = d ** 2
-        factor = 1.0 / (self.delta * (c2 + d2)) ** 2
+        factor = 2.0 / (self.delta * (c2 + d2)) ** 2
         cos_term = tt.cosh(cd) * tt.cos(dd) - 1
         sin_term = tt.sinh(cd) * tt.sin(dd)
 
-        C1 = 2 * (a * (c2 - d2) + 2 * b * c * d)
-        C2 = 2 * (b * (c2 - d2) - 2 * a * c * d)
+        C1 = a * (c2 - d2) + 2 * b * c * d
+        C2 = b * (c2 - d2) - 2 * a * c * d
 
         coeffs += [
             factor * (C1 * cos_term - C2 * sin_term),

paper/exoplanet.tex

@@ -807,22 +807,23 @@ \section{Exposure time integration for celerite models}
 \end{align}
 This results in a time-integrated kernel for $\tau \ge \Delta$ given by
 \begin{equation}
-k_\Delta(\tau) = A\exp\left(-c\,\tau\right)\cos\left(d\,\tau\right) + B\exp\left(-c\,\tau\right)\sin\left(d\,\tau\right).
+k_{\Delta\le\tau}(\tau) = A\exp\left(-c\,\tau\right)\cos\left(d\,\tau\right) + B\exp\left(-c\,\tau\right)\sin\left(d\,\tau\right).
 \end{equation}
 while for $0 \le \tau < \Delta$,
-\begin{align}
-k_\Delta(\tau) &= A\exp\left(-c\,\tau\right)\cos\left(d\,\tau\right) + B\exp\left(-c\,\tau\right)\sin\left(d\,\tau\right) \nonumber\\
-&+  \frac{2}{ \Delta^2(c^2+d^2)^2} \left\{(ac+bd)(c^2+d^2)(\Delta-\tau)\right. \nonumber \\
-&+\left.  \sinh(c(\tau-\Delta)) \left(C_1\cos(d(\Delta-\tau)) +C_2 \sin(d(\Delta-\tau))\right) \right\}.
-\end{align}
+\begin{eqnarray}
+k_{\Delta>\tau}(\tau) &=& k_{\Delta\le\tau}(\tau) + 2\,\left[
+  - C_1\,\sinh(c\,\Delta-c\,\tau)\,\cos(d\,\Delta-d\,\tau) \right. \nonumber\\
+&&  + C_2\,\cosh(c\,\Delta-c\,\tau)\,\sin(d\,\Delta-d\,\tau) \nonumber \\
+&&  \left. + (a\,c + b\,d)\,(c^2 + d^2)\,(\Delta-\tau)\right] / [\Delta\,(c^2+d^2)]^2
+\end{eqnarray}
 Note that in the limit of $b = d = 0$, this kernel reverts to the real kernel.
 As with the real case, the kernel is only semi-separable when exposures
 do not overlap in time, i.e. $\vert t_i - t_j \vert > \Delta$ for all $i,j$, while the diagonal
 components of the covariance matrix require adding the extra term:
-\begin{equation}
-\frac{2}{ \Delta^2(c^2+d^2)^2} \left\{(ac+bd)(c^2+d^2)\Delta
--\sinh(c\Delta) \left(C_1\cos(d\Delta) +C_2 \sin(d\Delta)\right) \right\}.
-\end{equation}
+\begin{eqnarray}
+  \frac{2\,\left[
+    - C_1\,\sinh(c\,\Delta)\,\cos(d\,\Delta) + C_2\,\cosh(c\,\Delta)\,\sin(d\,\Delta) + (a\,c + b\,d)\,(c^2 + d^2)\,\Delta\right]}{[\Delta\,(c^2+d^2)]^2}
+\end{eqnarray}
 As with the real case, the power spectrum of the time-integrated kernel simply
 requires multiplication by the sinc-squared function:
 \begin{equation}

paper/proofs/celerite-integral.ipynb

@@ -689,6 +689,42 @@
     "sm.simplify(ktest - res)"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "This corresponds to adding a function to the kernel:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 27,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "0"
+      ]
+     },
+     "execution_count": 27,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "delta = sm.simplify(ktest.expand(trig=True) - (A * sm.exp(-c*tau) * sm.cos(d*tau) + B * sm.exp(-c*tau) * sm.sin(d*tau)))\n",
+    "\n",
+    "C1 = (a*c**2 - a*d**2 + 2*b*c*d)\n",
+    "C2 = (b*c**2 - b*d**2 - 2*a*c*d)\n",
+    "sinh = (sm.exp(c*(dt - tau)) - sm.exp(-c*(dt - tau)))/2\n",
+    "cosh = (sm.exp(c*(dt - tau)) + sm.exp(-c*(dt - tau)))/2\n",
+    "norm = dt**2 * (c**2 + d**2)**2\n",
+    "\n",
+    "delta_test = 2 * (C2 * cosh * sm.sin(d*(dt-tau)) - C1 * sinh * sm.cos(d*(dt-tau)) + (a*c + b*d) * (dt-tau) * (c**2 + d**2)) / norm\n",
+    "sm.simplify(delta.expand(trig=True) - delta_test.expand(trig=True))"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -709,7 +745,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 27,
+   "execution_count": 28,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -718,7 +754,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 28,
+   "execution_count": 29,
    "metadata": {},
    "outputs": [
     {
@@ -727,7 +763,7 @@
        "(-a*c**2*exp(2*c*dt) + a*c**2 + a*d**2*exp(2*c*dt) - a*d**2 - 2*b*c*d*exp(2*c*dt) + 2*b*c*d)*exp(-c*dt)/(dt**2*(c**4 + 2*c**2*d**2 + d**4))"
       ]
      },
-     "execution_count": 28,
+     "execution_count": 29,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -738,7 +774,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 29,
+   "execution_count": 30,
    "metadata": {},
    "outputs": [
     {
@@ -747,7 +783,7 @@
        "(-2*a*c*d*exp(2*c*dt) - 2*a*c*d + b*c**2*exp(2*c*dt) + b*c**2 - b*d**2*exp(2*c*dt) - b*d**2)*exp(-c*dt)/(dt**2*(c**4 + 2*c**2*d**2 + d**4))"
       ]
      },
-     "execution_count": 29,
+     "execution_count": 30,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -758,7 +794,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 30,
+   "execution_count": 31,
    "metadata": {},
    "outputs": [
     {
@@ -767,7 +803,7 @@
        "2*(a*c + b*d)/(dt*(c**2 + d**2))"
       ]
      },
-     "execution_count": 30,
+     "execution_count": 31,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -778,7 +814,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 31,
+   "execution_count": 32,
    "metadata": {},
    "outputs": [
     {
@@ -787,7 +823,7 @@
        "0"
       ]
      },
-     "execution_count": 31,
+     "execution_count": 32,
      "metadata": {},
      "output_type": "execute_result"
     }