-Interpretable gaussian processes for stellar light curves using starry. +Interpretable Gaussian processes for stellar light curves using starry.
+ +# A Gaussian Process for Stellar Variability + +The `starry_process` code implements an interpretable Gaussian process (GP) +for modeling stellar light curves. Whether your goal is to marginalize +over the stellar variability signal (if you think of it as noise) +or to understand the surface features that generated it (if you +think of it as data), this code is for you. The GP implemented here works +just like any other GP you might already use in your analysis, except that +its hyperparameters are *physically interpretable*. These are (among others) +the **radius of the spots**, the +**mean and variance of the latitude distribution**, +the **spot contrast**, and the **number of spots**. Users can also specify +things like the rotational period of the star, the limb darkening parameters, +and the inclination (or marginalize over the inclination if it is not known). + +The code is written in Python and relies on the +[Theano package](https://theano-pymc.readthedocs.io/en/stable/index.html), +so a little familiarity with that is recommended. Check out the crash +course [here](https://luger.dev/starry_process/notebooks/Quickstart.html#Compiling-theano-functions). + +# Quickstart + +Import the main interface: + +```python +from starry_process import StarryProcess +``` + +Draw samples from a Gaussian process with small mid-latitude spots: + +```python +# Instantiate the GP +sp = StarryProcess( + r=10, # spot radius in degrees + mu=30, # central spot latitude in degrees + sigma=5, # latitude std. dev. in degrees + c=0.1, # fractional spot contrast + n=10 # number of spots +) + +# Draw & visualize a spherical harmonic sample +y = sp.draw_ylm().eval() +sp.visualize(y) + +# Compute & plot the flux at some inclination +t = np.linspace(0, 4, 1000) +flux = sp.flux(y, t, i=60).eval() +plt.plot(t, flux) +``` + +
+
+Same as above, but for high-latitude spots:
+
+```python
+sp = StarryProcess(r=10, mu=0, sigma=10, c=0.1, n=10)
+```
+
+
+
+Large equatorial spots:
+
+```python
+sp = StarryProcess(r=30, mu=0, sigma=10, c=0.1, n=10)
+```
+
+
+
+Small, approximately isotropic spots:
+
+```python
+sp = StarryProcess(r=10, mu=0, sigma=40, c=0.1, n=10)
+```
+
+
+
+For more information check out the full
+[Quickstart tutorial](https://luger.dev/starry_process/notebooks/Quickstart.html) and
+the complete [documentation](https://luger.dev).
+
+# References & Attribution
+
+The code is described in this
+[JOSS paper](https://github.com/rodluger/starry_process/raw/joss-paper/joss/paper.pdf).
+It is the backbone of the
+[Mapping Stellar Surfaces](https://github.com/rodluger/mapping_stellar_surfaces)
+paper series, including:
+
+ - [Degeneracies in the rotational light curve problem](https://github.com/rodluger/mapping_stellar_surfaces/raw/paper1-pdf/ms.pdf)
+ - [An interpretable Gaussian process model for stellar light curves](https://github.com/rodluger/mapping_stellar_surfaces/raw/paper2-pdf/ms.pdf)
+
+If you make use of this code in your research, please cite
+
+```
+@article{Luger2021a,
+ author = {{Luger}, Rodrigo and {Foreman-Mackey}, Daniel and
+ {Hedges}, Christina},
+ title = {{Mapping stellar surfaces I: Degeneracies in the rotational light curve problem}},
+ journal = {in preparation},
+ year = {2021},
+ url = {https://github.com/rodluger/mapping_stellar_surfaces/raw/paper1-pdf/ms.pdf}
+}
+```
+
+```
+@article{Luger2021b,
+ author = {{Luger}, Rodrigo and {Foreman-Mackey}, Daniel and
+ {Hedges}, Christina},
+ title = {{Mapping stellar surfaces II: An interpretable Gaussian process model for light curves}},
+ journal = {in preparation},
+ year = {2021},
+ url = {https://github.com/rodluger/mapping_stellar_surfaces/raw/paper2-pdf/ms.pdf}
+}
+```
+
+```
+@article{Luger2021c,
+ author = {{Luger}, Rodrigo and {Foreman-Mackey}, Daniel and
+ {Hedges}, Christina},
+ title = {{starry\_process: Interpretable Gaussian processes for stellar light curves}},
+ journal = {in preparation},
+ year = {2021},
+ month = {Jan},
+ url = {https://github.com/rodluger/starry_process/raw/joss-paper/joss/paper.pdf}
+}
+```
diff --git a/docs/generate_readme_plots.py b/docs/generate_readme_plots.py
new file mode 100644
index 0000000..a4d2d95
--- /dev/null
+++ b/docs/generate_readme_plots.py
@@ -0,0 +1,107 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib.colors import Normalize
+from mpl_toolkits.axes_grid1.inset_locator import inset_axes
+from starry_process import StarryProcess
+import starry
+import os
+
+
+def plot_samples():
+
+ # Settings
+ nsamples = 5
+ norm = Normalize(vmin=0.5, vmax=1.1)
+ incs = [15, 30, 45, 60, 75, 90]
+ t = np.linspace(0, 4, 1000)
+ cmap = plt.get_cmap("plasma_r")
+ color = lambda i: cmap(0.1 + 0.8 * i / (len(incs) - 1))
+ map = starry.Map(15, lazy=False)
+ kwargs = [
+ dict(r=10, a=0.40, b=0.27, c=0.1, n=10, seed=0),
+ dict(r=10, a=0.31, b=0.36, c=0.1, n=10, seed=1),
+ dict(r=30, a=0.28, b=0.0, c=0.05, n=10, seed=2),
+ dict(r=10, a=0.06, b=0.0, c=0.1, n=20, seed=3),
+ ]
+
+ # One figure per `kwargs`
+ for n in range(len(kwargs)):
+
+ # Set up the plot
+ fig, ax = plt.subplots(
+ 2,
+ nsamples + 1,
+ figsize=(12, 2.5),
+ gridspec_kw={
+ "height_ratios": np.array([1, 0.5]),
+ "width_ratios": np.append(np.ones(nsamples), 0.1),
+ },
+ )
+
+ # Draw samples
+ sp = StarryProcess(marginalize_over_inclination=False, **kwargs[n])
+ y = sp.sample_ylm(nsamples=nsamples).eval()
+
+ # Normalize so that the background photosphere
+ # has unit intensity (for plotting)
+ y[:, 0] += 1
+ y *= np.pi
+
+ # Visualize each sample
+ for k in range(nsamples):
+ map[:, :] = y[k]
+ map.show(ax=ax[0, k], projection="moll", norm=norm)
+ ax[0, k].set_ylim(-1.5, 2.25)
+ ax[0, k].set_rasterization_zorder(1)
+ for i, inc in enumerate(incs):
+ map.inc = inc
+ flux = map.flux(theta=360.0 * t)
+ flux -= np.mean(flux)
+ flux *= 1e3
+ ax[1, k].plot(t, flux, color=color(i), lw=0.75)
+ if k == 0:
+ ax[1, k].spines["top"].set_visible(False)
+ ax[1, k].spines["right"].set_visible(False)
+ ax[1, k].set_xlabel("rotations", fontsize=8)
+ ax[1, k].set_ylabel("flux [ppt]", fontsize=8)
+ ax[1, k].set_xticks([0, 1, 2, 3, 4])
+ for tick in (
+ ax[1, k].xaxis.get_major_ticks()
+ + ax[1, k].yaxis.get_major_ticks()
+ ):
+ tick.label.set_fontsize(6)
+ ax[1, k].tick_params(direction="in")
+ else:
+ ax[1, k].axis("off")
+
+ # Appearance tweaks
+ cax = inset_axes(
+ ax[0, -1], width="70%", height="50%", loc="lower center"
+ )
+ cbar = fig.colorbar(
+ ax[0, k].images[0], cax=cax, orientation="vertical"
+ )
+ cbar.set_label("intensity", fontsize=8)
+ cbar.set_ticks([0.5, 0.75, 1])
+ cbar.ax.tick_params(labelsize=6)
+ ax[0, -1].axis("off")
+ lax = inset_axes(
+ ax[1, -1], width="80%", height="100%", loc="center right"
+ )
+ for i, inc in enumerate(incs):
+ lax.plot(
+ 0, 0, color=color(i), lw=1, label=r"{}$^\circ$".format(inc)
+ )
+ lax.legend(loc="center left", fontsize=5, frameon=False)
+ lax.axis("off")
+ ax[1, -1].axis("off")
+ dy = max([max(np.abs(ax[1, k].get_ylim())) for k in range(nsamples)])
+ for k in range(nsamples):
+ ax[1, 0].set_ylim(-dy, dy)
+
+ # We're done
+ fig.savefig("samples_{}.png".format(n), bbox_inches="tight", dpi=100)
+
+
+if __name__ == "__main__":
+ plot_samples()
diff --git a/joss/figures/samples.py b/joss/figures/samples.py
index 4b89963..165a99e 100644
--- a/joss/figures/samples.py
+++ b/joss/figures/samples.py
@@ -14,8 +14,8 @@
cmap = plt.get_cmap("plasma_r")
color = lambda i: cmap(0.1 + 0.8 * i / (len(incs) - 1))
kwargs = [
- dict(r=10, a=0.40, b=0.27, c=0.15, N=10, seed=7),
- dict(r=15, a=0.31, b=0.36, c=0.075, N=10, seed=8),
+ dict(r=10, a=0.40, b=0.27, c=0.15, n=10, seed=7),
+ dict(r=15, a=0.31, b=0.36, c=0.075, n=10, seed=8),
]
map = starry.Map(15, lazy=False)
diff --git a/starry_process/sp.py b/starry_process/sp.py
index 56c66fb..1ea4193 100644
--- a/starry_process/sp.py
+++ b/starry_process/sp.py
@@ -1009,8 +1009,8 @@ def flux(
as ``sample_ylm`` can thus be directly passed to this method.
t (vector): The time array in arbitrary units.
i (scalar, optional): The inclination of the star in degrees.
- Default is %%defaults["i"]%%. If ``marginalize_over_inclination``
- is set, this argument is ignored.
+ Default is %%defaults["i"]%%. This option is accepted even
+ if ``marginalize_over_inclination`` is ``True``.
p (scalar, optional): The rotational period of the star in the same
units as ``t``. Default is %%defaults["p"]%%.
u (vector, optional): The limb darkening coefficients for the