transit tutorial and defaults

docs/_static/notebooks/pymc3-extras.ipynb

@@ -107,10 +107,7 @@
     "    pm.MvNormal(\"x\", mu=np.zeros(ndim), chol=L, shape=(ndim,))\n",
     "    \n",
     "    # Run the burn-in and learn the mass matrix\n",
-    "    # For this specific problem, we get better performance if\n",
-    "    # we don't regularize the variance estimate:\n",
-    "    step_kwargs = dict(regular_window=0,\n",
-    "                       target_accept=0.9)\n",
+    "    step_kwargs = dict(target_accept=0.9)\n",
     "    sampler.tune(tune=2000, step_kwargs=step_kwargs)\n",
     "    \n",
     "    # Run the production chain\n",

docs/_static/notebooks/stellar-variability.ipynb

@@ -16,7 +16,7 @@
     "# Gaussian process models for stellar variability\n",
     "\n",
     "When fitting exoplanets, we also need to fit for the stellar variability and Gaussian Processes (GPs) are often a good descriptive model for this variation.\n",
-    "[*PyMC3* has support for all sorts of general GP models](https://docs.pymc.io/gp.html), but *exoplanet* includes support for scalable 1D GPs (see :ref:`gp` for more info) that can work with large datasets.\n",
+    "[PyMC3 has support for all sorts of general GP models](https://docs.pymc.io/gp.html), but *exoplanet* includes support for scalable 1D GPs (see :ref:`gp` for more info) that can work with large datasets.\n",
     "In this tutorial, we go through the process of modeling the light curve of a rotating star observed by Kepler using *exoplanet*.\n",
     "\n",
     "First, let's download and plot the data:"

docs/index.rst

@@ -5,8 +5,8 @@ User guide
 
 .. toctree::
    :maxdepth: 2
-   :caption: User Documentation
 
+   user/what
    user/install
    tutorials/citation
    user/api
@@ -21,6 +21,7 @@ Tutorials
    tutorials/intro-to-pymc3
    tutorials/pymc3-extras
    tutorials/rv
+   tutorials/transit
    tutorials/gp
    tutorials/stellar-variability
 

docs/user/what.rst

@@ -0,0 +1,8 @@
+.. _what:
+
+What is exoplanet?
+==================
+
+*exoplanet* is a set of tools that provide the foundation for using
+gradient-based inference algorithms for fitting exoplanets observed using
+radial velocities or transits.

exoplanet/light_curves.py

@@ -41,7 +41,7 @@ def __init__(self, u, r_star=1.0, model=None):
 
     def get_light_curve(self, orbit=None, r=None, t=None,
                         texp=None, oversample=7, order=0,
-                        use_approx_in_transit=True, duration_factor=3):
+                        use_in_transit=True):
         """Get the light curve for an orbit at a set of times
 
         Args:
@@ -73,10 +73,9 @@ def get_light_curve(self, orbit=None, r=None, t=None,
                 centered Riemann sum (equivalent to the "resampling" procedure
                 suggested by Kipping 2010), ``1`` for the trapezoid rule, or
                 ``2`` for Simpson's rule.
-            use_approx_in_transit (Optional[bool]): If ``True``, the model will
-                only be evaluated for the data points expected to be in transit
-                as computed using the ``approx_in_transit`` method on
-                ``orbit``.
+            use_in_transit (Optional[bool]): If ``True``, the model will only
+                be evaluated for the data points expected to be in transit
+                as computed using the ``in_transit`` method on ``orbit``.
 
         """
         if orbit is None:
