Note

This tutorial was generated from a Jupyter notebook that can be downloaded here.

Linear solve

In this notebook we discuss how to linearly solve for the posterior over spherical harmonic coefficients of a map given a light curve. This is similar to what we did in the Eclipsing Binary notebook. The idea is to take advantage of the linearity of the starry solution to analytically compute the posterior over maps consistent with the data.

import matplotlib.pyplot as plt
import numpy as np
import starry

np.random.seed(12)
starry.config.lazy = False
starry.config.quiet = True

We’re going to demonstrate the linear solve feature using a map in reflected light, since the presence of a day/night terminator breaks many degeneracies and makes the mapping problem much less ill-posed. Let’s begin by instantiating a reflected light map of the Earth. We’ll give it the same obliquity as the Earth and observe it at an inclination of 60 degrees. We’ll also give it an amplitude of 1.3: this is the value that will scale the light curve (and which we will infer below).

map = starry.Map(ydeg=10, reflected=True)
map.load("earth")
map.amp = 1.3
map.obl = 23.5
map.inc = 60
map.show(projection="rect", illuminate=False)

Now we generate a dataset. We’ll assume we have 10,000 observations over the course of a full orbit of the planet. We further take the planet’s rotation period to be one-tenth of its orbital period. This will give us good coverage during all seasons, maximizing the amount of data we have for all the different regions of the planet:

# Make the planet rotate 10 times over one full orbit
npts = 10000
nrot = 10
time = np.linspace(0, 1, npts)
theta = np.linspace(0, 360 * nrot, npts)

# Position of the star relative to the planet in the orbital plane
t = np.reshape(time, (1, -1))
p = np.vstack((np.cos(2 * np.pi * t), np.sin(2 * np.pi * t), 0 * t))

# Rotate to an observer inclination of 60 degrees
ci = np.cos(map.inc * np.pi / 180)
si = np.sin(map.inc * np.pi / 180)
R = np.array([[1, 0, 0], [0, ci, -si], [0, si, ci]])
xs, ys, zs = R.dot(p)

# Keywords to the `flux` method
kwargs = dict(theta=theta, xs=xs, ys=ys, zs=zs)

# Compute the flux
flux0 = map.flux(**kwargs)
sigma = 0.005
flux = flux0 + sigma * np.random.randn(npts)

fig, ax = plt.subplots(1, figsize=(12, 4))
ax.plot(time, flux)
ax.set_xlabel("orbital phase", fontsize=18)
ax.set_ylabel("flux", fontsize=18);

Now the fun part. Let’s instantiate a new map so we can do inference on this dataset:

map = starry.Map(ydeg=10, reflected=True)
map.obl = 23.5
map.inc = 60

We now set the data vector (the flux and the covariance matrix):

map.set_data(flux, C=sigma ** 2)

And now the prior, which is also described as a multivariate gaussian with a mean \(\mu\) and covariance \(\Lambda\). This is the prior on the amplitude-weighted spherical harmonic coefficients. In other words, if \(\alpha\) is the map amplitude (map.amp) and \(y\) is the vector of spherical harmonic coefficients (map.y), we are placing a prior on the quantity \(x \equiv \alpha y\). While this may be confusing at first, recall that the coefficient of the \(Y_{0,0}\) harmonic is fixed at unity in starry, so we can’t really solve for it. But we can solve for all elements of the vector \(x\). Once we have the posterior for \(x\), we can easily obtain both the amplitude (equal to \(x_0\)) and the spherical harmonic coefficient vector (equal to \(x / x_0\)). This allows us to simultaneously obtain both the amplitude and the coefficients using a single efficient linear solve.

For the mean \(\mu\), we’ll set the first coefficient to unity (since we expect \(\alpha\) to be distributed about one, if the light curve is properly normalized) and all the others to zero (since our prior is isotropic and we want to enforce some regularization to keep the coefficients small).

For the covariance \(\Lambda\) (L in the code), we’ll make it a diagonal matrix. The first diagonal entry is the prior variance on the amplitude of the map, and we’ll set that to \(10^{-1}\). The remaining entries are the prior variance on the amplitude-weighted coefficients. We’ll pick something small—\(5 \times 10^{-4}\)—to keep things well regularized. In practice, this prior is related to our beliefs about the angular power spectrum of the map.

mu = np.empty(map.Ny)
mu[0] = 1
mu[1:] = 0
L = np.empty(map.Ny)
L[0] = 1e-1
L[1:] = 5e-4
map.set_prior(L=L)

Finally, we call solve, passing in the kwargs from before. In this case, we’re assuming we know the orbital information exactly. (When this is not the case, we need to do sampling for the orbital parameters; we cover this in more detail in the Eclipsing Binary tutorial).

%%time
x, cho_cov = map.solve(**kwargs)

CPU times: user 503 ms, sys: 75.4 ms, total: 578 ms
Wall time: 3.18 s

The values returned are the posterior mean of the amplitude-weighted spherical harmonic coefficients mu and the Cholesky factorization of the posterior covariance matrix cho_cov. To view our mean map, we can do what we described above:

map.amp = x[0]
map[1:, :] = x[1:] / x[0]

map.show(projection="rect", illuminate=False)

map.amp

array(1.30265515)

(Note that we recovered the correct amplitude, with small error). We can also draw a random sample from the posterior (which automatically sets the map amplitude and coefficients) by calling

map.draw()

map.show(projection="rect", illuminate=False)

map.amp

array(1.31245808)

We can verify that we got a good fit to the data from this random sample:

fig, ax = plt.subplots(1, figsize=(12, 4))
ax.plot(time, flux)
plt.plot(time, map.flux(**kwargs))
ax.set_xlabel("Orbital phase", fontsize=18)
ax.set_ylabel("Normalized flux", fontsize=18);