better readme

.github/workflows/docs.yml

@@ -57,6 +57,12 @@ jobs:
         run: |
           make -C docs html
 
+      - name: Generate plots for README
+        shell: bash -l {0}
+        run: |
+          python docs/generate_readme_plots.py
+          mv docs/samples_*.png docs/_build/html/_images
+
       - name: Push
         if: ${{ github.event_name != 'pull_request' }}
         run: |

README.md

@@ -23,5 +23,132 @@
 </p>
 
 <p align="center">
-Interpretable gaussian processes for stellar light curves using <a href="https://github.com/rodluger/starry">starry</a>.
+Interpretable Gaussian processes for stellar light curves using <a href="https://github.com/rodluger/starry">starry</a>.
 </p>
+
+# 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)
+```
+
+<img src="https://github.com/rodluger/starry_process/raw/gh-pages/_images/samples_0.png"/>
+
+Same as above, but for high-latitude spots:
+
+```python
+sp = StarryProcess(r=10, mu=0, sigma=10, c=0.1, n=10)
+```
+
+<img src="https://github.com/rodluger/starry_process/raw/gh-pages/_images/samples_1.png"/>
+
+Large equatorial spots:
+
+```python
+sp = StarryProcess(r=30, mu=0, sigma=10, c=0.1, n=10)
+```
+
+<img src="https://github.com/rodluger/starry_process/raw/gh-pages/_images/samples_2.png"/>
+
+Small, approximately isotropic spots:
+
+```python
+sp = StarryProcess(r=10, mu=0, sigma=40, c=0.1, n=10)
+```
+
+<img src="https://github.com/rodluger/starry_process/raw/gh-pages/_images/samples_3.png"/>
+
+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}
+}
+```

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()

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)

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