@@ -90,10 +89,10 @@ def get_light_curve(self, orbit=None, r=None, t=None,
         r = tt.reshape(r, (r.size,))
         t = tt.as_tensor_variable(t)
 
-        if use_approx_in_transit:
-            transit_model = tt.zeros_like(t)
-            inds = orbit.approx_in_transit(t, r=r, texp=texp,
-                                           duration_factor=duration_factor)
+        if use_in_transit:
+            transit_model = tt.shape_padleft(tt.zeros_like(r), t.ndim) \
+                + tt.shape_padright(tt.zeros_like(t), r.ndim)
+            inds = orbit.in_transit(t, r=r, texp=texp)
             t = t[inds]
 
         if texp is None:
@@ -142,9 +141,8 @@ def get_light_curve(self, orbit=None, r=None, t=None,
         if texp is not None:
             stencil = tt.shape_padright(tt.shape_padleft(stencil, t.ndim), 1)
             lc = tt.sum(stencil * lc, axis=t.ndim)
-        lc = tt.sum(lc, axis=-1)
 
-        if use_approx_in_transit:
+        if use_in_transit:
             transit_model = tt.set_subtensor(transit_model[inds], lc)
             return transit_model
         else:

exoplanet/light_curves_test.py

@@ -44,7 +44,7 @@ def test_light_curve_grad():
     utt.verify_grad(lc, [u_val, b_val, r_val])
 
 
-def test_approx_in_transit():
+def test_in_transit():
     t = np.linspace(-20, 20, 1000)
     m_planet = np.array([0.3, 0.5])
     m_star = 1.45
@@ -62,13 +62,12 @@ def test_approx_in_transit():
 
     lc = StarryLightCurve(u)
     model1 = lc.get_light_curve(r=r, orbit=orbit, t=t)
-    model2 = lc.get_light_curve(r=r, orbit=orbit, t=t,
-                                use_approx_in_transit=False)
+    model2 = lc.get_light_curve(r=r, orbit=orbit, t=t, use_in_transit=False)
     vals = theano.function([], [model1, model2])()
     utt.assert_allclose(*vals)
 
     model1 = lc.get_light_curve(r=r, orbit=orbit, t=t, texp=0.1)
     model2 = lc.get_light_curve(r=r, orbit=orbit, t=t, texp=0.1,
-                                use_approx_in_transit=False)
+                                use_in_transit=False)
     vals = theano.function([], [model1, model2])()
     utt.assert_allclose(*vals)

exoplanet/orbits/keplerian.py

@@ -353,20 +353,16 @@ def get_radial_velocity(self, t, K=None, output_units=None):
         v = self.get_star_velocity(t)
         return conv * v[2]
 
-    def approx_in_transit(self, t, r=0.0, texp=None, duration_factor=3):
-        """Get a list of timestamps that are expected to be in transit
+    def in_transit(self, t, r=0.0, texp=None):
+        """Get a list of timestamps that are in transit
 
         Args:
             t (vector): A vector of timestamps to be evaluated.
             r (Optional): The radii of the planets.
             texp (Optional[float]): The exposure time.
-            duration_factor (Optional[float]): The factor by which to multiply
-                the approximate duration when computing the in transit points.
-                Larger values will be more conservative and might be needed for
-                large planets or very eccentric orbits.
 
         Returns:
-            The indices of the timestamps that are expected to be in transit.
+            The indices of the timestamps that are in transit.
 
         """
         z = tt.zeros_like(self.a)
@@ -395,26 +391,6 @@ def approx_in_transit(self, t, r=0.0, texp=None, duration_factor=3):
 
         return tt.arange(t.size)[mask]
 
-        # Estimate the maximum duration of the transit using the equations
-        # from Winn (2010)
-        # arg = (self.r_star + r) / (-self.a_planet)
-        # max_dur = self.period * tt.arcsin(arg) / np.pi
-        # if self.ecc is not None:
-        #     max_dur *= tt.sqrt(1-self.ecc**2) / (1+self.ecc*self.sin_omega)
-
-        # # Wrap the times into time since transit
-        # hp = 0.5 * self.period
-        # dt = tt.mod(self._warp_times(t) - self.t0 + hp, self.period) - hp
-
-        # # Estimate the data points that are within the maximum duration of the
-        # # transit
-        # tol = 0.5 * duration_factor * max_dur
-        # if texp is not None:
-        #     tol += 0.5 * texp
-        # mask = tt.any(tt.abs_(dt) < tol, axis=-1)
-
-        # return tt.arange(t.size)[mask]
-
 
 def get_true_anomaly(M, e, **kwargs):
     """Get the true anomaly for a tensor of mean anomalies and eccentricities

exoplanet/orbits/keplerian_test.py

@@ -107,7 +107,7 @@ def test_velocity():
     utt.assert_allclose(planet_vel, planet_vel_expect)
 
 
-def test_approx_in_transit():
+def test_in_transit():
     t = np.linspace(-20, 20, 1000)
     m_planet = np.array([0.3, 0.5])
     m_star = 1.45
@@ -125,7 +125,7 @@ def test_approx_in_transit():
     r_pl = np.array([0.1, 0.03])
     coords = theano.function([], orbit.get_relative_position(t))()
     r2 = coords[0]**2 + coords[1]**2
-    inds = theano.function([], orbit.approx_in_transit(t, r=r_pl))()
+    inds = theano.function([], orbit.in_transit(t, r=r_pl))()
 
     m = np.isin(np.arange(len(t)), inds)
     in_ = r2[inds] <= ((r_star + r_pl)**2)[None, :]

exoplanet/orbits/simple.py

@@ -69,11 +69,9 @@ def get_star_velocity(self, t):
     def get_radial_velocity(self, t, output_units=None):
         raise NotImplementedError("a SimpleTransitOrbit has no velocity")
 
-    def approx_in_transit(self, t, r=None, texp=None, duration_factor=None):
+    def in_transit(self, t, r=None, texp=None):
         """Get a list of timestamps that are in transit
 
-        For this orbit, this function is exact.
-
         Args:
             t (vector): A vector of timestamps to be evaluated.
             r (Optional): The radii of the planets.

exoplanet/sampling.py

@@ -42,7 +42,7 @@ def __init__(self, start=75, finish=50, window=25, dense=True):
         self._current_trace = None
 
     def get_step_for_trace(self, trace=None, model=None,
-                           regular_window=5, regular_variance=1e-3,
+                           regular_window=0, regular_variance=1e-3,
                            **kwargs):
         """Get a PyMC3 NUTS step tuned for a given burn-in trace
 
@@ -77,14 +77,17 @@ def get_step_for_trace(self, trace=None, model=None,
             if self.dense:
                 # Compute the regularized sample covariance
                 cov = np.cov(samples, rowvar=0)
-                cov = cov * N / (N + regular_window)
-                cov[np.diag_indices_from(cov)] += \
-                    regular_variance * regular_window / (N + regular_window)
+                if regular_window > 0:
+                    cov = cov * N / (N + regular_window)
+                    cov[np.diag_indices_from(cov)] += \
+                        regular_variance * regular_window / (N+regular_window)
                 potential = quad.QuadPotentialFull(cov)
             else:
                 var = np.var(samples, axis=0)
-                var = var * N / (N + regular_window)
-                var += regular_variance * regular_window / (N + regular_window)
+                if regular_window > 0:
+                    var = var * N / (N + regular_window)
+                    var += \
+                        regular_variance * regular_window / (N+regular_window)
                 potential = quad.QuadPotentialDiag(cov)
 
         return pm.NUTS(potential=potential, **kwargs)